BBGKY Hierarchy and Generalised Hydrodynamics

We consider fermions defined on a continuous one-dimensional interval and subject to weak repulsive two-body interactions. We show that it is possible to perturbatively construct an extensive number of mutually compatible conserved charges for any interaction potential. However, the contributions to the densities of these charges at second order and higher are generally non-local and become spatially localized only if the potential fulfils certain compatibility conditions. We prove that the only solutions to the first of these conditions are the Cheon-Shigehara potential (fermionic dual to the Lieb-Liniger model) and the Calogero-Sutherland potentials. We use our construction to show how Generalized Hydrodynamics (GHD) emerges from the Bogoliubov--Born--Green--Kirkwood--Yvon hierarchy, and argue that GHD in the weak interaction regime is robust under non-integrable perturbations.

We consider fermions defined on a continuous one-dimensional interval and subject to weak repulsive two-body interactions. We show that it is possible to perturbatively construct an extensive number of mutually compatible conserved charges for any interaction potential. However, the contributions to the densities of these charges at second order and higher are generally non-local and become spatially localized only if the potential fulfils certain compatibility conditions. We prove that the only solutions to the first of these conditions are the Cheon-Shigehara potential (fermionic dual to the Lieb-Liniger model) and the Calogero-Sutherland potentials. We use our construction to show how Generalized Hydrodynamics (GHD) emerges from the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy, and argue that GHD in the weak interaction regime is robust under non-integrable perturbations.
Finding an efficient description for the non-equilibrium dynamics of quantum many-particle systems has been a key problem in theoretical physics since the birth of quantum mechanics [1]. During the last two decades there has been an upsurge in interest as a result of significant experimental advances [2] and potential technological applications [3]. In spite of remarkable progress in our understanding of out-of-equilibrium quantum matter [4][5][6][7][8][9][10], however, an efficient and accurate description of its dynamics remains out of reach.
In essence, what makes this problem so hard is the lack of general methods to tackle it. Exact methods are restricted to a small number of fine-tuned many-body systems [11][12][13][14][15] and there is currently no controlled approximation scheme applicable to generic interacting systems. Numerical methods are typically limited by the accessible system sizes [5,16] or by the class of initial states that can be accommodated [17]. Even in one dimension, where techniques based on matrix product states [18][19][20] give access to large systems, their applicability is limited to short times by the rapid growth of quantum entanglement. Furthermore, continuum models -which describe many relevant experiments -provide additional obstacles to numerical approaches.
For weakly interacting systems the situation is substantially simpler because, at least in principle, a general description of the dynamics can be attained by using the celebrated Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [21][22][23], which encodes the Heisenberg equations of motion for the reduced density matrices. To illustrate it let us consider for simplicity and concreteness a system of spinless fermions: in this case the reduced density matrices are written as ρ n (x, y; t) = Tr ρ(t)ψ † x1 . . . ψ † xn ψ yn . . . ψ y1 , (1) where ψ † x and ψ x are fermionic creation and annihilation operators and ρ(t) denotes the density matrix specifying the state of the system at time t. Assuming that the fermions interact via a two-body potential V (x − y) the BBGKY hierarchy takes the form i∂ t ρ n (x, y; t) − (H (n) x − H (n) y )ρ n (x, y; t) = (2) where H (n) x is the first-quantized n-particle Hamiltonian in the position representation. In principle the system (2) gives a complete account of the non-equilibrium dynamics of many-particle quantum systems. In practice it is necessary to truncate it, which can be justified e.g. for weak interactions. Retaining only two-particle cumulants gives rise [24] to the ubiquitous quantum Boltzmann equation (QBE) for the Wigner function (3) The QBE holds in the regime of weak spatial variations and large times, i.e. the Euler scaling limit [24], and under the assumption of local relaxation can be further reduced to a set of hydrodynamic equations, which are obtained from the continuity equations of particle number (or mass), energy, and momentum [22, 23]. The situation is very different, and significantly richer, in quantum integrable systems [25,26], which are characterized by having an extensive number of (mutually compatible) conservation laws with densities that are sufficiently local in space. These conservation laws give rise to the existence of stable quasi-particles over macro states at finite energy densities [27]. Integrable models do not thermalize but relax to a much wider class of equilibrium states known as Generalized Gibbs Ensembles [28,29], which can be fully characterized in terms of their respective quasi-particle density in momentum space ρ(k). In these systems the dynamics of local observables close to local equilibrium is described by GHD [30,31], which can be expressed as an evolution equation for a space and time dependent density ρ x,t (k) of the stable quasiparticles Here v x,t (k) is a (known) quasi-particle velocity that depends on ρ x,t (·). GHD is obtained from the system of continuity equations for the extensive number of conservation laws, postulating local relaxation to an equilibrium state with quasi-particle density ρ x,t (k), and then inferring its evolution equation (4). The GHD equation (4) can be viewed as a kinetic theory governed by a dissipationless Boltzmann Equation [30] where the velocity is a non-linear functional of the density itself [32,33].
Hence it suggests the existence of an operator n(x, k) expressed in terms of fermions ψ x localized near x, whose expectation value is ρ x,t (k), and for which the BBGKY hierarchy would re-organize into a dissipationless QBE. Such an operator, however, is only known in the noninteracting case, where it is given by the Wigner operator (whose expectation value is (3)).
In this letter we show that under certain conditions the BBGKY hierarchy can be reformulated as GHD equations. We establish that, surprisingly, one can construct conserved charges for arbitrary local but weak interaction potentials by appropriately dressing the non-interacting modes. Crucially, however, these conserved charges are generally non-local and acquire good locality properties only for integrable models. This means that even though it is possible to explicitly construct operators fulfilling a dissipationless QBE for generic interacting models, these operators are non-local and the equation cannot be used to infer the dynamics of physical observables. In contrast, in the integrable case the evolution equation reduces to the GHD equation for expectation values of local mode occupation numbers. As discussed in the following, this also gives a simple operational criterion to assess the integrability of a given interaction potential.
In the remainder of this letter we consider a system of interacting fermions on a ring of size L with Hamiltonian where all summations are on "free momenta" p = 2πn/L with n ∈ Z and where ψ † p and ψ p are fermionic creation and annihilation operators with canonical anti Here V (p) is a two-particle interaction potential, which by virtue of the anti-commutation relations can be cast in the form where A k1···kn is an operator that acts by antisymmetrising in {k 1 , . . . , k n }. Without loss of generality we take V (k) to be even in k and V (0) = 0 [34].
In order to obtain GHD we follow the logic of [30, 31], but aim to explicitly construct the conserved charges from the BBGKY hierarchy. These are derived by assuming that the potential V (k) can be treated perturbatively. We note that for the potential (7) perturbation theory is believed to have a finite radius of convergence, cf. [38][39][40].
Our first step is to find all the conserved charges of (5). We begin by considering an operator Q and assuming that it has a regular perturbative expansion where where we set Q −1 = 0. In addition, we restrict our discussion to cases where the particle number operator N = p ψ † p ψ p is part of the tower of conserved charges. This implies that all Q n are expressed as sums of monomials involving equal numbers of ψ † p and ψ q . It is easy to see that the zeroth order of (9) is solved by the following number-conserving charges where f (p) is any function, which we take to be smooth. We note that charges corresponding to different f s commute, but (10) is clearly not the most general choice for a conserved charge of H 0 . Any energy-diagonal operator (i.e. any operator that is diagonal in an eigenbasis of H 0 ) could be taken as a conserved charge of H 0 , with any V -dependent prefactor appearing at higher orders. The densities of such more complicated energy-diagonal operators, however, are generically non-local in real space. We therefore restrict our search to zeroth-order charges of the form (10) and leave out energy-diagonal contributions at higher orders. We expect this "minimal ansatz" to be sufficient for the perturbative construction of complete sets of conserved charges in all integrable models of the form (5) featuring a single species of quasi-particles, i.e. cases in which the stable quasi-particles can be thought of as being "adiabatically connected" to the free fermions ψ † p ψ p .
Starting from (10) we can directly use (9) and our minimal ansatz to recursively generate the higher orders of Q f . In particular, for m = 1 the equation can be solved for any potential giving [41] where the function g (4) f ;1 (k) is non-singular for any choice of V (k), and has the same regularity as V (k): its explicit expression is given in (sm-14).
In contrast, for m = 2 Eq. (9) does not always admit a solution in the framework of our minimal ansatz, i.e. starting from the zeroth order (10) and omitting energydiagonal contributions at higher orders. This is because [H 1 , Q f ;1 ] generically contains an energy-diagonal component and therefore cannot be expressed in the form [H 0 , Q f ;2 ]. Indeed, the energy-diagonal part of the latter commutator is always zero. In particular, defining S 6 := M 6 \ N 6 with we find the following solvability condition for every k ∈ S 6 .
This condition can be shown to be equivalent to the request that the connected 2nd order contribution to the three-particle S-matrix vanishes for different sets of incoming and outgoing momenta [44]. If Condition 1 is fulfilled the second order charge can be expressed as [34] f ;2 (k)} n=4,6 are regular and their explicit expressions are given in Eqs. (sm-33) and (sm-34).
Proceeding in this way we obtain a set of charges conserved up to O(V 3 ). Our construction can be readily extended to higher orders. Crucially the charges (16) are mutually compatible (it is shown in [34] that they commute up to O(V 3 )) and for smooth f (k) and V (k) they are quasi-local, i.e. their density is exponentially localised. The latter property can be shown by expressing Q f ;m in terms of real space fermions ψ(x) = p ψ p e ixp / √ L. For instance, considering (11) we have where In the thermodynamic limit C f ;1 (x) becomes the Fourier transform of g (4) f ;1 (k), which is smooth for f (k) and V (k) smooth. Therefore, quasi-locality is guaranteed by standard Fourier analysis [42]. If V (k) is regular but not smooth, the densities still decay in real space, but generically as power-laws.
If Condition 1 is not fulfilled one can construct charges at second order by either: (i) adding appropriate energydiagonal contributions to Q f,0 (or Q f ;1 ) to subtract the energy diagonal part of [H 1 , Q f,1 ] [43]; (ii) appropriately deforming the dispersion relation of the free model: k 2 → k 2 + ǫη(k), with ǫ ≪ 1 and η(·) a suitable function. Using this deformed dispersion in (12) one can ensure that M 6 coincides with N 6 in any finite volume L and hence Condition 1 is always fulfilled for finite L. Both strategies (i) and (ii), however, necessarily produce second-order charges where the coefficient of the six-fermion term (g f ;2 (k) in Eq. (15)) has singularities in momentum space. In real space, these singularities translate into a charge density that does not decay, i.e. the second-order charges are non-local.
In summary, restricting to zeroth order charges with a quadratic term, Condition 1 is necessary and sufficient for the system (5) to have a complete set of charges with densities that decay in space at second order. We remark that if one changes the dispersion relation, for instance by considering a system on the lattice, the situation becomes richer and will be discussed in a separate work [44].
Let us now characterise the solutions to Condition 1. We have the following [34] Property 1. The only potentials fulfilling Condition 1 and admitting a power-series expansion around 0 are Restricting to b > 0 (we seek potentials that are well defined for all k ∈ R) we see that that (19) correspond to the inverse-sinh-squared Calogero-Sutherland potential up to mass and momentum rescaling [26]. In particular we find the two following limiting cases The first is nothing but the Fourier transform of the integrable Cheon-Shigehara potential (7) at order O(β) when the regulator is removed, while the second is the inverse-squared Calogero-Sutherland potential [26]. We have conducted a numerical check of Condition 1 for several classes of singular potentials but have failed to find any additional solutions. We note that to the best of our knowledge, (19) and (20) correspond to the only known integrable potentials for (5) in the thermodynamic limit. Importantly both cases give rise to theories with a single species of quasiparticles (smoothly connected to free fermions for vanishing V ). For the k 2 potential one can use a standard result of Fourier analysis [42, 45]to construct a complete set of "ultra-local" charges [25] with density supported on a single point. To this end it is sufficient to take the set of charges constructed choosing f (k) ∈ {k 2n } ∞ n=1 . Given the set of conserved charges (16) we can now define the quasi-particle number operator n(k, x). Our starting point are the operatorial densities q f (x) of the charges (16), which can be chosen with at least a powerlaw decaying density if Condition 1 is fulfilled. A convenient choice is corrected by higher orders in V . Note that this choice of density operator is not Hermitian but one can straightforwardly recover an Hermitian density via q f (x) → (q f (x) + q † f (x))/2. To construct n(k, x) we consider linear combinations of {q f (x)} with smooth f 's and select the contributions of a single non-interacting mode k. The linear combination is then chosen such that we obtain the density associated with f (p) = δ k,p . Since (21) is linear in f this is always possible if {f (k)} is a complete set in L 2 (R). Specifically we propose the following definition where k obeys free quantisation conditions. This operator is characterised by the following properties: As shown below this implies that the thermodynamic limit of the expectation value of n(k, x) on a translationally invariant state gives the density of quasi-particles.
(ii) n(k, x) is conserved, i.e., it fulfils where j f (x) is the current associated with q f (x) via the continuity equation [34] (iii) An explicit expression of n(k, x) can be obtained from (21). We note that for V = 0 the resulting expression differs from the free fermion Wigner operator, but the latter can be easily recovered by modifying the definition of the densities of the conserved charges.
If the moments of n(k, x) are sufficiently local in space [46] then (i) and (ii) allow us to interpret Eq. (24) as an operatorial progenitor of the GHD equation (4). Property (iii) makes the relation between GHD and the BBGKY hierarchy explicit. The expectation value of n(k, x), which fulfils the GHD equation (4), is recovered as a specific sum of cumulants fulfilling the BBGKY hierarchy. The expansion of n(k, x) up to order O(V n ) involves cumulants of order 2n + 2. Therefore, in order to recover the full GHD equation one needs the entire BBGKY hierarchy.
Our previous discussion implies that, at second order, n(k, x) is sufficiently local only if V (k) is taken to be the Cheon-Shigehara or Calogero-Sutherland potential. Remarkably, however, at first order the moments are quasi-local for any smooth potential V (k), integrable or not [47]. This suggests that in the weak-interaction regime GHD physics is robust against non-integrable perturbations. Namely, having a weak non-integrable potential in (5) does not preclude the formulation of the BBGKY hierarchy in terms of a GHD equation at first order in perturbation theory, but only at second order. This means that GHD will be applicable on a longer time scale than naively expected and perhaps partially explains the relevance of GHD in modelling actual experiments [48][49][50].
In order to fully make contact with GHD, the thermodynamic limit of the expectation values of n(k, x) and j n (k, x) in an energy eigenstate should match the known formulae for Thermodynamic Bethe ansatz (TBA) integrable models. Namely, at order m in perturbation theory one should have where ρ(k) is the quasi-particle density in momentum space and v(k) is the group velocity of stable particle and  (27) are fulfilled for any potential. Namely, using ψ † p ψ q ≡ δ p,q ϑ(q) and ψ † p ψ † q = 0, we find at first order [43] where K(k) = ∂ k (V (k)/k). These agree with the first order expansion of ρ(k) and v(k) in an integrable model with a single species of quasi-particles and a scattering phase shift V (k)/k (perturbatively small). The latter is then interpreted as the scattering phase shift of the effective integrable model describing (5) at first order.
In particular, using the potentials (19) and (20) we recover the first order expansion of the known TBA expressions of ρ(k) and v(k) in the Lieb-Liniger and Calogero-Sutherland models [34]. The extension of these results to higher orders in the integrable case will be reported elsewhere [44].
We note that at first order in V the above programme can be generalized to the bosonic case, i.e. when the operators ψ entering the Hamiltonian (5) satisfy canonical commutation relations. We are again able to construct charges for any potential, but now the expectation value of these charges in an eigenstate is in general divergent at first order in V , as reported in Sec.VI of the SM.
Discussion. In this Letter we showed how to systematically derive a dissipationless Boltzmann equation for certain "dressed" quasi-particles from the BBGKY hierarchy for weakly interacting many-particle systems and derived an explicit expression for the density of the corresponding mode occupation operator. In order to enable a GHD description this density must have good locality properties, which we find to be the case precisely for the known integrable potentials. This suggests that our "integrability condition" provides an exact characterisation of all integrable systems (even those not solvable by Bethe ansatz) with a single species of quasi-particles. Our construction further suggests that in weakly interacting models the GHD description is robust against non-integrable perturbations, which is relevant for applications of GHD to cold-atom experiments. Our work can be extended in a number of directions. Firstly, the case of integrable models with several species of stable quasi-particles can be analysed in a similar way, but a number of interesting complications occur [44]. Secondly, from our "operatorial" GHD equation it should be possible to derive corrections to GHD [55-59] in a systematic fashion. Finally, our approach can be used to study the effects of general weak integrability breaking interactions [60-67]. An important goal is to investigate how, and over what time-scales, they render perturbed GHD descriptions invalid.
Acknowledgments: This work has been supported by the Royal Society through the University Research Fellowship No. 201101 (BB) and the EPSRC under grant EP/S020527/1 (EG and FHLE).   (11) and (15); (ii) a proof of Property 1; (iii) a proof that the charges are in involution; (iv) an explicit expression of the currentĵ f,x (t) at first order; (v) the derivation of (28) and (29)   Supplemental Material for "BBGKY Hierarchy and Generalised Hydrodynamics" Here we report some useful information complementing the main text. In particular -In Section I we explicitly compute the charges at the first two orders in perturbation theory.
-In Section II we prove Property 1.
-In Section III we show that the charges (16) commute up to second order.
-In Section IV we present an explicit expression for the current j f,x associated to the density q f,x in Eq. (21).
-In Section V we present an explicit Thermodynamic Bethe ansatz derivation of (28) and (29) in the cases of Lieb-Liniger and Calogero-Sutherland.

