Crossover between snake instability and Josephson instability of dark solitons in superfluid Fermi gases

TL;DR

This paper investigates the decay mechanisms of dark solitons in superfluid Fermi gases, revealing a crossover from snake instability in the BEC regime to Josephson instability in the BCS regime, with distinct underlying physical processes.

Contribution

It demonstrates a qualitative difference between soliton decay in BEC and BCS regimes using simulations and analysis, highlighting a crossover in decay mechanisms.

Findings

In BEC regime, amplitude fluctuations induce snaking deformations.

In BCS regime, phase fluctuations lead to Josephson currents destroying the soliton.

The decay mechanisms are experimentally distinguishable.

Abstract

Dark solitons in superfluid Bose gases decay through the snake instability mechanism, unless they are strongly confined. Recent experiments in superfluid Fermi gases have also interpreted soliton decay via this mechanism. However, we show using both an effective field numerical simulation and a perturbative analysis that there is a qualitative difference between soliton decay in the BEC- and BCS-regimes. On the BEC-side of the interaction domain, the characteristic snaking deformations are induced by fluctuations of the amplitude of the order parameter, while on the BCS-side, fluctuations of the phase destroy the soliton core through the formation of local Josephson currents. The latter mechanism is qualitatively different from the snaking instability and this difference should be experimentally detectable.

Equations

$$S[\Psi ]={\int }_0^{\beta }d\tau \int d\mathbf{r}$$

$$+Q\frac{\partial {\Psi }^{\ast }}{\partial \tau }\frac{\partial \Psi }{\partial \tau }-R{\left(\frac{\partial |\Psi {|}^2}{\partial \tau }\right)}^2]$$

$$i\tilde{D}(|\Psi {|}^2)\frac{\partial \Psi }{\partial t}=-C{\nabla }_{\mathbf{r}}^2\Psi +Q\frac{{\partial }^2\Psi }{\partial t^2}+\left(\mathcal{A}(|\Psi {|}^2)+2E{\nabla }_{\mathbf{r}}^2|\Psi {|}^2-2R\frac{{\partial }^2|\Psi {|}^2}{\partial t^2}\right)\Psi$$

$$\tilde{D}=\frac{\partial \left(|\Psi {|}^{2}D\right)}{\partial \left(|\Psi {|}^2\right)}\mathcal{A}$$

$$-C{\partial }_x^2{\Psi }_s+\left(\mathcal{A}(|{\Psi }_s{|}^2)+2E{\partial }_x^2|{\Psi }_s{|}^2\right){\Psi }_s=0$$

$$\underset{x\to \pm \infty }{\lim }\Psi (x)=\mp {\Psi }_{\infty }$$

$$P_{\pm }(\mathbf{r},t)=\left[\delta \Psi (\mathbf{r},t)\pm \delta {\Psi }^{\ast }(\mathbf{r},t)\right]/2$$

$$\displaystyle\alpha_{1}(x)\,P_{+}^{{}^{\prime\prime}}+\alpha_{2}(x)\,P_{+}^{{}^{\prime}}+\Big{(}\alpha_{3}(x)-\omega^{2}\,\alpha_{4}(x)\Big{)}\,P_{+}=\,\omega\,\gamma(x)\,P_{-}$$

$$\displaystyle\beta_{1}\,P_{-}^{{}^{\prime\prime}}+\Big{(}\beta_{2}(x)-\omega^{2}\,\beta_{3}\Big{)}\,P_{-}=\omega\,\gamma(x)\,P_{+}$$

$$N_{\pm }=\int |P_{\pm }(x){|}^{2}dx$$

$$\eta =\frac{N_{+}}{N_{+}+N_{-}}$$

$$S[\psi ]={\int }_0^{\beta }d\tau \int d\text{x}\left[\sum _{\sigma \in \{\uparrow ,\downarrow \}}{\bar{\psi }}_{\sigma }(\text{x},\tau )\left(\frac{\partial }{\partial \tau }-{\nabla }_{\text{x}}^2-{\mu }_{\sigma }\right){\psi }_{\sigma }(\text{x},\tau )+g{\bar{\psi }}_{\uparrow }(\text{x},\tau ){\bar{\psi }}_{\downarrow }(\text{x},\tau ){\psi }_{\downarrow }(\text{x},\tau ){\psi }_{\uparrow }(\text{x},\tau )\right]$$

$$S_{\text{EFT}}[\Psi ]={\int }_0^{\beta }d\tau \int d\mathbf{r}$$

