New classical integrable systems from generalized T ¯ T -deformations

.

Introduction.-Since the inception of generalized hydrodynamics (GHD) in 2016 [1,2], there has been a resurgence of interests in understanding the nature of the dynamics in integrable systems.GHD has proven to be a powerful and universal tool that captures it at large scales, and its predictions have been tested against coldatom experiments in different platforms [3][4][5].GHD has also been studied beyond the quantum realm and applied to classical systems including various types of hard rods [6][7][8], the Toda chain [9,10], the nonlinear Schrodinger equation [11,12], the Calogero-Moser model [13] and the sinh-and sine-Gordon models [14][15][16].It provides the statistical framework [17] for the theory of soliton gases [18,19].The structure of GHD is extremely general, and it requires only limited data from the underlying system, such as the two-body scattering shift.
Yet, a full understanding of how GHD emerges from the microscopic dynamics is still lacking.In the hardrods and box-ball systems, rigorous proofs from slowly varying ensembles are available [8,20,21]; in soliton gases they are obtained from finite-gap solutions [18,[22][23][24]; ab initio derivations exist from kinetic theory [25] and Bethe-ansatz semiclassical principles [26]; and the equations of state are well understood [27,28].But every model has specific properties; to understand the universality of GHD, it is important to construct new integrable systems that can access the full space of GHD equations.
In this letter, we do just that.We define a new class of classical many-particle systems with short-range interactions that are integrable, and that cover a very large space of scattering functions.The new systems are shown to arise from generalised T T -deformations.T T -deformations were introduced in relativistic quantum field theory [29][30][31] as integrability-preserving deformations based on local conserved currents, that modify the scattering matrix by "CDD factors" [30,31].Matrix elements of local fields have been recently studied [32][33][34], and T T -deformations have been adapted to systems of different kinds [35][36][37][38][39][40][41], see the review [42].Here "generalised T T -deformations" are those proposed in [43], based on the larger space of extensive conserved charges first studied in the context of the non-equilibrium dynamics of integrable systems [44][45][46].One admits conserved quantities measuring the density of asymptotic momenta, and generalised T T -deformations modify the scattering matrix by an arbitrary momentum function, although no explicit construction was made.Our models provide the explicit construction for classical Galilean particle systems.We confirm that they are Liouville integrable, that many-particle scattering is elastic and factorises into twobody shifts, and that the two-body shift can be chosen as an almost arbitrary function of momenta.
Constructing effective, fully integrable many-particle Hamiltonians for other objects such as solitons in nonlinear media is an old problem, see e.g.[47].Our models are the first to do that, and give in particular the potential for precise initial state preparation in inhomogeneous soliton gases [18,19], a problem of current relevance.
In [39] it was shown that "mass-momentum" T Tdeformations simply change the length of particles.This can also be interpreted as a local change of the effective space particles freely travel through, a special case of a geometric interpretation [48] much like in GHD [49], but the exact local properties of (generalised) T T -deformation are nebulous.We obtain an explicit microscopic dynamics implementing the geometry, generalising the change of particle lengths.Particles are tracers for their asymptotic momenta and "go through" each other: additional available space is implemented by a slowing down at particles' proximity, while reduced space, by a novel process of creation and annihilation of pairs of particles and antiparticles, which effectively gives an acceleration.
We further provide an expression for the free energy, showing at infinite volume the thermodynamic Bethe ansatz (TBA) with Boltzmann-Maxwell statistics (see e.g.[50]); and, remarkably, a similar structure at finite volumes generalising the recent result [51] in hard rods [52] (we are not aware of any other examples).We then show that the GHD equation emerges in the large spacetime limit, in the generality of arbitrary two-body shift.We confirm this by numerical simulations.
The models we introduce are different from most known classical integrable systems, whose dynamics are not of tracer type.They widely generalise hard-rod gases [6,8], and are closely related to the quantum Bethe ansatz and the gas of Bethe wave packets introduced recently [26].We refer to them as semiclassical Bethe systems.Some of the results presented here are proved rigorously in the separate paper [53].
The model.-Consider the N -particle classical phase space with canonical coordinates (y, θ) ∈ R 2N , {y i , θ j } = δ ij -which will be identified with asymptotic coordinates -and the free-particle Hamiltonian H(y, θ) induces a canonical transformation to the "real" coordinates (x, p) as y = ∇ θ Φ [ψ] , p = ∇ x Φ [ψ] : where ∂ x , resp.∂ θ , means derivative with respect to the first, resp.second, argument of ψ(x, θ).The Hamiltonian takes the form where θ i (x, p) are obtained by solving (3) (see below) and the "quasi-potential" V [ψ] (x, p) is defined by the second equation.The trajectories in phase space t → (x(t), p(t)) are induced from the free dynamics y → y(t) = y + θt via the change of coordinates Eqs. ( 2) and (3).Note that there could be multiple trajectories that are admissible when ∂ x ∂ θ ψ(x, θ) is negative.Three statements hold: (i) The Hamiltonian (4) is Liouville integrable: there are N independent conserved quantities, including the Hamiltonian, that Poisson commute with each other and are nice enough functions of phase space.Indeed a natural set is a , a = 1, . . .N .(ii) The quasipotentials are short-range.This is in the weak sense that whenever particles lie on well-separated intervals, a (x Ii , p Ii ).An important consequence is that one can define conserved densities q a (x) such that Q a = dx q a (x), and with {q a (x 1 ), q a ′ (x 2 )} → 0 for x i ∈ A i .Thus, in the limit of N large for physically sensible finitedensity distributions, conserved densities commute at large distances.(iii) The multi-particle scattering processes are elastic -the sets of incoming and outoing momenta are the same -and factorise into two-body scattering processes.The two-body scattering shift for incoming momenta θ 1 , θ 2 is given by φ(θ 1 − θ 2 ) where φ(θ) = lim x→∞ (∂ θ ψ(x, θ)−∂ θ ψ(−x, θ)).Factorised scattering means that in an N -body scattering event, outgoing particle θ is shifted, with respect to the straight trajectory of incoming particle θ, by the sum ω(θ) = θ ′ ̸ =θ sgn(θ ′ − θ)φ(θ − θ ′ ) of the two-body shifts with all particles that it has crossed.Factorised scattering is a fundamental property of many-body integrable systems [54].Here, because ψ(x, θ) becomes constant in x as x → ±∞, at long times p i ∼ θ i , and Thus, θ i are asymptotic momenta and y i are simply related to the impact parameters of the scattering process.Note that each particle i has the same incoming and outgoing momentum θ i : thus it is a tracer for where the asymptotic momentum lies at finite times.We call this a tracer dynamics.
x 2 +α 2 for some α ̸ = 0. Therefore, this is a new family of classical Liouville integrable, short-range, factorised-scattering models, covering a large class of two-body shift functions.Only a few shift functions are known to date to correspond to classical integrable models, hence this is a large extension.Note that fixing the scattering does not fix the dynamics: there is still freedom in ψ(x, θ).In [53], assuming that ψ is twice continuously differentiable, and that ∂ θ ∂ x ψ(x, θ) ≥ 0 (thus φ(θ) ≥ 0) along with a finite-range condition, we rigorously show the statements above, and we show that (2), (3) have unique continuously differentiable solutions.We believe much of this stays true under weaker conditions; but uniqueness can be broken, leading to important physical effects that we discuss below.
By a judicious choice of ψ(x, θ) one can reproduce the hard rod dynamics [6,8,39], see the Supplemental Material [55].The generating function (1) has the structure of a phase Φ for Bethe wave functions Ψ = e iΦ , with Bethe roots θ i and dynamics from semiclassical arguments: p = ∇ x Φ are physical momenta, and y = ∇ θ Φ evolve trivially.In the Lieb-Liniger model, gases of wave packets are indeed described by ψ ll (x, θ) = 1 2 sgn(x)ϕ(θ) with ϕ(θ) the quantum scattering phase [26], and we expect a similar relation for most quantum many-body integrable systems.
Microscopic dynamics.-Theeffect of the interaction can be seen as a particle-dependent, dynamical change of metric from x-to y-space where it is free: the change of infinitesimal length is dy i = K i (x i )dx i where K i (x) = 1 + j̸ =i ∂ θ ∂ x ψ(x − x j , θ i − θ j ) measures the effective "free" space.It can best be pictured in the twoparticle case, see Fig. 1 and [55].For φ(θ) > 0 (e.g.ψ ll ) particles slow down during scattering, giving an effective backwards displacement (Fig. 1a) interpreted as the presence of additional, hidden space where particles must travel.If ψ(x, θ) = 1 2 sgn(x)ϕ(θ) for some ϕ(θ), they "stick" and acquire an internal clock that accounts for this extra space at collisions [26].For φ(θ) < 0, (2) does not necessarily have a unique solution.In case it does, particles speed up, giving an effective forward displacement and a reduction of effective space; hard rods of positive lengths are a limiting case, where the displacementwhich traces the momentum being transferred -is instantaneous.But for most choices of ψ(x, θ) with φ(θ) < 0, solutions to (2) can become multivalued.Then trajectories appear to go backward in time, and the generating function (1), although locally inducing a canonical transformation, globally does not on the standard phase space.One may consider three "regularisations", without affecting the large-scale physics, all implementing a reduction of effective space.(i) Choose any branch; e.g.follow one branch until it disappears, then jump to another branch (Fig. 1b).This is similar to the flea gas [56] (we do not know if there exist choices of ψ(x, θ) and branch that would exactly reproduce the flea gas algorithm).But, like for the flea gas, this regularisation is not Hamiltonian nor time reversible.(ii) Using the "hard-core" picture [39], relabel particles at the first collision (Fig. 1c).This gives a time-reversible dynamics (no longer a tracer dynamics), as long as such collisions always appear before any timebackward parts of trajectories; this re-labelling gives the rods in the hard-rod case.(iii) Inspired by Feynman's picture, interpret time-backward parts of trajectories as antiparticles (Fig. 1d).The proximity of a particle (say orange in the figure) occasions a spontaneous particleantiparticle pair creation (blue and green); the antiparticle (green) later annihilates with the incoming particle (blue) leaving the created particle (also blue) as outgoing physical particle.This is time-symmetric, and we believe it might define a canonical flow on the 'Fock phase-space' F = g R 2N g which admits an arbitrary number g of solutions to (2); however this would need to be investigated further.For smooth ψ(x, θ), multivaluedness can always be interpreted in this way, as the equations ( 2), (3), for x, p, t, define smooth curves in R 2N +1 .We now argue that this picture naturally arises from T T -deformations.
T T -deformations.-GeneralizedT T -deformations, as proposed in [43], are obtained as flows of Hamiltonians H (λ) → H (λ) + δH (λ) parametrized by λ: with w(x, θ) some deformation function.Here q (λ) θ (x) and j (λ) θ (x) are the charge densities and currents, with continuity equation i ) that measures the density of asymptotic momenta at θ in the deformed system.It turns out that if H is an integrable tracer dynamics, then so is H (λ) , and q (λ) θ (x) and j (λ) θ (x) exist and are short-range.Remarkably, starting from a system of free particles H (0) = i p 2 i /2, the deformed Hamiltonian H (λ) is nothing else but H [ψ λ ] , Eq. ( 4), with ).This generalises the massmomentum deformation yielding hard rods [39], see [55].The semiclassical Bethe systems are the first concrete example of generalized T T -deformations.We have rigorous proofs of these statements [53] under conditions guaranteeing invertibility of (2), (3), where H [ψ λ ] can be constructed on the standard phase space.
Going further, here we make the crucial observation that the relation Eq. ( 2), be it invertible or not, is still the correct T T -deformation of the impact-parameterposition relation y i = x i .Indeed, T T -deformations (5) can be obtained as canonical flows [35,39,53,57], and Eq. ( 2) arises directly from applying this flow.Thus, multivaluedness may appear along the T T -flow, and pair creation/annihilation processes occur (Fig. 1d), as claimed; see [55].In all cases, the dynamics remains local.
Thermodynamics.-Let us consider the generalised Gibbs ensembles [58], with Boltzmann weights e − a βaQa .We take more generally (x, θ)-dependent Lagrange parameters varying on scale L, with e − dxdθ β(x/L,θ)q θ (x) = e − i β(xi/L,θi) .In [53] we show, using methods of graph theory [59] and under certain further assumptions, that the free energy density is finite and given by f where the pseudo-energy ε L (x, θ) satisfies (7) This is a TBA-like equation for Maxwell-Boltzmann statistics (see e.g.[50]).The TBA is well known from quantum [60][61][62] and classical [9,10,12,14,15] integrability, at infinite volumes.It is striking that even at finite volumes the free energy possesses a TBA structure; the only result we are aware about this is for (positivelength) hard rods [52,63].We postulate that the finitevolume free energy in interacting quantum and classical integrable systems is given by Eq. ( 7) under an appropriate choice of ∂ θ ∂ x ψ(x, θ) reproducing the scattering shift φ(θ).The infinite-volume limit of Eq. ( 7) yields the expected TBA equation ε , from which lim L→∞ f L follows.Thus the thermodynamics of the semiclassical Bethe systems is described by the standard machinery of TBA, including its "local density approximation".
The free energy gives thermodynamic averages and fluctuations of conserved quantities.Interestingly, we also have the exact thermodynamic average ρ phys (p) = n(θ(p))/(2π) for the physical momentum distribution Here n(θ) = e −ε(θ) is the occupation function and θ(p) is the inverse of the "Dressed momentum" function p Dr (θ).The latter is known to be the physical momentum of an excitation at Bethe root θ in Bethe ansatz systems, and is fixed by TBA equations.See [55].As physical momenta of particles change throughout their trajectories, ρ phys (p) is a quantity that is typically hard to access in integrable models; this is the first exact expression that we are aware of.
GHD.-We provide a heuristic argument for GHD to emerge in the hydrodynamic limit, paralleling Ref. [26]; other techniques [20] should give rigorous results.
We take macroscopic space and time, x = Lx, t = L t (x, t finite, L → ∞), with scaled coordinates xi ( t) = x i (t)/L and ȳi = y i /L.The empirical density ρ e (θ, x, t) = L −1 i δ(x − xi ( t))δ(θ − θ i ), is assumed to converge "weakly": ρ e (θ, x, t) → ρ p (θ, x, t).Clearly, , we assume that, for every i, there is a fraction of particles j that tend to 1 as L → ∞ such that xi ( t) − xj ( t) > 0.Then, we can replace ∂ θ ψ(Lx, θ) → ψ sgn(x) φ(θ) where (9) Making the ansatz ẋi = f (x i , θ i ), the second term on the right-hand side is This is exactly the equation for the effective velocity in GHD [1,2].Putting this into Eq.( 8) and taking the limit L → ∞ we obtain the GHD equation, Thus, we have shown that, at large scales, the semiclassical Bethe system satisfies the GHD equation.This is not rigorous -in particular, in Eq. ( 9) one would need to use a regularization of the delta-function.We give an alternative derivation of GHD in [55], based on the fact that the metric change dy i = K i (x i )dx i converges to the GHD change of metric K i (x) → 2πρ s (θ i , x) determined by the "space" or "total" density ρ s (θ, x) [49].We numerically demonstrate that the GHD equation correctly captures the large-scale behaviour in an explicit example, see Fig. 2. For illustration, we use the phase shift from the quantum Lieb-Liniger model, but with an initial state that breaks the maximal fermionic occupation allowed by quantum mechanics (the maximal density of particle per state is 6.264 > 1).This initial state is nevertheless realizable, and its hydrodynamics makes sense, as indeed it is realised by a semiclassical Bethe system.The details of the numerical simulations can be found in [55].Compared to the evolution of noninteracting particles the expansion of the interacting particles is much slower, which is in line with the intuitive meaning of a positive phase-shift as an effective timedelay during the scattering of two particles.
Conclusions.-Weintroduced a new class of classical integrable models, dubbed semiclassical Bethe systems for their relation with the quantum Bethe ansatz, obtained as generalized T T -deformations of classical noninteracting particles.In these systems, each particle is a "tracer": it has the same incoming and outgoing momentum.The class is parametrised by a function determining the microscopic dynamics, and displays factorised scattering with a (largely) arbitrary two-body shift, including those found in many quantum integrable models.The microscopic dynamics displays special features including pair creations / annihilations; the thermodynamics in finite volumes surprisingly takes a form akin to the thermodynamic Bethe ansatz (TBA), reducing to the TBA at infinite volume; the distribution in physical phase space can be evaluated exactly; and the large-scale dynamics is described by GHD, and therefore identical to that of any quantum/classical integrable system with the same chosen two-body shift.We conjecture that, with shortrange interaction, the agreement persists at higher orders: the models should encode the universal hydrodynamic expansion of classical many-body integrability, as corrections due to specific interactions should be exponentially subleading.For instance, particles' positions in our models should approximate well the spatial distribution of solitons in dense soliton gases, something that can be useful for initial state preparation.
It would be interesting to construct the full particlenon-conserving Hamiltonian description of trajectories (2), (3) with negative shifts φ(θ) < 0, Fig. 1d.Finding the full integrability structure of our models, perhaps connecting with sine-Gordon soliton trajectories [64][65][66], would be interesting, as would quantising our models, perhaps in the spirit of [26] (integrability of generalized T T -deformed systems is established [43]); the notion of pair creation / annihilation may play an important role.Finally, adding an external potential is possible, and we anticipate that rigorous proofs of the emergence of the GHD equation can be obtained following ideas in the hard-rod case [8,20].This might also shed light on GHD beyond the Euler scale.
In this supplemental material, we refer to equations in the main text as (xMT) where x is the equation number.
If α is too small the inverse function of f becomes multivalued.In this case the trajectory becomes nonunique during scattering (see Fig. 1).

