###### Abstract

We study quantum double dark-solitons by constructing corresponding quantum
states in the Lieb-Liniger model for the one-dimensional Bose gas. Here we
expect that the Gross-Pitaevskii (GP) equation should play a central role in
the long distance mean-field behavior of the 1D Bose gas. We first introduce
novel quantum states of a single dark soliton with a nonzero winding number.
We show them by exactly evaluating not only the density profile but also the
profiles of the square amplitude and phase of the matrix element of the field
operator between the $N$-particle and $(N-1)$-particle states. For elliptic
double dark-solitons, the density and phase profiles of the corresponding
states almost perfectly agree with those of the classical solutions,
respectively, in the weak coupling regime. We then show that the scheme of the
mean-field product state is quite effective for the quantum states of double
dark solitons. Assigning the ideal Gaussian weights to a sum of the excited
states with two particle-hole excitations we obtain double dark-solitons of
distinct narrow notches with different depths. We suggest that the mean-field
product state should be well approximated by the ideal Gaussian weighted sum
of the low excited states with a pair of particle-hole excitations. The
results of double dark-solitons should be fundamental and useful for
constructing quantum multiple dark-solitons.

###### keywords:

dark soliton; quantum soliton; 1D Bose gas

xx 1 5 Received: date; Accepted: date; Published: date Quantum dark solitons
in the 1D Bose gas: From single to double dark-solitons Kayo Kinjo1, Eriko
Kaminishi2, Takashi Mori3, Jun Sato4, Rina Kanamoto5 and Tetsuo Deguchi6 Kayo
Kinjo, Eriko Kaminishi, Takashi Mori, Jun Sato, Rina Kanamoto and Tetsuo
Deguchi Correspondence<EMAIL_ADDRESS>
##  1 Introduction

The Gross-Pitaevskii (GP) equation plays a fundamental role in the study of
interacting nonuniform gases realized in cold atoms such as the Bose Einstein
condensation (BEC) achieved in trapsPitaevskii and Stringari (2003). It is
expressed in terms of the classical complex scalar field $\psi$ as follows.

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

The non-linearity of the GP equation gives rise to a novel series of physical
phenomena such as grey solitons. The GP equation has an analogous form as the
Scrödinger equation with a nonlinear term of the fourth degree. In one
dimension it corresponds to the so-called nonlinear Scrödinger (NS) equation,
which is one of the fundamental integrable nonlinear equations characterized
by multi-soliton solutionsFaddeev and Takhtajan (1987). If the coupling
constant is negative the NS equation has localized soliton solutions, called
bright solitons, while if it is positive the NS equation has dark or grey
solitons exhibiting notches in the constant density solution Tsuzuki (1971).
Here we remark that bright solitons are actively studied in association with
optical fibersAbdullaev et al. (1993), while dark solitons are observed in
cold atomic experiments Burger et al. (1999); Weller et al. (2008); Becker et
al. (2008).

It should be quite an interesting but nontrivial problem to study how far the
GP equation can describe the behavior of quantum many-body systems such as the
BECs. The GP equation is a partial differential equation for the wave function
as a complex-valued order parameter, which is simply a classical complex
scalar field. It does not contain any quantum fluctuations in a
straightforward manner. In order to study the relevance or applicability of
the GP equation explicitly we may consider the Lieb-Liniger model. It
describes the Bose particles in one dimension interacting through the delta-
function interaction potentials. We briefly call it the one-dimensional (1D)
Bose gas. The GP equation corresponds to the Heisenberg equation of motion of
the field operator for the second-quantized Hamiltonian of the Lieb-Liniger
model where we assume that the system size is finite.

It is known that the NS equation plays an important role in the quantum
inverse scattering method, which can be regarded as a quantum version of the
inverse scattering method for classical integrable systems. However, it is not
trivial to construct such quantum states that can be related to solitons in
the NS equation even through a certain kind of limiting procedure. In the
attractive case, bright solitons are realized in terms of string-type bound
solutions of the Lieb-Liniger model Wadachi and Sakagami (1984). However, in
the repulsive case, the development of the studies was slightly more involved.
Ishikawa and Takayama argued that the dispersion relation of the type II
excitation of the 1D Bose gas in the repulsive case corresponds to that of a
dark (or grey) soliton solutionIshikawa and Takayama (1980). Here we remark
that some aspects of the connection were also addressed by several authors
previously. However, any quantum state which can be explicitly related to a
dark or grey soliton solution of the GP equation was not known for some period
of time. Recently, two different approaches to the quantum states associated
with dark soliton solutions have appeared. First, the quantum states which
have the same density and phase profiles with those of a dark soliton have
constructed Sato et al. (2012, 2016); Then, it was found that the density
profile of a single grey soliton appears in the quantum measurement process
Syrwid and Sacha (2015). Furthermore, the construction of quantum states of
grey solitons has been elaborated with the Gaussian weights Shamailov and
Brand (2019); Kaminishi et al. (2020). In particular, an ansatz on how to
assign the Gaussian weights to the yrast states was proposed Kaminishi et al.
(2020). The matrix elements of the bosonic field operator of the Lieb-Liniger
model are characterized by grey solitons Sato et al. (2016), and double
elliptic dark solitons are derived in association with the quantum measurement
process Golletz et al. (2020).

In this article we demonstrate that quantum states characterized by a double
dark or grey soliton are constructed in terms of the eigenstates with a pair
of particle-hole excitations in the Lieb-Liniger model. It thus gives
nontrivial examples where the description of the GP equation is effective in
the quantum system. We shall also show that the ansatz of Ref. Kaminishi et
al. (2020) is quite effective as far as the density profile is concerned. The
contents of the present article are given as follows. We first review in §2.1
that quantum states associated with the density and phase profiles of single
dark or grey solitons of the GP equation are constructed from the one-hole
excitations, namely, the type II excitations in the Lieb-Liniger model. They
are also called the yrast states. We then explain the ansatz Kaminishi et al.
(2020) for constructing quantum dark solitons in terms of the ideal Gaussian
weights in §2.2. We shall introduce the mean-field product state as the BEC
state corresponding to a macroscopic wave function. Here we assume that the
BEC state should exist for a given macroscopic wave function in the weak
coupling regime of a finite system, if the coupling constant is small enough.
In §3.1 we introduce a novel quantum state which leads to an elliptic dark-
soliton with a nonzero winding number in the phase profile of the matrix
element of the bosonic field operator between the quantum states. In §3.2 we
derive the elliptic multiple dark soliton solution of the GP equation. We then
argue analytically in §3.3 that the phase profile should be consistent with
that of the elliptic double dark-soliton in the limit of sending the system
size to infinity. We present in §4 the density profiles of quantum double
dark-solitons, explicitly. For an illustration we discuss in §4.1 the results
of the free fermions in one dimension. In §4.2 we show that the quantum state
given by the sum of the two particle-hole states with an equal weight leads to
the density and phase profiles, which almost perfectly agree with those of the
elliptic double dark-soliton solutions of the GP equation, respectively, in
the weak coupling regime. In §4.3, assigning the ideal Gaussian weights of the
ansatz Kaminishi et al. (2020) to the quantum double dark-soliton states we
obtain two distinct notches with different depths in the density profiles. We
also present the square amplitude and the phase profiles of the matrix element
of the bosonic field operator between the quantum double dark soliton states
with $N$ particles and $(N-1)$-particles. In the phase profiles of the quantum
double dark soliton states we observe that one of the two notches has a
nonzero winding number. Thus, the ansatz of the ideal Gaussian weights is also
important for constructing double dark soliton states. However, we should
point out that it does not completely lead to the mean-field product state. In
fact, one of the key properties of the mean-field product state holds only
approximately for all the examples investigated in this article. On the other
hand, we expect that one can improve the constructed quantum state with the
ansatz so that one obtain a complete mean-field product state by increasing
the number of particles $N$ and the system size $L$, while by decreasing the
coupling constant $c$.

We now give some remarks on mathematical aspects. The derivation of soliton
solutions is much harder in the repulsive case than in the attractive case for
the NS (or GP) equation as far as theoretical formulation is concerned. A dark
soliton solution of the GP equation was first derived by Tsuzuki for single
solitons Tsuzuki (1971). The multi-soliton solutions of the NS equation were
constructed by Zakharov-Shabat, in particular, in the case of attractive
nonlinear interaction Zakharov and Shabat (1973). For the repulsive case
multi-soliton solutions were systematically derived by Faddeev and Takhtajan
through the inverse scattering methodFaddeev and Takhtajan (1987). The multi
dark soliton solutions are called ”soliton solutions in the case of finite
density”, in the reference Faddeev and Takhtajan (1987). Here we should remark
that the gap-soliton solutions of the NS equation were derived through the
algebro-geometric approach, which generalize the multi solitons of the NS
equation into gap solitonsBelokolos et al. (1994). The formulation of the
inverse scattering method is quite general and useful, and the algebro-
geometric approach produces quite powerful results. However, both of them are
mathematically rather rich and it seems that they are not easily accessible
for standard theoretical physicists. Due to the mathematical richness they are
not often referred to in the physics literature, in particular, in the context
of cold atoms. It is thus one of the purposes of the present article as a
review to point out some explicit connections of some researches in
theoretical physics to the useful references in mathematics.

##  2 Quantum single dark soliton

###  2.1 Construction of quantum states of single dark solitons

In this section we review the method to construct quantum dark soliton states
in the case of a single soliton. Let us consider the one-dimensional Bose gas
interacting via the $\delta$-function potential Lieb and Liniger (1963)

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

where $N$ is the number of particle, $c$ is the coupling constant and $\mu$ is
the chemical potential. We impose a periodic boundary condition with length
$L$. In the second quantized form it reads

$\displaystyle\mathcal{\hat{H}}_{\text{LL}}=\int_{0}^{L}dx\left[\partial_{x}\hat{\psi}^{\dagger}\partial_{x}\hat{\psi}+c\hat{\psi}^{\dagger}\hat{\psi}^{\dagger}\hat{\psi}\hat{\psi}-\mu\hat{\psi}^{\dagger}\hat{\psi}\right],$
(3)

where the operators $\hat{\psi}(x)$ and $\hat{\psi}^{\dagger}(x)$ are the
Boson annihilation and creation operators satisfying the canonical commutation
relations

$\displaystyle[\hat{\psi}(x),\hat{\psi}^{\dagger}(x^{\prime})]=\delta(x-x^{\prime}),$
(4)
$\displaystyle[\hat{\psi}(x),\hat{\psi}(x)]=[\hat{\psi}^{\dagger}(x),\hat{\psi}^{\dagger}(x)]=0.$
(5)

The exact eigenstates and eigenenergies of this system are obtained through
the Bethe ansatz method Bethe (1931); Lieb and Liniger (1963); Lieb (1963).
The central object of this method is the Bethe ansatz equation

$\displaystyle e^{ik_{j}L}=\prod_{\ell\neq
j}^{N}\frac{k_{j}-k_{\ell}+ic}{k_{j}-k_{\ell}-ic}\quad\text{for}\quad
j=1,2,\ldots,N,$ (6)

the solution of which yields the exact wavefunction

$\displaystyle\varphi_{k_{1},\cdots,k_{N}}(x_{1},\cdots,x_{N})=\frac{c^{N/2}}{\sqrt{N!}}\left(\prod_{j>\ell}^{N}\frac{1}{k_{j}-k_{\ell}}\right)\sum_{\sigma\in
S_{N}}^{N!}A_{\sigma}\exp\left[i\sum_{j=1}^{N}k_{\sigma_{j}}x_{j}\right],$ (7)
$\displaystyle
A_{\sigma}=(-1)^{\sigma}\prod_{j>\ell}^{N}\left[k_{\sigma_{j}}-k_{\sigma_{\ell}}-ic\,\text{sign}(x_{j}-x_{\ell})\right]$
(8)

with the energy $E$ and momentum $P$ eigenvalues

$\displaystyle E=\sum_{j=1}^{N}k_{j}^{2},\qquad P=\sum_{j=1}^{N}k_{j}.$ (9)

Recasting the Bethe ansatz equation (6) into the logarithmic form, one has

$\displaystyle k_{j}L=2\pi I_{j}-2\sum_{\ell\neq
j}^{N}\arctan\left({\frac{k_{j}-k_{\ell}}{c}}\right)\quad\text{for}\quad
j=1,2,\ldots,N,$ (10)

where $I_{j}$’s are the Bethe quantum numbers, which are integers for odd $N$
and half-odd integers for even $N$. In the case of repulsive interaction
$c>0$, Bethe ansatz equation (10) has a unique real solution for a given set
of Bethe quantum numbers Korepin et al. (1997). Especially the ground state is
given by the set $I_{j}=-(N-1)/2+(j-1)$ for $j=1,2,\cdots,N$. Particle and
hole excitations are constructed by putting particles and holes in this Bethe
quantum numbers for the ground state.

In the attractive case $c<0$, the mean field equation has the bright soliton
solution and the corresponding quantum bound states are constructed in Wadachi
and Sakagami (1984). In this work we consider the repulsive case $c>0$, where
the mean field equation has the dark soliton solution Tsuzuki (1971); Zakharov
and Shabat (1973). In the refs. Sato et al. (2012, 2016), it is shown that the
quantum states corresponding to the dark soliton solutions are the
superpositions of the hole excitations

