Quantum Dark Solitons in the 1D Bose Gas: From Single to Double Dark-Solitons

Abstract

We study quantum double dark-solitons, which give pairs of notches in the density profiles, 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.

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 traps [1]. It is expressed in terms of the classical complex scalar field $\psi$ as follows:

$$\begin{array}{c}\hfill i\frac{\partial }{\partial t}\psi =-\frac{{\partial }^2}{\partial x^2}\psi +2c{\left|\psi \right|}^2\psi -\mu \psi .\end{array}$$

(1)

The nonlinearity 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 Schrödinger (NS) equation, which is one of the fundamental integrable nonlinear equations characterized by multi-soliton solutions [2]. 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 [3]. Here, we note that bright solitons are actively studied in association with optical fibers [4], while dark solitons are observed in cold atomic experiments [5,6,7].

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 [8]. 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 solution [9]. Here, we note 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 as those that a dark soliton have constructed [10,11]; Then, it was found that the density profile of a single grey soliton appears in the quantum measurement process [12,13]. Furthermore, the construction of quantum states of grey solitons has been elaborated with the Gaussian weights [14,15]. In particular, an ansatz on how to assign the Gaussian weights to the yrast states was proposed [15]. The matrix elements of the bosonic field operator of the Lieb–Liniger model are characterized by grey solitons [11], and double elliptic dark solitons are derived in association with the quantum measurement process [16]. Here, we note that solutions of the NS equation in the periodic system are expressed in terms of Jacobi’s elliptic functions. In this paper, we call such solutions elliptic solitons.

We define quantum solitons by the quantum states that have solitons in the density and phase profiles and appear in interacting or many-body quantum systems. It is nontrivial to construct such quantum states that have localized excitations in the density profiles in interacting quantum systems such as quantum dark solitons. In the case of interacting quantum systems, it should be hopeless to construct quantum states that have designed density profiles straightforwardly by controlling the dispersion and nonlinearity of the systems. However, in terms of the quantum dark solitons, we can construct quantum states which have localized structures in the density profiles. We thus suggest that we might control quantum states in interacting quantum systems such as the one-dimensional, bose gas through quantum dark solitons [17]. We note that classical solitons appear in the balance of dispersion and nonlinearity. We also note that there are interesting studies on more composite objects such as dark-bright and dark-antidark solitons by numerically solving the GP equation [18].

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. [15] 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 Section 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 [15] for constructing quantum dark solitons in terms of the ideal Gaussian weights in Section 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 Section 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 Section 3.2, we derive the elliptic multiple dark soliton solution of the GP equation. We then argue analytically in Section 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 Section 4 the density profiles of quantum double dark-solitons explicitly. For an illustration, we discuss in Section 4.1 the results of the free fermions in one dimension. In Section 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 Section 4.3, assigning the ideal Gaussian weights of the ansatz [15] 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 obtains a complete mean-field product state by increasing the number of particles N and the system size L, while decreasing the coupling constant c. Here, we have assumed that, if the coupling constant is very small, the condensate fraction of BEC should be nonzero finite so that the mean-field picture should hold in the quantum soliton state. We believe that, particularly in Section 3.1, Section 4.2 and Section 4.3, new results are presented.

We now give some notes 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 [3]. The multi-soliton solutions of the NS equation were constructed by Zakharov–Shabat, in particular, in the case of attractive nonlinear interaction [19]. For the repulsive case, multi-soliton solutions were systematically derived by Faddeev and Takhtajan through the inverse scattering method [2]. The multi dark soliton solutions are called ”soliton solutions in the case of finite density”, in the reference [2]. Here, we should note 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 solitons [20]. 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 research in theoretical physics to the useful references in mathematics.

Finally, we menton some potentially relevant results in the following. Bound states of dark solitons are numerically studied by solving the GP equation [21], dynamics of a bright soliton in BEC with time-dependent atomic scattering length in an repulsive parabolic potential [22], quantized quasi-two-dimensional Bose–Einstein condensates with spatially modulated nonlinearity [23], matter rogue wave in Bose–Einstein condensates with attractive atomic interaction [24], exact soliton solutions, and nonlinear modulation instability in spinor Bose–Einstein condensates [25].

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 [26]

$$\begin{array}{c}\hfill {\widehat{\mathcal{H}}}_{\mathrm{LL}}=-\sum _{j=1}^N\frac{{\partial }^2}{\partial x_j^2}+2c\sum _{j<k}^N\delta (x_j-x_k)-\mu N,\end{array}$$

(2)

where N is the number of particles, c is the coupling constant, and $\mu$ is the chemical potential. We note that the Lieb–Liniger model is often characterized by the parameter $\gamma =c/n$ with the particle density $n=N/L$. We express it in terms of the coupling constant c instead of the parameter $\gamma$ throughout this paper since we assume $n=1$ that is, $N=L$. We impose a periodic boundary condition with length L. In the second quantized form, it reads

$$\begin{array}{c}\hfill {\widehat{\mathcal{H}}}_{\mathrm{LL}}={\int }_0^{L}dx\left[{\partial }_x{\widehat{\psi }}^{†}{\partial }_x\widehat{\psi }+c{\widehat{\psi }}^{†}{\widehat{\psi }}^{†}\widehat{\psi }\widehat{\psi }-\mu {\widehat{\psi }}^{†}\widehat{\psi }\right],\end{array}$$

(3)

where the operators $\widehat{\psi }\left(x\right)$ and ${\widehat{\psi }}^{†}\left(x\right)$ are the bosonic annihilation and creation operators satisfying the canonical commutation relations

$$\begin{array}{cc}\hfill & [\widehat{\psi }\left(x\right),{\widehat{\psi }}^{†}\left(x^{\prime }\right)]=\delta (x-x^{\prime }),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & [\widehat{\psi }\left(x\right),\widehat{\psi }\left(x^{\prime }\right)]=[{\widehat{\psi }}^{†}\left(x\right),{\widehat{\psi }}^{†}\left(x^{\prime }\right)]=0.\hfill \end{array}$$

(5)

The exact eigenstates and eigenenergies of this system are obtained through the Bethe ansatz method [26,27,28]. The central object of this method is the Bethe ansatz equations

$$\begin{array}{c}\hfill e^{i{k}_{j}L}=\prod _{\ell \ne j}^N\frac{k_j-k_{\ell }+ic}{k_j-k_{\ell }-ic}\mathrm{for}j=1,2,\dots ,N.\end{array}$$

(6)

Each solution of the equations has its corresponding exact wavefunction expressed in terms of the summation over the symmetric group $S_N$ of N particles

$$\begin{array}{cc}\hfill & {\phi }_{k_1,\dots ,k_N}(x_1,\dots ,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],\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & A_{\sigma }={(-1)}^{\sigma }\prod _{j>\ell }^N\left[k_{{\sigma }_j}-k_{{\sigma }_{\ell }}-ic\mathrm{sign}(x_j-x_{\ell })\right],\hfill \end{array}$$

(8)

where the symbol ${\sigma }_j$ denotes the element to which the permutation $\sigma$ of N particles maps the integer j. The energy E and momentum P eigenvalues are given by

$$\begin{array}{c}\hfill E=\sum _{j=1}^N{k}_j^2,P=\sum _{j=1}^N{k}_j.\end{array}$$

