# Dynamics of quantum double dark-solitons and an exact finite-size scaling of
Bose-Einstein condensation

Kayo Kinjo1, Jun Sato2 and Tetsuo Deguchi3 1 Department of Physics, Center for
Soft Matter Physics, Ochanomizu University, Bunkyo-ku, Tokyo 112-8610, Japan
2 Faculty of Engineering, Tokyo Polytechnic University, Atsugi 243-0297, Japan
3 Department of Physics, Faculty of Core Research, Ochanomizu University,
Bunkyo-ku, Tokyo 112-8610, Japan
<EMAIL_ADDRESS>
###### Abstract

We show several novel aspects in the exact non-equilibrium dynamics of quantum
double dark-soliton states in the Lieb-Liniger model for the one-dimensional
Bose gas with repulsive interactions. We also show an exact finite-size
scaling of the fraction of the quasi-Bose-Einstein condensation (BEC) in the
ground state, which should characterize the quasi-BEC in quantum double dark-
soliton states that we assume to occur in the weak coupling regime. First, we
show the exact time evolution of the density profile in the quantum state
associated with a quantum double dark-soliton by the Bethe ansatz. Secondly,
we derive a kind of macroscopic quantum wave-function effectively by exactly
evaluating the square amplitude and phase profiles of the matrix element of
the field operator between the quantum double dark-soliton states. The
profiles are close to those of dark-solitons particularly in the weak-coupling
regime. Then, the scattering of two notches in the quantum double dark-soliton
state is exactly demonstrated. It is suggested from the above observations
that the quasi-BEC should play a significant role in the dynamics of quantum
double dark-soliton states. If the condensate fraction is close to 1, the
quantum state should be well approximated by the quasi-BEC state where the
mean-field picture is valid.

††: J. Phys. A: Math. Gen.

  * June 2022

Keywords: quantum dynamics, dark-soliton, BEC, finite-size scaling

##  1 Introduction

The experimental realization of trapped atomic gases in one dimension (1D) has
provided a new motivation for the study of strong correlations in fundamental
quantum mechanical systems of interacting particles [1, 2, 3, 4, 5].
Furthermore, the non-equilibrium dynamics of closed interacting quantum
systems is now extensively studied in 1D by experiments and theories [6, 7,
8]. In many 1D quantum interacting systems quantum fluctuations may play a key
role and often lead to subtle nontrivial effects. We thus expect that
fundamental many-body properties such as the quasi-Bose-Einstein condensation
(BEC) should play a key role in the nontrivial quantum dynamics such as
quantum dark-solitons. We shall define it shortly with the Penrose-Onsager
criterion.

Let us introduce a theoretical model for the 1D system of interacting bosons
with repulsive short-range potentials. Here we call it the 1D Bose gas. For
simplicity we assume that the interactions are given by the delta-function
potentials, since they give nontrivial effects in the 1D case although they
are simple. For instance, the scattering length depends on the strength of the
delta-function potential in 1D systems. We thus have the Lieb-Liniger model
(LL model) as the system of the 1D Bose gas. The Hamiltonian of the LL model
is given by [9, 10]

$\displaystyle\mathcal{H}_{\rm{LL}}=-\sum_{j=1}^{N}{\frac{\partial^{2}}{\partial
x_{j}^{2}}}+2c\sum_{j<k}^{N}\delta(x_{j}-x_{k}).$ (1)

Here $N$ denotes the number of bosons, and we assume the periodic boundary
conditions of the system size $L$ on the wave-functions. We employ a system of
units with $2m=\hbar=1$, where $m$ denotes the mass of the particle. We recall
that the coupling constant $c$ is positive. It is an exactly solvable model of
the 1D quantum many-body system. It is known that all the eigenvectors are
constructed by the Bethe-ansatz method [11]. Furthermore, the Gross-Pitaevskii
(GP) equation appears as the Heisenberg equation of motion for the second-
quantized Hamiltonian of the LL model. It is expressed in terms of the
classical complex scalar field $\psi$ as follows [12].

$\displaystyle i\frac{\partial}{\partial t}\psi=-\frac{\partial^{2}}{\partial
x^{2}}\psi+2c|\psi|^{2}\psi-\mu\psi.$ (2)

We expect that the GP equation should play a central role in the long-distance
mean-field behavior of the 1D Bose gas in some quantum state if the quasi-BEC
occurs in the quantum state of the LL model especially in the weak-coupling
regime. If it is the case, the solution of the GP equation should correspond
to the macroscopic wave-function of the quasi-BEC state, and describe the
quantum state well at least approximately.

We define the quasi-BEC by the criterion due to Penrose and Onsager [13, 14]
(see also Section 4.2). Suppose that particle number $N$ is very large but
finite. The density matrix at zero temperature is given by the ground state
$|\lambda\rangle$ of the system as
$\hat{\rho}=|\lambda\rangle\langle\lambda|$. Then, we define the one-particle
reduced density matrix by its partial trace with respect to all but one degree
of freedom: $\hat{\rho}_{1}=N\rm{tr}_{23\cdots N}\hat{\rho}$. Let $N_{0}$
denote the largest eigenvalue of the one-particle reduced density matrix
$\hat{\rho}_{1}$. If it is of order $N$, i.e., the ratio $n_{0}=N_{0}/N$ is
nonzero and finite for large $N$, then we say that the system exhibits the
quasi-BEC, and we call $n_{0}$ the condensate fraction.

If the quasi-BEC occurs in some quantum states of the LL model, we expect that
the GP equation should play a central role for characterizing the quantum
state, although it is only a partial differential equation for a complex
scalar variable. In the present research, we assume that the quasi-BEC should
occur if the coupling constant is small enough with respect to the system size
or the number of bosons, and hence some solutions of the GP equation such as
multiple dark-solitons can be compared with the density profiles of some
quantum states in the quasi BEC of the 1D Bose gas. In fact, we shall show a
finite-size scaling of the quasi BEC in the present research.

It should be emphasized that such quantum states whose density profiles
coincide with those of single dark-solitons of the GP equation have been
constructed explicitly in the form of superposition of the yrast states in the
Lieb-Liniger model [15]. The construction resolved a long standing problem
suggested by Ishikawa and Takayama almost forty years ago [16]. Here we remark
that it was shown through the strong coupling limit [17, 18] that the yrast
states and the mean-field solitons are closely related to each other with
respect to quantum numbers. Furthermore, several significant properties in the
non-equilibrium dynamics of a quantum single dark-soliton have been exactly
investigated [19] and the generic and the ideal Gaussian weights have been
introduced [20, 21]. Moreover, the density and phase profiles of quantum
states of double dark-solitons have been explicitly constructed [22], and the
phase shift has numerically been estimated in the scattering of two quantum
dark-solitons [23].

There is another aspect of quantum dark-soliton states. Successive
measurements of particle positions in the Lieb–Liniger model also leads to
observing quantum dark-solitons numerically [24, 25]. There is a question of
how the density profile of a superposition of yrast states is related to the
successive measurements of particle positions. When the coupling constant $c$
is equal to zero it was analytically shown that the construction of the
quantum dark-soliton state with the Gaussian weight [21] is related to the
particle position method [24] as shown in Ref. [21]. When the coupling
constant is small and nonzero: $c>0$, an ansatz was proposed to bridge between
the calculation of single-particle density and the particle position method
[26].

In the present paper we show various novel aspects in the exact non-
equilibrium dynamics of quantum double dark-solitons, which give pairs of
notches in the density profiles, by explicitly constructing corresponding
quantum states in the Lieb-Liniger model of the 1D Bose gas with the repulsive
interactions. For instance, we exhibit the time evolution of the density
profile of the double dark-soliton whose two notches are located at the same
position, and that of the phase profiles of the quantum double dark-solitons.
In particular, we give an example where the winding number of the phase
profile changes during the scattering process of two notches. Furthermore, we
also show an exact finite-size scaling of the fraction of the BEC for the
ground state. It should characterize the quasi-BEC which we assume to occur in
quantum double dark-soliton states in the weak coupling regime. We show that
if the coupling constant decreases as a power of the system size, condensate
fraction does not vanish and remains constant when we send the system size to
a very large value with fixed density. We recall that if the condensate
fraction is nonzero for a large particle number $N$, we call it the quasi-BEC
by employing the Penrose-Onsager criterion. It follows from it that the quasi-
BEC occurs only if the coupling constant is very small with respect to the
system size. Therefore quantum states of dark-solitons may appear particularly
in the weak coupling regime.

Based on the definition of the quasi-BEC we derive a kind of macroscopic
quantum wave-function by exactly deriving the amplitude and phase profiles of
the matrix element of the bosonic field operator, by making use of Slavnov’s
formula of form factors [27]. Here we recall that the bosonic field operator
is defined in the second-quantized Hamiltonian of the Lieb-Liniger model [28].

Let us briefly summarize the finite-size scaling of the quasi-BEC for the
ground state, which we shall show in detail in Section 4. The scaling behavior
of the quasi-BEC in the 1D Bose gas is fundamental when we send particle
number $N$ or system size $L$ to very large values. We define the interaction
parameter $\gamma$ by $\gamma=c/n$ with coupling constant $c$ in the delta-
function potentials and density $n=N/L$. We show that if $\gamma$ is given by
a negative power of $N$, i.e. $\gamma=A/N^{\eta}$, condensate fraction $n_{0}$
is nonzero and constant for any large value of $L$ or $N$. We also show that
exponent $\eta$ and amplitude $A$ are independent of density $n$, and evaluate
them as functions of $n_{0}$. Thus, the condensate fraction $n_{0}$ for the
ground state is given by a scaling function of variable $\gamma N^{\eta}$,
which corresponds to amplitude $A$. If the condensate fraction of a given
quantum state with large $N$ is nonzero in the 1D Bose gas, we suggest that
the classical mean-field approximation such as the GP equation should be valid
for the state [15]. Furthermore, we show that the 1D Bose gas of a finite
particle number may have the same condensate fraction for any large $L$ in the
case of the ground state.

Finally, we mention some potentially relevant results in the following. For
strong and intermediate interaction strengths, the Lieb-Liniger Gross-
Pitaevski equation is introduced, which is an extension of the GP equation
[29]. Associated with the quantum states of dark solitons, bound states of
dark solitons are numerically studied by solving the GP equation [30],
dynamics of a bright soliton in the quasi-BEC with time-dependent atomic
scattering length in a repulsive parabolic potential [31], quantized quasi-
two-dimensional Bose-Einstein condensates with spatially modulated
nonlinearity [32], matter rogue wave in Bose-Einstein condensates with
attractive atomic interaction [33], exact soliton solutions, and nonlinear
modulation instability in spinor Bose-Einstein condensates [34].

