Observation of universal quench dynamics and Townes soliton formation from modulational instability in two-dimensional Bose gases

TL;DR

This study experimentally investigates the universal nonequilibrium dynamics of two-dimensional Bose gases quenched from repulsive to attractive interactions, revealing modulational instability, Townes soliton formation, and universal behaviors in a many-body quantum system.

Contribution

It demonstrates the formation of Townes solitons from modulational instability in 2D Bose gases and provides detailed real-time measurements of the dynamical process.

Findings

Observation of modulational instability leading to soliton formation

Density distributions match Townes soliton profiles

Universal behaviors in collapse and collision dynamics

Abstract

We experimentally study universal nonequilibrium dynamics of two-dimensional atomic Bose gases quenched from repulsive to attractive interactions. We observe the manifestation of modulational instability that, instead of causing collapse, fragments a large two-dimensional superfluid into multiple wave packets universally around a threshold atom number necessary for the formation of Townes solitons. We confirm that the density distributions of quench-induced solitary waves are in excellent agreement with the stationary Townes profiles. Furthermore, our density measurements in the space and time domain reveal detailed information about this dynamical process, from the hyperbolic growth of density waves, the formation of solitons, to the subsequent collision and collapse dynamics, demonstrating multiple universal behaviors in an attractive many-body system in association with the formation of a quasistationary state.

Equations

$$\mathcal{S}(\tilde{\tau })\approx \frac{\gamma }{{\bar{\gamma }}_i}[\tilde{S}(k_p,\tilde{\tau })-1],$$

$$\tilde{S}(k_p,\tau )\approx 1+\alpha \frac{{\bar{\gamma }}_i}{\gamma }{\sinh }^2(\gamma \tau ),$$

$${\Gamma }_s=\eta \frac{\hslash \pi }{m},$$

$$S(\mathbf{k})=\frac{1}{N}\sum _{\mathbf{q},{\mathbf{q}}^{\prime }}⟨{\hat{a}}_{\mathbf{q}+\mathbf{k}}^{†}{\hat{a}}_{\mathbf{q}}{\hat{a}}_{{\mathbf{q}}^{\prime }-\mathbf{k}}^{†}{\hat{a}}_{{\mathbf{q}}^{\prime }}⟩,$$

$$S(\mathbf{k})=⟨{\hat{a}}_{\mathbf{k}}^{†}{\hat{a}}_{\mathbf{k}}⟩+⟨{\hat{a}}_{-\mathbf{k}}{\hat{a}}_{-\mathbf{k}}^{†}⟩+⟨{\hat{a}}_{\mathbf{k}}^{†}{\hat{a}}_{-\mathbf{k}}^{†}⟩+⟨{\hat{a}}_{-\mathbf{k}}{\hat{a}}_{\mathbf{k}}⟩.$$

$${\hat{a}}_{\mathbf{k}}$$

$${\hat{a}}_{-\mathbf{k}}^{†}$$

$$\hat{H}=\frac{N\mu }{2}+\sum _{\mathbf{k}\ne 0}ϵ(k)({\hat{b}}_{\mathbf{k}}^{†}{\hat{b}}_{-\mathbf{k}}^{†}+{\hat{b}}_{\mathbf{k}}{\hat{b}}_{-\mathbf{k}})-\sum _{\mathbf{k}\ne 0}(\frac{{\hslash }^2{k}^2}{2m}+\mu ),$$

$${\hat{b}}_{\mathbf{k}}$$

$${\hat{b}}_{-\mathbf{k}}^{†}$$

$$S(\mathbf{k},\tau )$$

$$⟨{\hat{b}}_{0,\mathbf{k}}^{†}{\hat{b}}_{0,\mathbf{k}}⟩+⟨{\hat{b}}_{0,-\mathbf{k}}{\hat{b}}_{0,-\mathbf{k}}^{†}⟩$$

$$⟨{\hat{b}}_{0,\mathbf{k}}^{†}{\hat{b}}_{0,-\mathbf{k}}^{†}⟩=⟨{\hat{b}}_{0,\mathbf{k}}{\hat{b}}_{0,-\mathbf{k}}⟩$$

$$\zeta =⟨{\hat{c}}_{0,\mathbf{k}}^{†}{\hat{c}}_{0,\mathbf{k}}⟩+⟨{\hat{c}}_{0,-\mathbf{k}}{\hat{c}}_{0,-\mathbf{k}}^{†}⟩.$$

$$⟨{\hat{c}}_{0,\mathbf{k}}^{†}{\hat{c}}_{0,\mathbf{k}}⟩+⟨{\hat{c}}_{0,-\mathbf{k}}{\hat{c}}_{0,-\mathbf{k}}^{†}⟩=\frac{2}{e^{{ϵ}_i(k)/k_{B}T}-1}+1=\coth \frac{{ϵ}_i(k)}{2k_{B}T}.$$

$$S(k,\tau )=S_0(k)\left[1+\frac{2{ϵ}_k\hslash ({\gamma }_i+\gamma )}{ϵ(k{)}^2}{\sinh }^2\frac{ϵ(k)\tau }{\hslash }\right],$$

$$S_0(k)=\frac{{\hslash }^2{k}^2}{2m{ϵ}_i(k)}\coth \frac{{ϵ}_i(k)}{2k_{B}T}$$

$$\mathcal{S}(\tilde{\tau })=\zeta \left[\tilde{S}(k_p,\tilde{\tau })-1\right],$$

$$\mathcal{S}(\tilde{\tau })=2{\sinh }^2(\tilde{\tau }).$$