(9)

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

$$\begin{array}{c}\hfill k_{j}L=2\pi I_j-2\sum _{\ell \ne j}^{N}arctan\left(\frac{k_j-k_{\ell }}{c}\right)\mathrm{for}j=1,2,\dots ,N,\end{array}$$

(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) have a unique real solution for a given set of Bethe quantum numbers [29]. In particular, the ground state is given by the set $I_j=-(N-1)/2+(j-1)$ for $j=1,2,\dots ,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 [8]. In this work, we consider the repulsive case $c>0$, where the mean field equation has the dark soliton solution [3,19]. In refs. [10,11], it is shown that the quantum states corresponding to the dark soliton solutions are the superpositions of the hole excitations

$$\begin{array}{c}\hfill |X\left(q\right),N\rangle :=\frac{1}{\sqrt{N}}\sum _{p=0}^{N-1}exp(-2\pi ipq/N)|P\left(p\right),N\rangle .\end{array}$$

(11)

Here, $|P\left(p\right),N\rangle$ is the one-hole excitation with the total momentum $P\left(p\right)=2\pi p/L$, whose quantum numbers are

$$\begin{gathered} I_j=-(N+1)/2+jfor1\le j\le N-p, \\ =-(N+1)/2+j+1forN-p+1\le j\le N. \end{gathered}$$

(12)

The density profile of this state $\langle X,N|{\widehat{\psi }}^{†}\left(x\right)\widehat{\psi }\left(x\right)|X,N\rangle$ shows a density notch at the position $X\left(q\right)=(q/N+1/2)L$, which coincides with the squared amplitude of the classical dark soliton solution [10,11]. Here, by the determinant formula for the norms of Bethe eigenstates [30,31], we can effectively evaluate the matrix element

$$\begin{gathered} \langle X,N|{\widehat{\psi }}^{†}\left(x\right)\widehat{\psi }\left(x\right)|X,N\rangle \\ =\frac{1}{N}\sum _{p=0}^{N-1}\sum _{p^{\prime }=0}^{N-1}exp\left[2\pi i(p-p^{\prime })\frac{(x-X)}{L}\right]\langle P^{\prime },N|{\widehat{\psi }}^{†}\left(0\right)\widehat{\psi }\left(0\right)|P,N\rangle \end{gathered}$$

(13)

and that of the form factors of the density operator [32,33,34]

$$\begin{gathered} \langle P^{\prime },N|{\widehat{\psi }}^{†}\left(0\right)\widehat{\psi }\left(0\right)|P,N\rangle \\ ={(-1)}^{N(N+1)/2}(P-P^{\prime })\left(\prod _{j,\ell =1}^N\frac{1}{k_j^{\prime }-k_{\ell }}\right)\left(\prod _{j>\ell }^N{k}_{j,\ell }{k}_{j,\ell }^{\prime }\sqrt{\frac{\widehat{K}\left(k_{j,\ell }^{\prime }\right)}{\widehat{K}\left(k_{j,\ell }\right)}}\right)\frac{detU(k,k^{\prime })}{\sqrt{detG\left(k\right)detG\left(k^{\prime }\right)}}, \end{gathered}$$

(14)

where the quasimomenta $\{k_1,\cdots ,k_N\}$ and $\{k_1^{\prime },\cdots ,k_N^{\prime }\}$ give the eigenstates $|P\rangle$ and $|P^{\prime }\rangle$, respectively. We use the abbreviations $k_{j,\ell }:=k_j-k_{\ell }$ and $k_{j,\ell }^{\prime }:=k_j^{\prime }-k_{\ell }^{\prime }$. The kernel $\widehat{K}\left(k\right)$ is defined by $\widehat{K}\left(k\right)=2c/(k^2+c^2)$. The matrix $G\left(k\right)$ is called the Gaudin matrix, whose $(j,\ell )$ th element is given by

$$G{\left(k\right)}_{j,\ell }={\delta }_{j,\ell }\left[L+\sum _{m=1}^N\widehat{K}\left(k_{j,m}\right)\right]-\widehat{K}\left(k_{j,\ell }\right)forj,\ell =1,2,\cdots ,N.$$

(15)

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

$$\begin{array}{c}\hfill U{(k,k^{\prime })}_{j,\ell }=2{\delta }_{j\ell }\mathrm{Im}\left[\prod _{a=1}^N\frac{k_a^{\prime }-k_j+ic}{k_a-k_j+ic}\right]+\frac{{\prod }_{a=1}^N(k_a^{\prime }-k_j)}{{\prod }_{a\ne j}^N(k_a-k_j)}\left(\widehat{K}\left(k_{j,\ell }\right)-\widehat{K}\left(k_{N,\ell }\right)\right).\end{array}$$

(16)

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

$$\begin{gathered} \begin{array}{c}\hfill \begin{array}{cc}\hfill & \langle X^{\prime },N-1|\widehat{\psi }\left(x\right)|X,N\rangle \hfill \end{array} \\ =\frac{1}{\sqrt{N(N-1)}}\sum _{p=0}^{N-1}\sum _{p^{\prime }=0}^{N-2}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]\langle P^{\prime },N-1|\widehat{\psi }\left(0\right)|P,N\rangle ,\end{array} \end{gathered}$$

(17)

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 $|P,N\rangle$ and $|P^{\prime },N\rangle$, respectively. The determinant formula is given by [30,31,32,33,34,35]

$$\begin{gathered} \langle P^{\prime },N-1|\widehat{\psi }\left(0\right)|P,N\rangle ={(-1)}^{N(N+1)/2+1} \\ \times \left(\prod _{j=1}^{N-1}\prod _{\ell =1}^N\frac{1}{k_j^{\prime }-k_{\ell }}\right)\left(\prod _{j>\ell }^N{k}_{j,\ell }\sqrt{k_{j,\ell }^2+c^2}\right)\left(\prod _{j>\ell }^{N-1}\frac{k_{j,\ell }^{\prime }}{\sqrt{{\left(k_{j,\ell }^{\prime }\right)}^2+c^2}}\right)\frac{det\widehat{U}(k,k^{\prime })}{\sqrt{detG\left(k\right)detG\left(k^{\prime }\right)}}, \end{gathered}$$

(18)

where the quasi-momenta $\{k_1,\cdots ,k_N\}$ and $\{k_1^{\prime },\cdots ,k_{N-1}^{\prime }\}$ give the eigenstates $|P,N\rangle$ and $|P^{\prime },N-1\rangle$, respectively. Here, we note that the symbol $X^{{}^{\prime }}$ corresponds to an integer $q^{{}^{\prime }}$. We recall that the matrix $G\left(k\right)$ denotes the Gaudin matrix, whose $(j,\ell )$th element is given in (15). The matrix elements of the $(N-1)$ by $(N-1)$ matrix $\widehat{U}(k,k^{\prime })$ are given by

$$\begin{array}{cc}\hfill \widehat{U}{(k,k^{\prime })}_{j,\ell }& =2{\delta }_{j\ell }\mathrm{Im}\left[\frac{{\prod }_{a=1}^{N-1}(k_a^{\prime }-k_j+ic)}{{\prod }_{a=1}^N(k_a-k_j+ic)}\right]+\frac{{\prod }_{a=1}^{N-1}(k_a^{\prime }-k_j)}{{\prod }_{a\ne j}^N(k_a-k_j)}\left(\widehat{K}\left(k_{j,\ell }\right)-\widehat{K}\left(k_{N,\ell }\right)\right).\hfill \end{array}$$

