Numerical Identification of Stationary States and Their Stability in a Model of Quantum Droplets

TL;DR

This work studies stationary states and their stability in an extended NLS model for quantum droplets in Bose mixtures with Lee–Huang–Yang corrections, across 1D and 2D settings. It develops and deploys three robust numerical strategies—a companion-based multi-level method, a homotopy grid expansion, and a dimension-by-dimension homotopy continuation—to systematically uncover diverse stationary states and complex bifurcation structures arising from competing nonlinearities. The results reveal unprecedented bifurcation phenomena, including continuous connections between vortex and dark soliton stripe branches and several pitchfork and saddle-center events, enriching the landscape relative to the cubic defocusing NLS. These insights highlight the potential for similar complex bifurcations in higher dimensions or in models with competing nonlinearities, and they set a framework for future experimental and theoretical exploration of quantum droplets in Bose mixtures.

Abstract

In this work, we are motivated by a recent variant of the nonlinear Schrodinger (NLS) equation describing cold, dilute atomic condensates with quantum fluctuation effects. Our goal is to develop robust numerical methods capable of uncovering diverse stationary solutions in such NLS models. Specifically, and in line with recent theoretical and experimental interest, we focus on ultracold quantum droplets in Bose mixtures influenced by the Lee Huang Yang quantum fluctuation correction and study these systems in one and two dimensional settings. To this end, we deploy several numerical techniques. The homotopy grid method allows systematic refinement from coarse to fine spatial discretizations in one dimension, while the dimension by dimension homotopy approach extends one-dimensional solutions to two-dimensional domains. These methods effectively detect broad families of stationary states, many of which have not been previously reported, to the best of our knowledge. Furthermore, they enable the monitoring of solution continuation and bifurcation phenomena. During our investigation, we encounter unusual bifurcation events, including nonstandard pitchforks and saddle-center bifurcations, which exhibit novel stability transitions. For example, we identify continuous pathways connecting vortex and dark soliton stripe branches, absent in the standard cubic defocusing model. Overall, the presence of competing mean-field and quantum fluctuation interactions leads to a richer bifurcation structure than in traditional cubic NLS systems. These findings suggest that similar complex bifurcation and stability phenomena may appear in other settings, including higher-dimensional systems or models with competing nonlinearities such as cubic-quintic interactions, highlighting the importance of further theoretical and numerical exploration.

Figures (19)

• Figure 1: Overview of the numerical methods employed for both 1D and 2D systems. For 1D systems, we first use the companion-based multi-level method and the homotopy grid expansion method to compute multiple solutions. For 2D systems, the solutions obtained from 1D serve as the initial guess, and we then apply the dimension-by-dimension homotopy method and the homotopy grid expansion method to compute the corresponding 2D solutions.
• Figure 2: Illustration of the homotopy grid expansion method. Starting from coarse-grid solutions $\tilde{\bm \Psi}^{\,l}_{\mathbb{R}}$ and ${\bm \Psi}^{\,l}_{\mathbb{R}}$, the solutions are combined and the homotopy system is tracked from $s=0$ to $s=1$ to obtain the fine-grid solution ${\bm \Psi}^{\,l+1}_{\mathbb{R}}$. Here, "B" denotes boundary points.
• Figure 3: A total of 17 numerical solutions of Eq. \ref{['1d_re']} computed on 1025 grid points over the domain $D = [-12,12]$ with parameters $\mu = 0.8$ and $\Omega = 0.2$, under homogeneous Dirichlet boundary conditions.
• Figure 4: Based on solutions in Fig. \ref{['1d_2']}, we employed arclength continuation to trace the different solution branches, including those that converge to the trivial state for $\mu = 0.8$. Notably, the first branch encountered corresponds to the ground state, while the second one yields the first excited state, the third gives the second excited state, and so on.
• Figure 5: The first figure visualizes the bifurcation with respect to $\mu$. The computational domain is chosen as $D = (-12,12)^2$ with $\Omega = 0.2$. The discretization uses $N_x = 129$ and $N_y = 129$ grid points. In the lower panels, the left figure shows the density $|\psi|^2 = \psi_{\mathbb{R}}^2 + \psi_{\mathbb{C}}^2$, while the right figure displays the corresponding phase angle. The phase angle is defined by $\theta = \tan^{-1}(\psi_{\mathbb{C}} / \psi_{\mathbb{R}})$.