Universal shock-wave propagation in one-dimensional Bose fluids

We propose a protocol for creating moving, long-lived shock waves in interacting one-dimensional Bose fluids. The fluid is prepared in a moving state by phase imprinting and sent against the walls of a box trap. We demonstrate that the shock wave front thus formed remains robust against dephasing and thermal fluctuations, over several oscillation periods. We show that this large amplitude dynamics is universal across the whole spectrum of the interatomic interaction strength, from weak to strong interactions, and it is fully controlled by the sound velocity inside the fluid. The shock waves we propose are within reach for ultracold atom experiments.

Large-amplitude moving perturbations are found in all types of fluids, and even in solids. As a response to a sudden change of parameters, a shock wave-a sharp jump in hydrodynamic variables capable of propagating without dispersion-may form. Even ideal fluids can support shock waves as long as the infinitely sharp discontinuities are consistent with the conservation laws. Dissipative effects, present in real-world fluids, give the shock layer a thickness and a shape [1]. Superfluids can host shock waves, within the corresponding hydrodynamic two-fluid theory, as, e.g., in the case of Helium-4 [2,3]. Shock waves were also experimentally observed in dilute, weakly interacting Bose-Einstein condensates of ultracold atoms [4][5][6][7].
One-dimensional (1D) Bose fluids constitute particularly suitable media for a study of shock waves. Only collective modes are possible in such reduced dimensionality, and the fluids belong to the Luttinger liquid universality class [8], thus opening a possibility for a unified theory. Furthermore, at strong interactions, onedimensional Bose gases display a statistical transmutation, i.e., some of their properties coincide with those of an ideal, and thus easily tractable Fermi gas [9]. In addition, several theoretical methods are available in the full spectrum of the interaction strength [10][11][12], thanks to the integrability of the underlying model [13].
In a strongly interacting 1D Bose gas, shock waves of a specific class were predicted to form as a result of the time evolution following a density bump in the density profile [14]. This type of shock waves maps to a solution breakdown in the nonlinear transport equation (also known as inviscid Burgers' equation). This type of solution is well supported by the hydrodynamic equations, but the effects beyond hydrodynamics-e.g., quantum pressure and microscopic correlations-are not capable of stabilizing this shock wave and it disappears after the collapse of the underlying nonlinear transport equation.
In our work, we propose a dynamical protocol for generating the more conventional type of shock waves that are well stabilized by the beyond-hydrodynamics effects, across all parameter regimes. We obtain them by a large perturbation on top of a moving quantum fluid in one dimension. The shock-wave front is created when the fluid, with an initially imprinted velocity, hits against the walls of a box trap. We then follow the large-amplitude dynamics and in particular the velocity and the shape of the wavefront, as well as the current across the sample.
By combining three theoretical methods, i.e., classical field theory, generalized hydrodynamics and exact solution we describe all interaction strengths from weak to strong repulsion. We observe a remarkably robust behavior of the shock wave propagation at all interaction regimes: we find a universal trend for the wavefront in form of a stable step-like flow and of the current, which displays a triangular-shape oscillation. Both features are robust under inclusion of thermal fluctuations. We also show why our shock waves are not collapsing, at difference from [14]: the underlying equation is not a transport equation, rather, it is a free-wave (or massless Klein-Gordon) equation. Our microscopic approaches evidence also non-universal features which depend on the interaction strength: at weak interactions, formation of density modulations due to emissions of phonons and soliton trains, and at large interactions density modulations associated to the Friedel-like oscillations in proximity of a wall, due to quantum fluctuations of the density.
Model We consider a one-dimensional Bose gas with repulsive interactions described by the Lieb-Liniger Hamiltonian: where m is the mass of the particles, g is the onedimensional interaction strength [15] describing the collisions in a tight atomic waveguide. V (z) is a box trap potential of size L with infinitely high walls which we model by imposing hard-wall boundary conditions andΨ(z), Ψ † (z) are bosonic field operators satisfying the commu-tation relations [Ψ(z),Ψ † (z )] = δ(z − z ). The trap contains a fixed number of particles N = L 0 dz Ψ †Ψ . We define the dimensionless coupling strength γ = gm/ 2 n 0 , n 0 = N/L being the average fluid density.
We study the dynamics following a quench in momentum space: starting from the equilibrium state, at time t = 0 we apply a phase imprinting to all particles, generated by the shift operatorÛ = e ik0ẑ , yielding a boost of all the particles with velocity v = k 0 /m. As will be shown below the relevant velocity scale is the speed of sound c(γ) [16]. We then follow the quantum dynamics of the particle density n(z, t) = Ψ † (z, t)Ψ(z, t) and of the spatial average of the current density The system under consideration is, in general, exactly solvable by Bethe-Ansatz [17], however a quench dynamics can be difficult to compute, requiring to evaluate overlaps of excited-state Bethe wavefunctions. Thus, in order to cover the whole interaction range we use three complementary theoretical approaches: the mean-field Gross-Pitaevskii (GP) equation for the weakly interacting gas, the Generalized Hydrodynamic (GHD) theory for intermediate interactions [11,12,18] and the time-dependent Bose-Fermi mapping [9,19] for the strongly interacting Tonks-Girardeau (TG) limit.
The Gross-Pitaevskii equation describes the time evolution of the condensate wavefunction ψ(z, t) [20] by the nonlinear Schrödinger equation We solve it numerically by time evolving the initial equilibrium solution ψ eq (z), satisfying the box boundary conditions ψ eq (0) = ψ eq (L) = 0, boosted by the phase imprinting ψ(z, t = 0) = e ik0z ψ eq (z). The Gross-Pitaevskii equation breaks down at intermediate interactions, when quantum fluctuations start to significantly affect the dynamics and modify the equation of state of the Bose fluid. In this regime, we describe the fluid at long wavelengths using the generalized hydrodynamic equations for the distribution function n(z, k, t) of the quasiparticles of the Lieb-Liniger model solved self-consistently with the equation for the dressed velocity v eff n (k) = m where the dressing operation is defined by h dr (k) − scale associated to the interaction strength. To implement the quench and impose hard wall boundary conditions we use a mirror image method [21]. Once the self-consistent solution n(z, k, t) is found, we compute the current density according to j(z, t) = dk 2π k m n(z, k, t) and the total current J(t) = L 0 dz j(z, t). Finally, in the Tonks-Girardeau regime of infinitely strongly interacting bosons, we describe the dynamics using an exact solution based on the time-dependent Bose-Fermi mapping [22][23][24], where the many-body wavefunction Ψ T G (z 1 , ...., z N ) reads where ψ (z, t) is the solution of the single-particle Schrödinger equation i ∂ t ψ = − 2 ∂ 2 z /2m + V (z) ψ with the initial conditions ψ (z, 0) = ψ 0 (z), where ψ 0 (z) is the eigenfunction of the Schrödinger equation at initial time, boosted by the phase imprinting [21]. This approach allows us to describe in an exact way the full dynamics after the quantum quench.
Results Figure 1 shows the universal behavior of the density dynamics in a one-dimensional Bose gas in a box trap. At early times, the density develops a double step profile, corresponding to the shock wave: as the particles are moving towards one side of the box and bounce on the boundary, a high-density plateau develops upstream, while a low-density plateau develops downstream as particles move away from the other boundary. In between the density remains unchanged, until the two plateaus meet and the total current vanishes. At this point the two propagating fronts cross each other (see FIG. 2. (Color online) (a) Particle density n(z) (normalized to N/L) and (b) local current j(z) (normalized to c(γ)N/L) in the box at different times (t = τ × L/c(γ)), using GP (blue solid lines), GHD (orange solid lines) and TG (yellow solid lines) calculations, for the same parameters as in Fig. 1. The dashed black lines are the predictions of the step density profile model (see (7) and (8)). The light gray triangle is a guide to the eye emphasizing the common propagation velocity of the fronts. Particle densities (local currents) at successive times are shifted downwards by 0.15N/L (0.15c(γ)N/L) for clarity.
role of the two boundaries being exchanged because the flow is now reversed, the two plateaus develop again, and so on. As shown in Fig. 2(a), focusing on the early time evolution, the two fronts separating the density plateaus propagate at the speed of sound in the fluid at that fluid density, such that if time is rescaled by L/c(γ) the density cuts as a function of time fall onto each other displaying a remarkable universal dynamical behavior. In addition to universal features, we notice also small differences among the three regimes [25]: at weak interactions, in addition to a shock wave, we observe the formation of soliton trains upstream of the flow, see, e.g., Fig. 1 at times t = 2.5L/c(γ), as we have checked by analyzing the phase of the condensate wavefunction [21], and also reported in [26,27]. At very large interactions, we observe modulations in the density profile due to the effect of the box confinement, which correspond to Friedel oscillations of the mapped Fermi gas. These are present generally in strongly correlated fluids and are due to the quantum fluctuations of the density [8,28].
In order to gain insight into this dynamics, we analyze the shock wave formation at weak interactions, where we show that the long-wavelength limit of the Gross-Pitaevskii equation supports shock waves, providing the observed wave speed (i.e., the speed of sound) and the correct relationship between the velocity and density discontinuities. Assuming a weak perturbation on top of a density n 0 and neglecting box boundaries we set ψ(z, t) = √ n 0 + δψ(z, t) e −iµ0t/ , with µ 0 ≡ gn 0 .
Using (2), the equation to first order in δψ reads: where c GP ≡ µ 0 /m is the speed of sound, and ξ ≡ /mc GP is the healing length. Note that equation (6) is nothing else but a convenient rewriting of the standard Bogoliubov equations for small excitations of a homogeneous stationary Bose condensate. In the longwavelength limit, the typical length scale of the variation of the excitation wavefunction δψ(z, t) is much greater than the healing length ξ. To the zeroth order in ξ, Eq. (6) reduces to the free-wave or zero-mass Klein-Gordon equation 1 ∂z 2 δψ = 0, which admits as general solution a linear combination of right and left moving wavepackets We search for a δψ(z, t) of a form (7), that on one hand corresponds to a sharp right-moving jump in the field variables, and on the other hand constitutes a small excitation limit of two distinct exact solutions of the Gross-Pitaevskii equation (2) (with V (z) = 0) on either sides of the discontinuity. Consider the following Ansatz, corresponding, for simplicity, to a single boundary on the left at z = 0: (8) The discontinuity is created at t = 0, and it propagates in the positive direction, with the speed of sound. To the left of the discontinuity, we have a stationary condensate featuring a density depression ∆n, relatively to the "reference" density n 0 ; to the right, the condensate is moving with a velocity v 0 , at a density restored to its reference value n 0 . Linearization of the above for small density depression ∆n, velocity v, and phase ∆φ allows to write δψ in the form of Eq. (7) provided that [21] We have verified that this condition is satisfied in our simulations for all interaction regimes. It is remarkable that despite being very different microscopically, the three models behave identically in the long wavelength domain. Indeed, as our analysis of the GP model in the infrared limit given above indicates, the shock waves, in spite of being sharp features, are long wavelength phenomena, and as such they replicate the universality that one-dimensional gases with sound enjoy [8]. There is no contradiction here: within a long wavelength theory, the shock waves appear as structureless discontinuities, themselves fully regulated by the observables residing within the infrared domain. Note that more detailed properties of the shock waves, such as the density and velocity profiles, are governed by the ultraviolet part of the momentum spectrum, and they remain model-specific. A similar phenomenon-an availability of a self-sufficient description of the shock-wave-induced discontinuities and their dynamics by a theory incapable of resolving their width-is observed in dissipationless hydrodynamics (see [1], §82), with the discontinuity width and shape being specific to particular values of the viscosity and thermal conductivity coefficients [1] ( §87 there). Figure 3 shows the oscillations of the current at longer times obtained from the three theoretical approaches: GP (γ 1), GHD at γ = 1, and TG (γ → ∞). We have also used GHD to investigate the dynamics at the hydrodynamic scale for the whole interaction strength range. We find a good agreement for the current dynamics with the TG exact solution at large γ and with the GPE at small γ. We also obtain the period of the current oscillations, as shown in Fig. 3(b). We find that the period is well accounted for by the expression T = L/c(γ) where c(γ) is the exact speed of sound obtained from the solution of the Lieb-Liniger model [16]: this provides another confirmation that even though the shock wave is generated by a large-amplitude oscillation, its hydrodynamic nature implies that the speed of sound sets its dynamics.
Our microscopic calculation finally allows us to address the robustness of the shock waves created by the proposed protocol. At long times, as illustrated in Fig. 1 and in Fig. 2(a), one wavefront (the mid-to-low density one, [21]) gradually broadens during the propagation, resulting in a loss of contrast between the two density regions. Correspondingly, the oscillations of the current progressively damp and change shape from triangular to sinusoidal, see Fig. 3(a). The damping of the current oscillations weakly increases with interaction strength, and can be estimated within GHD by the dephasing time τ d = L/(v high − v low ) ∼ (n 0 /∆n) × L/c(γ) and hence is faster for stronger quenches. Here v high (v low ) corresponds to the effective velocity of the fastest (slowest) quasi-particles involved in the dynamics [21]. Its microscopic origin depends on the interaction regime: at weak interactions, it is due to the mode-mode coupling induced by the nonlinearity in the GPE [29], at strong interactions it is due to the slightly different dispersion of each single-particle mode with time.
We have also explored the effect of thermal fluctuations in the propagation of the shock waves. We find that the phenomenon persists at finite temperature up to T ∼ µ/k B , with µ being the chemical potential, and that the damping of the current oscillations increases with temperature (see [21] for details).
Conclusions. We have proposed a protocol for generating shock waves in a 1D Bose fluid: we use phase imprinting to impart a velocity flow onto the gas, driving it against the walls of the container. By combining several theoretical techniques, we have shown that the formed wavefront is stable, and it is capable of propagating over several periods of oscillations in the box trap; the effect persists for any interaction strength, from weak to strong repulsion, and it is robust against thermal fluctuations. We have found that even under such a strong quench the wavefront follows a universal dynamics fixed by the hydrodynamic sound velocity. From the theoretical point of view this means that the underlying microscopic theory supports the universal features and keeps them stable: the large-amplitude dynamics is fully consistent with infrared hydrodynamic regime, and does not depend on short-distance cutoff except for the details of the shape of the wavefront.
Our work calls for further studies on the dynamics at long times, e.g., exploration of the emergence of grey solitons in the weakly interacting regime and their analogues at strong interactions, and of the origin of the damping mechanisms in one dimension. More generally, our work constitutes a new avenue towards the theoretical and experimental study of strongly driven one-dimensional quantum systems, allowing for an access to quantum turbulence. Finally, our result implies an existence of a new kind of universality in out-of-equilibrium dynamics.
Note added. After completing this work we became aware of Ref. [30] which studies dispersive shock waves of the type proposed in [14].
We acknowledge fruitful discussions with J. Dubail on generalized hydrodynamics and the zero entropy subspace method. JP acknowledges Okinawa Institute of Science and Technology Graduate University and also the JSPS KAKENHI Grant Number 20K14417. We acknowledge financial support from the ANR project SuperRing (Grant No. ANR-15-CE30). LPL is a member of DIM SIRTEQ (Science et Ingénierie en RégionÎle-de-France pour les Technologies Quantiques).

Methods
Here we provide details on the different methods and approaches used in the main text.
Details on the solution of the Gross-Pitaevskii equation To describe the dynamics in the γ 1 regime we solve the Gross-Pitaevskii equation (2 -main text) numerically, using a spectral method relying on the discrete sine transform embedding the hard wall boundary conditions ψ(0) = ψ(L) = 0. We first use imaginary time propagation to find the groundstate in the box, then we quench the state at t = 0 and compute the subsequent dynamics. To ensure that the system is in the mean-field hydrodynamic regime we choose a sufficiently large non linear coefficient gN = γN 2 = 20000. We have checked that the transition between the single particle and meanfield regime occurs at gN ∼ 500.
Details on the solution of the generalized hydrodynamics equations The GHD equations are given in Eqs. (3)(4) of the main text, which we recall here: where n is the occupation function of the Lieb-Liniger quasi-particles, and the dressed velocity is given by: At first the GHD formalism seems incompatible with the box boundary conditions, because it relies on the local density approximation. One method to naturally include the effect of the hard-wall boundaries is to double the system size (from [0, L] to [−L, L]), impose periodic boundary conditions with period 2L, and use an anti-symmetric initial state: the right part z ≥ 0 (resp. left part z < 0) is quenched with a positive (resp. negative) velocity boost: wheren(k) is the equilibrium occupation function obtained from the equation of state. This approach is well adapted to GHD and exact at the level of the initial Lieb-Liniger Hamiltonian.
To integrate the GHD equations at zero temperature we use the zero entropy subspace method [31]. In this case it is sufficient to compute the evolution of the edges of the Fermi sea, that are located initially at k = ±k F . After the quench described by (9), the edges are shifted to ±k F ± k 0 . Furthermore the box boundary condition imposes that a quasi-particle arriving at the right boundary with quasi-momentum k > 0 is reflected at quasi-momentum −k (particles at k < 0 are already moving away from the boundary). A symmetric condition occurs at the left boundary. Therefore, immediately after the quench, the dynamics of the front moving to the left is fixed by the quasi-particles lying in k ∈ [−k F − k 0 , −k F + k 0 ], while the front moving to the right corresponds to quasi-particles in k ∈ [k F − k 0 , k F + k 0 ]. The broadening of the fronts is then explained by the fact that these quasi-particles move at different effective velocities: for example, the width of the front moving to the right will evolve as: . This simple explanation indicates that both fronts broaden within GHD, as is seen in the simulation. Therefore, GHD is not able to reproduce the microscopic details of the exact GP and TG results, while giving an accurate prediction for global observables, see figures S.1 and S.2 below.
We have also benchmarked our results using an independent integration scheme, based on the iterative method of [32], which also allows for finite temperature calculations. To summarize, the occupation function at time t + dt is obtained by solving the implicit equation: This is done by iterating this formula starting with the initial guess n(z, k, t + dt) = n(z, k, t). During this process, periodic boundary conditions are enforced on the interval [−L, L]. To proceed numerically, we use a discrete rectangular grid to store the values of n(z, k, t) at time t and rely on a cubic interpolation formula on this grid to evaluate equation (10). Once satisfactory convergence is obtained the same method is repeated to compute the next time step, until the desired final time is achieved. Details on the Tonks-Girardeau exact solution In the infinitely strongly repulsive limit, γ → ∞, we focus on the exact Tonks-Girardeau (TG) model [22]. In particular, we make use of the time-dependent Bose-Fermi mapping [22][23][24], where the many-body wavefunction Ψ T G is writtn in Eq. (5 -main text).
Our specific protocol is the following: we write the initial wavefunction as the ground state of a hard-wall box potential, constructed by the first N single-particle orbitals χ (z), which we then multiply by a phase profile, induced by the phase imprinting, obtaining the wavefunction ψ 0 (z) = e ik0z χ (z), which is used as starting point for the time evolution. The evolution is then calculated by projecting this state in the eigenbasis of the unperturbed system ψ (z, t) = ∞ χ |ψ 0 χ (z)e −i t/ and where is the -th single-particle eigenenergy [33,34].
The current of a TG gas at finite temperature is then readily obtained in terms of the evolved single-particle orbitals according to with f ( ) being the Fermi-Dirac distribution. In our specific quench setup, the current density after the phase imprinting reads with an amplitude of the excitations being given by The sound velocity of a TG gas is readily obtained from its equation of state. In order to compare it with the generalized hydrodynamics predictions, we have included the first order correction due to the boundary [17], such that at zero temperature it reads: .   predictions also at strong interactions. At lower number of particles we attribute the discrepancies to finite size effects, that are not captured within GHD. Our analysis shows that the study of the shock wave dynamics provides a very accurate test of the validity of the GHD equations.
Oscillations of the current at finite temperature and strong interactions In the strongly interacting regime, at finite temperatures, bosonic particles can be described using the Bose-Fermi mapping, in which particles populate the eigenstates of the system following the Fermi-Dirac distribution. When considering a quench into such Fermi sphere, the Hilbert space over which the quenched state projects increases, leading to more low energy excitations during the quench [35]. In Fig. S.3 we calculate the total current, Eq. (11), at different temperatures. Note that at temperatures lower than the Fermi temperature, the current oscillations are still visible and follow a few full oscillations, which shows the robustness of the universal features discussed in the main text. For temperatures of the order of the Fermi temperature, the damping increases dramatically and shock waves diffuse rapidly.
Oscillations of the current at finite temperature using GHD In order to include finite temperature effects in the GHD we use the Thermodynamic Bethe-Ansatz [36]. The initial equilibrium occupation function is: where the pseudo-energy k is the solution of: The box boundary conditions and quench protocol are implemented in the initial state as in the zero temperature case and we use the iterative integration algorithm explained above. Figure S.4 shows the decay of the current oscillations as temperature increases. The phenomenon reported in the paper is robust up to T ∼ 0.25 × 2 n 2 0 /(mk B ) where n 0 = N/L is the one-dimensional density.
Identification of the density dips as a soliton train is seen as small density dips associated to well defined "steps" in the phase profile. Therefore it seems that for our scenario the density oscillations associated to the shock front propagation are mainly due to fast grey solitons. As the solitons propagate with slightly different speeds (the shallower the faster) and bounce back on the hard wall boundaries, the phase profile can be complicated to interpret at later times where solitons propagates in both directions and overlap. We have checked that the number of generated solitons increases with the quench strength. A simple fitting function for the density profiles We analyse the output of the GPE and TG simulations by using a trial density profile: n(z, t) = n 0 + ∆n 2 tanh x − x 1 (t) ξ 1 (t) + tanh x − x 2 (t) ξ 2 (t) .
(16) As shown in Fig. S.6 this simple model allows to fit reasonably well the density profiles in the GP and TG regimes, at different times, up to t = L/c GP . In particular the shape of the density steps is well captured, despite the density oscillations (soliton train in GP and Friedel oscillations in TG). By analyzing the fit results as a function of time we confirm to very good accuracy that: (i) the two fronts propagate in opposite directions but at the same velocityẋ 1 = −ẋ 2 = c GP , corresponding to the speed of sound; (ii) ξ 2 remains constant while ξ 1 grows approximately linearly with time; and (iii) the parameters n 0 and ∆n remain constant. Observation (i) is not a surprise as it follows from the long wavelength analysis discussed in the main text. Observations (ii) and (iii) suggest that it may be possible to find a simple analytical model describing the class of hydrodynamic shock waves studied in this work and suggest a common mechanism is at play throughout the whole interaction strength range. This analysis goes beyond the scopes of the current work.