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

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


Introduction
The experimental realization of trapped atomic gases in one dimension (1D) has provided a new motivation for the study of strong correlations in fundamental quantum mechanical systems of interacting particles [1,2,3,4,5]. Furthermore, the nonequilibrium dynamics of closed interacting quantum systems is now extensively studied in 1D by experiments and theories [6,7,8]. In many 1D quantum interacting systems quantum fluctuations may play a key role and often lead to subtle nontrivial effects. We thus expect that fundamental many-body properties such as the quasi-Bose-Einstein condensation (BEC) should play a key role in the nontrivial quantum dynamics such as quantum dark-solitons. We shall define it shortly with the Penrose-Onsager criterion.
Let us introduce a theoretical model for the 1D system of interacting bosons with repulsive short-range potentials. Here we call it the 1D Bose gas. For simplicity we assume that the interactions are given by the delta-function potentials, since they give nontrivial effects in the 1D case although they are simple. For instance, the scattering length depends on the strength of the delta-function potential in 1D systems. We thus have the Lieb-Liniger model (LL model) as the system of the 1D Bose gas. The Hamiltonian of the LL model is given by [9,10] Here N denotes the number of bosons, and we assume the periodic boundary conditions of the system size L on the wave-functions. We employ a system of units with 2m = = 1, where m denotes the mass of the particle. We recall that the coupling constant c is positive. It is an exactly solvable model of the 1D quantum manybody system. It is known that all the eigenvectors are constructed by the Betheansatz method [11]. Furthermore, the Gross-Pitaevskii (GP) equation appears as the Heisenberg equation of motion for the second-quantized Hamiltonian of the LL model. It is expressed in terms of the classical complex scalar field ψ as follows [12].
We expect that the GP equation should play a central role in the long-distance meanfield behavior of the 1D Bose gas in some quantum state if the quasi-BEC occurs in the quantum state of the LL model especially in the weak-coupling regime. If it is the case, the solution of the GP equation should correspond to the macroscopic wave-function of the quasi-BEC state, and describe the quantum state well at least approximately. We define the quasi-BEC by the criterion due to Penrose and Onsager [13,14] (see also Section 4.2). Suppose that particle number N is very large but finite. The density matrix at zero temperature is given by the ground state |λ of the system asρ = |λ λ|. Then, we define the one-particle reduced density matrix by its partial trace with respect to all but one degree of freedom:ρ 1 = N tr 23···Nρ . Let N 0 denote the largest eigenvalue of the one-particle reduced density matrixρ 1 . If it is of order N , i.e., the ratio n 0 = N 0 /N is nonzero and finite for large N , then we say that the system exhibits the quasi-BEC, and we call n 0 the condensate fraction.
If the quasi-BEC occurs in some quantum states of the LL model, we expect that the GP equation should play a central role for characterizing the quantum state, although it is only a partial differential equation for a complex scalar variable. In the present research, we assume that the quasi-BEC should occur if the coupling constant is small enough with respect to the system size or the number of bosons, and hence some solutions of the GP equation such as multiple dark-solitons can be compared with the density profiles of some quantum states in the quasi BEC of the 1D Bose gas. In fact, we shall show a finite-size scaling of the quasi BEC in the present research.
It should be emphasized that such quantum states whose density profiles coincide with those of single dark-solitons of the GP equation have been constructed explicitly in the form of superposition of the yrast states in the Lieb-Liniger model [15]. The construction resolved a long standing problem suggested by Ishikawa and Takayama almost forty years ago [16]. Here we remark that it was shown through the strong coupling limit [17,18] that the yrast states and the mean-field solitons are closely related to each other with respect to quantum numbers. Furthermore, several significant properties in the non-equilibrium dynamics of a quantum single dark-soliton have been exactly investigated [19] and the generic and the ideal Gaussian weights have been introduced [20,21]. Moreover, the density and phase profiles of quantum states of double dark-solitons have been explicitly constructed [22], and the phase shift has numerically been estimated in the scattering of two quantum dark-solitons [23].
There is another aspect of quantum dark-soliton states. Successive measurements of particle positions in the Lieb-Liniger model also leads to observing quantum darksolitons numerically [24,25]. There is a question of how the density profile of a superposition of yrast states is related to the successive measurements of particle positions. When the coupling constant c is equal to zero it was analytically shown that the construction of the quantum dark-soliton state with the Gaussian weight [21] is related to the particle position method [24] as shown in Ref. [21]. When the coupling constant is small and nonzero: c > 0, an ansatz was proposed to bridge between the calculation of single-particle density and the particle position method [26].
In the present paper we show various novel aspects in the exact non-equilibrium dynamics of quantum double dark-solitons, which give pairs of notches in the density profiles, by explicitly constructing corresponding quantum states in the Lieb-Liniger model of the 1D Bose gas with the repulsive interactions. For instance, we exhibit the time evolution of the density profile of the double dark-soliton whose two notches are located at the same position, and that of the phase profiles of the quantum double dark-solitons. In particular, we give an example where the winding number of the phase profile changes during the scattering process of two notches. Furthermore, we also show an exact finite-size scaling of the fraction of the BEC for the ground state. It should characterize the quasi-BEC which we assume to occur in quantum double dark-soliton states in the weak coupling regime. We show that if the coupling constant decreases as a power of the system size, condensate fraction does not vanish and remains constant when we send the system size to a very large value with fixed density. We recall that if the condensate fraction is nonzero for a large particle number N , we call it the quasi-BEC by employing the Penrose-Onsager criterion. It follows from it that the quasi-BEC occurs only if the coupling constant is very small with respect to the system size. Therefore quantum states of dark-solitons may appear particularly in the weak coupling regime.
Based on the definition of the quasi-BEC we derive a kind of macroscopic quantum wave-function by exactly deriving the amplitude and phase profiles of the matrix element of the bosonic field operator, by making use of Slavnov's formula of form factors [27]. Here we recall that the bosonic field operator is defined in the second-quantized Hamiltonian of the Lieb-Liniger model [28].
Let us briefly summarize the finite-size scaling of the quasi-BEC for the ground state, which we shall show in detail in Section 4. The scaling behavior of the quasi-BEC in the 1D Bose gas is fundamental when we send particle number N or system size L to very large values. We define the interaction parameter γ by γ = c/n with coupling constant c in the delta-function potentials and density n = N/L. We show that if γ is given by a negative power of N , i.e. γ = A/N η , condensate fraction n 0 is nonzero and constant for any large value of L or N . We also show that exponent η and amplitude A are independent of density n, and evaluate them as functions of n 0 . Thus, the condensate fraction n 0 for the ground state is given by a scaling function of variable γN η , which corresponds to amplitude A. If the condensate fraction of a given quantum state with large N is nonzero in the 1D Bose gas, we suggest that the classical mean-field approximation such as the GP equation should be valid for the state [15]. Furthermore, we show that the 1D Bose gas of a finite particle number may have the same condensate fraction for any large L in the case of the ground state.
Finally, we mention some potentially relevant results in the following. For strong and intermediate interaction strengths, the Lieb-Liniger Gross-Pitaevski equation is introduced, which is an extension of the GP equation [29]. Associated with the quantum states of dark solitons, bound states of dark solitons are numerically studied by solving the GP equation [30], dynamics of a bright soliton in the quasi-BEC with timedependent atomic scattering length in a repulsive parabolic potential [31], quantized quasi-two-dimensional Bose-Einstein condensates with spatially modulated nonlinearity [32], matter rogue wave in Bose-Einstein condensates with attractive atomic interaction [33], exact soliton solutions, and nonlinear modulation instability in spinor Bose-Einstein condensates [34].
The contents of the paper consist of the following. In Section 2 we explain the Bethe ansatz and useful formulas for evaluating the form factors of the field operator. We also define the winding number for solutions of the GP equation under the periodic boundary conditions. In Section 3 we show the time evolution of the quantum double dark-soliton state constructed with equal weight for the following two cases: (i) The soliton positions X 1 and X 2 are different: X 1 = L/4 and X 2 = 3L/4; (ii) the soliton positions are the same: X 1 = X 2 = 0. We also show the time evolution of the quantum double darksoliton state constructed with the Gaussian weights. Here, two notches have different speeds thanks to the Gaussian weights, and we evaluate the phase shift in the collision of the two dark solitons. We remark that two notches have mostly the same speed if the quantum double dark-soliton state is constructed with equal weight. In Section 4 we show the finite-size scaling behavior of the condensate fraction in the ground state for the 1D Bose gas with repulsive interactions at zero temperature. According to it, we can estimate that the fraction of the quasi-BEC condensate should be equal to 0.99 for the quantum double dark-soliton state with N = L = 20 and c = 0.05 studied in the present research.

