Dispersive Shock Waves in a 1D Quantum Liquid

Abstract

We implement a moving boundary condition in a 1D quantum liquid to study nonlinear wave breaking and its regularization by dispersion. Programmable optical potentials allow us to compress a weakly interacting ultra-cold Bose gas trapped on an atomchip at tunable speeds of up to three times the speed of sound and subsequently measure the quasi-in-situ density distribution to extract the shock wave edge dynamics. We resolve both leading and trailing edge velocities and observe a shock wave width that increases linearly in time, which are distinguishing features of dispersive shock waves, consistent with asymptotic predictions using Whitham’s method. Quantitative agreement is found with finite temperature non-polynomial Schrödinger equation simulations, taking into account the imaging process and the finite height of the piston potential. Our results constitute a controlled, quantitative test of dispersive shock dynamics in a 1D quantum fluid and demonstrate that the coarse-grained dispersive-shock phenomenology remains robust even as the microscopic dynamics depart from the strictly integrable 1D regime.

Introduction $-$ In nonlinear media steepening of wave fronts leads to a so-called gradient catastrophe, requiring a mechanism to regularize the ensuing wave breaking. Classical fluids exhibit viscosity, resulting in the formation of viscous shocks and entropy production, whereas in dissipationless media dispersion yields an expanding rank-ordered train of oscillations, namely a dispersive shock wave (DSW) 53; 24; 23; 52; 17; 27.

Ultracold dilute Bose gases present a powerful platform for investigating quantum many-body dynamics 5; 11; 4; 38. Confined to 1D, they become close to integrable 8 and effectively dissipationless over experimental timescales. This makes them an ideal testbed for dispersive phenomena, while offering precise experimental control. Following recent progress in the dynamics of 1D systems, investigation of higher-order hydrodynamic effects has attracted a lot of interest 9; 1; 13; 2; 46; 7; 47; 14; 37; 35; 50; 45. The ensuing rich nonequilibrium dynamics at the wave breaking offer insights beyond the ballistic (Euler) regime. To that end, the canonical piston protocol constitutes a direct, controlled, and quantitative test of dispersive shock dynamics.

Dispersive shock waves have been reported in nonlinear optics for an initial discontinuity 41; 19; 54 and during four-wave mixing 18. However, in systems of ultracold Bose gases, investigations of shock dynamics have been limited to 3D 16; 48; 26, albeit some highly anisotropic configurations 33; 12; 34. Thus, the ensuing dynamics exhibited dissipation in the form of formation of topological defects. While previous works present dynamics consistent with DSWs 26; 33; 12, they do not resolve both edge velocities, the shock width growth or the behavior at the critical velocity.

Here we implement the canonical piston problem in a strongly confined quasi-1D quantum liquid using a programmable optical wall and probe the shock dynamics for a wide range of piston speeds. We analyze coarse-grained density profiles after time-of-flight (TOF) to directly measure both shock edge velocities and the shock region width, a distinguishing feature of DSWs 53; 28; 52; 17; 29. We subject our effectively dissipationless system to large perturbations, pushing local densities beyond the strict 1D-condition 38; 11; 42, testing dispersive nonlinear wave dynamics. Further, our work constitutes a quantitative test of Whitham’s method 53; 51; 17; 52, which provides an asymptotic solution to the dispersive shock wave in the piston problem.

Experimental Setup $-$ Our experimental system is a one-dimensional ultracold Bose gas of $4-10\times 10^{3}$ ⁸⁷Rb atoms trapped below an Atom Chip 40. The chip produces a cigar-shaped magnetic trap with $\omega_{\perp}=2\pi\times 1.38~\mathrm{kHz}$ along the two transverse directions, satisfying the 1D condition $\hbar\omega_{\perp}>k_{B}T,\mu$, where $\mu$ is the chemical potential, $k_{B}$ Boltzmann’s constant, and $T$ the temperature, which is typically in the range $20-90$ nK. The interaction strength is quantified by the dimensionless Lieb-Liniger parameter $\gamma=mg_{1\mathrm{D}}/(\hbar^{2}n_{\text{1D}})$ 31. In our experiment, $\gamma\approx 0.002$, which, combined with the low temperature, places the gas deeply in the quasi-condensate regime. To implement arbitrary dynamic longitudinal potentials, we utilize a Digital Micromirror Device (DMD) combined with a blue detuned laser ($\lambda=766.7~\mathrm{nm}$) to project pre-engineered intensity patterns onto the gas along the weak longitudinal confinement in $z$-direction 49, see Fig. 1(a).

