Wavelet representation of hardcore bosons

We consider the one-dimensional Tonks–Girardeau gas with a space-dependent potential out of equilibrium. We derive the exact dynamics of the system when divided into n boxes and decomposed into energy eigenstates within each box. It is a representation of the wave function that is a mixture between real space and momentum space, with basis elements consisting of plane waves localized in a box, giving rise to the term ‘wavelet’. Using this representation, we derive the emergence of generalized hydrodynamics in appropriate limits without assuming local relaxation. We emphasize that a generalized hydrodynamic behaviour emerges in a high-momentum and short-time limit, in addition to the more common large-space and late-time limit, which is akin to a semi-classical expansion. In this limit, conserved charges do not require numerous particles to be described by generalized hydrodynamics. We also show that this wavelet representation provides an efficient numerical algorithm for a complete description of the out-of-equilibrium dynamics of hardcore bosons.


Introduction
Computing dynamics in many-body quantum systems is a difficult task even in systems that can be diagonalized exactly, since it involves going from a real space basis (to specify the initial state or to compute expectation values) to an energy eigenstate basis (to time-evolve the state).In several areas of physics, hydrodynamics is known to provide a good approximation of the dynamics of a many-body system in certain conditions [1]- [4].Hydrodynamics postulates that the system can be divided into "fluid cells" that are relaxed to an equilibrium state, which in a microcanonical ensemble (or generalized microcanonical if there are several conserved charges) can be taken to be an energy eigenstate.At a larger scale, the system is thus described by equilibrium parameters (particle density, energy density, conserved charges, etc) that vary with space and time.Hydrodynamics is thus an approximate representation of a state that is mixed between real space (through the division into fluid cells) and energy eigenstates (through the description in terms of equilibrium parameters within each fluid cell).This mixed representation is key to its efficiency, as the real space resolution allows for a simple encoding of the initial state, while the energy eigenstate resolution allows for describing complex states locally with few parameters.
Recently, this hydrodynamic description has proven to be particularly efficient to describe 1D quantum integrable models [5], [6], where due to the infinite number of conservation laws it is called "generalized" hydrodynamics (GHD).A paradigmatic example is a gas of bosons in δ interaction, called the Lieb-Liniger model [7], which is known to provide an excellent description of cold atoms confined in a 1D geometry [8]- [14].The limit of infinite coupling is particularly studied and is called the Tonks-Girardeau gas of hardcore bosons [15], [16].But although exactly solvable with the Bethe ansatz, computing dynamics analytically or numerically in the Lieb-Liniger model remained a challenge for a long time.Contrary to the Bethe ansatz solution, this hydrodynamics description is amenable to efficient numerical implementation [17]- [22], yielding useful comparison material that was unavailable before.Besides its convenient practical use, it has also been reported to match remarkably well experimental [14], [23], [24] and numerical data [6], [19], [25], [26], see [27] for a review.
In this paper, inspired by the efficiency of hydrodynamics, we introduce an exact representation of hardcore bosons that is mixed between real space and momentum space.We explicitly implement the division of the system into fluid cells, and decompose a state in the basis given by the tensor product of energy eigenstates within each fluid cell.As the basis elements of this decomposition are small plane waves localized in each fluid cell, we call this decomposition a "wavelet" representation.Contrary to a coarse-graining or a hydrodynamic description, we do not remove any degree of freedom.Using the Girardeau mapping of hardcore bosons to free fermions [16], we derive the exact dynamics of the model in this wavelet representation.
The output of this work is twofold.Firstly, (i) we show the emergence of (generalized) hydrodynamics in certain limits in absence of potentials, without assuming local relaxation of the gas within the fluid cells.The usual limit in which hydrodynamics holds is the Euler scale limit, namely large space and large time at fixed ratio space over time.In this scaling each fluid cell becomes infinitely large.But we show in particular that a hydrodynamic behaviour also emerges in a short-time, high-momentum limit, which does not require any space rescaling.This is akin to a semi-classical expansion → 0. In this limit, every particle with momentum α initialised in a fluid cell moves "classically" from fluid cell to fluid cell with velocity 2α.In this scaling each fluid cell becomes infinitely small.Importantly, several observables (like particle density) do not require a large number of particles to be described by GHD, which is compatible with observations that GHD works for very low numbers of particles [14].Despite this observation being very simple, it is appealing since much more physically applicable to coldatom gases, and does not seem to have been observed elsewhere.To the best of our knowledge, previous mentions of the limit → 0 in the context of GHD were always accompanied with the thermodynamic limit N ∼ 1/ → ∞ [28], [29].
Secondly, (ii) we show that this wavelet representation provides an efficient way of simulating hardcore bosons with inhomogeneous fields.Assuming that the initial state is a tensor product of energy eigenstates on each fluid cell, the state possesses a determinant structure at all times in the wavelet decomposition which can be computed efficiently.As said above, this representation is mixed between real space (through the division into fluid cells) and energy eigenstates space (through the decomposition into the local eigenstates in each fluid cell).This wavelet representation allows for a truncation of the higher momenta in each fluid cell, which is a physically meaningful truncation.We find that keeping a very small number of momenta in each fluid cell already gives very good results.Besides, this wavelet representation allows for a simple formula for the single-particle Green's function, which is usually difficult to compute [30]- [35].As an illustration, we present simulations of quantum-Newton-cradle-like protocols for hardcore bosons with different confining potentials [36]- [38].Some comments about GHD are in order, the hydrodynamic theory that describes the Lieb-Liniger model [5], [6].Hydrodynamics is here "generalized" in the sense that the Lieb-Liniger model being integrable, it has several other conserved charges beyond particle density, momentum and energy, encoded in a generating function ρ(λ) called root density.The derivation of hydrodynamics in physical systems from first principles is generally a difficult task.The standard approach to argue for GHD is to divide the system into a collection of fluid cells assumed to be at equilibrium [39], [40].Then, one computes the expectation value of the current of the conserved charges that parametrize the equilibrium state [5], [6], [41]- [43].Popular toy models for these currents between fluid cells are bi-partitioning protocols [44]- [49].This current leads to a change in the conserved quantities of the neighbouring fluid cells, yielding partial differential equations holding at a much larger scale than that of the fluid cells.Although intuitive, the drawback of this approach is that it entirely relies on the assumption of local equilibrium in each fluid cell, which is the difficult point to show.We also note that observations of large spatial correlations have challenged the fluid cell picture recently [50].
The hardcore boson limit of the Lieb-Liniger model is particularly studied for ist mapping to free fermions [15], [16].Finite couplings can then be studied perturbatively around free fermions [51]- [56].For free fermions, any correlation function can be expressed in terms of the Wigner function W (x, λ) through Wick's theorem.And this Wigner function straightforwardly satisfies a hydrodynamic-like equation holding at finite x, t [40], [57]- [60] Moreover analogous continuity equations can be obtained for the Lieb-Liniger model perturbatively at large coupling around the free fermion limit [61].This equation is reminiscent to the GHD equation for the space-time dependent root density in the hardcore boson limit [5], [6] Besides, the integration of W (x, λ) over x in an equilibrium state yields the root density ρ(λ) of this equilibrium state.But to go from W (x, λ; t) satisfying (1) to an equation on a local root density ρ x,t (λ) is more subtle and is where lies hydrodynamics.It would require again a coarse-graining, a separation of scale and an assumption of local relaxation [40], and would only hold in certain limits.For finite x, t, although (1) holds true, the expression for some observables at equilibrium in terms of ρ(λ) would not be given by the expression in terms of W (y, µ; t) obtained with Wick's theorem.This expression would have "unwanted" terms that vanish only in certain limits.
The paper is organized as follows.In Section 2 we introduce the model, the out-ofequilibrium setup and the wavelet basis.In Section 3 we derive the representation of the state of the system out-of-equilibrium in this wavelet basis.In Section 4 we explain how to compute expectation values in this representation.In Section 5 we define the hydrodynamic limit that we consider, and derive a hydrodynamic behaviour in our setup.In Section 6 we present a numerical algorithm based on this wavelet for simulating the dynamics of hardcore bosons with inhomogeneous potentials.