Bethe ansatz equations
In the LL model, the Bethe ansatz offers an exact eigenstate with an exact energy eigenvalue for a given set of quasi-momenta k 1 , k 2 , . . . , k N satisfying the Bethe ansatz equations (BAE) for j = 1, 2, . . . , N : Here I j 's are integers for odd N and half-odd integers for even N . We call them the Bethe quantum numbers. The total momentum P and the energy eigenvalue E are expressed in terms of the quasi-momenta as If we specify a set of Bethe quantum numbers I 1 < · · · < I N , the BAE in Equation (3) have a unique real solution k 1 < · · · < k N [28,11]. In particular, the sequence of the Bethe quantum numbers of the ground state is given by The Bethe quantum numbers for low lying excitations are systematically derived by putting holes or particles in the perfectly regular ground-state sequence.

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

Quantum double dark-soliton state
A quantum state that has two notches in both profiles of density and square amplitude of the matrix element of the field operator was proposed in [22]. We call it the quantum double dark-soliton state, and it is given by the superposition of "two-hole" excitation states as follows. with a normalization factor M N for N particles. The quantum state |p 1 , p 2 , N is characterized by a configuration of Bethe quantum numbers that has two vacancies located at p 1 and p 2 in the series of the Bethe quantum numbers, which is illustrated in Figure 1 (a). This configuration represents the Bethe quantum numbers of the ground state of N particles along with those of additional two particles. In Equation (6) Figure 1 some configurations with two holes p 1 and p 2 are exhibited. In the third configuration, two holes p 1 and p 2 are located in its middle part of the series which corresponds to the ground state of N particles. Here we remark that in order for two notches have positive velocities we derive two hole excitations derived from the configuration constructed by adding two particles to the right of the "Fermi momentum" as shown in Figure 1 (a). If we add the two particles to the right and left of the "Fermi momentum" symmetrically, then the sum of the momenta vanishes. The density profile of this state X 1 , X 2 , N |ψ † (x)ψ(x)|X 1 , X 2 , N shows the two density notches at the positions x = X 1 , X 2 , which coincides with the squared amplitude of the elliptic soliton [22]. Here, by the determinant formula for the norms of Bethe eigenstates [35,36] we can effectively evaluate the matrix element × e −i(p 1 X 1 +p 2 X 2 ) p 1 , p 2 , N |ψ † (0)ψ(0)|p 1 , p 2 , N .
Here, P and P in an exponential term denote the total momentum of the state |p 1 , p 2 , N and |p 1 , p 2 , N calculated through Equation (4), respectively. The sum in the above equation is taken over all pairs of p = {p 1 , p 2 } and p = {p 1 , p 2 } that belong to the set P N . The matrix element of the form factors of the density operator [27,37,38] is given by where the quasimomenta {k 1 , · · · , k N } and {k 1 , · · · , k N } give the eigenstates |p 1 , p 2 , N and |p 1 , p 2 , N , 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 The matrix elements of the (N − 1) by (N − 1) matrix U (k, k ) are given by We have also considered the matrix element of the single field operator where P and P denote the total momenta of the state |p 1 , p 2 , N and |p 1 , p 2 , N − 1 , respectively. The determinant formula is given by [35,36,27,39,37,38] where the quasi-momenta {k 1 , · · · , k N } and {k 1 , · · · , k N −1 } give the eigenstates |p 1 , p 2 , N and |p 1 , p 2 , N − 1 , respectively. We recall that the matrix G(k) denotes the Gaudin matrix, whose (j, )th element is given in Equation (9). The matrix elements of the (N − 1) by (N − 1) matrix U (k, k ) are given by 2.4. One-particle reduced density matrix The matrix element of the one-particle reduced density matrix, ρ 1 (x, y) := x|ρ 1 |y , for a quantum system is expressed as a correlation function in the ground state |λ : In the LL model we can numerically evaluate the correlation function by the form factor expansion. Inserting the complete system of eigenstates, µ |µ µ|, we have where P µ denotes the momentum eigenvalues of eigenstates |µ . Each form factor in the sum (15) is expressed as a product of determinants by making use of the determinant formula for the norms of Bethe eigenstates [35] and that for the form factors of the field operator [27,38,37]: where the quasi-momenta {k 1 , · · · , k N } and {k 1 , · · · , k N −1 } give the eigenstates |λ and |µ , respectively. Here we have employed the abbreviated symbols k j, := k j − k and k j, := k j − k . The matrix G(k) is the Gaudin matrix, whose (j, )th element is given . The matrix elements of the (N − 1) by (N − 1) matrix U (k, k ) are given by [37,27,38,35] For the ground state |λ we have shown that the sum of the form factor expansion is almost saturated for the one-particle and one-hole (1p1h) excitations together with twoparticles and two-holes (2p2h) excitations. The saturation rate is explicitly presented in Table 1 of Section 4.3. However, for excited states the saturation rate has not been evaluated. It should be technically nontrivial to evaluate it for excited states. For the quantum states of double dark-solitons, we suggest that the saturation rate should be close to one in the weak coupling case in the form factor expansion up to some excitations with relatively small numbers of particles and holes. It is based on the observation that the density profiles of quantum double dark-soliton states are similar to those of the double dark-solitons of the GP equation, as we shall show in Section 3.