The contents of the paper consist of the following. In Section 2 we explain
the Bethe ansatz and useful formulas for evaluating the form factors of the
field operator. We also define the winding number for solutions of the GP
equation under the periodic boundary conditions. In Section 3 we show the time
evolution of the quantum double dark-soliton state constructed with equal
weight for the following two cases: (i) The soliton positions $X_{1}$ and
$X_{2}$ are different: $X_{1}=L/4$ and $X_{2}=3L/4$; (ii) the soliton
positions are the same: $X_{1}=X_{2}=0$. We also show the time evolution of
the quantum double dark-soliton state constructed with the Gaussian weights.
Here, two notches have different speeds thanks to the Gaussian weights, and we
evaluate the phase shift in the collision of the two dark solitons. We remark
that two notches have mostly the same speed if the quantum double dark-soliton
state is constructed with equal weight. In Section 4 we show the finite-size
scaling behavior of the condensate fraction in the ground state for the 1D
Bose gas with repulsive interactions at zero temperature. According to it, we
can estimate that the fraction of the quasi-BEC condensate should be equal to
0.99 for the quantum double dark-soliton state with $N=L=20$ and $c=0.05$
studied in the present research.

##  2 Method

###  2.1 Bethe ansatz equations

In the LL model, the Bethe ansatz offers an exact eigenstate with an exact
energy eigenvalue for a given set of quasi-momenta $k_{1},k_{2},\ldots,k_{N}$
satisfying the Bethe ansatz equations (BAE) for $j=1,2,\ldots,N$:

$\displaystyle k_{j}L=2\pi I_{j}-2\sum_{\ell\neq
j}^{N}\arctan\left({\frac{k_{j}-k_{\ell}}{c}}\right).$ (3)

Here $I_{j}$’s are integers for odd $N$ and half-odd integers for even $N$. We
call them the Bethe quantum numbers. The total momentum $P$ and the energy
eigenvalue $E$ are expressed in terms of the quasi-momenta as

$\displaystyle P=\sum_{j=1}^{N}k_{j}=\frac{2\pi}{L}\sum_{j=1}^{N}I_{j},\quad
E=\sum_{j=1}^{N}k_{j}^{2}.$ (4)

If we specify a set of Bethe quantum numbers $I_{1}<\cdots<I_{N}$, the BAE in
Equation (3) have a unique real solution $k_{1}<\cdots<k_{N}$ [28, 11]. In
particular, the sequence of the Bethe quantum numbers of the ground state is
given by

$\displaystyle I_{j}=-(N+1)/2+j\qquad\textrm{for $j\in\mathbb{Z}$ with $1\leq
j\leq N$}.$ (5)

The Bethe quantum numbers for low lying excitations are systematically derived
by putting holes or particles in the perfectly regular ground-state sequence.

###  2.2 Coupling constant

In the thermodynamic limit several physical quantities of the LL model are
characterized by the single parameter $\gamma=c/n$, where $n=N/L$ is the
density of particle number $N$. We often fix the particle-number density as
$n=1$ throughout the present paper, and change coupling constant $c$ so that
we have different values of $\gamma$.

###  2.3 Quantum double dark-soliton state

A quantum state that has two notches in both profiles of density and square
amplitude of the matrix element of the field operator was proposed in [22]. We
call it the quantum double dark-soliton state, and it is given by the
superposition of “two-hole” excitation states as follows.

$\displaystyle|X_{1},X_{2},N\rangle:=\frac{1}{\sqrt{\mathcal{M}_{N}}}\sum_{\bm{p}\in\bm{P}_{N}}e^{i(p_{1}X_{1}+p_{2}X_{2})}|p_{1},p_{2},N\rangle$
(6)

with a normalization factor $\mathcal{M_{N}}$ for $N$ particles. The quantum
state $|p_{1},p_{2},N\rangle$ is characterized by a configuration of Bethe
quantum numbers that has two vacancies located at $p_{1}$ and $p_{2}$ in the
series of the Bethe quantum numbers, which is illustrated in Figure 1 (a).
This configuration represents the Bethe quantum numbers of the ground state of
$N$ particles along with those of additional two particles. In Equation (6),
the two holes are denoted as $\bm{p}:=\\{p_{1},p_{2}\\}$ hereinafter, and the
set of all allowed pairs of two holes is represented by $\bm{P}_{N}$. For
example, when the number of particles $N=5$,
$\bm{P}_{N}=\left\\{\\{-2,-1\\},\\{-2,0\\},\\{-2,1\\},\\{-2,2\\},\cdots,\\{2,4\\},\\{3,4\\}\right\\}$,
which has $|\bm{P}_{N}|=21$ elements. In the panel (b) of Figure 1 some
configurations with two holes $p_{1}$ and $p_{2}$ are exhibited. In the third
configuration, two holes $p_{1}$ and $p_{2}$ are located in its middle part of
the series which corresponds to the ground state of $N$ particles. Here we
remark that in order for two notches have positive velocities we derive two
hole excitations derived from the configuration constructed by adding two
particles to the right of the ”Fermi momentum” as shown in Figure 1 (a). If we
add the two particles to the right and left of the ”Fermi momentum”
symmetrically, then the sum of the momenta vanishes.

Figure 1:  (Color online) Configuration of Bethe quantum numbers for a quantum
double dark-soliton state. (a) Configuration derived by adding two particles
(yellow circles) to that of the ground state (red circles). (b) A series of
configurations derived by punching two holes, which correspond to
$\bm{p}=\\{p_{1},p_{2}\\}$, in the configuration of the Bethe quantum numbers
constructed in panel (a). The third configuration has two holes at $p_{1}$ and
$p_{2}$.

The density profile of this state $\langle
X_{1},X_{2},N|\hat{\psi}^{\dagger}(x)\hat{\psi}(x)|X_{1},X_{2},N\rangle$ shows
the two density notches at the positions $x=X_{1},X_{2}$, which coincides with
the squared amplitude of the elliptic soliton [22]. Here, by the determinant
formula for the norms of Bethe eigenstates [35, 36] we can effectively
evaluate the matrix element

$\displaystyle\langle
X_{1},X_{2},N|\hat{\psi}^{\dagger}(x)\hat{\psi}(x)|X_{1},X_{2},N\rangle=\frac{1}{\mathcal{M}_{N}}\sum_{\bm{p},\bm{p^{\prime}}\in\bm{P}_{N}}e^{i(P-P^{\prime})x}e^{i(p_{1}X_{1}+p_{2}X_{2})}$
$\displaystyle\times e^{-i(p^{\prime}_{1}X_{1}+p^{\prime}_{2}X_{2})}\langle
p^{\prime}_{1},p^{\prime}_{2},N|\hat{\psi}^{\dagger}(0)\hat{\psi}(0)|p_{1},p_{2},N\rangle.$
(7)

Here, $P$ and $P^{\prime}$ in an exponential term denote the total momentum of
the state $|p_{1},p_{2},N\rangle$ and
$|p^{\prime}_{1},p^{\prime}_{2},N\rangle$ calculated through Equation (4),
respectively. The sum in the above equation is taken over all pairs of
$\bm{p}=\\{p_{1},p_{2}\\}$ and
$\bm{p^{\prime}}=\\{p^{\prime}_{1},p^{\prime}_{2}\\}$ that belong to the set
$\bm{P}_{N}$. The matrix element of the form factors of the density operator
[27, 37, 38] is given by

$\displaystyle\langle
p^{\prime}_{1},p^{\prime}_{2},N|\hat{\psi}^{\dagger}(0)\hat{\psi}(0)|p_{1},p_{2},N\rangle$
$\displaystyle=(-1)^{N(N+1)/2}(P-P^{\prime})\left(\prod^{N}_{j,\ell=1}\frac{1}{k^{\prime}_{j}-k_{\ell}}\right)\left(\prod^{N}_{j>\ell}k_{j,\ell}k^{\prime}_{j,\ell}\sqrt{\frac{{\hat{K}}(k^{\prime}_{j,\ell})}{{\hat{K}}(k_{j,\ell})}}\right)$
$\displaystyle\times\frac{\det U(k,k^{\prime})}{\sqrt{\det G(k)\det
G(k^{\prime})}},$ (8)

where the quasimomenta $\\{k_{1},\cdots,k_{N}\\}$ and
$\\{k^{\prime}_{1},\cdots,k^{\prime}_{N}\\}$ give the eigenstates
$|p_{1},p_{2},N\rangle$ and $|p^{\prime}_{1},p^{\prime}_{2},N\rangle$,
respectively. We use the abbreviations $k_{j,\ell}:=k_{j}-k_{\ell}$ and
$k^{\prime}_{j,\ell}:=k^{\prime}_{j}-k^{\prime}_{\ell}$. The kernel
${\hat{K}}(k)$ is defined by ${\hat{K}}(k)=2c/(k^{2}+c^{2})$. The matrix
$G(k)$ is called the Gaudin matrix, whose $(j,\ell)$ th element is given by

$G(k)_{j,\ell}=\delta_{j,\ell}\left[L+\sum_{m=1}^{N}{\hat{K}}(k_{j,m})\right]-{\hat{K}}(k_{j,\ell})\quad\mbox{for}\quad
j,\ell=1,2,\cdots,N.$ (9)

The matrix elements of the $(N-1)$ by $(N-1)$ matrix $U(k,k^{\prime})$ are
given by

$\displaystyle U(k,k^{\prime})_{j,\ell}$
$\displaystyle=2\delta_{j\ell}\textrm{Im}\left[\prod^{N}_{a=1}\frac{k^{\prime}_{a}-k_{j}+ic}{k_{a}-k_{j}+ic}\right]$
(10)
$\displaystyle+\frac{\prod^{N}_{a=1}(k^{\prime}_{a}-k_{j})}{\prod^{N}_{a\neq
j}(k_{a}-k_{j})}\left({\hat{K}}(k_{j,\ell})-{\hat{K}}(k_{N,\ell})\right).$

We have also considered the matrix element of the single field operator