$$+Q\frac{\partial \bar{\Psi }}{\partial \tau }\frac{\partial \Psi }{\partial \tau }-R{\left(\frac{\partial |\Psi {|}^2}{\partial \tau }\right)}^2]$$

$${\Omega }_s$$

$$D$$

$$C$$

$$E$$

$$Q$$

$$R$$

$$+\frac{4({\xi }_{\mathbf{k}}^2-2E_{\mathbf{k}}^2)}{3}{f}_3(\beta ,E_{\mathbf{k}},\zeta )+2E_{\mathbf{k}}^2|\Psi {|}^2{f}_4(\beta ,E_{\mathbf{k}},\zeta )]$$

$$f_j(\beta ,ϵ,\zeta )=\frac{1}{\beta }\sum _n\frac{1}{{\left[{\left({\omega }_n-i\zeta \right)}^2+{ϵ}^2\right]}^j}$$

$$\frac{\partial {\Omega }_s}{\partial |\Psi {|}^2}\Psi =0$$

$$i\tilde{D}(|\Psi {|}^2)\frac{\partial \Psi }{\partial t}=-C{\nabla }_{\mathbf{r}}^2\Psi +Q\frac{{\partial }^2\Psi }{\partial t^2}+\left(\mathcal{A}(|\Psi {|}^2)+2E{\nabla }_{\mathbf{r}}^2|\Psi {|}^2-2R\frac{{\partial }^2|\Psi {|}^2}{\partial t^2}\right)\Psi$$

$$\tilde{D}=\frac{\partial \left(|\Psi {|}^{2}D\right)}{\partial \left(|\Psi {|}^2\right)}\mathcal{A}$$

$$\Psi (x)=\Delta \tanh \left(\frac{x}{\sqrt{2}\xi }\right)$$

$$F_{\text{EFT}}[\Psi ]={\int }_{-\infty }^{\infty }dx\left[X(|\Psi {|}^2)+C{\partial }_x\bar{\Psi }{\partial }_x\Psi -E{\left({\partial }_x|\Psi {|}^2\right)}^2\right]$$

$$X(|\Psi {|}^2)={\Omega }_s(|\Psi {|}^2)-{\Omega }_s(|{\Psi }_{\infty }{|}^2)$$

$$F=\sqrt{2}{\int }_{-\infty }^{\infty }du\left[\xi X(u)+\frac{\tilde{C}{\Delta }^2}{2\xi }{\text{sech}}^4\left(u\right)-\frac{\tilde{E}{\Delta }^4}{\xi }{\text{sech}}^4\left(u\right){\tanh }^2\left(u\right)\right]$$

$$\frac{dF}{d\xi }=\sqrt{2}\left[{\int }_{-\infty }^{\infty }X(u)du-\frac{2\tilde{C}{\Delta }^2}{3{\xi }^2}+\frac{4\tilde{E}{\Delta }^4}{15{\xi }^2}\right]$$

$$\xi =\sqrt{\frac{10\tilde{C}{\Delta }^2-4\tilde{E}{\Delta }^4}{15B}}$$

$$B={\int }_{-\infty }^{\infty }X(u)du$$

$$\varphi =\frac{\partial \Psi }{\partial t}$$

$$\bar{\varphi }=\overline{\frac{\partial \Psi }{\partial t}}=\frac{\partial \bar{\Psi }}{\partial t}$$

Full text

I Introduction

Arising from an interplay between nonlinear and dispersive effects, solitons are solitary matter waves which retain their shape while propagating at a constant velocity. They emerge in a wide variety of physical systems including optical fibers, classical fluids, plasmas and ultracold atomic gases. Due to their high experimental tunability, ultracold atom clouds in particular form an ideal environment for studying the physics of nonlinear matter waves in a well-controlled way. In these atomic clouds, solitons often manifest themselves as dark solitons, which are characterized by a localized density dip and a jump in the phase profile of the order parameter. Dark solitons have been theoretically and experimentally studied in Bose-Einstein condensates (BECs) Frantzeskakis (2010); Denschlag et al. (2000); Burger et al. (1999); Anderson et al. (2001); Becker et al. (2008), as well as superfluid Fermi gases Antezza et al. (2007); Scott et al. (2011); Liao and Brand (2011). In both systems, they are subject to an instability mechanism called the snake instability Ku et al. (2016); Muryshev et al. (1999); Brand and Reinhardt (2002); Cetoli et al. (2013), which makes the soliton decay into vortices if the radial width of the atom cloud is too large. The name “snake instability” comes from the observation that in the decay process the depletion plane of the soliton starts to oscillate until it breaks up into one or more vortex structures. The snaking mechanism and the subsequent decay of the dark soliton plane have been observed experimentally in both BECsDonadello et al. (2014) and superfluid Fermi gases Ku et al. (2016).