(19)

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 [15] 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 ${\varphi }_P\left(x\right)$ 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. [15]: 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 [36].

The proposed BEC state [15], i.e., the quantum state introduced with the wave function ${\varphi }_P$, which we shall call the mean-field product state, exactly reproduces the density and phase profiles of the classical dark-soliton ${\varphi }_P$. Furthermore, it was numerically studied [15] that the main components of the product state are given by a weighted sum of the yrast states. Here, we note 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 [15]. For a given macroscopic wave function $\varphi \left(x\right)$ 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 $\varphi \left(x\right)$, the N-particle state is given by

$$\frac{1}{\sqrt{N!}}{\left({\int }_{-L/2}^{L/2}dx\varphi \left(x\right){\widehat{\psi }}^{†}\left(x\right)\right)}^N|\Omega \rangle ,$$

where $|\Omega \rangle$ denotes the vacuum. Here, the number N may be macroscopically large. We call it the mean-field product state associated with the wave function $\varphi \left(x\right)$ in the paper. In ref. [15], as an N-particle quantum soliton state corresponding to a classical dark soliton solution ${\varphi }_P(x-X)$, the following BEC state is proposed:

$$|N,X;P\rangle =\frac{1}{\sqrt{N!}}{\left({\int }_{-L/2}^{L/2}dx\frac{{\varphi }_P(x-X)}{\sqrt{N}}{\widehat{\psi }}^{†}\left(x\right)\right)}^N|\Omega \rangle .$$

(20)

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

$$\langle N,X;P|{\widehat{\psi }}^{†}\left(x\right)\widehat{\psi }\left(x\right){|N,X;P\rangle }_{\mathrm{BEC}}={\left|{\varphi }_P(x-X)\right|}^2,$$

(21)

and satisfies

$$\widehat{\psi }\left(x\right)|N,X;P\rangle ={\varphi }_P(x-X)|N-1,X;P\rangle .$$

(22)

Here, we note that an explicit analytic expression of a dark soliton solution ${\varphi }_P\left(x\right)$ in a finite system will be given in Section 3.2 [11].

It should be nontrivial to express the mean-field product state $|N,X;P\rangle$ in terms of the energy eigenstates of the Lieb–Liniger model. However, an ansatz is given in ref. [15] 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 $|N,X;P\rangle$ in terms of the eigenstates $|{\left\{I_j\right\}}_N\rangle$ of the Lieb–Liniger model as follows:

$$|N,X;P\rangle =\sum _{\left\{I_j\right\}}{C}_{\left\{I_j\right\}}|{\left\{I_j\right\}}_N\rangle .$$

(23)

Quite an important numerical observation is made in ref. [15] in which 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 Equation (20) should be essentially expressed as a superposition of the yrast states.

The expansion coefficients of the yrast states in the weak-coupling regime are evaluated in ref. [15] in the large N limit. We now explain the method about 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 $\widehat{P}={\int }_{-L/2}^{L/2}dx{\widehat{\psi }}^{†}\left(x\right)(-i{\partial }_x)\widehat{\psi }\left(x\right)$ are given by

$$\underset{N\to \infty }{lim}\langle N,X;P|\widehat{P}|N,X;P\rangle =P$$

(24)

and

$$\underset{N\to \infty }{lim}\langle N,X;P|{(\widehat{P}-P)}^2|N,X;P\rangle =\frac{4}{3}{\beta }^{3}n\sqrt{nc}\equiv {\sigma }_P^2,$$

(25)

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_{c,\infty }^2$. We note that n and $\beta$ correspond to ${\rho }_0$ and $\gamma$ in ref. [15]. Associated with Equations (24) and (25), we assume that, in the BEC state $|N,X;P\rangle$, the momenta ${\left\{p_j\right\}}_{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:

$$\begin{array}{cc}\hfill |N,X;P\rangle & \approx {\mathcal{N}}^{-1/2}\sum _{P^{\prime }}{e}^{i{P}^{\prime }(X-L/2)}{e}^{-G_P\left(P^{\prime }\right)/2}{|N,P^{\prime }\rangle }_{\mathrm{yr}}\hfill \\ \hfill & \equiv {|N,X;P\rangle }_{\mathrm{yr}},\hfill \end{array}$$

(26)

where

$$G_P\left(P^{\prime }\right)=\frac{{(P^{\prime }-P)}^2}{2{\sigma }_P^2}$$

(27)

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

In Figure 2, it is clear that the density profiles of the quantum and classical solitons completely agree with each other. Here, we note that the density profiles for the quantum soliton state of Equation (26), $\langle N,0,P|{\widehat{\psi }}^{†}\left(x\right)\widehat{\psi }\left(x\right)|N,0,P\rangle$ with $P=\pi n$ (the black soliton) and $P=(\pi /2-1)n$ (a gray soliton) are plotted in red curves for $N=L=100$ and $c=0.01$ together with the corresponding classical solutions of the GP equation in blue dashed curves in Figure 2. Furthermore, it is also clear that the expression of Equation (26) is much more useful for numerical calculation than Equation (20). The range of the summation in Equation (26) 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.

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 Equation (26) with Equation (27). In fact, the ansatz with the ideal Gaussian weights [15] 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 Equation (26) compared to the uniform superposition introduced in ref. [10] 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 Equations (26) and (27) 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 Equation (20) and the analytic expression of ${\sigma }_P$ is obtained. Here, it should be noted that the idea of a Gaussian superposition of yrast states has also been proposed in a recent work [14], where the Gaussian form is assumed and the Gaussian width ${\sigma }_P$ is given by a fitting parameter.

We note that quantum depletion was observed in experiments [37]. It is suggested that the quantum depletion should lead to a less black profile than the corresponding classical one. However, for the one-dimensional system in weak interacting regime such as $c=0.01$, we expect that the effect of quantum depletion could be small since the scattering length should be small for $c=0.01$, so that the quantum profile may coincide with the classical one.

2.2.3. Numerical Study on the Validity of the Ansatz with the Ideal Gaussian Weights

The ansatz with the ideal Gaussian weights [15] 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 Equation (22) 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 note that the soliton momentum P and the variance ${\sigma }_P^2$ are determined by the target soliton depth d through Equations (93) and (94). We also note that the definition of target soliton depth will be given in Equation (92).

For the quantum state associated with soliton momentum P in Equation (26), the density profile is exactly evaluated by

$$\begin{array}{cc}\hfill & \langle N,X;P|{\widehat{\psi }}^{†}\left(x\right)\widehat{\psi }\left(x\right)|N,X;P\rangle \hfill \\ & =\frac{1}{N}\sum _{p=0}^N\sum _{p^{\prime }=0}^N{e}^{-G_P\left(P\right)/2}{e}^{-G_P\left(P^{\prime }\right)/2}exp\left[2\pi i(p-p^{\prime })\frac{(x-X)}{L}\right]\langle P^{\prime },N|{\widehat{\psi }}^{†}\left(0\right)\widehat{\psi }\left(0\right)|P,N\rangle ,\hfill \end{array}$$

(28)

and the matrix element of the field operator between the quantum soliton states with N and $N-1$ particles by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \langle N-1,X^{\prime };P|\widehat{\psi }\left(x\right)|N,X;P\rangle \hfill \\ \hfill & =\frac{1}{\sqrt{N(N-1)}}\sum _{p=0}^N\sum _{p^{\prime }=0}^{N-1}{e}^{-G_P\left(P\right)/2}{e}^{-G_P\left(P^{\prime }\right)/2}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]\hfill \\ \hfill & \times \langle P^{\prime },N-1|\widehat{\psi }\left(0\right)|P,N\rangle .\hfill \end{array}\end{array}$$