$\displaystyle\langle X_{1},X_{2},N-1|\hat{\psi}(x)|X_{1},X_{2},N\rangle$
$\displaystyle=\frac{1}{\sqrt{\mathcal{M}_{N}\mathcal{M}_{N-1}}}\sum_{\bm{p^{\prime}}\in\bm{P}_{N-1}}\sum_{\bm{p}\in\bm{P}_{N}}e^{i(P-P^{\prime})x}e^{i(p_{1}X_{1}+p_{2}X_{2})}e^{-i(p^{\prime}_{1}X_{1}+p^{\prime}_{2}X_{2})}$
$\displaystyle\times\langle
p^{\prime}_{1},p^{\prime}_{2},N-1|\hat{\psi}(0)|p_{1},p_{2},N\rangle,$ (11)

where $P$ and $P^{\prime}$ denote the total momenta of the state
$|p_{1},p_{2},N\rangle$ and $|p^{\prime}_{1},p^{\prime}_{2},N-1\rangle$,
respectively. The determinant formula is given by [35, 36, 27, 39, 37, 38]

$\displaystyle\langle
p^{\prime}_{1},p^{\prime}_{2},N|\hat{\psi}(0)|p_{1},p_{2},N\rangle=(-1)^{N(N+1)/2+1}\left(\prod^{N-1}_{j=1}\prod^{N}_{\ell=1}\frac{1}{k^{\prime}_{j}-k_{\ell}}\right)$
$\displaystyle\times$
$\displaystyle\left(\prod^{N}_{j>\ell}k_{j,\ell}\sqrt{k_{j,\ell}^{2}+c^{2}}\right)\left(\prod^{N-1}_{j>\ell}\frac{k^{\prime}_{j,\ell}}{\sqrt{(k^{\prime}_{j,\ell})^{2}+c^{2}}}\right)\frac{\det\widehat{U}(k,k^{\prime})}{\sqrt{\det
G(k)\det G(k^{\prime})}},$ (12)

where the quasi-momenta $\\{k_{1},\cdots,k_{N}\\}$ and
$\\{k^{\prime}_{1},\cdots,k^{\prime}_{N-1}\\}$ give the eigenstates
$|p_{1},p_{2},N\rangle$ and $|p^{\prime}_{1},p^{\prime}_{2},N-1\rangle$,
respectively. We recall that the matrix $G(k)$ denotes the Gaudin matrix,
whose $(j,\ell)$th element is given in Equation (9). The matrix elements of
the $(N-1)$ by $(N-1)$ matrix $\widehat{U}(k,k^{\prime})$ are given by

$\displaystyle\widehat{U}(k,k^{\prime})_{j,\ell}$
$\displaystyle=2\delta_{j\ell}\textrm{Im}\left[\frac{\prod^{N-1}_{a=1}(k^{\prime}_{a}-k_{j}+ic)}{\prod^{N}_{a=1}(k_{a}-k_{j}+ic)}\right]$
(13)
$\displaystyle+\frac{\prod^{N-1}_{a=1}(k^{\prime}_{a}-k_{j})}{\prod^{N}_{a\neq
j}(k_{a}-k_{j})}\left({\hat{K}}(k_{j,\ell})-{\hat{K}}(k_{N,\ell})\right).$

###  2.4 One-particle reduced density matrix

The matrix element of the one-particle reduced density matrix,
$\rho_{1}(x,y):=\langle x|\hat{\rho}_{1}|y\rangle$, for a quantum system is
expressed as a correlation function in the ground state $|\lambda\rangle$:

$\displaystyle\rho_{1}(x,y)=\langle\lambda|\hat{\psi}^{\dagger}(y)\hat{\psi}(x)|\lambda\rangle.$
(14)

In the LL model we can numerically evaluate the correlation function by the
form factor expansion. Inserting the complete system of eigenstates,
$\sum_{\mu}|\mu\rangle\langle\mu|$, we have

$\displaystyle\rho_{1}(x,y)=\sum_{\mu}e^{{\rm
i}(P_{\mu}-P_{\lambda})(y-x)}|\langle\mu|\hat{\psi}(0)|\lambda\rangle|^{2},$
(15)

where $P_{\mu}$ denotes the momentum eigenvalues of eigenstates $|\mu\rangle$.
Each form factor in the sum (15) is expressed as a product of determinants by
making use of the determinant formula for the norms of Bethe eigenstates [35]
and that for the form factors of the field operator [27, 38, 37]:

$\displaystyle\langle\mu|\hat{\psi}(0)|\lambda\rangle=(-1)^{N(N+1)/2+1}\left(\prod^{N-1}_{j=1}\prod^{N}_{\ell=1}\frac{1}{k^{\prime}_{j}-k_{\ell}}\right)\left(\prod^{N}_{j>\ell}k_{j,\ell}\sqrt{k_{j,\ell}^{2}+c^{2}}\right)$
$\displaystyle\times\left(\prod^{N-1}_{j>\ell}\frac{k^{\prime}_{j,\ell}}{\sqrt{(k^{\prime}_{j,\ell})^{2}+c^{2}}}\right)\frac{\det
U(k,k^{\prime})}{\sqrt{\det G(k)\det G(k^{\prime})}}\,,$ (16)

where the quasi-momenta $\\{k_{1},\cdots,k_{N}\\}$ and
$\\{k^{\prime}_{1},\cdots,k^{\prime}_{N-1}\\}$ give the eigenstates
$|\lambda\rangle$ and $|\mu\rangle$, respectively. Here we have employed the
abbreviated symbols $k_{j,\ell}:=k_{j}-k_{\ell}$ and
$k^{\prime}_{j,\ell}:=k^{\prime}_{j}-k^{\prime}_{\ell}$. The matrix $G(k)$ is
the Gaudin matrix, whose $(j,\ell)$th element is given by
$G(k)_{j,\ell}=\delta_{j,\ell}\left[L+\sum_{m=1}^{N}K(k_{j,m})\right]-K(k_{j,\ell})$
for $j,\ell=1,2,\cdots,N$, where the kernel $K(k)$ is defined by
$K(k)=2c/(k^{2}+c^{2})$. The matrix elements of the $(N-1)$ by $(N-1)$ matrix
$U(k,k^{\prime})$ are given by [37, 27, 38, 35]

$\displaystyle U(k,k^{\prime})_{j,\ell}$
$\displaystyle=2\delta_{j\ell}\textrm{Im}\left[\frac{\prod^{N-1}_{a=1}(k^{\prime}_{a}-k_{j}+ic)}{\prod^{N}_{a=1}(k_{a}-k_{j}+ic)}\right]$
(17)
$\displaystyle+\frac{\prod^{N-1}_{a=1}(k^{\prime}_{a}-k_{j})}{\prod^{N}_{a\neq
j}(k_{a}-k_{j})}\left(K(k_{j,\ell})-K(k_{N,\ell})\right).$

For the ground state $|\lambda\rangle$ we have shown that the sum of the form
factor expansion is almost saturated for the one-particle and one-hole (1p1h)
excitations together with two-particles and two-holes (2p2h) excitations. The
saturation rate is explicitly presented in Table 1 of Section 4.3. However,
for excited states the saturation rate has not been evaluated. It should be
technically nontrivial to evaluate it for excited states. For the quantum
states of double dark-solitons, we suggest that the saturation rate should be
close to one in the weak coupling case in the form factor expansion up to some
excitations with relatively small numbers of particles and holes. It is based
on the observation that the density profiles of quantum double dark-soliton
states are similar to those of the double dark-solitons of the GP equation, as
we shall show in Section 3.

###  2.5 Winding number

We introduce the winding number $J$ associated with solutions of the GP
equation under the periodic boundary conditions. Let us assume that a solution
of the GP equation $\phi(x)=\sqrt{\rho(z)}\exp[i\varphi(x)]$ satisfies the
periodic boundary conditions:

$\displaystyle\varphi(x+L)=\varphi(x)+2\pi J$ (18)

where $J$ is an arbitrary integer. The integer $J$ is called the winding
number [17, 18]. In the previous study, we constructed the quantum single
dark-soliton with a nonzero-winding number.

##  3 Dynamics of quantum double dark-soliton

###  3.1 Time evolution of quantum double dark-soliton state constructed with
equal weight

By making use of the time dependent field operator $\hat{\psi}(x,t)$, the
local density and the matrix element of the quantum state at a given time $t$
are expressed as follows.

$\displaystyle\rho_{Q}(x,t)$ $\displaystyle:=\langle
X_{1},X_{2},N|\hat{\rho}(x,t)|X_{1},X_{2},N\rangle$ (19)
$\displaystyle=\frac{1}{\mathcal{M}_{N}}\sum_{\bm{p},\bm{p^{\prime}}\in\bm{P}_{N}}e^{i(P-P^{\prime})x}e^{i(p_{1}X_{1}+p_{2}X_{2})}e^{-i(p^{\prime}_{1}X_{1}+p^{\prime}_{2}X_{2})}e^{\left[-\mathrm{i}\left(E-E^{\prime}\right)t\right]}$
$\displaystyle\times\langle
p_{1}^{\prime},p_{2}^{\prime},N|\hat{\rho}(0,0)|p_{1},p_{2},N\rangle,$
$\displaystyle\psi_{Q}(x,t)$ $\displaystyle:=\langle
X_{1},X_{2},N-1|\hat{\psi}(x,t)|X_{1},X_{2},N\rangle$ (20)
$\displaystyle=\frac{1}{\sqrt{\mathcal{M}_{N-1}\mathcal{M}_{N}}}\sum_{\bm{p^{\prime}}\in\bm{P}_{N-1}}\sum_{\bm{p}\in\bm{P}_{N}}e^{i(P-P^{\prime})x}e^{i(p_{1}X_{1}+p_{2}X_{2})}e^{-i(p^{\prime}_{1}X_{1}+p^{\prime}_{2}X_{2})}$
$\displaystyle\times e^{-i(E-E^{\prime})t}\langle
p_{1}^{\prime},p_{2}^{\prime},N-1|\hat{\psi}(0,0)|p_{1},p_{2},N\rangle,$

where $E$ is the energy of the state $|p_{1},p_{2},N\rangle$, and
$\hat{\rho}(x,t)=\hat{\psi}^{\dagger}(x,t)\hat{\psi}(x,t)$ denotes the local
density operator. We have obtained the exact expressions of the time evolution
in Equations (19) and (20) since the Bethe ansatz method gives the exact
energies for the quantum state $|X_{1},X_{2},N\rangle$.

####  3.1.1 Quantum dark-soliton located at $X_{1}=L/4$ and $X_{2}=3L/4$
initially