I. PERTURBATIVE DETERMINATION OF THE CHARGES
A. First Order The first order correction to the charge, Q f ;1 , is obtained by solving the following equation where Q f ;0 is given in Eq. (10). The commutator on the r.h.s. is immediately evaluated as Therefore making the ansatz The latter equation admits solutions only if the r.h.s vanishes when the l.h.s. does. This means that we should have where we set First we note that if the condition (sm-5) is immediately satisfied. Indeed V (k 1 , k 2 , k 3 , k 4 ) vanishes because of the momentum-conserving Kronecker Delta in Eq. (6). Therefore, the non-trivial condition is It is easy to see that (sm-8) is immediately satisfied. Indeed, by solving explicitly the constraints and writing k 2 1 + k 2 2 − k 2 3 − (k 1 + k 2 − k 3 ) 2 = −2(k 3 − k 1 )(k 3 − k 2 ), the manifold M 4 can be rewritten as (sm-10) and it is immediate to see that Therefore (sm-1) is solved by (11) with the choice Using the momentum-conserving delta function in V (k 1 , k 2 , k 3 , k 4 ) we can then write where we introduced for k 2 1 + k 2 2 − k 2 3 − (k 1 + k 2 − k 3 ) 2 = 0, extended by continuity at k 3 = k 1 and k 3 = k 2 .

