Dark solitons revealed in Lieb-Liniger eigenstates

TL;DR

This paper demonstrates how dark solitons can be derived from the eigenstates of the Lieb-Liniger model of interacting bosons, unifying previous approaches and revealing solitonic structures through measurement-based analysis.

Contribution

It introduces a simple Ansatz that unifies prior methods and enables the study of large systems, showing solitons emerge from many-body eigenstates via measurements.

Findings

Dark solitons emerge in the mean-field limit from Lieb-Liniger eigenstates.

The Ansatz allows analysis of systems with hundreds of atoms.

Post-measurement states are wave packets of yrast states.

Abstract

We study how dark solitons, i.e. solutions of one-dimensional single-particle nonlinear time-dependent Schr\"odinger equation, emerge from eigenstates of a linear many-body model of contact interacting bosons moving on a ring, the Lieb-Liniger model. This long-standing problem was addressed by various groups, which presented different, seemingly unrelated, procedures to reveal the solitonic waves directly from the many-body model. Here, we propose a unification of these results using a simple Ansatz for the many-body eigenstate of the Lieb-Liniger model, which gives us access to systems of hundreds of atoms. In this approach, mean-field solitons emerge in a single-particle density through repeated measurements of particle positions in the Ansatz state. The post-measurement state turns out to be a wave packet of yrast states of the reduced system.

Equations

$$i\hslash {\partial }_t{\varphi }_{\mathrm{MF}}(x,t)=\left(-\frac{{\hslash }^2{\partial }_x^2}{2m}+g(N-1)|{\varphi }_{\mathrm{MF}}(x,t){|}^2\right){\varphi }_{\mathrm{MF}}(x,t),$$

$$\hat{H}=-\frac{{\hslash }^2}{2m}\sum _{j=1}^N{\partial }_{x_j}^2+g\sum _{\begin{array}{c}j,l\\ j<l\end{array}}^N\delta (x_j-x_l),$$

$$\psi (x_1+\Delta x,...,x_N+\Delta x)=e^{i\frac{2\pi }{L}K\Delta x}\psi (x_1,...,x_N).$$

$$\psi (x_1,\dots ,x_N)=\mathcal{N}{\int }_0^{L}dy{e}^{i\frac{2\pi }{L}Ky}\prod _{j=1}^N\varphi (x_j-y),$$

$${\psi }_{\mathrm{Ansatz}}(x_1,\dots ,x_N)=\mathcal{N}{\int }_0^{L}dy{e}^{i\frac{2\pi }{L}Ky}\prod _{j=1}^N{\varphi }_{\mathrm{MF}}(x_j-y)$$

$$|{\psi }_{\mathrm{Ansatz}}⟩=\mathcal{N}{\int }_0^{L}dy{e}^{i\frac{2\pi }{L}Ky}{e}^{-i\hat{P}y/\hslash }|{\varphi }_{\mathrm{MF}}{⟩}^{\otimes N}.$$

$$|{\psi }_{\mathrm{Ansatz}}⟩={\alpha }_0|K⟩+\sum _{j=1}^{\infty }{\alpha }_j|K^j⟩,$$

$$E_{\mathrm{Ansatz}}$$

$$|⟨{\psi }_{\mathrm{Ansatz}}|K⟩{|}^2\ge \frac{E_{\mathrm{exc}}^1-E_{\mathrm{Ansatz}}}{E_{\mathrm{exc}}^1-E_{\mathrm{yrast}}}=:F_B,$$

$$\hat{H}=\int_{0}^{L}\text{d}x\,\bigg{[}\frac{-\hbar^{2}}{2m}\hat{\Psi}^{\dagger}(x)\partial_{x}^{2}\hat{\Psi}(x)+\frac{g}{2}\hat{\Psi}^{\dagger}(x)^{2}\hat{\Psi}(x)^{2}\bigg{]},$$

$${\rho }_m(x)\propto$$

$$⟨K|{\hat{\Psi }}^{†}(x_1)...{\hat{\Psi }}^{†}(x_{m-1}){\hat{\Psi }}^{†}(x)\hat{\Psi }(x)\hat{\Psi }(x_{m-1})...\hat{\Psi }(x_1)|K⟩,$$

$${\mathrm{g}}_1(x)\propto$$

$$⟨{\hat{\Psi }}^{†}(x_1)...{\hat{\Psi }}^{†}(x_{m-1}){\hat{\Psi }}^{†}(x)\hat{\Psi }(0)\hat{\Psi }(x_{m-1})...\hat{\Psi }(x_1)⟩,$$

$$|{\varphi }_{\mathrm{MF}}{⟩}^{\otimes N}\approx \sum _{K^{\prime }}{a}_{K^{\prime }}|K^{\prime }{⟩}_N$$

$$|{\tilde{\psi }}^m⟩\propto \hat{\Psi }(x_m)\hat{\Psi }(x_{m-1})\dots \hat{\Psi }(x_1)|\psi ⟩.$$

$$|{\tilde{\psi }}_{\mathrm{Ansatz}}^m⟩\propto \hat{\Psi }(x_m)\hat{\Psi }(x_{m-1})\dots \hat{\Psi }(x_1)|{\psi }_{\mathrm{Ansatz}}⟩\propto$$

$${\int }_0^{L}dy{e}^{i\frac{2\pi }{L}Ky}\left(\prod _{j=1}^m{\varphi }_{\mathrm{MF}}(x_j-y)\right)e^{-i\hat{P}y/\hslash }|{\varphi }_{\mathrm{MF}}{⟩}^{\otimes (N-m)}.$$

$$|{\tilde{\psi }}_{\mathrm{Ansatz}}^m⟩=\sum _{K^{\prime }}{a}_{K^{\prime }}|K^{\prime }{⟩}_{N-m}+\sum _j{b}_j|{\psi }_j{⟩}_{N-m}$$

$$\rho (x)=|{\varphi }_{\mathrm{MF}}(x){|}^2,$$

$$\phi (x)=\mathrm{Arg}\left\{{\varphi }_{\mathrm{MF}}(x)\right\},$$

$$\rho (x)$$

$$\phi (x)$$

$$v$$

$$a_1$$

$$a_2$$

$$a_3$$

$${\phi }_{\mathrm{black}}(x):=\pi \left(\mathrm{sgn}\left(x-L/2\right)-x/L\right),$$

$$\tilde{\varphi }(x):=e^{-\hat{T}\delta t}{e}^{-\hat{V}\delta t}\varphi (x),$$

$$\tilde{\varphi }(x)$$

$$\widetilde{\stackrel{~}{\varphi }}(x):=\frac{\tilde{\varphi }(x)}{\int \mathrm{d}x|\tilde{\varphi }(x){|}^2}.$$