Various theoretical methods have been used to analyze the instability of dark solitons in fermionic systems Cetoli et al. (2013); Muñoz Mateo and Brand (2014); Lombardi et al. (2017). In Ref. Lombardi et al. (2017), the critical length scale of the instability was calculated in the BEC-BCS crossover by means of a recently developed low-energy effective field theory (EFT) Klimin et al. (2015). This formalism, which is capable of describing Fermi superfluids in a wide range of temperatures and interaction strengths, is based on the assumption that the order parameter of the system changes slowly in both space and time. It has already been successfully employed for the description of the properties and dynamics of dark solitons as a function of temperature and population imbalance Klimin et al. (2014); Lombardi et al. (2016); Alphen et al. (2018).

In the present paper, numerical simulations based on the EFT reveal that the dynamics of the dark soliton decay change significantly across the interaction domain. A perturbative analysis of the amplitude and phase fluctuations of the order parameter demonstrates that this crossover in the instability dynamics is related to a change in the nature of the unstable modes, which shift from amplitude-like to phase-like when one moves from the BEC- to the BCS-regime. All calculations are preformed for the case of a Fermi superfluid with a uniform background. While traditionally ultracold gases are studied in set-ups with harmonic trapping potentials, the recent realization of box-like optical traps Gaunt et al. (2013) provides an incentive to investigate uniform superfluids and the opportunity to experimentally test the predictions of the present work.

II Theoretical model

The system under consideration is an ultracold Fermi gas in which particles of opposite pseudo-spin interact via an $s$-wave contact potential. In the context of effective field theories, this system can be described in terms of a superfluid order parameter $\Psi(\mathbf{r},t)$, representing the bosonic field of Cooper pairs. Under the assumption that this field varies slowly in both space and time, a gradient expansion of the Euclidean-time action functional of the fermionic system can be performed, resulting in an effective action for the bosonic field:

[TABLE]

where we use the natural units of $\hbar=1$, $2m=1$, $E_{F}=1$, and where $\beta$ is the inverse temperature. A more detailed explanation of this model can be found in Ref. Klimin et al. (2015) and in Appendix A, together with the analytical expressions for the thermodynamic potential $\Omega_{s}$ and the coefficients $C$, $D$, $E$, $Q$ and $R$ in terms of the average chemical potential $\mu$ and the interaction parameter $(k_{F}a_{s})^{-1}$. The coefficients $C$, $E$, $Q$ and $R$ are only a function of the amplitude of the order parameter in bulk $|\Psi_{\infty}|$ (i.e. the superfluid gap). The coefficient $D$ and the thermodynamic potential $\Omega_{s}$, on the other hand, depend fully upon the local value of the amplitude of the order parameter Klimin et al. (2014). In this work, we assign to $|\Psi_{\infty}|$ and $\mu$ the mean-field values that are obtained by simultaneously solving the saddle-point gap and number equations Tempere and Devreese (2012). These background values could be further improved upon by, for example, including fluctuations around the saddle point, using the values of quantum Monte-Carlo simulations Carlson et al. (2003); Astrakharchik et al. (2004); Carlson and Reddy (2005) or even using the values derived from experimental measurements Ku et al. (2012).

From the Euclidian action functional (1), the real-time three-dimensional (3D) equation of motion for the pair field $\Psi(\mathbf{r},t)$ can be derived:

[TABLE]

where the coefficients $\tilde{D}$ and $\mathcal{A}$ are defined as

[TABLE]

This equation is a type of non-linear Schrödinger equation which is closely related to both the Gross-Pitaevskii equation for Bose-Einstein condensates Klimin et al. (2016) and the Ginzburg-Landau equation for Fermi superfluids (Sá de Melo et al., 1993). We find an analytical solution $\Psi_{s}(x)$ for a one-dimensional (1D) stationary dark soliton (also called a black soliton) by solving the time-independent equation

[TABLE]

with boundary conditions

[TABLE]