Problem setting 2.1 Lieb-Liniger model and hardcore bosons
We consider the Lieb-Liniger model on a ring of size L with a Bose field φ(x) and periodic boundary conditions.We will be interested in the limit c → ∞ where the bosons become hardcore and map to free fermions [16].The normalized eigenfunctions in this limit c → ∞ are parametrized by a set of distinct quantized momenta and if N is odd For simplicity and without loss of generality, we will restrict to N even for the rest of the paper, and take K as defined in (4).The explicit expression of the normalized eigenstates is with when Their energy is

Wavelet representation
We divide the entire system irrespectively of the parity of the number of particles, and have identical representations as in (7).The particular choice of boundary conditions is purely for technical convenience, so that elements of K box can never be equal to elements of K in (4), and will not influence the physics.The important point is that they form a basis of the Hilbert space in box b, and that local observables can be expressed simply in terms of them.We note the factor n compared to (4), coming from the change of system size.As these states form a basis of the Hilbert space of wave functions on box b, a basis of the entire Hilbert space on [0, L] is obtained with tensor products of theses states in each box.Namely, given α (b)  α (b) α (b) ⊂ K box for b = 1, ..., n we introduce with by definition ᾱ ᾱ ᾱ ⊂ K ⊗ where We will denote ᾱ ∈ K ⊗ the couple ᾱ = ( 2πnm L , b), with α = 2πnm L ∈ K box the momentum and b(ᾱ) = b ∈ {1, ..., n} the box index.We will call |ᾱ ᾱ ᾱ a wavelet state.This term is motivated by the fact that ᾱ ᾱ ᾱ is a collection of "wavelets", i.e. small plane waves localized in each box.Explicitly, the wave function of the wavelet state |ᾱ ᾱ ᾱ is with N b denoting the number of particles in each : Schematic picture of the wavelet decomposition.The original system (top) is divided into boxes (bottom) and described by plane waves with different amplitudes and frequencies within each box.