$$\varphi (x)=|\widetilde{\stackrel{~}{\varphi }}(x)|e^{i{\phi }_{\mathrm{black}}(x)}.$$

$${\psi }_{\mathrm{Ansatz}}(x_1,\dots ,x_N)=\mathcal{N}{\int }_0^{L}dy{e}^{i\frac{2\pi }{L}Ky}\prod _{j=1}^N{\varphi }_{\mathrm{MF}}(x_j-y),$$

$$⟨\hat{p}⟩:=-i\hslash {\int }_0^{L}dx{\varphi }_{\mathrm{MF}}^{\ast }(x){\partial }_x{\varphi }_{\mathrm{MF}}(x)=\frac{2\pi \hslash }{L}K/N.$$

Full text

I Introduction

The famous Lieb-Liniger (LL) model Lieb (1963); Lieb and Liniger (1963) describes particles moving along a circle and interacting via delta inter-atomic potential. Such a simple interaction turns out to be a well-suited approximation for realistic interactions between neutral slow atoms. Thus, the LL model and its extensions, remain an active research topic in theoretical and experimental physics and mathematics Gogolin et al. (2004); Jiang et al. (2015); Lang et al. (2017); Cazalilla et al. (2011); Sowiński and García-March (2019).

The same system, of $N$ atoms with contact interaction, is often treated within a simple mean-field (MF) approximation, based on the non-linear Schrödinger equation (NLSE):

[TABLE]

where the wave function $\phi_{\rm MF}(x,t)$ is interpreted as an orbital occupied by a macroscopic number of atoms, $m$ is the particle mass, and $g$ is the interaction strength. The latter equation (1), is useful in many areas of physics ranging from quantum optics Kivshar and Luther-Davies (1998) to hydrodynamics Gross (1961); Johnson (1976). It is also a rare example of a model with physical applications supporting solitonic solutions Zakharov and Shabat (1973), observed in atomic gases Burger et al. (1999), plasma Heidemann et al. (2009), water waves Chabchoub et al. (2013), ferromagnetic materials Tong et al. (2010). In the case of gases, when their atoms repel each other, i.e. $g>0$, a soliton is a rarefaction in the atomic density, which moves with a constant speed, preserves its shape and is unusually robust thanks to the balance between dispersion and nonlinearity Frantzeskakis (2010). In this case ($g>0$) the soliton is called a dark soliton, which can be either black or gray. Black solitons are characterized by a zero-density dip i.e. a point where the atomic density is exactly zero and the phase of $\phi_{\rm MF}(x)$ undergoes a sharp $\pi$-jump. While gray solitons have a non-zero density dip with the phase jump strictly smaller than $\pi$.

There is a puzzling link between the MF dark solitons and the solutions of the underlying many-body LL model. More than a decade after the seminal paper by E. Lieb Lieb (1963), a coincidence between the dispersion relations of dark solitons and certain many-body eigenstates, the so-called type-II elementary excitations, was observed in the weak interaction limit Kulish et al. (1976); Ishikawa and Takayama (1980). The type-II excitations are simply the many-body eigenstates that minimize the energy for a fixed total momentum, sometimes called yrast states Mottelson (1999). Together with the type-I excitations (corresponding to Bogoliubov quasiparticles Lieb (1963)) they constitute two branches of elementary excitations, with dispersion relations sketched in the left panel of Fig. 1.

Why is the correspondence between the yrast states and dark solitons bizarre? Firstly, yrast states, as eigenstates of the system, are stationary solutions, whereas the MF dark solitons are solutions of time-dependent Eq. (1). Moreover, the LL model includes all correlations between particles in a linear Hamiltonian, while the MF approach gets rid of mutual correlations but introduces the non-linearity in the description. Finally, the eigenstates of the many-body model have to be translationally invariant. On the other hand, the dark solitons do not obey this symmetry. As the nature of the MF solitons and the yrast states is so different, the question arises whether they have something in common except the same dispersion relation.

There have been efforts to show, that the relation between these objects is deeper, and in particular that the MF soliton can be extracted directly from the yrast states. In Refs. Syrwid and Sacha (2015); Syrwid et al. (2016) it was indicated that MF solitons are already hidden in a single yrast state, and they will emerge in sufficiently high order correlation functions. In turn, in Refs. Sato et al. (2012, 2016); Kaminishi et al. (2018); Shamailov and Brand (2019) it was shown that MF soliton appears in a single-particle density, when calculated for an appropriate superposition of many yrast states. The other relations between the many-body states and solitons were also presented in Martin and Ruostekoski (2010); Katsimiga et al. (2017a, 2018); Kaminishi et al. (2018); Delande and Sacha (2014). Still, these interesting results leave the field in an unpleasant situation of many seemingly unrelated views on the connection between MF solitons and many-body yrast states.

Here, we unify different approaches by employing a simple but powerful Ansatz for the yrast state of the LL model in the MF regime of parameters. We use this Ansatz, to show that the state which appears after "measuring" many particles drawn from a high order correlation function is a random superposition of yrast states. The mutual unification between two approaches becomes apparent as we consider systems consisting of hundreds of atoms close to the MF regime.

The paper is organized as follows. In Section II we remind LL model, yrast state and define the parameter regime we are interested in. Our Ansatz for yrast states is introduced in Section III. Validity of this Ansatz is discussed in Section IV. In Sections V and VI we show how the different constructions of the MF solitons out of the yrast states presented in Ref. Syrwid and Sacha (2015) and Refs. Sato et al. (2012, 2016); Kaminishi et al. (2018) prove to be different views of the same object. To make this paper self-complete we describe in details all relevant analytical previous results and our numerical approaches in Appendices.

II The Lieb-Liniger model and yrast states

We will investigate eigenstates of the many-body system of $N$ bosons moving along a circle of length $L$ governed by the LL Hamiltonian 111We rescaled the original Lieb-Liniger model by a factor $2$, to have the form which is used more frequently now.:

[TABLE]

where $x_{j}$ denotes position of $j$-th particle. As this system is translationally invariant, the values of total momentum $\hat{P}=-i\hbar\,\sum_{j=1}^{N}\partial_{x_{j}}$ can be used to label the energy eigenstates, even in the case with interaction. The exact solutions for the eigenstates of the LL Hamiltonian (2) for repelling ($g>0$) particles are known since 1963 Lieb (1963). Among the eigenstates, there are special ones that are called elementary excitations, which can be divided into two families. Before Ref. Lieb (1963), the approximated theories was applied to find energies of weak perturbations of an atomic gas. This originated in a single family of the Bogoliubov elementary excitations (identified with the type-I excitations). Its dispersion relation given by $E_{\rm B}(P)=\sqrt{\frac{P^{2}}{2m}\left(\frac{P^{2}}{2m}+\frac{2gN}{L}\right)}$ DeWitt (1958).