$\displaystyle\ket{X(q),N}:=\frac{1}{\sqrt{N}}\sum_{p=0}^{N-1}\exp(-2\pi
ipq/N)\ket{P(p),N}.$ (11)

Here $\ket{P(p),N}$ is the one-hole excitation with the total momentum
$P(p)=2\pi p/L$, whose quantum numbers are

$\displaystyle I_{j}$ $\displaystyle=-(N+1)/2+j\quad\mbox{for}\quad 1\leq
j\leq N-p$ $\displaystyle=-(N+1)/2+j+1\quad\mbox{for}\quad N-p+1\leq j\leq N.$
(12)

The density profile of this state
$\bra{X,N}\hat{\psi}^{\dagger}(x)\hat{\psi}(x)\ket{X,N}$ shows a density notch
at the position $X(q)=(q/N+1/2)L$, which coincides with the squared amplitude
of the classical dark soliton solution Sato et al. (2012, 2016). Here the
matrix element $\bra{X,N}\hat{\psi}^{\dagger}(x)\hat{\psi}(x)\ket{X,N}$ can be
effectively evaluated by the determinant formula for the norms of Bethe
eigenstates Gaudin (1983); Korepin (1982) and that of the form factors of the
density operator Slavnov (1989, 1990); Calabrese and Caux (2007)

$\displaystyle\bra{P^{\prime},N}\hat{\psi}^{\dagger}(x)\hat{\psi}(x)\ket{P,N}$
$\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)\frac{\det
U(k,k^{\prime})}{\sqrt{\det G(k)\det G(k^{\prime})}},$ (13)

where the quasimomenta $\\{k_{1},\cdots,k_{N}\\}$ and
$\\{k^{\prime}_{1},\cdots,k^{\prime}_{N}\\}$ give the eigenstates $\ket{P}$
and $\ket{P^{\prime}}$, 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.$ (14)

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}=2\delta_{j\ell}\text{Im}\left[\prod^{N}_{a=1}\frac{k^{\prime}_{a}-k_{j}+ic}{k_{a}-k_{j}+ic}\right]+\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).$
(15)

We can also consider the matrix element of the single field operator

$\displaystyle\begin{split}&\bra{X^{\prime},N-1}\hat{\psi}(x)\ket{X,N}\\\
&=\frac{1}{\sqrt{N(N-1)}}\sum^{N-1}_{p=0}\sum^{N-2}_{p^{\prime}=0}\exp\left[2\pi
i(p-p^{\prime})\frac{x}{L}\right]\exp\left[-2\pi
i\left(\frac{pq}{N}-\frac{p^{\prime}q^{\prime}}{N}\right)\right]\bra{P^{\prime},N-1}\hat{\psi}(0)\ket{P,N},\end{split}$
(16)

where $P=2\pi p/L$ and $P^{\prime}=2\pi p^{\prime}/L$ denote the total momenta
of the normalized Bethe eigenstates in the type II branch $\ket{P,N}$ and
$\ket{P^{\prime},N}$, respectively. The determinant formula is given by Gaudin
(1983); Korepin (1982); Slavnov (1989, 1990); Calabrese and Caux (2007);
Kojima et al. (1997)

$\displaystyle\begin{split}&\bra{P^{\prime},N-1}\hat{\psi}(0)\ket{P,N}=(-1)^{N(N+1)/2+1}\\\
&\times\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)\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})}},\end{split}$ (17)

where the quasi-momenta $\\{k_{1},\cdots,k_{N}\\}$ and
$\\{k^{\prime}_{1},\cdots,k^{\prime}_{N-1}\\}$ give the eigenstates
$\ket{P,N}$ and $\ket{P^{\prime},N-1}$, respectively. Here we remark that the
symbol $X^{{}^{\prime}}$ corresponds to an integer $q^{{}^{\prime}}$. We
recall that the matrix $G(k)$ denotes the Gaudin matrix, whose $(j,\ell)$th
element is given in (14). 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}\text{Im}\left[\frac{\prod^{N-1}_{a=1}(k^{\prime}_{a}-k_{j}+ic)}{\prod^{N}_{a=1}(k_{a}-k_{j}+ic)}\right]+\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).$
(18)

###  2.2 Construction with Gaussian weights

####  2.2.1 Effective Bose-Einstein condensation in the weak coupling case

In this subsection we briefly review the approach Kaminishi et al. (2020) to
quantum states with a dark soliton. We shall also investigate its validity
numerically. In the approach a corresponding quantum state is defined by
making use of a dark soliton solution $\phi_{P}(x)$ of the classical nonlinear
Schrödinger equation, and constructed explicitly with the Gaussian weights.

Based on the concept of the Bose-Einstein condensation (BEC) a “single-
particle dark soliton state” is proposed as an ideal quantum state of a dark
soliton in Ref. Kaminishi et al. (2020): All the particles occupy the same
one-particle state in the BEC. In the most strict sense, the BEC never happens
in one dimension. However, in a weak-interaction regime of a large but finite
system, a quasi-condensate state may appear. Thus, we assume a macroscopic
wave function for the quantum dark-soliton state. Here we suggest that in a
quasi-condensate state the coupling constant should be very small. It should
be inversely proportional to the system size in order to keep the condensation
fraction constant as the system size increases. The condensation fraction is
evaluated in the ground-state of the Lieb-Liniger model for the repulsive
coupling by applying the Onsager-Penrose criterion of the BEC Sato et al.
(2013).

The proposed BEC state Kaminishi et al. (2020), i.e., the quantum state
introduced with the wave function $\phi_{P}$, which we shall call the mean-
field product state, exactly reproduces the density and phase profiles of the
classical dark-soliton $\phi_{P}$. Furthermore, it was numerically studied
Kaminishi et al. (2020) that the main components of the product state are
given by a weighted sum of the yrast states. Here we remark that it is not
trivial to express it in terms of the eigenstates of the Lieb-Liniger model.

####  2.2.2 Ansatz with the ideal Gaussian weights and the mean-field product
state

Let us explain the approach Kaminishi et al. (2020). For a given macroscopic
wave function $\phi(x)$ and a given positive integer $N$ we consider a
corresponding BEC state consisting of $N$ particles. If all the particles
occupy the same single-particle state associated with the wave function
$\phi(x)$, the $N$-particle state is given by

$\frac{1}{\sqrt{N!}}\left(\int_{-L/2}^{L/2}dx\,\phi(x)\hat{\psi}^{\dagger}(x)\right)^{N}\ket{\Omega},$

where $\ket{\Omega}$ denotes the vacuum. Here, the number $N$ may be
macroscopically large. We call it the mean-field product state associated with
the wave function $\phi(x)$ in the paper. In Ref. Kaminishi et al. (2020) as
an $N$-particle quantum soliton state corresponding to a classical dark
soliton solution $\phi_{P}(x-X)$ the following BEC state is proposed

$\ket{N,X;P}=\frac{1}{\sqrt{N!}}\left(\int_{-L/2}^{L/2}dx\,\frac{\phi_{P}(x-X)}{\sqrt{N}}\hat{\psi}^{\dagger}(x)\right)^{N}\ket{\Omega}.$
(19)

This quantum state exactly reproduces the classical dark soliton density
profile which is given by

$\bra{N,X;P}\hat{\psi}^{\dagger}(x)\hat{\psi}(x)\ket{N,X;P}_{\mathrm{BEC}}=|\phi_{P}(x-X)|^{2},$
(20)

and satisfies

$\hat{\psi}(x)\ket{N,X;P}=\phi_{P}(x-X)\ket{N-1,X;P}.$ (21)

Here we remark that an explicit analytic expression of a dark soliton solution
$\phi_{P}(x)$ in a finite system will be given in section 3.2 Sato et al.
(2016).

It should be nontrivial to express the mean-field product state $\ket{N,X;P}$
in terms of the energy eigenstates of the Lieb-Liniger model. However, an
ansatz is given in Ref. Kaminishi et al. (2020) that it should be given by the
Gaussian weighted sum of yrast states. Let us consider an expansion of the
$N$-particle quantum soliton state $\ket{N,X;P}$ in terms of the eigenstates
$\ket{\\{I_{j}\\}_{N}}$ of the Lieb-Liniger model as follows:

$\ket{N,X;P}=\sum_{\\{I_{j}\\}}C_{\\{I_{j}\\}}\ket{\\{I_{j}\\}_{N}}.$ (22)

Quite an important numerical observation is made in Ref. Kaminishi et al.
(2020) that the expansion coefficients of the yrast states are dominant among
all the coefficients. It is shown in Figure 1. In other words, it is strongly
suggested that a quantum soliton state defined by Eq. (19) should be
essentially expressed as a superposition of the yrast states.

Figure 1: The overlap
$|C_{\\{I_{j}\\}}|^{2}=|\braket{\\{I_{j}\\}_{N}}{N,X;P=\pi}|^{2}$ for $n=1$,
$c=0.1$, and $N=8$. Each circle represents an eigenstate $\\{I_{j}\\}$ with
the momentum $\sum_{j=1}^{N}k_{j}$ and the energy $\sum_{j=1}^{N}k_{j}^{2}$,
and deep color implies a large overlap. The eigenstate with the lowest energy
for a given momentum $P$ corresponds to an yrast state
$\ket{N,P}_{\mathrm{yr}}$.

The expansion coefficients of the yrast states in the weak-coupling regime are
evaluated in Ref. Kaminishi et al. (2020) in the large $N$ limit. We now
explain the method how to evaluate them. In the thermodynamic limit with a
fixed coupling constant $c$, the mean $P$ and the variance $\sigma_{P}^{2}$ of
the total momentum operator
$\hat{P}=\int_{-L/2}^{L/2}dx\,\hat{\psi}^{\dagger}(x)(-i\partial_{x})\hat{\psi}(x)$
are given by

$\lim_{N\to\infty}\braket{N,X;P}{\hat{P}}{N,X;P}=P$ (23)

and

$\lim_{N\to\infty}\braket{N,X;P}{(\hat{P}-P)^{2}}{N,X;P}=\frac{4}{3}\beta^{3}n\sqrt{nc}\equiv\sigma_{P}^{2},$
(24)

respectively. Here $n$ denotes the density of particles: $n=N/L$ and $\beta$
will be introduced in section 3.3.1 as $\beta=1-v^{2}/v^{2}_{c,\infty}$. We
remark that $n$ and $\beta$ correspond to $\rho_{0}$ and $\gamma$ in Ref.
Kaminishi et al. (2020). Associated with Eqs. (23) and (24) we assume that in
the BEC state $\ket{N,X;P}$ the momenta $\\{p_{j}\\}_{j=1}^{N}$ of the $N$
particles can be regarded as independent random variables. Hence, the central
limit theorem implies that the total momentum $\sum_{j=1}^{N}p_{j}$ obeys a
Gaussian distribution. It thus follows from the above assumption that if the
quantum soliton state consists of only the yrast states, the former is given
by a Gaussian superposition of the latter for any large system size:

$\displaystyle\ket{N,X;P}$
$\displaystyle\approx\mathcal{N}^{-1/2}\sum_{P^{\prime}}e^{iP^{\prime}(X-L/2)}e^{-G_{P}(P^{\prime})/2}\ket{N,P^{\prime}}_{\mathrm{yr}}$
$\displaystyle\equiv\ket{N,X;P}_{\mathrm{yr}},$ (25)

where

$G_{P}(P^{\prime})=\frac{(P^{\prime}-P)^{2}}{2\sigma_{P}^{2}}$ (26)

and $\mathcal{N}$ is a normalization constant, and the sum is taken over
$P^{\prime}=2\pi M/L$ with $M=0,1,\dots,N$.

In Fig. 2 it is clear that the density profiles of the quantum and classical
solitons completely agree with each other. Here we remark that the density
profiles for the quantum soliton state of Eq. (25),
$\braket{N,0,P}{\hat{\psi}^{\dagger}(x)\hat{\psi}(x)}{N,0,P}$ with $P=\pi n$
(the black soliton) and $P=(\pi/2-1)n$ (a gray soliton) are plotted for
$N=L=100$ and $c=0.01$ together with the corresponding classical solutions of
the GP equation in Fig. 2. Furthermore, it is also clear that the expression
of Eq. (25) is much more useful for numerical calculation than Eq. (19). The
range of the summation in Eq. (25) is restricted only to the yrast states and
the expression of expansion coefficients is explicitly given. Thus, we can
calculate almost any physical observables by making use of the Bethe ansatz
method even for very large system sizes.

Figure 2: Comparison of the density profiles of quantum soliton states
$\ket{N,0;P}_{\mathrm{yr}}$ (the Gaussian superposition of yrast states) and
those of classical dark solitons, for $P=\pi n$ (left), $P=0.4502\pi n$
(middle) and $P=(\pi/2-1)n$ (right). The system size is set as $N=L=100$.