(29)

Here, $\langle P^{\prime },N|{\widehat{\psi }}^{†}\left(0\right)\widehat{\psi }\left(0\right)|P,N\rangle$ and $\langle P^{\prime },N-1|\widehat{\psi }\left(0\right)|P,N\rangle$ are computed by the determinant formula as explained in Section 2.1. 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 note that the degree of overlap of the two profiles suggests the degree of how far the quantum soliton state constructed by the ansatz of the ideal Gaussian weights [15] 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). We compare the density profiles and the square-amplitude profiles for the quantum states represented in Equations (11) and (26). The profiles shown in the upper left panel of Figure 3 for the state in Equation (11) are evaluated by Equations (13) and (17), while those shown in the other panels of Figure 3 for the state in Equation (26) are evaluated by Equations (28) and (29). In Figure 3, we observe that, in each value of target soliton depth d, which will be defined by Equation (92), 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.

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 $\varphi \left(x\right)=\sqrt{\rho \left(z\right)}exp\left[i\phi \left(x\right)\right]$ satisfying the periodic boundary conditions

$$\phi (x+L)=\phi \left(x\right)+2\pi J,$$

(30)

where J is an arbitrary integer. It is called the phase-winding number [38]. 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 [11,15]. Let us consider a new quantum state in the following:

$$\begin{array}{c}\hfill {|X\left(q\right),N\rangle }_{J=1}:=\frac{1}{\sqrt{N}}\sum _{p=1}^{N}exp(-2\pi ipq/N)|P\left(p\right),N\rangle .\end{array}$$

(31)

The difference between the quantum states (11) and (31) 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 Equation (31) is obtained from that of Equation (11) by removing the ground state ($p=0$) and adding the quantum state with momentum $P=2\pi /L(p=N)$. We suggest that the shift increases the momentum of the quantum state so that the winding number J increases.

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

$${\psi }_Q\left(x\right):={\left[\langle X,N-1|\widehat{\psi }\left(x\right)|X,N\rangle \right]}_{J=1},$$

(32)

and write ${\psi }_Q\left(x\right)$ in the polar form: ${\psi }_Q\left(x\right)=\left|{\psi }_Q\left(x\right)\right|exp\left[iarg\left[{\psi }_Q\left(x\right)\right]\right]$. The matrix element is computed by the determinant formula such as Equation (18).

3.1.2. Profiles of the Square Amplitude and the Phase of the Matrix Element of the Field Operator

Figure 4 and Figure 5 present the profiles of the square amplitude $|{\psi }_Q{\left(x\right)|}^2$ and those of the phase field $arg\left[{\psi }_Q\left(x\right)\right]$ in red solid curves, respectively, for four different values of the coupling constant such as $c=0.01,0.1,10$ and 100. Although the system is in the strong coupling regime in the cases of $c=10,$ and 100, we have displayed these plots for comparison. 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.1. Here, $a_1$ denotes the smallest root of the equation: $U\left(\rho \right)=0$, where the potential $U\left(\rho \right)$ will be given by Equation (39). We can draw the density profile of the corresponding elliptic dark-soliton by applying Equation (49), as will be shown in Section 3.2. Here, we note 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 Equations (61) and (62), 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$.

In Figure 4, the plots of the square amplitudes $|{\psi }_Q{\left(x\right)|}^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{\left(x\right)|}^2$ is wider than that of the elliptic solution. Moreover, the oscillating behavior around the notch is observed as shown in Figure 4. This behavior is a feature of the density profile in free fermions indicated in Section 4.1.

In Figure 5, all the phase profiles $arg\left[{\psi }_Q\left(x\right)\right]$ exhibited satisfy the boundary condition

$$\begin{array}{c}\hfill arg\left[{\psi }_Q(x+L)\right]=arg\left[{\psi }_Q\left(x\right)\right]+2\pi .\end{array}$$

(33)

It corresponds to the case of $J=1$ in Equation (30). 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 Equation (37). 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 Figure 4 and Figure 5 are listed in

Table 1. Parameters of the profiles of the elliptic solutions in Figure 4 and Figure 5. The notch depth d is given by the minimum of the square amplitude profile for each value of c. Here, for $c=0.01$, $n=0.830840$ denotes $n=0.830840$.