$$S(k,\tau)=S_{0}(k)\Big{[}1+\frac{\epsilon_{i}(k)^{2}-\epsilon^{2}(k)}{\epsilon(k)^{2}}\sin^{2}\frac{\epsilon(k)\tau}{\hbar}\Big{]},$$

$$\frac{d{\bar{N}}_a}{d{\tau }_2}=-\frac{L_3}{9\sqrt{3}{\pi }^3\sigma ({\tau }_2{)}^4{l}_z^2}{\bar{N}}_a^3,$$

$${\bar{N}}_a({\tau }_2)=\frac{{\bar{N}}_a(0)}{\sqrt{1+\frac{2L_3{\bar{N}}_a(0{)}^2}{27\sqrt{3}{\pi }^3{l}_z^2|\dot{\sigma }|}\left[\frac{1}{({\sigma }_0+\dot{\sigma }{\tau }_2{)}^3}-\frac{1}{{\sigma }_0^3}\right]}},$$

Full text

Predicting the evolution of multidimensional nonlinear systems under attractive interactions is a challenging task, owing to the instability to collapse [1, 13, 3]. Bright solitons are remarkable stationary states, established when the self-focusing effect responsible for collapse is exactly compensated by wave dispersion [4, 5]. In uniform two-dimensional (2D) systems with standard cubic interactions, such as Kerr medium [4, 3] or matter waves formed by weakly interacting 2D Bose gases [7, 1], however, such intricate balance cannot be fulfilled except when a wave packet possesses a critical norm (or total atom number) known as the Townes threshold and a specific waveform known as the Townes profile [4, 9, 10] – only at which bright solitons can form. A Townes soliton is predicted to be unstable as long as the norm deviates from the Townes threshold [3, 9]. Despite extensive interest in multidimensional bright solitons [11, 12, 9, 10, 13], including recent advancements on 2D spatial solitons in various nonlinear optics settings [11, 3, 14], an experimental realization and characterization of Townes solitons has remained elusive.

In soliton formation dynamics, modulational instability (MI) is a ubiquitous mechanism that causes the amplification of initial wave disturbances and fragmentation into solitary waves [1, 15, 10, 11]. In one-dimensional (1D) systems, MI is responsible for the formation of stable soliton trains, for example, in nonlinear fiber optics [18, 19, 20], in 1D Bose gases [21, 22, 10, 11], or in Bose-Einstein condensates in optical lattices [23, 24, 25]. In higher spatial dimensions, transverse MI and wave fragmentation were studied in various types of bulk nonlinear optical media [26, 27, 14]. However, detailed dynamics of multidimensional MI and its possible connection to the universal formation of a quasistationary state, the unstable Townes solitons, have not been clearly demonstrated.

Using ultracold 2D Bose gases, we show that universal MI dynamics supports the critical formation of Townes solitons. Starting with a 2D superfluid of an initial density $n_{i}$ and quenched to an attractive interaction $g_{f}<0$, we show that MI causes collective modes with a wave number around $k_{\mathrm{MI}}=\pi/\xi$ associated with the interaction length $\xi=\pi/\sqrt{2n_{i}|g_{f}|}$ to grow predominantly [28, 29, 30], fragmenting the superfluid as illustrated in Fig. 1(a). Intriguingly, the characteristic atom number in each wave packet $\sim n_{i}\xi^{2}=\pi^{2}/2|g_{f}|$ well approaches the Townes threshold $N_{\mathrm{th}}=5.85/|g_{f}|$ [10], thus opening up a pathway for Townes soliton formation. This relation should apply universally for any $n_{i}$ and $g_{f}$ provided no other length scales set in. We note that, due to the scaling symmetry in 2D, a Townes soliton forms under a scale-invariant profile [4] and MI can set the physical length scale $\xi$ that depends only on the product of $n_{i}$ and $|g_{f}|$.

In this Letter, we report the observation of universal dynamics and Townes soliton formation in quenched 2D Bose gases. We observe MI that breaks up an otherwise large 2D sample into fragments universally around the Townes threshold. We clearly identify solitary waves whose density distributions agree well with the Townes profiles – the stationary state solution of the 2D nonlinear Schrödinger equation (NLSE) [4, 10]. Our measurement further reveals universal solitary wave dynamics governed by the MI timescale and a universal scaling behavior in the density power spectra, allowing us to clearly identify a distinct time period for the unstable growth of density waves (while conserving total atom number), followed by a short era of wave collapse and soliton formation.

We begin the experiment with a uniform 2D Bose gas deep in the superfluid regime [31], which is formed by $N\approx 1.5\times 10^{4}$ cesium atoms trapped inside a quasi-2D box potential. The atomic surface density is $n_{i}\approx 5/\mu$m2 and the surface area is controlled by a wall-like potential that is removed following the interaction quench [31]. The tight vertical ($z$) confinement of the box freezes all atoms in the harmonic ground state along $z$ axis. The trap vibrational level spacing ($\omega_{z}=2\pi\times 1750~{}$Hz) is more than 2 orders of magnitude larger than the attractive interaction energy studied, ensuring that the observed wave dynamics is effectively 2D [32]. The interaction strength $g=\sqrt{8\pi}a/l_{z}$ is controlled by the tunable s-wave scattering length $a$ via a magnetic Feshbach resonance [2], and $l_{z}=208~{}$nm is the vertical harmonic oscillator length; $g=g_{i}=0.115$ is the initial interaction strength.