Figure 1: Illustration of the experimental protocol and measured density carpet. (a): Ultracold ⁸⁷Rb atoms are trapped in a 1D confinement using an atomchip. Superimposing (blue detuned) optical dipole potentials created by a DMD enables the creation of arbitrary dynamic longitudinal potentials, including (b): the realization of a moving piston by successively turning on columns of pixels. The averaged density profile $n(z,t)$ is shown in red (shaded: standard error), with extracted shock edge positions $z_{\pm}$. (c): Density carpet measured during a piston protocol: $n(z,t)$ as a function of position and time. A piston with constant velocity $v_{\mathrm{p}}$ compresses the condensate from the left. The extracted shock edge positions $z_{\pm}$ are shown as red dots and the corresponding linear fits $v_{\pm}$ are shown as solid blue/red lines with confidence intervals as shaded regions. Theoretical predictions using Whitham’s method are shown as solid lines. Measurement time is limited by the finite system size.

The piston $-$ To implement the piston protocol we prepare a 1D thermal state directly into a box potential. At $t=0$ we start to successively turn on columns of DMD pixels next to the wall, which locally increases the light intensity, moves the wall, and pushes the atoms inward at a constant velocity, thus implementing the piston. In the resulting dynamics, atoms accumulate at the edge of the piston and a shock region forms, which is characterized by a rapidly oscillating density profile. In classical hydrodynamics, a balance of dissipation and nonlinearity leads to a viscous monotonic shock region of constant width 52. In contrast, in the absence of dissipation, dispersion gradually fans out the oscillations depending on their amplitude, resulting in an oscillating growing shock region 52. Towards the leading edge, short wavelength small-amplitude linear oscillations are predicted on the unperturbed background, while the trailing edge matches a (large-amplitude) soliton on the compressed region 52. The experimental protocol is illustrated in Fig. 1.

Measurement $-$ We measure the quasi-in-situ linear density profile $n(z,t)$ via absorption imaging after 2 ms TOF expansion, repeating the experiment for several evolution times $t$ and piston speeds $v_{\mathrm{p}}$ for 50-80 repetitions. We determine the leading and trailing edges of the shock region for intermediate evolution times $t$, that is, after the shock region has expanded to an observable size, but before the train of waves created by the piston has reached the opposite wall of the box. An example profile is shown in Fig. 1(b). We identify an unperturbed region matching the initial density distribution far from the action of the piston. The leading edge $z_{+}$ of the shock region is found where a fitted error-function of the relative density $n(z,t)-n(z,0)$ exceeds the standard error of the density distribution (details see Supplementary Material 44). The trailing edge $z_{-}$ is found at the edge of the error function, except close to the critical velocity $1.4v_{\mathrm{s}}<v_{\mathrm{p}}<2.2v_{\mathrm{s}}$, where $z_{-}$ coincides with the position of the piston. Notably, we find that the wave fronts $z_{\pm}$ propagate with constant velocity, see Fig. 1(c), enabling the straightforward extraction of the corresponding DSW velocities $v_{\pm}$ via a linear fit.

