# Swallow-tail dispersions of moving solitons in a two-dimensional fermionic
superfluid

Jan Major<EMAIL_ADDRESS>Dodd-Walls Centre for Photonic and Quantum
Technologies, Centre for Theoretical Chemistry and Physics, New Zealand
Institute for Advanced Study, Massey University, Private Bag 102904, North
Shore, Auckland 0745, New Zealand  Institute of Theoretical Physics, Wrocław
University of Science and Technology, 50-370 Wrocław, Poland  Joachim Brand
Dodd-Walls Centre for Photonic and Quantum Technologies, Centre for
Theoretical Chemistry and Physics, New Zealand Institute for Advanced Study,
Massey University, Private Bag 102904, North Shore, Auckland 0745, New Zealand

(August 28, 2024)

###### Abstract

Soliton-like localised wave solutions in a two-dimensional Fermi superfluid
are studied by solving the Bogoliubov-de Gennes equations in the BCS regime of
weak pairing interactions. The dispersion relations of these solitons are
found to exhibit a peculiar swallow-tail shape, with cusps and multiple
branches. The effective mass of the solitons is found to diverge and change
sign at the cusp. This behavior is in contrast to the smooth dispersion
relations and negative effective masses of solitons in the three-dimensional
Fermi superfluid. The swallow-tail dispersion relations are shown to be a
consequence of counterflow of the superfluid and sign-changing contributions
to the superfluid current from different transverse momenta in the Bogoliubov-
de Gennes formalism. The results are relevant for the understanding of
solitonic excitations in two-dimensional Fermi superfluids, such as ultracold
atomic gases and high-temperature superconductors.

##  I Introduction

Phenomena of superfluidity and superconductivity are known already for more
than a hundred years [1, 2, 3], yet they still inspire new research and
provide new riddles to solve. The same can be said about solitons with an even
longer research history [4] and connections to transport and critical
phenomena in a wide variety of systems: from DNA and polymers to cosmological
models [5, 6]. Solitonic excitations in superfluids appear in form of dark
solitons, which are characterised by an abrupt change of the complex phase of
the superfluid order parameter [7]. Dark solitons were confirmed and
characterised experimentally in Bose-Einstein condensates [8, 9, 10], but in
fermionic superfluids, including superconductors, their existence is still a
matter of debate. Only observations of transient phenomena attributed to
solitons have been reported so far [11]. Theoretical studies of solitons in
fermionic superfluids have been conducted in one dimension [12, 13, 14, 15]
and three dimensions [16, 17, 18, 19, 20, 21, 22], but to the best of our
knowledge, no studies of moving solitons in two-dimensional fermionic
superfluids have been reported.

The thorough understanding of soliton properties and their connection to
transport in fermionic superfluids (including superconductors) may be vital
for progress in such fields as high-$T_{c}$ superconductivity [23] or in the
analysis of novel topologically nontrivial superfluid systems [24]. This paper
focuses on the case of a balanced s-wave superfluid in two dimensions in the
crossover regime close to the Bardeen-Cooper-Schrieffer (BCS) limit. In this
regime we observe some previously unseen behavior that is directly connected
to the dynamic properties of solitons. For low binding energies, the
dispersion relation undergoes a qualitative change, from a smooth inverted
lobe into a swallow-tail shape with non-smooth cusps. From this change results
a range of phenomena such as the coexistence of several solitonic solutions
for the same velocity. We also find sign changes, zero crossings and
divergence of the soliton’s effective mass.

As the movement of solitons in many classes of systems is hindered by
impurities, the ideal medium for studying soliton dynamics are “clean” systems
of ultracold atoms [25]. Apart from the lack of defects (unless intentionally
included), those systems allow for unprecedented control over the parameters.
This makes it possible, for example, to choose the dimensionality, create
external static and dynamic potentials, and generate artificial gauge fields
such as spin-orbit coupling [26]. The most useful in the context of this paper
is the possibility of changing the interparticle interaction strength via
magnetic-field tunable Feshbach resonances [27]. This means that a range of
models can be simulated: from the weakly-coupled BCS phase through the
crossover to the Bose-Einstein condensate (BEC) of tightly bound fermionic
pairs.

There are several theoretical descriptions of solitons in fermionic
superfluids. In one dimensional superfluids, solitons have been studied both
in full many body framework using an exactly-solvable model [12] as well as
using an analytically solvable mean field approximation [13], showing some
significant differences. In three dimensions mean field theory generally
requires numerical solutions [16, 28] although scaling relations provide some
analytical insights for the unitary gas [17]. Our study aims at filling the
gap for two-dimensional case and uses numerical solutions of the Bogoliubov-de
Gennes (BdG) equations. We are focusing on the deep BCS regime near the BCS
limit. The reason for that is twofold: It is known that the BdG approach fails
in the BEC regime in two dimensions but is expected to be safer to use in
weakly interacting regime [29]. Moreover, we find very interesting phenomena
for systems in the deep BCS regime.

In Sec. II we describe the fermionic superfluid model studied in this paper as
well as the mean-field BdG approach used to simulate its excitations. In Sec.
III we show results of the self-consistent calculations hosting solitonic
excitations and examine their spatial profiles. The dispersion relations for
these solitons are presented in Sec. IV along with calculations of their
inertial and physical masses. In section V we focus on the phenomenon of
counterflow – the fact that the sign of contributions to superfluid current
may vary for different transverse momenta. Finally we conclude in VI.

Appendix A contains a discussion of how superfluid flow affects the chemical
potential and App. B discusses the change of the free energy of a soliton due
to a change of reference frames.

##  II Model and method

We are considering a balanced Fermi gas with contact $s$-wave interactions at
zero temperature. The grand canonical Hamiltonian of such a system reads:

$\displaystyle\mathcal{H}=\sum_{\sigma}\int\mathrm{d}\mathbf{r}\left\\{\psi^{\dagger}_{\sigma}(\mathbf{r})\hat{h}\psi_{\sigma}(\mathbf{r})+g\psi^{\dagger}_{\sigma}(\mathbf{r})\psi_{\bar{\sigma}}^{\dagger}(\mathbf{r})\psi_{\bar{\sigma}}(\mathbf{r})\psi_{\sigma}(\mathbf{r})\right\\},$
(1)

where $\hat{h}=-\frac{\hbar^{2}}{2m}\nabla^{2}-\mu$ is the one particle
Hamiltonian, with $\mu$ the chemical potential and $m$ the fermion mass. $g$
is a renormalized interaction constant, $\sigma$ denotes the spin degrees of
freedom ($\pm\nicefrac{{1}}{{2}}$) and $\bar{\sigma}$ is opposite to $\sigma$.
A standard way of deriving a mean-field theory for an inhomogeneous superfluid
system is the BdG approach [30]. As it is a mean-field method, quantum
correlations are neglected and it is expected to fail in low dimensions. In
two-dimensional systems it is considered to be a good approximation in the BCS
limit, while it is known to fail when approaching the BEC limit [29].

In this framework the system is described using the self-consistent BdG
equation:

$\displaystyle i\hbar\partial_{t}\left[\begin{array}[]{c}u(\mathbf{r},t)\\\
v(\mathbf{r},t)\end{array}\right]=\left[\begin{array}[]{cc}\hat{h}&\Delta(\mathbf{r},t)\\\
\Delta^{*}(\mathbf{r},t)&-\hat{h}\end{array}\right]\left[\begin{array}[]{c}u(\mathbf{r},t)\\\
v(\mathbf{r},t)\end{array}\right].$ (8)

The $u(\mathbf{r},t)$ and $v(\mathbf{r},t)$ are quasi-particle amplitudes
which act as weights for the approximated many body wave function. The
superconducting order parameter $\Delta(\mathbf{r},t)$ has to be determined
self-consistently from a complete and orthogonal set of solutions for the
quasi-particle amplitudes [25].

We will further restrict ourselves to two-dimensional systems and consider
solutions which are: 1) inhomogeneous only in the $x$ direction, 2) are moving
with velocity $v_{s}$ in direction $x$. Because translational invariance in
the transverse $y$ direction is maintained, the BdG matrix splits into
disconnected blocks according to different transverse momentum channels. We
rewrite it in the frame co-moving with the soliton [17]. Eventually we can
obtain a time independent version of Eq. (8) for each transverse momentum
$p_{\mathrm{t}}$:

$\displaystyle
H^{(p_{\mathrm{t}})}_{\mathrm{BdG}}\left[\begin{array}[]{c}u_{\nu}^{(p_{\mathrm{t}})}(x)\\\
v_{\nu}^{(p_{\mathrm{t}})}(x)\end{array}\right]=\varepsilon_{\nu}^{(p_{\mathrm{t}})}\left[\begin{array}[]{c}u_{\nu}^{(p_{\mathrm{t}})}(x)\\\
v_{\nu}^{(p_{\mathrm{t}})}(x)\end{array}\right],$ (13)

where $\nu$ is the index of the eigenstate and
$\varepsilon_{\nu}^{(p_{\mathrm{t}})}$ is the corresponding eigenenergy. The
BdG matrix for each $p_{\mathrm{t}}$ is given by