As far as the density profiles are concerned we might conclude that the
quantum soliton state should be completely constructed by the superposition of
the yrast states with the ideal Gaussian weights in Eq. (25) with Eq. (26). In
fact, the ansatz with the ideal Gaussian weights Kaminishi et al. (2020) has
the following two important properties: (i) We can construct a family of
quantum dark solitons by changing the value of $P$, where $P$ determines the
depth of the soliton. This is definitely an advantage of Eq. (25) compared to
the uniform superposition introduced in Ref. Sato et al. (2012) because the
uniform superposition of the yrast states cannot control the depth of the
soliton. Figure 2 shows that the Gaussian superposition indeed reproduces dark
solitons with an arbitrary depth (i.e., an arbitrary $P$), while the uniform
superposition does not; (ii) The Gaussian weights in Eqs. (25) and (26)
assigned to the yrast states are completely determined without any free
parameter. Here we recall that the Gaussian form is derived by applying the
central limit theorem to the quantum soliton state in Eq. (19) and the
analytic expression of $\sigma_{P}$ is obtained. Here it should be remarked
that the idea of a Gaussian superposition of yrast states has also been
proposed in a recent work Shamailov and Brand (2019), where the Gaussian form
is assumed and the Gaussian width $\sigma_{P}$ is given by a fitting
parameter.

####  2.2.3 Numerical study on the validity of the ansatz with the ideal
Gaussian weights

The ansatz with the ideal Gaussian weights Kaminishi et al. (2020) should thus
be quite useful and significant with respect to the density profile. However,
we now point out that it does not completely produce the mean-field product
state and also that Eq. (21) does not hold exactly but only approximately. In
Figure 3 the density profiles (red curves) are plotted together with the
corresponding profiles (blue curves) of the square amplitude of the matrix
element of the field operator between the quantum soliton states with $N$ and
$N-1$ particles. The two plots do not completely overlap: The former plots are
larger than the latter plots at any position of $x$. Here they are plotted for
$N=L=100$ and $c=0.1$. It is the case for all the examples investigated so
far. In the three panels of Figure 3 the density and square-amplitude profiles
are plotted for the quantum soliton states constructed by the ansatz of the
ideal Gaussian weights for three different values of soliton momentum $P$.
Here we remark that the soliton momentum $P$ and the variance $\sigma_{P}^{2}$
are determined by the target soliton depth $d$ through Eqs. (89) and (90). We
also remark that the definition of target soliton depth will be given in Eq.
(88).

For each value of the three values of depth $d$ the density profile of the
quantum soliton state constructed by the ansatz is clearly larger than the
corresponding square-amplitude profile. Here we remark that the degree of
overlap of the two profiles suggests the degree how far the quantum soliton
state constructed by the ansatz of the ideal Gaussian weights Kaminishi et al.
(2020) is consistent with the mean-field product state, as we shall argue in
section 4.3.3.

We suggest that the quantum states constructed by the ansatz of the ideal
Gaussian weights should be much closer to the mean-field product state than
those of the superposition of the yrast states with equal weight (11). In
Figure 3 we observe that in each value of target soliton depth $d$ which will
be defined by Eq. (88) the difference of the two profiles is smaller than that
of the profiles of the quantum states without the Gaussian weights.
Furthermore, we expect that one can improve the quantum soliton state
constructed by the ansatz so that we make it closer to an ideal mean-field
product state, by increasing the number of particles $N$ and the system size
$L$, while by decreasing the coupling constant $c$.

Figure 3: Density profiles (red curves) and the profiles of the square
amplitude of the matrix element of the field operator between the quantum
soliton states with $N$ and $N-1$ particles (blue curves) for $N=L=100$ and
the coupling constant $c=0.1$. The vertical axis represents the local density.
The profiles of the quantum soliton state without Gaussian weights are plotted
in the upper left panel. In other panels, the profiles for the quantum soliton
states with the ideal Gaussian weights are plotted for three values of target
soliton depth $d$ as $d=0,0.5$, and $0.6$, respectively. They correspond to
the parameters of the Gaussian weights,
$(P,\sigma^{2})=(\pi,0.421637),(0.18169\pi,0.149071)$, and
$(0.124027\pi,0.106667)$, respectively, through Eqs. (89) and (90).

##  3 Quantum state of a dark soliton with nonzero winding number and the
corresponding classical solution

We present the quantum state associated with an elliptic dark-soliton solution
of the GP equation with a nonzero winding number, in this section. Hereafter,
we call it the quantum state of a dark soliton with a nonzero winding number.
In order to define the winding number we briefly construct elliptic multiple
dark-solitons of GP equation and derive the single dark-soliton in the
thermodynamic limit.

###  3.1 Novel quantum state of a single dark soliton with a nonzero winding
number

####  3.1.1 Construction of quantum states with nonzero winding number

There exists a classical periodic solution of the NS equation
$\phi(x)=\sqrt{\rho(z)}\exp[i\varphi(x)]$ satisfying the periodic boundary
conditions

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

where $J$ is an arbitrary integer. It is called the phase-winding number
Kanamoto et al. (2009). We construct the quantum state of a single dark-
soliton with a nonzero winding number. When we construct the quantum dark-
soliton state as a superposition of Bethe eigenstates, we select a different
sequence of Bethe eigenstates, which have not been investigated in the
previous studies Sato et al. (2016); Kaminishi et al. (2020). Let us consider
a new quantum state in the following:

$\displaystyle\ket{X(q),N}_{J=1}:=\frac{1}{\sqrt{N}}\sum_{p=1}^{N}\exp(-2\pi
ipq/N)\ket{P(p),N}.$ (28)

The difference between the quantum states (11) and (28) is only one point that
the lower and upper bounds are shifted by 1 in the summation with respect to
$p$. That is, the summation of Eq. (28) is obtained from that of Eq. (11) by
removing the ground state ($p=0$) and adding the quantum state with momentum
$P=2\pi/L\ (p=N)$.

We express the matrix element of the field operator $\hat{\psi}(x)$ for a
quantum single dark soliton state $\ket{X,N}$ with a non-zero winding number
by

$\displaystyle\psi_{Q}(x):=[\braket{X,N-1}{\hat{\psi}(x)}{X,N}]_{J=1}$ (29)

and write $\psi_{Q}(x)$ in the polar form:
$\psi_{Q}(x)=|\psi_{Q}(x)|\exp[{i\arg[\psi_{Q}(x)]}]$. The matrix element is
computed by the determinant formula such as Eq. (17).

####  3.1.2 Profiles of the square amplitude and the phase of the matrix
element of the field operator

Figures 4 and 5 present the profiles of the square amplitude
$|\psi_{Q}(x)|^{2}$ and those of the phase field $\arg[\psi_{Q}(x)]$ in red
solid curves, respectively, for four different values of the coupling constant
such as $c=0.01,0.1,10$ and 100. The plots of the corresponding elliptic
solutions of the GP equation are exhibited in blue dashed curves for the four
values of the coupling constant $c$. Here the system size $L$ is given by
$L=20$ and the number $N$ of particles by $N=20$, so that the density $n$ is
given by 1: $n=N/L=1$.

When we search the corresponding classical solution of the GP equation to the
density profile of a given quantum dark-soliton state, we numerically evaluate
the soliton depth $d$ in the density profile, and then determine the value of
modulus $k$ through the relation: $a_{1}=d$, as will be shown in section
3.2.2. Here $a_{1}$ denotes the smallest root of the equation: $U(\rho)=0$,
where the potential $U(\rho)$ will be given by Eq. (36). We can draw the
density profile of the corresponding elliptic dark-soliton by applying Eq.
(46), as will be shown in section 3.2. Here we remark that an elliptic dark-
soliton has four parameters $c$, $L$, density $n$ and modulus $k$. We
calculate the chemical potential $\mu$ and the critical velocity $v_{c}$ by
Eqs. (58) and (59), respectively, as will be shown later. Moreover, when we
search the corresponding classical elliptic dark-soliton solution to the
matrix element of the field operator for a given quantum dark-soliton state,
we first evaluate the density $n$ by numerically integrating the whole profile
of the square amplitude of the matrix element. We secondlly evaluate the
soliton depth $d$ by finding out the minimum of the square amplitude profile.
We then determine the value of modulus $k$ through the relation: $a_{1}=d$.

Figure 4: Square amplitudes of the matrix elements $|\psi_{Q}(x)|^{2}$ for $N=L=20,$ and $c=0.01,\ 0.1,\ 10,\ \text{and}\ 100$ (red solid curves) and the corresponding elliptic solutions (blue dashed curves). The corresponding quantum state is given by Eq. (28), which has $J=1$. In the panel of $c=0.01$, the red solid curve agrees with the blue dashed curve. Figure 5: Phase profiles $\arg[\psi_{g}(x)]/\pi$ for $N=L=20,$ and $c=0.05,\ 0.1,\ 10,\ \text{and}\ 100$ (red solid curves) and the corresponding elliptic solutions (blue dashed curves). The corresponding quantum state is given by Eq. (28), which has $J=1$. In the panels of $c=0.01$ and 0.1, the red solid curves completely agree with the blue dashed curves, respectively. Table 1: Parameters of the profiles of the elliptic solutions in Figures 4 and 5. The notch depth $d$ is given by the minimum of the square amplitude profile for each value of $c$.  $c$ | $n$ | $d$ | $k$ | $\mu/2nc$ | $v_{c}/2\sqrt{nc}$ | $v/v_{c}$
---|---|---|---|---|---|---
0.01 | 0.830 840 | 9.670 81 $\times 10^{-7}$ | 0.680 75 | 1.446 70 | 1.753 96 | 0.982 297
0.1 | 0.903 971 | 0.00172 103 | 0.993 356 | 1.272 54 | 0.853 159 | 0.575 645
1.0 | 0.788 041 | 0.00192 928 | 1- 5.693 69$\times 10^{-8}$ | 1.111 45 | 0.899 118 | 0.147 633
10 | 0.461 185 | 0.00419 681 | 1- 7.864 09$\times 10^{-19}$ | 1.05 184 | 0.949 564 | -0.0188 539
100 | 0.300 536 | 0.00640 668 | 1- 2.2907 21$\times 10^{-47}$ | 1.022 293 | 0.977 873 | -0.117 299

In Figure 4, the plots of the square amplitudes $|\psi_{Q}(x)|^{2}$ agree with
those of the corresponding elliptic solution represented with the blue dashed
lines, in particular, for the weak coupling cases of $c=0.01$ and 0.1. In the
cases of $c=0.01,0.1$ and 10, the width of the notch decreases in the square
amplitude profile as the coupling constant increases. In the case of $c=100$,
however, the notch of the square amplitude $|\psi_{Q}(x)|^{2}$ is wider than
that of the elliptic solution, as shown in Figure 4.

In Figure 5 all the phase profiles $\arg[\psi_{Q}(x)]$ exhibited satisfy the
boundary condition

$\displaystyle\arg[\psi_{Q}(x+L)]=\arg[\psi_{Q}(x)]+2\pi\,.$ (30)

It corresponds to the case of $J=1$ in Eq. (27). Moreover, the phase profiles
of the quantum states almost perfectly agree with those of the classical
solutions. All the blue dotted curves for the phase profiles of the classical
solutions with $J=1$ are depicted by choosing the plus in the plus/minus sign
in the constant $W$, which will be defined in Eq. (34). The choice should be
consistent with the abruptly increasing behavior appearing at each notch
position, as we shall argue in section 3.3.2.

The parameters for the fitted curves in Figures 4 and 5 are listed in Table 1.
The estimates of the velocity parameter $v$ are given by making use of Eq. (52
) with $J=1$. For the strong coupling cases of $c=10$ and $100$ the velocity
parameter $v$ becomes negative, which will be explained with Eqs. (74) and
(75) in section 3.3.4.

Figure 6: Density profiles for $N=L=20,$ and $c=0.01,\ 0.1,\ 10,\ \text{and}\ 100$ (red solid curves), and the corresponding elliptic solutions (blue dashed curves). The corresponding quantum state is given by Eq. (28), which has $J=1$. In the panels of $c=0.01$ and 0.1, the red solid curves completely agree with the blue dashed curves, respectively. Table 2: Parameters of the density profiles of the elliptic dark-soliton solutions in Figure 6.  $c$ | $d$ | $k$ | $\mu/2nc$ | $v_{c}/2\sqrt{nc}$ | $v/v_{c}$
---|---|---|---|---|---
0.01 | 0.170 659 | 0.680 343 | 1.690 17 | 1.444 67 | 0.978 349
0.1 | 0.155 884 | 0.991 939 | 1.40 035 | 0.667 869 | 0.317 074
10 | 0.0406 314 | 1- 3.64 668$\times 10^{-17}$ | 1.04 120 | 0.959 205 | -0.151 773
100 | 0.0475 028 | 1- 4.973 27$\times 10^{-85}$ | 1.013 54 | 0.986 482 | -0.202 798

In Figure 6 the density profiles
$\left[\braket{X(q),N}{\hat{\psi}^{\dagger}(x)\hat{\psi}(x)}{X(q),N}\right]_{J=1}$
of the novel quantum state (28) are depicted for the four values of the
coupling constant $c=0.01,0.1,10$ and 100 with red solid curves. The blue
dotted curves are the density profiles of the corresponding classical
solutions, i.e., elliptic single dark-solitons. Here we recall that they are
depicted by evaluating only the soliton depth $d$.

It is interesting to observe that the density profiles of the quantum state
with $J=1$ in Figure 6 are similar to those of $J=0$ depicted in earlier
studies. It is consistent with the fact that the local density $\rho(x)$ of
the classical solution does not depend on the winding number $J$.