This solution for the order parameter possesses a phase jump $\pi$ and an amplitude equal to zero at the core of the soliton Klimin et al. (2014); Lombardi et al. (2016). When studying solitonic excitations across the BEC-BCS crossover, it is convenient to express the length scale in units of the healing length $\xi$, which is here defined as the width of the soliton. An analytic expression for $\xi$ can be derived through a variational ansatz for the stationary soliton solution and a minimization of the EFT free energy (see also Appendix B). Using the system parameters of the experiment in Ref. Ku et al. (2016), the variational result for $\xi$ yields a soliton width of about 550 nanometers at unitarity.

The main assumption of the EFT model is that the order parameter $\Psi(\mathbf{r},t)$ varies slowly in both space and time Klimin et al. (2015). In terms of spatial fluctuations, this assumption corresponds to the condition that the pair field should vary over a spatial region larger than the pair correlation length (also referred to as the Pippard correlation length in the context of superconducting systems). In terms of the frequency and energy of the collective excitations, the validity of the theory is mainly determined by the role of pair-breaking processes. Broken pairs are only present in the EFT in a local equilibrium state, there is no explicit pair-breaking dynamics. As a function of energy, pair-breaking processes become important at the bottom of the single-particle excitation spectrum ($2\Delta$ in the BCS-regime, $2\sqrt{\Delta^{2}+\mu^{2}}$ in the BEC-regime), which means the validity of the EFT can only be guaranteed for fluctuations of which the energy lies below this threshold value. A detailed study of the validity of the model reveals that the theory is reliable, except in some cases in the BCS-regime (Lombardi et al., 2016, Figure 5), where $\Delta$ becomes small and the pair correlation length becomes large. Earlier EFT calculations on the snake instability mechanism show a good agreement with results obtained from other theoretical formalisms Lombardi et al. (2017) across the whole BEC-BCS crossover. This is a consequence of the fact that the unstable mode is in general a long-wavelength mode, with energies sufficiently far below $2\Delta$. Accordingly, we also expect our current study on the nature of the instability mechanism to remain within the validity domain of the EFT.

III Results

To investigate the dynamics of the soliton instability across the BEC-BCS crossover, we perform numerical simulations of the decay of a stationary dark soliton in a uniformly trapped quasi-two-dimensional (2D) superfluid Fermi gas, using the EFT’s non-linear equation of motion (2). The initial state is constructed by extending the 1D stationary soliton solution $\Psi_{s}(x)$ into two dimensions (using periodic boundary conditions in the transverse direction) and adding a small amount of random noise to trigger the instability. The subsequent numerical time evolution is carried out by discretizing the space-time grid and applying a finite-difference fourth order Runge-Kutta (RK4) algorithm. For the present calculations, the longitudinal and transverse system size are chosen to be respectively $L_{x}=40\,\xi$ and $L_{y}=100\,\xi$, while the spatial and temporal resolution are taken to be respectively $5\%$ of $\xi$ and $2\%$ of $t_{F}=\omega_{F}^{-1}=(E_{F}/\hbar)^{-1}$. All calculations are carried out at a temperature $T=0.001\,T_{F}$, which for all practical purposes corresponds to zero temperature. A more detailed explanation of this procedure is given in Appendix C.

Figure 1 shows the evolution of the pair density $|\Psi(x,y)|^{2}$ of the superfluid during the decay of a stationary dark soliton. The upper row shows the evolution for $(k_{F}a_{s})^{-1}=1$ (BEC-regime), while the low row shows the evolution for $(k_{F}a_{s})^{-1}=-2$ (BCS-regime). To make the dynamics of the soliton core more apparent, the spatial profile of the pair density is displayed using only two colors: black if $|\Psi(x,y)|^{2}<0.1\,|\Psi_{\infty}|^{2}$ and white otherwise. In the BEC-regime, one can clearly observe the onset of snaking deformations of the depletion plane, which eventually break up the soliton into vortices. The manifestation of these snaking oscillations during the soliton decay is well-known from superfluid Bose gases. In the BCS-regime, on the other hand, the depletion plane does not exhibit snaking deformations during the decay process. Instead, the soliton core simply dissolves into vortices.

To further investigate the observed transition in the dynamics of the decay across the BEC-BCS crossover, we perform a perturbative study of the collective excitations around the stationary soliton solution. These excitations are represented by a small perturbation field $\delta\Psi(\mathbf{r},t)=\Psi(\mathbf{r},t)-\Psi_{s}(x)$. The equation of motion (2) can be linearized with respect to this perturbation field in order to describe small fluctuations of the order parameter (see Appendix D). To capture the physics of the instability more clearly, we also introduce the fields