$\displaystyle
H^{(p_{\mathrm{t}})}_{\mathrm{BdG}}=\left[\begin{array}[]{cc}\hat{h}_{x}+\frac{p_{\mathrm{t}}^{2}}{2m}&\Delta(x)e^{2imv_{\mathrm{s}}x/\hbar}\\\
\Delta^{*}(x)e^{-2imv_{\mathrm{s}}x/\hbar}&-\hat{h}_{x}+\frac{p_{\mathrm{t}}^{2}}{2m}\end{array}\right],$
(16)

where $\hat{h}_{x}=-\frac{\hbar^{2}}{2m}\partial_{x}^{2}-\mu$ is a
differential operator acting on $x$ only. The quasiparticle amplitudes
$u_{\nu}^{(p_{\mathrm{t}})}(x)$ and $v_{\nu}^{(p_{\mathrm{t}})}(x)$ are
eigenstates of the respective $H^{(p_{\mathrm{t}})}_{\mathrm{BdG}}$ matrix.
$\Delta(x)$ has to be calculated self-consistently as

$\displaystyle\Delta(x)=-g\sum_{\nu}^{\prime}\sum_{p_{\mathrm{t}}}u_{\nu}^{(p_{\mathrm{t}})}(x){v_{\nu}^{(p_{\mathrm{t}})}}^{*}(x),$
(17)

combining contributions from all transverse momentum channels. The particle
number density is given by

$\displaystyle
n(x)=\sum_{\nu}^{\prime}\sum_{p_{\mathrm{t}}}|v_{\nu}^{(p_{\mathrm{t}})}(x)|^{2},$
(18)

where the primed sum is done over states with positive eigenenergies
($\sum^{\prime}_{\nu}\equiv\sum_{\nu:\varepsilon_{\nu}^{(p_{\mathrm{t}})}>0}$).

Solutions resembling dark solitons on an infinite domain are expected to have
a phase jump in $\Delta(x)$,
$\phi_{\mathrm{s}}=\arg[\Delta(\infty)/\Delta(-\infty)]$, which is a signature
of the soliton and is expected to depend on $v_{\mathrm{s}}$ [7, 17, 18]. In
addition, the BdG matrix of Eq. (16) has a phase gradient in the off-diagonal
terms, originating from a frame transformation into the soliton’s rest frame.
In order to solve the BdG equations with periodic (or anti-periodic) boundary
conditions, it is therefore convenient to use a unitary transformation that,
in combination with a convenient boundary condition, allows us to keep only
the residual backflow in the off-diagonal terms. This transformation is given
by:

$\displaystyle U_{v_{\mathrm{fr}}}=\left[\begin{array}[]{c
c}\exp\left(-i\frac{m}{\hbar}v_{\mathrm{fr}}x\right)&0\\\
0&\exp\left(i\frac{m}{\hbar}v_{\mathrm{fr}}x\right)\end{array}\right],$ (21)

and using it leads to:

$\displaystyle\tilde{H}^{(p_{\mathrm{t}})}_{\mathrm{BdG}}=\left[\begin{array}[]{c
c}\hat{h}_{v_{\mathrm{fr}}}+\frac{p_{\mathrm{t}}^{2}}{2m}&\Delta(x)e^{-i\frac{2m}{\hbar}v_{\mathrm{bf}}x}\\\
\Delta^{*}(x)e^{+i\frac{2m}{\hbar}v_{\mathrm{bf}}x}&-\hat{h}_{-v_{\mathrm{fr}}}+\frac{p_{\mathrm{t}}^{2}}{2m}\end{array}\right].$
(24)

Here $v_{\mathrm{fr}}$ is associated with the unitary transformation,
$v_{\mathrm{bf}}=v_{\mathrm{s}}-v_{\mathrm{fr}}$ has the interpretation of a
backflow velocity, and
$\hat{h}_{v_{\mathrm{fr}}}=\frac{\hbar^{2}}{2m}(-i\partial_{x}+mv_{\mathrm{fr}}/\hbar)^{2}-\mu$.
In this way the soliton velocity $v_{\mathrm{s}}$ is split into the velocity
of the condensate backflow $v_{\mathrm{bf}}$ and $v_{\mathrm{fr}}$, which is a
free parameter to be set in the calculations. Enforcing periodic (anti-
periodic) boundary conditions in the $x$ direction, will ensure that
$\arg[\Delta(L_{x}/2)/\Delta(-L_{x}/2)]-2mv_{\mathrm{bf}}L_{x}/\hbar=0$
($\pi$) 111Due to $\Delta(x)$ being a sum of products of the eigenstates of
the BdG matrix, phase twists applied to either diagonal sector of BdG matrix
($\phi_{u}$ and $\phi_{v}$, respectively) add in contributing to $\Delta(x)$,
while the individual phases do not have physical significance. In our
calculations we usually set $\phi_{u}$ to $0$ and $\phi_{v}$ to $\pi$.. In
sufficiently large systems with $\pi$-twisted boundary conditions (as used
through most of our calculations) $v_{\mathrm{fr}}$ and $v_{\mathrm{s}}$
differ only sightly, as $v_{{\mathrm{bf}}}\ll v_{\mathrm{fr}}\approx
v_{\mathrm{s}}$.

Having found a self-consistent $\Delta(x)$ and thus knowing an approximate
eigenfunction of the many-body Hamiltonian $\mathcal{H}$ one can calculate
various observables of the system. The free energy of the soliton-bearing
system in the frame co-moving with the soliton is given by

$\displaystyle F^{\mathrm{SF}}_{\mathrm{s}}$
$\displaystyle=\langle\mathcal{H}\rangle_{\mathrm{s}}=\int\mathrm{d}x\sum_{\nu,p_{\mathrm{t}}}^{\prime}\varepsilon_{\nu}^{(p_{\mathrm{t}})}|u_{\nu}^{(p_{\mathrm{t}})}(x)|^{2},$
(25)

where the subscript s denotes system with a soliton, and the superscript SF
means the soliton frame of reference. In order to characterize the properties
of the soliton itself it is more interesting to consider the change in free
energy caused by adding a soliton rather than the (extensive) quantity
$F^{\mathrm{SF}}_{\mathrm{s}}$, which is the free energy of the whole system.
Moreover, to compare different cases, the free energy has to be brought back
to the laboratory frame of reference [32]. Thus we have to calculate the
difference between the free energy of system with ($F_{\mathrm{s}}$) and
without ($F_{0}$) a soliton in the lab frame, and make an appropriate
transformation to obtain a relation to the observables in the soliton frame
(see App. B). The soliton (free) energy is then defined as

$\displaystyle\Delta
F=F_{\mathrm{s}}-F_{0}=F_{\mathrm{s}}^{\mathrm{SF}}-F_{0}^{\mathrm{SF}}+N_{\mathrm{d}}mv_{\mathrm{s}}^{2}.$
(26)

Here, $N_{\mathrm{d}}<0$ is particle number depleted by the soliton

$\displaystyle
N_{\mathrm{d}}=\int[n_{\mathrm{s}}(x)-n_{0}]\mathrm{d}\mathbf{r}=N_{\mathrm{s}}-N_{0},$
(27)

where $N_{\mathrm{s}}$ is the number of particles in the soliton-bearing
system with density $n_{\mathrm{s}}(x)$, and $N_{0}$ is the number of
particles in the system without the soliton with the constant background
density $n_{0}$.

Next we can calculate the physical momentum of the system:

$\displaystyle P^{\mathrm{SF}}_{\mathrm{ph}}$
$\displaystyle=\sum_{\sigma}\int\mathrm{d}\mathbf{r}\left\langle\psi\left|\psi_{\sigma}^{\dagger}(\mathbf{r})(-i\hbar\partial_{x}+mv_{\mathrm{fr}})\psi_{\sigma}(\mathbf{r})\right|\psi\right\rangle$
$\displaystyle=4\hbar\int\mathrm{d}x\sum_{\nu,p_{\mathrm{t}}}^{\prime}{v_{\nu}^{(p_{\mathrm{t}})*}}(x)\partial_{x}v_{\nu}^{(p_{\mathrm{t}})}(x).$
(28)

In this frame of reference this is the sum of the physical momentum of the
soliton and a momentum contribution from the $v_{\mathrm{fr}}$. Because we are
using finite system sizes with an arbitrarily set phase twist
$\phi_{\mathrm{t}}$ at the boundary, we have to add a boundary contribution to
obtain the canonical momentum of the soliton

$\displaystyle
P^{\mathrm{SF}}_{\mathrm{c}}=P^{\mathrm{SF}}_{\mathrm{ph}}+P^{\mathrm{SF}}_{\phi_{\mathrm{t}}},$
(29)

where
$P^{\mathrm{SF}}_{\phi_{\mathrm{t}}}=-\frac{\hbar}{2}\frac{\phi_{\mathrm{t}}}{L_{x}}N_{\mathrm{s}}$.
After transforming the momentum to the lab frame we obtain a more natural
splitting:

$\displaystyle P_{c}=P_{\mathrm{ph}}+P_{\phi_{\mathrm{s}}},$ (30)