Framework $-$ The piston problem is typically studied using Whitham’s method applied to the zero temperature Gross-Pitaevskii equation (GPE). This yields an asymptotic description of the coarse-grained shock wave envelope in the form of a traveling wave, from which the oscillations can be reconstructed (see Supplementary Material44) 27; 17; 29; 25; 24; 52; 53. Because the shock wave microstructure lies below our imaging resolution (healing length $\xi_{\mathrm{h}}\approx 0.26~\mu\mathrm{m}$, imaging point-spread-function width $\sigma_{\mathrm{PSF}}\approx 3~\mu\mathrm{m}$), the highly oscillatory wave train of the DSW cannot be imaged directly. Additionally, finite temperature effects and experimental fluctuations decrease the contrast of the DSW, especially towards the leading edge 46; 37. Furthermore, free evolution during TOF expansion mixes phase and density (see Supplementary Material 44): large phase gradients at the trailing edge of the DSW lead to an increase in observed atomic density 37, whereas the leading edge is minimally affected, as phase gradients are small. These experimental constraints and factors are not included in Whitham’s treatment of the piston. Therefore, we compare experimental data to finite-temperature simulations of the quasi-1D non-polynomial Schrödinger equation (NPSE) 42, which is the conventional GPE including transverse broadening. Furthermore, we account for the measurement process. The NPSE, its parameters, and the simulation procedure are detailed in the End Matter.

Results $-$

Figure 2: Shock region edge speeds for the piston
protocol. The experimentally measured speed for the leading edge is shown in blue, trailing edge in orange. The corresponding theory values are shown as solid lines in blue and orange respectively. The blue shaded region is obtained via finite temperature NPSE simulations including TOF and imaging. The piston speed $v_{\mathrm{p}}$ is shown as a dashed black line. Experimental standard errors are derived from statistics of multiple measurements as well as errors extracted from the fitting procedure and are in the range of: $\bar{\sigma}_{v_{\mathrm{p}}/v_{s}}\approx 0.01-0.07~v_{\mathrm{s}}$, $\bar{\sigma}_{v_{\mathrm{\pm}}/v_{s}}\approx 0.1-0.4~v_{\mathrm{s}}$. Inset: Shock wave widths, $W(t)=z_{+}(t)-z_{-}(t)$ and linear fits as solid lines. The widths are colored in order of increasing $v_{\mathrm{p}}$ as : red, green, purple, yellow, black, gray

First, we extract and study the propagation velocities of the shock wave fronts across a wide range of piston velocities. The results are presented in Fig. 2. The measured velocities are consistent with the theoretical predictions of Whitham’s method for dispersive shock waves (see End Matter Eq. 3): The velocity of the leading edge $v_{+}$ asymptotically approaches $2v_{\mathrm{p}}+2v_{\mathrm{s}}$. Meanwhile, the trailing edge velocity $v_{-}$ experiences a kink, a discontinuity in the first derivative, at the critical velocity $v_{\mathrm{p}}=2v_{\mathrm{s}}$. Below, it scales as $v_{\mathrm{p}}/2$ and above $v_{-}$ approaches a linear growth scaling as $2v_{\mathrm{p}}$. Crucially, the observed kink and widening shock region are both distinguishing properties of DSWs. The measured leading edge speeds $v_{+}$ lie systematically below the theoretically predicted values, which we attribute to finite temperature and imaging effects. Comparison with simulations reveals a similar behavior. Error bars are derived from standard errors of the density distribution, as well as errors from the fitting process.

Next, we study the accumulation of atoms between the piston and the shock region. Whitham’s method predicts a region of constant density for $v_{\mathrm{p}}<2v_{\mathrm{s}}$, whereby the plateau density $n_{\mathrm{plat}}$ increases quadratically until the critical piston speed $v_{\mathrm{p}}=2v_{\mathrm{s}}$, where the piston catches the trailing edge. Beyond, $v_{\mathrm{p}}=2v_{\mathrm{s}}$ the plateau is replaced by large density oscillations with a maximum amplitude of $n_{\mathrm{plat}}=4n_{0}$ independent of $v_{\mathrm{p}}$ (see End Matter Eq. (4)). The results are shown in Fig. 3.

Figure 3: Maximum atom accumulations for the piston protocol. Measured plateau densities (black) are compared to Whitham theory predictions (dark blue line) and simulations of increasing complexity: zero temperature GPE in-situ (blue pluses), finite temperature NPSE in-situ (light green region) and finite temperature NPSE simulations including 2 ms TOF and imaging (NPSE TOF imaged - dark green region). Experimental standard errors are in the range $\approx 0.02-0.06~n_{0}$.