The unexpected second family, revelead by the exact solution presented in Lieb (1963), consists of the aforementioned yrast states. These are also called "one hole excitations" Lieb (1963), or the "lowest energy solutions for fixed total momentum" Mottelson (1999). The yrast state with the total momentum of $P=\frac{2\pi\hbar}{L}K$ is represented by $|K\rangle$, where $K$ is an integer.

As the subject of this paper is the relationship between the yrast states and solitons, we will restrict our considerations to MF regime as the NLSE should work, in principle, only there. That is, within the regime of weak interactions with only slightly correlated atoms, in which quantum phenomena, like the quantum depletion of the ground state, are small. On the other hand, it is desired to see the effects that are substantially different from the ideal gas case. Therefore, we require the healing length $\xi:=\sqrt{\hbar^{2}L/gmN}$, which is close to the soliton width, be much shorter than the system size $L$ Ołdziejewski et al. (2018). If we were to consider a small number of atoms, the latter condition would lead to a large $g$, resulting in strong interactions. Therefore, we are interested in the limit in which number of atoms $N$ converges to infinity, an interaction strength $g$ goes to [math], but the MF parameter $ng:=\frac{\hbar^{2}}{m\xi^{2}}$, with $n:=N/L$ denoting a gas density, is fixed. In this limit the LL coupling constant $\gamma\propto g/N$ Lieb (1963) decreases with the number of atoms as $(ng)/N^{2}$, indicating that indeed the system enters quickly a weak interaction regime.

The MF regime defined in such a way is very difficult to handle in the frame of many-body analysis, which is usually limited to systems with small number of atoms $N$. Apart from the few existing semi-analytical results Sato et al. (2012); Shamailov and Brand (2019); Kaminishi et al. (2018), the majority of approaches are devoted to small systems of $\approx 10-20$ atoms Kanamoto et al. (2008); Kaminishi et al. (2011); Fialko et al. (2012); Syrwid and Sacha (2015); Syrwid et al. (2016); Ołdziejewski et al. (2018) solved with brute force methods or around $\approx 100$ atoms solved with sophisticated and time consuming numerics Katsimiga et al. (2017b); Delande and Sacha (2014); Mistakidis et al. (2018); Katsimiga et al. (2018). Our way around these numerical difficulties is to use a natural and simple Ansatz for the yrast states in the MF regime.

III The Ansatz for yrast states

Here we shall discuss, step-by-step, our construction of the Ansatz for the yrast states in the MF regime. The main building block of the Ansatz consists of a product state of $N$ particles occupying a single orbital $\phi(x)$ represented by $\prod_{j=1}^{N}\phi(x_{j})$. Next, the Ansatz has to belong to the same momentum space as the yrast state $|K\rangle$, i.e. the translation of the $N$-particle wave function by $\Delta x$ needs to be equivalent to multiplying it by a factor $e^{i\frac{2\pi}{L}K\Delta x}$:

[TABLE]

For any orbital $\phi$ one may construct the states satisfying above condition by taking a continuous superposition of product states shifted by the translation operator $e^{-i\hat{P}y/\hbar}$ and multiplied by the phase factor $e^{i\frac{2\pi}{L}Ky}$ over all possible shifts $y$ (see Appendix B.1 for formal justifications):

[TABLE]

where $\mathcal{N}$ is a (real) normalization factor.

As a yrast state $|K\rangle$ is the lowest energy state for fixed total momentum equal exactly to $\frac{2\pi\hbar K}{L}$, therefore the task would be to find an orbital ${\phi}$, that minimizes the average energy of the state (4). Finding such orbital would be a difficult task, as the average energy of the state (4) is given by a complicated formula (see Appendix B.4). On the other hand it is known that energies of the yrast states $|K\rangle$ and MF solitons agree Kulish et al. (1976). Therefore, as the Ansatz for yrast state we choose the state (4) with ${\phi}(x)={\phi}_{\rm MF}(x)$, where ${\phi}_{\rm MF}(x-vt)$ is the solitonic solution of the NLSE with the average single particle momentum $\left\langle-i\hbar\,\partial_{x}\right\rangle$ equal to $2\pi\hbar\,K/(NL)$:

[TABLE]

We also use state defined in Eq. (5) in the Dirac notation:

[TABLE]

The exact form of solitonic solution on the circle ${\phi}_{\rm MF}$ is quite complicated – we give the appropriate formulas and our numerical methods for handling them, in Appendix A.

The construction of the Ansatz is sketched in Fig. 1 – the many-body Ansatz is understood as a continuous superposition of macroscopically occupied MF solitons. Each soliton in the superposition (5) appears with a phase prefactor $e^{i2\pi Ky/L}$ (distinguished in Fig. 1 with a color) depending on the position shift $y$.

We remark that the Ansatz follows the ideas partially spread in the community, that the yrast states should be somehow related to the product states of MF solitons but with unknown position of the density dips, i.e. smeared over the whole circle as shown in Fig. 1. Such Ansatz was presented in the context of Bose-Einstein condensation Castin and Herzog (2001). Condensate, as defined via Penrose-Onsager criterion Penrose and Onsager (1956), is supposed to appear in the system of bosons at very low temperature. As show in Ref. Castin and Herzog (2001), surprisingly, even at $T=0$, when the system is in a ground state, there may be no condensation at all. This happens when the system has some continuous symmetry, like the translational invariance in our case. Once the symmetry is broken, for instance by measuring positions of few bosons, condensation may emerge immediately. The situation presented in this paper is similar but the resulting condensate is (i) not a ground state of the system and (ii) is temporary as it may disappear in time Syrwid et al. (2016); Sato et al. (2016).

How accurately does the energy of the Ansatz agree with the energy of the exact LL solution, the yrast state, as the interaction strength increasing? In Fig. 2, we present the energies as functions of the MF parameter $ng$, where $n=N/L$ is the gas density, for the yrast state with $K=N/2$ (black soliton) and $K=N/4$ (gray soliton) and the corresponding Ansatz (6) (see Appendices B.4, B.6 and C for the details of computations). For the reference, we plot results obtained via the first order perturbation theory with interaction strength $g$ being a small parameter (dashed burgundy line). In this case the average energy is evaluated in the yrast state corresponding to $g=0$, which is a state with $N-K$ motionless particles and the remaining $K$ of them with momentum $2\pi\hbar/L$. Its average energy is a linear function of interaction strength $g$ equal to $K\frac{4\pi^{2}\hbar^{2}}{2mL^{2}}+\frac{g}{2L}\left(N(N-1)+2NK-2K^{2}\right)$. The second reference curve, is the average energy of the MF state, with $N$ atoms occupying a single orbital $\phi_{\rm MF}$ (solid gray line).