The parameters for the fitted curves in Figure 6 are listed in Table 2. The
estimates of the velocity parameter $v$ are given by making use of Eq. (52)
with $J=1$. Similarly as in Table 1, for the strong coupling cases of $c=10$
and $100$ the velocity parameter $v$ are negative. It will be explained with
Eqs. (74) and (75) in section 3.3.4.

Finally we remark that quantum states of a single dark soliton with zero
winding number, i.e., $J=0$, have been studied in the previous study Sato et
al. (2016). In all the phase profiles the argument $\arg[\psi_{Q}(x)]$
satisfies the periodic boundary condition:
$\arg[\psi_{Q}(x)]=\arg[\psi_{Q}(x+L)]$.

###  3.2 Construction of an elliptic multiple dark soliton for the GP equation

####  3.2.1 Local density $\rho(x)$ and the phase field $\varphi(x)$ of an
elliptic multiple dark-soliton

Let us consider stationary solutions of the GP equation in the region of
length $L$ under the periodic boundary conditions. We assume that they depend
on the coordinate $x$ and time $t$ through variable $z=x-vt$ as follows:

$\psi(x,t)=\phi(x-vt).$ (31)

Substituting it into Eq. (1) we have

$\phi^{{}^{\prime\prime}}-iv\phi^{{}^{\prime}}+\mu\phi-2c|\phi|^{2}\phi=0.$
(32)

If there are $N$ particles in the system of length $L$, the average density
$n$ is given by $n=N/L$. Here we remark that in the limit of sending $L$ to
infinity, the chemical potential $\mu$ approaches $2cn$ if the local density
$|\phi(x)|^{2}$ approaches $n$ for some large region. We express it as
$\mu_{\infty}=2cn$.

We express the complex scalar field $\phi(z)$ in terms of the local density
$\rho(z)$ and the phase field $\varphi(z)$ as follows.

$\phi(z)=\sqrt{\rho(z)}\exp\left(i\varphi(z)\right).$ (33)

By putting the ansatz (33) into Eq. (32), from the real part and the imaginary
part of the derived equation we can show that there are two first-integrals of
motion:

$\frac{d\varphi}{dx}=\frac{v}{2}+\frac{W}{\rho(x)},$ (34)

and

$\left(\frac{1}{2}\frac{d\rho}{dx}\right)^{2}+U(\rho)=0,$ (35)

where

$U(\rho)=-c\rho^{3}+\left(\mu+{\frac{v^{2}}{4}}\right)\rho^{2}-2V\rho+W^{2}.$
(36)

Here $W$ and $V$ are constants with respect to the first integrals (34) and
(35), respectively.

Let $a_{1}$, $a_{2}$ and $a_{3}$ denote the three roots of the equation:
$U(\rho)=0$. We assume that they are real. We put them in increasing order:
$a_{1}\leq a_{2}\leq a_{3}$. It follows from the first integral (35) that they
are non-negative. The coefficient of $\rho^{3}$ in the potential $U(\rho)$ of
Eq. (36) is negative since $c$ is positive: $c>0$, and hence the potential
$U(\rho)$ is negative or zero only in the regions $a_{1}\leq\rho\leq a_{2}$
and $a_{3}\leq\rho$. Physically the local density $\rho$ should be non-
negative and have an upper bound. It thus follows that we consider only the
region $a_{1}\leq\rho\leq a_{2}$ and $a_{1}$ is non-negative: $a_{1}\geq 0$.
In fact, the potential $U(\rho)$ is non-negative at the origin:
$U(\rho=0)=W^{2}\geq 0$. It is straightforward to show the following:

$W=\pm\sqrt{ca_{1}a_{2}a_{3}}.$ (37)

By solving Eq. (35) we express the local density $\rho(z)$ in terms of the
Jacobi’s elliptic function ${\rm sn}(u,k)$. Integrating the differential
$2dx=d\rho/\sqrt{-U(\rho)}$ with the initial condition: $\rho(x=0)=a_{1}$ and
applying the variable transformation: $\rho=a_{1}+(a_{2}-a_{1})z^{2}$, we have

$\rho(x)=a_{1}-(a_{1}-a_{2}){\rm
sn}^{2}\left(\sqrt{c(a_{3}-a_{1})}\,x,k\right)\,,$ (38)

where $k$ denotes the elliptic modulus

$k=\sqrt{\frac{a_{3}-a_{1}}{a_{2}-a_{1}}}.$ (39)

We now assume the period $L/m$ with a positive integer $m$ for stationary
solutions. we thus have

$\sqrt{c(a_{3}-a_{1})}={2mK}/L,$ (40)

where $K$ denotes the complete elliptic integrals of the first kind.

We now express the three roots $a_{1},a_{2}$ and $a_{3}$ in terms of the
average density $n$, $c$, $L$ and $k$. Since the average of the local density
over the period $L$ is equal to the average density $n$, we have

$\frac{1}{L}\int_{0}^{L}\rho(x)dx=n.$ (41)

Through Eq. (39) and the relation: ${\rm sn}^{2}(u,k)+k^{2}{\rm
dn}^{2}(u,k)=1$, the local density $\rho(z)$ is expressed in terms of the
Jacobi’s elliptic function ${\rm dn}(u,k)$ as

$\rho(x)=a_{3}-(a_{3}-a_{1}){\rm dn}^{2}(2mKx/L,\,k).$ (42)

It follows from the definition of the complete elliptic integral of the second
kind: $E=\int_{0}^{K}{\rm dn}^{2}(u,\,k)du$ and Eq. (41) that we have

$a_{3}=n+\frac{1}{c}\left(\frac{2mK}{L}\right)^{2}\,\frac{E}{K}.$ (43)

Let us introduce a parameter $\beta_{k}$ by

$\beta_{k}=\frac{1}{nc}\left(\frac{2mK}{L}\right)^{2}\,.$ (44)

It follows from Eqs. (39) and (41) that the three roots are expressed as

$\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle
n\left(1+\beta_{k}\left(\frac{E}{K}-1\right)\right)\,,$ $\displaystyle a_{2}$
$\displaystyle=$ $\displaystyle
n\left(1+\beta_{k}\left(\frac{E}{K}-1+k^{2}\right)\right)\,,$ $\displaystyle
a_{3}$ $\displaystyle=$ $\displaystyle
n\left(1+\beta_{k}\frac{E}{K}\right)\,.$ (45)

Thus, the three roots $a_{1},a_{2},a_{3}$ are expressed in terms of
$\beta_{k}$, the density $n$, and the modulus $k$. The parameter $\beta_{k}$
depends on $c$, $n$, $k$ and $L/m$. Here the integer $m$ denotes the number of
notches in the multiple dark-soliton, i.e., the multiplicity. The parameter
$\beta_{k}$ has to satisfy the condition: $\beta_{k}\leq 1/(1-E/K)$ since
$a_{1}$ is non-negative: $a_{1}\geq 0$.

The local density is expressed in the form of the average density $n$ and the
oscillating part

$\rho(x)=n+n\,\beta_{k}\,\left(\frac{E}{K}-{\rm
dn}^{2}\left(\frac{2mKx}{L},k\right)\right)\,.$ (46)

Here we remark that the second term in the right-hand side of Eq. (46) leads
to a logarithmic correction associated with the conservation of the number of
particles when we send $k$ to 1 and $L$ to infinity with $K/L$ fixed.

Let us now express the phase field in terms of elliptic functions. By
integrating Eq. (34) with respect to variable $x$ we have

$\varphi(x)-\varphi(0)={\frac{v}{2}}x+{\frac{W}{a_{1}}}\int_{0}^{x}\frac{dy}{1-n_{e}\,{\rm
sn}^{2}\left({2mKy}/L,k\right)}$ (47)

where $n_{e}$ is given by

$n_{e}=1-a_{2}/a_{1}\,.$ (48)

In terms of the elliptic integral of the third kind $\Pi(r,\phi,k)$ with
$\sin\phi={\rm sn}(u,k)$ and a given number $r$

$\Pi(r,\phi={\rm am}\,u,k)=\int_{0}^{u}\frac{dv}{1-r\,{\rm sn}^{2}(v,k)}$ (49)

we have

$\varphi(x)-\varphi(0)={\frac{v}{2}}x+{\frac{W}{a_{1}}}{\frac{L}{2m}}{\frac{1}{K}}\Pi\left(n_{e},\phi={\rm
am}\,\frac{2mKx}{L},k\right).$ (50)

In order to determine the velocity $v$, let us consider the shifted periodic
boundary condition (27) with a winding number $J$.

$\varphi(x=\frac{L}{2m})-\varphi(x=0)={\frac{2\pi J}{2m}}.$ (51)

Putting Eq. (50) into it, we express the velocity $v$ in terms of elliptic
functions as follows:

$v=-\frac{2W}{a_{1}}\frac{1}{K}\Pi(n_{e},k)+{\frac{4\pi}{L}}J.$ (52)

We now show that the phase field $\varphi(x)$ is expressed in terms of
Jacobi’s Theta function. Let $\Theta(u)$ and $Z(u)$ denote Jacobi’s Theta and
Zeta functions, respectively. Making use of the connection of the Theta
functions to the elliptic integral of the third kind Copson (1978):

$\int_{0}^{u}\frac{k^{2}{\rm sn}(a,k){\rm cn}(a,k){\rm dn}(a,k){\rm
sn}^{2}(u,k)}{1-k^{2}{\rm sn}^{2}(a,k){\rm
sn}^{2}(u,k)}\,du=\frac{1}{2}\log\left(\frac{\Theta(u-a)}{\Theta(u+a)}\right)+uZ(a)$
(53)

we derive the following formula

$\displaystyle\int_{0}^{u}\frac{dv}{1-r{\rm sn}^{2}(v,k)}$ $\displaystyle=$
$\displaystyle u\left(1+Z(a)\sqrt{\frac{r}{(k^{2}-r)(1-r)}}\right)+$ (54)
$\displaystyle\qquad+\frac{1}{2}\log\left(\frac{\Theta(u-a)}{\Theta(u+a)}\right)\sqrt{\frac{r}{(k^{2}-r)(1-r)}}$

where $r$ is any real number while the parameter $a$ is related to $r$ by
$k^{2}{\rm sn}^{2}(a,k)=r$. In order to apply Eq. (54) to the phase field, we
set $r=n_{e}=1-a_{2}/a_{1}$. We define the parameter $a$ by a pure imaginary
number $a=i\alpha$ with $\alpha>0$ satisfying the following relation:

${\rm sn}(i\alpha,k)=i\sqrt{\frac{\beta_{k}}{1-\beta_{k}(1-E/K)}}\,.$ (55)

Then, we can numerically evaluate the real number $\alpha$ by making use of
the relations: ${\rm sn}u=\Theta(K)H(u)/H(K)\Theta(u)$ and
$\sqrt{k}=H(K)/\Theta(K)$.

Thus, making use of the relation between the elliptic integral of the third
kind and Jacobi’s Zeta and Theta functions, we express the phase field
$\varphi(x)$ as the logarithm of a ratio of Theta functions

$\varphi(x)-\varphi(0)={\frac{2\pi
J}{L}}x\pm{\frac{i}{2}}\log\left(\frac{\Theta({2mKx}/L-i\alpha)}{\Theta({2mKx}/L+i\alpha)}\right)\,,$
(56)

for $W=\pm\sqrt{ca_{1}a_{2}a_{3}}$. Here the plus/minus signs are in
correspondence. In the derivation of Eq. (56), through Eq. (54) we have
expressed the complete elliptic integral of the third kind in terms of the
Zeta function

$\Pi(r,k)=K\left(1+Z(a,k)\sqrt{\frac{r}{(k^{2}-r)(1-r)}}\right).$ (57)

####  3.2.2 Chemical potential and the critical velocity

When the system size $L$, the coupling constant $c$, and the density $n$ are
fixed, once the modulus $k$ is specified, the corresponding classical solution
of the multiple dark-soliton with multiplicity $m$ is determined. Here, the
chemical potential $\mu$ is related to the parameters $c$, $n$, $L$ and $k$ as
follows:

$\mu=3nc-{\frac{nc}{4}}v^{2}+{\frac{4}{\left(L/m\right)^{2}}}\,K(3E-(2-k^{2})K).$
(58)

Here we remark that velocity $v$ is given by Eq. (52) with a winding number
$J$.

For a given quantum dark-soliton state, we search the corresponding elliptic
solution of the GP equation, as follows. Suppose that the system size $L$, the
coupling constant $c$, and the density $n$ are also given. We then determine
the modulus $k$ by equating the first root $a_{1}$ with the depth $d$ at the
notch, i.e., the minimum value of the density profile. Furthermore, we
evaluate the chemical potential $\mu$ by Eq. (58). Here we recall that
velocity $v$ is evaluated by Eq. (52).

Let us derive the range of velocity $v$. When the system size $L$, the
coupling constant $c$, and the density $n$ are fixed, and the chemical
potential $\mu$ is also fixed, then the value of velocity $v$ specified by Eq.
(58) is a function of modulus $k$. We thus have the range of velocity $v$ as
we change the value of modulus $k$ from 0 to 1. We can show that the function
$f(k)=K(3E-(2-k^{2})K)$ is monotonically decreasing with respect to modulus
$k$. The largest value of $f(k)$ is therefore given by $f(0)=\pi^{2}/4$. It
follows that the largest value $v_{c}$ of the velocity in the range is given
by

