# Motion of dark solitons in a non-uniform flow of Bose-Einstein condensate

S. K. Ivanov  A. M. Kamchatnov  Institute of Spectroscopy, Russian Academy of
Sciences, Troitsk, Moscow, 108840, Russia<EMAIL_ADDRESS>
###### Abstract

We study motion of dark solitons in a non-uniform one-dimensional flow of
Bose-Einstein condensate. Our approach is based on Hamiltonian mechanics
applied to the particle-like behavior of dark solitons in a slightly non-
uniform and slowly changing surrounding. In one-dimensional geometry, the
condensate’s wave function undergoes the jump-like behavior across the soliton
and this leads to generation of the counterflow in the background condensate.
For correct description of soliton’s dynamics, the contributions of this
counterflow to the momentum and energy of the soliton are taken into account.
The resulting Hamilton equations are reduced to the Newton-like equation for
the soliton’s path and this Newton equation is solved in several typical
situations. The analytical results are confirmed by numerical calculations.

###### pacs:

05.45.Yv, 47.35.Fg, 03.75.-b, 45.20.Jj

> The idea that localized solitary waves behave in external fields like point
> particles of the classical mechanics is very old and it was studied and used
> in a number of articles. However, its application to dynamics of dark
> solitons has some peculiarities which are not still fully understood. For
> example, from naive point of view a dark soliton in a cloud of Bose-Einstein
> condensate (BEC) is just a dip in the density distribution and,
> consequently, its motion must result in the motion of the center of mass in
> the opposite direction. In reality both experiments and numerical solutions
> of the Gross-Pitaevskii equation demonstrate different laws of motion of the
> dark soliton and the center of mass of the cloud of a quasi-one-dimensional
> condensate confined in a harmonic trap. This means that the motion of dark
> solitons is accompanied by some counterflow in the BEC cloud which
> contributes to the location of its center of mass. Correspondingly, the
> contributions of the counterflow to the momentum and the energy of dark
> solitons must be taken into account. Previously, this idea was applied to
> description of dark solitons motion in a quiescent condensate’s cloud. In
> this paper, we extend it to the general case of slightly non-uniform and
> slowly changing background distributions. The resulting Hamilton equations
> for the soliton’s motion are reduced to the Newton-like equation and this
> equation is solved for several typical situations of the background
> evolutions, such as a rarefaction wave, a self-similar expansion of the BEC
> cloud, a hydraulic flow of BEC past a wide obstacle. Solutions of the Newton
> equation agree very well with numerical solutions of the full Gross-
> Pitaevskii equation.

##  I Introduction

The name of solitons was coined by Zabusky and Kruskal [1] due to their
particle-like behavior: two solitons interact elastically and pass through one
another without losing their identity. This concept of solitons as particle-
like nonlinear wave excitations was confirmed and amplified by the properties
of solitons propagation in a non-uniform external field (see, e.g., [2, 3, 4,
5]). In particular, bright soliton clouds of attractive BEC oscillate in the
external trap potential $U=m\omega_{0}^{2}x^{2}/2$ ($m$ is an atomic mass)
with the trap frequency $\omega_{0}$. However, behavior of dark solitons in
repulsive BEC has some peculiarities. For example, in case of BEC confined in
a quasi-one-dimensional (1D) harmonic potential trap, they propagate along a
stationary Thomas-Fermi profile of the density and oscillate with the
frequency $\omega_{0}/\sqrt{2}$, although the center of mass of the condensate
still oscillates with the trap frequency $\omega_{0}$ (see [6]). This means
that the counterflow of the condensate should be taken into account. Such a
counterflow is caused by the condensate’s phase jump across a dark soliton [7,
8] and it leads to the difference between the canonical momentum $p$ of the
soliton quasi-particle and its ‘naive’ mechanical momentum [9, 10]. As a
result, the canonical momentum $p$ of a dark soliton, its energy $\varepsilon$
and its velocity $V$ are related by the classical formula

$V=\frac{\partial\varepsilon}{\partial p}$ (1)

of Hamiltonian mechanics, and this explicitly demonstrates that a dark soliton
behaves like an effective particle (quasi-particle). In particular,
conservation of the energy $\varepsilon$ in a stationary field leads
immediately [11] to the frequency $\omega_{0}/\sqrt{2}$ of dark solitons
oscillations in a harmonic potential. Similar approach yields easily the
theory of evolution of ring solitons in cylindrically symmetrical traps [12]
and this theory agrees very well with numerical simulations.

The aim of this paper is to extend the above approach to situations with a
non-uniform and stationary or time-dependent flow of BEC. The Hamilton
equations are reduced to the Newton-like equation and in this way we reproduce
the results obtained in Ref. [6] by the perturbation method as well as the
results derived in Ref. [13] from the soliton limit of the Whitham modulation
equations. In a stationary case, the Hamiltonian approach provides the non-
trivial energy conservation law.

We illustrate our approach by solutions of the dark soliton motion equation
for such typical situations as its motion along a rarefaction wave [13], along
a self-similar expanding cloud of BEC [14, 15], along a stationary 1D flow of
BEC past an obstacle in the subcritical regime [16, 17]. Solutions of the
Newton equation agree very well with numerical solutions of the Gross-
Pitaevskii (GP) equation.

##  II Hamiltonian approach to dark solitons motion

One-dimensional dynamics of BEC is described with high accuracy by the GP
equation which can be written in the non-dimensional form as

$i\psi_{t}+\frac{1}{2}\psi_{xx}-|\psi|^{2}\psi=U(x)\psi,$ (2)

where $U(x)$ is the external potential. Transition from the BEC wave function
$\psi$ to the more convenient hydrodynamic variables, namely, the condensate’s
density $\rho$ and its flow velocity $u$, is performed by means of the
substitution

