Variational Approach to the Snake Instability of a Bose-Einstein Condensate Soliton

TL;DR

This paper develops a variational framework to analyze the snake instability of a dark soliton in a 2D BEC under anisotropic confinement. By coupling soliton bending with vortex-like perturbations along the nodal line in a carefully constructed ansatz, it derives an effective two-parameter dynamical system whose linearized form reveals growth rates, most unstable modes, and a critical trap anisotropy for stabilization. The TF background and phase-offset extensions yield quantitative predictions that closely match full GP simulations, demonstrating both stable oscillations and instability-induced vortex necklaces. The approach provides a compact, efficient analytical tool for predicting soliton stability in experimentally relevant trap geometries and suggests routes for generalizations to other nonlinear wave phenomena.

Abstract

Solitons are striking manifestations of nonlinearity, encountered in diverse physical systems such as water waves, nonlinear optics, and Bose-Einstein condensates (BECs). In BECs, dark solitons emerge as exact stationary solutions of the one-dimensional Gross-Pitaevskii equation. While they can be long-lived in elongated traps, their stability is compromised in higher dimensions due to the snake instability, which leads to the decay of the soliton into vortex structures among other excitations. We investigate the dynamics of a dark soliton in a Bose-Einstein condensate confined in an anisotropic harmonic trap. Using a variational ansatz that incorporates both the transverse bending of the soliton plane and the emergence of vortices along the nodal line, we derive equations of motion governing the soliton's evolution. This approach allows us to identify stable oscillation modes as well as the growth rates of the unstable perturbations. In particular, we determine the critical trap anisotropy required to suppress the snake instability. Our analytical predictions are in good agreement with full numerical simulations of the Gross-Pitaevskii equation.

Figures (9)

• Figure 1: The maximum density $n_\text{max}$ as a function of the anisotropy ratio $\omega$ for $gN=500$. Numerically obtained GP simulations (dots) follow the 2DTF (dashed) approximation for smaller anisotropies and the 1DTF (dash-dotted) in the large anisotropy limit. The insets display cross-sections of the density profile along the $y$-axis for selected values of $\omega$, demonstrating the crossover from the 2DTF to the 1DTF regime in large anisotropies.
• Figure 2: Time evolution of a dark soliton in a trapped Bose-Einstein condensate for two different trap anisotropies, both with $g = 1$ and $N = 400$. The colorbars show the density of BECs. Top row: For $\omega = 3$, the soliton is dynamically unstable and decays via the snake instability, leading to the formation of vortex structures. Bottom row: For $\omega = 35$, the soliton remains stable and persists without significant deformation. Note that the vertical axis ($y$) is scaled differently in the two panels to reflect the change in trap geometry.
• Figure 3: Density profiles of variational wavefunction defined in Eq. \ref{['eqn:ansatz-crossover']} with constant density $n_a$ for various choices of $\alpha$ and $\beta$. The colorbar shows the density of the wavefunction. The parameter $\alpha$ controls the nodal line breaking up into vortices, while $\beta$ determines the bending amplitude of the soliton plane.
• Figure 4: Density profiles of variational wavefunction defined in Eq. \ref{['eqn:ansatz-crossover']} inside a TF-profile background density for $\alpha=0.7$ and $\beta=0.4$. The wavevector $k$ controls the spacing between vortices, which is related to the number of vortices inside the condensate. We show two modes corresponding to (a) four vortices and (b) a single vortex as a decay product of snake instability inside an anisotropic condensate with $gN=400$ and $\omega=3$, as used in the first row of Fig. \ref{['fig:SI']}.
• Figure 5: Phase space of the variational dynamical system defined by Eqs.\ref{['eqn:const_density-eom']} in the complex $A=\alpha+\mathrm i\beta$ plane, illustrating the dependence of the flow topology on the perturbation wavenumber $k$. The colorbars indicate the potential energy $U(\alpha,\beta)$ calculated with Eq. \ref{['eqn:U']}, and the arrows show the direction of the flow in the phase space. For large $k=2.0\sqrt{\mu}$,(top), the origin is a stable center, and all nearby trajectories form closed orbits, indicating dynamically stable soliton configurations. For small $k=0.5\sqrt{\mu}$(bottom), the origin becomes a saddle point and two stable fixed points appear along the $\beta=0$ axis, signaling stable necklace of vortices and antivortices.