$v_{c}=\sqrt{(2m\pi/L)^{2}+4(3nc-\mu)}.$ (59)

We call it the critical velocity.

###  3.3 Reduction of elliptic multiple dark-soliton

####  3.3.1 Reduction of an elliptic multiple dark soliton into a series of
single dark-solitons

We shall show that an elliptic multiple dark-soliton reduces into separate
single dark-solitons if we increase the system size in a certain way. Let us
consider the limit of sending the system size $L$ to infinity ($L\to\infty$)
while the modulus $k$ to 1 ($k\to 1$) simultaneously so that the ratio
$K(k)/L$ is kept constant. We define the ratio $b$ by $b=K/L$. Here we remark
that we keep the density $n$ constant. We define $\beta$ by the simultaneous
limit of sending $L$ to infinity and $k$ to 1 with $K/L=b$:

$\beta=\lim_{k\to 1,L\to\infty}\beta_{k}=\frac{(2mb)^{2}}{cn}\,.$ (60)

We introduce the complementary modulus $k^{{}^{\prime}}$ by
$k^{{}^{\prime}}=\sqrt{1-k^{2}}$. We define $K^{{}^{\prime}}(k)$ and
$E^{{}^{\prime}}(k)$ by $K^{{}^{\prime}}(k)=K(k^{{}^{\prime}})$ and
$E^{{}^{\prime}}(k)=E(k^{{}^{\prime}})$, respectively. We remark that $K$ and
$iK^{{}^{\prime}}$ are the quarter-periods of Jacobi’s elliptic functions. We
define an expansion parameter $p$ by

$p=\exp(-\pi K/K^{{}^{\prime}}),$ (61)

which is small when modulus $k$ is close to 1.

We now consider the density profile of the stationary solution of an elliptic
$m$-fold multiple dark-soliton at time $t=0$. We assume that at time $t=0$ the
position of the $j$-th notch $x_{j}$ is given by $x_{j}=jL/m$ for
$j=0,1,\ldots,m-1$, in the interval $[0,L]$. Near the position of the $j$th
notch, we express the $x$-coordinate in terms of the position of the $j$-th
notch and the deviation from it:

$x=\Delta x+x_{j}.$ (62)

Here we assume that $|\Delta x|$ is small: $|\Delta x|\ll L/2m$. Near the
$j$-th notch the square amplitude $\rho(x)$ is expanded as an asymptotic
series of parameter $p$ as follows:

$\rho(x)/n=1-\beta\,\text{sech}^{2}(2b\,\Delta x)+2\beta/\log(1/p)+\cdots\,.$
(63)

It is easy to show that the asymptotic expansion of the square amplitude up to
the order of $1/\log p$, i.e., up to the order of $1/L$, satisfies the
normalization condition. Thus, near the $j$-th notch, the density profile of
the elliptic multiple dark-soliton reduces to a single dark-soliton in the
limit of sending $L$ to infinity and the elliptic modulus $k$ to 1 with the
ratio $K/L$ fixed as $K/L=b$:

$\rho(x)=\left(\frac{\mu_{\infty}}{2c}\right)\left\\{1-\beta\,{\rm
sech}^{2}\left[\left(\frac{\beta\mu_{\infty}}{2}\right)^{1/2}\,(x-x_{j})\right]\right\\}.$
(64)

Here $\beta$ is expressed as

$\beta=1-\left(\frac{v}{v_{c,\infty}}\right)^{2}$ (65)

where $v_{c,\infty}=\sqrt{2\mu_{\infty}}$. Here we remark that we have derived
the square amplitude of the single dark-soliton due to Tsuzuki Tsuzuki (1971)
in the limit of sending $p$ to zero. Hereafter we express $\mu_{\infty}$
simply as $\mu$ for convenience.

We now consider the phase profile near a notch of an elliptic multiple dark-
soliton. We recall that near the $j$th notch at time $t=0$ we have expressed
the $x$-coordinate as $x=\Delta x+x_{j}$ where $x_{j}=jL/m$. In the limit of
sending $L$ to infinity and $k$ to 1 with $K/L=b$ fixed the phase field is
expressed in terms of the logarithmic functions as follows:

$\displaystyle\varphi(x)-\varphi(0)$ $\displaystyle=$
$\displaystyle{\frac{2\pi
J}{L}}\,x\,\mp\frac{1}{2i}\log\left(1-i\sqrt{\frac{\beta}{1-\beta}}\tanh\left(\sqrt{\frac{\mu\beta}{2}}\,\Delta
x\right)\right)$ (66)
$\displaystyle\qquad\pm\frac{1}{2i}\log\left(1+i\sqrt{\frac{\beta}{1-\beta}}\tanh\left(\sqrt{\frac{\mu}{2}\beta}\,\Delta
x\right)\right)\quad+O(1/\log(1/p))$

for $W=\pm\sqrt{ca_{1}a_{2}a_{3}}$. Here the plus/minus signs are in
correspondence.

For every $j$-th notch we have the expression (64) of the local density as a
function of the deviation $\Delta x$ from the $j$-th notch at time $t=0$.
Thus, in the limit of sending $L$ to infinity and $k$ to 1 where ratio $K/L$
is fixed as $K/L=b$ an elliptic multiple dark-soliton reduces into a series of
separated single dark-solitons.

####  3.3.2 Abrupt change of the phase profile near a notch

We express the position $x$ as $x=\Delta x+x_{j}$ near the $j$th notch at time
$t=0$. It follows from Eq. (66) that for very small $|\Delta x|$ the phase
field $\varphi(x)$ is approximately given by the following if the system size
$L$ is very large

$\varphi(x)-\varphi(0)={\frac{2\pi J}{L}}\left(\Delta
x+x_{j}\right)\pm\sqrt{\frac{\beta}{1-\beta}}\,\sqrt{\frac{\mu\beta}{2}}\,\Delta
x+o(|\Delta x|)\quad(|\Delta x|\ll 1).$ (67)

Thus, if $L$ is very large, the phase change near the $j$-th notch is
approximately given by

$\varphi(x)-\varphi(0)-{\frac{2\pi
J}{L}}\,x_{j}\approx\left(\pm\sqrt{\frac{\mu(1-\beta)}{2}}\cdot\frac{\beta}{1-\beta}+{\frac{2\pi
J}{L}}\right)\,\Delta x\quad(|\Delta x|\ll 1)$ (68)

for $W=\pm\sqrt{ca_{1}a_{2}a_{3}}$. Here we recall that the plus/minus signs
are in correspondence.

The gradient of the phase near the $j$-th notch is therefore approximated by

$\left.\frac{d\varphi}{dx}\right|_{\Delta
x=0}=\pm\sqrt{\frac{\mu(1-\beta)}{2}}{\frac{\beta}{1-\beta}}+{\frac{2\pi
J}{L}}$ (69)

for $W=\pm\sqrt{ca_{1}a_{2}a_{3}}$. Thus, the phase field $\varphi(z)$ may
abruptly increase ($W>0$) or decrease ($W<0$) near the $j$th notch
$x_{j}=jL/m$ if $|J|$ is not large, in particular, when $\beta$ is close to 1,
since $1/\sqrt{1-\beta}$ can be much larger than 1.

When the winding number is zero: $J=0$, the velocity parameter $v$ is positive
if $W$ is negative from Eq. (52) or as we shall show in Eq. (70). The phase
field abruptly decreases at the $j$th notch, i.e., we have the minus sign in
Eq. (68).

The abrupt change of the phase at the notches of multiple dark-soliton will
appear in the phase profile of quantum states corresponding to double dark-
solitons in §4.

####  3.3.3 Evaluation of velocity in an infinite limit

We first recall that through Eq. (57) the velocity is expressed in terms of
the Zeta function:

$\frac{v}{2}=-\frac{W}{a_{1}}\,\left(1+Z(a,k)\sqrt{\frac{n_{e}}{(k^{2}-{n_{e}})(1-{n_{e}})}}\right)+{\frac{2\pi
J}{L}}\,$ (70)

where $W$ is given by $W=\pm\sqrt{ca_{1}a_{2}a_{3}}$.

Let us evaluate the velocity $v$ in the limit of sending $L$ to infinity and
$k$ to 1 where $K/L$ is kept constant as $K/L=b$. Through the Jacobi’s
imaginary transformation we can show

$Z(i\alpha,k)=i\sqrt{\frac{\beta}{1-\beta}}+O(1/\log(1/p))\quad(p\rightarrow
0).$ (71)

Putting it into Eq. (70) we have

$\frac{v}{2}\approx\mp\sqrt{\frac{\mu(1-\beta)}{2}}+{\frac{2\pi
J}{L}}\qquad(p\ll 1)$ (72)

for $W=\pm\sqrt{ca_{1}a_{2}a_{3}}$.

We recall that for $J=0$ the absolute value of the phase gradient at a notch
is given by

$\left.\frac{d\varphi}{dx}\right|_{\Delta
x=0}={\frac{v}{2}}\cdot{\frac{\beta}{1-\beta}}.$ (73)

It is larger than the value of velocity $v$ if $\beta$ is close to 1. Here we
remark that the factor $\beta/(1-\beta)$ is larger than 1 if $\beta>1/2$.

####  3.3.4 Velocity for a nonzero winding number

Let us consider the velocity $v$ in the case of $J=1$. For
$W=\pm\sqrt{ca_{1}a_{2}a_{3}}$ we have the following approximation for $v$ if
the system size $L$ is very large

$\frac{v}{2}=\mp\sqrt{\frac{\mu(1-\beta)}{2}}+\frac{2\pi}{L}.$ (74)

When $W>0$, the condition that $v$ is positive is given by the following:

$cn(1-\beta)<\left(\frac{2\pi}{L}\right)^{2}.$ (75)

It holds if the coupling constant $c$ is small or $\beta$ is close to 1.

If the system size $L$ is very large, the phase change near the $j$th notch of
an elliptic $m$-fold multiple dark-soliton is approximately given by

$\varphi(x)-\varphi(0)-\frac{2\pi}{L}x_{j}\approx\left(\pm\sqrt{\frac{\mu(1-\beta)}{2}}\,\frac{\beta}{1-\beta}\,+{\frac{2\pi}{L}}\right)\,\Delta
x\qquad(|\Delta x|\ll 1)$ (76)

for $W=\pm\sqrt{ca_{1}a_{2}a_{3}}$.

It follows from Eqs. (74) and (76) that if $\beta$ is close to 1 it happens
very rarely that both the phase gradient $d\varphi/dx$ and the velocity $v$
are positive simultaneously. It occurs only if the coupling constant $c$ is
very small.

For instance, if $W$ is negative, the sign in Eq. (74) is positive and hence
$v$ is positive, while that of Eq. (76) is negative. If $\beta$ is very close
to 1 then the factor $1/(1-\beta)$ is large, so that the phase gradient is
negative unless the coupling constant $c$ is extremely small such as it is
proportional to $1/L^{2}$.

On the other hand, if $W$ is positive, the sign in Eq. (74) is negative and
hence $v$ is positive only if the inequality in Eq. (75) holds, while that of
Eq. (76) is positive and hence the phase gradient is positive. Thus, for
$W>0$, both the phase gradient $d\varphi/dx$ and the velocity $v$ are positive
simultaneously if the inequality (75) holds.

##  4 Quantum double dark soliton

###  4.1 Density profile in the free fermions

In order to explore a quantum state in the Lieb-Liniger model whose density
profile has two notches, we now consider the one-dimensional free fermionic
system, which corresponds to the Lieb-Liniger model in the limit of sending
the coupling constant to infinity Girardeau (1960). First, we construct the
quantum state in the free fermions which possesses one notch in its density
profile shown in Figure 7 (a). Then, we derive a quantum state whose density
profile has two notches.

Let us introduce the second quantized free fermions under the periodic
boundary conditions with period $L$. The creation (annihilation) operator
$a_{\ell}^{\dagger}$ ($a_{\ell}$) creates (annihilates) a particle with a
momentum $k_{\ell}=2\pi n_{\ell}/L$, where $n_{\ell}$ is an integer. The
operators $a_{\ell}$ and $a_{\ell}^{\dagger}$ satisfy the anti-commutation
relations

$\displaystyle\\{a_{\ell},a^{\dagger}_{m}\\}=\delta_{\ell,m},$
$\displaystyle\\{a_{\ell},a_{m}\\}=\\{a^{\dagger}_{\ell},a^{\dagger}_{m}\\}=0.$
(77)

The field operators for the one dimensional free fermions $\psi(x)$ and
$\psi^{\dagger}(x)$ are defined by the creation and annihilation operators
respectively

$\displaystyle\psi(x)=\frac{1}{\sqrt{L}}\sum_{\alpha=1}^{\infty}\mathrm{e}^{\mathrm{i}k_{\alpha}x}a_{\alpha},$
$\displaystyle\psi^{\dagger}(x)=\frac{1}{\sqrt{L}}\sum_{\alpha=1}^{\infty}\mathrm{e}^{-\mathrm{i}k_{\alpha}x}a_{\alpha}^{\dagger}.$
(78)