[TABLE]

In the linearized theory, these fields correspond to fluctuations of respectively the amplitude and phase of the order parameter. Since we are interested in transverse unstable modes which propagate along the soliton plane in the $y$-direction and are localized around the soliton in the $x$-direction, we assume $P_{\pm}(\mathbf{r},t)=P_{\pm}(x)\,e^{i(ky-\omega t)}$. Then, the linearization of equation (2) results in two coupled equations for the amplitude and phase field:

[TABLE]

The expressions for the position-dependent coefficients $\alpha_{i}(x)$, $\beta_{i}(x)$ and $\gamma(x)$ are given in Appendix D. In general, the two equations are coupled by a coupling coefficient $\gamma(x)$. However, this coupling coefficient becomes very small with respect to the other coefficients on the BCS-side of the interaction domain, causing the amplitude and phase modes to decouple for $(k_{F}a_{s})^{-1}\ll-1$. On the BEC-side, on the other hand, $\gamma(x)$ becomes much larger, resulting in a strong coupling between the equations for $(k_{F}a_{s})^{-1}\gg 1$. Consequently, tuning the interaction parameter $(k_{F}a_{s})^{-1}$ from the BCS- to the BEC-regime increases the coupling between the amplitude and phase modes.

The spectrum of eigenmodes $\omega(k)$ is obtained by numerically solving the coupled equations (7) and (8). In particular, the soliton will be unstable if there is at least one mode for which $\text{Im}(\omega)>0$, since such a mode will grow exponentially in time. Additionally, we can determine for every mode whether it is more amplitude-like or phase-like in nature by calculating the norms of the eigenfields as

[TABLE]

and then using these quantities to define a mixing parameter

[TABLE]

with $0\leq\eta\leq 1$. For $\eta=1$, the excitation is a pure amplitude mode, while for $\eta=0$, it is purely a phase mode.

Figure 2 shows the real frequencies of several low-lying energy modes (right graph) and the imaginary part of the frequency of the unstable mode (left graph) in function of $(k_{F}a_{s})^{-1}$ for $k\approx 0$. The inset of the left graph shows an example of the profiles of the localized eigenfunctions $P_{+}(x)$ and $P_{-}(x)$ for the unstable mode at $(k_{F}a_{s})^{-1}=-0.9$. Each point on the main graphs is assigned a color based on the value of $\eta$ for the associated eigenmode: as $\eta$ goes from 0 to 1, the assigned color shifts from blue to red.

We observe that as we tune the interactions from the BCS-regime ( $(k_{F}a_{s})^{-1}<0$) to the BEC-regime ($(k_{F}a_{s})^{-1}>0$), the unstable mode changes from a phase-like mode to an amplitude-like mode. This shift in nature also appears to affect the behavior of the stable collective excitations across the interaction domain, inducing a steeply increasing amplitude mode which runs through the continuum of phase modes across the BEC-BCS crossover.

Figure 3 shows the imaginary part $\text{Im}(\omega)$ of the unstable modes for several values of $(k_{F}a_{s})^{-1}$. The insets of the top left and bottom right graphs show the localized eigenfunctions $P_{+}(x)$ and $P_{-}(x)$ for the modes which have been observed to manifest during the decay processes in Figure 1.

It is clear that the unstable mode is a long-wavelength mode which only exists up to a maximum wave number $k_{c}$. We observe that, in the deep BCS-regime, the unstable mode is a phase mode across the whole $k$-domain, while closer to unitarity and in the BEC-regime, its nature strongly depends on the value of $k$: for lower values of $k$ it is amplitude-like, while for higher $k$ it is phase-like. We can interpret these observations in terms of the width of the system, as for a given transverse system size $L_{y}$, the largest transverse mode which fits into the system corresponds to the wave number $k=2\pi/L_{y}$. If it is assumed that the mode with the largest growth rate will be the one to manifest itself, one would need a considerably small transverse width $L_{y}$ at unitarity or in the BEC-regime to observe the phase-like modes close to $k_{c}$. In the deep BCS-regime, on the other hand, the nature of the unstable mode will always be phase-like, independent of the system width. Considering the typical system sizes which are currently employed in dark soliton experiments in superfluid Fermi gases Yefsah et al. (2013); Ku et al. (2014, 2016), we expect the phase-like character of the unstable mode will only be observable in the BCS-regime.

