Dissipative shock waves generated by a quantum-mechanical piston

The piston shock problem is a prototypical example of strongly nonlinear fluid flow that enables the experimental exploration of fluid dynamics in extreme regimes. Here we investigate this problem for a nominally dissipationless, superfluid Bose-Einstein condensate and observe rich dynamics including the formation of a plateau region, a non-expanding shock front, and rarefaction waves. Many aspects of the observed dynamics follow predictions of classical dissipative—rather than superfluid dispersive—shock theory. The emergence of dissipative-like dynamics is attributed to the decay of large amplitude excitations at the shock front into turbulent vortex excitations, which allow us to invoke an eddy viscosity hypothesis. Our experimental observations are accompanied by numerical simulations of the mean-field, Gross-Pitaevskii equation that exhibit quantitative agreement with no fitting parameters. This work provides an avenue for the investigation of quantum shock waves and turbulence in channel geometries, which are currently the focus of intense research efforts.

F rom the generation of localized solitons and quantized vortices 1 to the extended coherence of dispersive shock waves 2 , quantum hydrodynamics exhibit an intriguingly rich phenomenology. While many pioneering observations have been made in superfluid helium 3,4 , dilute-gas Bose-Einstein condensates (BECs) provide an exceptionally versatile medium in which to access quantum hydrodynamics 5 . The experimental control and theoretical understanding of BECs enable novel techniques for entering new quantum hydrodynamic regimes. A central focus of superfluid helium studies has been the investigation of quantum turbulence, including the origin of dissipation within the system 6 . Despite strong experimental and theoretical research efforts spanning many decades (see ref. 7 and references therein), quantum turbulence still poses many open questions. For example, while many aspects of quantum turbulence in a homogeneous system have been clarified, the nature of quantum turbulence in channel geometries is now under intense investigation in superfluid helium systems. In dilute-gas BECs, quantum turbulence has been experimentally observed only in a very limited number of settings so far. Those include the generation of vortex turbulence in stirred BECs [8][9][10] , the observation of weakwave turbulence in a shaken BEC confined in a box potential 11 , and the observation of spinor turbulence 12,13 . Aspects involving the definition of a superfluid Reynolds number 14 , or energy and enstrophy cascades [15][16][17] , remain under very active theoretical investigation. A discussion of relevant experimental realizations but also of theoretical attempts to study the problem has recently been compiled in ref. 18 .
Here, we introduce a setting for the observation of rich quantum hydrodynamics by studying a BEC piston shock in a channel geometry. The piston shock is a paradigmatic example of strongly nonlinear flow, probing hydrodynamics in an extreme regime. For a one-dimensional (1D) channel, the BEC piston shock is theoretically predicted to be an expanding, coherent dispersive shock wave (DSW) with rank-ordered, nonlinear oscillations 19 . A related setting, the collision of two BECs with strong transverse confinement in a channel, has been shown experimentally and theoretically to give rise to similar dynamics, such as the continuous transformation of a sinusoidal interference pattern into a train of dark solitons that was interpreted as two adjacent dispersive shock waves 20 . However, in the presence of weaker transverse confinement, collision experiments in BECs 21 (and also in dilute Fermi gases 22 ) cannot be described by DSWs but by two counterpropagating, viscous or dissipative shock waves (VSWs) 23 . Two features that distinguish a VSWhowever weak the dissipation-from a DSW are: (i) its shock width is independent of time and proportional to the medium's dissipation, (ii) its speed is uniquely determined by the Rankine-Hugoniot jump relations 24 . In contrast, a DSW exhibits an expanding series of rank-ordered oscillations with two edge speeds, each of which satisfies DSW closure relations that are entirely different from the Rankine-Hugoniot relations 2 .
The fully three-dimensional (3D) BEC piston shock problem explored here provides a clean setting for the quantitative study of the roles of dispersion and dissipation in nonlinear quantum hydrodynamics. We observe the generation of non-expanding, large-scale shocks that satisfy the Rankine-Hugoniot jump relations, image features indicative of vortex turbulence at small, healing length scales, and perform numerical simulations of the conservative, mean-field Gross-Pitaevskii (GP) equation that quantitatively agree with experiment. Piston compression is found to continually generate two distinct wave-field components: sound waves and solitons. The solitons breakup into vortices via the well-known snake instability. We argue that the decay of large amplitude excitations generated at the shock front into vortex excitations lead to an emergent dissipative-like behavior in a coarse-grained description of the fluid. Consequently, piston compression provides both the generation mechanism for turbulence into which large amplitude soliton like excitations dissipate and the sustenance of a subsonic to supersonic shock front. Thus, the piston shock problem also opens a pathway to the study of quantum turbulence with BECs in channel geometries.