Winding number
We introduce the winding number J associated with solutions of the GP equation under the periodic boundary conditions. Let us assume that a solution of the GP equation φ(x) = ρ(z) exp[iϕ(x)] satisfies the periodic boundary conditions: where J is an arbitrary integer. The integer J is called the winding number [17,18].
In the previous study, we constructed the quantum single dark-soliton with a nonzerowinding number.

Time evolution of quantum double dark-soliton state constructed with equal weight
By making use of the time dependent field operatorψ(x, t), the local density and the matrix element of the quantum state at a given time t are expressed as follows.
where E is the energy of the state |p 1 , p 2 , N , andρ(x, t) =ψ † (x, t)ψ(x, t) denotes the local density operator. We have obtained the exact expressions of the time evolution in Equations (19) and (20) since the Bethe ansatz method gives the exact energies for the quantum state |X 1 , X 2 , N .
3.1.1. Quantum dark-soliton located at X 1 = L/4 and X 2 = 3L/4 initially Figure 2 shows the time evolution of the density profile , i.e., the graph of ρ Q (x, t) versus x at a given time t, for the quantum double dark-soliton state with initial soliton positions X 1 = L 4 and X 2 = 3L 4 under the periodic boundary conditions. We call the plot in the left panel of Figure 2 the two-dimensional (2D) density plot of the local density. Here, the value of the local density ρ Q (x, t) at position x and time t is expressed by the brightness of the point at (x, t) in the space-time diagram, where the horizontal axis corresponds to the x coordinate, while the vertical axis to time t. In the right panels of Figure 2 snapshots of the density profile of ρ Q (x, t) at t = 0, 2, 4, and 11 are plotted.
We note that the density profile shown in panel (a) of Figure 2 is identical to the upper-left panel of Figure 9 for c = 0.05 in Ref. [22]. In the latter panel it was shown that the density profile of the quantum double dark-soliton state completely coincides with the density profile of the elliptic double dark-soliton solution of the GP equation. Thus, at t = 0, the density profile of the quantum double dark-soliton state coincides with that of the elliptic soliton solution of the GP equation.
The positions of notches are expressed by the areas of the darker color in the 2D density plot at the left panel of Figure 2. The trajectories of the positions of the two notches in the density profile are given by two parallel linearly elongated regions in the diagram of time t and coordinate x, as shown in the left panel of Figure 2. Thus, the two notches moves at the same velocity in the positive x direction. In the snapshots of the density profiles, the soliton notches are gradually filled, i.e., they become shallower in time evolution, as shown in panels (a), (b), (c), and (d) of Figure 2. That is, the distance between the bottoms of the notches is kept constant through the time evolution, while the depths of the notches become smaller. Here we have defined the depth of a notch by the difference between the largest and smallest values in the density profile. For example, at t = 11, the notches are located at x 1 = 1.9115 and x 2 = 11.9115, and the distance between the two notches is given by ∆x = x 1 − x 2 = 10 = L/2, which is equal to that of t = 0.
It was reported in Ref. [40] that quantum double dark-solitons with notches of almost the same depths can appear again after their depths of notches become much smaller over a time scale of 1/c. However, the quantum double dark-soliton states constructed in the present research do not show this reappearing or recurrent behavior in time evolution. Once the soliton notches in the density profile are completely filled, i.e., their depths vanish, the density profile remains flat and uniform in time evolution, as illustrated in Figure 2. We note that the construction of the quantum soliton in Ref. [40] is different from that of the present research, and also that the number of particles is equal to N = 8 in Ref. [40], which is smaller than N = 20 for the system in Figure 2. The notches in the density profile of ρ Q (x, t) and those in the profile of the square amplitude |ψ Q (x, t)| 2 of matrix element ψ Q (x, t) exhibit different decaying behaviors in time evolution. Figure 3 shows the time evolution of the square amplitude profile of matrix element ψ Q (x, t) with initial soliton positions X 1 = L 4 and X 2 = 3L 4 under the periodic boundary conditions. The average density is decreasing in the time evolution of the profile of the square amplitude |ψ Q (x, t)| 2 in Figure 3, while the notches in Figure 2 are filled gradually. In the density profile, the average density is kept constant as time t increases, since the density is conserved as a whole for any time t: On the other hand, we suggest that the amplitude of the matrix element between the two different quantum states of double dark-soliton should gradually decrease and finally vanish in time evolution, since they have different energies and particle numbers.
In the 2D density plot at the left panel of Figure 3, the trajectories of notches in the space-time diagram are depicted by linearly elongated parallel regions with darker color . The values at the bottoms of the notches are almost equal to zero constantly in time evolution in panels (a), (b), (c), and (d) of Figure 3. Consequently, Figure 3 shows the trajectories of the notches more clearly than Figure 2, as depicted in the 2D density plot at the left panel.
The snapshots of the phase profile at different times in time evolution are shown in Figure 4. Here we remark that the phase is given by the argument of the matrix element of Equation (20) as a complex number. In Figure 4 the abrupt jumps of the phase profile are located at the positions of the notches in Figure 3. The abrupt jumps of the phase profile move with the same constant velocity as the notches in the square amplitude profile. Furthermore, the whole phase profile is gradually shifted toward the negative direction in time evolution. Moreover, the shape of the phase profile as a whole remains the same at least up to t = 40. At the initial time t = 0, the profiles of the square amplitude and the phase of the matrix element ψ Q (x, t) shown in panels (a) of Figure 3 and Figure 4 are identical to those of Figures 10 and 11 for c = 0.05 in Ref. [22], respectively. Panel (a) of Figure 3, the square amplitude profile of the matrix element, corresponds to the panel of c = 0.05 in Figure 10 of Ref. [22], where it was shown that the square amplitude profile of the classical and quantum double dark-soliton overlap completely. Panel (a) of Figure 4, the phase profile of the matrix element, corresponds to the panel of c = 0.05 in Figure 11 of Ref. [22], where the phase profiles of the classical and quantum double dark-solitons overlap completely.
However, the time evolution of the phase profile in the quantum double darksoliton state is different from that of the elliptic dark-soliton solution, which is given by the travelling wave solution of the GP equation.
We recall that the phase profile in the quantum double dark-soliton is gradually shifted toward the negative direction in time evolution in Figure 4, while the phase profile of the travelling wave solution is not shifted. Thus, the time evolution of the quantum dark-solitons that we have constructed is slightly different from that of the classical elliptic soliton solution.
We remark that two notches have mostly the same velocity as shown in Figures 2  and 3 for the quantum double dark-soliton constructed with equal weight. In Section 3.2 we shall show that two notches have different velocities for the quantum double dark-soliton state constructed with the Gaussian weights.
3.1.2. Quantum dark-soliton positions located at X 1 = X 2 = 0 initially By placing the positions of the notches for the quantum dark-solitons X 1 and X 2 at the same point , the profiles of the density and square amplitude derived in time evolution are plotted in Figures 5 and 6, respectively. In both profiles of the density and the square amplitude it seems as if the two notches repel each other in time evolution. The quantum double dark-soliton state with overlapping positions of two notches has different properties in the profiles of the density and the square amplitude from the quantum double dark-soliton state in Equation (6) with different initial positions of two notches as X 1 = L 4 and X 2 = 3L 4 . In the density profile of Figure 5 the notches are much deeper than those in Figure 2, similarly as the notches in Figure 6 of the square amplitude profile. In the profile of the square amplitude , the values at the bottoms of the notches increase in time evolution: The values of the square amplitude at the bottoms of notches are not close to zero at t = 11 in Figure 6. Thus, for the quantum double dark-soliton with two overlapping positions of notches the difference between the density profile and the square amplitude profile is smaller than in Figure 2 and Figure 3.
The snapshots of the phase profile in time evolution are exhibited in Figure 7 for  . At t = 0, the winding number was given by J = 2, while it suddenly changed to J = 1 at t = 0.05. After the change of the winding number, the phase profile became smoother in shape gradually in time evolution. Furthermore, we observe in Figure 7 that the whole phase profile was shifted toward the negative direction step-by-step in time evolution. It is also the case in Figure 4: The whole phase profile was shifted in the negative direction for the quantum double dark-soliton state with initial positions of notches placed at X 1 = L 4 and X 2 = 3L 4 . The abrupt change of the winding number may occur in time evolution for the phase profile associated with the quantum states, i.e., the phase profile of the matrix element of the field operator between the quantum double dark-soliton states in Equation (20). The boundary condition of the phase is given by the form of Equation (18) for solutions of the GP equations and also for the phase profile associated with the quantum states in Equation (20). However, the quantum states do not depend on the boundary conditions of classical solutions. It is sufficient if the phase profile associated with the quantum states satisfies one of the boundary conditions of Equation (18) specified by an integer J, which we have called the winding number. Thus, the winding number J may change abruptly in time evolution in the phase profile associated with the quantum double dark-soliton states in Equation (20).