as $P_{\mathrm{ph}}$ is just physical momentum of the soliton:
$P_{\mathrm{ph}}=N_{\mathrm{d}}mv_{\mathrm{s}}$, while
$P_{\phi_{\mathrm{s}}}=-\frac{\hbar}{2}\frac{\phi_{\mathrm{s}}}{L_{x}}N_{\mathrm{s}}$
is proportional to the phase step of the soliton.

The dispersion relation of $\Delta F$ vs. $P_{\mathrm{c}}$ of the soliton
solutions is parameterized by the soliton velocity $v_{\mathrm{s}}$. Besides
the interaction strength, these quantities also depend on the chemical
potential $\mu$, which parameterizes the background particle number density.
In the presence of a slowly-varying external potential (on the length scale of
the soliton), the dispersion relations can be used to predict the motion of
the soliton by treating the soliton like a quasiparticle [33]. The free energy
of the soliton at position $z(t)$ then depends on its velocity
$v_{\mathrm{s}}=\dot{z}$ and the chemical potential at local equilibrium
$\mu(z)$, which is modified by the external potential according to the local-
density approximation. This leads to an equation of motion that resembles
Newton’s equation [33, 34]:

$\displaystyle\frac{1}{v_{\mathrm{s}}}\left.\frac{\partial\Delta F}{\partial
v_{\mathrm{s}}}\right|_{\mu}\ddot{z_{\mathrm{s}}}=-m\left.\frac{\partial\Delta
F}{\partial\mu}\right|_{v_{\mathrm{s}}}\frac{\partial\mu}{\partial z}.$ (31)

The coefficients on the left and right sides of the equation can be
interpreted as the inertial and the physical mass, respectively:

$\displaystyle M_{\mathrm{P}}$ $\displaystyle=m\left.\frac{\partial\Delta
F}{\partial\mu}\right|_{v_{\mathrm{s}}},$ (32) $\displaystyle M_{\mathrm{I}}$
$\displaystyle=\frac{1}{v_{\mathrm{s}}}\left.\frac{\partial\Delta F}{\partial
v_{\mathrm{s}}}\right|_{\mu}.$ (33)

Near the minima or maxima of the local chemical potential $\mu(z)$, the masses
can be approximated as constants and Eq. (31) describes harmonic oscillations.
The frequency of these small-amplitude oscillations is then related to the
mass ratio by $\omega_{z}^{2}/\Omega^{2}=M_{\mathrm{I}}/M_{\mathrm{P}}$, where
$\Omega$ is the harmonic oscillator frequency related to the extremum of the
external potential [34].

##  III Self consistent calculations

The calculations were conducted for $N$ fermions confined in a two-dimensional
box with dimensions $L_{x}\times L_{y}$ and periodic (torus) boundary
conditions for a wide range of binding energies $E_{\mathrm{b}}$. To obtain a
soliton solution we often apply an additional phase change of $\pi$ in the $x$
direction, realizing antiperiodic boundary condition, in order to offset the
main contribution from the soliton phase step, as described in Sec. II.

The renormalized interaction strength $g$ for a two-dimensional Fermi gas
depends on $E_{\mathrm{b}}$ as

$\displaystyle\frac{1}{g}=\sum_{p_{x}=-p_{x}^{\mathrm{cutoff}}}^{p_{x}^{\mathrm{cutoff}}}\sum_{p_{y}=-p_{y}^{\mathrm{cutoff}}}^{p_{y}^{\mathrm{cutoff}}}\frac{1}{\frac{1}{m}(p_{x}^{2}+p_{y}^{2})+E_{\mathrm{b}}}.$
(34)

where the values of $p_{x}$ and $p_{y}$ are integer multiples of
${2\pi\hbar}/{L_{x}}$ and ${2\pi\hbar}/{L_{y}}$, and $p_{x}^{\mathrm{cutoff}}$
and $p_{y}^{\mathrm{cutoff}}$ are the respective momentum cutoffs. By changing
$E_{\mathrm{b}}$ from $0$ to $+\infty$ we can obtain systems in both BCS and
BEC regimes. Additional calculations in one-dimensional systems were made for
comparison.

While presenting results we are using Fermi units based on the Fermi energy
$\varepsilon_{\mathrm{F}}=\pi\hbar^{2}N/mL_{x}L_{y}=\hbar^{2}k_{\mathrm{F}}^{2}/2m$,
from which one can further derive the length scale
$l_{\mathrm{F}}=2\pi/k_{\mathrm{F}}$ and the velocity
$v_{\mathrm{F}}=k_{\mathrm{F}}/m$. Other important quantities used throughout
the paper are BCS correlation length $\xi_{\mathrm{BCS}}=\hbar
v_{\mathrm{F}}/\Delta_{0}$ and the pair-breaking velocity
$v_{\mathrm{pb}}\equiv\Delta_{0}/\hbar
k_{\mathrm{F}}=v_{\mathrm{F}}\sqrt{E_{\mathrm{b}}/2\varepsilon_{\mathrm{F}}}$
for the homogeneous two-dimensional Fermi superfluid with
$\Delta_{0}=\sqrt{2E_{\mathrm{b}}\varepsilon_{\mathrm{F}}}$ and
$\mu_{0}=\varepsilon_{\mathrm{F}}-E_{\mathrm{b}}/2$ [35].

In the two dimensional system the transverse direction $y$ is assumed to be
homogeneous. For the simulations we have used a square lattice discretization.
Lattice spacings were set to $dL_{x}=dL_{y}=0.5l_{\mathrm{F}}$ for
$E_{\mathrm{b}}<\varepsilon_{\mathrm{F}}$. For higher $E_{\mathrm{b}}$ and for
few other cases when additional precision were needed the lattice spacings
were set to $0.25l_{\mathrm{F}}$. The lattice sizes were set to at least
$L_{x}=15\xi_{BCS}$ for the inhomogeneous direction and $L_{y}=10\xi_{BCS}$
for homogeneous one in the BCS regime. In the BEC regime we used at least
$L_{x}=40l_{\mathrm{F}}$ and $L_{y}=30l_{\mathrm{F}}$, respectively. Gradients
were approximated using finite differences with a highly-accurate 9-point
stencil. Transverse momenta were summed in the range
$p_{\mathrm{t}}\in(-\pi\hbar/dL_{y},\pi\hbar/dL_{y}]$, with discretization
$dp_{\mathrm{t}}=2\pi\hbar/L_{y}$.

In order to obtain a self-consistent solution one has to start from some
assumed form of the order parameter $\Delta_{0}(x)$ in order to construct an
initial BdG matrix (24). The respective eigenstates are used to calculate the
new order parameter using (17) – which we denote $\Delta_{0}^{\prime}(x)$. In
general $\Delta_{0}(x)\neq\Delta_{0}^{\prime}(x)$ and thus a new iterate
$\Delta_{1}(x)$ has to be constructed leading to an iterative procedure.
Convergence is achieved when $||\Delta_{n}-\Delta^{\prime}_{n}||<\epsilon$
with some prescribed accuracy $\epsilon$. In our case the simplest self-
consistent iteration scheme with $\Delta_{n}(x)=\Delta^{\prime}_{n-1}(x)$ does
not work, as it is unstable when applied to moving solitonic solutions (i.e.
when $v_{\mathrm{s}}\neq 0$). Instead, we use Broyden’s method [36], which is
a quasi-Newton method where the Jacobian is not explicitly calculated at every
step. As a starting point for the search we have used a step function or
hyperbolic tangent (with the length scale $\xi_{\mathrm{BCS}}$). The initial
Jacobians for Broyden’s method were calculated explicitly using finite
differences with a 2-point stencil. Broyden iterations were continued until
the difference between the 2-norms and $\infty$-norms of $\Delta(x)$ and
$\Delta^{\prime}(x)$ were smaller than the threshold values, typically set-up
around $\epsilon=10^{-6}\varepsilon_{\mathrm{F}}$, which provided well
converged results in most cases. All calculations were done using the Julia
programming language [37] and its standard linear algebra solvers based on
LAPACK [38] were used for matrix diagonalization.

Figure 1: Absolute value of the order parameter $|\Delta(x)|$ in a system with
stationary dark soliton solution in the BCS regime of a two-dimensional Fermi
superfluid. The self consistent result is shown for
$E_{\mathrm{b}}=0.02\varepsilon_{\mathrm{F}}$, $v_{s}\approx 0$ and
$\phi_{\mathrm{s}}\approx\pi$ (blue solid line). The red dashed line is
$|\tanh(x/\xi_{\mathrm{BCS}})|$.

A typical example of a self-consistent solution in the BCS regime with a
soliton can be seen in Fig. 1. For the absolute value $|\Delta(x)|$ we can
observe a wide envelope given approximately by $|\tanh(x/\xi_{\mathrm{BCS}})|$
with additional oscillations at the Fermi length scale $l_{\mathrm{F}}$ and a
central through, which is almost unaffected by the binding energy. We
attribute the oscillations to Friedel oscillations and note that they appear
also in quasi-particle eigenstates (not shown). The general shape of the
soliton order parameter is reminiscent of the dark soliton in the 3D Fermi
superfluid at unitarity [16].