II. SPECIALISATION TO HARD RODS
For hard rods of positive length λ [1], one may take the continuous piecewise linear The particles then trace the hard rods' momenta.Indeed, this guarantees that whenever particles' positions are separated by distances greater then λ, they move freely, while the first pair of particles coming to a distance λ experiences an instantaneous motion pass each other, effectively exchanging their positions.As almost surely (with respect to physically sensible distributions of positions and momenta) only at most one pair of particles at a time come to a distance λ, only at most one two-body process occurs at a time.This indeed reproduces the hard-rods dynamics.This should be seen as the limit ϵ → 0 of the choice This choice gives invertible Eqs.(2MT), (3MT) on all configurations where at most one pair of particles are a distance smaller than (1 + ϵ)λ, and as ϵ → 0, at a collision of a pair of particles the speeds of the particles tend to infinity, giving, in the limit, the instantaneous exchange of their positions.
of the particles stuck with each other simply takes the average momentum [4], in fact, this model represents the center-of-mass position of any group of rods in the wrong order; the microscopic dynamics is therefore slightly different.
In the case of positive-length rods, in fact there is one subtlety concerning the dynamics with ψ hr,+ .Let us refer to the "normal sector" as the set of configurations under the constraint that all but at most one pair of particles are a distance strictly greater than λ from each other (or (1 + ϵ)λ in the regularised version).Although the solutions to Eqs. (2MT), (3MT) in are indeed unique in the normal sector, there are in general other solutions where many particles are near to each other.These "ghost solutions" are not required for a well-defined dynamics when starting with configurations where all particles are far enough from each other, because, as explained above, almost surely the normal sector only is explored over time.However, if admitting all solutions, as per our general particle-production process, these solutions should be considered; they produce "vacuum" loops.Such loops in fact allow us to define the hard-rod dynamics, for positive rod lengths λ, with rods being closer than a distance λ, thus in particular at densities higher than 1/λ.Therefore, we see that the hard-rod model with positive rod lengths is that given by the restriction to the normal sector of the semiclassical Bethe systems with ψ hr,+ ; this sector is almost surely invariant under the dynamics, hence a consistent restriction.
We note that the choice w(x, θ) = 1 2 δ(x)θ in (5MT), for the generalised T T -deformation, corresponds to the mass-momentum deformation, argued in [2] to give rise to the hard rods.From our analysis, we see that a number of subtleties arise.
First, in [2], it is the analytical continuation in λ, from positive rod lengths, that was argued to give negativelength rods.Here, we see that, as per our discussion above, the direct deformation towards negative rod lengths gives ψ hr,− , the "center-of-mass" version where for groups of rods in the wrong order, we only trace the momentum of the center-of-mass position.Second, the choice w(x, θ) = 1 2 δ(x)θ, in the direction of the T Tflow giving positive rod lengths, in fact only gives the true positive-length hard rods under the hard-core regularisation where exchanges are made at first collisions (see the main text); this was indeed the picture taken in [2].Here, we see that one must choose a λ-dependent w hr,+ (x, θ) = 1 4λ Θ(|x| − λ) in order for T T -deformation to give the semiclassical Bethe system with ψ hr,+ that directly represent, when restricting to the normal sector, the tracers of hard rods' momenta (the go-through picture).
Finally, the specialisation of the finite-volume TBA equation (7MT) to the center-of-mass negative-length hard rods, with ψ hr,− , is It is interesting that we find no correction from the local structure at all: the same TBA equation arises at infinite volume.As (7MT) is derived in [5] under assumptions that are not satisfied for the positive-length hard rods, ψ hr,+ , its specialisation to this case is not expected to give the correct answer (this turns out to be of a similar form to that of [6,7], but sligthly different).