Combining the results of the numerical simulations and the perturbative analysis now allows us to explain the change in the behavior of the soliton instability observed in Figure 1. In the BEC-regime, the mode with the largest growth rate has a finite wave number and is a mix between an amplitude-like and a phase-like mode (right bottom graph of Figure 3). The amplitude-like contribution causes the characteristic snaking oscillations of the soliton core 111The $k=0$ amplitude mode $\delta\Psi\propto\delta x\frac{d\Psi_{s}}{dx}$ corresponds to a uniform shift of the soliton in the $x$-direction over a distance $\delta x$. For finite $k$, the magnitude and direction of this translation vary along the soliton plane, leading to the creation of the characteristic snaking pattern., while the phase-like contribution creates local Josephson currents, leading to the manifestation of vortices. In the deep BCS-regime, on the other hand, unstable modes are purely phase-like (top left graph of Figure 3). As such, they cause vortex nucleation without the occurrence of snaking deformations. Therefore, we suggest that, on the BCS-side of the interaction domain, “Josephson instability” might be a more suitable name than snake instability to describe the unstable mode, as no observable snaking pattern is present. Moreover, since in the deep BCS-regime the imaginary part of the frequency reaches its maximum value at $k=0$, the number of generated vortices is expected to be smaller than in the BEC-regime, as is indeed the case in Figure 1.

IV Conclusions

The combined observations of the nonlinear dynamics and the perturbative analysis of the soliton decay indicate that the crossover in the behavior of the soliton instability is caused by the fact that the nature of the unstable mode changes from amplitude-like to phase-like when one moves from the BEC- to the BCS-side of the interaction domain. In the BCS-regime, where the unstable mode is purely phase-like in nature, the creation of local Josephson currents leads to the disintegration of the soliton core into vortices without changing the soliton’s position. In the BEC-regime, on the other hand, where the nature of the unstable mode is a mix between amplitude-like and phase-like, the nucleation of vortices is preceded by the onset of snaking deformations of the soliton core. We proposed to name the soliton instability in the BCS-regime the “Josephson instability”, in order to distinguish it from the snake instability in the BEC-regime.

We can compare our predictions with the results of the experiment carried out in Ref. Ku et al. (2016), in which the authors imaged the decay of a dark soliton in a 3D harmonically trapped superfluid Fermi gas at unitarity. By applying the system parameters of this experiment to the case of a uniformly trapped superfluid, we obtain that $t_{F}\approx 5.2\times 10^{-5}$ seconds, which yields a typical decay time on the order of (tens of) milliseconds for our numerical simulations across the BEC-BCS crossover. This time scale is in good agreement with the decay time observed in Ref. Ku et al. (2016), even though the experiment was performed in a harmonically trapped system. We furthermore observe that, on some of the experimental images in Ref. Ku et al. (2016), the decay process actually looks like a crossover between the snake instability mechanism and the aforementioned Josephson instability mechanism, indicating that, in practice, the predicted transition in the dynamics of the decay might already be observable around unitarity. This hypothesis will have to be investigated more thoroughly in future experiments. We suspect that large, quasi-2D systems, like the one employed in the experiment in Ref. Park et al. (2018), might be most suitable to test the current predictions, as they closely resemble the theoretical system described in this work.

Acknowledgements.

W. Van Alphen acknowledges financial support in the form of a Ph. D. fellowship of the Research Foundation - Flanders (FWO). This research was supported by the University Research Fund (BOF) of the University of Antwerp, by the Flemish Research Foundation (FWO-Vl) project nr G.0429.15.N., by JSPS KAKENHI Grants No. JP17K05549 and JP17H02938, and in part by the Osaka City University (OCU) Strategic Research Grant 2017 for young researchers.

Appendix A Overview of the EFT

In this section we provide a brief overview of the EFT model and the expressions for the EFT expansion coefficients. More detailed derivations and explanations can be found in Ref. Tempere and Devreese (2012); Klimin et al. (2015); Lombardi (2017).

The system of interest is an ultracold Fermi gas, in which particles of opposite pseudo-spin interact via an $s$-wave contact potential. The Euclidian-time action functional of this system can be written down in terms of the fermionic (Grassmann) fields $\psi_{\sigma}(\mathbf{x},\tau)$ and $\bar{\psi}_{\sigma}(\mathbf{x},\tau)$:

[TABLE]