Figure 2 shows the time evolution of the density profile , i.e., the graph of
$\rho_{Q}(x,t)$ versus $x$ at a given time $t$, for the quantum double dark-
soliton state with initial soliton positions $X_{1}=\frac{L}{4}$ and
$X_{2}=\frac{3L}{4}$ under the periodic boundary conditions. We call the plot
in the left panel of Figure 2 the two-dimensional (2D) density plot of the
local density. Here, the value of the local density $\rho_{Q}(x,t)$ at
position $x$ and time $t$ is expressed by the brightness of the point at
$(x,t)$ in the space-time diagram, where the horizontal axis corresponds to
the $x$ coordinate, while the vertical axis to time $t$. In the right panels
of Figure 2 snapshots of the density profile of $\rho_{Q}(x,t)$ at $t=0,2,4$,
and 11 are plotted.

Figure 2:  (Color online) 2D density plot of local density $\rho_{Q}(x,t)$ in
the space-time diagram (left panel) and snapshots of the density profile of
$\rho_{Q}(x,t)$ at $t=0,\ 2,\ 4$, and 11, in panels (a), (b), (c), and (d),
respectively (right panels), for the quantum double dark-soliton state with
$N=L=20$ and $c=0.05$. The corresponding quantum state is given by Equation
(6) with $X_{1}=\frac{L}{4}$ and $X_{2}=\frac{3L}{4}$. In the left panel, the
vertical (horizontal) axis corresponds to time $t$ (coordinate $x$). Panel (a)
shows the density profile at $t=0$, which corresponds to the green dotted line
in the 2D density plot of the left panel. Panels (b), (c), and (d) show the
density profile at $t=2,4$, and 11, which correspond to the blue, pink, and
ocher dotted lines in the 2D density plot (left panel) , respectively. The
frames of the legends of panels (a), (b), (c), and (d) correspond to the
green, blue, pink, and ocher colors, respectively.

We note that the density profile shown in panel (a) of Figure 2 is identical
to the upper-left panel of Figure 9 for $c=0.05$ in Ref. [22]. In the latter
panel it was shown that the density profile of the quantum double dark-soliton
state completely coincides with the density profile of the elliptic double
dark-soliton solution of the GP equation. Thus, at $t=0$, the density profile
of the quantum double dark-soliton state coincides with that of the elliptic
soliton solution of the GP equation.

The positions of notches are expressed by the areas of the darker color in the
2D density plot at the left panel of Figure 2. The trajectories of the
positions of the two notches in the density profile are given by two parallel
linearly elongated regions in the diagram of time $t$ and coordinate $x$, as
shown in the left panel of Figure 2. Thus, the two notches moves at the same
velocity in the positive $x$ direction. In the snapshots of the density
profiles, the soliton notches are gradually filled, i.e., they become
shallower in time evolution, as shown in panels (a), (b), (c), and (d) of
Figure 2. That is, the distance between the bottoms of the notches is kept
constant through the time evolution, while the depths of the notches become
smaller. Here we have defined the depth of a notch by the difference between
the largest and smallest values in the density profile. For example, at
$t=11$, the notches are located at $x_{1}=1.9115$ and $x_{2}=11.9115$, and the
distance between the two notches is given by $\Delta x=x_{1}-x_{2}=10=L/2$,
which is equal to that of $t=0$.

It was reported in Ref. [40] that quantum double dark-solitons with notches of
almost the same depths can appear again after their depths of notches become
much smaller over a time scale of $1/c$. However, the quantum double dark-
soliton states constructed in the present research do not show this
reappearing or recurrent behavior in time evolution. Once the soliton notches
in the density profile are completely filled, i.e., their depths vanish, the
density profile remains flat and uniform in time evolution, as illustrated in
Figure 2. We note that the construction of the quantum soliton in Ref. [40] is
different from that of the present research, and also that the number of
particles is equal to $N=8$ in Ref. [40], which is smaller than $N=20$ for the
system in Figure 2.

Figure 3:  (Color online) 2D density plot of the square amplitude
$|\psi_{Q}(x,t)|^{2}$ of the matrix element $\psi_{Q}(x,t)$ in the space-time
diagram (left panel) and snapshots of the square amplitude profile (i.e.,
$|\psi_{Q}(x,t)|^{2}$ versus $x$) at $t=0,2,4$ and 11, in panels (a), (b),
(c), and (d), respectively, (right panels), for $N=L=20$ and $c=0.05$. The
corresponding quantum states are given by Equation (6) with
$X_{1}=\frac{L}{4}$ and $X_{2}=\frac{3L}{4}$ for $N-1$ and $N$ particles. In
the left panel, the vertical (horizontal) axis corresponds to time $t$
(coordinate $x$). Panels (a), (b), (c), and (d) exhibit the square amplitude
profile at $t=0,2,4$, and 11, which correspond to the green, blue, pink, and
ocher dotted lines in the 2D density plot (left panel).

The notches in the density profile of $\rho_{Q}(x,t)$ and those in the profile
of the square amplitude $|\psi_{Q}(x,t)|^{2}$ of matrix element
$\psi_{Q}(x,t)$ exhibit different decaying behaviors in time evolution. Figure
3 shows the time evolution of the square amplitude profile of matrix element
$\psi_{Q}(x,t)$ with initial soliton positions $X_{1}=\frac{L}{4}$ and
$X_{2}=\frac{3L}{4}$ under the periodic boundary conditions. The average
density is decreasing in the time evolution of the profile of the square
amplitude $|\psi_{Q}(x,t)|^{2}$ in Figure 3, while the notches in Figure 2 are
filled gradually. In the density profile, the average density is kept constant
as time $t$ increases, since the density is conserved as a whole for any time
$t$: $\int_{0}^{L}dx\rho_{Q}(x,t)=N$. On the other hand, we suggest that the
amplitude of the matrix element between the two different quantum states of
double dark-soliton should gradually decrease and finally vanish in time
evolution, since they have different energies and particle numbers.

In the 2D density plot at the left panel of Figure 3, the trajectories of
notches in the space-time diagram are depicted by linearly elongated parallel
regions with darker color . The values at the bottoms of the notches are
almost equal to zero constantly in time evolution in panels (a), (b), (c), and
(d) of Figure 3. Consequently, Figure 3 shows the trajectories of the notches
more clearly than Figure 2, as depicted in the 2D density plot at the left
panel.

The snapshots of the phase profile at different times in time evolution are
shown in Figure 4. Here we remark that the phase is given by the argument of
the matrix element of Equation (20) as a complex number. In Figure 4 the
abrupt jumps of the phase profile are located at the positions of the notches
in Figure 3. The abrupt jumps of the phase profile move with the same constant
velocity as the notches in the square amplitude profile. Furthermore, the
whole phase profile is gradually shifted toward the negative direction in time
evolution. Moreover, the shape of the phase profile as a whole remains the
same at least up to $t=40$.

Figure 4:  (Color online) Snapshots of the phase profile for the matrix
elements $\psi_{Q}(x,t)$ at $t=0,2,4$, and 11, for $N=L=20$ and $c=0.05$. The
corresponding quantum states are given by Equation (6) with
$X_{1}=\frac{L}{4}$ and $X_{2}=\frac{3L}{4}$ for $N-1$ and $N$ particles.

At the initial time $t=0$, the profiles of the square amplitude and the phase
of the matrix element $\psi_{Q}(x,t)$ shown in panels (a) of Figure 3 and
Figure 4 are identical to those of Figures 10 and 11 for $c=0.05$ in Ref.
[22], respectively. Panel (a) of Figure 3, the square amplitude profile of the
matrix element, corresponds to the panel of $c=0.05$ in Figure 10 of Ref.
[22], where it was shown that the square amplitude profile of the classical
and quantum double dark-soliton overlap completely. Panel (a) of Figure 4, the
phase profile of the matrix element, corresponds to the panel of $c=0.05$ in
Figure 11 of Ref. [22], where the phase profiles of the classical and quantum
double dark-solitons overlap completely.

However, the time evolution of the phase profile in the quantum double dark-
soliton state is different from that of the elliptic dark-soliton solution,
which is given by the travelling wave solution of the GP equation. We recall
that the phase profile in the quantum double dark-soliton is gradually shifted
toward the negative direction in time evolution in Figure 4, while the phase
profile of the travelling wave solution is not shifted. Thus, the time
evolution of the quantum dark-solitons that we have constructed is slightly
different from that of the classical elliptic soliton solution.

We remark that two notches have mostly the same velocity as shown in Figures 2
and 3 for the quantum double dark-soliton constructed with equal weight. In
Section 3.2 we shall show that two notches have different velocities for the
quantum double dark-soliton state constructed with the Gaussian weights.

####  3.1.2 Quantum dark-soliton positions located at $X_{1}=X_{2}=0$
initially

By placing the positions of the notches for the quantum dark-solitons $X_{1}$
and $X_{2}$ at the same point , the profiles of the density and square
amplitude derived in time evolution are plotted in Figures 5 and 6,
respectively. In both profiles of the density and the square amplitude it
seems as if the two notches repel each other in time evolution.

Figure 5:  (Color online) 2D density plot of local density $\rho_{Q}(x,t)$ in
the space-time diagram (left panel) and snapshots of the density profile of
$\rho_{Q}(x,t)$ at $t=0,\ 2,\ 4$ and 11 (right panels), for $N=L=20$ and
$c=0.05$. The corresponding quantum state is given by Equation (6) with
$X_{1}=X_{2}=0$. Panels (a), (b), (c), and (d) show the density profile at
$t=0,\ 2,\ 4$, and 11, which correspond to the green, blue, pink, and ocher
dotted lines in the 2D density plot (left panel), respectively. The frames of
the legends of panels (a), (b), (c), and (d) correspond to the green, blue,
pink, and ocher colors, respectively.

The quantum double dark-soliton state with overlapping positions of two
notches has different properties in the profiles of the density and the square
amplitude from the quantum double dark-soliton state in Equation (6) with
different initial positions of two notches as $X_{1}=\frac{L}{4}$ and
$X_{2}=\frac{3L}{4}$. In the density profile of Figure 5 the notches are much
deeper than those in Figure 2, similarly as the notches in Figure 6 of the
square amplitude profile. In the profile of the square amplitude , the values
at the bottoms of the notches increase in time evolution: The values of the
square amplitude at the bottoms of notches are not close to zero at $t=11$ in
Figure 6. Thus, for the quantum double dark-soliton with two overlapping
positions of notches the difference between the density profile and the square
amplitude profile is smaller than in Figure 2 and Figure 3.