III. DETAILS ON THE NUMERICAL SIMULATIONS
In the letter we simulate the Euler scale time evolution of a macrosopic initial state described by: using both the underlying microscopic particle model and the GHD equation.In the following we provide details on both simulations.
A. Particle simulations

Generating the initial state
In order to generate an initial state whose quasiparticle density corresponds to Eq. ( 3), we use an adaption of an algorithm that is used to generate hard rods initial states [8].The algorithm first generates the positions of particles using their total density ρ0 (x) = dθρ 0 (x, θ) = 10 √ 2π e −x 2 /2 and then subsequently randomly chooses their initial θ according to the double Gaussian momentum distribution f (θ) = ρ0(x,p) ρ0(x) = 1 5 √ 2π e −25(θ−1) 2 /2 + e −25(θ+1) 2 /2 .The procedure for generating the particle positions is as follows: Fix the location of the first particle (we choose x 1 = 0).The locations of all the particles for x > 0 are generated consecutively using the following rule: where η n ∼ Exp(Lρ 0 (x n )) are exponentially distributed random numbers.Note that this can be viewed as a discrete time random walk.As L → ∞ the η n = O(1/L) become very small, meaning one can ignore locally the x dependence of ρ0 (x n ).This stochastic process then locally generates a Poisson point process with the average density ρ0 (x).The slow dependence of ρ0 (x) on x then creates variations of the density on the macroscopic scale.
For practical purposes the procedure has to be stopped eventually, which we do once the density ρ0 (x) < 0.01 is below a threshold.The above procedure generates the particle locations for x > 0. In order to generate the negative particle locations x < 0 as well, we simply repeat the algorithm, but run it in the negative direction.