Figure 2: The phase jump for families of dark/grey soliton solutions. Panel
(a): Phase step of the soliton in the two-dimensional superfluid as a function
of the soliton velocity for various binding energies:
$E_{\mathrm{b}}=0.005\varepsilon_{\mathrm{F}}$ – solid red,
$E_{\mathrm{b}}=0.009\varepsilon_{\mathrm{F}}$ – blue dashed,
$E_{\mathrm{b}}=0.03\varepsilon_{\mathrm{F}}$ – yellow dotted,
$E_{\mathrm{b}}=0.05\varepsilon_{\mathrm{F}}$ – green dash-dotted. Panel (b):
The same relation for a one-dimensional BdG model (setting $p_{\mathrm{t}}=0$)
with $g=-2\varepsilon_{\mathrm{F}}$ – blue dashed and
$g=-3\varepsilon_{\mathrm{F}}$ – yellow dotted lines. The black solid line is
a theoretical prediction for a 1D model [13].

In Fig. 2 we present phase-steps of the solitons as a function of the soliton
velocity $v_{\mathrm{s}}$ (scaled with pair-breaking velocity
$v_{\mathrm{pb}}$) for various interaction strengths. The soliton phase step
is calculated as a difference between phases on the edges of the system
transformed to a frame in which there is no backflow:

$\displaystyle\phi_{\mathrm{s}}$
$\displaystyle={\arg(\Delta(x\rightarrow-\infty)/\Delta(x\rightarrow+\infty))}|_{v_{\mathrm{bf}}=0}$
$\displaystyle\approx\left.{\arg\left(\Delta\left(-{L_{x}}/{2}\right)/\Delta\left({L_{x}}/{2}\right)\right)}\right|_{v_{\mathrm{bf}}=0}.$
(35)

The phase step of a soliton is a characteristic feature of solitons in
superfluids and is expected to depend on the soliton velocity. The top panel
of Fig. 2 shows results for two-dimensional systems for several binding
energies. Interestingly, for low binding energies we observe a dynamic change
of the shape of the relation $\phi_{\mathrm{s}}(v_{\mathrm{s}})$, which
eventually leads to the appearance of three distinct solitonic solutions for
the same velocity. It is also remarkable that for sufficiently high velocities
all phase-step values seem to lie on one universal curve.

The bottom panel of Fig. 2 shows analogous results for one-dimensional
systems, obtained by restricting our calculation to a single transverse
momentum $p_{\mathrm{t}}=0$. The numerical results are compared with the
theoretical prediction from Ref. [13],
$\phi_{\mathrm{s}}=2\arccos(v_{\mathrm{s}}/v_{\mathrm{pb}})$ and show a
qualitative agreement with this result. Note that Ref. [13] applied additional
approximations to the BdG equations like a linearization of the fermion
dispersion relation, which are not used in our calculations.

The differences between the situation in two and one dimensions as shown in
Fig. 2 are striking. A further analysis of the change in the shape of the
$\phi_{\mathrm{s}}(v_{\mathrm{s}})$ relation will be given in the next
section.

##  IV Dispersion relations and inertial masses

Further we focus on the soliton masses and their dependence on the binding
energy and velocity. The physical mass defined by Eq. (32) can be also
obtained from the depletion number [34]:

$\displaystyle M_{\mathrm{P}}=N_{\mathrm{d}}m.$ (36)

The inertial mass can be calculated using Eq. (33). However, this equation
introduces significant numerical errors for $v_{\mathrm{s}}\approx 0$. To
alleviate this we are using the equivalent relation
$M_{\mathrm{I}}=\nicefrac{{\partial P_{\mathrm{c}}}}{{\partial
v_{\mathrm{s}}}}$ (derived using: $v_{\mathrm{s}}=\nicefrac{{\partial\Delta
F}}{{\partial P_{\mathrm{c}}}}$). This equation can be further transformed
into

$\displaystyle M_{\mathrm{I}}=\frac{\partial P_{\mathrm{c}}}{\partial
v_{\mathrm{s}}}=\frac{\partial P_{\mathrm{ph}}}{\partial
v_{\mathrm{s}}}+\frac{\partial P_{\phi}}{\partial
v_{\mathrm{s}}}=M_{\mathrm{P}}+mL_{y}\frac{\partial\phi_{\mathrm{s}}}{\partial
v_{\mathrm{s}}}.$ (37)

It shows that the difference between the inertial and the physical mass is
proportional to the slope of the $\phi_{\mathrm{s}}(v_{\mathrm{s}})$ relation
presented in Fig. 2.

The soliton energy – momentum dispersion relations are presented for a range
of binding energies in Fig. 3, using the same color scheme as in Fig. 2. The
quantities $\Delta F_{\mathrm{s}}$, $P_{\mathrm{c}}$, $M_{\mathrm{I}}$ and
$M_{\mathrm{P}}$ are extensive in $y$ due to the homogeneity of the system in
this direction. Thus we are plotting the scaled, intensive values $\Delta
F_{\mathrm{s}}/L_{y}$, $P_{\mathrm{c}}/L_{y}$, $M_{\mathrm{I}}/L_{y}$ and
$M_{\mathrm{P}}/L_{y}$.

Figure 3: Soliton dispersion relations. Free energy of the soliton as a
function of canonical momentum (scaled with transverse system length). Results
for
$E_{\mathrm{b}}=0.009\varepsilon_{\mathrm{F}},0.03\varepsilon_{\mathrm{F}},0.05\varepsilon_{\mathrm{F}}$
are shown with the dashed blue, dotted yellow and dash-dotted green lines,
respectively. Insets present closeups on the features of the
$E_{\mathrm{b}}=0.03\varepsilon_{\mathrm{F}}$ and
$E_{\mathrm{b}}=0.009\varepsilon_{\mathrm{F}}$ curves.

As one can see in Fig. 3, there are qualitative differences between the
presented dispersion curves. For $E_{\mathrm{b}}=0.05\varepsilon_{\mathrm{F}}$
we have a smooth inverted lobe with negative curvature. For
$E_{\mathrm{b}}=0.03\varepsilon_{\mathrm{F}}$ the inverted lobe is split in
the center by two cusps. The outer arms retains a negative curvature but the
central part is convex. For even smaller
$E_{\mathrm{b}}=0.009\varepsilon_{\mathrm{F}}$, the curvature of the central
part changes. In the very center the curvature is negative while it is
positive at the sides. There are two inflection points where the curvature
vanishes and the inertial mass diverges according to (33).

To better understand the latter, most complicated, case we plot the inertial
mass as a function of the soliton velocity in Fig. 4. The central branch
(depicted with solid blue line) is symmetric around $v_{\mathrm{s}}=0$. It
always stays negative and diverges to $-\infty$ for
$|v_{\mathrm{s}}|\rightarrow v_{\mathrm{c}}$, where $v_{\mathrm{c}}$ is the
critical velocity at which the curvature of the dispersion relation vanishes.

Figure 4: Inertial mass for $E_{\mathrm{b}}=0.009\varepsilon_{\mathrm{F}}$ as
function of the soliton velocity. The blue solid line is the central branch,
which contains the $\pi$-phase stationary soliton. The green and purple dashed
lines are the outer branches containing non-$\pi$-phase stationary solutions.
Yellow vertical dotted lines denotes the critical velocity $v_{\mathrm{c}}$ at
which the inertial mass diverges.

One of the outer branches (green line) starts at $-v_{\mathrm{c}}$ where it
has $+\infty$ singularity, then it monotonically decreases, eventually
crossing $0$ (which is not well visible on the plot). The other is a mirror
image of the first with opposite velocities.

As there is a particular interest in slowly moving solitons, we consider the
stationary ($v_{\mathrm{s}}=0$) solutions specifically, and plot their masses
for a range of binding energies in Fig. 5. The inertial mass is shown in panel
Fig. 5(a) and the physical mass in panel Fig. 5(b).

Figure 5: Soliton mass as function of the pair binding energy. (a) Inertial
and (b) physical mass per transverse length for stationary
($v_{\mathrm{s}}=0$) solitons for a range of binding energies. The blue solid
line shows data for the soliton solution with $\pi$-phase step, while the
orange dashed line shows the mass for the two degenerate soliton solutions
with non-$\pi$ phase step. The dotted magenta line marks the value of
divergence $E_{\mathrm{b}}^{\infty}=0.01025\pm
0.0005\varepsilon_{\mathrm{F}}$. The inset presents a closeup of the two upper
curves to show that they both are crossing $M_{\mathrm{I}}=0$ at some point.

Analyzing the data presented in Fig. 5 we can mark three transition points: At
$E_{\mathrm{b}}^{0}=0.034\pm 0.0005\varepsilon_{\mathrm{F}}$ the inertial mass
vanishes for the sole $v_{\mathrm{s}}=0$ soliton solution (with
$\phi_{\mathrm{s}}=\pi$). At the smaller binding energy
$E_{\mathrm{b}}^{\infty}=0.01025\pm 0.0005\varepsilon_{\mathrm{F}}$ signifying
weaker interactions, the effective mass of all $v_{\mathrm{s}}=0$ soliton
solutions diverge. Finally at even weaker interactions with
$E_{\mathrm{b}}^{0\prime}=0.006\pm 0.001\varepsilon_{\mathrm{F}}$ a third
transition point is found where the inertial mass vanishes for two $v_{s}=0$
solutions with non-$\pi$ phase step.