Figure 6:  (Color online) 2D density plot of the square amplitude
$|\psi_{Q}(x,t)|^{2}$ of the matrix element $\psi_{Q}(x,t)$ (left panel) and
snapshots of the profile of the square amplitude (i.e., $|\psi_{Q}(x,t)|^{2}$
versus $x$) at $t=0,2,4$ and 11 (right panels), for $N=L=20$ and $c=0.05$. The
corresponding quantum states are given by Equation (6) with $X_{1}=X_{2}=0$
for $N-1$ and $N$ particles. Panel (a) shows the density profile at $t=0$,
which corresponds to the green dotted line of the 2D density plot in the left
panel. Panels (b), (c), and (d) show the density profile at $t=2,4$, and 11,
which correspond to the blue, pink, and ocher dotted lines of the 2D density
plot in the left panel, respectively. The frames of the legends of panels (a),
(b), (c), and (d) correspond to the green, blue, pink, and ocher colors,
respectively. Figure 7:  (Color online) Snapshots in the time evolution of the
phase profile for the matrix element $\psi_{Q}(x,t)$ at $t=0,2,4$ and 11 for
$N=L=20$ and $c=0.05$. The corresponding quantum states are given by Equation
(6) with $X_{1}=X_{2}=0$ for $N-1$ and $N$ particles.

The snapshots of the phase profile in time evolution are exhibited in Figure 7
for the quantum double dark-soliton state with two overlapping positions of
notches. In each panel, the phase profile satisfies the boundary condition
specified by a nonzero winding number $J$:
$\textrm{Arg}[\psi_{Q}(x+L,t)]=\textrm{Arg}[\psi_{Q}(x,t)]+2\pi J$
($J\in\mathbb{Z}$). At $t=0$, the winding number was given by $J=2$, while it
suddenly changed to $J=1$ at $t=0.05$. After the change of the winding number,
the phase profile became smoother in shape gradually in time evolution.
Furthermore, we observe in Figure 7 that the whole phase profile was shifted
toward the negative direction step-by-step in time evolution. It is also the
case in Figure 4: The whole phase profile was shifted in the negative
direction for the quantum double dark-soliton state with initial positions of
notches placed at $X_{1}=\frac{L}{4}$ and $X_{2}=\frac{3L}{4}$.

The abrupt change of the winding number may occur in time evolution for the
phase profile associated with the quantum states, i.e., the phase profile of
the matrix element of the field operator between the quantum double dark-
soliton states in Equation (20). The boundary condition of the phase is given
by the form of Equation (18) for solutions of the GP equations and also for
the phase profile associated with the quantum states in Equation (20).
However, the quantum states do not depend on the boundary conditions of
classical solutions. It is sufficient if the phase profile associated with the
quantum states satisfies one of the boundary conditions of Equation (18)
specified by an integer $J$, which we have called the winding number. Thus,
the winding number $J$ may change abruptly in time evolution in the phase
profile associated with the quantum double dark-soliton states in Equation
(20).

###  3.2 Time evolution of quantum double dark-soliton state with the ideal
Gaussian weights

Let us consider the Gaussian weighted superposition of the excited states
consisting of two particle-hole excitations which are determined by a pair of
holes $\bm{p}=\\{p_{1},p_{2}\\}$ in the set $\bm{P}$:

$\displaystyle|X_{1},X_{2},N\rangle_{G}=\frac{1}{\sqrt{\mathcal{N}}}\sum_{\bm{p}\in\bm{P}}G_{P_{0}\sigma_{0}}(p_{1})G_{P^{\prime}_{0}\sigma^{\prime}_{0}}(p_{2})e^{i(p_{1}X_{1}+p_{2}X_{2})}\left|p_{1},p_{2},N\right\rangle.$
(21)

Here, $\mathcal{N}$ is a normalization factor and the set $\bm{P}$ is the same
as given in Section 2.3. The Gaussian function is given by

$\displaystyle
G_{P\sigma}(q)=\exp\left[-\frac{\left(q-P\right)^{2}}{4\sigma^{2}}\right]$
(22)

with two Gaussian parameters $(P,\sigma)$ [21]. The parameters $P$ and
$\sigma$ are determined by the target soliton depth $d$ and the density
$n=N/L$:

$\displaystyle P(d)$
$\displaystyle=2n\left\\{\frac{\pi}{2}-\left[\sqrt{\frac{d}{n}\left(1-\frac{d}{n}\right)}+\arcsin{\left(\sqrt{\frac{d}{n}}\right)}\right]\right\\}\,,$
(23) $\displaystyle\sigma^{2}(d)$
$\displaystyle=\frac{4}{3}n\sqrt{nc}\left(1-\frac{d}{n}\right)^{\frac{3}{2}}.$
(24)

Here we have defined the soliton depth $d$ by the smallest value in the
density profile of a single dark-soliton. It is different from the “depth of a
notch” defined in Section 3.1.1. The target soliton depth $d$ is expressed
with the dark soliton solution to the GP equation moving with velocity $v$ in
the thermodynamic limit $\phi_{P}^{\infty}(x)$ [21]:

$\displaystyle\sqrt{d}=|\phi_{P}^{\infty}(x=0)|=\sqrt{n}\left|\frac{v}{v_{c,\infty}}\right|.$
(25)

Here $|\phi_{P}^{\infty}(x=0)|$ denotes the square root of the local density
at the origin, which is the position of the notch in the thermodynamic limit,
and $v_{c,\infty}$ is called the critical velocity of the infinite system.
When the system size $L$ is finite, the largest velocity of the elliptic dark-
soliton solution of the GP equation is denoted by the critical velocity
$v_{c}$[22]. It approaches the critical value $v_{c,\infty}$ in the limit of
sending the system size $L$ to $\infty$.

The exact profiles in time evolution are numerically derived for the local
density $\rho_{Q}(x,t)$ and the square amplitude of the matrix element
$\psi_{Q}(x,t)$ of the field operator by calculating the time-dependent matrix
elements of the field operator between the Gaussian weighted quantum states of
Equation (21), similarly as we have demonstrated in Equation (19) and Equation
(20) of Section 3.1 for the quantum double dark-soliton state constructed with
equal weight. For the Gaussian weighted quantum double dark-soliton state, by
assigning a pair of proper values of the target soliton depth $d$ to the two
notches of a given superposition of quantum states of Equation (21), we can
construct a quantum double dark-soliton state such that its density profile
has two distinct notches with different depths.

We have constructed several quantum double dark-soliton states in which the
density profile has two distinct notches with different depths. In Figures 8,
9, 10, 11, 12 and 13 we set the target soliton depths as $d=0.6$ and $d=0.0$
to the two notches, respectively, and we generated the Gaussian weights by
making use of Equation (22). Here, the corresponding Gaussian parameters are
given by $(P_{0},\sigma)=(0.124027\pi,0.106667)$ and
$(P^{\prime}_{0},\sigma^{\prime})=(\pi,0.421637)$, respectively, which are
derived by making use of Equations (23) and (24). We have thus obtained the
quantum double dark-soliton state of distinct narrow notches with different
depths. Here we recall that single dark-solitons with different depths have
different speeds in the same direction for the GP equation.

We observe the scattering of two notches in the density and phase profiles of
the quantum double dark-soliton state. It exhibits the phase shift which is a
characteristic property in soliton-soliton collisions [41, 23], as shown in
the density profile. We remark that the 2D density plot of the local density
in the space-time diagram and the snapshots of the density profile at
different times are presented in Figure 8 for the quantum double dark-soliton
state constructed with the Gaussian weights for $c=0.05$. As the two notches
of the double dark-soliton approached each other, they moved along
approximately straight and linear trajectories with different constant
velocities. The collision occurred around at a time interval including $t=11$
(see panel (c), which corresponds to the pink dotted line in the left panel of
Figure 8). After the collision, each of the dark solitons travelled at the
same velocity before the collision. Furthermore, we confirm that the phase
shift occurred after the collision in the left panel of Figure 8.

Figure 8:  (Color online) 2D density plot of local density $\rho_{Q}(x,t)$ in
the space-time diagram (left panel) and the snapshots of the density profile
for the quantum double dark-soliton state constructed with the Gaussian
weights for $N=L=20$ and $c=0.05$. The corresponding quantum state is given by
Equation (21) with $X_{1}=L/4$ and $X_{2}=3L/4$. Panels (a), (b), (c), and (d)
show the density profile at $t=0,7,11$, and 21, which correspond to the green,
blue, pink, and ocher dotted lines in the left panel of the 2D density plot,
respectively. The frames of the legends of panels (a), (b), (c), and (d)
correspond to the green, blue, pink, and ocher colors, respectively.

Let us investigate the phase shift explicitly. By applying the Galilean
transformation, that is, in Figure 9 we observe the scattering process in the
inertial frame of reference moving with the left-hand-side notch of the
quantum double dark-soliton in Figure 8. We clearly confirm the phase shift
after the collision as shown in Figure 9. The position of the deeper notch is
shifted from $x=5$ to $x=2$ as shown in panel (d) of Figure 9. The position
shift of the deeper notch corresponds to the phase shift due to the collision
of the two dark-solitons.

Figure 9:  (Color online) 2D density plot of local density $\rho_{Q}(x,t)$ in
the space-time diagram derived by the Galilean transformation (left panel) and
corresponding snapshots of the density profile at $t=6,12,14$, and 20 in the
quantum double dark-soliton state constructed with the Gaussian weights for
$N=L=20$ and $c=0.05$. The corresponding quantum state is given by Equation
(21) with $X_{1}=-L/4$ and $X_{2}=L/4$. Panels (a), (b), (c), and (d) show the
density profile at $t=6,12,14$, and 22, which correspond to the green, blue,
pink, and ocher dotted lines in the left panel of the 2D density plot,
respectively. The frames of the legends of panels (a), (b), (c), and (d)
correspond to the green, blue, pink, and ocher colors, respectively.

Figure 10 shows the time evolution of the square amplitude profile of the
matrix element of the field operator for the quantum double dark-soliton
states constructed with the Gaussian weights for $c=0.05$. The quantum state
is the same as that of Figure 8. We have constructed the double dark-solitons
of distinct narrow notches with different depths not only in the density
profile but also in the square amplitude profile, i.e., the graph of
$|\psi_{Q}(x,t)|^{2}$ versus $x$. We observe the scattering of two notches in
the quantum double dark-soliton states. As the two notches of the double dark-
soliton states approached each other, they moved along approximately straight
and linear trajectories with different constant velocities, as shown in Figure
10. The collision occurred in the time interval including $t=11$ (see panel
(c), which corresponds to the pink dotted line in the left panel of Figure
10). After the collision, each of the dark solitons travelled with the same
velocity before the collision. We observe at least approximately the same
phase shift as shown in Figure 8. We remark that panel (a) of Figure 10, the
square amplitude profile of the matrix element for the quantum states
constructed with the Gaussian weights, corresponds to the panel of $c=0.05$ in
Figure 13 of Ref. [22].