We perform an interaction quench (in 1 ms) to various $g=g_{f}$ in samples with a large surface area $A\approx(60~{}\mu$m$)^{2}$. Following a variable hold time $\tau$, we perform single-shot absorption imaging to record the sample density distribution as shown in Fig. 1(b). Around 30 samples are imaged for ensemble analyses. In a short hold time, we observe density blobs randomly clumping up throughout a sample. The sizes of the blobs are smaller with larger $|g_{f}|$. At a longer hold time, $\tau\geq 30~{}$ms, the number of observed blobs reduces, becoming more isolated, although the mean size and atom number of the surviving blobs remain nearly unchanged for $\tau\leq 200~{}$ms (see Fig. S2 in Ref. [31]). The same quench protocol applied to samples in a three-dimensional trap ($\omega_{z}\approx 2\pi\times 100~{}$Hz with a weak radial trap frequency $\omega_{r}\approx 2\pi\times 13$ Hz) results in no observed density blobs.

We characterize these isolated blobs (solitary waves) and compare their atom number with the Townes threshold. We approximate their density profiles by 2D Gaussians [10] and fit the mean size $\sigma$ and atom number $\bar{N}_{a}$ [31]. Within the interaction range studied, $-0.004\geq g_{f}\geq-0.034$, in Fig. 1(c) we find that all rescaled atom numbers $\bar{N}_{a}|g_{f}|$ fall around $N_{\mathrm{th}}|g_{f}|=5.85$, giving a mean $\overline{\bar{N}_{a}|g_{f}|}=6(1)$. Interestingly, the standard deviation of the atom number $\delta N_{a}$ scales with $g_{f}$ accordingly, giving a mean fluctuation $\overline{\delta N_{a}|g_{f}|}=$3.2(5)$\sim 0.5\overline{\bar{N}_{a}|g_{f}|}$. We believe that the number variation around the threshold results from the energy-time uncertainty relation, as we later show that these blobs form in a timescale $\sim\gamma^{-1}$, where $\hbar\gamma=\hbar^{2}n_{i}|g_{f}|/m$ is the interaction energy, $m$ is the atomic mass, and $\hbar$ is the reduced Planck constant. Moreover, the size of these solitary waves also agrees with the prediction $\sigma\approx\xi$ in Fig. 1(d), indicating that MI provides the length scale for the formation of blobs. A size discrepancy at $\gamma\approx 2\pi\times 1.2$ Hz ($g_{f}=-0.004$) could likely be attributed to the influence of a very weak horizontal corrugation in the vertical confining potential.

To confirm that the quench-induced solitary waves indeed form Townes solitons, we compare their density distributions with the scale-invariant, isotropic Townes profile $n(r)=\alpha^{2}|\phi(\alpha r)|^{2}/(2|g_{f}|)$, where $\alpha=\sqrt{2n_{0}|g_{f}|/|\phi(0)|^{2}}$ is a scale factor. Given the peak density $n_{0}$, the characteristic size and density profile of a Townes soliton are uniquely determined. Here, the radial function $\phi(\tilde{r})$ is the stationary solution of the scale-invariant 2D NLSE [4, 10]. The stationary profile $\phi(\tilde{r})$ has been evaluated numerically [4], giving $|\phi(0)|\approx 2.207$, and the norm $\int|\phi(\tilde{r})|^{2}d\tilde{\mathbf{r}}\approx 11.7$ sets the Townes threshold.

Since MI sets the length scale during the soliton formation process, $\alpha\sim\xi^{-1}$, the peak density of most solitons should be comparable to the initial sample density $n_{0}\approx n_{i}$. Figure 2 plots three isolated solitary waves of similar peak density that are randomly chosen from quenched samples ($g_{f}=-0.034$) at a long hold time $\tau=100~{}$ms. Their individual density distributions, as well as the averaged radial density profile, indeed agree fairly well with the expected Townes profiles with no fitting parameters. More agreement with the Townes profiles is discussed in Ref. [31], where a single array of well-isolated solitons can form in an elongated sample following an interaction quench. Our observation confirms that Townes solitons can prevail from MI and explains why many randomly formed solitary waves, as observed in Fig. 1, are long-lived. We note that, in the absence of MI and fragmentation, a wave was observed to collapse only partially to a Townes profile [3].

We now turn to study the universal dynamics during the soliton formation process. We focus on the time evolution of density power spectrum [6] $S(k,\tau)=\langle|n(k,\tau)|^{2}\rangle/N$ in the spatial frequency domain (momentum space), using the Fourier transform of the sample density distribution. Here $k=|\mathbf{k}|$ is the momentum wave number of the azimuthally averaged spectrum, $N$ is the total atom number, and $\langle\cdots\rangle$ denotes ensemble averaging. In the power spectra [Fig. 3(a)], we clearly observe rapid growth of a nonzero momentum peak at short hold time (marked by an arrow), indicating the development of density waves at a dominant length scale throughout the sample. The nonzero momentum peak then dissipates at longer hold time until the power spectrum finally becomes monotonic and stationary, which signifies the collapse of density waves and fragmentation of the sample into solitons that later becomes uncorrelated in coordinate space.

The evidence of MI-induced wave amplification at different interaction strengths is best illustrated when we plot the relative growth spectra $\tilde{S}(k,\tau)=S(k,\tau)/S(k,\tau_{0})$, normalized by the initial power spectrum at $\tau_{0}=2~{}$ms. This allows us to determine which mode has the largest growth rate. For different samples in Figs. 3(b) and 3(c), the momentum peak is clearly visible within $0.2/\mu$m$<k<1/\mu$m at short hold time. The growth patterns look similar for samples with different $g_{f}$, although the peak location, height, and the evolution timescale vary. We identify the peak wave number $k_{p}$ and find consistency with the prediction from MI in Fig. 3(d).