Looking from the right-hand side, for high binding energies
$E_{\mathrm{b}}>E_{\mathrm{b}}^{0}$ the inertial mass is negative. Below
$E_{\mathrm{b}}^{0}$ it becomes positive and grows towards infinity
approaching $E_{\mathrm{b}}\rightarrow E_{\mathrm{b}}^{\infty(+)}$. On the
other side of this singularity for $E_{\mathrm{b}}<E_{\mathrm{b}}^{\infty}$ we
have three distinct solutions with stationary solitons. One branch associated
with $\pi$-phase solutions diverges to $-\infty$ for
$E_{\mathrm{b}}^{\infty(-)}$, while for lower $E_{\mathrm{b}}$ values it is
finite and negative. The other two degenerate curves are associated with
non-$\pi$-phase solutions. These solutions only exist for
$E_{\mathrm{b}}<E_{\mathrm{b}}^{\infty}$ and their inertial mass diverges to
$+\infty$ at the termination point $E_{\mathrm{b}}^{\infty}$. Further
decreasing $E_{\mathrm{b}}$, the mass of the non-$\pi$-phase solitons
decreases monotonically eventually becoming negative for
$E_{\mathrm{b}}<E_{\mathrm{b}}^{0\prime}$.

The physical mass shown on the lower panel does not exhibit as dramatic
changes as the inertial mass. For all binding energies in our numerical data
it remains negative and decreases monotonically with increasing
$E_{\mathrm{b}}$. On the BEC side of the crossover it decreases, while towards
the BCS regime it increases towards zero, which is consistent with the
observation that we do not observe any change in the particle number density
caused by the presence of the soliton in that limit. No noticeable abrupt
change of the physical mass can be seen upon crossing the threshold
$E_{\mathrm{b}}^{0}$. Below $E_{\mathrm{b}}^{\infty}$ three different branches
appear but differences in $M_{\mathrm{P}}$ between the $\pi$-phase solution
and the others are orders of magnitude smaller than the actual
$M_{\mathrm{P}}$ values.

In Fig. 6 we present the ratio of the inertial and the physical mass as a
function of the binding energy. Due to the qualitatively simple dependence of
$M_{\mathrm{P}}$ on $E_{\mathrm{b}}$, the behavior of the ratio
$M_{\mathrm{I}}/M_{\mathrm{P}}$ is governed mostly by
$M_{\mathrm{I}}(v_{\mathrm{s}})$. An interesting takeaway is that the ratio
appears to correctly approach the value of 2 for large $E_{\mathrm{b}}$ in
agreement with predictions from Gross-Pitaevskii theory for a BEC of strongly-
bound bosons [39, 33], while we know that BdG mean-field theory in two
dimensions fails to correctly predict the equation of state in the BEC regime
[29].

Figure 6: Ratio of inertial and physical mass for stationary
($v_{\mathrm{s}}=0$) solitons for a range of binding energies. The inset shows
the only branch for $E_{\mathrm{b}}>E_{\mathrm{b}}^{0}$ for a wider range of
binding energies. The black horizontal line is the known value for BEC limit
$M^{\mathrm{BEC}}_{\mathrm{I}}/M^{\mathrm{BEC}}_{\mathrm{P}}=2$ [39, 33].
Otherwise the line styles and color schemes are as in Fig. 5.

Before summarizing the three different physical regimes, we list a few general
rules that are useful when interpreting the presented data:

  * •

The value of the inertial mass $M_{\mathrm{I}}$ is largely determined by the
slope of
$P_{\mathrm{c}}(v_{\mathrm{s}})\propto-\phi_{\mathrm{s}}(v_{\mathrm{s}})$ and
can be alternatively obtained from the inverse curvature of the dispersion
relation $\Delta F_{\mathrm{s}}(P_{\mathrm{c}})$.

  * •

The soliton velocity $v_{\mathrm{s}}$ is determined by the slope of the
dispersion relation through $v_{\mathrm{s}}=d\Delta
F_{\mathrm{s}}(P_{\mathrm{c}})/dP_{\mathrm{c}}$.

  * •

Therefore, smooth extrema of the dispersion relation mark points where
$v_{\mathrm{s}}=0$, i.e. stationary soliton solutions exist.

  * •

When the dispersion relation has a cusp the slope of the relation
$P_{\mathrm{c}}(v_{\mathrm{s}})$ vanishes and consequently $M_{\mathrm{I}}=0$.

  * •

When the dispersion relation has an inflection point the curve
$P_{\mathrm{c}}(v_{\mathrm{s}})$ has a vertical asymptote and the inertial
mass $M_{\mathrm{I}}$ diverges.

_Tightly-bound fermion pairs:_ For $E_{\mathrm{b}}>E_{\mathrm{b}}^{0}$ (green
dash-dotted line in Figs 2 and 3) the canonical momentum is monotonous as a
function of the soliton velocity ($v_{\mathrm{s}}$). The dispersion relation
is an inverted lobe and the inertial mass is always negative. Qualitatively,
the behavior of systems in this regime resembles the one in the BEC limit.

When $E_{\mathrm{b}}=E_{\mathrm{b}}^{0}$, the relation
$P_{\mathrm{c}}(v_{\mathrm{s}})$ flattens for $v_{\mathrm{s}}=0$ and a cusp
appears in the center of the dispersion relation, which means that
$M_{\mathrm{I}}(v_{\mathrm{s}}=0)=0$.

_Simple swallow-tail regime with positive inertial mass:_ For
$E_{\mathrm{b}}^{0}>E_{\mathrm{b}}>E_{\mathrm{b}}^{\infty}$ (yellow dotted
line in Figs 2 and 3), the dispersion relation has a swallow-tail shape. It is
split into three branches by two cusps. The outer branches are symmetric and
have negative curvature while the central one has positive curvature. The
slope at the cusp defines a threshold velocity $v_{\mathrm{t}}$ at which the
sign of $M_{\mathrm{I}}$ changes. The canonical momentum decreases with
growing $v_{\mathrm{s}}$ in the central part (corresponding to the central
branch of the dispersion relation) and increases for
$|v_{\mathrm{s}}|>v_{\mathrm{t}}$. The inertial mass is positive for the
central branch and negative for the rest.

At the point $E_{\mathrm{b}}=E_{\mathrm{b}}^{\infty}$ the dispersion relation
retains a swallow-tail shape but its curvature drops to $0$ at
$v_{\mathrm{s}}=0$. The $P_{\mathrm{c}}(v_{\mathrm{s}})$ relation acquires a
vertical asymptote. The inertial mass is positive for
$|v_{\mathrm{s}}|<v_{\mathrm{t}}$ and grows to indefinitely as
$v_{\mathrm{s}}$ approaches $0$.

_Complex swallow-tail regime:_ For $E_{\mathrm{b}}<E_{\mathrm{b}}^{\infty}$
(blue dashed line in Figs 2 and 3), the central branch in the dispersion
relation has a maximum in the center, indicating a stationary soliton solution
with negative inertial mass. Two other symmetric extrema exist. Depending on
the value of $E_{\mathrm{b}}$ they appear in the central or the outer branches
and are minima or maxima respectively. The $P_{\mathrm{c}}(v_{\mathrm{s}})$
relation becomes three-valued in the central part. This leads to the
appearance of a range of velocities around $v_{\mathrm{s}}=0$ for which there
are three distinct solutions. Only one $v_{\mathrm{s}}=0$ solution has
$\phi_{\mathrm{s}}=\pi$ while the other two have phase steps that are
symmetric around $\pi$. This leads also to existence of three separated
branches of $M_{\mathrm{I}}(v_{\mathrm{s}})$ as seen in Fig. 4.

The rich scenario for the inertial mass with zero crossings and divergences in
the same dispersion relation as seen in Fig. 4 has interesting consequences
for the soliton dynamics in a trapped Fermi gas. It is important to note,
though, that the Newton-like equations for quasiparticle-like soliton dynamics
(31) with the mass parameters of Eqs. (32) and (33) only hold for a slowly
varying potential and assumes the adiabatic adjustment of the soliton to its
environment. In the case of rapid changes of the inertial mass under small
changes of velocity these conditions may be compromised and the quasiparticle
picture of the soliton dynamics may break down.

An important physical consequence of the sign of the effective mass is that it
predicts the presence (or absence) of the snaking instability. The snaking
instability is a dynamical instability of a soliton solution in two or more
dimensions, where an initially homogeneous solution (a soliton stripe) in the
transverse direction spontaneously disintegrates through a snaking process
[40, 41, 42, 43]. The snaking instability can occur when the inertial mass
$M_{\mathrm{I}}$ is negative, which leads to unstable small-amplitude
oscillations of sinusoidal modulations of the soliton shape in the presence of
a restoring force provided by the surface tension coming from the energy
density of the soliton $\Delta F_{\mathrm{s}}/L_{y}$ [44]. The snaking
instability is not present when the inertial mass is positive, in which case
stable small amplitude oscillations are obtained. Thus, the regions of
negative inertial mass in Fig. 4 indicate unstable solutions under the snaking
instability, while positive $M_{\mathrm{I}}$ solutions are stable.