We now demonstrate that the winding number changed during the scattering
process in the time evolution of the Gaussian weighted quantum double dark-
soliton states. Figure 11 shows the time evolution of the phase profile. In
each panel, the phase profile satisfies the boundary condition:
$\textrm{Arg}[\phi_{Q}(x+L,t)]=\textrm{Arg}[\phi_{Q}(x,t)]+2\pi J$ with a
winding number $J$. At the initial time $t=0$, the two notches of the quantum
dark-soliton were located at the most distant points from each other such as
$X_{1}=L/4$ and $X_{2}=3L/4$, and the winding number is given by $J=1$. When
the two notches of the quantum dark-soliton states became very close in space,
the winding number was suddenly changed to $J=0$, in the time interval
including $t=11$, as shown in panel (c) of Figure 11. After the collision, the
winding number was recovered: The winding number at $t=21$ was given by $J=1$,
as shown in panel (d) of Figure 11. We remark that panel (a) of Figure 11, the
phase profile of the matrix element between the Gaussian weighted quantum
double dark-soliton states, corresponds to the panel of $c=0.05$ in Figure 14
of Ref. [22].

Figure 10:  (Color online) 2D density plot of the square amplitude of the
matrix element $\psi_{Q}(x,t)$ and the snapshots of the square amplitude
profile at $t=0,7,11$, and 21, in the Gaussian weighted quantum double dark-
soliton state for $N=L=20$ and $c=0.05$. The corresponding quantum states are
given by Equation (21) with $X_{1}=\frac{L}{4}$ and $X_{2}=\frac{3L}{4}$ for
$N-1$ and $N$ particles. Panels (a), (b), (c), and (d) show the density
profile at $t=0,7,11$, and 21, which correspond to the green, blue, pink, and
ocher dotted lines in the left panel of the 2D density plot, respectively. The
frames of the legends of panels (a), (b), (c), and (d) correspond to the
green, blue, pink, and ocher colors, respectively. Figure 11:  (Color online)
Snapshots of the phase profile of the matrix element $\psi_{Q}(x,t)$ at
$t=0,7,11$, and 21 for the Gaussian weighted quantum double dark-soliton
states with $N=L=20$ and $c=0.05$. The corresponding quantum states are given
by Equation (21) with $X_{1}=\frac{L}{4}$ and $X_{2}=\frac{3L}{4}$ for $N-1$
and $N$ particles.

We explicitly evaluate the phase shift due to the scattering of two notches in
the quantum double dark-soliton state. The left panel of Figure 12 shows the
square amplitude profile of the matrix element $\psi_{Q}(x,t)$ in time
evolution observed in the inertial frame of reference moving together with the
deeper notch of the quantum double dark-soliton. The abrupt increase (or
decrease) in the phase profile, which we call a phase jump, was located at the
position of the deeper notch of the double dark-soliton, as shown in panels
(a), (b), and (d) of Figure 12: It was located at $x=5$ in panels (a) and (b),
and at $x=2$ in panel (d). Thus, the position of the deeper notch in the
double dark-soliton was shifted after the collision in the inertial frame of
reference. It corresponds to the phase shift due to the scattering of the two
notches.

Let us investigate the changes of the winding number in time evolution in
detail. The winding number $J$ was equal to zero when the two notches of the
quantum double dark-soliton were close to each other in space, as shown in
panel (c) of Figure 12. Figure 13 exhibits that the abrupt changes of the
winding number $J$ from 1 to 0 and from 0 to 1 occurred before and after the
process of the soliton collision. At time $t=7.5$ the value of the phase
abruptly increases in space with respect to the $x$ coordinate at the position
$x=5$ of the phase jump, while at time $t=7.6$ it suddenly decreases in space
at the position $x=5$ of the phase jump. Similarly, the corresponding change
of the winding number occurred during the time interval between $t=20.6$ and
$20.7$ after the collision.

We recall that it is not necessary for the winding number in the phase profile
of a quantum state to be conserved during the time evolution of the quantum
system. The winding number is defined for the corresponding classical system,
i.e., the GP equation, or for the phase profile of the quantum system. The
dynamics of the quantum system can be much more complex than the solutions of
the GP equation. When the two notches are far from each other in space, the
phase profile of the quantum system is similar to that of the classical
solution, while it is not the case when they collide with each other since
they are very close in space.

In summary, the Gaussian weighted superposition of the two-hole excited states
has lead to the quantum double dark-soliton states in which two notches have
different depths [22]. It follows that the notches of the quantum double dark-
soliton state have different velocities, and hence we have observed the
scattering of two notches in the quantum double dark-soliton state exactly. We
have also shown that the winding number of a quantum double dark-soliton state
changed when the two notches approach each other, explicitly for the Gaussian
weighted quantum double dark-soliton states.

We remark that one can make the quantum single dark-soliton black by making
use of the Gaussian weights, as shown in Ref. [21]. However, for the quantum
double dark-soliton, it seems that it is difficult to construct the double
black-soliton only by applying the Gaussian weights to the superposition of a
set of two-hole excitations.

Figure 12:  (Color online) 2D density plot of the square amplitude profile of
the matrix element $\psi_{Q}(x,t)$ (left panel), and the snapshots of the
phase profile for the matrix element $\psi_{Q}(x,t)$ at $t=0,7,14$, and 21
(right panels), in the Gaussian weighted quantum double dark-soliton state for
$N=L=20$ and $c=0.05$, observed in the inertial frame of reference moving with
the deeper notch, derived by the Galilean transformation. The corresponding
quantum state is given by Equation (21) with $X_{1}=-L/4$ and $X_{2}=L/4$.
Panels (a), (b), (c), and (d) show the phase profile at $t=0,7,14,21$, which
correspond to the green, blue, pink, and ocher dotted lines in the left panel
of the 2D density plot, respectively. The frames of the legends of panels (a),
(b), (c), and (d) correspond to the green, blue, pink, and ocher colors,
respectively. Figure 13:  (Color online) Snapshots in the time evolution of
the phase profile associated with the matrix element $\psi_{Q}(x,t)$ for the
Gaussian weighted quantum double dark-soliton state, observed in the inertial
frame of reference derived by the Galilean transformation: (a) before
collision and (b) after collision at $t=7.5$ and 7.6, and $t=20.6$ and 20.7,
respectively. We set $N=L=20$ and $c=0.05$. The corresponding quantum state is
given by Equation (21) with $X_{1}=-L/4$ and $X_{2}=L/4$.

##  4 Finite-Size Scaling of BEC

###  4.1 Motivation to study the quasi-BEC in 1D for the ground state

In 1D systems quantum fluctuations play a key role and often give subtle and
nontrivial effects. It is known that BEC occurs even for bosons with repulsive
interactions due to the quantum statistical effect among identical particles
[13]. In fact, the existence of BEC has been proven rigorously for interacting
bosons confined in dimensions greater than one [42]. In 1D case there is no
BEC for bosons with repulsive interactions due to strong quantum fluctuations
if we assume the standard thermodynamic limit with fixed coupling constant
[43]. On the other hand, if the coupling constant is very weak, we may expect
that even the 1D bosons with a large but finite number of particles undergo a
quasi-condensation in which “a macroscopic number of particles occupy a single
one-particle state” [13]. We call it a quasi-BEC by following the Penrose and
Onsager criterion.

However, it has not been shown explicitly how such a quasi-condensation occurs
in interacting bosons in one dimension. Furthermore, it is nontrivial to
expect it for the 1D Bose gas that is solvable by the Bethe ansatz. No pair of
particles can have the same quasi-momentum in common for a Bethe-ansatz
solution. Here we recall that we call the 1D system of bosons interacting with
repulsive delta-function potentials the 1D Bose gas. For the impenetrable 1D
Bose gas where the coupling constant is taken to infinity, condensate
fractions are analytically and numerically studied [44], while in the weak
coupling case it is nontrivial to evaluate the fractions in the 1D Bose gas.

We thus study in section 4 how the condensation fraction $n_{0}$, i.e., the
degree of the quasi-BEC, explicitly depends on the system size $L$, the number
of particles $N$ and the coupling constant $c$ in the ground state of the LL
model and particularly in the weak coupling case. It will be an illustrative
example.

###  4.2 Onsager-Penrose criterion of BEC

Let us review the definition of BEC through the one-particle reduced density
matrix for a quantum system [13, 14]. We assume that the number of particles
$N$ is very large but finite. At zero temperature, the density matrix is given
by $\hat{\rho}=|\lambda\rangle\langle\lambda|$, where $|\lambda\rangle$
denotes the ground state of the quantum system. We define the one-particle
reduced density matrix by the partial trace of the density matrix with respect
to other degrees of freedom: $\hat{\rho}_{1}=N\rm{tr}_{23\cdots N}\hat{\rho}$.
This matrix is positive definite and hence it is diagonalized as

$\displaystyle\hat{\rho}_{1}=N_{0}|\Psi_{0}\rangle\langle\Psi_{0}|+N_{1}|\Psi_{1}\rangle\langle\Psi_{1}|+\cdots.$
(26)

Here we put eigenvalues $N_{j}$ in descending order: $N_{0}\geq N_{1}\geq
N_{2}\geq\cdots>0$. The sum of all the eigenvalues is given by the number of
particles: $\sum_{j}N_{j}=N$. Here we recall ${\rm{tr}}_{1}\hat{\rho}_{1}=N$
due to the normalization: $\rm{tr}_{123\cdots N}\hat{\rho}=1$. Let us denote
by $n_{0}$ the ratio of the largest eigenvalue $N_{0}$ to particle number $N$:

$\displaystyle n_{0}:=N_{0}/N.$ (27)

The criterion of BEC due to Penrose and Onsager [14] is given as follows: If
the largest eigenvalue $N_{0}$ is of order $N$, i.e., the ratio $n_{0}$ is
nonzero and finite for large $N$, then we say that the system exhibits BEC,
and we call $n_{0}$ the condensate fraction. Here we also define fractions
$n_{j}$ by $n_{j}=N_{j}/N$ for $j=1,2,\ldots$.