Results
Overview of shock dynamics. The general setting for our study is comprised of an elongated BEC confined by a cigar-shaped harmonic trap. A repulsive barrier, the height of which exceeds the chemical potential of the BEC by a factor of 10, is created by the dipole force of a far detuned laser beam. Initially, the barrier is placed outside the BEC and is then moved through the BEC at a constant speed denoted v p . Our axis convention is such that the weakly confined z-axis is oriented horizontally [ Fig. 1(a)]. For reference, the bulk speed of sound in the center of the initial, unperturbed BEC is c s,bulk ≈ 2.47 mm s −1 , and the speed with which sound pulses travel along the long axis of the initial, unperturbed BEC is calculated to be c s;1d ¼ c s;bulk = ffiffi ffi 2 p % 1:75 mm s À1 in the Thomas-Fermi regime 25 . Absorption images are taken at sequential times during the piston sweep to analyze the resulting dynamics.
For sweep speeds near or exceeding the speed of sound, c s,1d , such as the case shown in Fig. 1 (b), the dynamics are intriguingly rich. At t = 0 ms, the piston is located just to the right side of the BEC. As the piston enters the BEC, a pronounced density spike forms near the piston front. As the piston sweeps through the cloud, a plateau region of high density develops in front of the piston ( Fig. 1(b), 60 ms). The leading (left) edge of the plateau forms a steep shock front, where the density rapidly drops from the plateau density to the initial BEC density. Once the plateau reaches the left wing of the BEC, the shock front deforms and approaches the shape of a rarefaction wave (see discussion below). The onset of this behavior can be seen in Fig. 1(b) at t = 140 ms.
The experimentally observed dynamics are in excellent agreement with our 3D numerical simulations based on the Gross-Pitaevskii equation (GPE) with no fitting parameters. Parameters utilized in the numerics are taken directly from experiment. For more details, see Methods.
Characterization of the shock dynamics. For a quantitative analysis of the piston shock wave dynamics, we consider three characteristic features: the shock propagation speed, the shock width, and the plateau density. Their behavior as a function of the piston speed is shown in Fig. 2. The shock propagation speed is determined by tracking the position of the shock front as it moves through the central part of the BEC where it exhibits an approximately constant speed. For both experimental and numerical data, the shock front edge is obtained by calculating integrated cross sections and subtracting the cross section of an unperturbed BEC in the absence of a piston sweep. In the subtracted plots, the shock front is fit with a straight line. The zero intercept of the fit is recorded as the position of the shock front, and its speed is the shock propagation speed. Experiment and corresponding numerics indicate an approximately linear increase in the shock speed with increasing piston speed (Fig. 2a). The inset of Fig. 2a shows the evolution of the shock front width during sweeps with v p = 2 mm s −1 (blue solid) and v p = 3 mm s −1 (red dashed). The width is determined by following the shock front as it is driven through the BEC. Similar to finding the shock front position, we calculate integrated cross sections, subtract an unperturbed cross section and record the spacing between the zero intercept and the edge of the plateau region. Both piston speeds result in a time-averaged width of Δx w = 15.5 μm. This constant width over the course of a sweep is indicative of VSWs. The plateau density, as a function of piston speed, is determined by averaging the density of the plateau region that has formed in front of the barrier when the piston reaches the center of the BEC. This is where the initial density of the BEC is nominally uniform. The observed dependence on the piston speed is approximately linear for medium piston speeds, but undershoots a linear trend at large piston speeds (Fig. 2b). In all cases, we see remarkable agreement between experiment and numerics based on the GPE.
Because dilute-gas BECs are typically modeled as inviscid, dispersive superfluids, an obvious approach is to compare the observed behavior with that of BEC dispersive shock wave (DSW) theory [26][27][28] (see also the review ref. 2 ). The 1D DSW theory for the BEC piston problem was presented in ref. 19 (see also ref. 29 ). A dispersive shock wave is characterized by an expanding, coherent nonlinear waveform with a trailing large amplitude soliton edge, rank-ordered interior oscillations, and a leading small amplitude, harmonically oscillatory edge. However, the best fit to the predicted DSW soliton edge speed and plateau height is poor (see Supplementary Fig. 1). On the other hand, VSW theory, where a narrow, planar shock front is assumed to propagate through a uniform, weakly viscous medium (see, for example ref. 30 ), leads to a consistent description of both the shock speed data and the plateau height data if an effective speed of sound of c s,eff = 1.35 mm s −1 is assumed in the calculation. The resulting VSW predictions are shown as the red solid lines in Fig. 2. VSW speed and plateau height are independent of viscosity, and are determined by the flow conditions in front of and behind the shock front according to the Rankine-Hugoniot jump conditions 30 . Relevant details of this analysis are given in Supplementary Note 1. We note that the obtained effective speed of sound is lower than the bulk speed of sound in the BEC center (c s,bulk ≈ 2.47 mm s −1 ) and the speed of sound for effectively 1D waves along the BEC's long axis (c s,1d ≈ 1.75 mm s −1 ). Utilizing numerical simulations below, we argue that this difference is due to a decrease in the effective hydrodynamic density and hence pressure jump across the shock front.
Numerical simulations and vortex turbulence. The applicability of viscous shock theory may seem surprising for a nominally inviscid superfluid. Further insight can be gained from our comparative (3 + 1)D numerical simulations of the Gross-Pitaevskii equation, modeling the piston by a moving Gaussian potential that is swept through the BEC (see Methods and Supplementary Note 1 for details). While our numerics are performed in three spatial dimensions, the expected dimensionality of the dynamics can be classified in terms of the nondimensional where N is the number of trapped atoms, λ = ω z /ω ⊥ is the cylindrically symmetric trap aspect ratio, a s is the scattering length, and a ho ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h=mω ? p is the transverse harmonic trap length scale. When d ( 1 the transverse dynamics are significantly constrained such that the BEC exhibits quasi-1D behavior. When d ) 1, the BEC is no longer geometrically constrained and is truly described by 3D dynamics. A simulation in the 1D regime with corresponding dimensional spatial, temporal, and density scales of L = a ho = 0.713 μm, T = 1/ω ⊥ = 0.695 ms, Γ ¼ 1=ð4πa s a 2 ho Þ ¼ 29:1 μm À3 , with ω z = 2π × 0.24 Hz and N = 8732 is shown in Fig. 3a with d = 0.068. A coherent soliton train or DSW is generated whose soliton edge speed quantitatively agrees with DSW theory 19 , as shown in the Supplementary Fig. 2. For more details on this analysis, see Supplementary Note 2.
To match experimental parameters, the simulations in Fig. 3b consider a large number of atoms, N = 410,005, with a tighter trap geometry (a tenfold increase in the longitudinal trap frequency) and the same scales L, T, and Γ as in the low atom number case above. The dimensionality parameter is now d ≈ 32, which places the dynamics well into the 3D regime.
Isosurface plots corresponding to Fig. 3b panels are shown in Fig. 4. Based on these simulations, the dynamics can be characterized as follows. At short times, a soliton train is initially formed, as would be expected in DSW theory (Fig. 4, t = 17 ms). Due to the large atom number-or, equivalently, weak transverse confinement-the soliton train rapidly undergoes a transverse, snake instability 32 , and vorticity emerges in the neighborhood of the shock front (Fig. 4, t = 23 ms, t = 44 ms). The snake instability manifests itself when the transverse harmonic oscillator length a ho exceeds the healing length by a factor of order one 33 . As seen at t = 58 ms and experimentally in Fig. 5, the plateau region hosts a variety of topological defects including vortex rings, lines, and vortex interactions. Experimental evidence for vorticity and soliton dynamics in the plateau regions obtained in absorption images after 10.1 ms expansion can be found in Fig. 5. Since the images in Fig. 5 are integrated along the x-axis, vortex rings appear as two dark dots with a faint connection between them 34 . Transformation of dark solitons into vortex rings has been argued to be responsible for an apparent inelasticity of collision events 35 . Furthermore, such rings have also been recently identified in bosonic 36 and fermionic 37 systems.
Both tight harmonic transverse confinement and quantum turbulence contribute to the highly inhomogeneous background through which the shock propagates. The transverse trap width a ho and the healing length in the condensate center are less than 1 μm, thus these characteristic length scales of the system are significantly smaller than the measured shock width Δx w ≈ 15.5 μm. Consequently, we argue that the shock experiences an effectively averaged or filtered turbulent field and an associated decrease in the effective density (and hence pressure). To quantify this, we use a technique from the large eddy simulation (LES) framework 38 to filter by convolving the hydrodynamic field from the simulation in Fig. 3 at t = 83 ms. The filtering is done with a kernel, K(r) = w −3 , for r in a cube with side length w centered at the origin, and zero otherwise. This leads to a reduction in the central, filtered condensate density, ρðr; tÞ ¼ ðK Ã ρÞðr; tÞ ¼ R Kðr À r′Þρðr′; tÞ dr′, when compared to the unfiltered, initial, unperturbed density. By choosing a cube side length of w = 5.35 μm, we obtain excellent agreement with our observed reduced speed of sound, c s,eff , independently fit in Fig. 2. This kernel size leads to a 3.35-fold reduction in the central, filtered condensate density. In this case, the effective sound speed satisfies The filtered hydrodynamic field is well-described by an exact viscous shock profile of the viscous shallow water equations as is shown in Fig. 6 where u ¼ ðK Ã ρuÞ= ρ is the Favre-averaged velocity 38 (see Methods). Support for this interpretation can be found in the literature on shocks in turbulent gases 38,39 . Turbulent dynamics modeled as dissipation at larger scales is the basis for the effective eddy viscosity in LES 38 . We view the quantum turbulence here in a similar fashion. In a coarse-grained description, the vortex excitations act as a reservoir for energy dissipated from the large amplitude wave excitations generated at the large-scale shock front. We compare an exact viscous shock profile of the viscous shallow water equations with our filtered 3D numerical simulations in Fig. 6. This reveals a shock structure consisting of a smooth transition from a nonzero subsonic flow (j uj<c s;eff ) ahead of the shock to supersonic flow (j uj ¼ v p >c s;eff ) behind the shock. The width of this transition is proportional to the effective dissipation experienced by the shock. The numerically observed nonzero mean flow ahead of and in the same direction as shock propagation is additional evidence for a decrease in the shock's pressure (hence density) as observed in numerical simulations of shock waves in a turbulent gas 39 . We stress that the only fitting parameters in Fig. 6 are the filtering length scale, the shock width, and the mean flow ahead of the shock. All remaining quantitiesmean density/velocity behind the shock, the shock speed, and the viscous shock profile-are completely determined by the Rankine-Hugoniot jump conditions for a piston in a viscous fluid. The quantitative agreement between viscous shock theory, the filtered, conservative BEC simulation, and experiment (Fig. 2) constitutes strong support for our conclusions. For additional analysis of the turbulence and associated energy production, see Supplementary Note 3 along with corresponding Supplementary  Figs. 3 and 4. Rarefaction waves. The emergence of a shock front as described above crucially depends on the presence of a finite background density through which the front propagates. In the absence of such a background density, the phenomenology is completely different and rarefaction waves emerge ( Fig. 1(b), 140 ms). Rarefaction waves, which are commonly discussed in the context of the shock tube problem in gas dynamics 24 or the dam break problem in shallow water 30 , form when a region of high density expands into a region of zero density (vacuum).
For dilute-gas BECs, we can study the formation of rarefaction waves in a clean, unperturbed setting by starting with 4.1 × 10 5 atoms confined in the left half of the harmonic trap. A repulsive barrier initially prevents the atoms from spreading out into the empty right half of the trap. When the barrier is suddenly removed, the BEC begins to spread out to the right and the front edge of the BEC density assumes a parabolic shape. For further details and experimental cross sections, see Supplementary Note 4. The self-similar expansion of the rarefaction front is in stark contrast to the dynamics of a sharp shock front. The parabolic shape is consistent with the predictions of these types of waves in a Bose-condensed gas 20,26,28,40 . The theory also predicts that the expanding edge propagates at twice the local speed of sound 20,26 .
We test this prediction in experiment by fitting the parabolic expanding edge of the BEC integrated cross section at various times during the expansion, and deduce from this the edge propagation speed (See Fig. 7 and Supplementary Fig. 5). The edge speed determined from our experiment (5.50 ± 0.33 mm s −1 ) and from matching numerics (5.16 ± 0.01 mm s −1 ) is in good agreement with the predicted behavior based on the 3D speed of sound (i.e., 2 × c s,bulk = 5.6 mm s −1 ). The data and numerics in Fig. 7 are also compared to 2 c s;bulk = ffiffi ffi 2 p ¼ 2 c s;1d ¼ 4:02 mm s −1 , which is the expected rarefaction edge propagation speed in a 1D channel (d ( 1). We see a clear deviation from this behavior, further indicating the fully 3D structure of our system (d ) 1) and the inapplicability of 1D DSW theory.