For slow piston speeds, $v_{\mathrm{p}}<1.2v_{\mathrm{s}}$, the measured $n_{\mathrm{plat}}$ agrees closely with the predictions of Whitham theory. At larger $v_{\mathrm{p}}$, however, the observations lie considerably below the predicted values. We assess this discrepancy using simulations of progressively increasing complexity. First, finite-temperature in-situ NPSE simulations account for thermal effects and transverse broadening, isolating deviations arising from the quasi-1D dynamics. The reduced speed of sound leads to a slightly higher plateau density for $v_{\mathrm{p}}<2v_{\mathrm{s}}$. We then include the measurement process by accounting for 2 ms TOF expansion and the finite imaging resolution, thereby incorporating phase–density mixing and the coarse-graining of large density oscillations (more details in the Supplementary Material 44). This reproduces the observed reduction in maximum density and captures the experimental trend well. At the highest piston speeds, $v_{\mathrm{p}}>1.6v_{\mathrm{s}}$, residual deviations in both the edge velocities and maximum densities remain, which we trace to the finite piston height in the next paragraph.

Finite piston $-$ To better understand the experimental dynamics we need to accurately model the implemented piston. Unlike previous works that sweep a laser beam across the condensate 34, we shape the full beam with the DMD and overlap it with the cloud, which reduces the intensity, see Fig. 1(a) (more details in the End Matter). For high piston speeds $v_{\mathrm{p}}\geq 1.6v_{\mathrm{s}}$, we observe that the local chemical potential starts to exceed the finite dipole potential of the piston, whereby excess atoms are lost from the box. Importantly, these are local spillovers into the larger harmonic confinement, with the atoms still detectable, and thus differ fundamentally from genuine losses in their effect on 1D gases 6. We call this piston leaky.

Figure 4: Edge velocities for a leaky piston. Experimentally measured leading and trailing edge shock region speeds $v_{\pm}$ (blue/orange markers) are compared to theory predictions (solid blue/orange lines) and finite temperature NPSE simulations (shaded area $v_{+}$), including TOF and imaging, with varying piston heights $U_{\mathrm{p}}$ in relation to the initial chemical potential $\mu$. The shading is intended as a guide to the eye. Standard errors of $v_{+}/v_{\mathrm{s}}$ are in the range $0.05-0.1$. Inset: For the same datasets: $\gamma_{p}$: Percent of total atoms that leave the box per time as a function of $v_{\mathrm{p}}$. Experimental errors are obtained from standard errors of the density distribution as well as errors extracted from the fit. Standard errors of $\gamma_{\mathrm{p}}$ are in the range $0.05-0.30\%$.

The measured velocities of the shock wave edges for a leaky piston are presented in Fig. 4. We compare datasets with the same approximate initial system size $L\approx 70~\mu\mathrm{m}$ and density $n_{\mathrm{1D}}\approx 70~\mu\mathrm{m}^{-1}$. These datasets are driven by the same finite piston height (for more details, see the Supplementary Material 44), and thus by fixing the ratio of initial chemical potential to $U_{\mathrm{p}}$ we can observe the dynamics for varying $v_{\mathrm{p}}$. We find that both shock edge velocities, $v_{-}$ and $v_{+}$ exhibit a saturation at $v_{\mathrm{p}}\approx 2.3v_{\mathrm{s}}$. The effect of the finite piston height is most evident at the leading edge of the shock wave, since the highest-momentum particles that would form it instead leave the box.

To confirm that the observed dynamics are caused by a finite height piston, we compare with NPSE simulations for varying strengths of the piston potential. Here we find a good agreement. The piston loss rate: $\gamma_{\mathrm{p}}=(1/N)~\mathrm{d}N_{\mathrm{p}}/\mathrm{d}t$, where $N_{\mathrm{p}}$ is the number of atoms measured outside of the box, increases approximately linearly with $v_{\mathrm{p}}$ from $\approx 2.1~\%/$ms to $6~\%/$ms, similar to simulations, albeit systematically higher. We attribute the deviations to systematic imperfections and non-uniformity of the piston potential (for more details, see Supplementary Material 44). The genuine atomic loss rates are very low at $\leq 0.25\%$/ms and do not depend on $v_{\mathrm{p}}$. For very high $v_{\mathrm{p}}>3v_{\mathrm{s}}$, the measured $v_{+}$ and $v_{-}$ collapse to the piston speed and we observe a very narrow, almost constant shock region propagating rigidly with the piston, in agreement with simulations (for more details, see the Supplementary Material 44).