Inhomogeneities and out-of-equilibrium protocol
Given the division into boxes defined above, we consider the following Hamiltonian where v b are potentials that are uniform in each box but that can differ from box to box.As it does not complicate the discussion, we will consider a general setting where the potentials are given an arbitrary time dependence v b (t).For an initial wave function |Ψ(0) , we propose to study the dynamics induced by this Hamiltonian.Namely, we define where the evolution operator is given by the time-ordered exponential We will consider a particular class of initial states |Ψ(0) in ( 14).We will assume that at time t = 0 the initial state is a wavelet state with ᾱ ᾱ ᾱ ⊂ K ⊗ fixed.

Change of basis
Given an observable O, its expectation value at time t can be expressed the following way in the wavelet basis The form factors of the time evolution operator ᾱ ᾱ ᾱ|U (t)| β β β are non-trivial, as a tensor product of eigenstates of H 0 in each box is not an eigenstate of H 0 in the entire system (and even less so with inhomogeneous potentials).To compute these form factors ᾱ ᾱ ᾱ|U (t)| β β β , we first compute the overlap between eigenstates of H 0 and wavelet states.
Proof.We have The last line imposes a matching x j = x κj for κ ∈ S N , and the two products of signs give a factor (−1) κ .Then we write N j=1 and set π = π κ and then π → π to obtain We now reparametrize , and π into πσ 1 ...σ n .We obtain We now perform the integrals, and using α This can be expressed as a determinant with ϕ given in the Lemma.

Form factors of the time evolution operator in the wavelet basis
We now prove the following result.
Proof.Using Trotter's expansion we have Let us fix T , introduce δt = t T and define Let us show that we can always write with some function ψ (m) t (ᾱ, β).For m = 0, this is true with Let us assume it is true for m − 1 and show it is true for m.We have We now use Lemma 1 to express the form factors ᾱ ᾱ ᾱ|λ λ λ as determinants.We then use the following Lemma to sum over λ λ λ, proven e.g. in [63], [64].
To obtain a well-defined differential equation in the Trotter limit δt → 0, we go to the interaction picture by introducing the Green's function for ᾱ, and defining ψ(m) We note that we have From (39) we obtain the recurrence relation on ψ(m) (43) Performing the sum over κ using we obtain This form allows for a well-defined Trotter limit δt → 0. We thus obtain the representation with ψ t (ᾱ, β) given by with ψt (ᾱ, β) satisfying the differential equation This concludes our proof.
The reader could notice that in Theorem 1 we could have written a differential equation directly in terms of ψ t rather than ψt .In fact, the matrix appearing in the differential equation for ψ t would be defined in terms of a divergent series.This reflects the fact that ψ t itself is not differentiable at t = 0.This can also be seen at the level of G t , whose expression in terms of a series (40) does not allow for term-by-term time differentiation.

Expectation values of local operators
where the trace is over all the boxes c = b.As we will see, another related quantity has a simpler expression in terms of ψ t (ᾱ, β).This is the cumulative reduced density matrix that we define by where by b(x x x) = b we mean that for all x ∈ x x x we have b(x) = b.We note that from this expression, the definition of the reduced density matrix ρ b would be obtained by imposing that β β β \ x x x and γ γ γ \ ȳ ȳ ȳ have no particles in box b.The terms for which β β β \ x x x and γ γ γ \ ȳ ȳ ȳ have exactly one particle in box b can be expressed in terms of the reduced density matrix as T [ρ b ] with the linear operator applying on density matrices with φ † (λ) the bosonic creation operator for the mode λ ∈ K box in box b.Similarly, terms with higher numbers of particles can be expressed through powers of T .One thus has The inverse of this relation is then seen to be where |α b denotes a single particle eigenstate with momentum α in box b.

Expression of
This cumulative reduced density matrix admits the following simple expression in terms of ψ t (ᾱ, β).
Theorem 2. For x x x, ȳ ȳ ȳ ⊂ K ⊗ localized in box b with M particles we have where we recall that |ᾱ ᾱ ᾱ denotes the initial wavelet state.If x x x, ȳ ȳ ȳ have a different number of particles, then x x x|ρ cumul b |ȳ ȳ ȳ = 0.
Proof.Firstly, we see in the definition of ρ cumul b in (50) that β β β and γ γ γ must have the same number of particles N as in |Ψ(t) .Since β β β \ x x x = γ γ γ \ ȳ ȳ ȳ, there has to be the same number of particles in x x x and ȳ ȳ ȳ.
From the definition we have Using Lemma 1 to express the form factors of U (t) as determinants we have with β β β = β β β ∪ x x x and β β β = β β β ∪ ȳ ȳ ȳ imposed to have N particles.In this expression we can include in the sum the terms where β β β contains elements of x x x or ȳ ȳ ȳ, since the determinants vanish in this case.Hence, setting x x x and ȳ ȳ ȳ to be the first M indices of β β β , β β β and writing β β β = { βM+1 , ..., βN } we have Then we write that the product of determinants is the determinant of the product det i,j and expand the determinant as a sum over permutations Expanding the product, we see that if we pick twice the same q for i and j, then the product of the two terms is invariant under swapping σ(i) and σ(j), whereas this changes the sign of (−1) σ , making the contribution vanish.Hence we have The sum over βM+1 , ..., βN factorizes into products of single sums that read Then the sum over τ does not depend on the ordering of τ −1 (M + 1), ..., τ −1 (N ).It can be converted into a sum over subsets {q 1 , ..., q M } ⊂ {1, ..., N } with a factor (N − M )! and with a remaining permutation over τ −1 (1), ..., τ −1 (M ).It yields We used that σ leaves invariant the set {q 1 , ..., q M } because of (66), to write q σ(i) instead of σ(q i ) after a reparametrisation of σ.The sum over σ is exactly a M × M determinant.Then it vanishes if two q i 's are equal, so the sum over the subset {q 1 , ..., q M } can be converted into M sums over q ∈ K with a factor 1/M !, and the sum over τ gives a factor M !. Hence which concludes the proof.

