Solitary waves in a quantum droplet-bearing system

We unravel the existence and stability properties of dark soliton solutions as they extend from the regime of trapped quantum droplets towards the Thomas-Fermi limit in homonuclear symmetric Bose mixtures. Leveraging a phase-plane analysis, we identify the regimes of existence of different types of quantum droplets and subsequently examine the possibility of black and gray solitons and kink-type structures in this system. Moreover, we employ the Landau dynamics approach to extract an analytical estimate of the oscillation frequency of a single dark soliton in the relevant extended Gross–Pitaevskii model. Within this framework, we also find that the single soliton immersed in a droplet is stable, while multisoliton configurations exhibit parametric windows of oscillatory instabilities. Our results pave the way for studying dynamical features of nonlinear multisoliton excitations in a droplet environment in contemporary experimental settings.

Figure 1

Schematic representation of the underlying bifurcation diagram of the particle number $N\left(\mu \right)$ as a function of the chemical potential $\mu$ for the droplet ground state as well as the single-, two-, and three-droplet dark solitons for $\delta =1$. Vertical dashed lines indicate the turning point ${\mu }_{\mathrm{cr}}$ and the corresponding linear limit ${\mu }_{\mathrm{lin}}=(n+1/2)\mathrm{\Omega }$ of the adimensionalized model, where $\mathrm{\Omega }=0.1$ denotes the trap frequency. The solution (Hermite-Gauss state) departs from the linear limit and for decreasing chemical potential acquires a droplet background since the LHY dominates until it reaches ${\mu }_{\mathrm{cr}}$. Afterward, the configuration gradually deforms from a Gaussian-type backdrop towards a TF-type backdrop bearing one or more strongly localized dark soliton(s) for increasing $\mu$ where the cubic nonlinearity prevails. Characteristic density profiles of the ensuing localized waveforms at specific chemical potentials (see ellipse markers) are also provided. The latter correspond to $\mu$ $=-0.175$ and 0.115 for the ground states, $\mu$ $=-0.06$ and 0.2 for the single, $\mu$ $=0.02$ and 0.25 for the two, and $\mu$ $=0.17$ and 0.3 for the three dark soliton solutions, respectively. Note that for visualization purposes the curves are slightly shifted upwards from $N=0$. Namely, all the lower curves asymptote to $N=0$ as $\mu$ $\to \infty$.

Figure 2

Potential $V\left(q\right)$ (top panels), respective phase planes (middle panels), and bounded orbits (bottom panels) for $C_1=0$ and $\delta =1$ for which ${\mu }_{★}=-2/9$. The different columns of panels correspond to representative cases, for the values of $\mu$ indicated, for all the qualitative different scenarios that produce bounded orbits. Bounded (unbounded) orbits are depicted with solid (dashed) lines in the middle panels. For clarity of exposition, when symmetric ($q\leftrightarrow -q$) solutions exist, only positive bounded orbits are shown in the lower panels. Case (a) corresponds to $-{\delta }^2/4<\mu <{\mu }_{★}$ which gives rise to homoclinic orbits supported by a nonzero background. Case (b) corresponds to $\mu$ $={\mu }_{★}$ where the two homoclinic orbits of case (a) collide and give rise to heteroclinic orbits connecting zero and nonzero backgrounds. Case (c) corresponds to ${\mu }_{★}<\mu <0$ where the heteroclinic orbits of case (b) merge into heteroclinic orbits connecting nonzero backgrounds. Case (d) corresponds to $\mu >0$ where the homoclinic orbits of case (c) disappear. Note that, in order to better relate our results to experimental conditions, we opt to display in this figure (as it is the case for all subsequent figures) all of the relevant quantities in full dimensional units as indicated in the labels.

Figure 3

“Zoo” of nonlinear waves in the crossover from the droplet to the TF regime. Top and bottom sets of panels depict, respectively, the wave function, $q$, and the relevant density profiles, ${\left|q\right|}^2$. (a) Bubble wave function for $\mu$ $=-0.23$ referring to a homoclinic solution anchored at $q_{+}$. Bubbles are found within the interval $-{\delta }^2/4<\mu <{\mu }_{★}$ where $q_{+}$ exists and $V\left(0\right)>V(\pm q_{+})$, see also Fig. 2. (b) Kink solution corresponding to a heteroclinic connection between 0 and $q_{+}$ (or $-q_{+}$), occurring solely for $\mu$ $={\mu }_{★}=-2{\delta }^2/9$. Here, the maxima of the effective potential have the same height (energy), i.e., $V\left(0\right)=V(\pm q_{+})$ [see Fig. 2]. (c) Bright droplet for $\mu$ $=-0.2$. These states form for ${\mu }_{★}<\mu <0$ when $V\left(0\right)<V(\pm q_{+})$ and are homoclinic solutions at $q=0$ [see also Fig. 2]. (d) [(e)] Dark soliton configuration for $\mu$ $=-0.222$ $(\mu =0)$ being a heteroclinic orbit connecting $-q_{+}$ and $q_{+}$. These structures arise when $\mu >{\mu }_{★}$ and, as $\mu$ $\to {\mu }_{★}^{+}$, they become wider due to the bottleneck induced by the potential maximum at $q=0$ as evinced in panel (d). In all cases, $C_1=0$ and $\delta =1$.