where $g$ is the strength of the contact interaction and the label $\sigma$ denotes the spin degree of freedom. The quartic interaction term can be decoupled through the Hubbard-Stratonovich (HS) transformation, which introduces the bosonic pair field $\Psi(\mathbf{x},\tau)$ (the HS field is often also denoted as $\Delta$, but here we use $\Psi$ to emphasize its interpretation as a position- and time-dependent order parameter for the system) Tempere and Devreese (2012). The fermionic degrees of freedom can then be integrated out. If we assume that the pair field $\Psi(\mathbf{x},\tau)$ only varies slowly around its constant background value $\Psi_{\infty}$, we can perform a gradient expansion around $\Psi_{\infty}$ up to second order in the spatial and temporal derivatives of $\Psi(\mathbf{x},\tau)$ Klimin et al. (2015). This results in the following Euclidian-time effective action functional for the bosonic pair field:

[TABLE]

This effective action functional forms the starting point (1) for our study of the snake instability in the main work. The thermodynamic potential $\Omega_{s}$ is given by:

[TABLE]

while the gradient expansion coefficients $D$, $C$, $E$, $Q$ and $R$ are defined as

[TABLE]

The functions $f_{j}(\beta,\epsilon,\zeta)$ in the above expressions are defined by

[TABLE]

with the fermionic Matsubara frequencies $\omega_{n}=(2n+1)\pi/\beta$. In this treatment, the chemical potentials of the two pseudo-spin species $\mu_{\uparrow}$ and $\mu_{\downarrow}$ are combined into the average chemical potential $\mu=(\mu_{\uparrow}+\mu_{\downarrow})/2$ and the imbalance chemical potential $\zeta=(\mu_{\uparrow}-\mu_{\downarrow})/2$, the latter determining the difference between the number of particles in each spin-population. The quantity $\xi_{\mathbf{k}}=\frac{k^{2}}{2m}-\mu$ is the dispersion relation for a free fermion, $E_{\mathbf{k}}=(\xi_{\mathbf{k}}^{2}+|\Psi_{\mathbf{x},\tau}|^{2})^{1/2}$ is the local Bogoliubov excitation energy, and $a_{s}$ is the $s$-wave scattering length that determines the strength and the sign of the contact interaction. In absence of spatial and temporal variations, the thermodynamic potential $\Omega_{s}$ determines the value of the pair-breaking gap $|\Psi_{\infty}|$ for the uniform system through the saddle-point gap equation

[TABLE]

This equation is solved self-consistently together with the number equation to obtain the correct values of $|\Psi_{\infty}|$ and $\mu$ for a given set of system parameters.

In principle, all expansion coefficients (A.3)–(A.8) fully depend upon the order parameter $\Psi(\mathbf{x},\tau)$, but in practice, we assume that the coefficients associated with the second order derivatives of the pair field can be kept constant and equal to their bulk value, since retaining their full space-time dependence would lead us beyond the second-order approximation of the gradient expansion. This means that in expressions (A.5), (A.6), (A.7) and (A.8) for the coefficients $C$, $E$, $Q$ and $R$, we set $|\Psi(\mathbf{x},\tau)|^{2}\rightarrow|\Psi_{\infty}|^{2}$ and $E_{\mathbf{k}}\rightarrow E_{\mathbf{k},\infty}=(\xi_{\mathbf{k}}^{2}+|\Psi_{\infty}|^{2})^{1/2}$. For the thermodynamic potential $\Omega_{s}$ and the coefficient $D$, on the other hand, the full space-time dependence of the order parameter is preserved.

From the Euclidian-time action functional (A.2), the EFT equation of motion for the pair field $\Psi(\mathbf{r},t)$ is found to be

[TABLE]

where the coefficients $\tilde{D}$ and $\mathcal{A}$ are defined as

[TABLE]

The first term on the right-hand side of the equation can be identified as a kinetic energy term, while the non-linear term represents a system-inherent potential for the field. The ratio $\tilde{D}/C$ can be interpreted as a renormalization factor for the mass of the fermion pairs Klimin et al. (2016) and the coefficient $\mathcal{A}$ determines the uniform background value of the system, since $\mathcal{A}(\Psi)\,\Psi=0$ is nothing but the aforementioned gap equation (A.10). It has been verified that in the deep BEC-limit $\left(1/k_{F}a_{S}\gg 1\right)$, the equation correctly tends to the Gross-Pitaevskii equation for bosons with a mass $M=2m$ and an s-wave boson-boson scattering length $a_{B}=2\,a_{s}$ (Lombardi, 2017).