Discussion
In conclusion, we have observed and analyzed intriguingly rich dynamics in a quantum-mechanical piston shock. For our typical experimental parameters, the dynamics are described by dissipative, rather than dispersive, shock waves. The piston provides both the source of the shock front and the generation of superfluid quantum turbulence, manifested through the development of vortical patterns and dispersive waves, via a transverse, snake instability of a planar soliton train. We argue that an effective dissipation arises in a nominally inviscid superfluid as a consequence of the dissipation of large amplitude excitations from the large-scale shock front into small-scale vortex excitations.
Our experiments provide a versatile platform for the investigation of quantum turbulence, which is currently an area of intense research efforts in both cold atom and superfluid helium systems. Further studies into lower dimensional systems with similar geometries may also be of interest to determine the effect of dimensionality on quantum turbulence.

Methods
Experimental procedure. Our experimental setting consists of an elongated BEC of 4.1 × 10 5 87 Rb atoms confined in an optical dipole trap with trap frequencies {ω x , ω y , ω z } = 2π × {229, 222, 2.4} Hz. We estimate the atom number in the BEC by fitting integrated cross sections of the absorption images to the numerical ground state. Temperature is estimated to be ( T c , the critical temperature of the BEC, given no observable thermal cloud. We note that for velocities of 4 mm s −1 and lower, we see no noticeable atom loss in experiments. Agreement between experiment and zero-temperature GP numerics throughout the complete barrier sweeps further corroborates this point. The piston is generated by a repulsive laser beam of wavelength λ piston = 660 nm and an elliptical cross section of Gaussian waists {w x , w z } ≈ {68,11} μm. The barrier is swept from right to left (negative z-direction) using a high-speed mirror galvanometer. Effects of the initial acceleration of the galvanometer are avoided by initializing the barrier sweep far from the right edge of the BEC. This is done such that the acceleration range occurs outside of the BEC. Absorption images are taken after 2 ms expansion time to avoid image saturation.
Three-dimensional simulations. We perform (3 + 1)D numerical simulations of the Gross-Pitaevskii equation 31 . The initial condition consists of the ground state solution in the form ψ(r, t) = f(r)e −iμt in the presence of the harmonic trap without the barrier. μ is the chemical potential, determined by the total number of atoms N in the condensate according to where a ho ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi h=ðmω ? Þ p and Γ ¼ 1=ð4πa s a 2 ho Þ ¼ 29:1 μm À3 . We utilize a Fourier spatial discretization and a second order split-step method, exactly integrating the linear and nonlinear/potential terms separately. For the simulation of the experiment, we use a grid spacing of 0.071 μm on a box of size 8. Viscous shock fitting. Following Favre averaging 38 , filtering is performed by convolving the 3D density and momentum in the longitudinal direction with a cube of side length 5.35 μm. This cube size results in the downstream speed of sound c s,eff = 1.35 mm s −1 , used in the main text to explain the experimental and numerical observations. This cube size is a plausible filter in the context of large eddy simulation modeling as it is an intermediate length scale to the measured shock width (15.5 μm) and the healing length (~0.2 μm). The velocity is recovered by dividing the filtered momentum by the filtered density (Favre filtering). The shock profile is obtained from an exact, traveling wave solution of the 1D shallow water equations (the dispersionless limit of the Gross-Pitaevskii equation without a potential) with an additional, phenomenological viscous term for the nondimensional, filtered density ρ, velocity u, and viscosity parameter ν > 0. The traveling wave speed, plateau (rightmost) density and velocity are obtained from the Rankine-Hugoniot viscous shock conditions for a piston. The fitting parameters are the nondimensional viscosity ν = 25 and the traveling wave center z 0 = −25 μm, obtained by minimizing the sum of the absolute differences between the filtered and traveling wave superfluid velocity.
Code availability. All relevant code used for numerical studies in this work is available from the corresponding authors on reasonable request.

Data availability
All relevant experimental and numerical datasets in this work are available from the corresponding authors on reasonable request.  Rarefaction waves. A BEC initially confined to the left half of the trap is suddenly allowed to spread to the right. The plot shows the edge position vs. time, where experiment (green dots) and numerics (blue triangles) are plotted overlaid with expected results for 2c s,bulk (orange dashed) and 2c s,1d (red dot-dashed). Experimental data are mean ± s.d. for five runs at each measured time, where the green and blue solid lines are best fits to experiment and numerics, respectively. See text for more information