Quantum states of dark solitons in the 1D Bose gas

We present a series of quantum states that are characterized by dark solitons of the nonlinear Schr\"{o}dinger equation (i.e. the Gross-Pitaevskii equation) for the one-dimensional (1D) Bose gas interacting through the repulsive delta-function potentials. The classical solutions satisfy the periodic boundary conditions and we call them periodic dark solitons. Through exact solutions we show corresponding aspects between the states and the solitons in the weak coupling case: the quantum and classical density profiles completely overlap with each other not only at an initial time but also at later times over a long period of time, and they move together with the same speed in time; the matrix element of the bosonic field operator between the quantum states has exactly the same profiles of the square amplitude and the phase as the classical complex scalar field of a periodic dark soliton not only at the initial time but also at later times, and the corresponding profiles move together for a long period of time. We suggest that the corresponding properties hold rigorously in the weak coupling limit. Furthermore, we argue that the lifetime of the dark soliton-like density profile in the quantum state becomes infinitely long as the coupling constant approaches zero, by comparing it with the quantum speed limit time. Thus, we call the quantum states quantum dark soliton states.


Introduction
The experimental technique of trapped one-dimensional atomic gases [1,2] has now become a fundamental tool for revealing nontrivial effects in quantum many-body systems [3,4]. For the interacting Bose gas in one dimension (1D), the first set of exact results goes back to the pioneering work of Girardeau [5] on the impenetrable Bose gas where the strong interacting limit is considered. The 1D Bose gas interacting with the delta-function potentials, i.e. the Lieb-Liniger (LL) model, gives a solvable model for interacting bosons in 1D [6], where it is integrable even when the interaction parameter is generic. For the impenetrable Bose gas which corresponds to the Tonks-Girardeau (TG) limit, i.e. the strong interacting limit of the LL model, the one-body reduced density matrix is derived and successfully expressed in terms of the determinant of a Fredholm operator [7]. The exact result is followed by several important developments in mathematical physics [8,9]. Furthermore, dynamical correlation functions of the LL model are now systematically derived [10].
Soliton-like localized excitations in a periodic 1D Bose gas have attracted much interest recently and have been studied theoretically [11,12,13]. Here we remark that dark solitons are created experimentally in cold atoms [14], for instance, by the phase-imprinting method [15,16] (See also [17]). Localized quantum states are important for investigating dynamical responses of interacting quantum systems. Quantum dark solitons in confining potentials are studied by semiclassical quantization [18], and those in thermal equilibrium of a quasi-1D Bose gas by generating classical field ensembles [19]. However, it is not clear even at zero temperature how we can construct quantum states associated with dark solitons in the many-body system of the LL model.
Let us consider the Gross-Pitaevskii (GP) equation, which describes Bose-Einstein condensation (BEC) in the mean-field approximation [20]. We also call it the nonlinear Schrödinger equation. The GP equation has dark soliton solutions for the repulsive interactions, while it has bright soliton solutions for the attractive interactions [21]. It was conjectured that dark solitons are closely related to Lieb's type-II excitations, i.e. onehole excitations, by carefully studying the dispersion relations [22]. The dispersion relations of the LL model are briefly compared with those of the classical nonlinear Schrödinger equation in the weak coupling limit [23]. However, it has not been shown how one can construct such a quantum state that leads to a dark soliton in the classical limit or what kind of physical quantity can show a property of a dark soliton for some quantum state.
Here we remark that each of the type-II eigenstates has a flat density profile since the Bethe ansatz eigenstates are translationally invariant. Moreover, we consider dark solitons under the periodic boundary conditions, which are expressed in terms of elliptic functions [11,12,13].
In this paper we demonstrate that a quantum state derived from the Bethe ansatz eigenvectors of the LL model by taking the Fourier transform of the type-II excitations over a branch [24] has many aspects closely related to classical dark solitons of the GP equation under the periodic boundary conditions. We call the state the quantum dark soliton state and a classical dark soliton under the periodic boundary conditions simply a classical dark soliton. Through the Bethe ansatz we show several corresponding aspects in the weak coupling regime. Firstly, the density profile of the quantum dark soliton state at an initial time is consistent with that of a classical dark soliton. Here we call the graph of the expectation value of the local density operator for a given state versus the position coordinate the density profile of the state, and for a quantum dark soliton state we simply call it the density profile of quantum dark soliton; we call the graphs of the square amplitude and phase in the complex scalar field of a classical dark soliton versus the position coordinate the density profile and phase profile of classical dark soliton, respectively. Secondly, in time evolution the density profile of quantum dark soliton coincides with that of the corresponding classical dark soliton over the whole graph and they move together with the same velocity for a long period of time. Thirdly, for the matrix element of the bosonic field operator between two quantum dark soliton states where one state has N − 1 particles and another N particles, the profiles of the square amplitude and phase at the initial time graphically agree with those of classical dark soliton, respectively. In time evolution the profiles of square amplitude and phase overlap with those of classical dark soliton, respectively, over the whole region and the corresponding profiles move together in time for a long period of time. Here we remark that a classical dark soliton parametrized by elliptic functions becomes a standard dark soliton with open boundaries by simultaneously sending the system size to infinity and the elliptic modulus to unity. Furthermore, in order to illustrate the method for constructing quantum dark solitons, in the 1D free fermions we show from the anti-commutation relations that a notch appears in the density profile of some superposition of one-hole excitations. Interestingly, the density profile of the fermionic state coincides with that of quantum dark soliton for the 1D Bose gas in the large coupling case, i.e. near the TG limit, not only at an initial time but also during the quantum dynamics for some period of time.
The time evolution of the expectation value of the local density operator in the 1D Bose gas should be important also from the renewed interest on the fundamentals of quantum statistical mechanics [25,26,27,28,29]. The density profile of a quantum dark soliton state has initially a localized notch but collapses slowly to a flat profile in time evolution [24]. The relaxation behavior is consistent with the viewpoints of equilibration of an isolated quantum system [30,31] and thermalization due to the typicality of quantum states [32].
We now argue that the density profile of quantum dark soliton has a finite lifetime but it is much longer than the quantum speed limit time of the quantum state in the weak coupling case. Here we remark that the lifetime of a generic state is given by the quantum speed limit time for it [33,34,35]. By observing the exact time evolution of the density profile of quantum dark soliton we estimate its lifetime. We shall show in section 5 that the observed life time is inversely proportional to the coupling constant. We thus suggest that the localized density profile of a quantum dark soliton is much more stable than the density profile of a generic quantum state in the weak coupling case. Here we remark that Girardeau and Wright discussed a permanent quantum soliton for impenetrable bosons in 1D [36], which corresponds to the infinite coupling case of the LL model.
We also argue that the behavior of the wavefunctions is nontrivial in the weak coupling limit for the LL model. The wavefunctions do not simply become close to such wavefunctions of the 1D free bosons that could be consistent with the mean-field picture, when the coupling constant approaches zero but it takes a nonzero value. The exact Bethe ansatz wavefunctions consist of a large number of terms such as N !. Furthermore, we can show that many-body correlations increase if we increase the particle number N while we keep the coupling constant small but fixed and the particle density constant. For instance, the zero mode fraction in the ground state of the 1D Bose gas, which we regard as the condensate fraction of BEC, becomes small and decreases to zero if the number of particles increases while the coupling constant and the particle density are fixed [37]. We suggest that it is also the case for quantum dark solitons i.e., that some properties of quantum dark solitons corresponding to classical dark solitons do not necessarily hold if we increase the particle number N while we keep the coupling constant small but fixed and the density of particles fixed. In fact, it is not a priori clear how valid the mean-field approximation is for quantum dark solitons even in the weak coupling case. Here we remark that the breakdown of mean-field theory is addressed for impenetrable bosons [38], and a long-wavelength theory beyond the mean-field approximation is discussed for low-dimensional bose liquids [39].
Let us consider the Hamiltonian of the LL model for N bosons with coupling constant c Here we impose periodic boundary conditions with length L. We employ a system of units with 2m = = 1, where m is the particle mass. We introduce the canonical Bose field operatorψ(x, t) with the commutation relations [9] [ψ( The second-quantized Hamiltonian of the LL model is written in terms of the field operatorψ(x, t) as where µ is the chemical potential. The Heisenberg equation of motion of this system has the form In the classical limit where the quantum field operatorψ(x, t) is replaced by a complex c-number field ψ C (x, t), equation (5) becomes the following partial differential equation We call it the nonlinear Schrödinger equation. The equation (6) can be solved by the inverse scattering method and has soliton solutions [21] (see also [40]). We recall that it has dark soliton solutions for the repulsive interactions c > 0, while it has bright soliton solutions for the attractive interactions c < 0. For the attractive case, bright solitons are analytically derived from the quantum wave packets constructed through the Fourier transform of the bound states [41]. The construction was extended for optical fibers [42]. Time evolution of quantum states of bright solitons was investigated theoretically [43] and experimentally [44]. The effect of quantum noises on the quantum states of bright solitons was also studied [45,46,47]. The contents of the paper consists of the following. In section 2 we briefly introduce the notation for the Bethe ansatz, characterize the type II excitation branch, and then construct quantum dark soliton states. We argue that the structure of quantum dark soliton states is analogous to the Dirichlet kernel, and show for free fermions that the density profile of a superposition of one-hole excited states has a notch. It illustrates the construction of quantum dark soliton states. In section 3 we derive one-soliton solutions of the nonlinear Schrödinger equation (6) under the periodic boundary conditions. We express the classical dark solitons in terms of elliptic functions similarly as in Ref. [11,13]. By sending the elliptic modulus k to 1, the classical dark solitons under the periodic boundary conditions become the standard classical dark solitons of the nonlinear Schrödinger equation with open boundaries. Furthermore, we calculate the logarithmic corrections to them, in particular, to the chemical potential and the soliton velocity, by which we can confirm the numerical estimates of the soliton parameters shown in section 4 at least approximately. In section 4 we show several corresponding aspects between the quantum dark soliton states and the corresponding classical dark solitons, which hold in the weak coupling case and should be rigorously valid when the coupling constant approaches zero with the density and the particle number N fixed. Firstly, we show that the density profile of quantum dark soliton at an initial time is consistent with that of classical dark soliton. Secondly, we show that in exact time evolution the density profile of quantum dark soliton as a whole coincides with that of classical dark soliton for a long period of time.
Thirdly, we show that both the square amplitude and phase profiles for the matrix element of the quantum field operatorψ(x, t) in the LL model between the quantum dark soliton states with N − 1 and N particles overlap with those of classical dark soliton, respectively, all over the profiles. The square amplitude and phase profiles move in time evolution with exactly the same speed as those of classical dark soliton, respectively. Furthermore, we remark that the density profile of a sum of one-hole excitations in the 1D free fermions overlap with that of quantum dark soliton in the 1D Bose gas with the large coupling constant, i.e. in a regime close to the TG limit. For the matrix element of the fermionic field operator between two such fermionic states the square amplitude and phase profiles have several features in common with those of quantum dark soliton in the 1D Bose gas for the large coupling case, respectively, although they do not overlap each other. We thus obtain useful tools to describe approximately the profiles of quantum dark soliton in the large coupling case. However, it is still quite nontrivial that the profiles of quantum dark soliton in the weak coupling case are consistent with those of classical dark soliton. In section 5 we argue that the lifetime of a notch in the density profile of quantum dark soliton is much longer than that of a generic state when the coupling constant c is very small. We compare the lifetime of a notch in the density profile of quantum dark soliton with the quantum speed limit time, and show that the ratio of the former to the latter increases as the coupling constant decreases. It becomes infinitely large as the coupling constant approaches zero. Thus, although the localized density profile collapses in time evolution, we call the quantum states quantum dark soliton states.
2 Excited states of the 1D Bose gas with delta-function potentials

Bethe ansatz eigenwavefunctions
In the framework of the Bethe ansatz method, the Bethe ansatz wavefunction which is expressed in terms of quasi-momenta gives an eigenstate of the LL model (1), if the quasi-momenta satisfy the Bethe ansatz equations In the logarithmic form, we have Here I j 's are integers for odd N and half-odd integers for even N . We call them the Bethe quantum numbers. The total momentum P and energy eigenvalue E are written in terms of the quasi-momenta as In the case of repulsive interaction c > 0, the Bethe equations (9) have a unique real solution k 1 < · · · < k N if we specify a set of Bethe quantum numbers I 1 < · · · < I N [9]. The configuration of the Bethe quantum numbers for the ground state is given by the regular array around the center: I j = −(N + 1)/2 + j for j = 1, 2, . . . , N.
It is depicted in the top panel of Figure 1 (a) in the case of N = 5. We can construct low-lying excited states systematically by putting holes and placing particles in the regular array of the Bethe quantum numbers for the case of the ground state. The arrays of the Bethe quantum numbers Figure 1: (a) Configurations of Bethe quantum numbers for the ground state, 1-hole excitations and 1-particle excitations in the case of N = 5. (b) Dispersion relations of the 1-hole excitations and 1-particle excitations in the case of N = 5. A triangle represents the ground state, big circles represent 1-particle excitations, and squares represent 1-hole excitations. Small dots are the 1-particle 1-hole excitations.
for one-hole excitations are shown in the second panel of Figure 1 (a). In each sequence there exists an unoccupied number in the middle of the sequence of the five red circles, which we call the one-hole. We remark that the one-hole excitations are also called the type II excitations [6]. Similarly as in the second panel, the arrays of the Bethe quantum numbers for one-particle excitations are shown in the third panel of Figure 1 (a). The right most red particle corresponds to the one-particle in each array. The one-particle excitations are also called the type I excitations [6].
The dispersion relation of the one-particle excitations and that of the one-hole excitations are plotted with filled circles and filled squares, respectively, in Figure 1 (b). The former branch corresponds to the type I excitations and the latter branch the type II excitations.

Construction of quantum dark soliton states through type II excitations
We now explain how we construct a series of quantum states [24] through the type II excitations. Furthermore, we call them quantum dark soliton states.
In the type II branch, for each integer p in the set {0, 1, . . . , N − 1}, we consider momentum P = 2πp/L. We denote the normalized Bethe eigenstate of N particles with total momentum P by |P, N , which we call a one-hole excited state. The Bethe quantum numbers I j 's of the one-hole excitation |P, N are given by For any value of X satisfying 0 ≤ X < L we define the coordinate state |X, N by the discrete Fourier transformation: The density profile of this state X, N |ψ † (x)ψ(x)|X, N shows a density notch as will be shown in section 4.

Quantum state analogue of a truncated series of the delta function
The formal structure of quantum dark soliton states (13) is analogous to a truncated series of the delta function. For periodic functions with period L the delta function is expressed in the form of an infnite series We now truncate it by keeping only such terms with integers n satisfying |n| ≤ N c It is called the Dirichlet kernel in the Fourier series analysis. It has the peak around at x = 0 with peak height (2N c + 1)/L and oscillating behavior in the region for |x| > L/(2N c + 1) with period L/(2N c + 1). Here we assume that N c is independent of L. For N c ≫ L, the oscillations in the graph of the Dirichlet kernel become of very short wavelength, and we may call them ripples. The integral of the squared amplitude of the Dirichlet kernel over the period L is given by We therefore normalize the function D Nc (x) by dividing it by the square root ( In terms of the normalized momentum eigenfunctions x|P = exp(2πipx/L)/ √ L where P = 2πp/L the normalized functionD Nc (x) is expressed as It is analogous to the structure of the states of quantum dark solitons (13), where the number 2N c + 1 in (18) corresponds to the number of states N in (13), and the single-particle momentum eigenstates |P to the one-hole excitations |P, N consisting of N particles. By considering the phase factors exp(−2πipX/L) in (13) the graphs are shifted by X in the x-direction. We have D Nc (x − X) for the truncated series of the delta function modified with the phase factors.

Density notch in the sum over one-hole excitations for free fermions
Through the anti-commutation relation of the field operators we now show that a notch appears in the density profile of a linear combination of one-hole excited states in the 1D free fermions.
Under the periodic conditions of length L the momenta are given by k α = 2πn α /L with some integers n α in the 1D free fermions. We denote by a α and a † α for some integers α the creation and annihilation operators of the one-dimensional free fermions with momenta k α = 2πn α /L, respectively. We assume that an infinite series of integers n α for α = 1, 2, . . . , cover all integers. The field operator and its conjugate are given by Let us take a set of arbitrary integers n 1 , n 2 , . . . , n M and consider the state of M fermions with momenta k α 's For simplicity, we may assume that |M is the ground state. By applying the field operator ψ(0) to it, we have Here we have assumed the ordering of fermionic operators as We define the state |Φ by the sum over all one-hole excitations derived from |M Let us define the local density ρ Φ (x) for the state |Φ of the 1D free fermions by We now show that the local density vanishes at the origin: ρ Φ (0) = 0 for the state. We first recall that the square of the field operator at the origin vanishes, ψ(0) 2 = 0, due to the anti-commutation relations of the field operators. Thus, we have We can show that the integral of the local density over the entire interval L is given by a large positive value such as M − 1 It follows that the local density cannot be always equal to zero, although it vanishes at the origin. Therefore, for the state |Φ which is given by the sum over all the one-hole excited states of |M , the local density ρ Φ (x) has a notch at least at the origin. In section 4 we shall construct the free-fermionic analog of quantum dark soliton states for the 1D free fermions and show that the density profile coincides with those of the LL model if the coupling constant is large enough such as c = 100.

Traveling-wave solutions expressed in terms of elliptic functions
Let us assume the traveling wave solution with velocity v: ψ C (x, t) = ψ C (x − vt) to the nonlinear Schrödinger equation (6). Hereafter in section 3 we denote ψ C simply by ψ. We denote the time derivative by ∂ t ψ = −vψ ′ and the second spatial derivative by ∂ 2 x ψ = ψ ′′ , where ψ ′ denotes dψ/dθ with θ = x − vt. We therefore have We denote the solution of eq. (26) by ψ(x) as a function of x, for simplicity. We express the complex scalar field ψ(x) in terms of the amplitude ρ(x) and the phase ϕ(x), both of which are real, as By substituting (27) into equation (26) we derive a pair of coupled equations from the real and imaginary parts, respectively.
It follows from equation (29) that we have (ρϕ ′ − vρ/2) ′ = 0. By integrating it once with respect to x, we have Here W denotes a constant of integration. Substituting (30) into (28) and multiplying it by Integrating this equation with respect to x, we have where V is another constant of integration. Thus, we have where the potential function U (ρ) is given by Let us assume that the cubic equation: U (ρ) = 0 has three distinct real roots, a 1 , a 2 , and a 3 . We put them in increasing order: a 1 < a 2 < a 3 . The coefficients of the cubic equation are expressed in terms of the roots as follows.
We now derive a solution ρ(x) to equation (33) such that it satisfies a 1 ≤ ρ(x) ≤ a 2 in the interval of x with −L/2 ≤ x ≤ L/2. We set the initial condition: ρ(x = 0) = a 1 , for simplicity. It thus follows from equation (33) that we have Let us define the modulus of Jacobi's elliptic functions, k, by It is clear that 0 < k < 1. Through the transformation of variables from r to z we express the integral (36) in terms of the elliptic integral of the first kind: We therefore have the solution to eq. (33) in terms of Jacobi's elliptic function sn(u, k) We now integrate equation (30). Here we remark that the constant W is given by W = ± √ ca 1 a 2 a 3 . In terms of the elliptic integral of the third kind (A.15) we have where ϕ(0) corresponds to the constant of integration.

Classical dark solitons expressed in terms of elliptic integrals
We now express a 1 , a 2 and a 3 in terms of the complete elliptic integrals. Here we consider the periodic boundary conditions (PBC) for the square amplitude ρ(x) and the phase ϕ(x) of the classical complex scalar field ψ(x, t), and the normalization condition of the square amplitude ρ(x). Then, the elliptic modulus k is related to the depth of the classical dark soliton, i.e. the minimum value a 1 of the square amplitude ρ(x). We express the soliton velocity v in terms of the elliptic integral of the second kind, and determine the chemical potential µ.

PBCs for the square amplitude and the phase of classical dark soliton
Suppose that the square amplitude ρ(x) of a classical dark soliton increases from a 1 to a 2 and then returns back to a 1 with period L. We therefore assign the conditions at the initial and the middle points as It follows from (39) that in terms of the complete elliptic integral of the first kind, K(k), we have We have a 3 − a 1 = 4K 2 /(cL 2 ). Here and hereafter we often denote K(k) by K, suppressing the modulus k.
We assume the periodic boundary condition for the phase: ϕ(x = L) = ϕ(0). Suppose that the soliton velocity is positive: v > 0. It follows that we have W < 0 and the velocity v is expressed in terms of the complete elliptic integral of the third kind as

Normalization condition for the square amplitude
We impose the normalization condition to the square amplitude ρ(x) where n is the density of the number of particles N , i.e. we have n = N/L. By the formula (A.8) we transform equation (40) in terms of Jacobi's dn function into the following.
Thus, by the formula (A.11) in terms of the complete elliptic integral of the second kind E(k) we have 3.2.3 Parametrization of a 1 , a 2 and a 3 in terms of complete elliptic integrals By solving equations (37), (43) and (47) with respect to the roots a 1 , a 2 and a 3 , we express them in terms of the complete elliptic integrals as follows: Since the minimum value of the amplitude must be non-negative, i.e. a 1 ≥ 0, we assign the following condition Thus, if the coupling constant c, the system size L, the density n and the elliptic modulus k satisfy condition (51), the classical dark soliton exists as a solution of the nonlinear Schrödinger equation (6)

Numerical determination of the parameters of classical dark solitons
Let us explain how we evaluate the parameters of a dark soliton such as n and k in order to compare the profiles of the dark soliton with those of a given quantum dark soliton. Here we remark that when a quantum state is given, parameters c, L and the particle number N have been specified. We determine the elliptic modulus k from the depth of the dark soliton, i.e. a 1 : We solve equation (48) as an equation for modulus k. We assign the value N/L on the parameter of density n for the density profile of a quantum state. However, we determine the density n numerically for the squared amplitude profile of the matrix element of the bosonic field operator between different quantum states, since it may take a value smaller than N/L. For the density profile of a quantum state |Φ , the integral of the expectation value Φ|ρ(x)|Φ of the local density operator over the whole space is equal to N/L. However, for the matrix element Ψ|ψ(x)|Φ of the quantum field operatorψ(x) between two states |Φ and |Ψ , the integral of the squared amplitude | Ψ|ψ(x)|Φ | 2 over the whole space may be smaller than the value of N/L. In this case we determine the parameter n numerically by taking the integral of the squared amplitude | Ψ|ψ(x)|Φ | 2 over the whole space. Here we remark that the symbols will be defined in section 4.1.1.
In summary, after we fix the value of density n, we solve equation (48) as a function of modulus k and estimate its numerical value very precisely. In section 4 we shall make the profiles of classical dark soliton for given density profiles of quantum dark soliton by the method explained in the above.

Critical velocity
We define the function f (k) by f (k) := K(k) 3E(k) − (2 − k 2 )K(k) . It appears in the third term of the expression (53) for the chemical potential µ. We can show that it is a monotonically decreasing function of modulus k, which yields the inequality f (k) ≤ f (0) = π 2 /4. This gives an upper bound for the absolute value of soliton velocity v: Here we remark that the critical velocity is defined for the system of a finite size. We shall denote by v c,∞ the critical velocity for an infinite system.

Square amplitude ρ(x) and the phase ϕ(x) expressed in terms of elliptic functions
We now express a classical dark soliton explicitly in terms of the complete elliptic integrals. The square amplitude (46) is given by We shall show that the second term in the right-hand side of equation (55) leads to the logarithmic correction associated with the conservation of the number of particles. Let us introduce a parameter β k by We denote the argument 2Kx/L in (40) by u. Then, the square amplitude is expressed as follows: We now show that the phase ϕ(x) is expressed in terms of Jacobi's Theta function. Let us define a pure imaginary number a = iα with α > 0 by By formula (A.21) we evaluate the elliptic integral of the third kind with Jacobi's Zeta and Theta functions [48], and we express the phase ϕ(x) as the logarithm of a ratio of Theta functions 3.3.2 Parameters in the asymptotic expansion with respect to the system size L Let us consider the limit of sending the system size L to infinity (L → ∞) and the modulus k to 1 (k → 1) simultaneously so that the ratio K(k)/L is kept constant. We denote the ratio by b: b = K/L. Here we remark that we keep the density n constant. We define β by the simultaneous limit of sending L to infinity and k to 1 Let us introduce the complementary modulus k , respectively. We remark that K and iK ′ give quarter periods of Jacobi's elliptic functions. We introduce a parameter p by p = exp(−πK/K ′ ), which is small when modulus k is close to 1. Making use of formulae (A.5), (A.13), (A.14) and Legendre's relation (A.12) we expand K's and E's with respect to p as follows: We have p ≃ exp(−2bL), since we have fixed the ratio K/L = b.
3.3.3 Asymptotic expansion of the velocity and the chemical potential when k approaches 1 We recall that the soliton velocity v is expressed in terms of the complete elliptic integral of the third kind (52). Through (A.22) we express it in terms of Jacobi's Zeta function as follows: Here we recall that parameter α has been defined by (58).
Let us now evaluate Z(iα, k). We first express the left-hand side of (58) in terms of β k as When the modulus k is close to 1, by making use of Jacobi's imaginary transformation of sn function: sn(iα, k) = i sn(α, k ′ )/cn(α, k ′ ) and then by expressing it in terms of Jacobi's Eta and Theta functions, we show tan πα Through (A.30) we express the value of Zeta function Z(a, k) in terms of a = iα as follows.
Applying the expansions (64) and (65) to (62), we evaluate the velocity v upto the order of logarithmic terms and through (53) we evaluate the chemical potential up to the order of logarithmic terms It follows that chemical potential µ approaches 2nc when we send system size L to infinity. We denote the limiting value by µ ∞ . Here we remark that the logarithmic corrections correspond to the finite size corrections since we have log(1/p) = 2bL. Moreover, the chemical potential µ is larger than µ ∞ = 2nc when the system size is finite since we have x > arctan x for x > 0. We can confirm such behavior by checking the numerical estimates for the parameters of classical dark solitons listed in section 4.

Asymptotic expansion of the square amplitude of the classical complex scalar field
Through Jacobi's imaginary transformation for dn function we have We therefore have It is easy to show that the asymptotic expansion of the square amplitude up to the order of 1/ log p, i.e. the order of 1/L, satisfies the normalization condition (45). At k = 1, the square amplitude ρ(x) is expressed in terms of β = 4b 2 /(nc) by Thus, the critical velocity v c with the infinite system size (L = ∞) is given by v c,∞ = 2 √ 3nc − µ ∞ = 2 √ nc. Here we recall µ ∞ = 2nc. It follows that we have Here β is expressed as β = 1 − (v/v c,∞ ) 2 . We have thus derived the square amplitude of one-soliton solution due to Tsuzuki [21].

Asymptotic expansion of the phase of the classical complex scalar field
Applying (A.29) to (59) we derive the following expansion in terms of p = exp(−2bL) for any system size L: Making use of (64) we derive the expansion of ϕ(x) up to the order of the inverse of log(1/p) as follows: Therefore, in the limit of sending k to 1, we obtain the phase of the dark soliton solution by Tsuzuki [21]: Thus, we have shown that the classical dark soliton under the PBCs approaches the dark soliton solution by Tsuzuki [21] through the simultaneous limit of sending the system size to infinity and the modulus k to 1. We define the local density operatorρ(x) byρ(x) :=ψ † (x)ψ(x), in the second-quantized system of the 1D Bose gas interacting through the delta-function potentials (4). Here we denote byψ(x) the field operator at the initial time t = 0:ψ(x, t = 0). We now consider the graph of the expectation value of the local density operator for a quantum state |Φ (i.e. Φ|ρ(x)|Φ ) versus the position coordinate x. Here we recall that we call it the density profile of the state |Φ . We evaluate the expectation value of the local density operator for the state |X, N by the form factor expansion. We express the expectation value as the sum over the form factors between the Bethe eigenstates Here |P, N and |P ′ , N are the normalized Bethe eigenstates of N particles in the type II branch (i.e. they are one-hole excitations) and have total momentum P = 2πp/L and P ′ = 2πp ′ /L, respectively. Here the form factors P ′ , N |ρ(0)|P, N are effectively calculated [24] by the determinant formula for the norms of Bethe eigenstates [49] and that of the form factors of the density operator [50,51] where the quasimomenta {k 1 , · · · , k N } and {k ′ 1 , · · · , k ′ N } lead to the eigenstates |P and |P ′ , respectively. We use the abbreviations k j,ℓ := k j − k ℓ and k ′ j,ℓ := k ′ j − k ′ ℓ . The kernelK(k) is defined byK(k) = 2c/(k 2 + c 2 ). The matrix G(k) is called the Gaudin matrix, whose (j, ℓ) th element is given by (k j,m ) −K(k j,ℓ ) for j, ℓ = 1, 2, · · · , N.  Table 1 is plotted with a blue broken line to each value of c with c = 0.01, 1, 10 and 100.  Table 1: Parameters for the density profiles of the classical dark solitons with L = 20 plotted by blue broken lines in Fig. 2 and in the first and the second columns of Fig. 4.
The matrix elements of the (N − 1) by (N − 1) matrix U (k, k ′ ) are given by

Density profiles of quantum and classical dark solitons at the initial time
We now show that the density profile of quantum dark soliton is in good agreement with that of classical dark soliton, i.e. the profile of the square amplitude of the complex scalar field for the corresponding classical dark soliton: |ψ C (x)| 2 versus coordinate x. In particular, the quantum and classical density profiles are consistent with each other in the weak coupling case such as c = 0.01. In   Table 2 is plotted with a blue broken line to each of c with c = 0.01, 1, 10 and 100.  Table 2: Parameters for the classical dark solitons with L = 500 plotted by blue broken lines in Fig. 3 and in the third and the fourth columns of Fig. 4.
We depict the density profile of classical dark soliton by shifting the profile of square amplitude ρ(x) derived in section 3 by L/2 in the x direction, as shown in Figs. 2 and 3. Here we remark that in section 3 we have put the density notch of any classical dark soliton at the origin by assuming relations (42). On the other hand, we observe that if we set the parameter X as X = 0 in (13), it follows from formula (75) that the notch is located at x = L/2 in the density profile of quantum dark soliton. We therefore shift the profiles of classical dark soliton by L/2 as |ψ C (x)| 2 = ρ(x − L/2).
The numerical estimates for the parameters of the classical dark solitons, i.e. elliptic modulus k, chemical potential µ, soliton velocity v and critical velocity v c for the four values of the coupling constant c (i.e. c = 0.01, 1, 10 and 100) are listed in Table 1 and Table 2 for the cases of N = L = 20 and N = L = 500, respectively. They are obtained by applying the method of section 3.2.4 to the density profiles of quantum dark soliton.

Time evolution of the density profile of quantum dark soliton
Let us now study the dynamics of a quantum dark soliton state and compare it with that of a classical dark soliton. We shall show that the time evolution of the density profile of quantum dark soliton is consistent with that of the corresponding classical dark soliton in the weak coupling case.
Making use of the determinant formula of the form factors (i.e. the matrix elements) (76) we evaluate numerically the expectation values of the local density operator for the quantum dark soliton state at time t and position x by Here we consider the quantum field operatorψ(x, t) at position and time coordinates x and t, respectively, and denote the local density operator byρ(x, t) =ψ † (x, t)ψ(x, t).
We now show the time evolution of the density profile of quantum dark soliton and that of classical dark soliton, explicitly in Fig. 4. Here we remark that the density profile at time t is given by the graph of the expectation value of the local density operatorρ(x, t) in the quantum dark soliton state (i.e. X, N |ρ(x, t)|X, N ) at given time t plotted against the position coordinate x.
When N = L = 20 the snapshots with c = 0.01 and c = 1 at times t are plotted for the density profile of quantum dark soliton by red solid lines in the first and second columns of panels of Fig. 4 from the left, respectively. Here we remark that each column consists of five panels in Fig. 4. In the first column the snapshots are taken at t = 0, 10, 20, 30, and 40, while in the second column at t = 0, 1, 2, 3 and 4. The density profiles of classical dark soliton are plotted by blue broken lines for the cases of c = 0.01 and c = 1 in the first and second columns of Fig. 4 from the left, respectively. Here we recall that the density profile of classical dark soliton at time t is given by the square amplitude of the classical scalar field |ψ C (x − vt)| 2 = ρ(x − vt − L/2). The numerical estimates of the parameters are given in Table 1.
It is clear that in the weak coupling case of c = 0.01 and N = L = 20 the density profile of quantum dark soliton at time t and that of classical dark soliton |ψ C (x − vt)| 2 move together dynamically for a rather long period of time at least such as up to t = 40 or 50. It still keeps the initial form of the density profile in the snapshot at t = 40, as shown in Fig. 4. On the other hand, the density profile with c = 1 collapses in a much shorter period of time such as t = 5.
When N = L = 500 the snapshots with c = 0.01 and c = 1 at times t are plotted for the density profile of quantum dark soliton by red solid lines in the third and fourth columns of panels of Fig. 4 from the left, respectively. In the third column the snapshots are taken at t = 0, 100, 200, 300 and 400, while in the fourth column at t = 0, 1, 2, 3, and 4. The density profiles of classical dark soliton at time t, |ψ C (x − vt)| 2 = ρ(x − vt − L/2), are plotted by blue broken lines for c = 0.01 and c = 1 in the third and fourth columns of panels of Fig. 4 from the left, respectively. Here, numerical estimates of the parameters are given in Table 2.
We observe in the third column of Fig. 4 from the left that the density profile of quantum dark soliton for c = 0.01 moves extremely slowly in time evolution and keeps overlapping with the density profile of classical dark soliton up to t = 100. Then, the density notch collapses slowly over a long period of time, and the density profile relaxes to a nonzero flat profile.

Density profile of the free-fermionic quantum-dark-soliton state
In the large coupling case such as c = 100, the density profile of quantum dark soliton shows ripples at the shoulders of the notch, as shown in Figs. 2 and 3. We now show that they are described in terms of the Dirichlet kernel. They correspond to the Friedel oscillations in the strong coupling limit.  Table  1 and Table 2, respectively.
Let us recall the 1D free fermions in section 2.4. We define a state |Ψ which is analogous to the quantum dark soliton state (13) by The number of fermions in |Ψ is given by N . Here we remark that the state |Ψ is given by the sum over N one-hole states, while the state |Φ defined in (20) consists of the sum over N + 1 one-hole states. By applying the fermionic field operator ψ(x) (19) to the state |Ψ we have (81) By taking the Hermitian conjugate of (81) with x replaced by y and then by multiplying the original one (81) by it, we derive the one-particle density matrix for the state |Ψ , rigorously. The result is given by For an illustration, let us consider the odd-N case with N = 2N c + 1 for an integer N c and assume the set of momenta {k 1 , k 2 , . . . , k N } = {2πn/L| n = 0, ±1, ±2, . . . , ±N c }. It is for the ground state of N free fermions. Here, we put k N +1 = 2π(N c + 1)/L. In terms of the Dirichlet kernel D Nc (x) in (15) we have from (82) In terms of N , both for the odd and even N cases, we can show that the density profile of |Ψ is given by It is easy to show that the local density becomes very small but nonzero at x = 0.
Interestingly, in the large coupling case such as c = 100, the density profiles of quantum dark soliton shown in Figs. 2 and 3 coincide with those of the 1D free fermions (84) over the whole region. We thus have a conjecture that when we send the coupling constant c to infinity, i.e. in the TG limit, the density profile of quantum dark soliton in the 1D Bose gas approaches that of the 1D free fermions.
Let us now derive the time evolution of the density profile for the free-fermionic quantum-dark-soliton state |Ψ . The Hamiltonian of the 1D free fermions is given by When we make connection to the 1D Bose gas, we assume the quadratic dispersion relation: ω α = k 2 α . We now define the initial state |Φ(0) by |Φ(0) = |Φ and introduce its time evolution by applying the time evolution operator e −iHt |Φ(t) = e −iHt |Φ(0) .
Similarly as the case of the one-paticle density matrix at the initial time, we can show that the one-particle density matrix for the state |Φ(t) constructed by the sum over one-hole states is given by Here F N +1 (x, t) = F N +1 (x, t; {k 1 , . . . , k N +1 }) denotes the exponential sum for N + 1 momenta. For an illustration, we consider the odd N case with the set of momenta for the local density (83). We have In the large coupling case such as c = 100, the time evolution of the density profile of quantum dark soliton looks quite similar to that of the free-fermionic quantum-dark-soliton state |Ψ calculated by (89) with (90).

Gaussian weighted sum leading to a smooth density profile
We recall that in the large coupling case such as c = 100, the density profiles of quantum dark soliton show ripples at the shoulders of the notches, as shown in Figs. 2 and 3. They correspond to the oscillations in the Dirichlet kernel. By modifying the construction (13) of quantum dark soliton with the Gaussian weights we can remove the ripples in the density profile of quantum dark soliton so that it consists of only smooth curves. In Appendix B we shall show that a smooth Gaussian profile is obtained by introducing the Gaussian weight to the sum.
However, in the large coupling case, the width of the corresponding classical dark soliton can be much smaller than the width of the notch in the density profile of quantum dark soliton. Let us recall relation (48) in section 3. Here, the parameter a 1 gives the smallest value of the soliton density ρ(x). In the large coupling case such as c = 100, it is given by 1/L, as shown in (85). Here, we set n = N/L = 1. When the modulus k is close to 1, the complete elliptic integrals K and E are given by K = O(log(1/p)) and E = O(1) as shown in the asymptotic behavior (61) in section 3.3. It follows from relation (48)  Let us consider the matrix element of the quantum field operatorψ(x) between the quantum dark soliton states with N particles, |X, N , and that of N − 1 particles, |X ′ , N − 1 . We define the symbol ψ Q (x) by where P = 2πp/L and P ′ = 2πp ′ /L denote the total momenta of the normalized Bethe eigenstates in the type II branch |P, N and |P ′ , N , respectively. Here we remark that the two quantum states have different numbers of particles such as N and N − 1 but they share the same system size L. Hereafter we consider only the case of X = X ′ = 0. The matrix element P ′ , N − 1|ψ(0)|P, N are evaluated effectively by the determinant formula for the norms of Bethe eigenstates [49] and that for the form factors of the field operator [50,52,53] as  where the quasi-momenta {k 1 , · · · , k N } and {k ′ 1 , · · · , k ′ N −1 } lead to the eigenstates |P, N and |P ′ , N − 1 , respectively. Here we have employed the abbreviated symbols k j,ℓ := k j − k ℓ and k ′ j,ℓ := k ′ j − k ′ ℓ . Here we recall that the kernelK(k) is given byK(k) = 2c/(k 2 + c 2 ) and the matrix G(k) denotes the Gaudin matrix, whose (j, ℓ)th element is given in (77). The matrix elements of the (N − 1) by (N − 1) matrix U (k, k ′ ) are given by

Square amplitude and phase profiles for the matrix element of the quantum field operator
Let us denote the phase argument of a given complex number z = r exp(iθ) (0 ≤ r, −π < θ ≤ π) by Arg[z] = θ. For instance, we shall denote the phase of the matrix element ψ Q (x) by Arg[ψ Q (x)].
For N = L = 20, the profiles of square amplitude |ψ Q (x)| 2 and phase Arg[ψ Q (x)]/π are plotted by red solid lines in Figs. 5 and 6, respectively. We recall that they are defined for the matrix element ψ Q (x) between the quantum dark soliton states. The density profiles |ψ C (x)| 2 = ρ(x − L/2) and phase profiles Arg[ψ C (x)]/π = ϕ(x − L/2)/π of classical dark soliton are shown by broken blue lines in Figs. 5 and 6, respectively. The numerical estimates of the soliton parameters are given in Table 3. They are obtained by the numerical method of section 3.2.4.    The square amplitude profile |ψ Q (x)| 2 of quantum dark soliton is consistent with the density profile of classical dark soliton in the weak coupling case of c = 0.01, as shown in the upper left panel of Fig. 5. Moreover, the phase profiles Arg[ψ Q (x)] of quantum dark soliton perfectly agree with those of classical dark soliton for c = 0.01 and c = 0.1, as graphs. Here, no deviation is found between the quantum and classical profiles. We thus suggest that both the profiles of square amplitude |ψ Q (x)| 2 and phase Arg[ψ Q (x)] for the matrix element ψ Q (x) between the quantum dark soliton states should be consistent with those of classical dark soliton, respectively, in the weak coupling limit: c → 0.
For N = 500 the profiles of square amplitude |ψ Q (x)| 2 and phase Arg[ψ Q (x)]/π are plotted in Figs 7 and 8, respectively. The density profile |ψ C (x)| 2 = ρ(x − L/2) and phase profile Arg[ψ C (x)]/π = ϕ(x − L/2)/π of classical dark soliton are shown by broken blue lines. The numerical estimates of the soliton parameters are given in Table 4. They are obtained by the numerical method of section 3.2.4 .
The phase profiles Arg[ψ Q (x)]/π of the matrix element ψ Q (x) perfectly agree with those of the corresponding classical dark solitons for the cases of c = 0.01, 0.1 and 1, respectively. On the other hand, the square amplitude profile |ψ Q (x)| 2 of the matrix element ψ Q (x) is not necessarily in very good agreement with the density profile of classical dark soliton even for the case of c = 0.01. We suggest that it is due to the effect of manybody correlations in the 1D Bose gas such as observed in the study of the BEC fraction [37] as mentioned in Introduction. We also suggest that if the coupling constant c becomes much smaller, the square amplitude profiles |ψ Q (x)| 2 should be in much better agreement with the density profiles of classical dark soliton |ψ C (x)| 2 .

Time evolution of the square amplitude and phase profiles for the matrix element of the quantum field operator
We now show the time evolution of the square amplitude profile and the phase profile of the matrix element of the bosonic quantum field operator with respect to the quantum dark soliton state with N particles and that of N − 1 particles. We denote by ψ Q (x, t) the matrix element of the bosonic quantum field operatorψ(x, t) between the quantum dark soliton states |X, N and |X ′ , N − 1 as follows.
Here we recall that the two quantum states have different particle numbers but have the same system size L. The snapshots at times t for the profiles of the square amplitude and the phase of the matrix element ψ Q (x, t) of the bosonic quantum field operatorψ(x, t) between the quantum dark soliton states |X, N and |X ′ , N − 1 are plotted by red solid lines for the cases of N = L = 20 and N = L = 500 in Figs. 9 and 10, respectively. Here, we consider two values of the coupling constant such as c = 0.01 and c = 1 to each of the particle number N in the cases of N = 20 and 500.
The quantum and classical profiles are in good agreement in many aspects of the time evolution as shown in Figs. 9 and 10. We now describe three aspects as follows. Firstly, in the weak coupling case of c = 0.01, both the time evolution of the square amplitude profile |ψ Q (x, t)| 2 and that of the phase profile Arg[ψ Q (x, t)]/π are consistent with those of the corresponding classical dark soliton, |ψ C (x − vt)| 2 and Arg[ψ C (x − vt)]/π, respectively. The consistency of the quantum profiles of the square amplitude and the phase with the classical profiles of the corresponding classical dark soliton holds over a long period of time such as at least up to t = 40 as shown in the first and second columns of panels in Fig. 9 from the left. Secondly, the quantum profiles move with the same speed as the corresponding classical profiles, as shown in all the eight columns of panels in Figs Table 3. (ii) Snapshots at times t for the profiles of square amplitude |ψ Q (x, t)| 2 and phase Arg[ψ Q (x, t)]/π of quantum dark soliton with N = L = 20 for c = 1 are shown by red solid lines in the third and fourth columns, respectively. The snapshots of the density profile |ψ C (x − vt)| 2 and phase profile Arg[ψ C (x − vt)]/π of classical dark soliton are shown by blue broken lines in the third and fourth columns from the left, respectively. Here, the soliton parameters are shown in Table 3. soliton are shown by blue broken lines in the first and second columns, respectively. The numerical estimates of the soliton parameters are given in Table 3.
(ii) Case of N = 20 for c = 1 With c = 1 and N = L = 20, for the profiles of square amplitude |ψ Q (x, t)| 2 and phase Arg[ψ Q (x, t)]/π of the matrix element ψ Q (x, t), the snapshots at times t are plotted by red solid lines in the third and fourth columns of panels of Fig. 9 from the left, respectively. Here, five snapshots are shown in each column of panels for t = 0, 2, 4, 6 and 8. The density profile |ψ C (x − vt)| 2 = ρ(x − vt − L/2) and phase profile Arg[ψ C (x − vt)]/π = ϕ(x − vt − L/2)/π of classical dark soliton are shown by blue broken lines in the third and fourth columns, respectively. The numerical estimates of the soliton parameters are given in Table 3.
We observe that by comparing the snapshots of c = 0.01 with those of c = 1, the density notch in the square amplitude profile |ψ Q (x, t)| 2 at t = 0 decays in a shorter period of time as the coupling constant becomes larger, as shown in the first and third columns of Fig. 9 for c = 0.01 and c = 1, respectively. However, the position of the density notch moves together with that of classical dark soliton, although the notch itself becomes wider rather rapidly in time. In the phase profiles, the abrut step becomes softer and milder gradually in time, but the change in the phase profile is much slower than in the square amplitude profile.  Table 4.
We observe that the square amplitude |ψ Q (x, t)| 2 for the matrix element of the quantum field operatorψ(x, t) approaches zero in time evolution. The decaying behavior is clear for c = 1 as shown in the third column from the left of Figs. 9 and 10. In the first column of Fig. 10 from the left, we also observe it for c = 0.01 and N = L = 500 over a very long period of time up to t = 4000. However, for the density profiles of quantum dark soliton, the density approaches a nonzero constant value in time evolution, as shown in Fig. 4. We suggest that the quantum dark soliton states with different particle numbers such as N and N − 1 do not have the same set of energy eigenvalues in common, so that the off-diagonal matrix element ψ Q (x, t) vanishes in time, finally. (iv) Case of N = 500 for c = 1 With c = 1 and N = L = 500, for the profiles of square amplitude |ψ Q (x, t)| 2 and phase Arg[ψ Q (x, t)]/π in the matrix element ψ Q (x, t), the snapshots at times t are plotted by red solid lines in the third and fourth columns of panels of Fig. 10 from the left, respectively. Here, five snapshots are taken for almost the same intervals of time as in the case of Fig. 9 such as for t = 0, 3, 6, 9 and 12 in both cases of the profiles of square amplitude |ψ Q (x, t)| 2 and phase Arg[ψ Q (x, t)]. The density profile |ψ C (x − vt)| 2 = ρ(x − vt − L/2) and phase profile Arg[ψ C (x − vt)]/π = ϕ(x − vt − L/2)/π of classical dark soliton are shown by blue broken lines in the third and the fourth columns of Fig. 10 from the left, respectively. The numerical estimates of the soliton parameters are given in Table 4.
Here we remark that when we evaluate the matrix element ψ Q (x, t) straightforwardly by making use of the form factor formula, the phase of the matrix element ψ Q (x, t) increases with a constant velocity in time, and the phase profile moves downwards with a constant velocity in the animation. We suggest that it is related to the difference between the zero point energy of the eigenstates with N − 1 particles and that of N particles. In Figs. 9 and 10 we plotted the phase profiles which are subtracted by the constant shift in time.

Matrix element of the fermionic field operator between the free-fermionic quantum dark soliton states
As the quantum dark soliton states for the 1D free fermions with N and N −1 particles we consider the following  Table 4. (ii) Snapshots at times t for the profiles of square amplitude |ψ Q (x, t)| 2 and phase Arg[ψ Q (x, t)]/π of quantum dark soliton with N = L = 500 for c = 1 are shown by red solid lines in the third and fourth columns, respectively. The snapshots at times t for the density profile |ψ C (x − vt)| 2 and phase profile Arg[ψ C (x − vt)]/π of classical dark soliton are shown in the third and fourth columns from the left, respectively. Here, the soliton parameters are shown in Table 4.
We calculate the matrix element of the field operator between the two states as Here, the exponential sum is given by the Dirichlet kernel. For instance, when we take k N +1 = 2π/L and momenta 2πn/L for n = 0, -1, and ±2, ±3, . . . , ±N c where N = 2N c + 1 we have The square-amplitude and phase profiles for the matrix element of the fermionic field operator (99) have several features in common with those of the 1D Bose gas in the large coupling case shown in Figs. 5, 6, 7, and 8, although they do not overlap with those of the 1D Bose gas even in the large coupling case. For instance, the square-amplitude profile vanishes at the origin, and it has a density notch around at the origin and small ripples at the shoulders of the notch; the phase profiles change signs at the origin, etc.. 5 Stability of density profiles of quantum dark soliton 5.1 Quantum speed limit time in the weak coupling limit For a given quantum state we denote the expectation value of the Hamiltonian measured from the ground state energy by E, and the square root of the variance of the expectation value of the Hamiltonian by ∆E. We define the quantum speed limit time T 0 [33,34,35] by We remark that the quantum speed limit time gives a typical lifetime of a generic quantum state. Let us evaluate the quantum speed limit time T 0 for the quantum dark soliton state |X, N in the weak coupling limit where we send the coupling constant c to zero. We solve the Bethe ansatz equations (9) in the weak coupling limit. Here we recall that the Bethe quantum numbers for one-hole excitations are defined by eq. (12). Let us consider an expansion of function arctan x which is valid if the absolute value of x is larger than 1.
Here sign(x) = 1 for x > 0 while sign(x) = −1 for x < 0. Applying this expansion we show in the weak coupling case (0 < c ≪ 1) We have the solution of one-hole excited state |P, N approximately as follows.
Therefore, the one-hole excitation |P, N has the energy eigenvalue E p = p (2π/L) 2 . Let us denote by E the expectation value of the energy eigenvalues for the quantum dark soliton state |X, N . We have E = N −1 p=0 E p /N and hence Since the ground-state energy is given by 0 in the weak coupling limit, we denote it also by E. We calculate the variance of the energy eigenvalues |X, N as (E − E ) 2 = (N 2 − 1) (2π/L) 4 /12 . We thus have the square root of the variance, ∆E, as It follows that the quantum speed limit time for the quantum dark soliton state in the weak coupling limit is given by that for the inverse of the square root of variance π /2∆E Here we remark that the quantum speed limit time T 0 is constant in the weak coupling limit where only the coupling constant c approaches zero while the density n and the particle number N are kept constant.

Observed collapse time of a notch in the density profile
The estimates of the relaxation time (or collapse time) for the density notch in the density profile of quantum dark soliton for the case of L = N = 20 are plotted against coupling constant c in the double logarithmic scale in Fig. 11. By observing the time evolution of the density profile of quantum dark soliton for a given value of the coupling constant c we estimated the relaxation time by the time when the smallest value of the density notch in the density profile reaches the value of 0.5. Here we recall that at the initial time t = 0 the density profile of quantum dark soliton is consistent with that of classical dark soliton, and then the bottom level of the density profile of quantum dark soliton increases gradually in time evolution. In Fig. 11 we observe that the collapse time T col increases as the coupling constant c approaches zero. It is almost inversely proportional to the coupling constant c in the weak coupling regime: Hence we suggest that it becomes infinitely large in the weak coupling limit.

Stability of a notch in the density profile of quantum dark soliton
The quantum speed limit time T 0 for quantum dark soliton states approaches a non-zero constant value when we send the coupling constant c to zero, if we keep the system size L and the number of particles N fixed. In fact, in the weak coupling limit, all the solutions of the Bethe ansatz equations are given by some integral multiples of 2π/L. Hence, the quantum speed limit time of any quantum state is given by some finite value.
On the other hand, the collapse time of the density notch in the density profiles of quantum dark soliton is approximately given by the inverse of the coupling constant, T col ≃ 1/c, in the very weak coupling case. Therefore, the collapse time of the density notch in the density profile of quantum dark soliton is much longer than the quantum speed limit time when the coupling constant c is very small.
For the quantum dark soliton state with N particles in the system size L we have approximately In the case of N = L = 20 we have a rough estimate: T 0 ≃ 2.75. Here we recall that we have set = 1. In Fig. 11 the collapse time is given by approximately 10.0 for c = 0.1. Thus, for c < 0.1 the collapse time of the density profile in the quantum dark soliton state is much longer than the quantum speed limit time.
We conclude that the density notch in the density profile of quantum dark soliton, i.e. that of state |X, N , has a much longer lifetime than a generic quantum state when the coupling constant c is very small but nonzero in the 1D Bose gas interacting through the repulsive delta-function potentials.

A.3 Integral formulas
(A.12) Let us define τ and q by τ = iK ′ /K and q = exp(πiτ ), respectively [48]. We have (1 − q 2n ) 2 (1 + q 2n−1 ) 4 . Let us take a complex number a satisfying the following relation: √ n = k sn(a, k) . (A.20) We can express the elliptic integral of the third kind in terms of Jacobi's Zeta function as follows.
Π(n, k) = K 1 + Z(a, k) n (k 2 − n)(1 − n) . In order to reduce the short-wavelength oscillations or ripples such as shown in the graph of the Dirichlet kernel D Nc (x) we may modify the truncated sum (15) with some weights. Let us denote by k the wavenumbers k = 2πn/L for integers n. We take an integer n 0 and set k 0 = 2πn 0 /L as a reference wavenumber. We define ∆n by ∆n = n − n 0 and denote the difference of wavenumbers by ∆k = k − k 0 = 2π∆n/L. With the Gaussian weights exp(−σ 2 (k − k 0 ) 2 ) we modify the truncated sum (15)  We can show that the sum over ∆n in (B.2) does not become very large in the region of x satisfying |x|/σ < 4πσN c /L. If |x|/σ > 4πσN c /L, each term in the sum (B.2) has its absolute value larger than 1, and hence the sum can be much larger than 2N c + 1 which is the number of the terms in the sum. Here we remark that the bound 4πσN c /L is given by π if we set σ = L/(4N c ). We may therefore assume that the sum over ∆n in (B.2) is approximately constant in some region near the origin such as |x| < 4σ. Evaluating the sum by an integral we approximate the Gaussian weighted sum D σ Nc (x) as Here the error function erf(z) is given by Typically, by putting σ = L/(4N c + 2), i.e. half the period of the short-distance oscillations or ripples in the graph of the function D Nc (x), the ripples disappear, and the Gaussian weighted sum D σ Nc (x) approximately becomes a smooth Gaussian wave packet with the standard deviation 2σ.