Time evolution of quantum double dark-soliton state with the ideal Gaussian weights
Let us consider the Gaussian weighted superposition of the excited states consisting of two particle-hole excitations which are determined by a pair of holes p = {p 1 , p 2 } in the set P : Here, N is a normalization factor and the set P is the same as given in Section 2.3. The Gaussian function is given by with two Gaussian parameters (P, σ) [21]. The parameters P and σ are determined by the target soliton depth d and the density n = N/L: Here we have defined the soliton depth d by the smallest value in the density profile of a single dark-soliton. It is different from the "depth of a notch" defined in Section 3.1.1. The target soliton depth d is expressed with the dark soliton solution to the GP equation moving with velocity v in the thermodynamic limit φ ∞ P (x) [21]: Here |φ ∞ P (x = 0)| denotes the square root of the local density at the origin, which is the position of the notch in the thermodynamic limit, and v c,∞ is called the critical velocity of the infinite system. When the system size L is finite, the largest velocity of the elliptic dark-soliton solution of the GP equation is denoted by the critical velocity v c [22]. It approaches the critical value v c,∞ in the limit of sending the system size L to ∞.
The exact profiles in time evolution are numerically derived for the local density ρ Q (x, t) and the square amplitude of the matrix element ψ Q (x, t) of the field operator by calculating the time-dependent matrix elements of the field operator between the Gaussian weighted quantum states of Equation (21), similarly as we have demonstrated in Equation (19) and Equation (20) of Section 3.1 for the quantum double dark-soliton state constructed with equal weight. For the Gaussian weighted quantum double darksoliton state, by assigning a pair of proper values of the target soliton depth d to the two notches of a given superposition of quantum states of Equation (21), we can construct a quantum double dark-soliton state such that its density profile has two distinct notches with different depths.
We have constructed several quantum double dark-soliton states in which the density profile has two distinct notches with different depths. In Figures 8, 9, 10, 11, 12 and 13 we set the target soliton depths as d = 0.6 and d = 0.0 to the two notches, respectively, and we generated the Gaussian weights by making use of Equation (22). Here, the corresponding Gaussian parameters are given by (P 0 , σ) = (0.124027π, 0.106667) and (P 0 , σ ) = (π, 0.421637), respectively, which are derived by making use of Equations (23) and (24). We have thus obtained the quantum double dark-soliton state of distinct narrow notches with different depths. Here we recall that single dark-solitons with different depths have different speeds in the same direction for the GP equation.
We observe the scattering of two notches in the density and phase profiles of the quantum double dark-soliton state. It exhibits the phase shift which is a characteristic property in soliton-soliton collisions [41,23], as shown in the density profile. We remark that the 2D density plot of the local density in the space-time diagram and the snapshots of the density profile at different times are presented in Figure 8 for the quantum double dark-soliton state constructed with the Gaussian weights for c = 0.05. As the two notches of the double dark-soliton approached each other, they moved along approximately straight and linear trajectories with different constant velocities. The collision occurred around at a time interval including t = 11 (see panel (c), which corresponds to the pink dotted line in the left panel of Figure 8). After the collision, each of the dark solitons travelled at the same velocity before the collision. Furthermore, we confirm that the phase shift occurred after the collision in the left panel of Figure 8.
Let us investigate the phase shift explicitly. By applying the Galilean transformation, that is, in Figure 9 we observe the scattering process in the inertial frame of reference moving with the left-hand-side notch of the quantum double darksoliton in Figure 8. We clearly confirm the phase shift after the collision as shown in    Figure 10 shows the time evolution of the square amplitude profile of the matrix element of the field operator for the quantum double dark-soliton states constructed with the Gaussian weights for c = 0.05. The quantum state is the same as that of Figure 8. We have constructed the double dark-solitons of distinct narrow notches with different depths not only in the density profile but also in the square amplitude profile, i.e., the graph of |ψ Q (x, t)| 2 versus x. We observe the scattering of two notches in the quantum double dark-soliton states. As the two notches of the double dark-soliton states approached each other, they moved along approximately straight and linear trajectories with different constant velocities, as shown in Figure 10. The collision occurred in the time interval including t = 11 (see panel (c), which corresponds to the pink dotted line in the left panel of Figure 10). After the collision, each of the dark solitons travelled with the same velocity before the collision. We observe at least approximately the same phase shift as shown in Figure 8. We remark that panel (a) of Figure 10, the square amplitude profile of the matrix element for the quantum states constructed with the Gaussian weights, corresponds to the panel of c = 0.05 in Figure 13 of Ref. [22].
We now demonstrate that the winding number changed during the scattering process in the time evolution of the Gaussian weighted quantum double dark-soliton states. Figure 11 shows the time evolution of the phase profile. In each panel, the phase profile satisfies the boundary condition: with a winding number J. At the initial time t = 0, the two notches of the quantum dark-soliton were located at the most distant points from each other such as X 1 = L/4 and X 2 = 3L/4, and the winding number is given by J = 1. When the two notches of the quantum dark-soliton states became very close in space, the winding number was suddenly changed to J = 0, in the time interval including t = 11, as shown in panel (c) of Figure 11. After the collision, the winding number was recovered: The winding number at t = 21 was given by J = 1, as shown in panel (d) of Figure 11. We remark that panel (a) of Figure 11, the phase profile of the matrix element between the Gaussian weighted quantum double dark-soliton states, corresponds to the panel of c = 0.05 in Figure 14 of Ref. [22].
We explicitly evaluate the phase shift due to the scattering of two notches in the quantum double dark-soliton state. The left panel of Figure 12 shows the square amplitude profile of the matrix element ψ Q (x, t) in time evolution observed in the inertial frame of reference moving together with the deeper notch of the quantum double darksoliton. The abrupt increase (or decrease) in the phase profile, which we call a phase jump, was located at the position of the deeper notch of the double dark-soliton, as shown in panels (a), (b), and (d) of Figure 12: It was located at x = 5 in panels (a) and (b), and at x = 2 in panel (d). Thus, the position of the deeper notch in the double dark-soliton was shifted after the collision in the inertial frame of reference. It corresponds to the phase shift due to the scattering of the two notches.
Let us investigate the changes of the winding number in time evolution in detail. The winding number J was equal to zero when the two notches of the quantum double dark-soliton were close to each other in space, as shown in panel (c) of Figure 12. Figure  13 exhibits that the abrupt changes of the winding number J from 1 to 0 and from 0   We recall that it is not necessary for the winding number in the phase profile of a quantum state to be conserved during the time evolution of the quantum system. The winding number is defined for the corresponding classical system, i.e., the GP equation, or for the phase profile of the quantum system. The dynamics of the quantum system can be much more complex than the solutions of the GP equation. When the two notches are far from each other in space, the phase profile of the quantum system is similar to that of the classical solution, while it is not the case when they collide with each other since they are very close in space.
In summary, the Gaussian weighted superposition of the two-hole excited states has lead to the quantum double dark-soliton states in which two notches have different depths [22].
It follows that the notches of the quantum double dark-soliton state have different velocities, and hence we have observed the scattering of two notches in the quantum double dark-soliton state exactly. We have also shown that the winding number of a quantum double dark-soliton state changed when the two notches approach each other, explicitly for the Gaussian weighted quantum double dark-soliton states.
We remark that one can make the quantum single dark-soliton black by making use of the Gaussian weights, as shown in Ref. [21]. However, for the quantum double dark-soliton, it seems that it is difficult to construct the double black-soliton only by applying the Gaussian weights to the superposition of a set of two-hole excitations.