Transverse effects $-$ For an ideal piston, increasing $v_{\mathrm{p}}$, raises $\mu\gtrsim\hbar\omega_{\perp}$, populates transverse modes and eventually leads to topological excitations, producing dissipative-like dynamics as observed in prior piston experiments with lower transverse confinement 34. Notably, we do not observe topological excitations or dissipative dynamics. For $v_{\mathrm{p}}>2v_{\mathrm{s}}$ the measured peak densities yield $\mu/\hbar\omega_{\perp}$ in the range of $1.13-1.39$, while in-situ estimates are higher (simulations predict values in the range of $2.17-2.85$), indicating non-negligible transverse effects (see Supplementary Material 44). Nevertheless, the longitudinal dynamics remain well-reproduced by quasi-1D NPSE, explicitly including transverse broadening, which underscores the importance of transverse effects 36; 10; 22.

In brief, the coarse-grained DSW phenomenology remains observable for three reasons. First, the finite piston height limits the local atom accumulation at high piston speeds, as excess atoms simply exit the box. Second, the over-compressed zone is spatially localized next to the piston and the characteristics that form the leading edge quickly exit this region. Third, the observed time-scales are not long enough for the injected transverse energy, visible as transverse broadening (see Supplementary Material 44), to develop into collective 3D dynamics or topological excitations. Especially for large $v_{\mathrm{p}}$, finite size effects limit the measurement time considerably.

Conclusion $-$ We presented a systematic experimental study of the piston problem in a 1D Bose gas. Our observations demonstrate the presence of dispersive shock waves for piston speeds up to three times the speed of sound, authenticated by measuring linearly growing shock widths and two distinct shock region speeds that follow the predictions obtained by Whitham’s method. Notably, the trailing edge speed displays a kink at the critical piston velocity $v_{\mathrm{p}}=2v_{\mathrm{s}}$. Finite-temperature NPSE simulations that include TOF and imaging effects reproduce the measured shock dynamics and peak densities, while behavior at high piston speeds are consistent with a finite-height potential. Remarkably, despite local densities exceeding the strict 1D regime, the coarse-grained DSW phenomenology remains observable, highlighting that the dispersive-shock description survives even as the microscopic dynamics depart from the strictly integrable 1D regime. Our results constitute a controlled, quantitative test of dispersive effects in a quasi-1D quantum fluid.

This letter opens further experimental avenues: Increased optical power would allow to approach the ideal piston, enabling investigation into the expected breakdown of integrability at large piston speeds. Sub-micron optical resolution and interferometric phase readout could reveal the shock microstructure, further progressing the investigation into Whitham’s asymptotic predictions and the behavior of dispersive shock waves 21. Multi-piston protocols enable the study of multiple shock regions and shock-shock interactions 3. The shock wave can be characterized by measuring the local quasi-momentum distribution 15. The introduction of weak disorder combined with entropy measurements can provide insights into prethermalization and relaxation dynamics of nearly integrable quantum systems 20; 30.

Acknowledgments $-$ The research at TU Wien was supported by the European Research Council under Horizon Europe (project EmQ, Grant No. 101097858), and by the Austrian Science Fund (FWF) (Grant DOI: 10.55776/P36656) and the Austrian Science Fund (FWF) through the Cluster of Excellence ‘Quantum Science Austria’ (quantA) (Grant DOI: 10.55776/COE1).

F.M. and P.S. conceived the idea. P.S. carried out the experiments with contributions by F.M., M.T., F.C., S.-C.J. and N.B.. P.S. performed the data analysis, and simulations. P.S. wrote the paper with contributions from all authors. F.M. and J.S. supervised the project.

Data availability $-$ The data supporting the findings of this study are openly available in the Zenodo repository [DOI 10.5281/zenodo.21991533] 43.

End Matter