$\psi(x,t)=\sqrt{\rho(x,t)}\exp\left({i}\int^{x}u(x^{\prime},t)dx^{\prime}\right),$
(3)

so that the GP equation is cast to the system

$\begin{split}&\rho_{t}+(\rho u)_{x}=0,\\\
&u_{t}+uu_{x}+\rho_{x}+\left[\frac{\rho_{x}^{2}}{8\rho^{2}}-\frac{\rho_{xx}}{4\rho}\right]_{x}=-U_{x}.\end{split}$
(4)

We assume that the external potential $U(x)$ forms a non-uniform distribution
of the background density which changes at the characteristic distances of the
order of magnitude much greater then the healing length of BEC. Since the
typical soliton’s width has the order of magnitude of the healing length, in
derivation of the soliton solution with account of the flow velocity we can
neglect the potential term $U(x)$ in the above equations. It is well known
(see, e.g., [18]) that the traveling wave solution of Eqs. (4) without the
potential term can be written in the form ($\rho=\rho(\xi)$, $u=u(\xi)$,
$\xi=x-Vt$)

$\begin{split}\rho&=\rho_{1}+(\rho_{2}-\rho_{1})\mathrm{sn}^{2}\left(\sqrt{\rho_{3}-\rho_{1}}(x-Vt),m\right),\\\
u&=V\pm\frac{\sqrt{\rho_{1}\rho_{2}\rho_{3}}}{\rho},\quad\rho_{1}\leq\rho_{2}\leq\rho_{3},\end{split}$
(5)

where

$m=\frac{\rho_{2}-\rho_{1}}{\rho_{3}-\rho_{1}}$ (6)

and $\rho(\xi)$ oscillates in the interval

$\quad\rho_{1}\leq\rho\leq\rho_{2}.$ (7)

Two signs in the formula for $u(\xi)$ correspond to two directions of the
velocity of the background flow along which the wave propagates. We have the
same density profiles for both directions of the flow.

We are interested in the soliton solution, so we turn to the soliton limit
$\rho_{3}\to\rho_{2}$ ($m\to 1$) and obtain

$\begin{split}\rho&=\rho_{2}-\frac{\rho_{2}-\rho_{1}}{\cosh^{2}\left[\sqrt{\rho_{2}-\rho_{1}}(x-Vt)\right]},\\\
u&=V\pm\frac{\rho_{2}\sqrt{\rho_{1}}}{\rho}.\end{split}$ (8)

Far enough from the soliton’s center the distributions tend to their
asymptotic background values $\rho\to\rho_{0}=\rho_{2}$, $u\to u_{0}$, that is
$\rho_{1}=(V-u_{0})^{2}$, so we get the soliton solution in terms of the
physical variables

$\begin{split}\rho&=\rho_{0}-\frac{\rho_{0}-(V-u_{0})^{2}}{\cosh^{2}\left[\sqrt{\rho_{0}-(V-u_{0})^{2}}(x-Vt)\right]},\\\
u&=V-\frac{\rho_{0}(V-u_{0})}{\rho}.\end{split}$ (9)

Consequently, the absolute value of the soliton’s speed relative to the BEC
flow is smaller then the sound velocity $c_{0}=\sqrt{\rho_{0}}$ in a uniform
condensate with the density $\rho_{0}$:

$|u_{0}-c_{0}|\leq V\leq u_{0}+c_{0},$ (10)

where $V$ is defined in the ‘laboratory’ reference system. As is clear from
Eq. (9), the amplitude of the dark soliton is related with its velocity by the
formula

$a=\rho_{0}-(V-u_{0})^{2}.$ (11)

We consider a dark soliton as a quasi-particle [11], that is a localized
excitation propagating through the moving condensate. To describe dynamics of
such a soliton, we need to find expressions for its energy and canonical
momentum. As is known [7], the excitation of a dark soliton in BEC is
accompanied by generation of the counterflow in the background condensate:
this counterflow compensates the jump of the phase in BEC $\Delta\phi=\int
udx$ corresponding to the solution (9). Consequently, the canonical momentum
of a dark soliton must include the term resulting from the phase jump
contribution. In case of the background at rest the expression for the
canonical momentum reads

$p=-2V\sqrt{\rho_{0}-V^{2}}+2\rho_{0}\arccos\frac{V}{\sqrt{\rho_{0}}}.$ (12)

If the background moves with velocity $u_{0}$, then we have to subtract
$u_{0}$ from the soliton’s velocity $V$ to obtain the expression

$\begin{split}p=&-2(V-u_{0})\sqrt{\rho_{0}-(V-u_{0})^{2}}\\\
&+2\rho_{0}\arccos\frac{V-u_{0}}{\sqrt{\rho_{0}}}\end{split}$ (13)

for the canonical momentum of dark solitons in the moving background. Then
from Eq. (1) we have

$\varepsilon=\int Vdp=\int V\frac{dp}{dV}dV$

and substitution of (13) followed by simple integration yields

$\varepsilon=\frac{4}{3}\left[\rho_{0}-(V-u_{0})^{2}\right]^{3/2}+u_{0}p\equiv\varepsilon^{(0)}+u_{0}p,$
(14)

where $\varepsilon^{(0)}(V)=\frac{4}{3}(\rho_{0}-V^{2})^{3/2}$ is the well-
known expression for the dark soliton’s energy in the quiescent BEC. Formula
(14) corresponds to the Galileo transformation of the quasi-particle’s energy
in agreement with the Landau approach to the superfluidity theory [8]. The
last term in Eq. (14) can also be considered as the Doppler shift of the
frequency due to motion of the medium.