###  4.3 Saturation rate in the form factor expansion at the ground state

$c$ | 0.01 | 1 | 100
---|---|---|---
1p1h | 0.999984 | 0.971538 | 0.693620
2p2h | $1.59454\times 10^{-5}$ | 0.0280102 | 0.289056
$n_{\rm sat}$ | 1.00000 | 0.999548 | 0.982676
Table 1:  Fraction $n_{\rm sat}$ of the reduced density operator at the
origin, $\rho_{1}(0,0)$, to the density $n$, evaluated by taking the sum over
a large number of eigenstates $|\mu\rangle$ with one particle and one hole
(1p1h) or with two particles and two holes (2p2h) for $N=L=50$ ($n=1$):
$n_{\rm
sat}=\left(\sum_{\mu}^{1p1h}+\sum_{\mu}^{2p2h}\right)|\langle\mu|\hat{\psi}(0)|\lambda\rangle|^{2}/n.$

Numerically we calculate correlation function in Equation (15) by taking the
sum over a large number of eigenstates with one particle and one hole (1p1h)
and those with two particles and two holes (2p2h). In order to confirm the
validity of the restricted sum, we have estimated the ratio of the one-
particle reduced density operator at the origin to density $n$,
$\rho_{1}(0,0)/n$, through the form factor expansion in Equation (15) for the
excitations with 1p1h or 2p2h. We express it by $n_{\rm sat}$. The estimates
of $n_{\rm sat}$ are listed in Table 1. The graph of $n_{\rm sat}$ approaches
1 for small coupling constant $c$, while it is larger than 0.98 for any value
of $c$ in the case of $N=50$.

###  4.4 Evaluation of the one-particle reduced density matrix of the ground
state

For the LL model, the eigenfunctions of the one-particle reduced density
matrix are given by plane waves for any nonzero and finite value of $c$. It is
a consequence of the translational invariance of the Hamiltonian of the LL
model. We thus have

$\displaystyle\rho_{1}(x,y)$
$\displaystyle=\frac{N_{0}}{L}+\sum_{j=1}^{\infty}\frac{2N_{j}}{L}\cos\left[2\pi
j(x-y)/L\right]\,.$ (28)

The eigenvalues of the one-particle reduced density matrix, $N_{j}$, are
expressed in terms of the form factor expansion. We consider the sum over all
the form factors between the ground state, $|\lambda\rangle$, and such
eigenstates, $|\mu\rangle$, that have a given momentum $P_{j}$ as

$\displaystyle
N_{j}=L\sum_{\mu:P_{\mu}=P_{j}}|\langle\mu|\hat{\psi}(0)|\lambda\rangle|^{2}\,.$
(29)

In the LL model we have $P_{j}:=(2\pi/L)j$.

Solving the Bethe ansatz equations for a large number of eigenstates we
observe numerically that eigenvalues $N_{j}$ are given in decreasing order
with respect to integer $j$: $N_{0}>N_{1}>N_{2}>\cdots$. It thus follows that
condensate fraction which corresponds to the largest eigenvalue of the one-
particle reduced density matrix ${\hat{\rho}}_{1}$ is indeed given by
$n_{0}=N_{0}/N$, where $N_{0}$ has been defined by the sum of Equation (29)
over all eigenstates with zero momentum.

###  4.5 Condensate fraction in the weak coupling regime

Figure 14: (Color online) Dependence of fractions $n_{j}$ on coupling constant
$c$. In the upper panel: condensate fraction $n_{0}$ is plotted against
coupling constant $c$ for $N=4$, 10, 20, 40, 100, 200 and 400, from the top to
the bottom, in red, green, blue, black, orange, purple and cyan lines,
respectively. In the lower panel: condensate fraction $n_{0}$, fractions
$n_{1}$ and $n_{2}$ are shown against $c$ from the top to the bottom in blue,
red and green lines, respectively, for $N=20$. We recall $n=N/L=1$.

The estimates of condensate fraction $n_{0}$ are plotted against coupling
constant $c$ in the upper panel of Figure 14 over a wide range of $c$ such as
from $c=10^{-3}$ to $c=10^{3}$ for different values of particle number $N$
such as $N=4$, 10, …, 400. For each $N$, condensate fraction $n_{0}$ becomes
1.0 for small $c$ such as $c<0.01$, while it decreases with respect to $c$ and
approaches an asymptotic value in the large $c$ region such as $c>100$ or
1000. The asymptotic values depend on particle number $N$ for $N=4$, 10, …,
400, and they are consistent with the numerical estimates of occupation
numbers for the impenetrable 1D Bose gas (see Equation (56) of Ref. [44]). In
the lower panel of Figure 14, we plot fractions $n_{j}$ for $j=0,1$ and 2
against coupling constant $c$ from $c=10^{-3}$ to $c=10^{3}$ with $N=20$. The
asymptotic values of $n_{j}$ for large $c$ (i.e. $c=1000$) are consistent with
the numerical estimates for the impenetrable 1D Bose gas (for $n_{1}$ and
$n_{2}$, see Equations (57) and (58) of Ref. [44], respectively).

Figure 15: (Color online) Condensate fraction $n_{0}$ as a function of $1/N$
for $c=0.01$. Here $n=N/L=1.0$.

We observe that condensate fraction $n_{0}$ decreases as particle number $N$
increases where density $n=N/L$ is fixed. It is the case for $c<0.1$ in the
upper panel of Figure 14. Condensate fraction $n_{0}$ decreases as $N$
increases even for small $c$ such as $c=0.01$, as shown in Figure 15. Thus, it
is necessary for coupling constant $c$ to decrease with respect to $N$ so that
condensate fraction $n_{0}$ remains constant as $N$ increases with fixed
density $n$.

###  4.6 Exact finite-size scaling

Figure 16: (Color online) Contours of condensate fraction $n_{0}$ are plotted
for various values of $n_{0}$ in the $\gamma$ versus $1/N$ plane. Each contour
is approximated by (30): $\gamma$ as a function of $1/N$.

We now show the finite-size scaling of condensate fraction $n_{0}$. In Figure
16 each contour line gives the graph of interaction parameter $\gamma$ as a
function of the inverse of particle number $N$ for a fixed value of condensate
fraction $n_{0}$. They are plotted for various values of $n_{0}$ from
$n_{0}=0.6$ to 0.99, and are obtained by solving the Bethe-ansatz equations
numerically.

For different values of density such as $n=1$, 2 and 5, we have plotted
contour lines with fixed values of condensate fraction $n_{0}$ in the plane of
interaction parameter $\gamma$ versus inverse particle number $1/N$. We have
observed that the contours with the same condensate fraction $n_{0}$ but for
the different densities coincided with each other in the $\gamma$ versus $1/N$
plane. Furthermore, they are well approximated by

$\gamma=A/N^{\eta}.$ (30)

Thus, condensate fraction $n_{0}$ is constant as particle number $N$ becomes
very large if interaction parameter $\gamma$ is given by the power of particle
number $N$ as in Equation (30).

Applying the finite-size scaling arguments, we suggest from Equation (30) that
condensation fraction $n_{0}$ is given by a scaling function $\phi(\cdot)$ of
a single variable $\gamma N^{\eta}$: $n_{0}=\phi(\gamma N^{\eta})$. Here we
recall the coincidence of contours for the different values of density $n$ in
Figure 16. We thus observe that exponent $\eta$ and amplitude $A$ of Equation
(30) are determined only by condensate fraction $n_{0}$ and are independent of
density $n$.

Let us consider amplitude $A$ as a function of $n_{0}$. We denote it by
$A=f(n_{0})$. Then, the scaling function $\phi(\cdot)$ is given by the inverse
function: $n_{0}=f^{-1}(A)$. In Figure 17, exponent $\eta$ increases with
respect to $n_{0}$, and amplitude $A$ decreases monotonically with respect to
$n_{0}$.

###  4.7 Quasi-BEC according to the Onsager-Penrose criterion

It follows from (30) that BEC does not occur in the 1D Bose gas if we fix
parameter $\gamma$ and density $n$ as system size $L$ goes to infinity.
However, if $\gamma$ is small enough so that it satisfies Equation (30) for a
given value of condensate fraction $n_{0}$, the 1D Bose gas shows the quasi-
BEC from the viewpoint of the Penrose and Onsager criterion. We suggest that
if condensate fraction $n_{0}$ of a quantum state is nonzero and finite for
large $N$, the mean-field approximation is valid for the quantum state. For
instance, there exist such quantum states that correspond to classical dark-
solitons of the GP equation [15], if parameter $\gamma$ is small enough so
that it satisfies Equation (30).

Figure 17: (Color online) Exponent $\eta$ and amplitude $A$ as functions of
condensate fraction $n_{0}$.

###  4.8 Various limiting procedures

With the scaling behavior expressed in Equation (30) we derive various ways of
the thermodynamic limit such that condensate fraction $n_{0}$ is constant. For
instance, we consider the case of a finite particle number, $N=N_{\rm f}$.
Choosing a value of $n_{0}$, we determine $\gamma$ by Equation (30) as
$\gamma=A(n_{0})/N_{\rm f}^{\eta(n_{0})}$. Then, the 1D Bose gas with
$N=N_{\rm f}$ has the same condensate fraction $n_{0}$ for any large value of
$L$ if coupling constant $c$ is given by $c=A(n_{0})N_{\rm f}^{1-\eta}/L$. Let
us set $\eta=1$ and $N_{\rm f}=10$, for simplicity. We have $n_{0}=0.97$ in
Figure 17, and $\gamma=0.3$ at $1/N=0.1$ in the contour of $n_{0}=0.97$ in
Figure 16. By assuming $n=1$, it corresponds to the case of $L=10$ and
$c=0.3$, and we have $A=cL=3$, which is consistent with Figure 17. Therefore,
the 1D Bose gas with $N_{\rm f}=10$ has $n_{0}=0.97$ for any large $L$ if $c$
is given by $c=0.3/L$. Moreover, we may consider other types of thermodynamic
limits. When density $n$ is proportional to a power of $L$ as $L^{\alpha}$,
condensate fraction $n_{0}$ is constant as $L$ goes to infinity if we set
$c\propto L^{(1-\eta)(1+\alpha)-1}$.