Motivation to study the quasi-BEC in 1D for the ground state
In 1D systems quantum fluctuations play a key role and often give subtle and nontrivial effects. It is known that BEC occurs even for bosons with repulsive interactions due to the quantum statistical effect among identical particles [13]. In fact, the existence of BEC has been proven rigorously for interacting bosons confined in dimensions greater than one [42]. In 1D case there is no BEC for bosons with repulsive interactions due to strong quantum fluctuations if we assume the standard thermodynamic limit with fixed coupling constant [43]. On the other hand, if the coupling constant is very weak, we may expect that even the 1D bosons with a large but finite number of particles undergo a quasi-condensation in which "a macroscopic number of particles occupy a single one-particle state" [13]. We call it a quasi-BEC by following the Penrose and Onsager criterion.
However, it has not been shown explicitly how such a quasi-condensation occurs in interacting bosons in one dimension. Furthermore, it is nontrivial to expect it for the 1D Bose gas that is solvable by the Bethe ansatz. No pair of particles can have the same quasi-momentum in common for a Bethe-ansatz solution. Here we recall that we call the 1D system of bosons interacting with repulsive delta-function potentials the 1D Bose gas. For the impenetrable 1D Bose gas where the coupling constant is taken to infinity, condensate fractions are analytically and numerically studied [44], while in the weak coupling case it is nontrivial to evaluate the fractions in the 1D Bose gas.
We thus study in section 4 how the condensation fraction n 0 , i.e., the degree of the quasi-BEC, explicitly depends on the system size L, the number of particles N and the coupling constant c in the ground state of the LL model and particularly in the weak coupling case. It will be an illustrative example.