So far we have considered dark solitons moving along a uniform condensate with
the density $\rho_{0}$ and the flow velocity $u_{0}$. Now we assume that the
distributions of the density $\rho=\rho(x,t)$ and the flow velocity $u=u(x,t)$
are slow functions of the space coordinate $x$ and time $t$, that is they
change little along the distance of the order of magnitude of the soliton’s
width. Then we can introduce the coordinate $x=x(t)$ of the soliton and in the
main approximation $\rho(x(t),t)$, $u(x(t),t)$ represent the density and the
flow velocity of the background distributions at the soliton’s location. In
this approximation, if we make the replacements

$\rho_{0}\mapsto\rho(x,t),\quad u_{0}\mapsto u(x,t),\quad
V\mapsto\dot{x}(t)=\frac{dx}{dt}$

in Eqs. (13) and (14), then we arrive at the expressions for the canonical
momentum

$p=-2(\dot{x}-u)\sqrt{\rho-(\dot{x}-u)^{2}}+2\rho\arccos\frac{\dot{x}-u}{\sqrt{\rho}}$
(15)

and the energy

$\begin{split}\varepsilon=&\frac{4}{3}\left[\rho-(\dot{x}-u)^{2}\right]^{3/2}-2u(\dot{x}-u)\sqrt{\rho-(\dot{x}-u)^{2}}\\\
&+2u\rho\arccos\frac{\dot{x}-u}{\sqrt{\rho}}\end{split}$ (16)

of the particle-like dark soliton located at the moment of time $t$ at the
point $x$, where $\rho=\rho(x,t)$, $u=u(x,t)$. Hamilton equations for the dark
soliton motion have the standard form

$\frac{dx}{dt}=\left(\frac{\partial\varepsilon}{\partial
p}\right)_{x},\qquad\frac{dp}{dt}=-\left(\frac{\partial\varepsilon}{\partial
x}\right)_{p},$ (17)

where it is implied that the velocity $\dot{x}$ is excluded from the energy
(16) with the use of Eq. (15).

Since it is impossible to express the energy (16) as a function of $p$ and $x$
in an explicit form, it is convenient to transform the Hamilton equations (17)
to the Newton-like equation for the soliton’s path $x(t)$. To this end, we
differentiate (15) with respect to time $t$ and find the expression for the
left-hand side of the second equation (17),

$\begin{split}\frac{dp}{dt}=&-4(\ddot{x}-u_{x}\dot{x}-u_{t})\sqrt{\rho-(\dot{x}-u)^{2}}\\\
&+2(\rho_{x}\dot{x}+\rho_{t})\arccos\frac{\dot{x}-u}{\sqrt{\rho}}.\end{split}$
(18)

In Hamiltonian mechanics the velocity $\dot{x}=V$ is considered as a function
$V=V(p,x)$ of the momentum $p$ and coordinate $x$ defined implicitly by Eq.
(15). Therefore the derivative of the right-hand side of the second equation
(17) can be written as

$\begin{split}\left(\frac{\partial\varepsilon}{\partial
x}\right)_{p}=&\left.\frac{\partial\varepsilon^{(0)}}{\partial\rho}\right|_{V,u}\rho_{x}+\left.\frac{\partial\varepsilon^{(0)}}{\partial
u}\right|_{V,\rho}u_{x}\\\ &+\left.\frac{\partial\varepsilon^{(0)}}{\partial
V}\right|_{\rho,u}\cdot\left.\frac{\partial V}{\partial
x}\right|_{p}+u_{x}p,\end{split}$ (19)

where $\varepsilon^{(0)}=\frac{4}{3}[\rho-(V-u)^{2}]^{3/2}$. The derivatives
of $\varepsilon^{(0)}$ with respect to $\rho,u,V$ are calculated without any
difficulty and the derivative $\left.\frac{\partial V}{\partial x}\right|_{p}$
is to be calculated by differentiation of Eq. (15) with $\dot{x}=V(p,x)$ under
the condition that $p=\mathrm{const}$. After simple calculation we obtain

$\left(\frac{\partial V}{\partial
x}\right)_{p}=u_{x}+\frac{\rho_{x}}{2\sqrt{\rho-(V-u)^{2}}}\arccos\frac{V-u}{\sqrt{\rho}}.$
(20)

Substitution of all derivatives into Eq. (19) followed by substitution of Eqs.
(18) and (19) into the second equation (17) yields the equation

$\begin{split}&2\ddot{x}=\rho_{x}+(u+\dot{x})u_{x}+2u_{t}\\\
&+\frac{\rho_{t}+(\rho
u)_{x}}{\sqrt{\rho-(V-u)^{2}}}\arccos\frac{V-u}{\sqrt{\rho}}.\end{split}$ (21)

This equation describes the dynamics of a soliton with account the source
(pumping) and absorbtion of the condensate. Since in our case we assume that a
dark soliton propagates along a smooth background whose evolution obeys the
dispersionless limit of Eqs. (4),

$\rho_{t}+(\rho u)_{x}=0,\qquad u_{t}+uu_{x}+\rho_{x}=-U_{x},$ (22)

we can cast Eq. (21) to the form

$2\ddot{x}=-2U_{x}-\rho_{x}+(\dot{x}-u)u_{x}$ (23)

or

$2\ddot{x}=-U_{x}+u_{t}+u_{x}\dot{x}=-U_{x}+\dot{u}.$ (24)

Equation (24) was derived in framework of the perturbation theory in Ref. [6].

Let us list several important particular cases.

(i) If there is no external potential and a dark soliton propagates along a
large scale wave obeying the Euler equations (22) with $U=0$, then Eq. (24)
gives at once the integral of motion

$2\dot{x}-u=\mathrm{const},$ (25)

