Snake instability of dark solitons in fermionic superfluids

We present numerical calculations of the snake instability in a Fermi superfluid within the Bogoliubov-de Gennes theory of the Bose-Einstein condensate (BEC) to BCS crossover using the random-phase approximation complemented by time-dependent simulations. We examine the snaking behavior across the crossover and quantify the time scale and length scale of the instability. While the dynamics shows extensive snaking before eventually producing vortices and sound on the BEC side of the crossover, the snaking dynamics is preempted by decay into sound due to pair breaking in the deep BCS regime. At the unitarity limit, hydrodynamic arguments allow us to link the rate of snaking to the experimentally observable ratio of inertial to physical mass of the soliton. In this limit we witness an unresolved discrepancy between our numerical estimates for the critical wave number of suppression of the snake instability and recent experimental observations with an ultracold Fermi gas.

Figure 1

Rate of snake instability of a dark soliton as a function of the wave vector of the snaking oscillation. Stars joined by solid lines correspond to the resonant poles calculated with RPA for $1/\left(k_{F}a\right)$ equal to (a) 0.2, (b) 0, (c) $-$0.5, and (d) $-$0.75. The numerical error for these points is about $0.001E_F$. Triangles correspond to the growth rate obtained in the time-dependent BdG simulations. The solid straight line in each panel is the hydrodynamic prediction (4), valid in the small-$q$ limit, with $E_s$ and $m_s$ obtained from the stationary BdG equations. At unitarity [panel (b)] this line coincides with Eq. (8), while the dashed line represents Eq. (9), which is the same hydrodynamic relation, but using the experimental value for the chemical potential and the period of oscillation of solitons as measured in [8].

Figure 2

Decay of the soliton for various values of the interaction strength $1/k_{F}a$ in the BCS regime. The column on the left depicts the evolution of the soliton at unitarity, the middle column is about $1/k_{F}a=-0.5$, while the rightmost column describes the case $1/k_{F}a=-1$. The time evolution at unitarity is presented at the times $0\hslash /E_F$ (a), $34\hslash /E_F$ (d), $43\hslash /E_F$ (g), and $54\hslash /E_F$ (j); for the middle column the times are $0\hslash /E_F$ (b), $67\hslash /E_F$ (e), $79\hslash /E_F$ (h), and $90\hslash /E_F$ (k), while at $1/k_{F}a=-1$ the evolution is shown for $0\hslash /E_F$ (c), $157\hslash /E_F$ (f), $180\hslash /E_F$ (i), and $200\hslash /E_F$ (l). The horizontal arrow in between the first and second row indicates the length scale in both $x$ and $y$ directions.

Figure 3

Characteristic lengths in the BCS regime (${\xi }_c$) and BEC regime (${\xi }_h$) as a function of $1/\left(k_{F}a\right)$, as computed with Eqs. (B4) and (B1). The values of ${\xi }_{\max }^{\mathrm{time}}$ and ${\xi }_{\max }^{\mathrm{RPA}}$ are taken from the time-dependent approach by looking at the longest perturbation wavelength that gives an imaginary spectrum.