Appendix B Variational derivation of the healing length

We can derive an analytic expression for the healing length $\xi$ associated with a stationary (black) soliton in a superfluid by considering a variational ansatz for the wavefunction and minimizing the free energy of the system. Since the wavefunction of a stationary soliton is expected to be real and anti-symmetric, we assume a hyperbolic tangent model for the pair field:

[TABLE]

The EFT free energy for a one-dimensional Fermi superfluid is given by

[TABLE]

with

[TABLE]

The subtraction of the term $\Omega_{s}(|\Psi_{\infty}|^{2})$ indicates that the energy is measured with respect to the energy of the uniform system. By substituting the ansatz (A.13) for the pair field into the free energy and making a change of integration variable $u=x/(\sqrt{2}\xi)$, we obtain

[TABLE]

The integrals in the second and third term can be computed analytically. The derivative of $F$ with respect to the variational parameter $\xi$ then becomes:

[TABLE]

By setting the above equation equal to zero, we find the following variational expression for the healing length:

[TABLE]

with

[TABLE]

Figure 4 shows the behavior of this quantity in function of the interaction parameter $(k_{F}a_{s})^{-1}$. A similar expression was derived in the context of the EFT for the width of a vortex core in Ref. Verhelst et al. (2017). A more extensive study on the healing length of a fermionic superfluid across the BEC-BCS crossover can be found in Ref. Palestini and Strinati (2014).

Appendix C Discretization and evolution of the equation of motion

In this section we elaborate on how the EFT equation of motion (2) is discretized and solved numerically using the explicit RK4 algorithm. We introduce a field $\phi(\mathbf{r},t)$ such that

[TABLE]

and

[TABLE]

Substituting this into the equation of motion and making use of the fact that

[TABLE]

we have

[TABLE]

In order to get an equation of the form $\partial_{t}\phi=...\,$, we take the complex conjugate of (A.23), find an expression for $\partial_{t}\bar{\phi}$ in function of $\partial_{t}\phi$ and substitute this back into (A.23), yielding

[TABLE]

Equations (A.20) and (A.24) form a system of two coupled partial differential equations of the form:

[TABLE]

where $f(\phi)=\phi$ and $g(\Psi,\phi)$ is given by (A.24). In the case of a 2D system, we use finite mesh widths $\Delta x$ and $\Delta y$ and a finite time step $\Delta t$ to discretize space-time into a grid of $L\times M\times N$ points by writing $x_{l}=l\Delta x$ with $l=1,...,L$, $y_{m}=m\Delta y$ with $m=1,...,M$ and $t_{n}=n\Delta t$ with $n=1,...,N$. This allows us to approximate the spatial derivatives by central finite difference formulas:

[TABLE]

where we use the notation $\Psi_{l,m,n}=\Psi(x_{l},y_{m},t_{n})$. Since we expect the superfluid to assume its uniform bulk value sufficiently far from the soliton, we require the derivatives of he fields to be zero at the $x$-boundaries of the grid. In the $y$-direction, we apply periodic boundary conditions. If we now know the values $\Psi_{l,m,n}$ and $\phi_{l,m,n}$ at a certain time step $t_{n}$ for all positions $x_{l}$ and $y_{m}$, the explicit RK4 method allows us to calculate for every position the values $\Psi_{l,m,n+1}$ and $\phi_{l,m,n+1}$ of the next time step by using the following algorithm (Süli and Mayers, 2003):

[TABLE]

This scheme can be repeated until the solution has been evolved up to the desired point in time.

Appendix D Linearization of the equation of motion

To describe small fluctuations of the pair field, we add a perturbation field $\delta\Psi(\mathbf{r},t)$ to the stable soliton solution $\Psi_{s}(x)$:

[TABLE]

We plug this perturbed solution into the equation of motion (2) and expand the coefficients $\tilde{D}$ and $\mathcal{A}$ (which depend on the local value of the order parameter) up to first order in the perturbation field:

[TABLE]

In order to study explicitly the character of amplitude and phase modes, we introduce the fields

[TABLE]

If we then collect all terms of first order in the fluctuation fields, we find two linearized equations

[TABLE]

where the coefficients $\alpha_{i}(x)$, $\beta_{i}(x)$ and $\gamma(x)$ are given by

[TABLE]

Here, we have used the notations

[TABLE]