which constant value is determined by the initial velocity $V_{0}$ and the
value of the flow velocity $u(x_{0})$ at the initial point $x_{0}$ of the
soliton’s path. This integral follows also from the soliton limit of the
Whitham equations for the GP equation periodic solutions with slowly changing
parameters (see Ref. [13]).

(ii) If the flow is stationary, that is the distributions $\rho=\rho(x)$,
$u=u(x)$ do not depend on time $t$, then the motion equation

$2\ddot{x}=-U_{x}+u_{x}\dot{x}$ (26)

has the integral of energy $\varepsilon(\dot{x},x)=\mathrm{const}$, see Eq.
(16), as it follows immediately from the Hamilton equations (17).

(iii) If there is no flow ($u=0$) and a dark soliton propagates along the
stationary Thomas-Fermi distribution $\rho(x)+U(x)=\mathrm{const}$, then its
motion obeys the Newton equation

$2\ddot{x}=-U_{x}.$ (27)

In case of condensate confined in a harmonic trap $U(x)=\omega_{0}^{2}x^{2}/2$
such a soliton oscillates with the frequency $\omega_{0}/\sqrt{2}$ as was
predicted in Ref. [6] and confirmed in the experiments [19, 20].

Now we turn to consideration of several typical examples a dark soliton
propagation along a non-uniform background.

Figure 1:  Dynamics of a dark soliton moving along expanding BEC. At $t=0$,
the condensate is trapped by a harmonic potential $U(x)=\omega_{0}^{2}x^{2}/2$
(dotted curve). Here $\mu=1$, $\omega_{0}=0.1$ and $V_{0}=0.5$.  Figure 2:
(a) Dark soliton’s path $x(t)$ during its motion in the expanding condensate.
(b) Velocity of the dark soliton $V=\dot{x}(t)$ as a function of time. The red
curves show the results of numerical calculations and the dash-dotted blue
curves show the analytical solution. Dotted lines correspond to $V=V_{0}$ and
dashed lines correspond to $V=2V_{0}$. Here $\mu=1$, $\omega_{0}=0.1$ and
$V_{0}=0.5$.

##  III Motion of dark solitons along a rarefaction wave

One of the simplest non-uniform time-dependent solutions of equations (22) for
the background flow with $U=0$ is a rarefaction wave. Although the motion of
dark solitons along such a wave has already been considered in Ref. [13] by
the Whitham method, we shall consider it briefly for completeness by our
method with some additions and simplifications.

First, we introduce the Riemann invariants

$r_{+}=\frac{u}{2}+\sqrt{\rho},\qquad r_{-}=\frac{u}{2}-\sqrt{\rho},$ (28)

and transform Eqs. (22) with $U=0$ to the diagonal Riemann form

$\frac{\partial r_{-}}{\partial t}+v_{-}\frac{\partial r_{-}}{\partial
x}=0,\qquad\frac{\partial r_{+}}{\partial t}+v_{+}\frac{\partial
r_{+}}{\partial x}=0,$ (29)

where

$v_{-}=\frac{1}{2}(3r_{-}+r_{+}),\qquad v_{+}=\frac{1}{2}(r_{-}+3r_{+}).$ (30)

The physical variables are related with the Riemann invariants by the formulas

$u=r_{-}+r_{+},\qquad\rho=\frac{1}{4}(t_{+}-r_{-})^{2}.$ (31)

A rarefaction wave is a particular self-similar solution of Eqs. (29) with one
of the Riemann invariants constant. To be definite, we assume that
$r_{-}=\mathrm{const}$ and the self-similar rarefaction wave evolves from the
initial discontinuity, located at $x=0$, of the Riemann invariant
$r_{+}=u/2+\sqrt{\rho}$. Let the density equal to $\rho_{-}$ and the flow
velocity equal to zero, $u_{-}=0$ on the left side of the initial
discontinuity. We denote the density on the right side of the discontinuity as
$\rho_{+}$, $\rho_{+}>\rho_{-}$, and then from the condition that $r_{-}$ is
constant across the discontinuity, $r_{-}=u/2-\sqrt{\rho}=-\sqrt{\rho_{-}}$,
we find the flow velocity on the right from the discontinuity is equal to

$u_{+}=2(\sqrt{\rho_{+}}-\sqrt{\rho_{-}}).$ (32)

The first equation (29) is already satisfied and the self-similar solution
$r_{+}=r_{+}(x/t)$ of the second equation gives $v_{+}=(3r_{+}+r_{-})/2=x/t$,
or, with account of $r_{-}=-\sqrt{\rho_{-}}$, we get
$r_{+}=u+\sqrt{\rho_{-}}$. As a result we obtain the distribution of the flow
velocity in the rarefaction wave