##  V Counterflow in momentum space

The swallow-tail scenario for the dispersion relations of a soliton in a two-
dimensional Fermi superfluid presents an unexpected behavior and has no
analogy in bosonic superfluids. While swallow-tail dispersions are known to
occur for dark solitons in a BEC in the presence of a periodic potential [45,
46], or for transversely modulated vortex-ring solutions [47, 34], they are
not present for transversely homogeneous solutions in a flat potential [7].
Previous studies of fermionic superfluids with the BdG equation in one [13] or
three dimensions [17, 18, 19] also have not shown any such effects. The
numerical studies in three dimension have not explored the deep BCS regime,
however, due to the numerical challenges posed by resolving the diverging
length scales of the soliton envelope (the coherence length) and the Fermi
length scale.

In order to properly understand what is happening it would be convenient to
find analytic solutions of solitons in the BdG equations. Such solitons are
available in a one-dimensional setting [13], where, however, smooth dispersion
relations with a very simple structure were predicted (as shown in Fig. 2).
The one-dimensional analytical solutions, which rely on linearized single-
particle dispersions, cannot be simply extended to two spatial dimensions due
to the contributions from transverse momenta.

To better understand source of this difference and its connection with higher
dimensionality of our problem we will consider the case of
$E_{\mathrm{b}}=0.005\varepsilon_{\mathrm{F}}$ where three stationary (i.e.
$v_{\mathrm{s}}=0$) solutions exist with different phase steps.

Figure 7: Complex phases of $\Delta(x)$ for solitonic self-consistent
solutions. All three results are stationary and are obtained for the same
value of $E_{\mathrm{b}}=0.005\varepsilon_{\mathrm{F}}$. The blue solid line
is a soliton with $\phi_{\mathrm{s}}=\pi$. The red dashed and green dotted
lines are results with non-$\pi$ phases.

In Fig. 7 we show the complex phase of $\Delta(x)$ of the solitonic solutions.
Two qualitatively different behaviors can be observed. The $\pi$-phase
solution simply has a constant phase with a discontinuous phase jump of
exactly $\pi$. The solution can thus be presented as a purely real function.
The other two solutions, however, are irreducibly complex with a continuous
change of phase. One of these solutions is the complex conjugate of the other.

It is known that for bosonic superfluids the gradient of the phase is directly
connected with the superfluid velocity, which is proportional to the current
by the particle number density [48]. The same connection can also be made for
fermionic superfluids in case of a homogeneous flow. It is obviously not the
case here, as we observe nonzero phase gradients for the stationary soliton
solutions.

Figure 8: Current densities for different transverse momenta for
$E_{\mathrm{b}}=0.005\varepsilon_{\mathrm{F}}$, $v_{\mathrm{s}}=0$,
$\phi_{\mathrm{s}}=0.82\pi$. In the inset current values integrated over $x$
are presented as function of $p_{\mathrm{t}}$ are plotted as dashed blue line.
The solid red line is a cumulative sum for momenta from $-p_{\mathrm{t}}$ to
$p_{\mathrm{t}}$ (currents are symmetric in transverse momenta).

In order to understand this discrepancy let us recall that the soliton
solutions are inhomogeneous only in one dimension but a second, homogeneous,
transverse direction is present as well. In the BdG formalism the transverse
direction is included by summing multiple one-dimensional solutions with
different transverse momenta $p_{\mathrm{t}}$, as in Eq. (17) to obtain the
complex order parameter
$\Delta(x)=\sum_{p_{\mathrm{t}}}\Delta_{p_{\mathrm{t}}}(x)$. One can calculate
the partial contributions $\Delta_{p_{\mathrm{t}}}(x)$ as well as the
transverse contributions $j_{p_{\mathrm{t}}}(x)$ to the current density, where
$j(x)=\sum_{p_{\mathrm{t}}}j_{p_{\mathrm{t}}}(x)$ for each transverse momentum
channel:

$\displaystyle\Delta_{p_{\mathrm{t}}}(x)$
$\displaystyle=-g\sum_{\nu}^{\prime}u_{\nu}^{(p_{\mathrm{t}})}(x)v_{\nu}^{(p_{\mathrm{t}})*}(x),$
(38) $\displaystyle j_{p_{\mathrm{t}}}(x)$
$\displaystyle=-2\hbar\sum_{\nu}v_{\nu}^{(p_{\mathrm{t}})*}(x)\partial_{x}v_{\nu}^{(p_{\mathrm{t}})}(x).$
(39)

Inspecting the different transverse contributions reveals that the partial
current densities $j_{p_{\mathrm{t}}}(x)$ have different signs depending on
$p_{\mathrm{t}}$. Consequently, the current contributions have different
directions. This can be seen in Fig. 8. In the special case we are considering
now – the stationary soliton – all the flows add up to zero. This can be seen
in the inset of Fig. 8. The blue dashed line represents the integrated current
densities as a function of the transverse momentum. One can see that it
changes sign – for small $|p_{\mathrm{t}}|$ we observe flow in one direction
while for larger $|p_{\mathrm{t}}|$ the direction reverses. The red line
represents a cumulative sum. It goes to zero for large transverse momentum,
which demonstrates that eventually the total current vanishes as we should
expect for a stationary soliton.

The example presented above showed that stationary soliton solutions with
finite phase gradients and a non-$\pi$ phase-step can exist and are supported
by a cancellation of the partial flows for different transverse momenta.
However, the coexistence of different partial flow directions is not limited
to the specific soliton solution that was presented. Such coexistence of flows
with different directions is found also for moving solitons and for a much
wider range of binding energies, up to $E_{\mathrm{b}}\approx
0.2\varepsilon_{\mathrm{F}}$. In Fig. 9 we show the approximate threshold of
the transverse momentum $p^{(c)}_{\mathrm{t}}$ at which the transition between
different flow directions occurs for a range of velocities and for three
different binding energies.

Figure 9: Sign changing current contributions in the soliton solutions to the
BdG equations. Shown is the threshold transverse momentum
$p^{(c)}_{\mathrm{t}}$ (in units of the Fermi momentum) above which the
current contribution changes sign. It is presented as a function of the
soliton velocity (scaled with the pair breaking velocity) for three different
values of the binding energies. The inset shows a closeup for
$E_{\mathrm{b}}=0.009\varepsilon_{\mathrm{F}}$ and for small $v_{s}$.

##  VI Conclusions

We have presented numerical simulations of solitons in a two-dimensional Fermi
superfluid. The results are obtained as self-consistent solutions from the
Bogoliubov-de Gennes method. For the small binding energies on the far BCS
side of the crossover regime we have observed unexpected and previously unseen
behavior of the solitons.

Below a threshold binding energy of $E_{\mathrm{b}}^{0}=0.034\pm
0.0005\varepsilon_{\mathrm{F}}$ the dispersion relations change shape from
smooth lobes to swallow tail type, which also leads to a change of sign of the
inertial mass of the soliton. Below the even lower threshold of
$E_{\mathrm{b}}^{\infty}=0.01025\pm 0.0005\varepsilon_{\mathrm{F}}$ we observe
the coexistence of several solutions for the same velocity, as well as the
existence of velocities for which the inertial mass diverges to infinity.
Those phenomena are connected with the fact that different one-dimensional
solutions contributing to the full soliton shape contribute partial superfluid
currents flowing in opposite directions. This allows, for example, a situation
where a soliton with a nontrivial complex phase profile can be stationary, as
flows for different transverse momenta cancel each other.

The soliton solutions studied in this work, if realized experimentally, can
show a range of interesting properties. An example is the coexistence of
multiple solutions with the same velocity that react differently to an
external potential gradient. Another example are non-standard oscillations in
a harmonic external potential due to a change of sign of the inertial mass, or
a sudden breakdown of solitons when they approach a velocity where the
inertial mass vanishes.

In the model studied in this paper the transverse direction is assumed to be
homogeneous. However, a similar setup could be obtained for ultracold atoms in
a flat-bottom trap. This might allow for an experimental detection of the
current components with different signs for different transverse momenta
predicted by our numerical simulations. The different transverse momentum
components can be separated in an experiment by time-of-flight position
measurements after turning off the trap in the $y$ direction. Previously
overlapping components with different $p_{\mathrm{t}}$ will be separated. Then
different directions of flow of the different components could be observed,
which may be a method to experimentally confirm this phenomenon.

## Acknowledgments

The authors thank Uli Zülicke for stimulating discussions. J.M. wants to
acknowledge support of National Science Centre (Poland) via Sonatina grant
2020/36/C/ST2/00494. J.B.’s work was supported by the Marsden Fund of New
Zealand, Contract No. MFP-MAU2007 from government funding managed by the Royal
Society of New Zealand Te Apārangi.

##  Appendix A Chemical potential of moving superfluid