We consider a quantum state possessing one notch in its density profile

$\displaystyle\ket{\Phi_{1}}$
$\displaystyle=\psi(X_{1})\ket{gs}=\frac{1}{\sqrt{\mathcal{N}_{1}}}\sum_{\alpha=1}^{N}(-1)^{\alpha}e^{ik_{\alpha}X_{1}}\prod_{\ell=1;\ell\neq\alpha}^{N}a_{\ell}^{\dagger}\ket{0}.$
(79)

It is obtained by applying a field operator $\psi(X_{1})$ at the position
$x=X_{1}$ to the $N$-particle ground state $\ket{gs}$ defined by

$\displaystyle\ket{gs}\equiv\prod_{\ell=1}^{N}a_{\ell}^{\dagger}\ket{0}.$ (80)

Here, $\mathcal{N}_{1}$ in Eq. (79) is the normalization factor. The state
$\ket{\Phi_{1}}$ is the superposition of such states obtained by removing a
particle with momentum $k=2\pi n_{\alpha}/L$ from the $N$-particle ground
state where the index $\alpha$ runs from 1 to $N$. The phase factor
$e^{-ik_{\alpha}X_{1}}$ in Eq. (79) has a position $X_{1}$, implying the
position of the notch in the expectation value of the local density operator
for its state Sato et al. (2016).

Generalizing the state $\ket{\Phi_{1}}$, we assume that the fermionic state
having two notches is a state such that the ground state with two field
operators at $x=X_{1},\ X_{2}$, that is,

$\displaystyle\ket{\Phi_{2}}=\psi(X_{2})\psi(X_{1})\ket{gs}=\frac{1}{\sqrt{\mathcal{N}_{2}}}\sum_{\alpha<\beta}^{N}(-1)^{\alpha+\beta}\left(e^{i\left(k_{\alpha}X_{1}+k_{\beta}X_{2}\right)}-e^{i\left(k_{\beta}X_{1}+k_{\alpha}X_{2}\right)}\right)\times\prod_{\ell=1;\
\ell\neq\alpha,\beta}^{N}a_{\ell}^{\dagger}|0\rangle.$ (81)

The state $\ket{\Phi_{2}}$ is constructed by the double summation of the
$N$-particle states which are given by removing two particles momenta
$k_{\alpha}=2\pi n_{\alpha}/L$ and $k_{\beta}=2\pi n_{\beta}/L$ in the ground
state. Here, the sum is taken over $\alpha$ and $\beta$. Figure 7 (b) shows
the density profile for the state (81). It has two notches at the positions
$x=X_{1}$ and $X_{2}$.

Figure 7: Density profiles in the free fermions. (a) Density profile of the
state of Eq. (79) with $N=L=21,\ $ and $X_{1}=0$. (b) Density profile of the
state of Eq. (81) with $N=L=21,\ X_{1}=-5,\ $and $X_{2}=5$.

###  4.2 Quantum double dark soliton with equal weight

####  4.2.1 Construction of the quantum state for a double dark soliton with
equal weights

By analogy with the fermionic case in the previous section, we now construct
quantum double dark solitons associated with a superposition of two-hole
excitation states in the Lieb-Liniger model. The choice of the Bethe quantum
numbers is the key to construct the quantum double dark soliton state.

A quantum state for a double dark soliton for $N$ particles is given by

$\displaystyle\ket{X_{1},X_{2},N}=\frac{1}{\sqrt{\mathcal{M}}}\sum_{\bm{p}\in\bm{P}}e^{i(p_{1}X_{1}+p_{2}X_{2})}\ket{p_{1},p_{2},N}$
(82)

with the normalization constant $\mathcal{M}$, and the two holes’ positions
$\bm{p}:=\\{p_{1},p_{2}\\}$ in the set $\bm{P}$ indicating all possible two
holes’ positions for the set of Bethe quantum numbers $I^{+}$. The set $I^{+}$
consists of the ground state of $N$ particles and additional two particles, as
shown in Figure 8. For example, in the case of $N=5$, $I^{+}$ is the set
$\\{-2,-1,0,1,2,3,4\\}$, and
$\bm{P}=\left\\{\\{-2,-1\\},\\{-2,0\\},\\{-2,1\\},\\{-2,2\\},\cdots,\\{2,4\\},\\{3,4\\}\right\\}$,
which has $|\bm{P}|=\binom{7}{2}=21$ elements. The quantum double dark soliton
state is a superposition of $|\bm{P}|=\binom{N+2}{2}$ two-hole excitation
states.

With respect to the state (82), the quantum double dark soliton with two
notches at $X_{1}$ and $X_{2}$ is expressed by the expectation value of the
local density operator

$\displaystyle\begin{split}\rho_{Q}(x):=&\braket{X_{1},X_{2},N}{\hat{\rho}(x)}{X_{1},X_{2},N}\\\
=&\frac{1}{\mathcal{M}}\sum_{\bm{p},\bm{p^{\prime}}\in\bm{P}}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})}\braket{p_{1}^{\prime},p_{2}^{\prime},N}{\hat{\rho}(0,0)}{p_{1},p_{2},N},\end{split}$
(83)

where $P$ is the total momentum of the state with the system length $L$, which
is given by

$\displaystyle P=\frac{2\pi}{L}\sum_{m\in I}m,\quad I\equiv I^{+}-\bm{p}.$
(84)

Again, the form factors
$\braket{P_{1},P_{2},N}{\hat{\rho}(0,0)}{P_{1}^{\prime},P_{2}^{\prime},N}$ are
computed by the determinant formula. We recall that the notch positions of the
quantum dark soliton are at $x=X_{1}$ and $X_{2}$, as in the case of the
quantum single dark solitons constructed by quantum states such as
$\ket{X(q),N}$ in section 2.1 and $\ket{N,X;P}$ in section 2.2, both of which
have a notch at the position $X$.

####  4.2.2 Density profile of the quantum double dark soliton state

It is quite marvelous how completely the density profiles $\rho_{Q}(x)$ of the
quantum double dark solitons in the weak coupling regime agree with those of
the corresponding elliptic multiple dark solitons (46), as shown in Figure 9.
It shows the quantum double dark solitons (83) and their corresponding
elliptic soliton solutions with different coupling constants. The quantum
double dark solitons are symmetrical about the line $x=0$ in all cases,
reflecting the fact that their quantum states are constructed by a
superposition of two-hole excitations with equal weight.

However, when the coupling constant is not small, each notch of the quantum
double dark-soliton is not as deep as the single soliton is. Furthermore, it
has two minima for the case $c=10$ and $100$, as shown in Figure 9. An
undulating line connects two quantum solitons, which tend to localize more in
a larger interaction regime.

Figure 8: Configurations of the Bethe quantum numbers. (i) Configuration $I^{+}$, which is obtained by adding two particles (yellow circles) to the ground state (red circles). (ii) Two holes, which correspond to a subset $\bm{p}=\\{p_{1},p_{2}\\}$ of $I^{+}$, are punched in the Bethe quantum numbers $I^{+}$ with the value of $P$ defined in (84).  Figure 9: Density profiles of the quantum double dark soliton $\rho_{Q}(x)$ (red solid curves) and the corresponding elliptic solutions (blue broken curves) for $N=L=20,\ X_{1}=L/4,\ X_{2}=-L/4,$ and coupling constants $c=0.05,\ 0.1,\ 10,\ \text{and }100$. The parameters of the elliptic solutions are shown in Table 3. Table 3: Parameters of the elliptic solutions for Figure 8. $c$ | $k$ | $\mu/2nc$ | $v_{c}/2\sqrt{nc}$ | $v/v_{c}$
---|---|---|---|---
0.05 | 0.681 750 | 1.112 62 | 1.657 91 | 0.986 753
0.1 | 0.840 623 | 1.115 19 | 1.325 36 | 0.959 493
10 | 1- 4.489 79$\times 10^{-9}$ | 1.006 60 | 0.998 333 | 0.814 300
100 | 1- 1.484 64$\times 10^{-24}$ | 1.001 34 | 0.999 156 | 0.841 665

####  4.2.3 Matrix elements of the quantum double dark soliton state

We have calculated the matrix element by using the determinant formula (17).
The matrix element $\psi_{Q}(x)$ is given by

$\displaystyle\begin{split}&\psi_{Q}(x):=\braket{X_{1},X_{2},N-1}{\hat{\psi}(x)}{X_{1},X_{2},N}\\\
=&\frac{1}{\mathcal{M}}\sum_{\bm{p},\bm{p^{\prime}}\in\bm{P}}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})}\braket{P_{1},P_{2},N-1}{\hat{\psi}(0,0)}{P_{1}^{\prime},P_{2}^{\prime},N},\end{split}$
(85)

which takes complex value for position $x$. The classical counterpart of the
matrix element $\psi_{Q}(x)$ is the relevant elliptic solution which is given
by Eq. (33), with its absolute part and argument. Comparisons of the matrix
element $\psi_{Q}(x)$ and the classical soliton $\psi_{C}(x)$ are made by
plotting the square amplitudes and the arguments. Note that the density $n$ is
not always $N/L$, since the number of particles in the matrix element is not
conserved.

Figure 10 shows the profile of the square amplitude $|\psi_{Q}(x)|^{2}$ for
the matrix element (85). The smallest values of the square amplitude of the
notches are smaller than those of the density profiles in Figure 8. In the
case of $c=0.05$ and 0.1, the square amplitudes shown in Figure 10 almost
completely agree with the corresponding elliptic solutions. In the case of
$c=10$ and 100, however, the widths of the dips still remain wider than those
of the elliptic solutions.

Figure 11 shows the phase profile $\arg[\psi_{Q}(x)]$ for the matrix element
(85). The phase profiles in Figure 11 agree with those of the elliptic
solutions for all coupling constants. It is noted that the phase profiles with
four different coupling constants $c$ shown in Figure 11 are similar for both
quantum and classical cases. We can confirm that the phase profile
$\arg[\psi_{Q}(x)]$ is satisfied with the periodic boundary condition such as
$\arg[\psi_{Q}(x+L)]=\arg[\psi_{Q}(x)]$ for the period $L$. The abrupt jumps
at the solitons positions in Figure 11 are consistent with our previous result
which says that the deep notch in a soliton amplitude gives the abrupt jump in
its phase profile Sato et al. (2016).

The phase profile of the notch located at $x=5$ is completely the same as that
of $x=-5$. The phase profile in total is point symmetric with respect to the
origin $x=0$.

Figure 10: Square amplitudes $|\psi_{Q}(x)|^{2}$ for the matrix elements (85) (red solid curves) and their classical counterparts (blue broken curves) for coupling constants $c=0.05,\ 0.1,\ 10,$ and 100. The parameters of the elliptic solutions are shown in Table 4. Figure 11: Phase profiles $\arg[\psi_{Q}(x)]/\pi$ for the matrix elements (85) (red solid curves) and their classical counterparts (blue broken curves) for coupling constants $c=0.05,\ 0.1,\ 10,$ and 100. The parameters of the elliptic solutions are shown in Table 4. In the all panels, the red solid curves completely agree with the blue broken curves. Table 4: Parameters of the elliptic solutions for Figures 10 and 11. $c$ | $k$ | $n$ | $\mu/2nc$ | $v_{c}/2\sqrt{nc}$ | $v/v_{c}$
---|---|---|---|---|---
0.05 | 0.681 520 | 0.666918 | 1.44 586 | 1.75 158 | 0.982 197
0.1 | 0.838 694 | 0.699671 | 1.40214 | 1.267 41 | 0.937 102
10 | 1- 4.444 99 $\times 10^{-9}$ | 0.411534 | 1.091 57 | 0.916 966 | 0.168 885
100 | 1- 9.964 63$\times 10^{-23}$ | 0.267601 | 1.037 57 | 0.963 609 | 0.0630 238

###  4.3 Quantum states of a double dark soliton with Gaussian weights

####  4.3.1 Gaussian weights in terms of target soliton depth

We now demonstrate that by assigning the Gaussian weights to the sum of
eigenstates in Eq. (82) quantum states for double dark-solitons are quite
effectively constructed. The method is similar to that of constructing quantum
single dark soliton states Shamailov and Brand (2019); Kaminishi et al.
(2020): We assign the ideal Gaussian weights Kaminishi et al. (2020) to the
sum of the eigenstates with two particle-hole excitations given in Eq. (82).
Then, the constructed quantum states give such density profiles that have
distinct notches with different depths.

We employ the ideal Gaussian weights Kaminishi et al. (2020) explained in
section 2.2. Hereafter in this subsection we express them in terms of the
target soliton depth denoted by $d$. We shall define it shortly. By assigning
a target soliton depth to each of the two notches, we construct a quantum
double soliton state such that its density profile has two distinct notches
with different depths.

Let us consider the Gaussian 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\ket{X_{1},X_{2},N}_{G}=\frac{1}{\sqrt{\mathcal{N}}}\sum_{\\{p_{1},p_{2}\\}\in\bm{P}}G_{P_{0}}(p_{1})G_{P^{\prime}_{0}}(p_{2})e^{i(p_{1}X_{1}+p_{2}X_{2})}\left|p_{1},p_{2},N\right\rangle.$
(86)

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

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