$u(x,t)=\left\\{\begin{array}[]{ll}0,&x<c_{-}t,\\\
\frac{2}{3}\left(\frac{x}{t}-c_{-}\right),&c_{-}t<x<(3c_{+}-2c_{-})t,\\\
2(c_{+}-c_{-}),&x>(3c_{+}-2c_{-})t,\end{array}\right.$ (33)

where $c_{\pm}=\sqrt{\rho_{\pm}}$.

Let at the initial moment of time $t=0$ a dark soliton be located on the right
from the discontinuity at the point $x_{+}>0$, let it move with the velocity
$V_{+}$, so that its amplitude equals to (see Eq. (9))

$a_{+}=\rho_{+}-(V_{+}-u_{+})^{2}.$ (34)

Since the right edge of the rarefaction wave propagates with velocity
$3c_{+}-2c_{-}$, it catches up the soliton at the moment

$t_{+}=\frac{x_{+}}{3c_{+}-2c_{-}-V_{+}}$ (35)

provided the denominator of this formula is positive what will be assumed in
what follows.

For $t>t_{+}$ the soliton moves along the rarefaction wave and we can use the
conservation law (25), which in our case becomes the linear differential
equation

$\frac{dx}{dt}-\frac{x}{3t}=V_{+}-c_{+}+\frac{2}{3}c_{-}.$ (36)

It can be readily solved with the initial condition
$x(t_{+})=x_{+}+V_{+}t_{+}$, so we get the formula for the soliton’s path,

$x(t)=\left(V_{+}-\frac{x_{+}}{2t_{+}}\right)t+\frac{3x_{+}}{2}\left(\frac{t}{t_{+}}\right)^{1/3},\qquad
t>t_{+}.$ (37)

The soliton reaches the left edge of the rarefaction wave at the moment
$t_{-}$ defined by the condition $x(t_{-})=c_{-}t_{-}$ which gives

$t_{-}=\frac{1}{\sqrt{t_{+}}}\left(\frac{x_{+}}{c_{+}-V_{+}}\right)^{3/2}.$
(38)

Consequently, the soliton passes through the rarefaction wave if only
$c_{+}>V_{+}$ and in this case Eq. (25) gives at once that for $t>t_{-}$ the
soliton’s velocity becomes equal to $V_{-}=c_{-}-(c_{+}-V_{+})$. Its amplitude
$a_{-}=\rho_{-}-V_{-}^{2}$ can be transformed with the help of Eq. (34) to

$a_{-}=a_{+}-2(c_{+}-c_{-})\left(c_{+}\pm\sqrt{c_{+}^{2}-a_{+}}\right),$ (39)

where the choice of the sign is determined by the sign of the relative
velocity $V_{+}-u_{+}$.

If $c_{+}=V_{+}$, then the soliton’s path is given by

$x(t)=c_{-}t+\frac{3x_{+}}{2}\left(\frac{t}{t_{+}}\right)^{1/3},\quad
t>t_{+}=\frac{x_{+}}{2(c_{+}-c_{-})},$ (40)

that is it propagates ahead of the rear (left) edge of the rarefaction wave.

If $c_{+}<V_{+}$, the soliton also remains forever inside the rarefaction
wave, but this time its velocity tends according to Eq. (36) to the constant
value

$V_{\infty}=c_{-}+\frac{3}{2}(V_{+}-c_{+}),$ (41)

which is greater that the velocity $c_{-}$ of the left edge of the rarefaction
wave. Along the limiting trajectory $x=V_{\infty}t$ the background flow has
the values $u_{\infty}=V_{+}-c_{+}$,
$\rho_{\infty}=[c_{-}+(V_{+}-c_{+})/2]^{2}$. Hence, the relative soliton’s
velocity $V_{\infty}-u_{\infty}$ is equal to the local sound velocity
$c_{\infty}=\sqrt{\rho_{\infty}}$, that is the soliton’s amplitude tends to
zero in the limit $t\to\infty$.

##  IV Dark soliton’s motion along expanding condensate

Now we assume that the initial background state is formed by a harmonic
potential $U(x)=\omega_{0}^{2}x^{2}/2$ (in non-dimensional variables $m=1$)
and the distribution of density is given by the Thomas-Fermi formula
$\rho(x,0)=(2\mu-\omega_{0}^{2}x^{2})/2,|x|\leq\sqrt{2\mu}/\omega_{0}$, where
$\mu$ is the chemical potential of BEC, and the initial flow velocity is equal
to zero everywhere. At the initial moment of time $t=0$ the trap is switched
off and the condensate starts its expansion. The dispersionless equations (22)
can be solved exactly in this case, the solution was found in Ref. [14] for
the focusing NLS equation and it was adapted to the present repulsive BEC
situation in Ref. [15]. It is given by the formulas

$\begin{split}\rho(x,t)&=\frac{1}{f(t)}\left(\mu-\frac{x^{2}}{\phi_{0}^{2}f^{2}(t)}\right),\qquad
u(x,t)=xg(t)\\\
g(t)&=\frac{f^{\prime}(t)}{f(t)}=\frac{2}{\phi_{0}}\sqrt{\frac{f(t)-1}{f^{3}(t)}},\qquad\phi_{0}=\frac{\sqrt{2}}{\omega_{0}},\end{split}$
(42)

where the function $f(t)$ satisfies the differential equation

$\frac{df}{dt}=\frac{2}{\phi_{0}}\sqrt{\frac{f-1}{f}},\qquad f(0)=1$ (43)

and it can be defined in implicit form by the equation

$t(f)=\frac{\phi_{0}}{2}\left[\sqrt{f(f-1)}+\ln(\sqrt{f-1}+\sqrt{f})\right].$
(44)

Let the dark soliton be formed at the initial moment of time $t=0$ at the
point $x=0$ with the velocity $V_{0}>0$. Then its path can be found by solving
the equation (25) which in our case reads

$\frac{dx}{dt}-\frac{1}{2}g(t)x=V_{0},$ (45)

with the initial condition $x(0)=0$. This linear differential equation can be
easily solved and with the use of the expressions for $g(t)$ and $df/dt$ the
solution can be cast to the form

$x(f)=V_{0}\phi_{0}\sqrt{f(f-1)}.$ (46)

This formula together with Eq. (44) define the soliton’s path $x=x(t)$ in
parametric form. The soliton’s velocity is to be found from Eq. (45),

$\dot{x}(f)=V_{0}\left(2-\frac{1}{f}\right).$ (47)

Typical condensate profiles at different moments of time $t$ are shown in Fig.
1 for $\mu=1$, $\omega_{0}=0.1$ and $V_{0}=0.5$. Plots of $x(t)$ and
$V=\dot{x}(t)$ are shown in Fig. 2. In the limit $t\to\infty$, when
$f\to\infty$, the velocity tends to its limiting value
$\dot{x}(\infty)=2V_{0}$ equal to the value of the flow velocity
$u=x/t=2V_{0}$ along the limiting path $x(t)=2V_{0}t$. Consequently, at
asymptotically large times the soliton becomes black and it is transferred by
the flow.

##  V Dark soliton’s motion along a hydraulic flow of BEC through a penetrable
barrier

Figure 3: Soliton’s paths for fixed values of the parameters $\rho_{0}$,
$u_{0}$ of the flow at infinity $x\to\pm\infty$, the fixed potential barrier
shape and different values of the soliton’s initial velocity $V_{0}$. At some
critical value of the velocity $V_{0}$ the paths passing through the barrier
change to the paths reflected from it. The plot corresponds to the values
$\rho_{0}=1.0$, $u_{0}=0.1$, $U_{m}=0.4$, and $\sigma=300$ (see Eq. (55)).
Figure 4: Downstream propagation of a dark soliton through the non-uniform
flowing condensate. Red curve shows the density $\rho(x)$ and the dotted line
is the potential $U(x)$. The background density at infinity is taken equal to
$\rho_{0}=1$ and flow velocity is $u_{0}=0.1$, the initial soliton’s velocity
is $V_{s}=0.8$. The potential has a shape (55) with the maximal value
$U_{m}=0.4$ and the width $\sigma=300$. The arrows show direction of the
background flow. Figure 5: The dependence of the soliton position on time for
the dynamics shown in Fig. 4. The dash-dotted blue line corresponds to the
solution of Eq. (26) and the red line to the numerical solution of the GP
equation (2). The dotted line shows the soliton’s path through a uniform
condensate when $U_{m}=0$.

Now we shall consider a stationary flow of BEC past a wide barrier whose
action is described by the potential $U(x)$. The large scale dependence of the
density $\rho$ and the flow velocity $u$ is determined by the stationary
version of Eqs. (22),

$(\rho u)_{x}=0,\quad uu_{x}+\rho_{x}+U_{x}=0,$ (48)

which are readily integrated to give

$\rho
u=\rho_{0}u_{0},\quad\frac{1}{2}u^{2}+\rho+U=\frac{1}{2}u_{0}^{2}+\rho_{0},$
(49)

where $\rho_{0}$ and $u_{0}$ are the density and the flow velocity of the
condensate at $x\to\pm\infty$. We imply here that the potential is localized
around the point $x=0$ and its distance of action $\sigma$ is much greater
than the healing length which is equal to unity in our dimensionless notation.
Elimination of $\rho=\rho_{0}u_{0}/u$ yields the equation

$\frac{1}{2}\left[u_{0}^{2}-u^{2}(x)\right]+\rho_{0}\left(1-\frac{u_{0}}{u(x)}\right)=U(x),$
(50)

which defines in implicit form the function $u=u(x)$ for the space dependence
of the background flow velocity. When it is found, the distribution of the
background density is given by the formula

$\rho(x)=\frac{\rho_{0}u_{0}}{u(x)}.$ (51)

Equation (50) is cubic and its solution should be chosen in such a way that it
satisfies the boundary condition $u(x)\to u_{0}$ at both infinities
$x\to\pm\infty$. As was noticed in Refs. [16, 17], this imposes a very
important condition on possible values of $u_{0}$. Indeed, at the point where
the potential $U(x)$ reaches its maximal value
$U_{m}=\max\\{U(x):-\infty<x<\infty\\}$, the function in the left-hand side of
Eq. (50) has the maximum equal to
$\frac{1}{2}u_{0}^{2}-\frac{3}{2}(\rho_{0}u_{0})^{2/3}+\rho_{0}$ at
$u(x)=u_{m}=(u_{0}\rho_{0})^{1/3}$. Consequently, the solution $u(x)$ for all
values of $x$, $-\infty<x<\infty$, exists if only

$U_{m}\leq\frac{1}{2}u_{0}^{2}-\frac{3}{2}\left(\rho_{0}u_{0}\right)^{2/3}+\rho_{0}.$
(52)

This inequality becomes equality for $u_{0}$ equal to two roots $u_{\pm}$ of
the equation

$\frac{1}{2}u_{0}^{2}-\frac{3}{2}(\rho_{0}u_{0})^{2/3}+\rho_{0}=U_{m}.$ (53)

Then it is easy to find that the smooth solution $u=u(x)$ exists for all $x$
in the subcritical $u_{0}<u_{-}$ and supercritical $u_{0}>u_{+}$ regimes of
the flow past an obstacle described by the potential $U(x)$. We assume that
$u_{0}>0$, so the root $u_{-}$ should be put equal to zero for
$U_{m}\geq\rho_{0}$. If $U_{m}\ll\rho_{0}$, then $u_{0}$ must be close to the
asymptotic value of the sound velocity $\sqrt{\rho_{0}}$, so one can easily
find that up to the second order in the small parameter
$(U_{m}/\rho_{0})^{1/2}$ we have

$u_{\pm}\approx\sqrt{\rho_{0}}\left(1\pm\sqrt{\frac{3U_{m}}{2\rho_{0}}}+\frac{U_{m}}{12\rho_{0}}\right).$
(54)

In the transcritical regime $u_{-}<u_{0}<u_{+}$, the condition (52) does not
hold. In this case, the flow past a wide barrier leads to generation of
dispersive shock waves in this finite interval of the flow velocities [16,
17]. We are only interested in the situations when dispersive shock waves are
not formed and the soliton passes through a stationary fairly smooth profile
of the condensate.

In our concrete examples we shall use the potential

$U(x)=\frac{U_{m}}{\cosh(x/\sigma)},$ (55)

where $\sigma=300$ and $U_{m}=0.4$. We take the initial coordinate of the dark
soliton $x_{0}$ upstream the flow far away from the barrier and launch the
dark soliton with the initial velocity $V_{0}$. Soliton’s motion obeys the
Newton equation (26) which can be easily solved numerically with the initial
conditions $x(0)=x_{0}$, $\dot{x}(0)=V_{0}$.

Figure 6: Reflection of a dark soliton from the potential barrier. Red curve
shows the density $\rho(x)$, and the dotted line depicts the potential $U(x)$.
Background density is equal to $\rho_{0}=1$ and flow velocity to $u_{0}=0.1$,
initial soliton velocity is $V_{s}=0.6$. The potential has the form (55) with
$U_{m}=0.4$ and $\sigma=300$. Figure 7: The path of the soliton for its
propagation for the dynamics shown in Fig. 6. Figure 8: Soliton amplitude for
the dynamics shown in Fig. 6. The green dotted lines correspond to the
asymptotic amplitudes of the soliton, which are in agreement with Eq. (63).

First of all, we notice that there exist two types of solitons trajectories:
(i) trajectories which pass through the potential region from one infinity to
the other one, and (ii) trajectories with reflection of solitons from the
potential region. This is illustrated in Fig. 3. We fix the parameters
$\rho_{0}$, $u_{0}$ of the flow at infinity $x\to\pm\infty$ and calculate
solitons paths for different values of the soliton’s initial velocity $V_{0}$
in the range from $0.5$ to $0.8$. At some critical value of the soliton’s
velocity $V_{0}\approx 0.61$ the paths switch from those passing through the
barrier to the paths reflected from it. At the critical value $a_{0}=a_{cr}$,
which separates these two regimes, the trajectory reaches the point $x=0$ of
the potential maximum with zero velocity, that is $\dot{x}(x=0)=0$ at the
point with $u(0)=u_{m},\rho(0)=\rho_{m}$. Substitution of these values of into
Eq. (16) gives the soliton’s critical energy

$\begin{split}\varepsilon_{cr}=&\frac{2}{3\rho_{m}^{3}}(2\rho_{m}^{3}+\rho_{0}^{2}u_{0}^{2})\sqrt{\rho_{m}^{3}-\rho_{0}^{2}u_{0}^{2}}\\\
&+2u_{0}\rho_{0}\left(\frac{\pi}{2}+\arcsin\frac{\rho_{0}u_{0}}{\rho_{m}^{3/2}}\right).\end{split}$
(56)

Equating this soliton’s energy at the turning point to its initial energy, we
arrive at the equation

$\begin{split}\frac{4}{3}a_{cr}^{3/2}-2u_{0}\sqrt{a_{cr}(\rho_{0}-a_{cr})}+2u_{0}\rho_{0}\arcsin\sqrt{\frac{a_{cr}}{\rho_{0}}}=\varepsilon_{cr}.\end{split}$
(57)

If we define the variables

$A=\frac{a_{cr}}{\rho_{0}},\qquad M=\frac{u_{0}}{\sqrt{\rho_{0}}},$ (58)

then we can write this equation in the form

$M(A)=\frac{4A^{3/2}-3\varepsilon_{cr}/\rho_{0}^{3/2}}{6\left(\sqrt{A(1-A)}-\arcsin\sqrt{A}\right)}.$
(59)

It defines in implicit form the dependence of the non-dimensional critical
amplitude on the Mach number of the flow far from the obstacle. In particular,
for the values $\rho_{0}=1.0$, $u_{0}=0.1$, $U_{m}=0.4$, we obtain
$\rho_{min}=0.591$ and Eq. (59) yields $A_{cr}=0.738$ ($V_{0}=0.612$) in
reasonable agreement with the value obtained in numerical solutions of the GP
equation.

If we turn to the trajectories with reflection of solitons from the obstacle,
then we see that the initial and final amplitudes of such a soliton are
different due to change of sign of the velocity:

$a_{+}=\rho_{0}-(V_{+}-u_{0})^{2},\quad a_{-}=\rho_{0}-(V_{-}-u_{0})^{2},$
(60)

where both velocities are defined far enough from the obstacle. Let the
soliton start its downstream motion at $x=-\infty$ with the velocity
$V_{+}>0$, so its initial energy is given by the expression

$\varepsilon=\frac{4}{3}a_{+}^{3/2}-2u_{0}\sqrt{a_{+}(\rho_{0}-a_{+})}+2u_{0}\rho_{0}\arcsin\sqrt{\frac{a_{+}}{\rho_{0}}},$
(61)

and it must be equal to its energy after reflection when it moves upstream
with the velocity $V_{-}<0$ so that its energy is given by

$\begin{split}\varepsilon=\frac{4}{3}a_{-}^{3/2}+&2u_{0}\sqrt{a_{-}(\rho_{0}-a_{-})}\\\
+&2u_{0}\rho_{0}\left(\pi-\arcsin\sqrt{\frac{a_{-}}{\rho_{0}}}\right).\end{split}$
(62)

Equating these two expressions and introducing again the variables
$A_{+}=a_{+}/\rho_{0}$, $A_{-}=a_{-}/\rho_{0}$, $M=u_{0}/\sqrt{\rho_{0}}$, we
get the equation

$\begin{split}\frac{2}{3}A_{-}^{3/2}&+M\sqrt{A_{-}(1-A_{-})}+M\left(\pi-\arcsin\sqrt{A_{-}}\right)\\\
&=\frac{2}{3}A_{+}^{3/2}-M\sqrt{A_{+}(1-A_{+})}+M\arcsin\sqrt{A_{+}},\end{split}$
(63)

which determines the final amplitude of the soliton in terms of its initial
amplitude and the Mach number of the flow.

We compared our analytical findings with numerical solutions of the GP
equation. An example of soliton’s propagation for a subcritical
($u_{0}<u_{-}$) background flow velocity is shown in Fig. 4. The soliton moves
downstream the condensate flow. Fig. 5 shows the soliton’s trajectory
calculated according to Eq. (26) and by means of numerical solution of the GP
equation (2). As one can see, the approximate analytical theory (dash-dotted
blue line) agrees perfectly well with the exact numerical solution (red line).
Similar agreement between two approaches was found also for supercritical flow
of the background BEC and for the upstream initial motion of solitons.

Another case, when the soliton is reflected from the barrier, is illustrated
in Fig. 6. The solution of Eq. (26) agrees perfectly well with the exact
numerical solution of the GP equation, as is shown in Fig. 7. One can clearly
see that the soliton’s amplitude changes after reflection. Fig. 8 shows the
dependence of the amplitude on $t$. As one can see, the numerical amplitude of
the soliton after the collision coincides with the analytical one obtained
from Eq. (63).

##  VI Conclusion

The motion of dark solitons is an interesting problem of nonlinear physics
because such a motion cannot be completely separated from the motion of the
background medium. Most spectacularly this phenomenon is demonstrated by the
difference in the motion of the center of mass of BEC confined in a trap and
the motion of the soliton itself described by its mean coordinate.
Consequently, the soliton’s motion is accompanied by some counterflow in the
background condensate which leads to a proper redistribution of the density.
In 1D geometry such a counterflow is necessary for topological reasons: the
flow velocity of the condensate is a gradient of the phase of the condensate’s
wave function, so formation of a dark soliton leads to the jump of phase
across the soliton and, to keep the wave function single-valued, this jump
must be compensated by the flow outside the soliton [7]. We show that the
contribution of the counterflow into the momentum and the energy of a dark
soliton is crucially important for correct description of its motion as a
localized quasi-particle by the methods of Hamiltonian mechanics. In this way
we have reproduced the results obtained earlier by the methods of perturbation
theory [6] and in framework of the Whitham theory [13] (in the last case only
for situations without external potential). The advantage of the Hamiltonian
approach is that it provides immediately the quite non-trivial energy
conservation law for the soliton’s motion in stationary background flows and
it would be difficult to find this law by other methods. We have illustrated
our approach by several examples and confirmed the analytical results by their
comparison with numerical solutions of the GP equation.

We believe that the suggested here Hamiltonian approach can be applied to
other interesting situations, such as, for example, interaction of solitons
with non-convex flows [21] or “magnetic” solitons [22, 23] in BEC, when the
main contribution into the soliton’s canonical momentum is made by the
counterflow.

###### Acknowledgements.

We thank D. V. Shaykin for useful discussions. This work was partially
supported by the Foundation for the Advancement of Theoretical Physics and
Mathematics “BASIS.”

## References

  * [1] N. J. Zabusky, M. D. Kruskal, Phys. Rev. Lett. 15, 240 (1965).
  * [2] V. I. Karpman, E. M. Maslov, Sov. Phys. JETP, 46, 281 (1977).
  * [3] V. I. Karpman, E. M. Maslov, Sov. Phys. JETP, 48, 252 (1978).
  * [4] D. J. Kaup and A. C. Newell, Proc. Roy. Soc. London, A 361, 413 (1978).
  * [5] A. M. Kosevich, Physica D 41, 253 (1990).
  * [6] Th. Busch, J. R. Anglin, Phys. Rev. Lett. 84, 2298-2301 (2000).
  * [7] S. I. Shevchenko, Sov. J. Low Temp. Phys. 14, 553 (1988).
  * [8] L. P. Pitaevskii, S. Stringari, Bose-Einstain Condensation, (Clarendon, Oxfors, 2003).
  * [9] Yu. S. Kivshar, X. Yang, Phys. Rev. E 49, 1657 (1994).
  * [10] Yu. S. Kivshar, W. Królikowski, Optics Commun., 114, 353-362 (1995).
  * [11] V. V. Konotop, L. P. Pitaevskii, Phys. Rev. Lett. 93, 240403 (2004).
  * [12] A. M. Kamchatnov, S. V. Korneev, Phys. Lett. A, 374, 4625 (2010).
  * [13] P. Sprenger, M. A. Hoefer, G. A. El, Phys. Rev. E 97, 032218 (2018).
  * [14] V. I. Talanov, JETP Lett., 2, 138 (1965).
  * [15] V. A. Brazhnyi, A. M. Kamchatnov, V. V. Konotop, Phys. Rev. A 68, 035603 (2003).
  * [16] V. Hakim, Phys. Rev. E 55, 2835 (1997).
  * [17] A. M. Leszczyszyn, G. A. El, Yu. G. Gladush, A. M. Kamchatnov, Phys. Rev. A 79, 063607 (2009).
  * [18] A. M. Kamchatnov, Nonlinear Periodic Waves and Their Modulations—An Introductory Course, (World Scientific, Singapore, 2000).
  * [19] Ch. Becker, S. Stellmer, P. Soltan-Panahi, S. Dörscher, M. Baumert, E.-M. Richter, J. Kronjäger, K. Bongs, K. Sengstock, Nature Phys., 4, 496-501 (2008).
  * [20] A. Weller, J. P. Ronzheimer, C. Gross, J. Esteve, M. K. Oberthaler, D. J. Frantzeskakis, G. Theocharis, P. G. Kevrekidis, Phys. Rev. Lett., 101, 130401 (2008).
  * [21] K. van der Sande, G. A. El, M. A. Hoefer, J. Fluid Mech., 928, A21 (2021).
  * [22] C. Qu, L. P. Pitaevskii, and S. Stringari, Phys. Rev. Lett. 116, 160402 (2016).
  * [23] L. P. Pitaevskii, Phys.-Uspekhi, 59, 1028 (2016).