B. Second Order
Repeating the procedure we find that the second order correction to the charge, Q f ;2 , has to fulfil the following condition (sm- 15) In particular, using the quartic form of the potential and the expression (11) where we set f,2 (k 1 , · · · , k 6 ) : and A abc [·] is the operator that antisymmetrises with respect to a, b, and c. This means that Q f ;2 has to have the same structure, namely In order for (sm-19) to solve (sm-15) we have to impose f,2 (k 1 , · · · , k 6 ) . (sm-20) As before, it is easy to see that these equations admit solutions only if f,2 (k 1 , · · · , k 6 ) = 0 (k 1 , · · · , k 6 ) ∈ M 6 , (sm-22) where M 4 is defined in (sm-9) and we introduced It is easy to see that (sm-21) is satisfied for any two-body potential. Indeed, first we note that M 4 can be expressed as in Eq. (sm-10) and then we observe where we used V (k 1 , k 2 , k 3 , k 4 ) = V (k 4 , k 3 , k 2 , k 1 ) * . (sm-25) To treat (sm-22) we note c (6) f,2 (k 1 , · · · , k 6 ) = − A k1k2k3 A k4k5k6 V (k 4 − k 1 )V (k 3 − k 5 ) ∆ǫ(k 1 , k 2 , k 4 ) ∆ 3 f (k 1 , . . . , k 5 ) , k j ∈ M 6 (sm-26) where we used ∆ǫ(k 5 , k 6 , k 3 ) = ∆ǫ(k 1 , k 2 , k 4 ) , k j ∈ M 6 , (sm-27) and defined Furthermore, we note that the latter quantity vanishes identically on the sub-manifold Indeed, it is easy to see that where σ is a permutation and S 3 is the group of permutations of three objects.
This means that in order for (sm-22) to hold for generic functions f one has to impose This is nothing but Condition 1 of the main text. If this condition holds we can write extended by continuity for (k 1 , k 2 , k 3 , k 1 +k 2 −k 3 ) ∈ M 4 and (k 1 , k 2 , k 3 , k 4 , k 5 , k 1 +k 2 +k 3 −k 4 −k 5 ) ∈ M 6 respectively.
To conclude, we have to show that (sm-43) fulfil also the more general condition (sm-39). This is indeed the case as can be seen, for instance, using the symbolic simplification of Mathematica.
Restricting the family (sm-43) to real potentials we have that the only solutions of (sm-39) with regular power series expansion around 0, and hence decaying rapidly enough in real space, are V ab (k) = a 1 − √ bk coth( √ bk) , a, b ∈ R , (sm-47) as claimed in the main text.