Onsager-Penrose criterion of BEC
Let us review the definition of BEC through the one-particle reduced density matrix for a quantum system [13,14]. We assume that the number of particles N is very large but finite. At zero temperature, the density matrix is given byρ = |λ λ|, where |λ denotes the ground state of the quantum system. We define the one-particle reduced density matrix by the partial trace of the density matrix with respect to other degrees of freedom:ρ 1 = N tr 23···Nρ . This matrix is positive definite and hence it is diagonalized asρ Here we put eigenvalues N j in descending order: N 0 ≥ N 1 ≥ N 2 ≥ · · · > 0. The sum of all the eigenvalues is given by the number of particles: j N j = N . Here we recall tr 1ρ1 = N due to the normalization: tr 123···Nρ = 1. Let us denote by n 0 the ratio of the largest eigenvalue N 0 to particle number N : The criterion of BEC due to Penrose and Onsager [14] is given as follows: If the largest eigenvalue N 0 is of order N , i.e., the ratio n 0 is nonzero and finite for large N , then we say that the system exhibits BEC, and we call n 0 the condensate fraction. Here we also define fractions n j by n j = N j /N for j = 1, 2, . . ..  Table 1. Fraction n sat of the reduced density operator at the origin, ρ 1 (0, 0), to the density n, evaluated by taking the sum over a large number of eigenstates |µ with one particle and one hole (1p1h) or with two particles and two holes (2p2h) for N = L = 50 (n = 1): Numerically we calculate correlation function in Equation (15) by taking the sum over a large number of eigenstates with one particle and one hole (1p1h) and those with two particles and two holes (2p2h). In order to confirm the validity of the restricted sum, we have estimated the ratio of the one-particle reduced density operator at the origin to density n, ρ 1 (0, 0)/n, through the form factor expansion in Equation (15) for the excitations with 1p1h or 2p2h. We express it by n sat . The estimates of n sat are listed in Table 1. The graph of n sat approaches 1 for small coupling constant c, while it is larger than 0.98 for any value of c in the case of N = 50.