with two Gaussian parameters $(P,\sigma)$. They are determined by target
soliton depth $d$, as we shall show shortly.

We now define the target soliton depth $d$. It is related to the ratio
${v}/{v_{c,\infty}}$ through the dark soliton solution in the thermodynamic
limit Kaminishi et al. (2020) and defined by the following:

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

Here $|\phi_{P}^{\infty}(0)|$ denotes the square root of the local density
$\rho(x)$ at the origin, i.e., at the position of a notch, in the
thermodynamic limit. We recall that the local density $\rho(x)$ near a notch
has been given in Eq. (64). We also recall that the critical velocity
$v_{c,\infty}$ has been defined in section 3.3.1.

In terms of the target soliton depth $d$ and density $n=N/L$, the Gaussian
parameters $P$ and $\sigma$ are given by

$\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\\}\,,$
(89) $\displaystyle\sigma^{2}$
$\displaystyle=\frac{4}{3}n\sqrt{nc}\left(1-\frac{d}{n}\right)^{\frac{3}{2}}.$
(90)

####  4.3.2 Profiles of density $\rho_{Q}(x)$ and square amplitude
$|\psi_{Q}(x)|^{2}$

For the quantum state of Eq. (86) we exactly compute the density $\rho_{Q}(x)$
and the matrix element $\psi_{Q}(x)$ and plot the profiles of the density
$\rho_{Q}(x)$, the square amplitude $|\psi_{Q}(x)|^{2}$ and the phase
$\arg[\psi_{Q}(x)]$ in Figures 12,13, and 14, respectively, for several values
of the coupling constants. For all the figures in this subsection, the system
length $L$ and the number of particles $N$ are given by $N=L=20$. The target
soliton depths are set as $d=0.6$ and $d^{\prime}=0$ for each notch, so that
the pairs of Gaussian parameters in the quantum state of Eq. (86) are given by
$(P_{0},\sigma)=(0.124027\pi,0.106667)$ and
$(P^{\prime}_{0},\sigma^{\prime})=(\pi,0.421637)$.

In the density profiles of Figure 12 we observe double dark-solitons
consisting of two distinct narrow notches with different depths in the two
panels of the weak coupling regime for $c\leq 0.1$. In the two panels of
$c\leq 0.1$ it is clear that the depths of the two notches are different. We
consider it a consequence of assigning the Gaussian weights to the sum of
eigenstates with two particle-hole excitations. Thus, the plots in Figure 12
confirm that the construction with the ideal Gaussian weights reviewed in
section 2.2 is effective in the weak interaction regime such as $c\leq 0.1$.

We thus suggest from Figures 12 and 13 that the Gaussian superposition of the
eigenstates with two particle-hole excitations in Eq. (86) should produce a
quantum state rather close to the mean-field product state for a given double
dark-soliton, when the coupling constant $c$ is small enough such as $c\leq
0.1$. In particular, the almost complete agreement between the density
profiles in Figure 12 and the profiles of the square amplitude in Figure 13
for $c\leq 0.1$ suggests that the quantum state with the ideal Gaussian
weights in Eq. (86) should be close to the corresponding mean-field product
state. In the panels of $c=0.05$ and $0.1$, every notch in Figures 12 and 13
is distinct and narrow. Furthermore, the value of $\rho_{Q}(x)$ or
$|\psi_{Q}(x)|^{2}$ at the starting points of the dips, i.e., the maximal
points in the profiles, one of which is located around at $x=0$ for instance,
is given by 1.4 in both Figures 12 and 13.

Here we remark that if the Gaussian superposition of the eigenstates in the
Lieb-Liniger model of Eq. (86) completely produce the mean-field product state
satisfying Eq. (21), then the square amplitude of the matrix element
$\psi_{Q}(x)={}_{G}\braket{X_{1},X_{2},N-1}{\hat{\psi}(x)}{X_{1},X_{2},N}_{G}$
is equal to the density profile
$\rho_{Q}(x)={}_{G}\braket{X_{1},X_{2},N}{\hat{\psi}^{\dagger}(x)\hat{\psi}(x)}{X_{1},X_{2},N}_{G}$:
We have $\rho_{Q}(x)=|\psi_{Q}(x)|^{2}$.

Figure 12: Density profiles $\rho_{Q}(x)$ for $N=L=20$ and $c=0.05,\ 0.1,\
1.0,\ \text{and}\ 10$. The target soliton depths $d$ for the Gaussian weights
are given by $0.6$ and $0$, and the Gaussian parameters are given by
$(P_{0},\sigma)=(0.124027\pi,0.106667)$ and
$(P^{\prime}_{0},\sigma^{\prime})=(\pi,0.421637)$, which are calculated by Eqs
(89) and (90). The density profile is calculated in the same way as in Eq.
(83). Figure 13: Square amplitudes of the matrix element, $|\psi_{Q}(x)|^{2}$,
for $N=L=20,$ and $c=0.05,\ 0.1,\ 10,\ \text{and}\ 100$. The pairs of the
Gaussian parameters are given by the same sets of values as in Figure 12.

####  4.3.3 Difference between density $\rho_{Q}(x)$ and square amplitude
$|\psi_{Q}(x)|^{2}$

We now estimate how far the quantum soliton state constructed by Eq. (86)
satisfies the characteristic property of Eq. (21) for the mean-field product
state. Let us subtract the square amplitude of the matrix element,
$|\psi_{Q}(x)|^{2}$, from the density, $\rho_{Q}(x)$. Inserting the identity
operator between $\hat{\psi}^{\dagger}(x)$ and $\hat{\psi}(x)$ in the density
operator, we obtain the square amplitude of the matrix element and other
terms:

$\displaystyle\begin{split}&{}_{G}\braket{X_{1},X_{2},N}{\hat{\psi}^{\dagger}(x)\hat{\psi}(x)}{X_{1},X_{2},N}_{G}\\\
=&\sum_{\chi}{}_{G}\braket{X_{1},X_{2},N}{\hat{\psi}^{\dagger}(x)\ket{\chi}\bra{\chi}\hat{\psi}(x)}{X_{1},X_{2},N}_{G}\\\
=&|{}_{G}\braket{X_{1},X_{2},N-1}{\hat{\psi}(x)}{X_{1},X_{2},N}_{G}|^{2}+|\braket{\chi^{\prime}}{\hat{\psi}(x)}{X_{1},X_{2},N}_{G}|^{2}+\cdots.\end{split}$
(91)

We define the difference $D(x)$ by

$\displaystyle D(x):=$
$\displaystyle{}_{G}\braket{X_{1},X_{2},N}{\hat{\psi}^{\dagger}(x)\hat{\psi}(x)}{X_{1},X_{2},N}_{G}-|{}_{G}\braket{X_{1},X_{2},N-1}{\hat{\psi}(x)}{X_{1},X_{2},N}_{G}|^{2}.$
(92)

Since each term on the right hand side of Eq. (91) is non-negative, we have
the relation

$D(x)\geq 0.$ (93)

The difference $D(x)$ expresses how far the quantum soliton state is different
from the corresponding mean-field product state which should satisfy the
property of Eq. (21). Moreover, we remark that $D(x)$ depends on the position
$x$ and the coupling constant $c$.

We now argue that the quantum double dark soliton state of Eq. (86) should be
close to the corresponding mean-field product state in the weak coupling
regime such as $c\leq 0.1$. Let us evaluate $D(x)$ for weak and strong
coupling constants. When the coupling constant $c=0.05$, the difference
$D(x=0)$ is almost zero, while $D(x=5)=0.3$, from Figure 12 and 13 (see also
Table 5). On the other hand, for a strong coupling constant $c=10$, the
difference $D(x)$ is about 0.6 except for the notches.

Here we recall that in both Figures 12 and 13 around at $x=0$ one of the two
notches starts to dip, i.e., the value of $\rho_{Q}(x)$ or $|\psi_{Q}(x)|^{2}$
decreases, and it is approximately given by 1.4, as shown in Table 5.

Table 5: The values of the density profile$\rho_{Q}(x)$, those of the square amplitudes $|\psi_{Q}(x)|^{2}$, and their differences $D(x)$ at $x=0$ and $x=5$.  $c$ | $x=0$ |  | $x=5$
---|---|---|---
| $\rho_{Q}(x)$ | $|\psi_{Q}(x)|^{2}$ | $D(x)$ |  | $\rho_{Q}(x)$ | $|\psi_{Q}(x)|^{2}$ | $D(x)$
0.05 | 1.41865 | 1.38161 | 0.0370444 |  | 0.306279 | $9.80563\times 10^{-3}$ | 0.296473
0.1 | 1.44277 | 1.40292 | 0.0398479 |  | 0.283943 | $4.46456\times 10^{-3}$ | 0.279478
1.0 | 1.23462 | 1.03089 | 0.203732 |  | 0.334114 | $4.24589\times 10^{-6}$ | 0.334109
10 | 1.11822 | 0.528072 | 0.59015 |  | 0.528875 | $1.96588\times 10^{-4}$ | 0.528678

####  4.3.4 Phase profile of the quantum double soliton state

The phase profiles in Figure 14 are different from those of no Gaussian
weights plotted in Figure 11 with respect to the point symmetry and the
boundary conditions. The phase profiles do not have the point symmetry at the
origin in all the panels of the different values of the coupling constant in
Figure 14. Furthermore, except for the panel of $c=1.0$, the phase profile
satisfies the boundary condition:
$\arg[\psi_{g}(x+L)]=\arg[\psi_{g}(x)]+2\pi$. The phase profile is shifted by
$2\pi$ at the ends, while it is continuous with respect to the coordinate $x$.

Assigning the Gaussian weights leads to the change of the winding numbers for
the corresponding solitons. For an illustration, let us consider the upper-
right panel for $c=0.1$ of Figure 14. We can show that the winding number of
one soliton with its notch located at $x=-5$ is given by $J=0$, while that of
another located at $x=5$ is given by $J=1$. We shall explain it later in
detail.

Figure 14: Phase profiles $\arg[\psi_{g}(x)]/\pi$ for $N=L=20,$ and $c=0.05,\
0.1,\ 1.0,\ \text{and}\ 10$. The pairs of the parameters of the Gaussian
function are the same as Figure 12.

We decompose the phase profile in the upper-right panel for $c=0.1$ of Figure
14 around at $x=0$ into two phase sub-profiles of single dark solitons: one
sub-profile in the region of $x<0$ with an N-shaped figure centered at $x=-5$,
and another one in the region of $x>0$ with an abrupt jump at $x=5$. In Figure
15 the phase profile is divided by the vertical line at $x=0$ to two sub-
profiles in the regions of $x<0$ and $x>0$, respectively. The black dots
depict the centers of the phase sub-profiles. Here, they are located at $x=-5$
and 5, respectively.

We shall argue that the winding numbers of the two sub-profiles are given by
$J=0$ and $J=1$, respectively. Here we recall that the phase field is
approximately expressed by Eq. (68) if the system size $L$ is very large. We
also recall that elliptic multiple dark-soliton solutions are described in
section 3.2. The phase field $\varphi(x)$ is expressed in terms of the
elliptic integral of the third kind together with velocity $v$ and parameter
$W$ in Eq. (50) or more compactly in Eq. (56) in terms of the elliptic theta
function. The parameter $v$ is determined by the boundary condition (51) with
the winding number $J$. For the real parameter $W$, we have the choice of
signs, as shown in Eq. (37).

The left-hand side ($x<0$) of the phase profile in the upper-right panel
($c=0.1$) of Figure 14 has the same shape as the N shaped profile in the
leftmost panel of Figure 16. The unwinding phase profile ($J=0$) for $W<0$ and
that of $W>0$ are exhibited in the two panels from the left in Figure 16,
respectively, expressed as blue curves. Thus, we suggest that the left-hand
side phase profile in the region of $x<0$ in the upper-right panel ($c=0.1$)
of Figure 14 should correspond to the phase profile of a single elliptic dark-
soliton with $J=0$ and $W<0$. Here we recall that the abrupt decrease ($W<0$)
(or increase ($W>0$)) near the center ($x=0$) of the profile is explained in
section 3.3.2.

On the other hand, the right-hand side phase profile ($x>0$) in the upper-
right panel ($c=0.1$) of Figure 14 should correspond to the phase profile of a
single elliptic dark-soliton with $J=1$ and $W>0$. When the winding number
$J=1$, the term $4\pi/L$ is added to the velocity for $J=0$ due to the
boundary condition. It follows that the $J=1$ phase profile is given by the
superposition of the linear graph of $2\pi x/L$ and the $J=0$ phase profile.
In the middle panel of Figure 16, the line graph $\frac{2\pi}{L}x$ (black
broken line) and the phase profile with $J=0$ and $W>0$ (blue curve) are
plotted. The resultant winding phase profile with $J=1$ is shown in the
rightmost panel of Figure 16. We suggest that it should correspond to the
right-hand side phase profile ($x>0$) in the upper-right panel of Figure 14.

Thus, we have observed that by assigning the Gaussian weights with certain
parameters we obtain not only the density profiles consisting of two notches
with different soliton depths but also the phase profiles consisting of two
dark-solitons with different winding numbers.

