Quantum-Fluctuation-Driven Dynamics of Droplet Splashing, Recoiling and Deposition in Ultracold Binary Bose Gases

Given the wide practical relevance to both nature and industry [1, 2, 3, 4, 5, 6, 7], droplet impact dynamics on a surface has attracted much attention ever since the first study by Washington [8, 9]. Various impact outcomes-including splashing, receding/recoiling, rebound, and deposition-have been observed successfully in experiments [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. In general, these dynamics were characterized by two physical observables, namely, the maximum spreading factor ($\beta$) [12, 16, 23] and the splashing threshold(${\cal K}$) [11, 13, 14, 15, 17, 18, 19, 20, 21, 22, 25], which were shown not only to depend on the properties of the droplet itself (size, density, surface tension, viscosity, impact velocity), but also to be strongly influenced by the surface condition(roughness, wettability) and surrounding gas (pressure, composition). Because of the complexities associated with various unmanageable factors, it is extremely challenging to deterministically parametrize $\beta,\ {\cal K}$ and fully control the impact outcomes. In this situation, a common practice is to assume an ideal droplet impact (on a smooth solid surface at atmospheric condition) and then to quantify the actual dynamics by the Weber and Reynolds numbers, which, respectively, describe the relative strength of droplet inertia with respect to capillary and viscous forces [26, 27]. Various scaling laws between $\beta,\ {\cal K}$ and these numbers have been proposed in the literature [26, 27], based on different models or empirical fitting from experimental data.

In the past few decades, ultracold atoms have emerged as an ideal platform for quantum simulation, given their extremely clean environment and the high controllability on the species, number, dimension, interaction strength etc [28, 29]. In particular, a recent important achievement in this field was the realization of a quantum droplet in both dipolar gas [30, 31, 32, 33, 34, 35, 36] and alkali bosonic mixtures [37, 38, 39, 40, 41, 42], with extremely dilute densities ($\sim 10^{14}-10^{15}{\rm cm}^{-3}$) that can be eight orders of magnitude lower than water. In forming these gaseous droplets, quantum fluctuations play an essential role in providing the repulsive force for their stabilization, for which they are called quantum droplets [43]. To date, the idea of a quantum droplet has been successfully extended to various atomic systems, including low dimensional ones [44, 45, 46, 47, 48, 49], Bose-Fermi mixtures [50, 51, 52, 53, 54, 55] and multi-component dipolar or alkali atomic mixtures [56, 57, 58]. The non-equilibrium properties of quantum droplets have also been investigated in terms of their dynamical formations [59, 60, 61] and collisions [62, 63]. These developments offer an unprecedented opportunity for simulating droplet impact dynamics in ultracold atoms, particularly, at the microscopic quantum level and in a highly controllable manner.

In this work, we demonstrate the capability of using ultracold Bose gases to simulate the droplet impact outcomes in a fully controlled quantum setting. Contrary to the conventional droplet impact setup, here there is no impact surface for the droplet, and the driving force of its dynamics is purely from its intrinsic energy due to quantum fluctuations. Specifically, we study the dynamical property of a two-dimensional (2D) Bose-Bose mixture with repulsive intraspecies and attractive inter-species couplings, whose ground state is a quantum droplet. To highlight the quantum effect in the dynamics, we have chosen the initial state as the Townes soliton generalized from the single-species case [64], which features a zero mean-field energy with continuous scale invariance as recently confirmed in experiments [65, 66, 67, 68]. In this way, the dynamics here is purely driven by quantum fluctuations and thus can be called quantum-fluctuation-driven dynamics (QFDD). It is found that by tuning the fluctuation energy of the initial Townes state through its size($\sigma_{0}$) and number($N$), the subsequent QFDD can produce various outcomes, including splashing, recoiling and deposition, similar to those in droplet impact dynamics. We have mapped out the dynamical phase diagram in the ($\sigma_{0},N$) plane and employed the Weber number to characterize different phases. The splashing and recoiling thresholds are identified, and a universal scaling law is established between the maximum spreading factor and the Weber number for the recoiling dynamics, which is applicable for a considerably large parameter regime. Finally, we show that the long-time QFDD in the deposition regime can be used to extract the collective breathing modes of quantum droplet. These results can be directly tested in the current cold atoms experiments. In the Appendixes, we provide more details on the derivation of the generalized Townes soliton with unequal masses, as well as on the numerical simulations.

We start from the energy functional of two-species bosons in 2D ($\hbar=1$ for brevity):

here $\bm{\bm{\rho}}=(x,y)$ is the 2D coordinate; $\phi_{i}$ is the wavefunction of the $i$-th species and $n_{i}=|\phi_{i}|^{2}$ is its density; $g_{ij}$ is the bare coupling between $i$- and $j$-species, which can be expressed as $g_{ij}=4\pi a_{ij}/(ml_{z})$ in quasi-2D geometry, with $a_{ij}$ the s-wave scattering length and $l_{z}$ the characteristic length along the confined ($z$) direction; and ${\cal E}_{\rm LHY}$ is the Lee-Huang-Yang(LHY) correction from quantum fluctuations, and for quasi-2D bosons near the mean-field instability point ($\delta a\equiv a_{12}+\sqrt{a_{11}a_{22}}\sim 0$) it reads [44, 47]

It has been shown that this LHY term can balance with the mean-field force and result in a self-bound droplet as the ground state [44]. Given the $n^{2}\ln(n)$ dependence of LHY energy, the 2D droplet can exist at both mean-field collapse and stable regimes [44] and with any infinitesimal atom number [49]. At a sufficiently large number, the 2D droplet develops a flat-top structure in its density profile, similar to the 3D case [43].

Given Eq.(1), the dynamics of $\{\phi_{i}\}$ is governed by the time-dependent Gross-Pitaevskii(GP) equations:

In this work, we will focus on the solution with zero angular momentum, since the ground state and the Townes soliton both stay within this sector, and different angular momentum sectors are decoupled from each other. In this sector, we can replace the coordinate $\bm{\bm{\rho}}$ simply by its magnitude $\rho=|\bm{\rho}|$. More details of solving Eq.(3) in the discretized coordinate and time space have been given in the Appendixes.

Throughout the paper, we specifically consider the two hyperfine states of ³⁹K atoms, $|1\rangle\equiv|F=1,m_{F}=-1\rangle,\ |2\rangle\equiv|F=1,m_{F}=0\rangle$, as has been well studied in the droplet experiments [37, 38, 39]. In this case, $a_{11}=33.5a_{B},\ a_{12}=-53a_{B}$ ($a_{B}$ is the Bohr radius), and $a_{22}$ is highly tunable by magnetic field. The confinement length is chosen as $l_{z}=0.08\mu m$.

For spinless bosons in 2D, it is known that the kinetic term and mean-field attraction can support a special stationary solution called the Townes soliton [64]. This state features zero energy and continuous scale invariance, and it can only exist when the boson number and coupling strength satisfy $N|g|=5.85/m$. Such a special solution has been successfully observed in both non-linear optics [65] and ultracold atoms [66, 67, 68]. In these experiments, the LHY correction takes little effect as it is much smaller than the mean-field part. In the following, we will show that the Townes soliton can be generalized to two-species bosons if equally neglecting the LHY correction.

By omitting the LHY term in (3), we can see that the two GP equations for $\{\phi_{1},\phi_{2}\}$ can support a single-mode solution $\phi_{i}=\sqrt{N_{i}}\phi\exp(-i\mu_{i}t)$ as long as

where $\mu_{1}=\mu_{2}\equiv\mu$, and the single mode $\phi$ satisfies:

Here $N=N_{1}+N_{2}$ is the total number, and the effective interaction $g_{\rm eff}$ is given by

Apparently we have $g_{\rm eff}<0$ in the mean-field collapse regime ($g_{12}<-\sqrt{g_{11}g_{22}}$). It is then straightforward to check that under the condition

there exists a sequence of zero-energy eigenstates with continuous scale invariance, i.e., the eigenstate nature and zero-energy property will not change under an arbitrary scaling transformation $\phi(\bm{\rho})\rightarrow\lambda\phi(\lambda\bm{\rho})$ (accordingly $E\rightarrow\lambda^{2}E$). These stationary solutions are the generalized Townes soliton for two-species bosons. Note that a similar generalization also works for the case of unequal masses, with slight modifications in Eqs.(4,6) as shown in the Appendixes.

In figure 1(a), we confirm the stationary Townes soliton for two-species bosons once the total number $N$ satisfies (7) and the number ratio $N_{1}/N_{2}$ satisfies (4). In comparison, if we change $N$ to be smaller or larger, the original profile will shrink (figure 1(b)) or expand (figure 1(c)) with time. The profile is also unstable if $N_{1}/N_{2}$ deviates from (4), see figure 1(d). In a word, both conditions (7) and (4) are required in supporting a stationary two-species Townes solution.

A crucial difference between the single- and two-species bosons is that quantum fluctuations play an important role in the latter, which can lead to droplet formation as a ground state. It then follows that starting from the generalized Townes soliton, which is the mean-field stationary solution for two-species bosons, the quantum fluctuation can destabilize it strongly and drive its time evolution towards the droplet formation. Such dynamics can be called the quantum-fluctuation driven dynamics (QFDD), also in light of the fact that the total energy of Townes soliton is purely given by the LHY part, $E(t=0)=E_{\rm LHY}$.

In figure 2(a), we show that the total energy ($=E_{\rm LHY}$) of two-species Townes soliton can be conveniently tuned by its size $\sigma$, taking a typical combination of $\{N,g_{ij}\}$ that satisfies (4,7). In particular, $E_{\rm LHY}$ varies non-monotonically with $\sigma$ and shows a minimum at certain finite $\sigma$. To understand this behavior, we employ a Gaussian ansatz to approximate the single mode $\phi(\rho)=\exp[-\rho^{2}/(2\sigma^{2})]/(\sqrt{\pi}\sigma)$, which leads to the total energy $E=\int 2\pi\rho d\rho{\cal E}(\rho)$ as

with $\bar{a}\equiv\sqrt{a_{11}a_{22}}$. We can see that the first two terms can support zero-energy states with arbitrary $\sigma$ under the condition $Nm|g_{\rm eff}|=2\pi$, which are just the simplified Gaussian version of the Townes soliton. However, when including the third LHY term, the total energy($=E_{\rm LHY}$) will deviate from zero, and the Townes profile is no longer stationary. Given the expression $E_{\rm LHY}\sim-\ln(\sigma^{2})/\sigma^{2}$, we can easily arrive at a non-monotonic $E_{\rm LHY}\sim\sigma$ dependence with energy minimum $E_{\rm min}=-N\bar{a}/(4l_{z}^{3}e)$ at $\sigma_{\rm min}=\sqrt{4N\bar{a}l_{z}e}$. As shown in figure 2(a), $E_{\rm LHY}$ from a Gaussian ansatz provides a qualitatively good prediction to the $E\sim\sigma$ lineshape of real Townes solutions.

Given the easily tunable fluctuation energy of the initial Townes state, the subsequent QFDD can exhibit rich dynamical outcomes. In figure 2(b1,b2,b3), we show the time evolution of density profiles for three typical QFDDs:

(I) Deposition. When the initial size $\sigma_{0}$ is close to $\sigma_{\rm min}$ and $E_{\rm LHY}$ is small, the QFDD shows a typical deposition behavior (figure 2(b1)). Specifically, as time passes the system repels a small proportion of atoms outside, and the rest automatically follows the profile of a ground-state droplet with additional periodic oscillations. As discussed later, such residual oscillation can be used to probe the collective breathing modes of a quantum droplet.

(II) Recoiling. As $\sigma_{0}$ deviates more from $\sigma_{\rm min}$ and $E_{\rm LHY}$ gets larger, the system enters the recoiling regime. As shown in figure 2(b2), at early times a considerable portion of atoms are repelled outside, while at some point they stop to spread and flow back to merge with the central part. Such back-flow (or recoiling) can be attributed to the competition between the surface tension and the kinetic energy of the cloud during the dynamics.

(III) Splashing. When $\sigma_{0}$ deviates significantly from $\sigma_{\rm min}$ and $E_{\rm LHY}$ is large enough, the system shows a rapid splashing dynamics, see figure 2(b3). In this case, the large $E_{\rm LHY}$ causes a drastic change of the initial Townes profile in a short time, i.e., the cloud quickly splits into two pieces. In this process, $E_{\rm LHY}$ converts to the large kinetic energy of the outgoing part, such that it completely separates from the central part and flows away forever.

The above dynamics can be well distinguished by monitoring the mean size of the dynamical system, $\sigma(t)=\sqrt{\langle\bm{\rho}^{2}\rangle_{t}}$. As shown in figure 2(c), at longer times $\sigma(t)$ is almost static in the deposition regime, while it shows a non-monotonic behavior in the recoiling regime and a continuous increase in the splashing regime.

To this end, we have shown that the QFDD starting from the Townes states can produce rich dynamical phases, including the deposition, recoiling, and splashing, which perfectly mimic the droplet impact outcomes as studied in the literature [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. However, different from these existing studies of droplet impact on a surface, in our case the driving force of the dynamics is purely from the intrinsic energy contributed by quantum fluctuations. This introduces a mechanism for these fluid dynamics. Meanwhile, in QFDD there are no complexities caused by the impact surface or environment, and the dynamical outcome can be fully controlled by adjusting the size $(\sigma_{0})$ and the number $(N)$ of the initial state. Therefore, the ultracold atoms provide an extremely clean and convenient platform to simulate droplet dynamics, where the quantum effect can be well manipulated and the fluid mechanics can be understood in a more deterministic way.

We now quantify various dynamical phases in QFDD by the Weber number. Note that the Reynolds number is irrelevant here because of the zero viscosity of Bose condensates. In conventional droplet impact dynamics [26, 27], the Weber number is defined as $W=\rho_{d}Dv^{2}/\gamma$, where $\rho_{d},\ D,\ v,\ \gamma$ respectively denote the droplet density, diameter, impact velocity, and surface tension. It measures the relative strength of droplet inertia with respect to capillary force. Here, we generalize the definition of $W$ to describe the quantum dynamics in general:

where $\Delta E$ is the energy difference between the initial state (here the Townes soliton) and the true ground state given the same initial parameters ($N,\ g_{ij}$); $D$ and $\gamma$ are the same as before, i.e., the droplet diameter and surface tension. Specifically, we have $\gamma=\int d\bm{\bm{\rho}}\left({\cal E}(\bm{\bm{\rho}})-\mu_{1}|\phi_{1}|^{2}-\mu_{2}|\phi_{2}|^{2}\right)$, with $\mu_{i}$ the chemical potential for the i-th species.

As shown in the color plot in figure 2(a), $W$ defined in (9) can well characterize different dynamical phases in QFDD. Namely, a low $W$ corresponds to the deposition dynamics, where the small $\Delta E$ can be well absorbed by the surface change of the droplet; as $W$ increases, the system enters the recoiling regime and finally end up at splashing, where the large $\Delta E$ overwhelms the capacity of the droplet surface and causes it to change drastically. Hereafter, we refer to the critical $W$ at the recoiling-splashing boundary as the splashing threshold (${\cal K}_{s}$), and that at the deposition-recoiling boundary as the recoiling threshold (${\cal K}_{r}$).

In figure 3(a), we map out the dynamical phase diagram in the ($\sigma_{0},N$) parameter plane. Distinct dynamical outcomes of QFDD, including deposition (’D’), recoiling (’R’) and splashing (’S’), are identified by monitoring the mean size $\sigma(t)$ of the cloud during expansion (see figure 2(c)). In addition, we show the contour plot of $W$ in the ($\sigma_{0},N$) plane, and one can see clearly that the ’D’, ’R’ and ’S’ phases respectively correspond to the small, intermediate and large $W$ regions. Due to the non-monotonic dependence of $W$ on $\sigma_{0}$ (as shown in figure 2(a)), there are two recoiling regions in the diagram, as marked by ’R₁’ and ’R₂’. In figure 3(b), we extract the two recoiling thresholds(${\cal K}_{r1},{\cal K}_{r2}$) and the splashing threshold(${\cal K}_{s}$) along the phase boundaries as varying $N$. One can see that ${\cal K}_{r1}\approx{\cal K}_{r2}$ are given by a constant $\sim 0.12$, while ${\cal K}_{s}$ is a much larger value and continuously increases with $N$.

Another important physical quantity to characterize the droplet impact dynamics is the maximum spreading factor $\beta=\sigma_{\rm max}/\sigma_{0}$, as defined by the ratio between the maximum spreading radius ($\sigma_{\rm max}$) and the initial one ($\sigma_{0}$). Clearly one has $\beta\sim 1$ for the deposition dynamics and $\beta\rightarrow\infty$ for splashing. An interesting behavior of $\beta$ shows up in the recoiling regime, where $\beta$ is finite and varies sensitively with $W$.

In figure 4, we extract $\beta$ as a function of $W$ along the horizontal lines in figure 3(a), i.e., by varying $\sigma_{0}$ at several fixed $N$ in the recoiling regime. Apart from the region near the phase boundaries ($\beta\rightarrow 1$ or large enough), the data of $\beta$ for any given $N$ well follow the scaling relation:

As shown in figure 4, such scaling works well for a wide parameter regime with $N/10^{4}\in(0.5,2)$, $\beta\in(1.5,10)$ and $W\in(0.2,1)$. Therefore, we expect the scaling law in (10) to reflect a very robust intrinsic property of the ultracold fluid during the recoiling QFDD.

Finally, we demonstrate that the deposition regime of QFDD can be used to probe the breathing modes of quantum droplet. The breathing mode can be theoretically obtained as follows. Assuming a small fluctuation mode $\delta\phi_{i}$ for the i-th species boson, and only keeping the lowest-order fluctuations in the GP equations(3), we obtain the following equations for $\{\delta\phi_{i}\}$:

According to the standard Bogoliubov analysis, we search for solutions of the form:

Here $\omega_{j}$ is the $j$-th collective (eigen-)mode of the system. These modes can be extracted from the following coupled equations for $u_{ij}(\bm{\rho})$ and $v_{ij}(\bm{\rho})$:

where

In Fig.5, we show the lowest two breathing modes ($\omega_{1},\omega_{2}$) below the atom emission threshold ($-\mu_{1},-\mu_{2}$), see solid lines. In our numerics, we have obtained $\omega_{j}$ by discretizing the coordinate space and exactly diagonalizing the resulted large matrix in the left-hand-side of (14). Note that the breathing modes, by definition, stay in the same zero angular momentum sector as the Townes state and the ground state droplet. More details on the numerics are given in the Appendixes.

In figure 5, we show how to extract collective breathing modes from the QFDD in the deposition regime for a given set of coupling strengths $\{a_{ij}\}$ and initial atom number $N$. Figure 5(a) shows the typical density oscillations with time at different locations. One can see that such oscillations are synchronic for both species and at different locations, and therefore their Fourier transformations give the same peak frequencies; see figure 5(b). In figure 5(c), we plot out two sets of extracted peak frequencies ($\omega_{1},\omega_{2}$) with different residual droplet numbers $N_{d}$. Here $N_{d}$ can be tuned effectively by setting different sizes of the initial Townes state. For all the extracted data we have collected 30 samples by varying the real-space locations or the time intervals in Fourier transformation, from which we obtain both the mean value and the variance. One can see that the extracted results match very well with theoretical predictions of collective breathing modes from Eq.14 (solid lines).

We note that a previous study has extracted one branch of a collective mode from the formation dynamics of a quantum droplet [59], while the dynamics there is not QFDD. Here, we have shown in figure 5 that two branches of collective modes can be simultaneously extracted from the residual oscillations in QFDD. Moreover, we have checked that by taking different $\{a_{ij}\}$ and initial $N$, one can also get access to the other regime with only a single breathing mode. In comparison to the single-mode case, here the main complexity brought by the presence of two modes is the appearance of additional frequency peaks at $\omega_{2}\pm\omega_{1}$ due to their interference. To correctly identify the breathing modes $\{\omega_{j}\}$ from the multiple frequency peaks, it is important to note that these modes are associated with the most pronounced peaks among all others. This can be understood by analyzing the spin densities $n_{i}=|\phi_{i}+\delta\phi_{i}|^{2}$, with $\delta\phi_{i}$ given in Eq.(13). Obviously, one can see that the Fourier transform of $n_{i}$ results in multiple peaks at $\{\omega_{j}\}$, $\{2\omega_{j}\}$ and $\{\omega_{i}\pm\omega_{j}\}$, while the ones at $\{\omega_{j}\}$ have the largest magnitudes, which depend linearly on the strength of the fluctuation mode (the others all have a quadratic dependence). Based on this principle, one can easily identify the correct breathing modes ($\omega_{1},\omega_{2}$) as marked in figure 5(b).

Our results can be readily tested in the current cold-atom experiments. The initial Townes profile with a given amplitude and size can be imprinted by properly designing the optical potential applied to the atoms, as successfully implemented in previous experiments on two-component Bose gases [68, 69]. The subsequent dynamics of the system can then be measured through the in situ image, and various dynamical outcomes can be distinguished typically within tends of milliseconds, as shown in Fig.2. Within this time scale severe atom loss in the ³⁹K droplet can be effectively avoided [37]. On the other hand, for long time dynamics, the atom loss will play an essential role and may even impede the breathing mode detection, which requires a typical time scale of hundreds of milliseconds, as shown in Fig.5(a1,a2). For that, one can resort to the heteronuclear mixtures of ⁴¹K and ⁸⁷Rb, whose lifetime has been shown to be much longer ($\sim 1s$) due to the low droplet density therein to suppress the atom loss [40, 41]. Since the Townes soliton can be well extended to boson mixtures with unequal masses (see the Appendixes), and the controllability of LHY energy by the size and number generally applies to these systems, we expect that various dynamical outcomes and the breathing mode extraction can be tested equally in heteronuclear Bose gases.

In summary, we have demonstrated the quantum simulation of droplet impact outcomes in ultracold boson mixtures. Various dynamical phases including splashing, recoiling, and deposition, have been revealed. A remarkable difference here is that these dynamics are purely driven by quantum fluctuations, instead of the mechanical impact force in previous studies. Given the easy manipulation of initial states with tunable fluctuation energies, the current cold-atom platform provides complete control of these dynamics in the microscopic quantum level. To characterize different dynamics, we have introduced the Weber number and examined two important physical quantities, namely, the splashing/recoiling threshold and the maximum spreading factor. We have also proposed to extract the collective breathing modes from the residual dynamics in the deposition regime. These results are directly relevant to ongoing cold-atom experiments.

Furthermore, we remark that our work is in distinct contrast to the previous studies of QFDD in ultracold atoms, where a visible dynamics can only be achieved for small condensates [70, 71, 72, 73]. This is because usually the energy difference per particle between the mean-field and the true quantum ground states decays rapidly as the number $N$ increases [74, 75, 76, 77] , and thus for large $N$ the quantum fluctuation energy is negligibly small as compared to the mean-field one. However, this is not the case for a quantum droplet. For a static droplet, quantum fluctuations have been shown to provide an indispensable repulsive force for its stabilization, without which the whole system will collapse [43]. In parallel, our current work reveals the equally significant role of quantum fluctuations played in the non-equilibrium dynamics of a quantum dropet, regardless of whether the droplet is small or large.

Finally, we anticipate that the dynamical outcomes of splashing, recoiling, and deposition revealed in this work may be equivalent with other initial conditions when the dynamics is not purely driven by quantum fluctuations. In fact, the generalized Weber number as defined in Eq.(9) can be applied to any initial state, and the key issue is whether the additional energy of such a state (as compared to the true ground state) can be absorbed by the surface tension of a droplet during dynamics. In this sense, Eq.(9) may serve as a unified quantity to understand the different dynamics of quantum droplets with various initial conditions.

The work is supported by the National Key Research and Development Program of China (2018YFA0307600), the National Natural Science Foundation of China (12074419, 12134015, 12205365), and the Strategic Priority Research Program of Chinese Academy of Sciences (No. XDB33000000).