Another remarkable prediction from MI is that, regardless of the dimensionality of the system, the power spectrum at $k_{p}$ exhibits a universal time and amplitude scaling behavior with respect to the interaction timescale $\gamma^{-1}$. This is based on extending the Bogoliubov phonon picture to the regime under attractive interactions, which predicts that collective modes of opposite momenta are generated in pairs, as seeded from initial density perturbations, and subsequently form density waves along the associative directions while being amplified at a rate $\gamma$ until significant depletion of the ground state atoms occurs [31].

In Fig. 3(e), we experimentally test the scaling relation in the peak growth spectra $\tilde{S}(k_{p},\tau)$, covering the entire time period. We summarize the scaling relation as follows:

[TABLE]

where $\tilde{\tau}=\gamma\tau$ is the scaled time and $\mathcal{S}(\tilde{\tau})$ is the scaled spectrum. In the amplitude scaling, we normalize $\gamma$ with the mean initial interaction energy unit $\bar{\gamma}_{i}=306~{}$s*-1* before the quench [31]. We show that six power spectra, each with different $\gamma$, can collapse onto a single curve over a surprisingly long scaled time $\tilde{\tau}<10$. The only exception is the spectrum at $g_{f}=-0.004$, where we have used $\gamma^{*}=3.2\gamma$ to force its collapse within $\tilde{\tau}\leq 0.8$. This different behavior is likely due to a weak horizontal trap corrugation influencing the dynamics, as noted in Fig. 1(d).

From this universal spectrum, we identify two distinct regimes of dynamics, divided by a critical time $\tilde{\tau}_{c}\approx 0.8$ as shown in Fig. 3(e). We label the time period $\tilde{\tau}\leq\tilde{\tau}_{c}$ for MI with an amplified (hyperbolic) growth of density waves [31],

[TABLE]

where $\alpha=1.5(1)$ is determined from a fit to $\mathcal{S}(\tilde{\tau})$ for $\tilde{\tau}\leq 0.5$; $\alpha=2$ is obtained from the theory calculation for $\tilde{\tau}\ll 1$, neglecting the depletion of ground state atoms, dissipation, or interaction between the collective modes. Beyond $\tilde{\tau}\geq\tilde{\tau}_{c}$ after $\mathcal{S}(\tilde{\tau})$ reaches the order of unity, dynamics enters the second phase, decaying with a time constant $\Delta\tilde{\tau}\sim 0.8$ and transitioning to a slowly evolving, quasistationary behavior.

Our data suggest the existence of a universal time and amplitude scaling behavior and a limit for the amplified density wave, followed by a universal dynamics for the wave collapse and soliton formation. For $g_{f}=-0.004$, however, $\tilde{S}(k_{p},\tau)$ remains slowly growing within $1\leq\gamma^{*}\tau\leq 10$, suggesting a less severe wave collapse.

Following the observed density wave collapse dynamics, a 2D sample fragments into many solitary wave packets of characteristic size $\sigma\sim\xi$. As seen in Fig. 1, due to the close proximity of many wave packets (also with characteristic distances $\sim\xi$), collisions between them may induce collapses that lead to rapid atom number loss. Here, we show that the MI timescale continues to dominate the collision dynamics and the total atom number loss in a quenched 2D Bose gas.

In a recent study of 1D soliton collisions [9], it is shown that merger occurs when solitons of similar phases collide, and two solitons of opposite phases appear to repel each other. In 2D, the merger of two solitonlike wave packets should lead to a new atom number $N_{a}>N_{\mathrm{th}}$. This can induce collapse that quickly removes the merged soliton. For solitons or density blobs formed by MI with randomly seeded density waves in a large 2D sample, one expects no fixed phase relationship between neighbors. Merger can thus occur randomly throughout the sample and the total atom number loss may reveal the scaling of soliton binary collision losses.

In Fig. 4, we examine the total atom number loss for the samples quenched to different $g_{f}$ (Fig. 1). We observe the onset of loss in each sample at a time $\tau$ corresponding roughly to the critical time $\tilde{\tau}_{c}$, a behavior similarly observed for MI in 1D [10]. Beyond the critical time, we confirm that the loss curves can be well-captured by a simple two-body loss model, $\dot{N}/N=-\Gamma_{\mathrm{2body}}N/A$. We attribute this behavior to the dominance of binary collisions between solitons or density blobs that trigger collapse and atom number loss; without triggered collapse, the usual three-body recombination loss should be negligible [3]. In Fig. 4(b), we find a linear dependence on $|g_{f}|$ in the measured loss coefficients $\Gamma_{\mathrm{2body}}$. This in fact suggests a constant binary loss coefficient $\Gamma_{s}$ for colliding wave packets, where $\Gamma_{s}=\Gamma_{\mathrm{2body}}\bar{N}_{a}$ [31] and $\bar{N}_{a}|g_{f}|\approx 6$ is the measured universal number for solitary waves formed by MI.

We suspect this universal loss behavior results from MI scaling and 2D scale invariance, which suggests a constant binary loss coefficient,

[TABLE]