B. Time evolution of the particles
The particles can then be time-evolved to any time t by solving Eq. (2MT) for given y i (t) = y i − θ i t numerically (the initial y i can be computed from Eq. (2MT)).Numerically solving non-linear equations is a standard, but still non-trivial task.Appropriate methods could vary depending on the model, but here we choose to use the gradient descent algorithm.Its validity is guaranteed for models with positive phase-shifts, like the Lieb-Liniger phase shift -more precisely, whenever ∂ x ∂ θ ψ(x, θ ≥ 0. It is based on the observation that Eq. (2MT) is the minimization condition of the following convex function: where ∂ x γ(x, θ) = ∂ θ ψ(x, θ).In our case γ(x, θ) = √ x 2 + α 2 φ(θ).
Since S(x 1 , . . ., x n ) is a convex function the gradient descent algorithm for finding its minimizer is guaranteed to converge.

IV. CHANGE OF METRIC AND PHYSICAL MOMENTUM DISTRIBUTION
Here we describe a different method to derive the GHD equation which is based on the change of coordinates Eq. (2MT).We can rewrite Eq. (2MT) directly in terms of the empirical density as (omitting the time argument) where in particular we use the fact that ψ θ (0, 0) = 0 by anti-symmetry.Assuming that xi = x(θ i , ȳi ) for some well-behaved function x(θ, ȳ) and taking the large- The differential at fixed θ is dȳ = (1+ dθ ′ ρ p (θ ′ , x)φ(θ−θ ′ ))dx = 2πρ s (θ, x)dx (8) where ρ s (θ, x) is the spatial, or total, density (with the conventional normalisation used in quantum systems).Eq. ( 8) is the transformation of coordinates x → ȳ that is known to trivialise the GHD equation [11].This trivialised GHD equation can be solved immediately and its solution can be mapped to the actual particle density by inverting Eq. (7).Thus, as the ȳi coordinates indeed evolve trivially, this can form a basis for an alternative way of proving the emergence of the GHD equation.
In fact, a similar calculation may be made to obtain the physical momentum density.Using where ϕ(θ) = θ 0 dθ ′ φ(θ ′ ), we get from (3MT) from which we deduce p i = p Dr (θ i , xi ) where p Dr (θ, x) is the Dressed momentum (see e.g.[12]); the latter is the physical momentum variation upon adding a particle of rapidity θ in Bethe ansatz systems (here in the state at macroscopic position x).Therefore, changing variable, the physical momentum distribution of our classical particle system is ) where we used 2πρ s (θ, x) = ∂p Dr (θ, x)/∂θ, and θ(•, x) is the inverse of p Dr (•, x), i.e. p Dr (θ(p, x), x) = p.Using the expression for the occupation function n(θ, x) = ρ p (θ, x)/ρ s (θ, x), this reproduces the formula quoted in the main text.