Let us consider a frame transformation where we move with velocity
$v_{\mathrm{fr}}$ in the $x$ direction. It is important to note that we have
to change the value of the chemical potential $\mu$, if we want to preserve
the density in the moving frame of reference. The reason for this can be
easily understood considering the simpler model of a homogeneous fermionic
superfluid in a box with periodic boundary conditions carrying a ring-current.
For given binding energy $E_{\mathrm{b}}$ (parameterizing the coupling
strength) we can self consistently find $\Delta_{0}$. Constraining the density
to be constant fixes the chemical potential $\mu$, which also has to be
determined self consistently. For a stationary superfluid in two dimensions
the values are given by
$\Delta_{0}=\sqrt{2E_{\mathrm{b}}\varepsilon_{\mathrm{F}}}$ and
$\mu_{0}=\varepsilon_{\mathrm{F}}-E_{\mathrm{b}}/2$ [35].

A ring current can be introduced by adding a linear phase change to
$\Delta(\mathbf{r})=\Delta e^{ik_{\mathrm{fr}}x}$ (where
$2mv_{\mathrm{fr}}/\hbar=k_{\mathrm{fr}}=2\pi n/L_{x},\;n\in\mathbb{N}$, and
$L_{x}$ is the length of the box in $x$ direction). After absorbing the phase
gradient in equal portions in the $u$ and $v$ quasiparticle amplitudes, the
BdG matrix of such a system can be represented in momentum space as

$\displaystyle H_{k}=\left[\begin{array}[]{c
c}\xi_{{k}+k_{\mathrm{fr}}/2}&\Delta\\\
\Delta^{*}&-\xi_{{k}-k_{\mathrm{fr}}/2}\end{array}\right],$ (42)

where we have suppressed the transverse momentum indices (formally set
$p_{\mathrm{t}}=0$) for simplicity. As the one particle energy reads
$\xi_{{k}}=\frac{\hbar^{2}}{2m}k^{2}-\mu$, we have

$\displaystyle\xi_{{k\pm
k_{\mathrm{fr}}/2}}=\frac{\hbar^{2}}{2m}k^{2}\pm\frac{\hbar^{2}}{2m}kk_{\mathrm{fr}}+\frac{\hbar^{2}}{2m}\frac{k_{\mathrm{fr}}^{2}}{4}-\mu.$
(43)

The quadratic term can be grouped together with chemical potential

$\displaystyle\tilde{\mu}=\mu-\frac{\hbar}{2m}\frac{{k_{\mathrm{fr}}}^{2}}{4},$
(44)

while mixed term $kk_{\mathrm{fr}}$ has the same sign in both particle and
hole segments. As an energy shift it affects only the quasiparticle energies
but not eigenstates. As long as this energy shift is smaller than the half gap
$\Delta/2$ it will not change quasiparticle level occupation numbers. Thus we
have

$\displaystyle\varepsilon_{{k}}$
$\displaystyle=\pm\left(\frac{\hbar^{2}}{2m}{k}k_{\mathrm{fr}}+\sqrt{\tilde{\xi}_{{k}}^{2}+\Delta^{2}}\right),$
(45) $\displaystyle u_{{k}+k_{\mathrm{fr}}/2}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(1+\frac{\tilde{\xi}_{{k}}}{\sqrt{\tilde{\xi}_{{k}}^{2}+|\Delta|^{2}}}\right)^{1/2},$
(46) $\displaystyle v_{{k}-k_{\mathrm{fr}}/2}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(1-\frac{\tilde{\xi}_{{k}}}{\sqrt{\tilde{\xi}_{{k}}^{2}+|\Delta|^{2}}}\right)^{1/2},$
(47)

where

$\displaystyle\tilde{\xi}_{{k}}=\frac{\hbar^{2}}{2m}{k}^{2}-\tilde{\mu}.$ (48)

Note that the particle number density is given by Eq. (18). We choose to keep
it fixed and determine the chemical potential $\mu$ self-consistently. As the
equations (46) and (47) have the same form as equations without the current, a
self-consistent solution will yield the same values:

$\displaystyle\Delta$ $\displaystyle=\Delta_{0},$ (49)
$\displaystyle\tilde{\mu}$ $\displaystyle=\mu_{0}.$ (50)

Going back to the original variables:

$\displaystyle\mu$
$\displaystyle=\tilde{\mu}+\frac{\hbar^{2}}{2m}{k_{\mathrm{fr}}}^{2},$ (51)

where $\mu$ is the chemical potential in the frame co-moving with the ring
current, i.e. where no explicit current is present. $\tilde{\mu}$ is the
chemical potential in the frame where a current with the superfluid velocity
$v_{\mathrm{fr}}=\hbar k_{\mathrm{fr}}/2m$ is present. So if we want to keep
the same density in cases with different ring-currents (but also with backflow
caused by the presence of soliton) we have to change chemical potential in
accordance with Eq. (51). Another way to look at this phenomenon is to
consider Fermi spheres for a moving superfluid. Due to a velocity-induced
shift between particle and hole lobes, if the chemical potential remained
constant, size of the Fermi sphere would be smaller, which would lead to a
smaller density. To preserve the density, the Fermi sphere size has to be
increased by adjusting the chemical potential.

##  Appendix B Observable transformation under Galilean boost

In the main part of this article self-consistent solutions are obtained in the
frame co-moving with the soliton. This stems from the construction of the
time-independent method. The observables calculated using those results are
also in the co-moving frame. However, to analyze and compare those results,
this data have to be transformed back to the lab frame. For Galilean boost
energy and momentum transform as follows:

$\displaystyle P^{\mathrm{SF}}$ $\displaystyle=P-Mv_{\mathrm{s}},$ (52)
$\displaystyle E^{\mathrm{SF}}$
$\displaystyle=E+\frac{1}{2}Mv_{\mathrm{s}}^{2}-Pv_{\mathrm{s}},$ (53)

where we are using the same notation as in the main article: superscript
$\mathrm{(}SF)$ for soliton frame and no superscript for lab frame. As we are
working in the grand canonical ensemble instead of energy we are using free
energy of the system:

$\displaystyle F^{\mathrm{SF}}=E^{\mathrm{SF}}-\mu_{v_{\mathrm{s}}}N$ (54)

While density is kept constant, the chemical potential will change as shown in
App. A:

$\displaystyle\mu_{v_{\mathrm{s}}}=\mu_{0}+\frac{1}{2}mv_{\mathrm{s}}^{2},$
(55)

if we assume that in the lab frame we do not have backflow. Further:

$\displaystyle F^{\mathrm{SF}}+\mu_{v_{\mathrm{s}}}N$
$\displaystyle=F+\mu_{0}N+\frac{1}{2}Mv_{\mathrm{s}}^{2}-P_{c}v_{v_{\mathrm{s}}},$
$\displaystyle F^{\mathrm{SF}}+\mu_{0}N+\frac{1}{2}Mv_{v_{\mathrm{s}}}^{2}$
$\displaystyle=F_{0}+\mu_{0}N+\frac{1}{2}Mv_{\mathrm{s}}^{2}-P_{c}v_{\mathrm{s}},$
$\displaystyle F^{\mathrm{SF}}$ $\displaystyle=F-P_{c}v_{\mathrm{s}}.$ (56)

Free energy considered above is free energy of the whole system, with or
without the soliton ($F_{\mathrm{s}}$ or $F_{0}$). We are interested in the
free energy of the soliton itself:

$\displaystyle\Delta F=F_{\mathrm{s}}-F_{0}.$ (57)

Knowing that $P_{0}=0$ and $P_{\mathrm{s}}=M_{\mathrm{d}}v_{\mathrm{B}}$ where
$M_{\mathrm{d}}$ is depletion mass ($M_{\mathrm{d}}=m(N_{\mathrm{s}}-N_{0})$),
we get:

$\displaystyle\Delta F=\Delta
F_{v_{\mathrm{s}}}+M_{\mathrm{d}}v_{\mathrm{s}}^{2}.$ (58)