c	n	d	k	$\mathit{\mu }/2\mathit{nc}$	${\mathit{v}}_{\mathit{c}}/2\sqrt{\mathit{nc}}$	$\mathit{v}/{\mathit{v}}_{\mathit{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

. The estimates of the velocity parameter v are given by making use of Equation (55) with $J=1$. For the strong coupling cases of $c=10$ and 100, the velocity parameter v becomes negative, which will be explained with Equations (77) and (78) in Section 3.3.4.

In Figure 6, the density profiles ${[\langle X\left(q\right),N|{\widehat{\psi }}^{†}\left(x\right)\widehat{\psi }\left(x\right)|X\left(q\right),N\rangle ]}_{J=1}$ of the novel quantum state (31) 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 \left(x\right)$ of the classical solution does not depend on the winding number J. In addition, the ripples around the notch in Figure 6 for $c=100$ are already reported in earlier studies.

The parameters for the fitted curves in Figure 6 are listed in

Table 2. Parameters of the density profiles of the elliptic dark-soliton solutions in Figure 6.

c	d	k	$\mathit{\mu }/2\mathit{nc}$	${\mathit{v}}_{\mathit{c}}/2\sqrt{\mathit{nc}}$	$\mathit{v}/{\mathit{v}}_{\mathit{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

. The estimates of the velocity parameter v are given by making use of Equation (55) with $J=1$. Similarly as in Table 1, for the strong coupling cases of $c=10$ and 100, the velocity parameter v is negative. It will be explained with Equations (77) and (78) in Section 3.3.4.

Finally, we note that quantum states of a single dark soliton with a zero winding number, i.e., $J=0$, have been studied in the previous study [11]. In all the phase profiles, the argument $arg\left[{\psi }_Q\left(x\right)\right]$ satisfies the periodic boundary condition: $arg\left[{\psi }_Q\left(x\right)\right]=arg\left[{\psi }_Q(x+L)\right]$.

3.2. Construction of an Elliptic Multiple Dark Soliton for the GP Equation

3.2.1. Local Density $\rho \left(x\right)$ and the Phase Field $\phi \left(x\right)$ 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)=\varphi (x-vt).$$

(34)

Substituting it into Equation (1), we have

$${\varphi }^{″}-iv{\varphi }^{\prime }+\mu \varphi -2c{\left|\varphi \right|}^2\varphi =0.$$

(35)

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

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

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

(36)

By putting the ansatz (36) into Equation (35), 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\phi }{dx}=\frac{v}{2}+\frac{W}{\rho \left(x\right)},$$

(37)

and

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

(38)

where

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

(39)

Here, W and V are constants with respect to the first integrals (37) and (38), respectively.

Let $a_1$, $a_2$ and $a_3$ denote the three roots of the equation: $U\left(\rho \right)=0$. We assume that they are real. We put them in increasing order: $a_1\le a_2\le a_3$. It follows from the first integral (38) that they are non-negative. The coefficient of ${\rho }^3$ in the potential $U\left(\rho \right)$ of Equation (39) is negative since c is positive: $c>0$, and hence the potential $U\left(\rho \right)$ is negative or zero only in the regions $a_1\le \rho \le a_2$ and $a_3\le \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\le \rho \le a_2$ and $a_1$ is non-negative: $a_1\ge 0$. In fact, the potential $U\left(\rho \right)$ is non-negative at the origin: $U(\rho =0)=W^2\ge 0$. It is straightforward to show the following:

$$W=\pm \sqrt{c{a}_1{a}_2{a}_3}.$$

(40)

By solving Equation (38), we express the local density $\rho \left(z\right)$ in terms of the Jacobi’s elliptic function $\mathrm{sn}(u,k)$. Integrating the differential $2dx=d\rho /\sqrt{-U\left(\rho \right)}$ 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 \left(x\right)=a_1-(a_1-a_2){\mathrm{sn}}^2\left(\sqrt{c(a_3-a_1)}x,k\right),$$

(41)

where k denotes the elliptic modulus

$$k=\sqrt{\frac{a_3-a_1}{a_2-a_1}}.$$

(42)

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,$$

(43)

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 \left(x\right)dx=n.$$

(44)

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

$$\rho \left(x\right)=a_3-(a_3-a_1){\mathrm{dn}}^2(2mKx/L,k).$$

(45)

It follows from the definition of the complete elliptic integral of the second kind: $E={\int }_0^K{\mathrm{dn}}^2(u,k)du$ and Equation (44) that we have

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

(46)

Let us introduce a parameter ${\beta }_k$ by

$${\beta }_k=\frac{1}{nc}{\left(\frac{2mK}{L}\right)}^2.$$

(47)

It follows from Equations (42) and (44) that the three roots are expressed as

$$\begin{array}{ccc}\hfill a_1& =& n\left(1+{\beta }_k\left(\frac{E}{K}-1\right)\right),\hfill \\ \hfill a_2& =& n\left(1+{\beta }_k\left(\frac{E}{K}-1+k^2\right)\right),\hfill \\ \hfill a_3& =& n\left(1+{\beta }_k\frac{E}{K}\right).\hfill \end{array}$$

(48)

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\le 1/(1-E/K)$ since $a_1$ is non-negative: $a_1\ge 0$.

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

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

(49)

Here, we note that the second term on the right-hand side of Equation (49) 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 Equation (37) with respect to variable x, we have

$$\phi \left(x\right)-\phi \left(0\right)=\frac{v}{2}x+\frac{W}{a_1}{\int }_0^x\frac{dy}{1-n_e{\mathrm{sn}}^2\left(2mKy/L,k\right)},$$

(50)

where $n_e$ is given by

$$n_e=1-a_2/a_1.$$

(51)

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

$$\mathsf{\Pi }(r,\varphi =\mathrm{am}u,k)={\int }_0^u\frac{dv}{1-r{\mathrm{sn}}^2(v,k)},$$

(52)

we have

$$\phi \left(x\right)-\phi \left(0\right)=\frac{v}{2}x+\frac{W}{a_1}\frac{L}{2m}\frac{1}{K}\mathsf{\Pi }\left(n_e,\varphi =\mathrm{am}\frac{2mKx}{L},k\right).$$

(53)

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

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

(54)

Putting Equation (53) into it, we express the velocity v in terms of elliptic functions as follows:

$$v=-\frac{2W}{a_1}\frac{1}{K}\mathsf{\Pi }(n_e,k)+\frac{4\pi }{L}J.$$

(55)

We now show that the phase field $\phi \left(x\right)$ is expressed in terms of Jacobi’s Theta function. Let $\mathsf{\Theta }\left(u\right)$ and $Z\left(u\right)$ 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 [39]:

$${\int }_0^u\frac{k^2\mathrm{sn}(a,k)\mathrm{cn}(a,k)\mathrm{dn}(a,k){\mathrm{sn}}^2(u,k)}{1-k^2{\mathrm{sn}}^2(a,k){\mathrm{sn}}^2(u,k)}du=\frac{1}{2}log\left(\frac{\mathsf{\Theta }(u-a)}{\mathsf{\Theta }(u+a)}\right)+uZ\left(a\right),$$

(56)

we derive the following formula:

$$\begin{array}{ccc}\hfill {\int }_0^u\frac{dv}{1-r{\mathrm{sn}}^2(v,k)}& =& u\left(1+Z\left(a\right)\sqrt{\frac{r}{(k^2-r)(1-r)}}\right)+\hfill \\ & & +\frac{1}{2}log\left(\frac{\mathsf{\Theta }(u-a)}{\mathsf{\Theta }(u+a)}\right)\sqrt{\frac{r}{(k^2-r)(1-r)}},\hfill \end{array}$$

(57)

where r is any real number while the parameter a is related to r by $k^2{\mathrm{sn}}^2(a,k)=r$. In order to apply Equation (57) 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:

$$\mathrm{sn}(i\alpha ,k)=i\sqrt{\frac{{\beta }_k}{1-{\beta }_k(1-E/K)}}.$$

(58)

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

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 $\phi \left(x\right)$ as the logarithm of a ratio of Theta functions

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

(59)

for $W=\pm \sqrt{c{a}_1{a}_2{a}_3}$. Here, the plus/minus signs are in correspondence. In the derivation of Equation (59), through Equation (57), we have expressed the complete elliptic integral of the third kind in terms of the Zeta function

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

(60)

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).$$

(61)

Here, we note that velocity v is given by Equation (55) 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 Equation (61). Here, we recall that velocity v is evaluated by Equation (55).

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 Equation (61) 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\left(k\right)=K(3E-(2-k^2)K)$ is monotonically decreasing with respect to modulus k. The largest value of $f\left(k\right)$ is therefore given by $f\left(0\right)={\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 )}.$$

(62)

We call it the critical velocity.

3.3. Reduction of the 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\left(k\right)/L$ is kept constant. We define the ratio b by $b=K/L$. Here, we note 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 =\underset{k\to 1,L\to \infty }{lim}{\beta }_k=\frac{{\left(2mb\right)}^2}{cn}.$$

(63)

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

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

(64)

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,\dots ,m-1$, in the interval $[0,L]$. Near the position of the jth 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.$$

(65)

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

$$\rho \left(x\right)/n=1-\beta {\sec \mathrm{h}}^2(2b\Delta x)+2\beta /log(1/p)+\cdots .$$

(66)

It is easy to show that the asymptotic expansion of the square amplitude up to the order of $1/logp$, 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 \left(x\right)=\left(\frac{{\mu }_{\infty }}{2c}\right)\left\{1-\beta {\mathrm{sech}}^2\left[{\left(\frac{\beta {\mu }_{\infty }}{2}\right)}^{1/2}(x-x_j)\right]\right\}.$$

(67)