As expected the energy of yrast state, MF state, and Ansatz are close to each other, even for such strong interactions that first order perturbation theory fails. Moreover, in the ideal gas limit, our Ansatz (6) is exactly equal to the yrast state $\lim_{g\to 0}|\psi_{\rm Ansatz}\rangle\equiv\lim_{g\to 0}|K\rangle$ for any $K\neq 0$ Ołdziejewski et al. (2018). In the same limit, the energy of the dark MF soliton is actually slightly smaller than the energy of the yrast state (see Appendix B.2).

The intuitive definition of an Ansatz, together with the apparent agreement between its energy and the energy of yrast state, motivate us to use the Ansatz (6) instead of yrast state in our study on the correspondence between LL model and NLSE. Before we shall do it, we discuss in more detail the validity range of the Ansatz.

IV Validity range

We analyse the validity range by determination how well our Ansatz approximates the true yrast state for a given strength of interaction. The most objective and unambiguous measure of the similarity between two states is fidelity $|\langle\psi_{\rm Ansatz}|K\rangle|^{2}$. However, to calculate it directly one would need to express the yrast state in the position representation (or to express the Ansatz in terms of quasi momenta from solutions of the LL model), which would be an extremely demanding task. Therefore we propose a simple (but very rough) lower bound for the value of fidelity, based on the values of energies calculated independently in each of the formalisms mentioned.

As the Ansatz has a well-defined total momentum, therefore it can be decomposed into the basis of many-body eigenstates of the same total momentum:

[TABLE]

where $|K^{j}\rangle$ is the $j$-th excited eigenstate with total momentum $\frac{2\pi\hbar}{L}K$ and energy $E_{\rm exc}^{j}$, and $\sum_{j=0}^{\infty}|\alpha_{j}|^{2}=1$. The excited eigenstates are listed in the ascending order with respect to their energy, i.e. $E_{\rm exc}^{1}\leq E_{\rm exc}^{2}\leq E_{\rm exc}^{3}\leq\ldots$.

The average energy $\langle\psi_{\rm Ansatz}|\hat{H}|\psi_{\rm Ansatz}\rangle$ can be expressed with the help of Eq. (7) as:

[TABLE]

where we use an inequality $E_{\rm exc}^{j}>E_{\rm exc}^{1}$ and $E_{\rm yrast}$ is the energy of the yrast state with the total momentum $\frac{2\pi\hbar}{L}K$. Therefore, the fidelity between the Ansatz (6) and the corresponding yrast state, i.e. $|\alpha_{0}|^{2}$, obeys the inequality

[TABLE]

where we introduced a fidelity bound $F_{B}$.

In Fig. 3 we show how $F_{B}$ decreases with increasing interaction strength. First of all, we observe, that data points calculated for a different number of atoms, but for the same value of mean-field parameter $ng$, are close to each other. Second, we see that for small values of $ng$ the fidelity between the Ansatz and the yrast state has to be very high (close to one). Here, it is worth stressing the fact, that the relatively small value of $F_{B}$ for much bigger $ng$ does not automatically implicate uselessness of the Ansatz for stronger interactions – the fidelity may be still close to one, but calculating it directly would require much more advanced numerical methods and would be unfeasible for large number of atoms.

Given that for $ng\lesssim 25$ the Ansatz is a good approximation for the yrast state, we proceed to use it as a replacement for the yrast state to unify the results of other groups. We will benefit from the fact, that for the Ansatz many calculations may be done analytically and the remaining necessary numerical analysis is feasible even for large number of atoms.

V Dark solitons revealed in high order correlation functions

In this section, we study the emergence of dark solitons out of a single yrast state as it was done in Ref. Syrwid and Sacha (2015). To make the comparisons with the literature results easier, we use the second quantization formalism, in which the energy operator (2) may be written as

[TABLE]

where $\hat{\Psi}(x)$ ($\hat{\Psi}^{\dagger}(x)$) is a annihilation (creation) field operator of a boson at position $x$ satisfying the commutation relations $[\hat{\Psi}(x),\hat{\Psi}^{\dagger}(x^{\prime})]=\delta(x-x^{\prime}),\,[\hat{\Psi}(x),\hat{\Psi}(x^{\prime})]=[\hat{\Psi}^{\dagger}(x),\hat{\Psi}^{\dagger}(x^{\prime})]=0$. The second quantization formalism is also very handy in performing any computations within the Ansatz (see Appendix B).

The object of interest in Ref. Syrwid and Sacha (2015) is a $m$-th order correlation function

[TABLE]

normalized to $1$. From its mathematical structure, Eq. (11) is identified as a probability density function (PDF) from which one draws a random position $x_{m}$ of the $m$-th particle to be measured. In Ref. Syrwid and Sacha (2015), function (11) is considered for increasing $m$ after subsequent "measurements" of particles. It was observed that $\rho_{m}$ resembles density $|\phi_{\rm MF}|^{2}$ of the MF soliton from the NLSE (1). However, in Ref. Syrwid and Sacha (2015) the agreement between MF solitons and the $m$-th order correlation function of a yrast state was demonstrated only for the healing length $\xi>L$ and for a few particles, namely for almost a non-interacting system. Contrarily, in the present paper, the Ansatz enables the study of healing lengths much smaller than the size of the system $L$ for a large number of particles. To that end, we employ the following procedure using Eq. (11) with the yrast state $|K\rangle$ replaced by our Ansatz $|\psi_{\rm Ansatz}\rangle$ (6). We begin with a computation of the single-particle reduced density matrix $\rho_{1}$ which is used as a PDF to draw a random position $x_{1}$. Subsequently, we compute second-order correlation function $\rho_{2}$ which again serves as the PDF for the next random position $x_{2}$ draw. Repeating this process $m-1$ times outputs the $m-1$ positions, parameters of the marginal distribution $\rho_{m}(x)$ (11), evaluated for the Ansatz (6) (for details of our calculations see Appendix B.3).