Generalities
The previous construction is straightforwardly generalized to observables in several boxes.Given two boxes b 1 , b 2 , we define the reduced density matrix and the cumulative reduced density matrix Then, an identical formula to Theorem 2 gives the expression for x x x|ρ b 1 b 2 |ȳ ȳ ȳ in terms of ψ t (ᾱ, β).

Single-particle Green's function
The wavelet representation allows for an efficient computation of certain observables which would be difficult to compute otherwise.Let us for example consider the single-particle bosonic Green's function Within an eigenstate of the entire system λ λ λ, this Green's function does not have a simple expression, because of the products of signs in (7) [30]- [35].Similarly, if x, y are two coordinates within the same box, no simple expression exists.However, if the two positions correspond to the beginning of two different boxes, i.e. if we set x = (b 1 − 1)n/L and y = (b 2 − 1)n/L, then we have from (12) for two wavelet states ᾱ ᾱ ᾱ, β β β means the set of all particles in ᾱ ᾱ ᾱ that belong to box b.
This yields the following expression in terms of the cumulative reduced density matrix in boxes In terms of ψ t (ᾱ, β), this is thus where we recall that the ᾱq 's are the particles in the initial state.

Definition
In order to define the hydrodynamic limit that we study in this paper, let us introduce typical scales for the parameters of the problem.We define τ the typical time scale in the problem and T the rescaled time by We define 2πκ the momentum spacing in each of the boxes, namely We define Λ the typical momentum of the particles in the wavelet initial state |ᾱ ᾱ ᾱ .We note that we necessarily have Λ κ.We finally define X the rescaled position by The typical scales τ and Λ are parameters of the problem that we can choose freely.They correspond to the time scale at which we study the problem and to the particle content of the initial state.The status of κ is more subtle.Since it directly depends on the number n of boxes that we divide the system into, it is not per se a physical parameter of the system, but rather parametrizes the precision at which we want to describe the system.However, through the fact that we impose the initial condition to be a wavelet state, it also constrains the minimal typical distance at which the initial state varies.
We define the hydrodynamic limit by the two conditions The first condition means that we study the system at a time scale at which the energy levels within each box appear continuous.The second condition means that we study the system at "large" times, the meaning of "large" being set by the initial state.These two conditions do not impose any behaviour of τ κΛ.As we will see, this quantity is the typical number of boxes that a particle moves through during a time interval of order τ .If τ κΛ → 0, the system will have no time dependence for the rescaled time T .If τ κΛ = O(1), the system is not divided into small enough boxes and will appear discontinuous at scales X, T .We will thus assume which will ensure that the system is smooth in terms of the rescaled space and time variable X, T .
The condition (79) can be met in quite different physical contexts.Let us mention in particular two situations.The first context (i) is that of the Euler scale, corresponding to large x, t at fixed x/t, which is the one mostly mentioned in the GHD literature.In our notations this means τ → ∞ and Λ = O(1).In this case the size of the boxes 1/κ becomes infinitely large.The second context (ii) is a short-time, high density limit.Because of the Pauli principle, an initial state with high density of particles will also have a large typical particle momentum Λ.This context translates thus into τ → 0 and Λ → ∞ at e.g.fixed τ Λ.The smoothness condition (80) imposes then κ → ∞, which is infinitely small boxes.No rescaling of space is needed in this context.It is similar to a semi-classical expansion → 0.

Generalized Hydrodynamics
Before taking the hydrodynamic limit of our model, let us very briefly summarize GHD, the hydrodynamic theory describing the Lieb-Liniger model [5], [6].As any hydrodynamic theory, it postulates that as long as local operators are concerned, the state of the system at rescaled coordinates X, T can be considered in an equilibrium state.In the hardcore boson gas, equilibrium states |λ λ λ are characterized by a root density ρ(λ), defined by the fact that there are Lρ(λ)dλ particles with momentum between λ and λ + dλ in the thermodynamic limit.GHD provides a partial differential equation for the root density ρ X,T (λ) describing the equilibrium state at rescaled coordinates X, T .In case of hardcore bosons and in absence of external potentials, it reads The solution can be readily expressed in terms of ρ X,0 (λ) the initial root densities In presence of an external potential V (X), the GHD equations are [19], [65] (83)