Here, $\beta$ is expressed as

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

(68)

where $v_{c,\infty }=\sqrt{2{\mu }_{\infty }}$. Here, we note that we have derived the square amplitude of the single dark-soliton due to Tsuzuki [3] 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 jth 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:

$$\begin{array}{ccc}\hfill \phi \left(x\right)-\phi \left(0\right)& =& \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)\hfill \\ & & \pm \frac{1}{2i}log\left(1+i\sqrt{\frac{\beta }{1-\beta }}tanh\left(\sqrt{\frac{\mu }{2}\beta }\Delta x\right)\right)+O(1/log(1/p))\hfill \end{array}$$

(69)

for $W=\pm \sqrt{c{a}_1{a}_2{a}_3}$. Here, the plus/minus signs are in correspondence.

For every j-th notch, we have the expression (67) 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 jth notch at time $t=0$. It follows from Equation (69) that, for very small $|\Delta x|$, the phase field $\phi \left(x\right)$ is approximately given by the following if the system size L is very large

$$\phi \left(x\right)-\phi \left(0\right)=\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\left(\right|\Delta x\left|\right)\left(\right|\Delta x|\ll 1).$$

(70)

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

$$\phi \left(x\right)-\phi \left(0\right)-\frac{2\pi J}{L}{x}_j\approx \left(\pm \sqrt{\frac{\mu (1-\beta )}{2}}·\frac{\beta }{1-\beta }+\frac{2\pi J}{L}\right)\Delta x\left(\right|\Delta x|\ll 1)$$

(71)

for $W=\pm \sqrt{c{a}_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\phi }{dx}\right|}_{\Delta x=0}=\pm \sqrt{\frac{\mu (1-\beta )}{2}}\frac{\beta }{1-\beta }+\frac{2\pi J}{L}$$

(72)

for $W=\pm \sqrt{c{a}_1{a}_2{a}_3}$. Thus, the phase field $\phi \left(z\right)$ may abruptly increase ($W>0$) or decrease ($W<0$) near the jth notch $x_j=jL/m$ if $\left|J\right|$ 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 Equation (55) or as we shall show in Equation (73). The phase field abruptly decreases at the jth notch, i.e., we have the minus sign in Equation (71).

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 Section 4.

3.3.3. Evaluation of Velocity in an Infinite Limit

We first recall that, through Equation (60), 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},$$

(73)

where W is given by $W=\pm \sqrt{c{a}_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))(p\to 0).$$

(74)

Putting it into Equation (73), we have

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

(75)

for $W=\pm \sqrt{c{a}_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\phi }{dx}\right|}_{\Delta x=0}=\frac{v}{2}·\frac{\beta }{1-\beta }.$$

(76)

It is larger than the value of velocity v if $\beta$ is close to 1. Here, we note 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{c{a}_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}.$$

(77)

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

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

(78)

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 jth notch of an elliptic m-fold multiple dark-soliton is approximately given by

$$\phi \left(x\right)-\phi \left(0\right)-\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\left(\right|\Delta x|\ll 1)$$

(79)

for $W=\pm \sqrt{c{a}_1{a}_2{a}_3}$.

It follows from Equations (77) and (79) that, if $\beta$ is close to 1, it happens very rarely that both the phase gradient $d\phi /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 Equation (77) is positive and hence v is positive, while that of Equation (79) 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 Equation (77) is negative and hence v is positive only if the inequality in Equation (78) holds, while that of Equation (79) is positive and hence the phase gradient is positive. Thus, for $W>0$, both the phase gradient $d\phi /dx$ and the velocity v are positive simultaneously if the inequality (78) 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 [40]. First, we construct the quantum state in the free fermions which possesses one notch in its density profile shown in Figure 7a. 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 }^{†}$ ($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 }^{†}$ satisfy the anti-commutation relations

$$\begin{array}{ccc}\hfill & \{a_{\ell },a_m^{†}\}={\delta }_{\ell ,m},\hfill & \hfill \{a_{\ell },a_m\}=\{a_{\ell }^{†},a_m^{†}\}=0.\end{array}$$

(80)

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

$$\begin{array}{ccc}\hfill & \psi \left(x\right)=\frac{1}{\sqrt{L}}\sum _{\alpha =1}^{\infty }{\mathrm{e}}^{\mathrm{i}{k}_{\alpha }x}{a}_{\alpha },\hfill & \hfill {\psi }^{†}\left(x\right)=\frac{1}{\sqrt{L}}\sum _{\alpha =1}^{\infty }{\mathrm{e}}^{-\mathrm{i}{k}_{\alpha }x}{a}_{\alpha }^{†}.\end{array}$$

(81)

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

$$\begin{array}{cc}\hfill |{\mathsf{\Phi }}_1\rangle & =\psi \left(X_1\right)|gs\rangle =\frac{1}{\sqrt{{\mathcal{N}}_1}}\sum _{\alpha =1}^N{(-1)}^{\alpha }{e}^{i{k}_{\alpha }{X}_1}\prod _{\ell =1;\ell \ne \alpha }^N{a}_{\ell }^{†}|0\rangle .\hfill \end{array}$$

(82)

It is obtained by applying a field operator $\psi \left(X_1\right)$ at the position $x=X_1$ to the N-particle ground state $|gs\rangle$ defined by

$$\begin{array}{c}\hfill |gs\rangle \equiv \prod _{\ell =1}^N{a}_{\ell }^{†}|0\rangle .\end{array}$$

(83)

Here, ${\mathcal{N}}_1$ in Equation (82) is the normalization factor. The state $|{\mathsf{\Phi }}_1\rangle$ 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^{-i{k}_{\alpha }{X}_1}$ in Equation (82) has a position $X_1$, implying the position of the notch in the expectation value of the local density operator for its state [11].

Generalizing the state $|{\mathsf{\Phi }}_1\rangle$, 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,

$$\begin{array}{c}\hfill |{\mathsf{\Phi }}_2\rangle =\psi \left(X_2\right)\psi \left(X_1\right)|gs\rangle =\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 \ne \alpha ,\beta }^N{a}_{\ell }^{†}|0\rangle .\end{array}$$

(84)

The state $|{\mathsf{\Phi }}_2\rangle$ is constructed by the double summation of the N-particle states which are given by removing two particles of 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 7b shows the density profile for the state (84). It has two notches at the positions $x=X_1$ and $X_2$.

The ripples appear in Figure 7. This oscillating behavior in fermions has been explained with an approximate delta function called the Dirichlet kernel [11].

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

$$\begin{array}{c}\hfill |X_1,X_2,N\rangle =\frac{1}{\sqrt{\mathcal{M}}}\sum _{\mathit{p}\in \mathit{P}}{e}^{i(p_1{X}_1+p_2{X}_2)}|p_1,p_2,N\rangle \end{array}$$

(85)

