Hydrodynamics of local excitations after an interaction quench in 1D cold atomic gases

We discuss the hydrodynamic approach to the study of the time evolution -induced by a quench- of local excitations in one dimension. We focus on interaction quenches: the considered protocol consists in creating a stable localized excitation propagating through the system, and then operating a sudden change of the interaction between the particles. To highlight the effect of the quench, we take the initial excitation to be a soliton. The quench splits the excitation into two packets moving in opposite directions, whose characteristics can be expressed in a universal way. Our treatment allows to describe the internal dynamics of these two packets in terms of the different velocities of their components. We confirm our analytical predictions through numerical simulations performed with the Gross-Pitaevskii equation and with the Calogero model (as an example of long range interactions and solvable with a parabolic confinement). Through the Calogero model we also discuss the effect of an external trapping on the protocol. The hydrodynamic approach shows that there is a difference between the bulk velocities of the propagating packets and the velocities of their peaks: it is possible to discriminate the two quantities, as we show through the comparison between numerical simulations and analytical estimates. In the realizations of the discussed quench protocol in a cold atom experiment, these different velocities are accessible through different measurement procedures.

The outstanding performances of modern experiments in preparing and controlling setups of quantum gases [1,2] have promoted the study of out-of-equilibrium systems in ultracold gases as one of the arguably most challenging topics in the field [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. One of the reasons for the complexity of this challenge is due to the variety of ways in which a system can be driven out of equilibrium and the difficulty in having a general guidance principle to relate their phenomenologies.
Thus, over the years, the community has concentrated on a few protocols that have emerged to be sufficiently clean and interesting. One question at the forefront has been under which conditions (and in which sense) a system is able to reach an equilibrium, and whether this equilibrium can show universal characteristics such as thermalization [19][20][21][22][23][24][25][26][27][28][29].
One of the experimentally most relevant protocol to study out-of-equilibrium dynamics and the issue of thermalization is the quench protocol, in which, typically, a system described by a Hamiltonian H is prepared in its ground state and then, at a given moment of time, let evolved using a different Hamiltonian H [37]. The questions usually asked with quantum quenches refer to late times properties of the systems (and on whether the unitary evolution can lead to something that can locally be described in a semi-classical way) [38][39][40][41][42].
Variations to the usual quench protocol have been considered mainly in two forms: the initial state was taken to be an excited state of the initial Hamiltonian [43], or the interaction was changed only locally [44][45][46][47][48]. In a recent work [35], we introduced a different protocol: we proposed to start with a local excitation and to let it evolve after a global interaction quench. In order to isolate the effect of the quench on the dynamics, the initial state is prepared as a solitonic excitation of the original Hamiltonian, that is, a state for which the one-particle density profile evolves without (almost) changing its shape. This setup has the merit of being experimentally feasible and of showing universal properties already at short times. In [35] it was also shown that, by changing the underlying interaction during the soliton motion, the excitation breaks into two profiles: one moving in the same direction as the initial excitation, the other in the opposite. In a cold atom setup, such a quench can be triggered through the trapping or with an external magnetic field [1] to induce a change in the scattering length and in the speed of sound. The system is then let evolve for short times after the quench and the velocities and shapes of the two chiral profiles created by the quench can be measured either by direct imaging or reconstructed by releasing the trap, through time-of-flight measurements.
As we are going to discuss in this paper, a convenient way to study the dynamics of a local excitation is to employ a (non-linear) hydrodynamics description of the system in terms of the density and velocity fields. Moreover, if the initial density profile is not a big perturbation over the background, velocities and shapes of the transmitted and reflected profiles for short times after the quench can be expressed in a universal way, that is independent from the details of the quench and of the microscopic interaction [35]. The hydrodynamic approach [49] is clearly complementary to a microscopic computation of the dynamics (when practically doable), and it has the advantage that microscopic details of the underlying model enter as parameters of the hydrodynamic equation (e.g., the sound velocity). It also provides a standard tool to study collective excitations and dynamical properties in cold atom setups [50,51] and it applies as well to higher dimension: however in this paper, in line with the topic of the Focus Issue and with our choice of consider solitonic solutions we limit ourself to one dimensional systems and soliton excitations (in higher dimensions, solitonic states are stable only for limited times and one needs to take into account the spreading of the wave packet). We consider in detail the Gross-Pitaevskii (GP) equation, with δ-like interactions (which is relevant for the 1D Bose gases in the limit of small interactions [50]), and the Calogero model, since it is exactly solvable also with in the presence of a parabolic trapping potential [52][53][54].
The same quench setup has also been considered in [55][56][57] with a specific interest on integrable PDEs describing the evolution of one dimensional systems and on the translation of the quench protocol in the corresponding quantum inverse scattering problem. These papers signal the interest of the community in the interaction quench. We would like to stress that the proposed quench protocol could be implemented relatively easily in cold atoms experiments and that, by focusing on the short time universal dynamics, our predictions can be tested directly in the laboratory, which should be contrasted with typical large time results produced by other approaches and protocols. Universal properties of short time out-of-equilibrium dynamics has been considered in [35,[58][59][60].
In this work our main goal is to investigate in detail the hydrodynamic approach for the study of the dynamics of solitonic excitations in one-dimensional systems, clarifying the hypothesis behind the derivations based on the hydrodynamic approach and to discuss on the possible experimental realization of the quench protocol in cold atom systems. We find that in the hydrodynamic approach one can naturally introduce both the "bulk" and "peak" velocities of the transmitted and reflected packets: we provide expressions for both quantities, and we compare them with our numerical simulations for the GP and Calogero models. We also discuss the time scales relevant for the experimental realization of the proposed protocol and the efficacy of our approximation within them.
The plan of the paper is the following. In Section II we review the hydrodynamic approach to the description of one dimensional cold atomic gases [51]. In Section III we analyze the quench protocol and derive the expressions for the velocities and heights of the chiral profiles generated by the quench.

II. THE HYDRODYNAMIC APPROACH
All quantum systems at sufficiently low temperature acquire a collective behavior. In many such cases, the expectation value of the particle density operatorρ( becomes a smooth function and the quantum fluctuations around it are negligible. When this happens, the system dynamics can be described with a hydrodynamic approach [49]. In the following, we will primarily be interested in one dimensional systems with one particle species without additional internal degrees of freedom, so that it is sufficient to introduce the scalar density ρ(x, t) and velocity v(x, t) fields.
Assuming Galilean invariance and locality, the most general hydrodynamic Hamiltonian reads [51] where (ρ) is the internal energy and A(ρ) is related to quantum pressure. We also assume that dissipative effects can be neglected. The equations of motions can be obtained by remembering that the density and velocity are conjugated fields [61,62] satisfying They are the continuity and Euler equations [49]: where ω(ρ) = ∂ ρ ρ (ρ) is the specific enthalpy (which at zero temperature reduces to the chemical potential).
Under the foretold conditions, these equations are completely general: the functions (ρ) and A(ρ) encode the specificity of the system under consideration. A discussion of how one can extract an integrable dynamics out of a generic hydrodynamic system and the derivation of the KdV (used in the next Section) is reported in Appendix A, and useful facts about the KdV are in Appendix B.
and combining the hydrodynamic fields into the complex field the dynamical equations (3,4) can be written as the single complex equation which can be recognized as the 1D non-linear Schrödinger equation, aka the 1D GP equation [without the external potential in (7)]. In Eq. (7) ρ 0 is the density for x → ∞. Details about the the non-linear Schrödinger equation and its KdV reduction are in Appendix C.
The GP equation has been routinely used in the last two decades to describe the dynamics of ultracold bosons at T = 0 [50]. In 3D its validity stems from the fact that there is a condensate, whose macroscopic wavefunction obeys the GP equation, from which one can derive hydrodynamic equations [50]. In 1D there is no condensate, since there is actually a "quasi-condensate" [63]: however, the GP gives a good description of experimental results for small interactions (i.e., for the coupling constant of the Lieb-Liniger model, defined in (C1), 1 [64]) and in particular of the -dark and bright -soliton dynamics observed in 1D experiments [65,66] and of shock waves [67]. When in 1D the coupling constant is not much smaller than 1, one can nevertheless use hydrodynamic classical equations (as the ones given above) to study small deviations from equilibrium [68]. For larger deviations one has to study the quantum dynamics using directly the Lieb-Liniger model (which might be at present rather challenging from the computational point of view) or resort to 1D mean-field effective equations [69][70][71][72] from which hydrodynamic equations may be derived as above producing (time-dependent) non-linear Schrödinger equations with suitable general non-linear terms (ρ), the GP equation corresponding to (ρ) ∝ ρ. As discussed in [35] a good agreement is found between the hydrodynamic results and the GP equation with a power-law (ρ). We observe that such mean-field equations may fail in correctly describing interference between wavepackets, as shown for the Tonks-Girardeau limit γ → ∞ in [73]: in general one expects anyway that the mean-field equations work better at short times with respect to long times. For instance, the ultimate fate of a soliton configuration due to classical to quantum crossover has been studied in [74].

III. THE QUENCH PROTOCOL
In this Section we discuss in detail the quench protocol. The starting point is to prepare the system with a localized excitation. In general, since such state cannot be an eigenstate of a translational invariant system, in time it would diffuse and disperse. However, it is an empirical observation that several systems abruptly taken away from equilibrium, settle back in configurations displaying localized excitations that propagates for long times without degrading significantly. Such behavior is that of a soliton (or trains of solitons) and it is a manifestation of the emergent collective hydrodynamics ensuing in the system. True solitons are a characteristic feature only of integrable differential equations [75]. Nonetheless, individual (or well separated) solitonic waves are commonly observed in a variety of systems, including of course classical systems [76] and ultracold systems [77].
For the quantum many-body Lieb-Liniger model of interacting bosons in 1D it is not obvious what is the quantum content of the microscopical state that results into a soliton.
It makes sense that it would contain a large number of eigenstates such that the coherent dynamics prevents their macroscopic spreading, similarly to what happens in quadratic theories for coherent states. Indeed, the construction of quantum states with one particle density evolving in a solitonic way has not provided conclusive results so far [78][79][80], in the sense that it is not clear what is the definition of a quantum dark soliton from the Bethe ansatz solution, and if the quantum dark soliton may be defined at all (see however the very recent discussion in [81]).
In this work we do not try to solve this problem, but we hope that our analysis can contribute in this direction. Our starting point is the empirical observation that solitonic configurations are commonly excited (and manipulated) in cold atomic Bose [65,66,77,82] and Fermi [83] systems and that they are best understood in terms of an effective semiclassical hydrodynamical description of the system.
Once a solitonic excitation is created, its density profile will remain (approximately) constant and will only translate in space at a constant velocity. In our quench protocol, at some time during the soliton evolution we change the parameter g governing the interaction strength of the system to g . This interaction quench modifies the hydrodynamic equations so that the initial soliton cannot evolve unperturbed anymore. We aim at describing the dynamics, and in particular the short time dynamics, of the soliton after the quench.
To qualitatively illustrate what happens after an interaction quench, we take as an example the GP equation (7) and we change the interaction parameter g to a value g > g at the time t Q . The sound velocity changes from c = √ g to c = √ g (in units where = m = ρ 0 = 1 where ρ 0 is the density for x → ∞). As initial condition we take the gray soliton with velocity V [50,84] reported in Appendix C where γ = 1 − V 2 /c 2 . As illustrated in Fig. 1 there is a time t 2 at which the soliton splits in two (more precisely: at which a second minimum is seen), and a time t 3 in which another splits occurs and one of the two packets in turn splits in two as well. We very close to zero. The black line is the density at t = 0, and the others are at the times t = t 2 (in which a second minimum emerges) -red; t = 2t 2 -green; t = t int -magenta; t = t 3 (in which another minimum is seen in the transmitted packet) -blue. Numerical values are t 2 ≈ 0.708, t int ≈ 8.219 and t 3 ≈ 15.729.
that between t 2 and t 3 the peak velocities are rather stable, while close to t 2 and t 3 numerical fluctuations may be present).
In Fig. 2 we plot the analytical prediction ±c /c obtained from a linear approximation (as discussed below), where the + (−) is for the transmitted (reflected) components. The hydrodynamic approach we introduce in the following reproduces the leading behavior and it is found to be able to reproduce deviations from it. It is also possible to use the hydrodynamic approach to discriminate the bulk and peak velocities: indeed, from the hydrodynamic theory one has that for very short times after the interaction (well before t 2 where a second To set up the hydrodynamic approach for the propagation of the localized excitations after the interaction quench, as we argued above we take the dynamics before the quench to be captured by the equations (3,4), where the hydrodynamic parameters ω(ρ; g) and A(ρ; g) depend on a microscopical parameter g, setting the strength of the inter-particle interaction. The shape of the initial soliton is set by the coupling g and by its velocity V .
The approach we develop in this section is independent from the microscopic details and it applies to general 1D systems with solitonic excitations, and it will be compared with numerical results in the next sections.
The existence of a soliton is a strong indication that it is possible to isolate out of (3,4) an integrable core, with the remaining terms being negligible for a relatively long time. If the initial soliton has velocity V close to the speed of sound c (so that its amplitude is also small compared to the background density ρ 0 ), it is known that the hydrodynamics can be reduced to that of the integrable KdV. We outline the reduction of equations (3,4) to KdV in Appendix A. This procedure, developed for cold atoms in [51], is at the heart of our analysis and thus in the following we assume V c 1.
The KdV is a chiral equation that readṡ where c ≡ ρ 0 ω 0 is the sound velocity (ω 0 ≡ ∂ ρ ω| ρ 0 ) and the nonlinear and dispersive coefficients are given by Details on the derivation of the KdV equation in a generic hydrodynamic system are given in Appendix A, while in Appendix B we collect some basic results on the KdV.
The ± in (9) refers to the two chiralities of the waves while u(x, t) is approximately the density fluctuation over the background u ρ − ρ 0 . The KdV reduction neglects the interaction between the left and right moving sectors. Due to locality, such approximation is violated only when the two chiral profiles overlap. However, since they move with relative velocity of approximately 2c, such effects exists only for short times and hence can be neglected [84].
Suddenly during the evolution of the soliton, we change the interparticle coupling to g .
Thus, after the quench the hydrodynamic equations will change and the effective dynamics will be given by the KdV (9) with modified parameters ζ , α and a new speed of sound c . This sudden change in the interaction is seen by the soliton profile as a perturbation to which it reacts by splitting into a transmitted and reflected profile, exactly as it would happen to a linear wave, in the presence of an obstacle, or if the sound velocity is suddenly changed.
Before the quench, the initial state can be approximated by the KdV soliton (B5) We make the ansatz that after the quench that is, we assume that the quench acts as an external force which splits the soliton into a transmitted and a reflected profile and that, for short times after the quench, these bumps evolve with a given velocity, without changing their shape significantly. If V > 0, the reflected velocity V r < 0 and the transmitted one is V t > 0. Imposing continuity of the solution and conservation of momentum at t = 0 we obtain with Note that, applying the same quench protocol to an initial sound wave moving at the pre-quench speed of sound c would yield transmitted and reflected waves moving at the post-quench sound speed ±c . The above formulae specialize for this case to The derivation of (13,14) is completely general, based only on the ansatz that after the quench the soliton splits into two chiral profiles, and does not depend on the dynamics. The laws governing the evolutions are needed to determine the velocities of the two profiles.
We can estimate the profiles velocities by looking at the motion of their center of mass 1 The velocities are thus where we used the fact that the denominator is the integral of motion I 0 (B8) and thus does not evolve with time. We can trade the time derivative for a spacial one through the post-quench equation of motion (9): 1 We thank the Referee for suggesting this approach.
where the + (−) sign applies to the left (right) moving profile. We can now integrate by part (remembering that at large distances the profiles decay to zero exponentially) to get Note that this expression has a natural interpretation as the ratio between the first two integrals of motion (B8, B9), which can be interpreted as the mass and momentum. Indeed, in (B11) we show that this is the case for a soliton: The integral of motions of the two profiles after the quench are determined immediately at their creation and, remembering (13), are simple rescaling of the original ones: We can now use (20) and (14) to get Solving this system we finally get where we introduced the universal parameter and where the T and R are found consistently to be We note that these expressions are completely universal. Eqs. (25)(26)(27)(28)(29) already appeared in [35] with a different, more qualitative, derivation.
To further highlight the universality of Eqs. (25)(26)(27)(28)(29), we can make everything dimensionless, by measuring velocities in units of sound velocity. Thus, we introduce the reduced velocities: and we write (25,26) as Hence, the solution immediately after the quench is The reflection and transmission coefficients are and the height of each chiral profile is where U is the height of the pre-quench soliton.
Since the chiral profiles are not solitons of the post quench dynamics, their shape will change during the evolution, because different parts of each profile will move at different speeds. For short times after the quench, we can take advantage of the fact that each of the reflected and transmitted profiles are similar to the original soliton, just with a reduced height. Thus, we can use (B15) to estimate the velocities of the different parts: where we used the fact that the width of both profiles is W = 2 α cν and their heights are U r = 3 ζ Rcν and U t = 3 ζ T cν. As a consistency check, we notice that the velocity at the average height of each profile ( u r,t = 1 3 U r,t ) coincides with (31,32). The velocities of the profile peaks are which, in terms of the reduced velocities (30) are where and The parameters β r,t characterizes whether the dynamics of the chiral profiles is dominated by the non-linear term of the KdV (β r,t 1), by the dispersive term (β r,t 1), or is in the solitonic regime of equilibrium between the two (β r,t 1). Eq. (42) is obtained considering that before the quench the dimensionless ratio Ω ζ Ωα in (B4) is close to unity, since we prepared the initial state in a solitonic state, and thus the soliton parameter satisfy Typically, large quenches g g gives β 1 (see, for instance, (52)): we see that the peaks have reduced velocities three times bigger than the profile center of mass. For smaller quenches, dispersive effects will reduce the peak speeds and we see that for β r,t < 1, ν peak r,t < ν r,m . Moreover, for β r,t < 2 3 , the peak starts moving supersonically and thus we expect the profile to become unstable.
Let us now consider the enthalpy ω(ρ) in (4) to be a simple monomial, that is then Hence Note that the reflection and transmission coefficients are the same of the linear process (15), although the velocities are not.

A. Quench Protocol for the Gross-Pitaevskii Equation
To test the predictions based on the KdV universality, we performed some numerical simulations on physically relevant systems. Cold atomic gases with local interaction are commonly described by the GP equation where ρ(x, t) = |ψ(x, t)| 2 and f (ρ) = gρ. To maintain our discussion more general we detail the derivation of the hydrodynamic results for a general f (ρ), but we present numerical results for the GP equation (49)  and in presence of a trapping potential presented in [35] confirms the general validity of the hydrodynamic results for general GP equation with local interactions.
In one dimension, the GP equation reduces to (3,4) with the ansatz (6). The hydrodynamic functions are By substituting the hydrodynamic parameters derived in Appendices C into III for the Furthermore, remembering the definition (42) of the parameter β we find Hence, quenching to stronger interactions (higher speed of sound) takes the dynamics to a non-linearity driven regime, as we anticipated. Using (47,48) we get For the peak velocities (40,41) we have 4 we plot, again for V = 0.96c and V = 0.9c, the analytical predictions for R and T from (53,54) and the numerical GP results: the agreement is excellent, even better than for the peak velocities.
The comparison of Figs. 3 and 4 shows that the hydrodynamic approach gives good results for both bulk and peak velocities, and that the theory can discriminate between them. The data in these figures were taken at time t int defined as halfway between the time t 2 were the two profile become distinguishable after the quench and the time t 3 where the transmitted profile further splits in two. The dynamics around t 3 is clearly beyond our approximation scheme and its treatment requires a different approach, for instance that of [55], which we will not pursue here. Nonetheless, we empirically notice from our numerical experiment that t 3 is more than one order of magnitude greater than t 2 . Moreover, as Using (13), we can estimate the time t 2 at which it is first possible to discern the existence of the two profiles. The comparison with numerical results helps to quantitatively address the validity of the hydrodynamical approach we are following. To do so, we search for the instant at which there is a flex point in (13), that is we look for (x 2 , t 2 ) at which the first and second derivative of (13) vanish: where we introduced with V r,t given by (53,54) and W by (C7).
Once we determine the p and q that solves (57,58), we have Linear combinations of (57), (58) yield the simplified system of equations We can solve (63), for instance, as a quadratic equation in q. One solution is the trivial q = p, which corresponds to the equilibrium reached at x → ±∞ and is thus not the one we are looking for. The other solution is Substituting in (62) we have which can be solved as a cubic equation in p 2 . The only real solution is where we used the definitions of R, T in (53,54) in terms of the quench strength c c . Substituting this solution into (64) and both of them in (61) and further simplifying the resulting expression we find where We tested this prediction against our numerics, finding an excellent agreement and further supporting the validity of our approximation scheme. The comparison is presented in Fig.   6. As ths inset of Fig. 6 shows, it is seen that the results are rather good also for V /c as low as 0.6 and 0.4, where the KdV approximation for the GPE is not supposed to be reliable.
To quantitatively assess deviations from the prediction (67) we observe that the quantity (67) itself is independent from the initial velocity: Numerical results for V /c = 0.96, 0.9, 0.75, 0.6, 0.4 are reported in Fig. 7: as expected, one sees that, while for shallow initial solitons the agreement is excellent, by decreasing the initial velocity the agreement becomes less satisfactory, although still acceptable.
To further show how well and for how long the approximations we used remain valid, we plot in Fig. 8 the density ρ 0 −ρ(x, t) for different times between t 2 and t 3 from the numerical solution of the GP equation and from the KdV hydrodynamic approach. It is seen that the analytical prediction is rather good almost all the way to t 3 , even though approaching t 3 it cannot reproduce the deformation needed to expel the third soliton.

B. Quench protocol for the Calogero Model
In the previous section we showed that numerical simulations strongly support our ana- Newtonian evolution for a system composed of a large number of particles and extract the emerging collective behavior. We will see that the initial soliton will split into a reflected and transmitted density profile, showing an emergent wave behavior out of the individual particle dynamics.
with a dimensionless coupling constant λ. We can perform the numerical evolution in the classical limit [in the quantum Calogero models, the coupling undergoes the quantum shift λ 2 → λ(λ − 1)], which is integrable even in the presence of an external parabolic potential.
This model has a long-range (power-law) interaction and thus it lies in a different universality class, compared to local models. At the same time, as was recognized in Ref. [53], the potential 1/r 2 is also relatively short ranged in one dimension (despite being power law).
Thus, although such interaction has not yet been realized in cold atomic gases, this model provides a good platform to study systems beyond contact interaction, with the additional benefit of being exactly solvable, even with an external confinement.
The above model (Eq. 70) has been shown to admit soliton solutions in [52]. These are very special set of positions and momenta for each particle that collectively evolve as a robust bump on top of a curved density background (see Appendix E).
Due to the power-law interaction, the hydrodynamic description of (70) is not the KdV, but a different integrable equation, in the family of the Benjamin-Ono (BO) equation [85].
It differs from the KdV by its dispersive term and by the fact that its solitons have longer (power-law) tails. Another difference is that this system supports supersonic bright solitons, instead of the subsonic dark ones of local models [52]. We collect in appendix D and E some useful information and new results on the model.
We can calculate the bulk velocities of the reflected and transmitted profiles by operating as before, through the velocity of the center of mass of each. Solving the system of equations (22) for this case yields the same results as for the KdV, with η = 1. To take into account the supersonic nature of the Calogero excitations we define the reduced velocities as so that ν and ν r,t are all positive. In terms of these, we have and which coincides with (53,54).
The solitons in Calogero are Lorentzian (D9) and we notice that their height and width do not depend on the coupling λ. Hence, if β was of the order of unity before the quench, after the quench β r,t will be less than unity, proportional to the coefficients R, T , due to the height reduction. Hence, for the Calogero model, after then quench we will be inevitably in the dispersive regime. We can calculate the peak velocities using (D23), with a r = Rρ 0 τ , a t = T ρ 0 τ : We see that the reflected peak always move slower than the speed of sound (ν r < 0 in our notations), while the transmitted one remains supersonic only for c < 2c. This results is seemingly counterintuitive, since the profile densities lie above the background that their bulk velocities are supersonics. While, in a non-linearity dominated regime, the profile peaks would move faster than the bulk, we see that effect of the dispersive regime effective after the quench in the Calogero model makes the peaks move even slower that speed of sound.
Thus, while the bulk behavior of the Calogero system after the quench follows the same universality of local models, the peak velocities belong to a different universality.
In our numerical simulation, we used the results of [52], where it was shown that a system on N Calogero particle is dual to a system of M interacting complex parameter, where M counts the number of solitons in the system. Thus, to a system of N particle lying on the real axis and interacting through (70) we add a dual particle z(t) (see Appendix E) which draws an ellipse on the complex plane. The interaction of the dual variable with the real particles induces a soliton in the latter and the soliton peak follows the projection of the z particle on the real axis. As explained in [52], a very useful property of this system is that, given an initial condition for the particle position and momenta, the configuration at any given time t can be found exactly by diagonalizing a certain matrix system (by exploiting the Lax pair formulation of the model). In this way, we have been able to set the initial conditions of a soliton, to let it evolve for some time and to follow its evolution after the interaction quench exactly.
In [54] the ground state of the quantum Calogero model in a harmonic potential was studied after a quench and it was observed that the one-particle density starts oscillating and breathing. Such behavior is natural, since after the quench the equilibrium particle distance increases and the trapping cannot compensate for this repulsion. To neutralize this effect, it is thus important to quench at the same time both the interaction and the external potential, so that the background density stays constant. We found that for the model (70) the trapping has to be increased by the same amount of the interaction: Once we are able to stabilize the background in this way, we perform the quench experiment and measure the characteristics of the reflected and transmitted profiles.
To produce the initial soliton configuration, we specify a given complex number z and the initial value problem for the position and momentum of each particle in the system (x j (0), p j (0)) is given by the following equations which can be solved numerically: Having determined the initial values of x j (0) and p j (0), the time dynamics of the system of Calogero particles after the quench can be computed by exploiting the Lax pair formalism, similarly to what was done in [52] to study the dynamics of a hCM soliton. One introduces the following N × N matrices: which depend on time through x j (t) and p j (t). It is straightforward to show that the equations of motion,ẋ are equivalent to the following matrix equationṡ or equivalentlyL written in terms of L and M matrices usually referred to as a Lax pair.
One can then write the solution of the hCM as an eigenvalue problem for a matrix which can be explicitly constructed from the initial positions and velocities of the Calogero particles. Namely, the particle trajectories are given by eigenvalues of the following matrix where the matrices X(0) and L(0) are constructed using the initial conditions x j (0), p j (0) from (77,78) inserted in the definitions (79,80). It is worth pointing out that the above technique is non-iterative in time and hence there is no numerical error accumulation.
In our numerical investigations, we used two different initial soliton velocities: one moving 4% faster than the speed of sound and one 7%. In Fig. 9 we present results for the reduced (reflected and transmitted) peak velocities for both cases. While the numerical results for the reflected peak velocities are clearly closer to the analytical peak estimates than to the bulk ones, the agreement is less satisfactory for the transmitted peak velocities (the 4% data are closer to the peak velocities, while the 7% data are almost in the middle between the peak and the bulk reduced velocities). Nevertheless, these data are clear evidence of the subsonic dynamics predicted by the analytics above, while evidently the peculiarities of the Calogero model renders the quantitative comparison more troublesome. In Fig. 10, we show the comparison between analytics (Eq. 73) and numerics for the reflected and transmitted heights (R, T ) for the cases of ν = 4% and 7%. A remarkable agreement between these analytical predictions and the numerical calculations performed on the Calogero model is evident.

IV. DISCUSSION & CONCLUSIONS
We have discussed a novel quench protocol, in which a moving, localized excitation reacts to a global interaction quench. We focused on a special class of excitations, namely, soliton solutions, which are stable and experimentally achievable in cold atomic systems as well as in other relevant interacting low-dimensional systems. We provided a general hydrodynamic framework describing the collective behavior of such systems. Using this hydrodynamic description, we showed that the dynamics immediately after the quench is universal, in that it does not depend on the details of the microscopic interaction, but only on macroscopic quantities, such a the speed of sound before and after the quench.
The quench protocol can be seen as a initial value problem for the non-linear PDE of hydrodynamic type. This has been the approach of [55][56][57] and is also common in mathematics. It is thus known that a generic initial condition will eventually break into several components, each moving with different velocities (see Fig. 1). Under certain conditions, some of these components will be stable, while other (most of them) will disperse. In [55] it was noted that, for the integrable non-linear Schrödinger, a quench that brings the new speed of sound to be an integer multiple of the original one will only generate a train of solitons (approximately half of them moving in the same direction as the pre-quench soliton and the other half moving in the opposite direction) and no dispersive sound waves.
In our approach we focused on short times after the quench and, making no assumptions on the integrability of the models, but using the structure of the initial condition provided by the quench protocol, we predict that the initial excitation will immediately break into two counter-propagating packets. While they will eventually break further, for a certain time the two chiral packets will retain the shape of the original excitation, but with amplitudes reduced by a reflection and transmission coefficient. The universal form of these coefficients are given in (28,29), while the bulk velocities of the two profiles appear in (31,32). It is also possible to express the latter in dimensionless quantities as in (31,32). These velocities are measurable in a time of flight experiments, but in our numerical simulation we employ a direct measurement scheme. Thus, it is simpler to measure the velocities of the peak of the profiles, which differ from the bulk (that is, center of mass) velocities because of the internal redistribution of energies of the packets. The universal expressions for the peak velocities are given in (40,41). In Figs. 3 and 4 we plot the comparison between the analytical expressions and the numerical results obtained through 1D Gross-Pitaevskii equation, which show a good agreement for the R and T coefficients and an acceptable agreement for the peak velocities. As reported in [35], when expressed in natural units, the agreement between the analytical expressions and the numerical data for the velocities is excellent, although it is hard to discriminate between peak and bulk predictions. We also addressed the experimental feasibility of the proposed protocol: in (67)  We also considered the Calogero model. Although we are not aware of cold atomic gases that realize the inverse square interaction of this integrable system, this model allows us to establish several relevant points. First of all, we see that, despite the non-local nature of this interaction, the short times dynamics remain similar to the local case and we found that the R and T coefficients, as well as the bulk velocities of the chiral profiles, are given by the same universal expression as before (see also [35]). The peak velocities, instead, follow a different universality (which is dominated by dispersive effects, instead of the non-linear at play for local interactions). While the solitons and the chiral profiles move supersonically in the Calogero case, the peak velocities are predicted to be subsonic. The agreement between the analytical expectations and the numerical simulation are less satisfactory in this case, compared to the local interactions, but still support the qualitative behavior derived analytically and the result for the subsonic peak velocities. In our numerical simulation of the Calogero dynamics we employ a classical, Newtonian evolution of the particles and extracted their collective behavior, whose agreement with the hydrodynamic prediction is a further, somewhat independent, proof of the solidity of our approach.
Finally, the Calogero model remains integrable even in the case of an external harmonic potential. Thus, by including the effect of the trap, we established that by simultaneously quenching the interaction between the particles and the external field, it is possible to keep the background density fixed and isolate the quench dynamics of the soliton. We believe that this observation is applicable to the experimental realization of the protocol and give us insights on how to deal with the external trapping. Clearly, we plan to compare our theory with experiments very soon and this will provide input to improve our modeling. A more ambitious future direction would be to study quenches of excited states of quantum systems in regimes where they do not admit a hydrodynamic description and to discuss the behaviour of quantum states that in the hydrodynamical regime behaves like a soliton.

V. ACKNOWLEDGEMENTS
We acknowledge discussions with A. G. Abanov The existence of solitons has been foremost an empirical observation [76]. It then took a while to realize that solitons appear as solutions of integrable differential equations. The reason for which solitonic waves can propagate in actual physical systems is that under general assumptions it is possible to isolate an integrable core in the dynamics, while the rest of the terms can often be neglected up to a given time scale. The typical appearance of this "non-linear universality" is the emergence of the KdV equation as the integrable core of local fluid for shallow waves [76].
In [51] it was reproduced the standard derivation of KdV in classical fluids to extend it to the hydrodynamic treatment of cold quantum systems. We review in this Appendix their approach, to set the notations that we used in the body of this work.
Our starting point are the continuity and Euler equations (3), which assume no dissipation. The results are not modified qualitatively, as long as dissipative effects enter linearly.
We introduce the following notations with respect to the equilibrium density ρ 0 . We want to describe long wave excitations on top of the constant background density ρ 0 . Thus we are looking for solutions in the form where α = ± indicates the wave chirality and where we anticipated the fact that small waves will move with velocities close to the sound speed.
The perturbative expansion is based on a counting scheme with formal counting parameter . We introduce the expansion of velocity and density fields as We substitute this ansatz into (3,4) and collect terms ∼ 3 and ∼ 5 . At cubic order we get and the consistency At the next order, we can combine the two hydrodynamical equations to geṫ where we introduced We also used the identity (A9) so that Shifting (A10) back from the reference frame moving with the sound velocity to the laboratory frame, we have (9).
In this derivation, we considered separately the two chiral sectors. A generic initial condition, however, will consist of both chiralities and, in principle, one should take into account the interaction between them. However, these effects can be neglected since, due to locality, the two sectors interact only when they are overlapping, but this happens for a short time, as they pass through each other with a relative velocity of approximately 2c.
We have thus shown that the dominant non-linear contributions to the dynamics of shallow, long waves in a generic one-component hydrodynamic system (3) is given by (9), which is known as the integrable KdV equation.

Appendix B: Generalities on KdV
The non-linear term ζ in the KdV (9,A10) pushes the different parts of u(x, t) to move with different velocities, so that small perturbations over the background move close to the speed of sound, while the parts more distant from the asymptotic equilibrium ρ 0 move with higher δV . This term tends to generate a shock-like profile, as the tip of u moves faster than the base. The dispersive, α, term, instead, redistributes the kinetic energy within the u profile and effectively lead to a broadening of the u profile.
We can introduce two time scales capturing the effective strength of these two terms where W is the typical width of the disturbance over the background and U is its typical size (note that U has the unit of an inverse length, as it describes the height of a density bump). There are three possible regimes Physically, they correspond to a regime of dominant non-linearity (true non-linear KdV dynamics), dominant dispersion with linear evolution (dispersive waves) and solitonic (equilibrium between non-linearity and dispersion). The dimensionless ratio distinguish the non-linear ( The KdV supports true solitons which keep a perfect equilibrium between the non-linear and dispersive effect and thus propagate without changing their shape. The solitons can be both bright (density higher than background) or gray (a depletion in density). In fact, (A10) is invariant under the simultaneous reversal of the sign and chirality of u ± → −u ∓ .
The single (dark) soliton solutions has the form where δV = c − V is the velocity of the soliton (in the sound velocity frame). The height and width of the depletion depend on the soliton velocity as The integrals of motions are where we rescaled everything by the right factors of ζ, which sets the scale of the amplitude wave, see (9).
We note that the soliton velocity can be found as the ratio between its momentum I 1 and its "mass" I 0 : as can be expected by averaging the nonlinear term in the KdV (A10), see also the comment after (B15).
In section III we considered the evolution of an initial condition which is functionally the same as the soliton (B5), but not necessarily with the correct parameters U, W of a soliton.
The different parts of the KdV equation give: We notice that for an inverse cosh square initial condition every term in the KdV is proportional to its space derivative. Hence, we have For a soliton, U W 2 = 12α ζ and the terms multiplying u x become a constant, which equals the soliton velocity. For generic U and W , the terms within parenthesis give the velocity of each part of the profile.
Note that the center of mass of a profile of the type (B5), that is the height at which the profile has equal area above and below it, is at a third of its total amplitude. In fact the velocity at u = U/3 from (B15) is c − ζ 3 U = V , that is, the bulk velocity of the profile.
The NLSE is also integrable and its single (dark) soliton solutions is [50,84] ψ( We have where the superscript H stands for the Hilbert transform: Due to the long-range nature of the Calogero interaction, its hydrodynamic description does not reduce to a KdV equation, but is given by the so-called double Benjamin-Ono [85], which, for small (chiral) profiles reduces to the usual Benjamin-Ono equation with parameters The one-soliton density profile [88] for the rational Calogero is where τ = c 2 V 2 − c 2 . (D10) Note that for this model, the soliton is bright (positive density displacement) and hence its velocity V > c is supersonic. The height and width of (D9) are given by Let us again consider the evolution of an initial profile like (D9), but with generic parameters, such as We have Hence, the velocity of each part of the profile can be read from the BO equation as We see that for a soliton (a = ρ 0 b) the term proportional to u vanishes and the soliton moves with velocity c 1 + 1 2b , or b = 1 2 Similarly to what we did for the KdV, we can find the velocity of the profile (D13) from the conserved quantities, as in (B11). We have Hence, the average height of the profile (D13) is and its bulk velocity is Eq. (D13) becomes a soliton of (D5) for a = ρ 0 b and in that case (D22) correctly reproduces (D18). The velocity (D22) coincides with the velocity at the center-of-mass, calculated from (D17) at the height (D21), since at this height the dispersive effects are perfectly balanced and cancel out.
The center of mass velocity can be compared with the velocity at the peak by setting u = a/τ 2 in (D17) We see that the initial value b uniquely characterizes the soliton. Combining Eqs. (E4) and (D10), the soliton velocity at the center is given by The soliton velocity changes as it moves away from the center (i.e, as a function of time).
Since the soliton follows the external complex parameter z(t), its velocity matches that of the real part of the dual variable z(t): By inspection of Eq. (E9) and Eq. (E10) we get, As one can notice from the above equation, the soliton velocity decreases as it moves away from the center. It is also worthwhile noticing that the soliton width, y 1 (t), also decreases.
Therefore, we have a scenario where a soliton is moving slower, as it becomes thinner. In flat background, thinner solitons move faster and thus we see that the interesting interplay between the non-constant background and the soliton moving on top it contradicts the common wisdom valid in constant background.