independent of the interaction parameters ($n_{i},g_{f})$. This is because the collision rate $\Gamma_{s}\sim\sigma\bar{v}$ and the dependences on length scales in the linear cross section $\sigma\approx\xi$ and relative velocity $\bar{v}\approx\sqrt{2}\hbar\pi/m\xi$ cancel each other; the constant $\eta\approx\sqrt{2}$ is estimated for $\sim$50% probability of merger per collision or, equivalently, on average one soliton or density-blob loss per collision event. If Eq. (3) holds, we expect a collision lifetime $1/n_{s}\Gamma_{s}\sim\gamma^{-1}$, where $n_{s}=n_{i}/\bar{N}_{a}$ is the initial soliton/blob density [31]. This suggests that wave collapse and binary collision likely take place at the same rate during the second phase of the density evolution.

To unambiguously confirm the universality of collision dynamics, in Fig. 4(c) we deduce $\Gamma_{s}$ independently using experimentally determined values ($\Gamma_{\mathrm{2body}}$, $\bar{N}_{a}$) at each $g_{f}$. Our results conform very well with the prediction by Eq. (3), giving a mean $\eta\approx 1.5$. We emphasize here that the loss coefficients universally depend only on the physical constants $\hbar/m$ is a remarkable manifestation of MI and scale-invariant symmetry in 2D. The observations in Figs. 3 and 4 together confirm that interaction quench dynamics leads to Townes soliton formation at $\tau\gtrsim\gamma^{-1}$, followed by collision (that induces collapse) also at the same timescale $\gamma^{-1}$ universally governed by MI.

In summary, we study the universal nonequilibrium dynamics in degenerate 2D Bose gases quenched from repulsive to attractive interactions, and observe the dynamical formation of Townes solitons from MI. Townes solitons are observed to be collisionally unstable. However, further stabilization and manipulation may be possible [37, 38, 39, 11]. We note that the initial shape and finite sample size can be further manipulated to invoke strong boundary effects in quench-induced MI [31]. Using slow interaction ramps may also induce MI dynamics deviating from the reported universal behavior. Soliton formation may be disrupted, leading to, for example, only partial collapse to the Townes profiles [3]. Lastly, we comment that controlled formation of 2D solitons via pair-production processes in MI may find future applications in matter-wave interferometry [40, 41], or even in the generation and distribution of many-body entanglement [42, 43, 44, 45].

We acknowledge discussions with Qi Zhou, Chih-Chun Chien, Sergei Khlebnikov, and Chris Greene. We thank Cheng Chin for discussions and many critical comments. We are grateful to H. J. Kimble for instrument loan since the early stage of this project. We thank May Kim, Yiyang Fang, and Wuxiucheng Wang for laboratory assistance. This project is supported in part by the Purdue Research Foundation, the W. M. Keck Foundation, the NSF Grant No. PHY-1848316, and the U.S. Department of Energy (Grant No. DE-SC0019202).

Appendix A Supplemental Material

A.1 Experimental procedures

A.1.1 Formation of uniform 2D Bose gases

Figure S1 shows the schematics of the experiment. Our uniform 2D Bose gas is confined in a quasi-2D box potential formed by all repulsive optical dipole beams. The vertical box confinement is provided by a single node of a repulsive standing-wave potential (2$~{}\mu$m periodicity) along the vertical ($z$-) direction. The node defines the horizontal 2D plane. The measured vertical trap frequency in the node is $\omega_{z}/2\pi=1.75(2)~{}$kHz $\gg(k_{B}T,|\mu|)/\hbar$, freezing all atoms in the harmonic ground state. Here $k_{B}$ is the Boltzmann constant, $T\leq 8~{}$nK is the temperature measured in time-of-flight, and $\mu$ is the chemical potential. The horizontal boundary of the box is formed by a tightly focused repulsive optical beam that scans the box boundary to form a time-averaged wall potential. The beam is controlled by a pair of acousto-optic deflector and is projected through the same microscope objective (numerical aperture N.A. = 0.37) used for imaging. Optical resolution of our imaging system is $\lesssim 2\mu$m, determined using the smallest density blobs identified in the experiment.

The 2D Bose gas is prepared deeply in the 2D superfluid regime with a phase space density $n_{i}\lambda_{\mathrm{dB}}^{2}\gtrsim 15$ [1]. Here $n_{i}\approx 5/\mu$m2 is the initial 2D density and $\lambda_{\mathrm{dB}}$ is thermal de Broglie wavelength. Finite temperature mainly contributes initial thermal phonon fluctuations in a superfluid sample that seed the growth of density waves from modulational instability.

We perform the interaction quench while simultaneously removing the wall potential to allow hot collision products (due to soliton collapses) to eject out of the experiment field of view. When the wall potential is removed, the residual global trap frequency in the horizontal plane, due to a weak magnetic trap, is experimentally determined to be $\omega_{r}/2\pi<1.5~{}$Hz through superfluid dipole oscillations.

A.1.2 Magnetic two-body interaction tuning

The scattering length is tuned using a magnetic Feshbach resonance [2, 3]. We identify zero scattering length at the magnetic field $B=17.120(6)$G, using the superfluid in-situ size as well as the expansion rate in a 2D time-of-flight. We adopt the formula [3] $a(B)=(1722+1.52B/G)\left(1-\frac{\Delta B}{B/G-B_{0}}\right)$ for the scattering length conversion, where $\Delta B=28.72$ and $B_{0}=-11.60$ is adjusted to shift the zero-crossing to the measured value. The interaction strength is determined as $g=\sqrt{8\pi}a/l_{z}$, where $l_{z}=208(1)~{}$nm is the vertical harmonic oscillator length. The uncertainty ($\pm\delta g$) in $g$ is primarily contributed by the uncertainty in the magnetic field at the zero-crossing of scattering length. Within the range of reported interaction strengths $-0.004\geq g_{f}\geq-0.034$, we have $0.0005\lesssim\delta g\lesssim 0.0007$.