with the normalization constant $\mathcal{M}$, and the two holes’ positions $\mathit{p}:=\{p_1,p_2\}$ in the set $\mathit{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 $\mathit{P}=\left\{\{-2,-1\},\{-2,0\},\{-2,1\},\{-2,2\},\cdots ,\{2,4\},\{3,4\}\right\}$, which has $\left|\mathit{P}\right|=\left(\genfrac{}{}{0pt}{}{7}{2}\right)=21$ elements. The quantum double dark soliton state is a superposition of $\left|\mathit{P}\right|=\left(\genfrac{}{}{0pt}{}{N+2}{2}\right)$ two-hole excitation states.

With respect to the state (85), 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

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\rho }_Q\left(x\right):=& \langle X_1,X_2,N\left|\widehat{\rho }\left(x\right)\right|X_1,X_2,N\rangle \hfill \\ \hfill =& \frac{1}{\mathcal{M}}\sum _{\mathit{p},{\mathit{p}}^{\prime }\in \mathit{P}}{e}^{i(P-P^{\prime })x}{e}^{i(p_1{X}_1+p_2{X}_2)}{e}^{-i(p_1^{\prime }{X}_1+p_2^{\prime }{X}_2)}\langle p_1^{\prime },p_2^{\prime },N\left|\widehat{\rho }(0,0)\right|p_1,p_2,N\rangle ,\hfill \end{array}\end{array}$$

(86)

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

$$\begin{array}{c}\hfill P=\frac{2\pi }{L}\sum _{m\in I}m,I\equiv I^{+}-\mathit{p}.\end{array}$$

(87)

Again, the form factors $\langle P_1,P_2,N\left|\widehat{\rho }(0,0)\right|P_1^{\prime },P_2^{\prime },N\rangle$ 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 $|X\left(q\right),N\rangle$ in Section 2.1 and $|N,X;P\rangle$ 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\left(x\right)$ of the quantum double dark solitons in the weak coupling regime agree with those of the corresponding elliptic multiple dark solitons (49), as shown in Figure 9. It shows the quantum double dark solitons (86) and their corresponding elliptic soliton solutions with different coupling constants. Although the system is in the strong coupling regime for $c=10$ and 100, we have displayed these plots for comparison. 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. It is the case also for quantum single dark solitons. We suggest that this undulating line is characteristic to fermions, as explained in Section 4.1.

By calculating the time dependent field operator $\widehat{\psi }(x,t)$, we obtain the density profile evolving in time

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \langle X_1,X_2,N\left|\widehat{\rho }(x,t)\right|X_1,X_2,N\rangle =& \frac{1}{\mathcal{M}}\sum _{\mathit{p},{\mathit{p}}^{\prime }\in \mathit{P}}{e}^{i(P-P^{\prime })x}{e}^{i(p_1{X}_1+p_2{X}_2)}{e}^{-i(p_1^{\prime }{X}_1+p_2^{\prime }{X}_2)}{e}^{\left[-\mathrm{i}\left(E-E^{\prime }\right)t\right]}\hfill \\ \hfill & \langle p_1^{\prime },p_2^{\prime },N\left|\widehat{\rho }(0,0)\right|p_1,p_2,N\rangle ,\hfill \end{array}\end{array}$$

(88)

where E is the energy of the state $|p_1,p_2,N\rangle$, and $\widehat{\rho }(x,t)={\widehat{\psi }}^{†}(x,t)\widehat{\psi }(x,t)$. Though the positions of the notches are located at equal distances from each other in this paper, we have observed the time evolution in the case where the notches’ positions are located at different distances such as $X_1=0$ and $X_2=L/4$. In this case, each notch moves with different speeds although their depths are the same.

4.2.3. Matrix Elements of the Quantum Double Dark Soliton State

We have calculated the matrix element by using the determinant formula in Equation (18). The matrix element ${\psi }_Q\left(x\right)$ is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\psi }_Q\left(x\right):=\langle X_1,X_2,N-1|\widehat{\psi }\left(x\right)|X_1,X_2,N\rangle \hfill \\ \hfill =& \frac{1}{\mathcal{M}}\sum _{\mathit{p},{\mathit{p}}^{\prime }\in \mathit{P}}{e}^{i(P-P^{\prime })x}{e}^{i(p_1{X}_1+p_2{X}_2)}{e}^{-i(p_1^{\prime }{X}_1+p_2^{\prime }{X}_2)}\langle P_1,P_2,N-1|\widehat{\psi }(0,0)|P_1^{\prime },P_2^{\prime },N\rangle ,\hfill \end{array}\end{array}$$

(89)

which takes a complex value for position x. The classical counterpart of the matrix element ${\psi }_Q\left(x\right)$ is the relevant elliptic solution, which is given by Equation (36), with its absolute part and argument. Comparisons of the matrix element ${\psi }_Q\left(x\right)$ and the classical soliton ${\psi }_C\left(x\right)$ 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{\left(x\right)|}^2$ for the matrix element (89). The smallest values of the square amplitude of the notches are smaller than those of the density profiles in Figure 9. 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. In the strong coupling regime, such as $c=100$, the line connecting the two notches undulates in Figure 10 as we have already observed in Figure 9. We suggest that this undulating line is characteristic to fermions, as explained in Section 4.1.

Figure 11 shows the phase profile $arg\left[{\psi }_Q\left(x\right)\right]$ for the matrix element (89). 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\left[{\psi }_Q\left(x\right)\right]$ is satisfied with the periodic boundary condition such as $arg\left[{\psi }_Q(x+L)\right]=arg\left[{\psi }_Q\left(x\right)\right]$ 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 [11].

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$.

4.3. Quantum State 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 Equation (85), quantum states for double dark-solitons are quite effectively constructed. The method is similar to that of constructing quantum single dark soliton states [14,15]: We assign the ideal Gaussian weights [15] to the sum of the eigenstates with two particle-hole excitations given in Equation (85). Then, the constructed quantum states give such density profiles that have distinct notches with different depths.

We employ the ideal Gaussian weights [15] 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 $\mathit{p}=\{p_1,p_2\}$ in the set $\mathit{P}$:

$$\begin{array}{c}\hfill {|X_1,X_2,N\rangle }_G=\frac{1}{\sqrt{\mathcal{N}}}\sum _{\{p_1,p_2\}\in \mathit{P}}{G}_{P_0}\left(p_1\right)G_{P_0^{\prime }}\left(p_2\right)e^{i(p_1{X}_1+p_2{X}_2)}\left|p_1,p_2,N\right.〉.\end{array}$$

(90)

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

$$\begin{array}{c}\hfill G_P\left(q\right)=exp\left[-\frac{{\left(q-P\right)}^2}{4{\sigma }^2}\right]\end{array}$$

(91)

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 [15] and defined by the following:

$$\begin{array}{c}\hfill \sqrt{d}=\left|{\varphi }_P^{\infty }\left(0\right)\right|=\sqrt{n}\left|\frac{v}{v_{c,\infty }}\right|.\end{array}$$

(92)

Here, $|{\varphi }_P^{\infty }\left(0\right)|$ denotes the square root of the local density $\rho \left(x\right)$ at the origin, i.e., at the position of a notch, in the thermodynamic limit. We recall that the local density $\rho \left(x\right)$ near a notch has been given in Equation (67). 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

$$\begin{array}{cc}\hfill P\left(d\right)& =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\},\hfill \end{array}$$

(93)

$$\begin{array}{cc}\hfill {\sigma }^2& =\frac{4}{3}n\sqrt{nc}{\left(1-\frac{d}{n}\right)}^{\frac{3}{2}}.\hfill \end{array}$$

(94)

4.3.2. Profiles of Density ${\rho }_Q\left(x\right)$ and Square Amplitude $|{\psi }_Q{\left(x\right)|}^2$