Hydrodynamic limit of the Green's function
In order to determine the hydrodynamic limit of the expectation values determined in Section 3, let us first determine the hydrodynamic limit of the Green's function G t (ᾱ, β).
Lemma 3.For α, β ∈ K box of order Λ, we have the following hydrodynamic limit We used the notation x for the integer part of x, and the discrete derivative defined by with the δ function modulo n δj = δ j mod n (−1) j/n .(88) Proof.We write with We have Let us first treat the case α = β.Then Using with x the integer part of x, we decompose it as with In the sum, in the hydrodynamic limit the values √ tλ are spaced by an order √ τ κ which goes to 0 from condition (i) of ( 79).Hence the sum can be seen as 1/κ times a Riemann sum for tuα is square integrable but not its derivative, the Riemann sum converges as O(κ √ τ ).This yields Moreover, since u → e −iu 2 −1 u is square integrable but not its derivative, its Fourier transform as a function of ω decays as 1/ω.Hence because √ t scales as √ τ and α as Λ in the hydrodynamic limit, we have Which of these two terms dominates is exactly given by the behaviour of τ κΛ, which is not constrained by the physical hydrodynamic behaviour of the system but by the precision with which we describe it.With our smoothness convention (80), the term O( √ τ ) dominates.In either case, R t (α) is subleading in (94) from both conditions of (79).Summing over q in (93 with the notation (87).Let us now consider the case α = β.In this case we have for Using we have Performing the sum over q yields when α = β These are the expressions given in the Lemma.

Hydrodynamic behaviour in absence of chemical potentials
We now would like to study the behaviour of the system out of equilibrium in the hydrodynamic limit.In absence of chemical potentials, we have ψ t (ᾱ, β) = G t ( ᾱ, β).

Particle content
Let us first probe the particle content of box b at time t.To that end, we consider the observable O b (w) introduced in ( 56) for a function w(β).Its expectation value in an equilibrium state in box b with root density ρ is given by w(β)ρ(β)dβ, allowing one to probe the value of the root density ρ in each box.With the expression (57), Theorem 2 and the fact that ψ t (ᾱ, β) = G t ( ᾱ, β), we have at all times t where we recall that ᾱ ᾱ ᾱ = {ᾱ 1 , ..., ᾱN } is the initial wavelet state.In the hydrodynamic limit, this suggests to introduce for a fixed function w(β) for ᾱ ∈ K ⊗ , with G hydro t defined in Lemma 3. Let us evaluate the hydrodynamic limit of this quantity.Importantly, we will assume that w(β) varies only at scale much larger than κ the typical spacing between energy levels in a box.This means we will assume We note that this means the Fourier transform of w(β) varies on a scale that is much smaller than the size of a box 1/κ.
Using the expression for G hydro t , we have Writing α and β in terms of κ times integers, we have that only finite integer difference (but that can be arbitrarily large) will contribute to the sum in the hydrodynamic limit, because of the factor But then we have 2tκα = 2tκβ , from assumption (ii) of (79).Since w /w is of order o(κ −1 ), approximating w(β) by w(α) comes with an error o(1) which goes to 0 in the hydrodynamic limit.Hence Using (100) in the limit θ → 0 we have This yields It can be rewritten as This contribution of a single initial particle ᾱ to the root density is plotted in the left panel of  117), right) as a function of 2tκα.The left panel represents the contribution of initial particle ᾱ to the root density in the hydrodynamic limit, and the right panel to an "anomalous" correlation function that vanishes at equilibrium.In presence of several initial particles, the function on the right will undergo destructive interference.

Local relaxation: general remarks
We showed in Section 5.4.1 that the particle content of box b in the hydrodynamic limit exactly follows the GHD equations.The effective root density ρ b,t in box b at time t is then given by where β ∈ K ⊗ with b( β) = b.We say it is "effective" since it is only defined by the fact that the expectation value of O b (w) at time t is β ρ b (β)w(β), as in an equilibrium state.But for it to really be the root density describing local observables in box b, one needs as well local relaxation.Let us take the example of the connected two-point function of the density operator defined in (54).Its expectation value is (113) This can be written as with w(β) = e iβ(y−x) and where we introduced where β + 2πκδ ∈ K ⊗ denotes the wavelet in same box as β but with momentum β + 2πκδ.In an equilibrium state |β β β , this expectation value would be Hence we see that this formula does not correspond to (114) with the effective root density (112).The effective root density would only capture the terms with δ = δ = 0. To show local relaxation to an equilibrium state in box b in the hydrodynamic limit, one thus needs to show that the terms for which δ = 0 or δ = 0 vanish in the hydrodynamic limit.

Local relaxation: simplest case
Let us define the hydrodynamic limit of P w (ᾱ; δ) for a fixed function w(β) that varies on a scale larger than 1/κ.We note that for the example of the density 2-point function of Section 5.4.2where we have w(β) = e iβ(y−x) , this condition means that the distance separation y − x is much smaller than the size of the box 1/κ.We write and use the expression for G hydro t ( ᾱ, β) in Lemma 3. Let us first focus on the sum, that we denote S. We have The terms with α − β much larger than κ go to zero, whereas the terms with α − β = O(κ) are of order O(1).Hence only the terms for which α − β = O(κ) contribute.But then, we have 2tκα = 2tκβ = 2tκ(β + 2πδκ) because tκ 2 → 0 from (80), and w(β) = w(α) + o(1).Hence Writing β and α in terms of integers, and using κ 2 t → 0 from (79), this yields Using (93) in the limit θ → 0 we get This yields Coming back to (118) and using (86) we obtain (124) This contribution of a single initial particle ᾱ is plotted in the right panel of Fig 2 .Contrary to P hydro w (ᾱ; 0), these contributions P hydro w (ᾱ; δ) vanish when integrated over ᾱ.Hence when several particles are included in the initial state with close momenta, these contributions vanish.