4.4.
Evaluation of the one-particle reduced density matrix of the ground state For the LL model, the eigenfunctions of the one-particle reduced density matrix are given by plane waves for any nonzero and finite value of c. It is a consequence of the translational invariance of the Hamiltonian of the LL model. We thus have The eigenvalues of the one-particle reduced density matrix, N j , are expressed in terms of the form factor expansion. We consider the sum over all the form factors between the ground state, |λ , and such eigenstates, |µ , that have a given momentum P j as In the LL model we have P j := (2π/L)j. Solving the Bethe ansatz equations for a large number of eigenstates we observe numerically that eigenvalues N j are given in decreasing order with respect to integer j: N 0 > N 1 > N 2 > · · ·. It thus follows that condensate fraction which corresponds to the largest eigenvalue of the one-particle reduced density matrixρ 1 is indeed given by n 0 = N 0 /N , where N 0 has been defined by the sum of Equation (29) over all eigenstates with zero momentum.

Condensate fraction in the weak coupling regime
The estimates of condensate fraction n 0 are plotted against coupling constant c in the upper panel of Figure 14 over a wide range of c such as from c = 10 −3 to c = 10 3 for different values of particle number N such as N = 4, 10, . . . , 400. For each N , condensate fraction n 0 becomes 1.0 for small c such as c < 0.01, while it decreases with respect to c and approaches an asymptotic value in the large c region such as c > 100 or 1000. The asymptotic values depend on particle number N for N = 4, 10, . . . , 400, and they are consistent with the numerical estimates of occupation numbers for the impenetrable 1D Bose gas (see Equation (56) of Ref. [44]). In the lower panel of Figure  14, we plot fractions n j for j = 0, 1 and 2 against coupling constant c from c = 10 −3 to c = 10 3 with N = 20. The asymptotic values of n j for large c (i.e. c = 1000) are consistent with the numerical estimates for the impenetrable 1D Bose gas (for n 1 and n 2 , see Equations (57) and (58) of Ref. [44], respectively).
We observe that condensate fraction n 0 decreases as particle number N increases where density n = N/L is fixed. It is the case for c < 0.1 in the upper panel of Figure  14. Condensate fraction n 0 decreases as N increases even for small c such as c = 0.01, as shown in Figure 15. Thus, it is necessary for coupling constant c to decrease with respect to N so that condensate fraction n 0 remains constant as N increases with fixed density n.