For the quantum state of Equation (90), we exactly compute the density ${\rho }_Q\left(x\right)$ and the matrix element ${\psi }_Q\left(x\right)$ and plot the profiles of the density ${\rho }_Q\left(x\right)$, the square amplitude $|{\psi }_Q{\left(x\right)|}^2$, and the phase $arg\left[{\psi }_Q\left(x\right)\right]$ in Figure 12, Figure 13 and Figure 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 Equation (90) are given by $(P_0,\sigma )=(0.124027\pi ,0.106667)$ and $(P_0^{\prime },{\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\le 0.1$. In the two panels of $c\le 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\le 0.1$.

We thus suggest from Figure 12 and Figure 13 that the Gaussian superposition of the eigenstates with two particle-hole excitations in Equation (90) 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\le 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\le 0.1$ suggests that the quantum state with the ideal Gaussian weights in Equation (90) should be close to the corresponding mean-field product state. In the panels of $c=0.05$ and $0.1$, every notch in Figure 12 and Figure 13 is distinct and narrow. Furthermore, the value of ${\rho }_Q\left(x\right)$ or $|{\psi }_Q{\left(x\right)|}^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 Figure 12 and Figure 13.

Here, we note that, if the Gaussian superposition of the eigenstates in the Lieb–Liniger model of Equation (90) completely produce the mean-field product state satisfying Equation (22), then the square amplitude of the matrix element ${\psi }_Q\left(x\right)={}_G{\langle X_1,X_2,N-1|\widehat{\psi }\left(x\right)|X_1,X_2,N\rangle }_G$ is equal to the density profile ${\rho }_Q\left(x\right)={}_G{\langle X_1,X_2,N|{\widehat{\psi }}^{†}\left(x\right)\widehat{\psi }\left(x\right)|X_1,X_2,N\rangle }_G$: We have ${\rho }_Q\left(x\right)={\left|{\psi }_Q\left(x\right)\right|}^2$.

4.3.3. Difference between Density ${\rho }_Q\left(x\right)$ and Square Amplitude $|{\psi }_Q{\left(x\right)|}^2$

We now estimate how far the quantum soliton state constructed by Equation (90) satisfies the characteristic property of Equation (22) for the mean-field product state. Let us subtract the square amplitude of the matrix element, $|{\psi }_Q{\left(x\right)|}^2$, from the density, ${\rho }_Q\left(x\right)$. Inserting the identity operator between ${\widehat{\psi }}^{†}\left(x\right)$ and $\widehat{\psi }\left(x\right)$ in the density operator, we obtain the square amplitude of the matrix element and other terms:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {}_G{\langle X_1,X_2,N|{\widehat{\psi }}^{†}\left(x\right)\widehat{\psi }\left(x\right)|X_1,X_2,N\rangle }_G\hfill \\ \hfill =& \sum _{\chi }{}_G{\langle X_1,X_2,N|{\widehat{\psi }}^{†}\left(x\right)|\chi \rangle \langle \chi |\widehat{\psi }\left(x\right)|X_1,X_2,N\rangle }_G\hfill \\ \hfill =& |{}_G{\langle X_1,X_2,N-1|\widehat{\psi }\left(x\right)|X_1,X_2,N\rangle }_G{|}^2+{|{\langle {\chi }^{\prime }|\widehat{\psi }\left(x\right)|X_1,X_2,N\rangle }_G|}^2+\cdots .\hfill \end{array}\end{array}$$

(95)

We define the difference $D\left(x\right)$ by

$$\begin{array}{cc}\hfill D\left(x\right):=& {}_G{\langle X_1,X_2,N|{\widehat{\psi }}^{†}\left(x\right)\widehat{\psi }\left(x\right)|X_1,X_2,N\rangle }_G-{|{}_G{\langle X_1,X_2,N-1|\widehat{\psi }\left(x\right)|X_1,X_2,N\rangle }_G|}^2.\hfill \end{array}$$

(96)

Since each term on the right-hand side of Equation (95) is non-negative, we have the relation

$$D\left(x\right)\ge 0.$$

(97)

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

We now argue that the quantum double dark soliton state of Equation (90) should be close to the corresponding mean-field product state in the weak coupling regime such as $c\le 0.1$. Let us evaluate $D\left(x\right)$ 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, Figure 13, Figure 14 and Figure 15 (see also

Table 5. The values of the density profile ${\rho }_Q\left(x\right)$, those of the square amplitudes $|{\psi }_Q{\left(x\right)|}^2$, and their differences $D\left(x\right)$ at $x=0$ and $x=5$.

c	$\mathit{x}=0$			$\mathit{x}=5$
	${\mathit{\rho }}_{\mathit{Q}}\left(\mathit{x}\right)$	$|{\mathit{\psi }}_{\mathit{Q}}{\left(\mathit{x}\right)|}^{\mathbf{2}}$	$\mathit{D}\left(\mathit{x}\right)$	${\mathit{\rho }}_{\mathit{Q}}\left(\mathit{x}\right)$	$|{\mathit{\psi }}_{\mathit{Q}}{\left(\mathit{x}\right)|}^{\mathbf{2}}$	$\mathit{D}\left(\mathit{x}\right)$
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

). On the other hand, for a strong coupling constant $c=10$, the difference $D\left(x\right)$ is about 0.6 except for the notches.

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

In Figure 15, density profiles (red curves) and square amplitude of matrix element (blue curves) are simultaneously plotted. The deviations of the two plots are similar to those of quantum single dark solitons in the weak coupling regime of $c=0.1$ as shown in Figure 3.

By explicitly calculating the time dependent field operator $\widehat{\psi }(x,t)$, we obtain the density profile evolving in time. All figures in this section correspond to the profiles at $t=0$. In the time evolution of the density profile for $c=0.1$ with $N=L=20$ such as in the upper right panel of Figure 12, we observe the following. First, as the solitons approach each other, they move along approximately straight-line trajectories with constant speeds. The collision occurs around $t=15$. After the collision, the profiles of solitons appear again with the same constant speeds as before the collision. We can estimate the phase shift after the collision, although the result should be described in future reports.

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\left[{\psi }_g(x+L)\right]=arg\left[{\psi }_g\left(x\right)\right]+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.

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 16, 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 Equation (71) 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 $\phi \left(x\right)$ is expressed in terms of the elliptic integral of the third kind together with velocity v and parameter W in Equation (53) or more compactly in Equation (59) in terms of the elliptic theta function. The parameter v is determined by the boundary condition (54) with the winding number J. For the real parameter W, we have the choice of signs, as shown in Equation (40).

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 17. 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 17, 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 17, 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 17. 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.

5. Conclusions

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\left(x\right)$. 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. [15] 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., Equation (22), 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 the 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. [12,41]. They showed quantum single/double dark-solitons by measuring particle positions in a single one-/two-hole excitation state. Even though their method is different from ours in Section 2.1 and Section 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 [15] and the particle position method [12] in ref. [15]. 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 recent research [16].

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 a 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 equation, obtained by the inverse scattering method does not satisfy the periodic boundary conditions. The finite-gap soliton [20] 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 method in the paper is directly extended to the study of time evolution as suggested in Equation (88). We note that some aspects of the beyond mean-field dynamics of quantum dark solitons are investigated [13,42].