Local relaxation: summary
Let us summarize our findings.In Section 5.4.1 we showed that in the hydrodynamic limit defined in (79), the effective root density (and so the expectation value of any conserved charges) follows from the GHD equations particle by particle.This means that as long as one is interested in conserved charges, no large number of particles is required for the GHD description to be valid.Now, the usual meaning of hydrodynamics is that a state can be considered to be in an equilibrium state as far as expectation values of local observables are concerned.This means more than just having the right expectation values of conserved charges, as explained in Section 5.4.2, since these do not fix all correlations functions in an out-of-equilibrium state (contrary to equilibrium).In Section 5.4.3, we showed that this local relaxation is not achieved for single particles, since there are some terms that do not vanish in the hydrodynamics limit (79) whereas they would vanish in an equilibrium state.However, these terms have zero integral and so will undergo destructive interference when several particles are included in the initial state.Hence, to summarize, local relaxation requires a large number of particles to occur, but expectation values of conserved charges do not require a large number of particles to be described by GHD in the limit (79).
6 Algorithm for simulation of hardcore bosons 6.1 The algorithm Theorem 1 allows for an efficient algorithm to simulate hardcorse bosons in 1D.According to this result, the function ψ t (ᾱ, β) fully characterizes the state during the out-of-equilibrium protocol described in Section 2. Expectation values of local observables can be expressed directly in terms of this function, according to Theorem 2. This function ψ t (ᾱ, β) represents the state in a way that is mixed between real space (through the dependence in the box indices b(ᾱ), b( β) ∈ {1, ..., n} of ᾱ, β ∈ K ⊗ ) and momentum space (through the dependence in the momenta α, β ∈ K box ).This mixed representation is key to the efficiency of the algorithm.
Its resolution in space is well adapted to inhomogeneous settings and the computation of local observables, while its resolution in momentum is well adapted to compute the time evolution of the state efficiently and with good precision.
In practice, we consider an even number of boxes n and define K and K (P ) a truncated version of K and see ψ t (ᾱ, β) at fixed t as a nP × nP matrix.For a given time step δt, we compute its time evolution iteratively through with the unitary matrix Here, we defined ϕ (P ) (λ, ᾱ) as the nP × nP unitary matrix obtained by applying a QR decomposition to the nP ×nP matrix (ϕ(λ, ᾱ)) λ∈K (P ) , ᾱ∈K (P )

⊗
. The iteration (127) can be implemented as a simple matrix multiplication.The exact time evolution is recovered in the limits δt → 0 and P → ∞.
Some comments on the QR decomposition are in order.Firstly (i), the unitary matrix ϕ (P ) (λ, ᾱ) obtained from a QR decomposition is unique only up to a right multiplication by a diagonal matrix with modulus 1 coefficients.But this corresponds to conjugating ψ t with this diagonal unitary matrix, which has no effects on the expectation values of observables.Secondly (ii), the QR decomposition while preserving the unitarity of U breaks parity symmetry, namely the resulting U is not invariant under simultaneous reversing of the box indices b → n + 1 − b and momenta α, β → −α, −β, whereas it would be so without the QR decomposition.This will introduce a slight left/right asymmetry in the results even if the initial condition is parity symmetric.This symmetry is however restored in the limit P → ∞.

Protocol
As an example of application of the algorithm, we consider a quantum Newton cradle setup [23], [36].At time t < 0, the system of N particles is initialized in the ground state of the Tonks-Girardeau model with a harmonic potential Then at time t = 0 a "Bragg pulse" is applied to the system, that is a potential given by with a large space frequency q and large amplitude A, for a short time window ∆t pulse .Then the system is time evolved with the original potential (129).