## References

  * Onnes [1911] H. K. Onnes, The resistance of pure mercury at helium temperatures, Commun. Phys. Lab. Univ. Leiden. Suppl. 29 (1911).
  * Kapitza [1938] P. Kapitza, Viscosity of Liquid Helium below the $\lambda$-Point, Nature 141, 74 (1938).
  * Allen and Misener [1938] J. F. Allen and A. D. Misener, Flow Phenomena in Liquid Helium II, Nature 142, 643 (1938).
  * Korteweg and de Vries [1895] D. D. J. Korteweg and D. G. de Vries, Xli. on the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 39, 422 (1895), https://doi.org/10.1080/14786449508620739 .
  * Muto  _et al._ [1988] V. Muto, J. Halding, P. Christiansen, and A. Scott, Solitons in DNA, J Biomol Struct Dyn. 5, 873 (1988).
  * Zurek [1996] W. Zurek, Cosmological experiments in condensed matter systems, Physics Reports 276, 177 (1996).
  * Tsuzuki [1971] T. Tsuzuki, Nonlinear waves in the Pitaevskii-Gross equation, J. Low Temp. Phys. 4, 441 (1971).
  * Burger  _et al._ [1999] S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, A. Sanpera, G. V. Shlyapnikov, and M. Lewenstein, Dark solitons in Bose-Einstein condensates, Phys. Rev. Lett. 83, 5198 (1999).
  * Denschlag  _et al._ [2000] J. Denschlag, J. E. Simsarian, D. L. Feder, C. W. Clark, L. A. Collins, J. Cubizolles, L. Deng, E. W. Hagley, K. Helmerson, W. P. Reinhardt, S. L. Rolston, B. I. Schneider, and W. D. Phillips, Generating Solitons by Phase Engineering of a Bose-Einstein Condensate, Science 287, 97 (2000).
  * Becker  _et al._ [2008] C. Becker, S. Stellmer, P. Soltan-Panahi, S. Dörscher, M. Baumert, E.-M. Richter, J. Kronjäger, K. Bongs, and K. Sengstock, Oscillations and interactions of dark and dark–bright solitons in Bose–Einstein condensates, Nat. Phys. 4, 496 (2008).
  * Ku  _et al._ [2016] M. J. H. Ku, B. Mukherjee, T. Yefsah, and M. W. Zwierlein, Cascade of Solitonic Excitations in a Superfluid Fermi gas: From Planar Solitons to Vortex Rings and Lines, Phys. Rev. Lett. 116, 045304 (2016), arxiv:1507.01047 .
  * Shamailov and Brand [2016] S. S. Shamailov and J. Brand, Dark-soliton-like excitations in the Yang–Gaudin gas of attractively interacting fermions, New Journal of Physics 18, 075004 (2016).
  * Efimkin and Galitski [2015] D. K. Efimkin and V. Galitski, Moving solitons in a one-dimensional fermionic superfluid, Phys. Rev. A 91, 023616 (2015).
  * Zou  _et al._ [2016] P. Zou, J. Brand, X.-J. Liu, and H. Hu, Traveling Majorana Solitons in a Low-Dimensional Spin-Orbit-Coupled Fermi Superfluid, Phys. Rev. Lett. 117, 225302 (2016), arxiv:1509.01803 .
  * Dziarmaga and Sacha [2004] J. Dziarmaga and K. Sacha, Soliton in BCS superfluid Fermi gas (2004), arxiv:cond-mat/0407585 .
  * Antezza  _et al._ [2007] M. Antezza, F. Dalfovo, L. P. Pitaevskii, and S. Stringari, Dark solitons in a superfluid Fermi gas, Phys. Rev. A 76, 043610 (2007).
  * Liao and Brand [2011] R. Liao and J. Brand, Traveling dark solitons in superfluid Fermi gases, Phys. Rev. A 83, 041604(R) (2011).
  * Scott  _et al._ [2011] R. Scott, F. Dalfovo, L. Pitaevskii, and S. Stringari, Dynamics of Dark Solitons in a Trapped Superfluid Fermi Gas, Phys. Rev. Lett. 106, 185301 (2011).
  * Spuntarelli  _et al._ [2011] A. Spuntarelli, L. D. Carr, P. Pieri, and G. C. Strinati, Gray solitons in a strongly interacting superfluid Fermi gas, New J. Phys. 13, 035010 (2011).
  * Scott  _et al._ [2012] R. G. Scott, F. Dalfovo, L. P. Pitaevskii, S. Stringari, O. Fialko, R. Liao, and J. Brand, The decay and collisions of dark solitons in superfluid Fermi gases, New J. Phys. 14, 023044 (2012), arxiv:1109.6444 .
  * Klimin  _et al._ [2014] S. N. Klimin, J. Tempere, and J. T. Devreese, Finite-temperature effective field theory for dark solitons in superfluid Fermi gases, Phys. Rev. A 90, 053613 (2014).
  * Van Alphen  _et al._ [2019] W. Van Alphen, H. Takeuchi, and J. Tempere, Crossover between snake instability and Josephson instability of dark solitons in superfluid Fermi gases, Phys. Rev. A 100, 023628 (2019).
  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, From quantum matter to high-temperature superconductivity in copper oxides, Nature 518, 179 (2015).
  * Sato and Ando [2017] M. Sato and Y. Ando, Topological superconductors: a review, Reports on Progress in Physics 80, 076501 (2017).
  * Challis  _et al._ [2007] K. J. Challis, R. J. Ballagh, and C. W. Gardiner, Bragg scattering of cooper pairs in an ultracold fermi gas, Phys. Rev. Lett. 98, 093002 (2007).
  * Bloch  _et al._ [2012] I. Bloch, J. Dalibard, and S. Nascimbène, Quantum simulations with ultracold quantum gases, Nature Physics 8, 267 (2012).
  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010).
  * Bulgac  _et al._ [2011] A. Bulgac, Y.-L. Luo, P. Magierski, K. J. Roche, and Y. Yu, Real-Time Dynamics of Quantized Vortices in a Unitary Fermi Superfluid, Science 332, 1288 (2011).
  * Bertaina and Giorgini [2011] G. Bertaina and S. Giorgini, Bcs-bec crossover in a two-dimensional fermi gas, Phys. Rev. Lett. 106, 110403 (2011).
  * De Gennes [1999] P. G. De Gennes,  _Superconductivity of Metals and Alloys_ (CRC Press, 1999).
  * Note [1] Due to $\Delta(x)$ being a sum of products of the eigenstates of the BdG matrix, phase twists applied to either diagonal sector of BdG matrix ($\phi_{u}$ and $\phi_{v}$, respectively) add in contributing to $\Delta(x)$, while the individual phases do not have physical significance. In our calculations we usually set $\phi_{u}$ to $0$ and $\phi_{v}$ to $\pi$.
  * Shamailov and Brand [2019] S. S. Shamailov and J. Brand, Quantum dark solitons in the one-dimensional Bose gas, Phys. Rev. A 99, 043632 (2019).
  * Konotop and Pitaevskii [2004] V. V. Konotop and L. Pitaevskii, Landau Dynamics of a Grey Soliton in a Trapped Condensate, Phys. Rev. Lett. 93, 240403 (2004).
  * Mateo and Brand [2015] A. M. Mateo and J. Brand, Stability and dispersion relations of three-dimensional solitary waves in trapped Bose–Einstein condensates, New Journal of Physics 17, 125013 (2015).
  * Randeria  _et al._ [1989] M. Randeria, J.-M. Duan, and L.-Y. Shieh, Bound states, cooper pairing, and bose condensation in two dimensions, Phys. Rev. Lett. 62, 981 (1989).
  * Broyden [1965] C. G. Broyden, A class of methods for solving nonlinear simultaneous equations, Math. Comp. 19, 577 (1965).
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing, SIAM review 59, 65 (2017).
  * Anderson  _et al._ [1999] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK Users’ Guide_ , 3rd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999).
  * Busch and Anglin [2000] Th. Busch and J. R. Anglin, Motion of Dark Solitons in Trapped Bose-Einstein Condensates, Phys. Rev. Lett. 84, 2298 (2000), 11018869 .
  * Kuznetsov and Turitsyn [1988] E. A. Kuznetsov and S. K. Turitsyn, Instability and collapse of solitons in media with a defocusing nonlinearity, Sov. Phys. JETP 67, 1583 (1988).
  * Muryshev  _et al._ [1999] A. E. Muryshev, H. B. van Linden van den Heuvell, and G. V. Shlyapnikov, Stability of standing matter waves in a trap, Phys. Rev. A 60, R2665 (1999).
  * Brand and Reinhardt [2002] J. Brand and W. P. Reinhardt, Solitonic vortices and the fundamental modes of the "snake instability": Possibility of observation in the gaseous Bose-Einstein condensate, Phys. Rev. A 65, 043612 (2002).
  * Cetoli  _et al._ [2013] A. Cetoli, J. Brand, R. G. Scott, F. Dalfovo, and L. P. Pitaevskii, Snake instability of dark solitons in fermionic superfluids, Phys. Rev. A 88, 043639 (2013), arxiv:1307.3717 .
  * Kamchatnov and Pitaevskii [2008] A. Kamchatnov and L. Pitaevskii, Stabilization of Solitons Generated by a Supersonic Flow of Bose-Einstein Condensate Past an Obstacle, Phys. Rev. Lett. 100, 160402 (2008).
  * Mueller [2002] E. J. Mueller, Superfluidity and mean-field energy loops: Hysteretic behavior in Bose-Einstein condensates, Phys. Rev. A 66, 063603 (2002).
  * Muñoz Mateo  _et al._ [2019] A. Muñoz Mateo, V. Delgado, M. Guilleumas, R. Mayol, and J. Brand, Nonlinear waves of Bose-Einstein condensates in rotating ring-lattice potentials, Phys. Rev. A 99, 023630 (2019).
  * Komineas and Papanicolaou [2003] S. Komineas and N. Papanicolaou, Solitons, solitonic vortices, and vortex rings in a confined Bose-Einstein condensate, Phys. Rev. A 68, 043617 (2003).
  * Lan [1965] 46 - The theory of superfluidity of helium II, in  _Collected Papers of L.D. Landau_, edited by D. ter Haar (Pergamon, 1965) pp. 301–330.