In the left panel of Fig. 4, we show samples of correlation functions $\rho_{m}(x)$ of different orders $m$ calculated for the parameters $ng=25$, $N=500$, and $K=50,\,125,\,250$. Densities $\left|\phi_{\rm MF}(x)\right|^{2}$ of the solitonic solution of Eq. (1) are shifted so that their notch positions overlap with that of $\rho_{150}(x)$ for direct comparison. We observed that the notch position of $\rho_{m}(x)$ is determined early i.e. for low order correlation functions and stabilizes as $m$ increases, with slight fluctuations dependent on the random particle position draws. Even highly disruptive particle measurements caused by unlikely draws, as that exemplified by $\rho_{50}(x)$ for $K=100$ in Fig. 4, do not prevent the formation of a dark soliton. It is important to mention that every curve presented in Fig. 4 is a result obtained for a single simulation. In all simulations we have made, the high order correlation functions always resemble the density of a MF soliton. Note, that our result corresponds to the short healing length $\xi=0.2L$. We observe that the MF soliton emerges from the correlation function as $m$ increases. We also notice that obtaining a very good agreement between the density emerging from the many-body calculation and the MF soliton requires calculation of high order correlation function.

The agreement between the MF approach and the Ansatz encourages us to investigate the spatial phase, which is peculiar for dark solitons. A phase ${\rm arg}\left\{\phi_{\rm MF}(x)\right\}$ of the soliton changes quickly within the density notch, but remains linear far from it (see the thick gray line in the right panel of Fig. 4). Therefore, when the system is in the solitonic state, the majority of atoms moves along the circle with a constant velocity, apart from the place of rarefaction where particles move quickly in the opposite direction. We extract the phase of the many-body wave function using

[TABLE]

evaluated for the Ansatz (6) (for details of our calculations see Appendix B.3). As shown in the right panel of Fig. 4, the phase of $\text{g}_{1}$ converges to $\arg\left\{\phi_{\rm MF}(x)\right\}$ for increasing $m$. The phase of MF soliton is shifted by the same amount as the corresponding density $\left|\phi_{\rm MF}(x)\right|^{2}$ in the left panel.

The above results prove that, indeed, the MF solitons emerge in the high order correlation function evaluated for the Ansatz (6) in the true MF regime with $\xi$ significantly smaller than $L$.

On the other hand, in Refs Sato et al. (2012, 2016); Kaminishi et al. (2018) the MF solitons were constructed from the many-body eigenstates of the LL Hamiltonian (2) in a completely different way, as explain in the next section.

VI Dark solitons as superpositions of yrast states

In the previous section, we have shown using our Ansatz that the MF solitons emerge in high order correlation functions. An opposite direction was taken in Refs. Sato et al. (2012, 2016); Kaminishi et al. (2018); Shamailov and Brand (2019) where the dark solitons are constructed as a specific superposition of yrast states. Namely the MF product state is expressed as:

[TABLE]

where $a_{K^{\prime}}$ are expansion coefficients ( drawn from a chosen distribution) and $|K^{\prime}\rangle_{N}$ denotes a yrast state of the system with $N$ particles and the total momentum $\frac{2\pi\hbar}{L}K^{\prime}$. A comprehensive discussion of different $a_{K^{\prime}}$ choices can be found in Refs. Sato et al. (2012, 2016); Kaminishi et al. (2018); Shamailov and Brand (2019).

An interesting question arises whether these two approaches of linking the yrast states of the LL model with the dark solitons from the NLSE complement each other or are completely unrelated. We shall answer this question by appealing to the definition of the high order correlation function $\rho_{m}$ (11) and using the Ansatz (6) for the yrast state.

Calculation of any correlation function by means of the second quantization requires the sequential action of the annihilation field operators at some points in space. One can say that such procedure conditions a system’s wave function. Physically, it corresponds to an instantaneous destructive measurement of certain particle positions. Thus, we introduce a conditional wave function $|\tilde{\psi}^{m}\rangle$ of the system in a state $|\psi\rangle$ after measuring positions of $m$ particles given by

[TABLE]

To maximize the reliability of such a measurement in any theoretical considerations it has to be performed according to a multivariate probability distribution determined by the wave function for a given state. Therefore, each position $x_{i}$ from Eq. (14) should be taken from the particular PDF $\rho_{i}$ defined in Eq. (11).

The average density in the conditional state (14) $\rho(x):=\left\langle\hat{\Psi}^{\dagger}(x)\hat{\Psi}(x)\right\rangle$ is equal to the $\left(m+1\right)$-th order correlation function $\rho_{m+1}$ (11) studied in the previous section. Therefore, to bridge the different views on the correspondence between the MF solitons and yrast states, one has to verify whether the conditional wave function (14) for $|\psi\rangle$ being a yrast state can be represented as a wave packet of yrast states, each with $N-m$ atoms and different total momentum. As calculations with the help of the exact many-body states would be limited to a small number of atoms, we again refer to the family of Ansatzes (6) as an approximation for the yrast states with different total momenta $\frac{2\pi\hbar}{L}K^{\prime}$.

For the Ansatz (6) one can easily find the conditional state 222For details of our calculations see Appendix B.3.

[TABLE]

Due to the factors $\phi_{\rm MF}(x_{j}-y)$, the solitons centered close to the positions $x_{j}$, where a measurement occurred, enter the conditional state with lower weights, as compared to solitons with a density dip far from $x_{j}$. The conditional wave function $|\tilde{\psi}^{m}_{\rm Ansatz}\rangle$ is no longer an eigenstate of LL system with $N-m$ particles as it is not translationally invariant. However, we can always decompose this state into the set of eigenstates of LL model for $N-m$ particles in the following way:

[TABLE]

where $|\psi_{j}\rangle_{N-m}$ is an eigenstate of the system, which is not the yrast state, and $\sum_{K^{\prime}}|a_{K^{\prime}}|^{2}+\sum_{j}|b_{j}|^{2}=1$. The question is whether the conditional state $|\tilde{\psi}^{m}_{\rm Ansatz}\rangle$ remains in the subspace of the yrast states in the form of a wave packet. To answer the question we calculate the overlap between the Ansatz (6) and conditional state (15) finding the weight of the yrast subspace given by $\sum_{K^{\prime}}|a_{K^{\prime}}|^{2}$ and the $a_{K}^{\prime}$ distribution 333 For details of our computations see Appendix B.5.

In the top panel of Fig. 5, we present sum of weights $\sum_{K^{\prime}}|a_{K^{\prime}}|^{2}$ of the yrast subspace as a function of a measured number of atoms $m$ for $ng=25$, $N=100,\,250,\,500,\ 1000$ and a fixed total momentum $K=N/10$ (left) and $K=N/4$ (right). We observe that the greater the number of atoms $N$ the closer to one the weight of the yrast subspace for a given value of $m$ is. It means that, indeed, the conditional wave function, which reveals the dark soliton, is approximated by a superposition of the yrast states. In the bottom panel of Fig. 5, we plot five representative distributions of $a_{K^{\prime}}$ as functions of total momentum $\frac{2\pi\hbar}{L}K^{\prime}$ for $ng=25$, $N=1000$, $K=N/10$, $m=10$ (left) and $m=40$ (right). The resulting distributions differ from shot-to-shot but they give the same single-particle density. Note that this agrees with the existing literature where different distribution models were considered.