A.1.3 Solitary wave characterization

To identify and characterize a solitary wave at long hold time in samples as shown in Fig. 1, we implement 2D peak searching algorithms to search for isolated blobs with peak density above a threshold around $80$% of the initial density $n_{i}$, and then fit each blob with a 2D Gaussian in approximation to the Townes profile [4]. From each fit, we obtain the root-mean-square (r.m.s.) diameters and extract the atom number $N_{a}$ under the fitted density profile. We discard fit results that fail to converge in size or have exaggerated aspect ratios that often result from two or more blobs being too close to each other. We then evaluate the mean size $\sigma$ and atom number $\bar{N}_{a}$. Examples at long hold time $\tau\geq 50~{}$ms can be found in Fig. S2, which shows no signs of rapid collapse or expansion in the mean size and atom number following a single interaction quench.

A.1.4 Formation of a Townes soliton array

To fully demonstrate that Townes solitons can generally form from our quench recipe, we induce MI in another set of 2D samples initially confined and then released from a narrow rectangular wall potential as shown in Fig. S3(a). We adjust the short side of the samples to be comparable to the MI length scale, so that a single array of solitons can form following the interaction quench. This avoids close proximity with many neighboring solitary waves or dispersing blobs, visible in large samples shown in Fig. 1(b). From these narrow samples, we clearly observe well-isolated solitary waves in almost every experiment repetition as shown in Fig. S3(a) at hold time $\tau\geq 30~{}$ms. In these solitary waves, we find ubiquitous agreement with the Townes profiles [Fig. S3(b-c)]. Together with Fig. 2, our observation confirms that Townes solitons can prevail from MI.

A.2 Quench-induced dynamics in the density power spectrum

In a 2D gas with uniform mean density distribution, the density power spectrum at finite $k$ is essentially the density static structure factor, which is the Fourier transform of the density-density correlation function [5]. In the following, we discuss the quench evolution of the static structure factor (density power spectrum) measured in our samples, while neglecting density perturbations due to boundary effects.

A.2.1 Theory of density structure factor after an interaction quench to $g_{f}<0$

For a Bose superfluid with initial density $n_{i}$ immediately following the interaction quench, we expect density waves with wavenumber $0<k<\sqrt{4|g_{f}|n_{i}}$ to growth unstably since the usual Bogoliubov dispersion becomes purely imaginary. There is no straightforward theory for evaluating quench evolution at all hold time $\tau$. To gain insights, here we analytically evaluate the quench dynamics only in the very early stage when most of the atoms still remain in the zero momentum state. We focus on the time-evolution of the static structure factor [6, 7]. Analytically, it can be evaluated as

[TABLE]

where $\hat{a}_{\mathbf{k}}$($\hat{a}^{\dagger}_{\mathbf{k}}$) stands for the annihilation (creation) operator of a momentum state $|\mathbf{k}\rangle$ and $N$ is the total atom number. At very short hold time $\tau\ll\gamma^{-1}$, the Bose gas is still primarily populated by ground state atoms ($\hat{a}^{({\dagger})}_{0}\approx\sqrt{N}$), where $\hbar\gamma=\hbar^{2}n_{i}|g_{f}|/m$ is the interaction energy, $m$ is the atomic mass, and $\hbar$ is the reduced Planck constant. The structure factor reduces to

[TABLE]

We perform the following transformation for momentum state within the range $0<|\mathbf{k}|<\sqrt{4|g_{f}|n_{i}}$, expressing the momentum state operator with a set of bosonic mode operators $\hat{b}_{\mathbf{k}}~{}(\hat{b}^{\dagger}_{-\mathbf{k}})$ as

[TABLE]

Here, we set the coefficients $u_{k}=\sqrt{\frac{\hbar\gamma}{2\epsilon(k)}+\frac{1}{2}}$ and $v_{k}=\sqrt{\frac{\hbar\gamma}{2\epsilon(k)}-\frac{1}{2}}$, $\epsilon(k)=\sqrt{|\epsilon_{k}^{2}-2\epsilon_{k}\hbar\gamma|}$ is the imaginary part of the Bogoliubov energy, and $\epsilon_{k}=\hbar^{2}k^{2}/2m$ is the single particle dispersion. $\hat{b}_{\mathbf{k}}~{}(\hat{b}^{\dagger}_{-\mathbf{k}})$ obeys the usual bosonic commutation relation. Using procedures similar to the standard Bogoliubov transformation, we can recast the weakly-interacting Hamiltonian into the following form

[TABLE]

where the summation runs over half of momentum space and $\mu=-\hbar\gamma$ is the chemical potential. Note that, under this transformation, new excitations are generated (and also annihilated) in pairs as hold time increases. In the Heisenbeg picture, these operators obey a set of coupled equations of motion, $\dot{\hat{b}}_{\mathbf{k}}=\frac{i}{\hbar}[\hat{H},\hat{b}_{\mathbf{k}}]=-\frac{i}{\hbar}\epsilon(k)\hat{b}^{\dagger}_{-\mathbf{k}}$ and its Hermitian conjugate, which lead to the following solution

[TABLE]