Numerical implementation
In our setting, the initial state is necessarily a tensor product of energy eigenstates in each of the boxes, which cannot encode the ground state of the Tonks-Girardeau model with the potential (129).We thus have to include a first initial stage to prepare this ground state, before applying the Bragg pulse.We proceed as follows.We initialize the system with one particle in the zero momentum space of each of the n/4 boxes from box 3n/8 to 5n/8.Numerically, we take n = 256 boxes (so N = 64 particles) and set the system size to L = 32.Then from t = 0 to t = t 0 we slowly change the potential from an infinite (or large) square well to the potential (129).At time t = 0, the wave function relaxes quickly to the ground state of the square well potential, which then according to the adiabatic theorem remains in the ground state of the time-dependent potential, provided the potential changes slowly enough.Specifically, we apply the following potential for 0 < t < t 0 2 where b ∈ {1, ..., n} is the box index, and with the numerical values t 0 = 2 and ω 2 = 16000.Then, at time t 0 and onward, we fix the potential to be (129).
To modelize the Bragg pulse, we consider the instantaneous limit with ∆t pulse → 0 and A → ∞ at fixed A∆t pulse , as in [36].Then the effect of the pulse is at t = t 0 to apply on the wave function the operator e −i∆t pulse dxV pulse (x)φ † (x)φ(x)dx .Since our setup requires the potential to be constant in each box, we approximate V pulse to be constant equal to A on even boxes, and to −A on odd boxes.Numerically, we take A∆t pulse = π/2.
Finally, for comparison purposes, we carry out the same simulation with a quartic potential V quartic (x) ∝ x 4 , with a proportionality constant chosen such that V quartic (1/6) = V (1/6).
We present the numerical results in Fig 3 and 4. In Fig 3 we measure the number of particles per unit length at the middle of the system, averaging over 8 boxes, and plot it as a function of time.We observe a large number of oscillations with a large amplitude for the harmonic potential.With the quartic potential, we see that the oscillations are immediately damped.
In Fig 4 we present the entire profile of expectation value of the particle number as a function of the box index b, at different times.We display as well the mode occupation number for the P box momenta β ∈ { 2πnp L , p = −(P − 1)/2, ..., (P − 1)/2}, averaged over all the n boxes.⊗ , averaged over all the boxes, displayed in increasing order from left to right.For better readability, the vertical scale of the green bars is arbitrary but consistent for different times.The algorithm parameters are δt = 0.0001 and P = 15.The first plot at t = 2 is just before the Bragg pulse, all the others correspond to times after the Bragg pulse.

Convergence towards hydrodynamics
We now study the convergence of the numerics with the hydrodynamic limit (79) in the quantum Newton cradle setup described in Section 6.2.2.We implement the hydrodynamic limit (79) the following way.If the system is divided into n boxes, we define the rescaled time with the rescaled Bragg pulse instant τ 0 = 512 n , and the rescaled potential amplitude The numbers are chosen so as to recover those of Section 6.2.2 for the value n = 256.The initial state is always chosen to be one particle in the zero momentum mode of the n/4 boxes at the middle of the system, as in Section 6.2.2 (there are thus N = n/4 particles).The number of boxes n plays the role here of the parameter to scale to reach the hydrodynamic limit (79), which is obtained in the limit n → ∞.We note that this exactly corresponds to a short-time, high-density limit mentioned below (79).We measure the expectation value of the particle number at the middle of the system N n/2 .We show in  6.4 Quench from double well to single well potential: convergence with δt and P We now present numerical results of the algorithm in a setup analogous to the quantum Newton cradle setup, which is a quench from a double well potential to a single well potential.At time t = 0, we take the initial wavelet state to be a tensor product of lowest energy eigenstates with m particles in each of the leftmost n/4 boxes and rightmost n/4 boxes (the total state has thus nm/2 particles).We choose this particular initial setup only to ease the first relaxation stage with the double well potential.Then for time 0 < t < t 0 we impose the following potential with x = b−1/2 n − 1 2 where b ∈ {1, ..., n} is the box index and with ω 2 an overall scale.Then for time t > t 0 we time evolve the system with the single potential well

Summary and discussion
In this paper we studied the emergence of hydrodynamics in a 1D gas of hardcore bosons in inhomogeneous potentials from first principles.We provided an exact implementation of the ubiquitous "fluid cell" view of hydrodynamics.Specifically, we divided the system into n different boxes and decomposed the out-of-equilibrium wave function of the system in the basis given by tensor products of energy eigenstates on each of the boxes.Such representation is mixed between real space and momentum space, and is composed of small plane waves localized in each box, called "wavelets".Contrary to a hydrodynamic or coarse-grained description, there is no removal of degrees of freedom in this wavelet representation.Because hardcore bosons can be mapped to free fermions, the dynamics in this wavelet basis can be exactly computed using form factor expansions.The output of this work is twofold.Firstly, we provide a derivation of (generalized) hydrodynamics for hardcore bosons from first principles, without the assumption of local equilibrium which was systematically required previously.In particular we emphasize that hydrodynamics emerges in a short-time, high-momentum limit.In this limit, the expectation values of conserved charges do not require a large number of particles in the initial state to be described by GHD equations.However, local relaxation (that is required to deduce all local correlations from the conserved charges) is obtained only when including a large number of particles to produce a destructive interference of "non-relaxed terms".Regarding the cold-atom experiments applications of GHD, this seems more realistic than the Euler limit of large time and large space.Secondly, we show that the decomposition of the system in this wavelet basis provides an efficient numerical algorithm to simulate the dynamics of hardcore bosons with inhomogeneous potentials.It is indeed a representation that is mixed in real space (through the division into boxes) and momentum space (through the decomposition into energy eigenstates in each box) that can be implemented efficiently and that gives access to several observables for a large number of particles and for a relatively long time, including the single-particle Green's function, otherwise difficult to compute.
The results presented in this paper suggest a number of directions.A natural direction would be to generalize this approach to finite coupling c between the bosons.The strategy would be to perform a strong coupling 1/c expansion, which has proven efficient before [54], [56], [61].At large finite coupling the bosons can be mapped to fermions with a weak coupling, allowing for perturbative expansions [51], [52], [56], [66].Another approach could be to generalize the results of geometric quenches to work directly at finite coupling c [67].Other possible directions include using the algorithm we introduced to simulate more closely cold-atom experiments, for example by including atom losses which play a significant role [68]- [72].
[0, L] into n boxes [(b−1) L n , b L n ] of equal size L/n, for b = 1, ..., n, as depicted in Fig 1.We can define in each box b the same Hamiltonian H 0 restricted to this box, and set anti-periodic boundary conditions for even number of particles, and periodic boundary conditions for odd number of particles.The eigenstates of this Hamiltonian are parametrized by sets of distinct momenta α (b)