The scaling law in Equation (30) and the estimates of condensate fraction in
the present paper should be useful for estimating conditions in experiments of
trapped cold atomic gases in one dimension [45]. For instance, we suggest from
Figure 14 that BEC may appear in 1D systems with a small number of bosons such
as $N=20$ or $40$ for $c=1$ or $10$.

##  5 Concluding remarks

In the first part, we have shown that the density profile and the square
amplitude evolved in time differently, in particular, for the equal weight
case. In the former the notches were filled progressively, while the amplitude
of the latter decreased gradually. Furthermore, the Gaussian weights led to
the different depths for quantum double dark-solitons [22]. This gave the two
notches of the quantum double dark-soliton the different speeds, and we
observed the scattering of the two notches in the quantum double dark-soliton
state exactly. Interestingly, the winding number of the quantum double dark-
soliton state has changed when the two notches approach. Here we recall that
it is not necessary for the winding number to be conserved in the time
evolution of the quantum system, since it is defined for the corresponding
classical system.

In the second part, we exactly calculated the condensate fraction of the 1D
Bose gas with repulsive interaction by the form factor expansion for the
ground state. We have shown the finite-size scaling behavior such that
condensate fraction $n_{0}$ is given by a scaling function of interaction
parameter $\gamma$ times some power of particle number $N$: $n_{0}=\phi(\gamma
N^{\eta})$. Consequently, if parameter $\gamma$ decrease as
$\gamma=A/N^{\eta}$, condensate fraction $n_{0}$ remains nonzero and constant
as particle number $N$ becomes very large. By modifying the thermodynamic
limit, the 1D Bose gas shows BEC from the viewpoint of the Penrose-Onsager
criterion.

##  6 Acknowledgements

The present research is partially supported by Grant-in-Aid for Scientific
Research No. 21K03398. K. K. is supported by the Japan Science Technology
Agency (CREST Grant Number JPMJCR 19T4).

## References

## References

  * [1] Görlitz A, Vogels J M, Leanhardt A E, Raman C, Gustavson T L, Abo-Shaeer J R, Chikkatur A P, Gupta S, Inouye S, Rosenband T and Ketterle W 2001 Phys. Rev. Lett. 87(13) 130402 URL https://link.aps.org/doi/10.1103/PhysRevLett.87.130402
  * [2] Greiner M, Bloch I, Mandel O, Hänsch T W and Esslinger T 2001 Phys. Rev. Lett. 87(16) 160405 URL https://link.aps.org/doi/10.1103/PhysRevLett.87.160405
  * [3] Kinoshita T, Wenger T and Weiss D S 2004 Science 305 1125–1128 URL https://www.science.org/doi/abs/10.1126/science.1100700
  * [4] Weller A, Ronzheimer J P, Gross C, Esteve J, Oberthaler M K, Frantzeskakis D J, Theocharis G and Kevrekidis P G 2008 Phys. Rev. Lett. 101(13) 130401 URL https://link.aps.org/doi/10.1103/PhysRevLett.101.130401
  * [5] Becker C, Stellmer S, Soltan-Panahi P, Dörscher S, Baumert M, Richter E M, Kronjäger J, Bongs K and Sengstock K 2008 Nature Physics 4 496–501 URL https://doi.org/10.1038/nphys962
  * [6] Burger S, Bongs K, Dettmer S, Ertmer W, Sengstock K, Sanpera A, Shlyapnikov G V and Lewenstein M 1999 Phys. Rev. Lett. 83(25) 5198–5201 URL https://link.aps.org/doi/10.1103/PhysRevLett.83.5198
  * [7] Katsimiga G, Mistakidis S, Koutentakis G, Kevrekidis P and Schmelcher P 2017 New J. Phys. 19 123012 URL https://dx.doi.org/10.1088/1361-6455/ab787f
  * [8] Delande D and Sacha K 2014 Phys. Rev. Lett. 112(4) 040402 URL https://link.aps.org/doi/10.1103/PhysRevLett.112.040402
  * [9] Lieb E H and Liniger W 1963 Phys. Rev. 130(4) 1605–1616 URL https://link.aps.org/doi/10.1103/PhysRev.130.1605
  * [10] Lieb E H 1963 Phys. Rev. 130(4) 1616–1624 URL https://link.aps.org/doi/10.1103/PhysRev.130.1616
  * [11] Dorlas T C 1993 Commun. Math. Phys. 154 347–376 URL https://doi.org/10.1007/BF02097001
  * [12] Pitaevskii L and Stringari S 2003 Bose-Einstein Condensation and Superfluidity International Series of Monographs on Physics (Oxford University Press) ISBN 9780191076688 URL https://doi.org/10.1093/acprof:oso/9780198758884.001.0001
  * [13] Leggett A J 2006 Quantum Liquids: Bose condensation and Cooper pairing in condensed-matter systems (Oxford University Press) ISBN 9780198526438 URL https://doi.org/10.1093/acprof:oso/9780198526438.001.0001
  * [14] Penrose O and Onsager L 1956 Phys. Rev. 104(3) 576–584 URL https://link.aps.org/doi/10.1103/PhysRev.104.576
  * [15] Sato J, Kanamoto R, Kaminishi E and Deguchi T 2012 Phys. Rev. Lett. 108 110401 URL https://link.aps.org/doi/10.1103/PhysRevLett.108.110401
  * [16] Ishikawa M and Takayama H 1980 Journal of the Physical Society of Japan 49 1242–1246 URL https://doi.org/10.1143/JPSJ.49.1242
  * [17] Kanamoto R, Carr L D and Ueda M 2008 Phys. Rev. Lett. 100(6) 060401 URL https://link.aps.org/doi/10.1103/PhysRevLett.100.060401
  * [18] Kanamoto R, Carr L D and Ueda M 2009 Phys. Rev. A 79(6) 063616 URL https://link.aps.org/doi/10.1103/PhysRevA.79.063616
  * [19] Sato J, Kanamoto R, Kaminishi E and Deguchi T 2016 New J. Phys. 18 075008 URL https://doi.org/10.1088%2F1367-2630%2F18%2F7%2F075008
  * [20] Shamailov S S and Brand J 2019 Phys. Rev. A 99(4) 043632 URL https://link.aps.org/doi/10.1103/PhysRevA.99.043632
  * [21] Kaminishi E, Mori T and Miyashita S 2020 J. Phys. B: At. Mol. Opt. Phys 53 095302 URL https://dx.doi.org/10.1088/1361-6455/ab787f
  * [22] Kinjo K, Kaminishi E, Mori T, Sato J, Kanamoto R and Deguchi T 2022 Universe 8 2 URL https://doi.org/10.3390/universe8010002
  * [23] Ishiguro Y, Sato J, Ezaki T and Nishinari K 2022 Phys. Rev. Res. 4(3) L032047
  * [24] Syrwid A and Sacha K 2015 Phys. Rev. A 92(3) 032110 URL https://link.aps.org/doi/10.1103/PhysRevA.92.032110
  * [25] Syrwid A 2021 J. Phys. B: At. Mol. Opt. Phys 54 103001 URL https://doi.org/10.1088/1361-6455/abd37f
  * [26] Golletz W, Górecki W, Ołdziejewski R and Pawłowski K 2020 Phys. Rev. Res. 2(3) 033368 URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.033368
  * [27] Slavnov N A 1989 Theor. Math. Phys. 79 502–508 URL https://doi.org/10.1007/BF01016531
  * [28] Korepin V E, Bogoliubov N M and Izergin A G 1993 Quantum Inverse Scattering Method and Correlation Functions Cambridge Monographs on Mathematical Physics (Cambridge University Press)
  * [29] Kopyciński J, Łebek M, Marciniak M, Ołdziejewski R, Górecki W and Pawłowski K 2022 SciPost Phys. 12 023 URL https://scipost.org/10.21468/SciPostPhys.12.1.023
  * [30] Morera Navarro I, Guilleumas M, Mayol R and Mateo A M n 2018 Phys. Rev. A 98(4) 043612 URL https://link.aps.org/doi/10.1103/PhysRevA.98.043612
  * [31] Liang Z, Zhang Z and Liu W 2005 Phys. Rev. Lett. 94 050402 URL https://link.aps.org/doi/10.1103/PhysRevLett.94.050402
  * [32] Wang D S, Hu X H, Hu J and Liu W M 2010 Phys. Rev. A 81(2) 025604 URL https://link.aps.org/doi/10.1103/PhysRevA.81.025604
  * [33] Wen L, Li L, Li Z D, Song S W, Zhang X F and Liu W 2011 Eur. Phys. J. D 64 473–478 URL https://doi.org/10.1140/epjd/e2011-20485-4
  * [34] Li L, Li Z, Malomed B A, Mihalache D and Liu W 2005 Phys. Rev. A 72 033611 URL https://link.aps.org/doi/10.1103/PhysRevA.72.033611
  * [35] Gaudin M 2014 The Bethe Wavefunction (Cambridge University Press)
  * [36] Korepin V E 1982 Commun. Math. Phys. 86 391–418 URL https://doi.org/10.1007/BF01212176
  * [37] Caux J S, Calabrese P and Slavnov N A 2007 J. Stat. Mech.: Theory Exp 2007 P01008–P01008 URL https://doi.org/10.1088/1742-5468/2007/01/p01008
  * [38] Kojima T, Korepin V E and Slavnov N 1997 Commun. Math. Phys. 188 657–689
  * [39] Slavnov N A 1990 Theor. Math. Phys. 82 273–282 URL https://doi.org/10.1007/BF01029221
  * [40] Syrwid A, Brewczyk M, Gajda M and Sacha K 2016 Phys. Rev. A 94(2) 023623 URL https://link.aps.org/doi/10.1103/PhysRevA.94.023623
  * [41] Kaminishi E, Mori T and Ueda M 2019 Dynamics of the quanutm dark soliton The 74th Annual Meeting of the Physical Society of Japan
  * [42] Lieb E H and Seiringer R 2002 Phys. Rev. Lett. 88(17) 170409 URL https://link.aps.org/doi/10.1103/PhysRevLett.88.170409
  * [43] Pitaevskii L and Stringari S 1991 J. Low Temp. Phys. 85 377–388 URL https://doi.org/10.1007/BF00682193
  * [44] Forrester P J, Frankel N E, Garoni T M and Witte N S 2003 Phys. Rev. A 67(4) 043607 URL https://link.aps.org/doi/10.1103/PhysRevA.67.043607
  * [45] Pitaevskii L and Stringari S 2003 Bose-Einstein Condensation International series of monographs on physics (Oxford University Press) ISBN 0 19 850719-4