and $\hat{b}_{0,\mathbf{k}}(\hat{b}^{\dagger}_{0,-\mathbf{k}})$ is the bosonic mode operator at time $\tau=0$ right after the interaction quench. Plugging this solution into Eq. (S3) and evaluate the structure factor Eq. (S2), we then find the following time-dependent evolution

[TABLE]

Here, the first line contains mode contributions that are seeded by the initial bosonic mode populations right after the quench. These modes grow ‘hyperbolically’ in the early stage of the quench dynamics. The second and third lines contain the contributions from mode populations that are generated or annihilated from the interaction quench, where the hyperbolic term is expected to vanish, leaving only the constant term (see below).

A.2.2 Amplification of density waves from density fluctuations prior to the quench

By using Eq. (S3) and the Bogoliubov transformation, we can further relate the expectation values of the bosonic modes in Eq. (S7) to those of the phonon modes before the interaction quench. We find

[TABLE]

and

[TABLE]

Here, $\hat{c}_{0,\mathbf{k}}(\hat{c}^{\dagger}_{0,\mathbf{k}})$ is the phonon annihilation (creation) operator, $\epsilon_{i}(k)=\sqrt{\epsilon_{k}^{2}+2\epsilon_{k}\hbar\gamma_{i}}$ is the Bogoliubov dispersion at interaction $g_{i}>0$ prior to the quench, $\gamma_{i}=\hbar n_{i}g_{i}/m$, and we have used $\langle\hat{c}^{\dagger}_{0,\mathbf{k}}\hat{c}^{\dagger}_{0,-\mathbf{k}}\rangle=\langle\hat{c}_{0,\mathbf{k}}\hat{c}_{0,-\mathbf{k}}\rangle=0$ in the above relation since there is no source or sink for phonons in our 2D gas. It is clear that the bosonic mode population Eqs. (S8-S9) is seeded by the initial thermal phonon population and zero-point fluctuations

[TABLE]

Using the above relations and keeping only the wavenumber $k$-dependence, Eq. (S7) can now be simplified as

[TABLE]

where the overall factor

[TABLE]

is exactly the equilibrium static structure factor [8] right before the interaction quench.

A.2.3 The scaling behavior

In Eq. (S12), the second term in the bracket represents contributions from the MI-amplified density waves, suggesting density waves at wavenumber $k_{p}=\sqrt{2|g_{f}|n_{i}}$ have the largest amplification rate $\epsilon(k_{p})/\hbar=\gamma$. Thus, at short hold time $\tau\ll\gamma^{-1}$, we expect the growth of density power spectrum $\tilde{S}(k_{p},\tau)\equiv S(k_{p},\tau)/S(k_{p},0)$ to obey the following scaling relation

[TABLE]

where $\tilde{\tau}=\gamma\tau$ is the scaled time, $\zeta=\frac{\gamma}{\gamma_{i}+\gamma}$ is a dimensionless amplitude scaling factor, and $\mathcal{S}(\tilde{\tau})$ is the scaled spectrum; $\zeta\approx\gamma/\gamma_{i}$ when $\gamma_{i}\gg\gamma$. The scaled spectrum should display a universal hyperbolic growth at short hold time

[TABLE]

A.2.4 Relationship with quench dynamics at $g_{f}>0$

We note that the time-dependent density power spectrum Eq. (S12) is essentially the analytical continuation of quench-induced Sakharov oscillations in the static structure factor at $g_{f}>0$ [6],

[TABLE]

where now a coherent, sinusoidal oscillation in the structure factor replaces the hyperbolic growth. The calculation is in principle valid for all hold time $\tau$ for $g_{f}>0$, in contrast to the case of $g_{f}<0$, provided no global trap effect or other effects set in.

A.2.5 Experimental test of the scaling behavior in the density power spectrum

We experimentally test the scaling behavior Eq. (S14) over an extended time period. We have empirically searched for the best amplitude scaling relation $\zeta$ for a larger time range $0<\tilde{\tau}<5$. We find that the spectra scale the best with $\zeta\propto\gamma$, keeping explicit dependence on $n_{i}|g_{f}|$. We adopt a simple amplitude scaling factor $\zeta=\gamma/\bar{\gamma}_{i}$, using mean $\bar{\gamma}_{i}=\hbar\bar{n}_{i}g_{i}/m=306$ s*-1* evaluated from all samples (mean $\bar{n}_{i}=6/\mu$m2 and standard deviation $\delta n_{i}=1/\mu$m2). We note that an alternative amplitude scaling factor $\zeta=\gamma/(\bar{\gamma}_{i}+\gamma)$ close to the exact form in Eq. (S14) gives a similar data collapse.

A.3 Soliton collision dynamics

A.3.1 Soliton binary collision coefficient $\Gamma_{s}$

Binary collisions between solitons can lead to merger [9], which makes soliton atom number $N_{a}>N_{\mathrm{th}}$ much greater than the Townes threshold and induces collapse and rapid atom loss. Soliton binary loss behavior can be effectively described by $\dot{N}_{s}/N_{s}=-\Gamma_{s}n_{s}$, where $N_{s}$ is the number of solitons, $n_{s}=N_{s}/A$ is the surface density, and $\Gamma_{s}$ is the 2D binary loss coefficient as discussed in the main text. For total atom number loss dominated by soliton binary collision loss, we should have $\dot{N}/N=-\Gamma_{\mathrm{2body}}n$, where $n$ is the atom number density. We convert the soliton number into total atom number $N$ using $N_{s}=N/\bar{N}_{a}$ [10, 11] and assuming $\bar{N}_{a}$ is approximately constant over the experiment time. We find the simple relation $\Gamma_{s}=\Gamma_{\mathrm{2body}}\bar{N}_{a}$, relating the measured loss coefficient $\Gamma_{\mathrm{2body}}$ to the binary collision coefficient $\Gamma_{s}$.