Our efforts in bridging the dark solitons and the yrast states also result in the unification of the previous attempts done in the literature Sato et al. (2012, 2016); Kaminishi et al. (2018); Shamailov and Brand (2019). In this section we have shown that the dark solitons hosted in the non-ideal gas of $N$ atoms and revealed by partial measurements can be almost exactly expressed as a wave packet of the yrast states of a gas of $N-m$ atoms.

VII Conclusions

We studied the correspondence between the yrast states of the Lieb-Liniger Hamiltonian and the mean-field solitons from the non-linear Schrödinger equation. To this end we proposed a simple construction for the yrast state (5) based on mean-field product states with appropriate phase factors. Using this Ansatz we were able to unify previous literature results and observations Syrwid and Sacha (2015); Sato et al. (2012, 2016); Kaminishi et al. (2018); Shamailov and Brand (2019) about the subject at hand.

The conditional wave function, which results from annihilation of $m$ particles in the Ansatz state at random positions, reveals the ultimate utility of our approach. The single-particle density evaluated in the conditional wave function is the $m$-th order correlation function resembling the mean-field soliton, as discussed in Ref. Syrwid and Sacha (2015). Moreover, the conditional wave function is found to be a wave packet of yrast states of $N-m$ atom system with different total momenta, as analysed in Refs. Sato et al. (2012, 2016); Kaminishi et al. (2018); Shamailov and Brand (2019). As can be readily seen, our proposal complements various preceding studies reproducing their results with a singular construction and thus tying them into a single picture.

The "measurements" needed to break the translational symmetry could be realized spontaneously, due to particle losses which are inevitable in the ultracold gases. We plan to study this in detail, using many-body methods. Another remaining question concerns dynamical stability of the conditional wave function. In all previous works, the emerging solitonic profiles were, unlike the mean-field solitons, blurred during evolution Syrwid et al. (2016); Sato et al. (2012); Kaminishi et al. (2018); Shamailov and Brand (2019). It was argued that the time of blurring should increase to infinity in the thermodynamic limit Syrwid et al. (2016); Sato et al. (2016). However, this hypothesis has to be verified. This can be done again employing our Ansatz.

Acknowledgements.

We acknowledge fruitful discussions with K. Rzążewski, M. Gajda and K. Sacha. In particular we thank to M. Gajda for convicing us, that the results concerning the unification of different views on relations between yrast states and MF solitons are the most interesting aspects of our work. This work was supported by the (Polish) National Science Center Grants 2016/21/N/ST2/03432 (R.O.), 2015/19/B/ST2/02820 (W. Górecki), 2014/13/D/ST2/01883 (W. Golletz and K.P.) and 2018/31/B/ST2/00349 (W. Golletz).

Appendix A Mean-field solitons

A.1 Mean-field gray solitons

The solitonic solution $\phi_{\rm MF}(x)$ of NLSE (1) was discussed several times in literature Zakharov and Shabat (1973); Carr et al. (2000); Sato et al. (2016). Here we briefly present the final formulas, in the form which was used to produce results of this paper.

We followed a procedure described in Sato et al. (2016). Solitonic solution of NLSE is a running wave, $\phi_{\rm MF}(x,\;t)=\phi_{\rm MF}(x-vt)$ with speed $v$. In what follows we will omit the time dependence and give separately a solitonic density:

[TABLE]

and its phase

[TABLE]

and its speed $v$. The approach presented in Sato et al. (2016) was devoted to the case of gray solitons, i.e. when $\varphi(x)$ is continuous and $\rho(x)$ is always larger than [math]. The density and the phase and the velocity of the dark soliton obeying NLSE, in terms of four parameters denoted with $a_{1}$, $a_{2}$, $a_{3}$ and $k$, are given by:

[TABLE]

where ${\rm sn}(u,\,k)$ is the Jacobi elliptic function, $\Pi$ is the incomplete elliptic integral of the third kind, with the modulus of Jacobi’s elliptic function $k$ and ${\rm am}$ is the Jacobi amplitude. Parameter $a_{1}$ has simple physical interpretation — it is the minimal density in the solitonic solution. Periodic boundary conditions for the phase and density and normalization condition $\int|\phi_{\rm MF}|^{2}=1$ lead to the following relations between parameters $a_{1}$, $a_{2}$, $a_{3}$ and the elliptic modulus $k$:

[TABLE]

where $K(k)$ and $E(k)$ are the elliptic integrals of the first and the second kind, respectively.

To compute phase and density of a soliton that has a desired minimum of the density $a_{1}$, we first solve numerically Eq. (20) for the elliptic modulus $k$, and then we use Eqs. (21) and (22) to find the remaining parameters $a_{2}$ and $a_{3}$. Having determined $a_{1}$, $a_{2}$, $a_{3}$ and $k$ we can compute the soliton wave function $\phi_{\rm MF}(x)$ at any position $x$, using Eqs. (19) and relation $\phi_{\rm MF}(x)=\sqrt{\rho(x)}e^{i\varphi(x)}$. The average momentum $\langle\hat{p}\rangle:=\frac{\hbar}{i}\int\,{\rm d}x\,\phi_{\rm MF}^{*}(x)\partial_{x}\phi_{\rm MF}^{*}(x)$ is then computed numerically.

To find the MF soliton with target momentum $p_{\rm target}=2\pi\hbar K/(NL)$, we repeat the steps described above varying $a_{1}$, until the numerically determined momentum matches $p_{\rm target}$. After a few bisection steps with respect to $a_{1}$ the target MF soliton is found – we stop bisection when relative discrepancy between numerically computed momentum and the target momentum is below $10^{-6}$.

A.2 Mean-field black solitons

The formulas presented in the previous section are derived under the assumption that the density $\rho(x)$ given in Eq. (19) is always grater than [math]. Therefore, they can not be used in the case of a black soliton. In fact, already for a gray but very deep solitons, the equation (20) becomes very demanding, as discussed in Sato et al. (2016) (see for instance Table 4 in Sato et al. (2016) for the values of the parameter $k$).

On the other hand, one can use the properties of the black soliton to quickly find it numerically with other method. The trick is to compute the lowest energy state of NLSE (1), but in the space of functions with a given phase. We have learned this trick from Tomasz Karpiuk 444Private communication., and used it successfully in Pawłowski and Rzażewski (2015). Precisely, we look for solutions with the phase:

[TABLE]

where the signum function ${\rm sgn}$ is equal to $1$ for positive arguments and [math] for negative ones. The signum function introduces a discontinuity in the phase, a $\pi$ jump, which is the characteristic feature of a black soliton.

The minimal energy state in the space of functions with phase $\varphi_{\rm black}(x)$ is found with the split-step imaginary time evolution, implemented as follows:

We start with an arbitrary function $\phi(x)$.

• 2)

We compute $\tilde{\phi}(x)$ according to the split-step formula:

[TABLE]

where $\hat{V}=g|\phi(x)|^{2}$ and $\hat{T}=-\frac{\hbar^{2}}{2m}\partial_{x}^{2}$. To act with the operator $e^{-\hat{T}\,\delta t}$ we apply Fourier transform $\mathcal{F}$ and its inverse $\mathcal{F}^{-1}$:

[TABLE]

to replace the cumbersome operator $e^{(\hbar^{2}\delta t/2m)\partial^{2}_{x}}$ with its Fourier representation.

• 3)

We normalize the output of the previous step:

[TABLE]

• 4)

We define a new function $\phi(x)$ as $\tilde{\tilde{\phi}}(x)$ but with the phase "overwritten" with $\varphi_{\rm black}$:

[TABLE]

We repeat steps (2)–(4) until the energy of $\phi(x)$ converges.

The procedure described above wasn’t proven to give the exact result, although it may be rooted in the relations between an yrast state and a ground state solution for the interacting bosons placed in a one-dimensional hard-wall box potential Reichert et al. (2019), found by Gaudin Gaudin (1971). We further verify numerically if the final state $\phi(x)$ is indeed the solution of the NLSE  (1). We use $\phi(x)$ as an initial state $\phi(x,\,t=0)$ for the Eq. (1) and check whether $\left|\phi(x,t>0)\right|^{2}$ preserves its shape during evolution and whether the density dip moves with the expected speed. Example of such verification is presented in Fig. 6. Additionally we compare $\phi(x)$ with the series of the deepest solitonic solution we were able to find with the methods for gray solitons described in the previous section to check if they converge to the black soliton found with the method describe in this section.

Appendix B The Ansatz

In the main text we use the following Ansatz for a yrast state $|K\rangle$:

[TABLE]

where $\mathcal{N}$ is a normalization factor, and $\phi_{\rm MF}(x)$ is the solitonic solution of the NLSE (1) which has the average momentum equal to $\frac{2\pi\hbar}{L}K/N$,

[TABLE]

The normalization factor $\mathcal{N}$ from Eq. (28) is evaluated from the normalization condition:

[TABLE]

where the overlap $\int_{0}^{L}dx\,\phi_{\rm MF}^{*}(x)\,\phi_{\rm MF}(x-y+y^{\prime})$ and integrals over $y$ and $y^{\prime}$ are evaluated numerically.

B.1 The Ansatz as an eigenstate of the total momentum operator

Let us start with a comment that MF product state

[TABLE]

does not have well defined momentum, i.e. it is a wave packet of eigenstates with different momenta. In contrast the Ansatz  (28) is an eigenstate of the momentum operator. To prove that we begin by showing that translation of all particles by an arbitrary shift $\Delta x$ is equivalent to multiplication by a global phase factor:

[TABLE]

In the second last equality we have shifted the integration limits with no impact on its value due to the periodicity of integrated function.

Using the relation above, one can explicitly check that $\psi_{\rm Ansatz}$ is indeed an eigenstate of $\hat{P}$ with corresponding eigenvalue $\frac{2\pi\hbar}{L}K$:

[TABLE]

Note that the Ansatz is an eigenstate of the total momentum operator, irrespectively of the choice of the orbital $\phi(x)$. As discussed in the main text, only with an appropriate choice of orbital does the Ansatz become good approximation of the yrast state.

B.2 Comparisons between the mean-field product state and Ansatz in the limit $g\to 0$

As stated before, the MF product state (31) is not a state with well-defined momentum and therefore it cannot be a good approximation of the exact yrast state. That is why we have decided to consider a properly weighted, by a phase factor, superposition of MF solitons. To get some intuition on how these two are related, it is convenient to discuss their properties in the limit $g\to 0$. Here we discuss the case with mean total momentum $\frac{2\pi\hbar}{L}\frac{N}{2}$, as for this one the analytical formulas are the simplest. The MF orbital is:

[TABLE]

which we denote symbolically as $|\phi^{g\to 0}_{\rm MF}\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|2\pi\hbar/L\rangle)$. Corresponding $N$-particle product state is a superposition of the yrast states with coefficients given by square roots of the binomial distribution coefficients:

[TABLE]

where $|n_{0}:N-k,n_{2\pi\hbar/L}:k\rangle$ denote the Fock state with $N-k$ atoms in orbital with momentum [math] and $k$ atoms in obital with momentum $2\pi\hbar/L$. Similar analysis may also be done for any gray soliton Kaminishi et al. (2018). The expectation value of kinetic energy is, as expected, the same as for the exact yrast state $|n_{0}:N/2,n_{2\pi\hbar/L}:N/2\rangle$:

[TABLE]

However, for the mean value of interaction energy situation turns out to be slightly more complicated. We get:

[TABLE]

where in the first step we have used the fact than interaction energy operator $\hat{E}_{\rm int}$ preserve the total momentum of the system and therefore $\forall_{k\neq k^{\prime}}\langle n_{0}:N-k,n_{2\pi\hbar/L}:k|\hat{E}_{int}|n_{0}:N-k^{\prime},n_{2\pi\hbar/L}:k^{\prime}\rangle=0$. On the other hand, interaction energy of the yrast state $|n_{0}:N/2,n_{2\pi\hbar/L}:N/2\rangle$ is:

[TABLE]

We see, that the expectation value of the mean-filed product state’s energy is smaller than the energy of the yrast state. It must be so, as the formula for interaction energy for the Fock state $|n_{0}:N-k,n_{2\pi\hbar/L}:k\rangle$: $\frac{g}{2L}(N(N-1)+2Nk-2k^{2})$ takes its maximum in $K=N/2$, therefore increasing the impact of other Focks in MF product state may only decrease the energy.

We want to stress again, that this result does not lead to contradiction with the definition of the yrast state (i.e. the state with the lowest energy for given total momentum), as the MF product state is not an eigenstate of total momentum operator.

B.3 Conditional states, single-particle densities and function $g_{1}$ after measurement of particle positions

The measurement of the particle positions is expressed by an action of the field operator $\hat{\Psi}(x)$ on the Ansatz. To some extent it can be evaluated analytically. If the particle has been measured at random position $x_{1}$, then the conditional wave function $|\tilde{\psi}^{1}_{\rm Ansatz}\rangle$ is given by:

[TABLE]

where we used the fact that $\hat{\Psi}(x)|f\rangle^{\otimes N}=f(x)|f\rangle^{\otimes(N-1)}$, and introduced the proportionality symbol $\propto$ because the state $\hat{\Psi}(x_{1})|\psi_{\rm Ansatz}\rangle$ is not normalized.

By repetitive action of the field operator we can write down a conditional state after $m$ subsequent measurements which occurred at random positions $x_{1},\,x_{2}\,\ldots,\,x_{m}$:

[TABLE]

Given the conditional wave function, we write down its single particle density which is equal to ($m+1$)-th order correlation functions used in the main text:

[TABLE]

B.4 Interaction and kinetic energy

Here we discuss our procedure of computing the interaction energy:

[TABLE]

where the square of the norm is in fact a double integral:

[TABLE]

where the last term under integral is overlap between product states of orbitals $\phi_{\rm MF}$ occupied by ($N-2$) atoms, shifted by $(y-y^{\prime})$. In the limit $N\to\infty$ the term $\mathcal{N}^{2}\left(\int_{0}^{L}dz\,\phi^{*}_{\rm MF}(z-y^{\prime})\phi_{\rm MF}(z-y)\right)^{N-2}$ quickly decay at points $y\neq y^{\prime}$. If this term was approximated by a delta function, precisely $L\delta(y-y^{\prime})$, then the interaction energy (43) would coincide with the interaction energy of a product of MF solitons. Yet, we keep this term and see that for our finite size system it makes a difference.

Similarly we evaluate the kinetic energy:

[TABLE]

where we used indistinguishability of bosons.

B.5 Overlap between the state after $m$ measurements and yrast state

In Fig. 5 we present the projection of the state $|\tilde{\psi}^{m}_{\rm Ansatz}\rangle$ on the subspace of yrast states. This projection is evaluated as $\sum_{K^{\prime}=0}^{N-m}|a_{K^{\prime}}|^{2}$, where

[TABLE]

is the overlap between a state after $m$ measurement and an yrast state $|K^{\prime}\rangle$ with $N-m$ atoms. We compute the values of $a_{K^{\prime}}$ under an assumption that the Ansatz (28) is a fair approximation of the yrast state. Then $a_{K^{\prime}}$ reads:

[TABLE]

where we introduce another notation for solitonic solution of NLSE (1), $\phi_{\rm MF;K}$ to indicate the momentum $2\pi\hbar K/(NL)$. The parameters $x_{j}$ are the random positions, at which particle were measured, drawn from probability density function (41), and $\mathcal{N}_{m}$, $\mathcal{N}_{K^{\prime}}$ are the normalization factors of $|\tilde{\psi}^{m}_{\rm Ansatz}\rangle$ and $|K^{\prime}\rangle$, respectively.

B.6 Numerical methods

We discuss our numerical methods on an example of the interaction energy Eq. (43). The computation requires evaluation of overlap integrals $A[w]:=\int dz\,\phi^{*}_{\rm MF}(z)\,\phi_{\rm MF}(z-w)$ and $B[w]:=\int dz\,\left(\phi^{*}_{\rm MF}(z)\,\phi_{\rm MF}(z-w)\right)^{2}$ which we evaluate for discrete values of shift $w$ and store in a computer memory before the main computation. Then the interaction energy is approximated by:

[TABLE]

where integrals $\int dy$ were discretized to $\sum_{y}$ with discrete values of $y$ separated by $\Delta y$. The function $D[y-y^{\prime}]$ is equal to $A[y-y^{\prime}]$ for $y>y^{\prime}$ and it is equal to $\left(A[y^{\prime}-y]\right)^{*}$ otherwise. The results presented in Figs. 2 and 3 were evaluated on a numerical grid $1000$ points, i.e. with $\Delta y=0.001$. Already for the simplest integration scheme possible, based on the rectangle rule, we achieved converging results.

The results presented in other figures were not so sensitive to the numerical grid (as we were not interested in small differences of energies evaluated in different approaches), therefore we used numerical grid with $80$ points only, i.e. with $\Delta y=0.0125$.

Appendix C The Lieb-Liniger equations

Any eigenstate of LL Hamiltonian (2) is clearly defined by a set of real numbers $p_{j}$, so-called quasi-momenta, satisfying (see Eq. (2.15) in Lieb and Liniger (1963)):

[TABLE]

where total momentum and total energy of such state may be expressed respectively as $\sum_{j=1}^{N}p_{j}$ and $\sum_{j=1}^{N}\frac{p_{j}^{2}}{2m}$. After taking the logarithm of both sides of (49) and multiplying by $i$ we get:

[TABLE]

where $I_{j}$’s are called Bethe quantum numbers which uniquely characterize the state of a system. It is worth noting that the total momentum can be expressed as $\frac{2\pi\hbar}{L}\sum_{j=1}^{N}I_{j}$. The ground state corresponds to $\{I_{j}\}_{j=1}^{N}$ satisfying $I_{j+1}-I_{j}=1$, $I_{1}=-I_{N}$, i.e.

[TABLE]

Excitations of the type-I (Bogoliubov states) and the type-II (yrast states), as well as the first excited state with total momentum $\frac{2\pi\hbar}{L}K$ (for $K\geq 2$), may be generated by increasing the value of appropriate $I_{j}$ by $K$:

• •

Bogoliubov states: $I^{\prime}_{N}=I_{N}+K$

• •

yrast states: $I^{\prime}_{N+1-K}=I_{N+1-K}+K$

• •

the first excited state: $I^{\prime}_{N+2-K}=I_{N+2-K}+K$

One may check that in the limit $g\rightarrow 0$ the formula for the first excited states correctly reproduces the Fock state with $N-K+1$ particles in orbital with momentum [math], $K-2$ particles with momentum $2\pi\hbar/L$ and one with $2\cdot 2\pi\hbar/L$, i.e. $|n_{0}:N-K+1,n_{2\pi\hbar/L}:K-2,n_{2\cdot 2\pi\hbar/L}:1\rangle$.

C.1 Numerical evaluation

Choosing $L$, $\hbar^{2}$/$mL^{2}$ and $\hbar/L$ as the units of length, energy and momentum respectively, we get:

[TABLE]

For numerical convenience, instead of solving the above system of $N$ equations, we translate it to the problem of minimizing the function of $N$ variables:

[TABLE]

which may be straightforwardly solved numerically even for the number of particles $N$ on the order of $500$. Total energy of such state is given as:

[TABLE]