Technical details of the experimental procedure $-$ We use a sequence of standard cooling and trapping techniques to prepare $N=4-10\times 10^{3}$ atoms at a temperature of $T=20-90$ nK. Temperature is measured via density ripples thermometry, utilizing the mixing of phase and density in TOF, whereby thermal in-situ phase fluctuations translate into density fluctuations (ripples) 32. Crucially, we cool the condensate directly into the box trap to prepare initially thermal states. At $t=0$ ms we start the protocol and after the desired evolution time the condensate expands freely for $2$ ms TOF after which we measure the transverse 2D atomic density distribution via absorption imaging, thereby destroying the system. We utilize an imaging system with magnification $12.39$, pixel size in object space: $1.05~\mu\mathrm{m}$ and an effective numerical aperture of NA$=0.2$. The estimated point spread function obtained from density ripples thermometry after $12.2$ ms TOF is $\sigma_{\mathrm{PSF}}\approx 3~\mu\mathrm{m}$. We obtain $50-80$ measurements per evolution time and finally post-select on total atom number, fitting a normal distribution and retaining shots within a $\pm 1\sigma$ window ($\pm 0.7\sigma$ for the leaky piston) around the mean and discarding approximately $30\%$ ($50\%$) of shots. Although the atom number determines the local speed of sound and in turn the shock dynamics, we observe no cross-dataset correlation of shock observables and post-selection window.

Creating arbitrary longitudinal potentials $-$ More details regarding the creation of arbitrary longitudinal potentials are found in Ref. 49. Here we summarize the setup for a better understanding of the resulting finite piston power. We utilize a $\lambda_{\mathrm{dip}}=766.8~\mathrm{nm}$ blue detuned (with respect to the D2 transition of ⁸⁷Rb at $\lambda_{\mathrm{D2}}=780.24~\mathrm{nm}$) Gaussian beam whose horizontal extension is elongated by a cylindrical telescope to illuminate the DMD, a SuperSpeed V-9501 Module from ViALUX, featuring 1080 by 1920 micromirrors each sized $10.8~\mu\mathrm{m}$. The light which is reflected by the DMD in the ”ON” position is superimposed on the transverse imaging path of the experimental setup via a polarizing beam splitter cube. The horizontal x-axis of the DMD matches the longitudinal z-axis of the condensate. Two demagnification stages, including a lens on a motorized stage, are used to image and focus the DMD pattern in the plane of the atoms, resulting in a final pixel size in the plane of the atoms of $0.45~\mu\mathrm{m}$. As the transverse harmonic oscillator ground state width in the trap is $l_{\mathrm{ho}}\approx 0.3~\mu\mathrm{m}$, smaller than the DMD pixel size, we use a horizontal slit in the Fourier plane of the first demagnification stage to implement grayscaling via spatial filtering, which reduces the total intensity. We optimize DMD patterns via iteratively comparing the in-situ atomic density to the desired potential and updating pixel states until convergence. To offset the weak longitudinal harmonic confinement, the box pattern displays the inverse profile with high pixel density at the center decreasing towards the edges of the box. The walls are implemented using the maximally available intensity, that is columns with all pixels on. The piston is then implemented by successively switching on columns of pixels atop this box, starting at the wall, resulting in a reduced effective power towards the center of the condensate, where many pixels are already active. Overall, the piston power is reduced by grayscaling and exhibits longitudinal non-uniformity due to the elongated Gaussian beam and the shared use of DMD for both box and piston.

Theoretical description (NPSE) $-$ The dynamics of a quasi-1D gas of $N$ bosons with repulsive contact interactions is well described by the non-polynomial Schrödinger equation (NPSE) 42, which differs from the conventional Gross-Pitaevskii equation (GPE) 39 by including transverse Gaussian broadening of the wave-function:

$$i\hbar\partial_{t}\Psi=\left[-\frac{\hbar^{2}}{2m}\partial_{z}^{2}+V+\hbar\omega_{\perp}\sqrt{1+2a_{\mathrm{s}}|\Psi|^{2}}\right]\Psi\text{.}$$ (1)