4. 1
Cumulative reduced density matrixWe now would like to express expectation values of observables localized in box b within state |Ψ(t) in terms of ψ t (ᾱ, β).To that end, we need to compute the reduced density matrix in box b ρ b = tr c=1,...,n c =b

)
Hence knowing ρ cumul b or ρ b is equivalent.

4. 2
Expectation values within a box in terms of ρ cumul b In fact, the expectation value of many local operators can be directly expressed in terms of ρ cumul b , without using ρ b .For example, let us take the example of the two-point function of the boson density operator

with |p 1
, p 2 b denoting the state with only two particles p 1 , p 2 in box b.Let us define observables that allow us to probe the particle content of a box b.Given a function w(β), we define the operator O b (w) as being localized in box b and being diagonal in the wavelet basis |ᾱ ᾱ ᾱ with eigenvalue O b (w)|ᾱ ᾱ ᾱ = ᾱ∈ᾱ ᾱ ᾱ b( ᾱ)=b w(ᾱ)|ᾱ ᾱ ᾱ .(56) The expectation value of O b (w) in state |Ψ(t) is readily expressed in terms of ρ cumul b O b (w) t = α∈K box w(α) b α|ρ cumul b |α b ,

where b(x x x) ∈ {b 1 , b 2 }
means that for all x ∈ x x x, we have b(x) ∈ {b 1 , b 2 }.A formula identical to (52) relates ρ cumul b 1 b 2 to ρ b 1 b 2 .Expectation values of operators localized in boxes b 1 , b 2 can also be written in terms of ρ cumul b 1 b 2 .For example, with O b (w) denoting the operator in (56), we have with x ∈ [0, L] the original position, related to the box index b by b = xκ .

Figure 3 :
Figure 3: Quantum Newton cradle setup.Number of particles per unit length at the middle of the system as a function of time, for the setup described in Section 6.2.2, with the quadratic potential (orange, top panel) and the quartic potential (green, bottom panel).The insets show the profile of the expectation value of the particle number in each box, as a function of the box index, at times indicated by the black line.The algorithm parameters are δt = 0.0002 and P = 15.The dotted vertical lines indicate the times at which the Bragg pulse is applied.

Figure 4 :
Figure 4: Quantum Newton cradle setup.Setup described in Section 6.2.2 with the harmonic potential, at different times t.Orange curve: expectation value of the particle number N b in each box, as a function of the box index b.Green bars: amplitude of each of the P box mode occupation number α ∈ K (P ) Fig 5 the results of the numerics.

Figure 5 :
Figure 5: Convergence towards the hydrodynamic limit.Expectation value of the particle number at the middle of the system N n/2 as a function of the rescaled time τ , for the setup described in Section 6.3, with N = 8 particles (pink), N = 32 particles (orange), N = 64 particles (green), N = 128 particles (purple).The algorithm parameters are δt = 0.0256/n and P = 7.

V ( 2 )
(x) = ω 2 x 2 .(135)We present numerical tests of the convergence speed of the algorithm with δt and P .To fix values, with set L = 32 the size of the system, n = 256 the number of boxes, m = 1 the number of particles in each occupied box initially (so N = 128 particles in total), V = 8000 the strength of the potential, and t 0 = 2 the duration of the first relaxation stage, for a total running time t fin = 16.Then we measure N b the number of particles in box b as a function of time t, for different algorithm parameters δt and P .The results are plotted in Fig6 and 7. We see that for these values of physical parameters, the algorithm parameters δt = 0.0005 and P = 15 show almost converged curves.The simulation with these parameters takes just a few minutes on a laptop.

Figure 6 :
Figure 6: Influence of the time step δt.Number of particles per unit length at the middle of the system as a function of time, for the physical parameters quoted in Section 6.4, with δt = 0.005 (orange), δt = 0.002 (green), δt = 0.0005 (purple), all for P = 15.The dotted vertical lines indicate the moment of the quench from the double well to the single well.

Figure 7 :
Figure 7: Influence of truncation parameter P .Expectation value of the number of particles (left panel) and of the momentum (right panel) in each box, as a function of the box index b, for the physical parameters quoted in Section 6.4.These are snapshots at times t = 2.75 (continuous line) and t = 3.15 (dotted) for the left panel, and at time t = 3 for the right panel.The algorithm parameters are δt = 0.0005 for all curves, with P = 3 (orange), P = 7 (green), P = 15 (purple).The momentum curve is averaged over a time window of size 0.03 to remove very fast oscillations.