Figure 15: Phase profile $\arg[\psi_{g}(x)]/\pi$ (red curve) for $N=L=20$ and
$c=0.1$ divided by the vertical line $x=0$ (blue line).  Figure 16: Phase
profiles of the elliptic soliton solution of Eq. (50). For $J=0$, $W<0$ and
hence $v>0$ (in the leftmost panel); for $J=0$, $W>0$ and hence $v<0$ (in the
middle panel); for $J=1$ and $W>0$ (in the rightmost panel). The black broken
line $\frac{2\pi x}{L}$ illustrates the first term of the phase profile in Eq.
(56) which is also contained in Eq. (50).

##  5 Conclusion

We have introduced novel quantum states of single dark solitons with a nonzero
winding number. Then, we have constructed the quantum double dark-soliton
states by taking the analogy to the case of the free fermions and the
construction of quantum single dark soliton states with and without Gaussian
weights. The superposition of the eigenstates with two particle-hole
excitations in the Lieb-Liniger model leads to such quantum states that have
two notches in the density profile as a double dark-soliton. Two notches with
the same depths appear in the superposition of two particle-hole excitations
with equal weight. Two notches with different depths are generated by those of
non-equal weights, in particular, the ideal Gaussian weights.

For quantum double dark solitons, we have controlled their depths by some
parameters. It has never been performed in the previous studies, as far as we
have investigated. We have obtained the quantum double dark soliton which
possesses distinct and narrow notches with controlled depths in both density
profiles and the profiles of the square amplitude of the matrix element of the
field operator, in the case when $N=L=20$ and $c=0.05$ and $0.1$.

Furthermore, we have confirmed that the quantum double dark soliton state
satisfies the assumption of the mean-field product state for small coupling
constant by evaluating the difference $D(x)$. That is, the density profiles
are close to the square amplitudes of the matrix elements of the field
operator if the coupling constant is small enough, such as $c<0.1$ for $N=20$.

As an important topic we have shown that the ansatz of Ref. Kaminishi et al.
(2020) is effective and useful as far as the density profile is concerned. The
ansatz was also quite useful for constructing double dark soliton states as we
have shown in §4. However, we should point out that it does not completely
lead to the mean-field product state. In fact, one of the key properties of
the mean-field product state, i.e., Eq. (21), holds only approximately for all
the examples investigated in this article. Furthermore, the phase field of the
constructed quantum state sometimes has a nonzero winding number, i.e., it
satisfies one of non-periodic boundary conditions. It is maybe due to the
finite-size effect. We expect that one can improve the construction so that
one obtain a complete mean-field product state by increasing the number of
particles $N$ and the system size $L$, while by decreasing the coupling
constant $c$.

##  6 Discussion

There is another approach to observe quantum dark-solitons numerically:
successive measurements of particle positions in the Lieb-Liniger model
studied by Syrwid et. al Syrwid and Sacha (2015); Syrwid et al. (2016). They
showed quantum single/double dark-solitons by measuring particle positions in
a single one-/two-hole excitation state. Even though their method are
different from ours in section 2.1 and 4.2, both studies share the results
that the phase profiles have the abrupt jumps at the solitons positions, and
the two notches in a double soliton have the same depths.

The question is how the density profile of a superposition of yrast states are
associated with the successive measurements of particle positions. For the
quantum single dark-soliton case, some groups worked on this question. When a
coupling constant $c=0$, it was analytically shown the relationship between
the construction with Gaussian weight Kaminishi et al. (2020) and the particle
position method Syrwid and Sacha (2015) in Ref. Kaminishi et al. (2020). In
the case of a small coupling constant $c>0$, an ansatz was proposed to bridge
between the calculation of single-particle density and the particle position
method, in a recent research Golletz et al. (2020).

On the other hand, for the quantum multi dark-soliton case, the connection
between the calculation of single-particle density and the particle position
method is still unclear. We suggest that mean-field state with the double dark
soliton solution of the NS equation could reveal the connection. The
connection is still unclear, but we showed some quantum states associated with
multiple dark-solitons and their properties explicitly.

We leave the comparison of the quantum double dark-soliton with Gaussian
weights and its classical counterpart for future work. The system we consider
here is under the periodic boundary conditions, while the soliton solution of
the NS eq., obtained by the inverse scattering method does not satisfy the
periodic boundary conditions. The finite-gap soliton Belokolos et al. (1994)
should be the one to compare with the quantum double dark-soliton with
Gaussian weights. For the appropriate parameters, the finite-gap soliton owns
two notches with the different depths. We will find the finite-gap solitons
which correspond to the quantum double dark-soliton with Gaussian weights.

Moreover, we shall study the dynamics of the quantum dark solitons, applying
the construction with Gaussian weights to them. We have not yet observed the
dynamics of the quantum double dark soliton even if both of the solitons are
black. The dynamics of the quantum dark solitons should reveal soliton
properties such as phase shifts.

The present research was partially supported by the Japan Society for the
Promotion of Science (KAKENHI Grants-in-Aid for Scientific Research: Grant No.
21K03398). E.K. was supported by JSPS Grant Number 20K14388.

###### Acknowledgements.

The authors declare no conflict of interest. References yes

## References

  * Pitaevskii and Stringari (2003) Pitaevskii, L.; Stringari, S.  Bose-Einstein Condensation; International series of monographs on physics, Clarendon Press, 2003.
  * Faddeev and Takhtajan (1987) Faddeev, L.D.; Takhtajan, L.A.  Hamiltonian methods in the theory of solitons; Classics in Mathematics, Springer-Verlag Berlin Heidelberg, 1987.
  * Tsuzuki (1971) Tsuzuki, T.  Nonlinear waves in the Pitaevskii-Gross equation.  Journal of Low Temperature Physics 1971, 4, 441–457.
  * Abdullaev et al. (1993) Abdullaev, F.; Darmanyan, S.; Khabidullaev, P.  Optical Solitons; Springer series in nonlinear dynamics, Springer-Verlag, 1993.
  * Burger et al. (1999) Burger, S.; Bongs, K.; Dettmer, S.; Ertmer, W.; Sengstock, K.; Sanpera, A.; Shlyapnikov, G.V.; Lewenstein, M.  Dark Solitons in Bose-Einstein Condensates.  Phys. Rev. Lett. 1999, 83, 5198–5201.  doi:black10.1103/PhysRevLett.83.5198.
  * Weller et al. (2008) Weller, A.; Ronzheimer, J.P.; Gross, C.; Esteve, J.; Oberthaler, M.K.; Frantzeskakis, D.J.; Theocharis, G.; Kevrekidis, P.G.  Experimental Observation of Oscillating and Interacting Matter Wave Dark Solitons.  Phys. Rev. Lett. 2008, 101, 130401.  doi:black10.1103/PhysRevLett.101.130401.
  * Becker et al. (2008) Becker, C.; Stellmer, S.; Soltan-Panahi, P.; Dörscher, S.; Baumert, M.; Richter, E.M.; Kronjäger, J.; Bongs, K.; Sengstock, K.  Oscillations and interactions of dark and dark–bright solitons in Bose–Einstein condensates.  Nature Physics 2008, 4, 496–501.
  * Wadachi and Sakagami (1984) Wadachi, M.; Sakagami, M.a.  Classical Soliton as a Limit of the Quantum Field Theory.  Journal of the Physical Society of Japan 1984, 53, 1933–1938, [https://doi.org/10.1143/JPSJ.53.1933].  doi:black10.1143/JPSJ.53.1933.
  * Ishikawa and Takayama (1980) Ishikawa, M.; Takayama, H.  Solitons in a One-Dimensional Bose System with the Repulsive Delta-Function Interaction.  Journal of the Physical Society of Japan 1980, 49, 1242–1246, [https://doi.org/10.1143/JPSJ.49.1242].  doi:black10.1143/JPSJ.49.1242.
  * Sato et al. (2012) Sato, J.; Kanamoto, R.; Kaminishi, E.; Deguchi, T.  Exact Relaxation Dynamics of a Localized Many-Body State in the 1D Bose Gas.  Phys. Rev. Lett. 2012, 108, 110401.
  * Sato et al. (2016) Sato, J.; Kanamoto, R.; Kaminishi, E.; Deguchi, T.  Quantum states of dark solitons in the 1D Bose gas.  New Journal of Physics 2016, 18, 075008.  doi:black10.1088/1367-2630/18/7/075008.
  * Syrwid and Sacha (2015) Syrwid, A.; Sacha, K.  Lieb-Liniger model: Emergence of dark solitons in the course of measurements of particle positions.  Phys. Rev. A 2015, 92, 032110.  doi:black10.1103/PhysRevA.92.032110.
  * Shamailov and Brand (2019) Shamailov, S.S.; Brand, J.  Quantum dark solitons in the one-dimensional Bose gas.  Phys. Rev. A 2019, 99, 043632.  doi:black10.1103/PhysRevA.99.043632.
  * Kaminishi et al. (2020) Kaminishi, E.; Mori, T.; Miyashita, S.  Construction of quantum dark soliton in one-dimensional Bose gas.  Journal of Physics B: Atomic, Molecular and Optical Physics 2020, 53, 095302.
  * Golletz et al. (2020) Golletz, W.; Górecki, W.; Ołdziejewski, R.; Pawłowski, K.  Dark solitons revealed in Lieb-Liniger eigenstates.  Phys. Rev. Research 2020, 2, 033368.  doi:black10.1103/PhysRevResearch.2.033368.
  * Zakharov and Shabat (1973) Zakharov, V.E.; Shabat, A.B.  Interaction between solitons in a stable medium.  Sov. Phys. JETP 1973, 37, 823–828.
  * Belokolos et al. (1994) Belokolos, E.; Bobenko, A.; Enol’skii, V.; Its, A.; Matveev, V.  Algebro-geometric approach to nonlinear integrable equations; Springer series in nonlinear dynamics, Springer-Verlag, 1994.
  * Lieb and Liniger (1963) Lieb, E.H.; Liniger, W.  Exact analysis of an interacting Bose gas. I. The general solution and the ground state.  Physical Review 1963, 130, 1605.
  * Bethe (1931) Bethe, H.  Eigenvalues and eigenfunctions of the linear atom chain.  Z. Phys 1931, 71, 205–226.
  * Lieb (1963) Lieb, E.H.  Exact analysis of an interacting Bose gas. I. The general solution and the ground state.  Physical Review 1963, 130, 1616.
  * Korepin et al. (1997) Korepin, V.E.; Bogoliubov, N.M.; Izergin, A.G.  Quantum inverse scattering method and correlation functions; Vol. 3, Cambridge university press, 1997.
  * Gaudin (1983) Gaudin, M.  La fonction d’onde de Bethe; Masson, 1983.
  * Korepin (1982) Korepin, V.E.  Calculation of norms of Bethe wave functions.  Communications in Mathematical Physics 1982, 86, 391–418.
  * Slavnov (1989) Slavnov, N.A.  Calculation of scalar products of wave functions and form factors in the framework of the algebraic Bethe ansatz.  Teoreticheskaya i Matematicheskaya Fizika 1989, 79, 232–240.
  * Slavnov (1990) Slavnov, N.A.  Nonequal-time current correlation function in a one-dimensional Bose gas.  Theoretical and Mathematical Physics 1990, 82, 273–282.
  * Calabrese and Caux (2007) Calabrese, P.; Caux, J.S.  Dynamics of the attractive 1D Bose gas: analytical treatment from integrability.  Journal of Statistical Mechanics: Theory and Experiment 2007, 2007, P08032.  doi:black10.1088/1742-5468/2007/08/p08032.
  * Kojima et al. (1997) Kojima, T.; Korepin, V.E.; Slavnov, N.  Determinant representation for dynamical correlation functions of the quantum nonlinear Schrödinger equation.  Communications in mathematical physics 1997, 188, 657–689.
  * Sato et al. (2013) Sato, J.; Kaminishi, E.; Deguchi, T.  Finite-size scaling behavior of Bose-Einstein condensation in the 1D Bose gas, 2013, [arXiv:cond-mat.quant-gas/1303.2775].
  * Kanamoto et al. (2009) Kanamoto, R.; Carr, L.D.; Ueda, M.  Metastable quantum phase transitions in a periodic one-dimensional Bose gas: Mean-field and Bogoliubov analyses.  Phys. Rev. A 2009, 79, 063616.  doi:black10.1103/PhysRevA.79.063616.
  * Copson (1978) Copson, E.  An Introduction to the Theory of Functions of a Complex Variable; Springer series in nonlinear dynamics, Oxford University Press, 1978\.
  * Girardeau (1960) Girardeau, M.  Relationship between Systems of Impenetrable Bosons and Fermions in One Dimension.  J. Math. Phys. 1960, 1, 516.  doi:black10.1063/1.1703687.
  * Syrwid et al. (2016) Syrwid, A.; Brewczyk, M.; Gajda, M.; Sacha, K.  Single-shot simulations of dynamics of quantum dark solitons.  Phys. Rev. A 2016, 94, 023623.  doi:black10.1103/PhysRevA.94.023623.
  * Golletz et al. (2020) Golletz, W.; Górecki, W.; Ołdziejewski, R.; Pawłowski, K.  Dark solitons revealed in Lieb-Liniger eigenstates.  Phys. Rev. Research 2020, 2, 033368.  doi:black10.1103/PhysRevResearch.2.033368.