Exact finite-size scaling
We now show the finite-size scaling of condensate fraction n 0 . In Figure 16 each contour line gives the graph of interaction parameter γ as a function of the inverse of particle number N for a fixed value of condensate fraction n 0 . They are plotted for various values of n 0 from n 0 = 0.6 to 0.99, and are obtained by solving the Bethe-ansatz equations numerically.
For different values of density such as n = 1, 2 and 5, we have plotted contour lines with fixed values of condensate fraction n 0 in the plane of interaction parameter γ  versus inverse particle number 1/N . We have observed that the contours with the same condensate fraction n 0 but for the different densities coincided with each other in the γ Thus, condensate fraction n 0 is constant as particle number N becomes very large if interaction parameter γ is given by the power of particle number N as in Equation (30). Applying the finite-size scaling arguments, we suggest from Equation (30) that condensation fraction n 0 is given by a scaling function φ(·) of a single variable γN η : n 0 = φ(γN η ). Here we recall the coincidence of contours for the different values of density n in Figure 16. We thus observe that exponent η and amplitude A of Equation (30) are determined only by condensate fraction n 0 and are independent of density n.
Let us consider amplitude A as a function of n 0 . We denote it by A = f (n 0 ). Then, the scaling function φ(·) is given by the inverse function: n 0 = f −1 (A). In Figure 17, exponent η increases with respect to n 0 , and amplitude A decreases monotonically with respect to n 0 .

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

Various limiting procedures
With the scaling behavior expressed in Equation (30) we derive various ways of the thermodynamic limit such that condensate fraction n 0 is constant. For instance, we consider the case of a finite particle number, N = N f . Choosing a value of n 0 , we determine γ by Equation (30) as γ = A(n 0 )/N η(n 0 ) f . Then, the 1D Bose gas with N = N f has the same condensate fraction n 0 for any large value of L if coupling constant c is given by c = A(n 0 )N 1−η f /L. Let us set η = 1 and N f = 10, for simplicity. We have n 0 = 0.97 in Figure 17, and γ = 0.3 at 1/N = 0.1 in the contour of n 0 = 0.97 in Figure 16. By assuming n = 1, it corresponds to the case of L = 10 and c = 0.3, and we have A = cL = 3, which is consistent with Figure 17. Therefore, the 1D Bose gas with N f = 10 has n 0 = 0.97 for any large L if c is given by c = 0.3/L. Moreover, we may consider other types of thermodynamic limits. When density n is proportional to a power of L as L α , condensate fraction n 0 is constant as L goes to infinity if we set c ∝ L (1−η)(1+α)−1 .
The scaling law in Equation (30) and the estimates of condensate fraction in the present paper should be useful for estimating conditions in experiments of trapped cold atomic gases in one dimension [45]. For instance, we suggest from Figure 14 that BEC may appear in 1D systems with a small number of bosons such as N = 20 or 40 for c = 1 or 10.

Concluding remarks
In the first part, we have shown that the density profile and the square amplitude evolved in time differently, in particular, for the equal weight case. In the former the notches were filled progressively, while the amplitude of the latter decreased gradually. Furthermore, the Gaussian weights led to the different depths for quantum double dark-solitons [22]. This gave the two notches of the quantum double dark-soliton the different speeds, and we observed the scattering of the two notches in the quantum double dark-soliton state exactly. Interestingly, the winding number of the quantum double dark-soliton state has changed when the two notches approach. Here we recall that it is not necessary for the winding number to be conserved in the time evolution of the quantum system, since it is defined for the corresponding classical system.
In the second part, we exactly calculated the condensate fraction of the 1D Bose gas with repulsive interaction by the form factor expansion for the ground state. We have shown the finite-size scaling behavior such that condensate fraction n 0 is given by a scaling function of interaction parameter γ times some power of particle number N : n 0 = φ(γN η ). Consequently, if parameter γ decrease as γ = A/N η , condensate fraction n 0 remains nonzero and constant as particle number N becomes very large. By modifying the thermodynamic limit, the 1D Bose gas shows BEC from the viewpoint of the Penrose-Onsager criterion.

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