where $\Psi(z,t)$ is the order parameter, $V(z,t)$ the external longitudinal potential (we leave out the $(z,t)$-dependence for readability), $\hbar$ the reduced Planck constant, $m$ the atomic mass, $\omega_{\perp}$ the transverse trapping frequency and $a_{\mathrm{s}}$ the scattering length. Writing the GPE in the Madelung representation (More Details see Supplementary Material 44) reveals a term $\propto\partial_{z}^{2}|\Psi|$, which is often called quantum pressure and becomes relevant on length scales smaller than the healing length $\xi_{\mathrm{h}}=\hbar/\sqrt{2mn_{\text{1D}}g_{\text{1D}}}$, regularizing the solution at the shock. The effective 1D interaction strength is denoted by $g_{\mathrm{1D}}=2\hbar\omega_{\perp}a_{\mathrm{s}}/\sqrt{1+2a_{\mathrm{s}}n_{\mathrm{1D}}}$. The homogeneous 1D (background) density $n_{\text{1D}}$ determines the speed of sound:

$$v_{\mathrm{s}}=\sqrt{\frac{\hbar\omega_{\perp}a_{\mathrm{s}}n_{\text{1D}}(2+3a_{\mathrm{s}}n_{\text{1D}})}{m(1+2a_{\mathrm{s}}n_{\text{1D}})^{3/2}}}\;.$$ (2)

Note that both $g_{\mathrm{1D}}$ and $v_{\mathrm{s}}$ are slightly lower compared to the conventional GPE definitions due to the effects of the transverse broadening.

Whitham DSW observables $-$ Whitham’s Modulation Theory applied to zero temperature GPE for the piston protocol results in a growing shock region where the shock edge velocities can be found by enforcing boundary conditions with the dispersionless region 53; 52; 27; 17:

$$v_{+} =\frac{2v_{\mathrm{p}}^{2}+4v_{\mathrm{p}}v_{\mathrm{s}}+v_{\mathrm{s}}^{2}}{v_{\mathrm{p}}+v_{\mathrm{s}}}\;,$$ (3)

$$v_{-} =0.5v_{\mathrm{p}}+v_{\mathrm{s}}~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\text{for}~v_{\mathrm{p}}<2v_{\mathrm{s}}\;$$

$$v_{-} =v_{\mathrm{p}}+\frac{2v_{\mathrm{s}}(v_{\mathrm{p}}-2v_{\mathrm{s}})K(M)}{v_{\mathrm{p}}E(M)-(v_{\mathrm{p}}-2v_{\mathrm{s}})K(M)}~~~\text{for}~v_{\mathrm{p}}\geq 2v_{\mathrm{s}}\;,$$

where K and E are the complete elliptic integrals of first and second kind respectively, of $M=4v_{\mathrm{s}}^{2}/v_{\mathrm{p}}^{2}$. The plateau density follows

$$n_{\text{plat}} =\left(\frac{v_{\mathrm{p}}}{2}\sqrt{\frac{m}{g_{\text{1D}}}}+\sqrt{n_{0}}\right)^{2}~~~~~~~~\text{for}~v_{\mathrm{p}}<2v_{\mathrm{s}}\;$$ (4)

$$n_{\text{plat}} =4n_{0}~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\text{for}~v_{\mathrm{p}}\geq 2v_{\mathrm{s}}\;.$$

Numerical simulations $-$ We use a standard Fourier-split-step method to simulate the NPSE and GPE dynamics. First, ground states are obtained via imaginary time propagation and used as initial guesses for further stochastic GPE/NPSE: We evolve the conventional GPE/NPSE with an additional noise term that simulates the system being in contact with a thermal bath, whose correlations are related to the final temperature, arriving at thermal initial states. The box potential is modeled via two finite-sized error function walls on both sides. The piston height is chosen at five times the chemical potential to exclude atoms leaking over the piston, except in the cases explicitly testing the effects of the finite piston height. The piston is implemented by moving one of the walls at constant speed starting at $t=0$. Finally we obtain statistics from 100 realizations. In order to account for TOF and the imaging process we allow for a free propagation and afterwards convolve with a Gaussian, simulating a point spread function of $\sigma_{\mathrm{PSF}}=3~\mu\mathrm{m}$ consistent with our density ripples measurements. Notably, our simulations do not account for systematic box imperfections, experimental drifts of atom number and temperature, or for longitudinally varying piston power.