A.3.2 Collision lifetime of solitons

Following the determination of universal collision dynamics, we also obtain the collision lifetime of solitons right after the formation process near $\tilde{\tau}\approx\tilde{\tau}_{c}$. We calculate the lifetime $\Delta\tau_{s}=\frac{1}{\Gamma_{s}n_{s}}\approx\frac{m\bar{N}_{a}}{\eta\pi\hbar n_{i}}\approx\gamma^{-1}$, where we have used $n_{s}=n_{i}/\bar{N}_{a}$, the measured universal threshold $\bar{N}_{a}|g_{f}|=6$ (Fig. 1), and the measured constant $\eta=1.5$ (Fig. 4).

A.4 Collapse and expansion dynamics of an unstable Townes soliton

This section is dedicated for gaining further information about the evolution of unstable Townes solitons. We determine their collapse and expansion time scales by driving the atom number more than three times away from the Townes threshold via a second interaction quench.

A.4.1 Second interaction quench

We apply a second interaction quench after solitons form at a sufficiently long hold time $\tau=50~{}$ms, hold for an additional time $\tau_{2}$, and perform imaging. In a quench path labeled (i) in Fig. S4(a) to a less attractive interaction, we induce immediate soliton expansion, during which $\bar{N}_{a}$ remains constant, as expected, but the number of solitons $N_{s}$ observed from the ensemble measurements greatly reduces; see Fig. S4(c) and inset. For a reversed quench path (ii) to a more attractive interaction, solitons collapse. In a short time scale $\gamma^{\prime-1}\sim 20~{}$ms corresponding to the interaction energy difference between the two quenches, we observe rapid reductions in both $\sigma$ and $\bar{N}_{a}$.

In the collapse dynamics, we tentatively attribute the atom number loss within a soliton to few-body inelastic collisions that are primarily due to three-body recombination. However, in Fig. S4(c), we obtain an unphysical loss coefficient that is five orders of magnitude larger than that measured in thermal samples [3]. A rapid three-body loss rate was also reported in the collapse within 1D solitons in a related experiment setting [10]. We believe that a more likely explanation for the rapid loss may be due to a combination of three-body loss under higher local collapse density and collective matter-wave ejection out of the solitons during the collapse [12, 13, 14, 15], which is challenging to detect given limited image resolution and signal-to-noise. Interestingly, the collapse seemingly slows down as $\bar{N}_{a}$ drops to the new threshold value. We suspect either the atom loss has regulated the collapse, perhaps due to collective wave emission, or solitons with initial $N_{a}$ close to the new threshold survive. From the inset of Fig. S4(c), the observed continuing decrease of soliton number beyond $\tau_{2}\geq 36~{}$ms may hint more of the latter case.

In either quench paths (i) and (ii), soliton evolution is clearly visible within the interaction time scale $\gamma^{\prime-1}$. For the dynamics after just a single quench, such rapid evolutions in mean atom number and soliton size are not observed. We therefore conclude that those surviving solitons formed by MI in a single quench are quasi-stationary within our experiment time $<200~{}$ms as their norm are sufficiently close to the Townes threshold.

A.4.2 Fitting atom number loss during rapid soliton collapse

Following the second interaction quench in path (ii) as shown in Fig. S4, mean atom number in a soliton $\bar{N}_{a}(\tau_{2}=0)$ exceeds the new Townes threshold by three-fold and solitons begin to collapse. From the linear fit in Fig. S4(b), we obtain an approximate linear time-dependent size $\sigma(\tau_{2})=\sigma_{0}+\dot{\sigma}\tau_{2}$, where $\sigma_{0}=8.0~{}\mu$m and $\dot{\sigma}=-0.11~{}$mm/s. Here we consider a simple case where the atom number loss within a soliton is fully due to three-body recombination. Because the vertical oscillator length $l_{z}=208~{}$nm far exceeds the magnitude of the 3D scattering length $|a|=1.1~{}$nm, we expect the three-body loss behavior to be 3D in nature. To quantify the three-body loss rate, we develop a model that takes into account the shrinking soliton size in the 2D plane, while assuming that the vertical wave packet remains in the harmonic ground state. This is justified because the 2D density $n$ needs to increase by $\sim 170$ times for the interaction energy to approach the vertical trap vibrational energy. Such dramatic increase in 2D density is not observed in our images, where the soliton size remains larger than the image resolution during the time of collapse. We therefore adopt the standard three-body recombination loss model

[TABLE]

where $L_{3}$ is the three-body loss coefficient and we have used a Gaussian form to approximate the soliton 3D density profile. From the above equations and the approximate linear time-dependence in $\sigma$, we derive an analytical formula

[TABLE]

where $L_{3}$ is the fit parameter. This effective three-body loss model appears to fit our data except beyond $\tau_{2}>36~{}$ms when the atom number reaches the new Townes threshold. The fitted $L_{3}=1.1(1)\times 10^{-23}~{}$cm${}^{6}/$s is nonetheless five orders of magnitude larger than that measured in a thermal sample [3], which we believe is not physical. To reconcile this discrepancy may require a different collapse dynamics that, for example, creates three-body loss under higher local collapse density beyond our image resolution and also rapidly ejects atoms out of a collapsing soliton [12, 13, 14, 15].