V. TRAJECTORIES FROM T T -DEFORMATIONS
In [5] we prove that the Hamiltonian Eq. (4MT) is well defined and arises from a generalised T T -deformation of the system of free particles as per Eq.(6MT), under certain conditions on ψ, including ∂ x ∂ θ ψ(x, θ) ≥ 0. There, we show that these conditions are sufficient to guarantee that Eqs.(2MT), (3MT) have unique solutions for x and θ, respectively.If the conditions is broken, there is in fact no guarantee that Eq. (3MT) can be inverted to define the Hamiltonian as a function of x, p.Likewise, there is no guarantee that Eq. (2MT) can be inverted to provide the trajectories.Nevertheless, we discussed in the main text how to make sense of trajectories in cases where invertibility does not hold: Eq. (2MT) still defines the trajectories, but one must interpret it correctly; one interpretation is that one must modify phase space in order to account for different particle numbers, and particle creation / annihilation occurs.
Here, we show that, indeed, T T -deformations of freeparticle trajectories give Eq.(2MT) without further conditions on ψ(x, θ): any solution to Eq. (2MT) is compatible with the corresponding generalised T T -deformation, even when Eq. (2MT) is not invertible.Thus, the interpretation above is sensible; in particular, one may augment the phase space, and at critical values of λ, where multiple solutions arise for some times, it is consistent to consider these multiple solutions at once, interpreted as the creation of particles and anti-particles.A similar construction from Eq. (3MT), to construct the full Hamiltonian on this larger phase space, is left for future works.
(12) This is in the sense that the Hamiltonian is deformed as H λ → H λ + δλ{H λ , X λ }.In fact, any conserved quantity transforms in this way, including the asymptotic momenta θ λ i (x, p).
As the flow in λ is a Poisson flow, it preserves the Poisson algebra, thus it is natural to define the asymptotic impact parameters as those obtained from his flow: y λ i (x, p) → y λ i (x, p) + δλ {y λ i (x, p), X λ (x, p)}.(13) We now show that: if Eq. (2MT) holds, in which we replace ψ → ψ λ = λψ, as well as y i , θ i → y λ i , θ λ i , and if θ λ i 's satisfy the flow equation, then the impact parameters y λ i 's satisfy the flow equation (13), where the flow is defined using w = ∂ x ψ/2 (this is in agreement with Eq. (6MT)).Therefore, any solution to Eq. (2MT) is compatible with the corresponding generalised T T -deformation.
Note that we do not need to define the Hamiltonian for the above statement to make sense: we directly use the flow on trajectories, much as was done in [2].Note also that under the above identification, we have the asymptotic condition |x|w(x, θ) → 0 (|x| → ∞), thus the function X λ (x, p) is well defined.In fact it takes the form where here and below we omit the superscript λ and use the condensed notation x ij = x i − x j , etc., for lightness of notation.The proof is as follows.Under the flow, the left-hand side of Eq. (2MT) gives (we use indices x, θ to denote derivatives) {y i , X} = jk {y i , x j }ψ x (x jk , θ jk )+ j ψ θ (x ij , θ ij ).(15) We may use Eq.(2MT) to obtain {y i , x j } = l ψ θθ (x il , θ il ){θ il , x j } (16) whence {y i , X} = jkl ψ θθ (x il , θ il )ψ x (x jk , θ jk ){θ il , x j }+ j ψ θ (x ij , θ ij ).

FIG. 1 .√ x 2
FIG. 1. Trajectories of two particles during scattering for (a) φ(θ) > 0 and (b), (c), (d) φ(θ) < 0. The thin dashed lines are the asymptotic trajectories.We show three interpretations of the multivalued solutions for φ(θ) < 0: (b) Particles jump instantaneously.(c) Particles are relabelled during collisions.(d) Spontaneous creation of a particle-antiparticle pairs: starting from the bottom, at a certain time (marked by a star) the blue and orange particles are close enough to allow for two pair creations (there are 3 solutions: two blue-orange (outer), one red-green (inner)).Later one particle of each pair (the red and green) each annihilates one of the original particles (marked by a circle), leaving the new blue and orange particles as outgoing.We used ψ(x, θ) = x √ x 2 +α 2 2 arctan θ c , with c = 0.5, α = 0.4 for (a) and c = −1, α = 1 for (b), (c), (d).