Figure 4

(a), (c) Bifurcation diagram of the particle number $N$ for the single dark soliton solutions as a function of the chemical potential $\mu$ for different values of the LHY correction as indicated. (b), (d) Corresponding dependence of the internal (anomalous) mode of eigenfrequency ${\mathrm{\Omega }}_{\mathrm{osc}}/\mathrm{\Omega }$ vs $\mu$. Notice that, in addition to the dipolar mode of frequency ${\mathrm{\Omega }}_{\mathrm{osc}}/\mathrm{\Omega }=1$ (see the horizontal black line), the anomalous mode starts from the linear limit ${\mu }_{\mathrm{lin}}=(3/2)\mathrm{\Omega }$. As $N$ is increased, for small enough $\delta$ values (see the cases of $\delta =0.01$ and $\delta =0.1$), the anomalous mode decreases and then asymptotes to ${\mathrm{\Omega }}_{\mathrm{osc}}/\mathrm{\Omega }=1/\sqrt{2}$ for large $\mu$ values [see the horizontal dashed black line in panel (d)]. The inset in panel (d) corresponds to a magnification at small chemical potentials where also the difference from the prediction of Eq. (17) is maximal. In contrast, for large enough values of $\delta$ (see the cases of $\delta =0.5$ and $\delta =1$), the anomalous mode frequency first increases and then decreases as $N$ increases while the chemical potential also reaches a turning point for a minimal critical value of $\mu$ at ${\mu }_{\mathrm{cr}}$ (${\mu }_{\mathrm{cr}}=-0.085$ for $\delta =1$ and ${\mu }_{\mathrm{cr}}=0.091$ for $\delta =0.5$). The trapping strength for all cases is fixed at $\mathrm{\Omega }=0.1$.

Figure 5

Critical value of the chemical potential ${\mu }_{\mathrm{cr}}$ versus the LHY strength $\delta$ for (a) single dark soliton solutions and different trap frequencies $\mathrm{\Omega }$ (see legend) and (b) distinct dark soliton complexes ranging from one to three (see legend) and for $\mathrm{\Omega }=0.1$. Evidently, for fixed $\delta$ the turning point gets shifted to lower chemical potentials, either by decreasing the trap frequency or the soliton number. In all cases larger $\delta$ reduces ${\mu }_{\mathrm{cr}}$. The dashed horizontal line at ${\mu }_{\mathrm{cr}}=0$ is added as a guide to the eye.

Figure 6

Selected density profiles of the droplet dark soliton for (a) $\mu$ $=-0.15$ and (b) $\mu$ $=0.2$ and distinct $\mathrm{\Omega }$ values (see legends) in the case of $\delta =1$. Reducing the trap strength leads to a flat-top droplet dark configuration. Panel (c) and its inset depict different magnifications of the density profile of panel (a) for $\mathrm{\Omega }=0.0004$ illustrating the spatial resolution used.

Figure 7

Dependence of the anomalous mode, ${\mathrm{\Omega }}_{\mathrm{osc}}/\mathrm{\Omega }$, for (a) two and (b) three dark soliton solutions as a function of $\mu$ for $\delta =1$. The relevant turning points, being shifted to more positive $\mu$ values for higher soliton complexes, appear at ${\mu }_{\mathrm{c}}=0.014$ and ${\mu }_{\mathrm{cr}}=0.115$. In both cases the collision of the second anomalous mode (green and purple lines) with the background ones (black lines) signals the occurrence of an oscillatory instability.

Figure 8

Density evolution of anomalous modes ${|\mathrm{\Psi }(x,t)}_{{\mathrm{AM}}_i}{|}^2$. In each case, the steady state is perturbed with the corresponding anomalous modes initialized with an amplitude of $\epsilon =0.3$ and then allowed to evolve in time. Panels (a) and (c) [(b) and (d)] correspond to stable modes while panel (e) [(f)] pertains to the unstable mode for the two- [three-] soliton configuration. The anomalous modes in panels (a) and (b) represent in-phase oscillations while panels (c) and (d) correspond to out-of-phase oscillations. The unstable cases depicted in panels (e) and (f) for two and three solitons, respectively, showcase the destabilization of the second anomalous mode (${\mathrm{AM}}_2$) through a complex eigenfrequency quartet. It is observed that the dynamical evolution of this instability is manifested as amplified out-of-phase vibrations that eventually saturate and are subsequently subject to recurrence.

Figure 9

The potential $V\left(q\right)$ (top panels) and respective phase planes (middle panels) and orbits (bottom panels) for $C_1=0.1$ and $\delta =1$. Same layout and notation as in Fig. 2. For (a) $\mu$ $=-0.5$ and (b) $\mu$ $=-0.2018$, all trajectories are unbounded. For values of $\mu$ larger than case (b) ($\mu$ $=-0.2018$) there exist bounded periodic and homoclinic orbits. For example, for (c) $\mu$ $=-0.15$ and (d) $\mu$ $=0.5$ there exist homoclinic orbits corresponding, respectively, to a shallow and a deep (almost black) gray soliton solution; see the relevant solid red curves in each case.