Dark solitons revealed in Lieb-Liniger eigenstates

We study how dark solitons, i.e., solutions of one-dimensional, single-particle, nonlinear, time-dependent Schrödinger 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 has been 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 uniﬁcation 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-ﬁeld solitons emerge in a single-particle density through repeated measurements of particle positions in the ansatz state. The postmeasurement state turns out to be a wave packet of yrast states of the reduced system. DOI: 10.1103/PhysRevResearch.2.033368


I. INTRODUCTION
The famous Lieb-Liniger (LL) model [1,2] describes particles moving along a circle and interacting via δ interatomic potential.Such a simple interaction turns out to be a wellsuited 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 [3][4][5][6][7].
The same system, of N atoms with contact interaction, is often treated within a simple mean-field (MF) approximation, based on the nonlinear Schrödinger equation (NLSE): where the wave function φ 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 nonlinear optics [8] to hydrodynamics [9,10].It is also a rare example of a model with physical applications supporting solitonic solutions [11], observed in atomic gases [12], plasma [13], water waves [14], and ferromagnetic materials [15].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 [16].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 φ MF (x) undergoes a sharp π jump, while gray solitons have a nonzero density dip with the phase jump strictly smaller than π .
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 Lieb [1], 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 [17,18].The type-II excitations are simply the manybody eigenstates that minimize the energy for a fixed total momentum, sometimes called yrast states [19].Together with the type-I excitations (corresponding to Bogoliubov quasiparticles [1]), 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?First, 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 nonlinearity in the FIG. 1.As shown in Ref. [1], among eigenstates of the Lieb-Liniger model there are special states, forming two branches of the many-body elementary excitations (left panel).In this paper, we show that states forming the type-II spectrum, i.e., yrast states, can be approximated by a superposition of the mean-field dark solitons with relative phases (marked with color) depending on the solitons' positions (right panel).We discuss the validity of a suitable ansatz and use it to unify different views on correspondence between the yrast states of the Lieb-Liniger model and the mean-field solitons from the nonlinear Schrödinger equation.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 natures of the MF solitons and the yrast states are so different, the question arises whether they have something in common other than 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.[20,21] 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.[22][23][24][25] 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 Refs.[24,[26][27][28][29].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 Sec.II, we remind readers of the LL model and yrast state and define the parameter regime we are interested in.Our ansatz for yrast states is introduced in Sec.III.Validity of this ansatz is discussed in Sec.IV.In Secs.V and VI, we show how the different constructions of the MF solitons out of the yrast states presented in Ref. [20] and Refs.[22][23][24] prove to be different views of the same object.To make this paper self-complete, we describe in detail all relevant analytical previous results and our numerical approaches in the Appendixes.

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 [30]: where x j denotes position of the jth particle.As this system is translationally invariant, the values of total momentum P = −i h N j=1 ∂ 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 have been known since 1963 [1].Among the eigenstates, there are special ones that are called elementary excitations, which can be divided into two families.Before Ref. [1], the approximated theories were 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 is given by [31].The unexpected second family, revealed by the exact solution presented in Ref. [1], consists of the aforementioned yrast states.These are also called "one-hole excitations" [1] and the lowest energy solutions for fixed total momentum [19].The yrast state with the total momentum of P = 2π h L K is represented by |K , 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 desirable to see the effects that are substantially different from the ideal gas case.Therefore, we require the healing length ξ := h2 L/gmN, which is close to the soliton width, be much shorter than the system size L [32].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 0, but the MF parameter ng := h2 mξ 2 , with n := N/L denoting a gas density, is fixed.In this limit, the LL coupling constant γ ∝ g/N [1] 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 semianalytical results [22,24,25], the majority of approaches are devoted to small systems of ≈10-20 atoms [20,21,[33][34][35][36] solved with brute force methods or around ≈100 atoms solved with sophisticated and time-consuming numerics [28,29,37,38].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 φ(x) represented by N j=1 φ(x j ).Next, the ansatz has to belong to the same momentum space as the yrast state |K ; i.e., the translation of the N-particle wave function by x needs to be equivalent to multiplying it by a factor e i 2π L K x : For any orbital φ, one may construct the states satisfying above condition by taking a continuous superposition of product states shifted by the translation operator e −i Py/ h and multiplied by the phase factor e i 2π L Ky over all possible shifts y (see Appendix B1 for formal justifications): where N is a (real) normalization factor.
As a yrast state |K is the lowest energy state for fixed total momentum equal exactly to 2π hK L , and therefore the task would be to find an orbital φ 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 B4).On the other hand, it is known that energies of the yrast states |K and MF solitons are related, at least in the weak interaction regime [17].Therefore, as the ansatz for yrast state we choose the state (4) with φ(x) = φ MF (x), where φ MF (x − vt ) is the solitonic solution of the NLSE with the average single-particle momentum −i h ∂ x equal to 2π h K/(NL): We also use state defined in Eq. ( 5) in the Dirac notation: L Ky e −i Py/ h|φ MF ⊗N .( The exact form of solitonic solution on the circle φ 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π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 an ansatz was presented in the context of Bose-Einstein condensation [39].Condensate, as defined via Penrose-Onsager criterion [40], is supposed to appear in the system of bosons at very low temperature.As show in Ref. [39], surprisingly, even at T = 0, when the FIG. 2. Expectation value of Lieb-Liniger Hamiltonian (10) evaluated as a function of the mean-field parameter ng in the meanfield approximation, solution of Eq. (1) (solid gray line), for the yrast state (Appendix C) (empty yellow circles), for the ansatz (6) (green crosses), and using perturbation theory (dashed burgundy line).In the top panels N = 100, while in the bottom ones N = 500.The parameter K is set to N/4 in the panels to the left and to N/2 in the panels to the right, as indicated in the upper left corner of each graph.The top panels share a common energy scale and so do the bottom ones.
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 [21,23].
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 Appendixes B4, B6, 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π h/L.Its average energy is a linear function of interaction strength g equal to The second reference curve is the average energy of the MF state, with N atoms occupying a single orbital φ MF (solid gray line).
As expected, the energies 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→0 |ψ ansatz ≡ lim g→0 |K for any K = 0 [32].In the same limit, the energy of the dark MF soliton is actually slightly smaller than the energy of the yrast state (see Appendix B2).
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 analyze 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 | ψ ansatz |K | 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 quasimomenta 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, it can be decomposed into the basis of many-body eigenstates of the same total momentum: where |K j is the jth excited eigenstate with total momentum 2π h L K and energy E j exc , and ∞ j=0 |α j | 2 = 1.The excited eigenstates are listed in the ascending order with respect to their energy, i.e., E 1 exc The average energy ψ ansatz | Ĥ|ψ ansatz can be expressed with the help of Eq. ( 7) as where we use an inequality E j exc > E 1 exc and E yrast is the energy of the yrast state with the total momentum 2π h L K. Therefore, the fidelity between the ansatz (6) and the corresponding yrast state, i.e., |α 0 | 2 , obeys the inequality 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 FIG. 3. A rough lower bound for the fidelity between the ansatz (6) and yrast state, F B (9), as a function of the mean-field parameter ng obtained for N = 20 (triangles), 100 (crosses), and 500 (circles).The parameter K is set to N/4 in the panel to the left and to N/2 in the panel to the right, as indicated in the upper right corner of each graph.The panels share a common vertical scale.implicate uselessness of the ansatz for stronger interactionsthe fidelity may be still close to one but calculating it directly would require much more advanced numerical methods and would be unfeasible for a large number of atoms.
Given that for ng 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 was done in Ref. [20].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 where ˆ (x) [ ˆ † (x)] is an annihilation (creation) field operator of a boson at position x satisfying the commu- The second quantization formalism is also very handy in performing any computations within the ansatz (see Appendix B).The object of interest in Ref. [20] is a mth-order correlation function 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 mth particle to be measured.In Ref. [20], function ( 11) is considered for increasing m after subsequent measurements of particles.It was observed that ρ m resembles density |φ MF | 2 of the MF soliton from the NLSE (1).However, in Ref. [20] the agreement between MF solitons and the mth-order correlation function of a yrast state was demonstrated only for the healing length ξ > L and for a few particles, namely for almost a noninteracting system.In contrast, 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 replaced by our ansatz |ψ ansatz (6).We begin with a computation of the singleparticle reduced density matrix ρ 1 which is used as a PDF to draw a random position x 1 .Subsequently, we compute second-order correlation function ρ 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 ρ m (x) (11), evaluated for the ansatz (6) (for details of our calculations, see Appendix B3).
In the left panel of Fig. 4, we show samples of correlation functions ρ m (x) of different orders m calculated for the parameters ng = 25, N = 500, and K = 50, 125, 250.Densities |φ MF (x)| 2 of the solitonic solution of Eq. ( 1) are shifted so that their notch positions overlap with that of ρ 150 (x) for direct comparison.We observed that the notch position of ρ 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, such as that exemplified by ρ 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 ξ = 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 arg{φ MF (x)} 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 evaluated for the ansatz (6) (for details of our calculations, see Appendix B3).As shown in the right panel of Fig. 4, the phase of g 1 converges to arg {φ MF (x)} for increasing m.The phase of MF soliton is shifted by the same amount as the corresponding density |φ MF (x)| 2 in the left panel.(11).Right: Phase arg{g 1 } of the correlation function (12).Both are plotted with respect to the position of mth particle x.Different solid lines correspond to different orders m from the shallowest to the deepest dip in the left panels m = 1 (violet), 25 (blue), 50 (green), 100 (yellow), and 150 (red).The thick gray line corresponds to the mean-field (MF) solitonic solution of Eq. (1).Parameters: ng = 25, N = 500, and K = 50, 125, and 250 for the top, middle and bottom rows, respectively, as indicated in the bottom-right corner of each graph.
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 ξ significantly smaller than L.
On the other hand, in Refs.[22][23][24], the MF solitons were constructed from the many-body eigenstates of the LL Hamiltonian (2) in a completely different way, as explained 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.[22][23][24][25] where the dark solitons are constructed as a specific superposition of yrast states.Namely, the MF product state is expressed as where a K are expansion coefficients (drawn from a chosen distribution) and |K N denotes a yrast state of the system with N particles and the total momentum 2π h L K .A comprehensive discussion of different a K choices can be found in Refs.[22][23][24][25].
An interesting question arises as to 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 ρ 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 | ψm of the system in a state |ψ after measuring positions of m particles given by 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 ρ i defined in Eq. (11).
The average density in the conditional state ( 14) ρ(x) := ˆ † (x) ˆ (x) is equal to the (m + 1)-th-order correlation function ρ 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 |ψ 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 2π h L K .For the ansatz (6), one can easily find the conditional state [41] Because of the factors φ 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 | ψm Ansatz 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 where |ψ j N−m is an eigenstate of the system, which is not the yrast state, and The question is whether the conditional state | ψm Ansatz 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 K |a K | 2 and the a K distribution [42].
In the top panel of Fig. 5, we present sum of weights K |a K | 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 correlates with the weight of the yrast subspace for a given value of m being closer to 1.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 as functions of total momentum 2π h L K 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 [22][23][24][25].In this section, we have shown that the dark solitons hosted in the nonideal 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 nonlinear 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 [20,[22][23][24][25] 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 mthorder correlation function resembling the mean-field soliton, as discussed in Ref. [20].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 analyzed in Refs.[22][23][24][25].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 [21,22,24,25].It was argued that the time of blurring should increase to infinity in the thermodynamic limit [21,23].However, this hypothesis has to be verified.This can be done again employing our ansatz.

ACKNOWLEDGMENTS
We acknowledge fruitful discussions with K. Rz ążewski, M. Gajda, and K. Sacha.In particular, we thank M. Gajda for convincing 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.

APPENDIX A: MEAN-FIELD SOLITONS 1. Mean-field gray solitons
The solitonic solution φ MF (x) of NLSE (1) was discussed several times in the literature [11,23,43].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 Ref. [23].Solitonic solution of NLSE is a running wave, φ MF (x, t ) = φ MF (x − vt ), with speed v.In what follows, we will omit the time dependence and give separately a solitonic density, and its phase, and its speed v.The approach presented in Ref. [23] was devoted to the case of gray solitons, i.e., when ϕ(x) is continuous and ρ(x) is always larger than 0. 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 where sn(u, k) is the Jacobi elliptic function, is the incomplete elliptic integral of the third kind, with the modulus of Jacobi's elliptic function k, and 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 |φ MF | 2 = 1 lead to the following relations between parameters a 1 , a 2 , a 3 , and the elliptic modulus k: ) ) where K (k) and E (k) are the elliptic integrals of the first and the second kinds, respectively.
To compute phase and density of a soliton that has a desired minimum of the density a 1 , we first solve numerically Eq. (A4) for the elliptic modulus k, and then we use Eqs.(A5) and (A6) 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 φ MF (x) at any position x, using Eqs.(A3) and relation φ MF (x) = √ ρ(x)e iϕ(x) .The average momentum p := h i dx φ * MF (x)∂ x φ * MF (x) is then computed numerically.To find the MF soliton with target momentum p target = 2π hK/(NL), we repeat the steps described above varying a 1 , until the numerically determined momentum matches p target .After a few bisection steps with respect to a 1 , the target MF  1).The mean-field black soliton φ(x, t = 0) with density dip at x = 0, used as an initial state, was found numerically according to steps (A8)-(A11).The black line is the reference line-a trajectory of a point moving with the speed v black = hπ/mL.soliton is found-we stop bisection when relative discrepancy between numerically computed momentum and the target momentum is less than 10 −6 .

Mean-field black solitons
The formulas presented in the previous section are derived under the assumption that the density ρ(x) given in Eq. ( A3) is always greater than 0. Therefore, they cannot be used in the case of a black soliton.In fact, already for a gray but very deep solitons, Eq. (A4) becomes very demanding, as discussed in Ref. [23] (see, for instance, Table 4 in Ref. [23] 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 method from Karpiuk [44] and used it successfully in Ref. [45].Precisely, we look for solutions with the phase where the signum function sgn is equal to 1 for positive arguments and 0 for negative ones.The signum function introduces a discontinuity in the phase, a π jump, which is the characteristic feature of a black soliton.
The minimal energy state in the space of functions with phase ϕ black (x) is found with the split-step imaginary time evolution, implemented as follows: (1) We start with an arbitrary function φ(x).
(3) We normalize the output of the previous step: We define a new function φ(x) as φ(x) but with the phase "overwritten" with ϕ black : We repeat steps ( 2)-( 4) until the energy of φ(x) converges.The procedure described above was not 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 [46], found by Gaudin [47].We further verify numerically if the final state φ(x) is indeed the solution of the NLSE (1).We use φ(x) as an initial state φ(x, t = 0) for the Eq. ( 1) and check whether |φ(x, t > 0)| 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 φ(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 , where N is a normalization factor, and φ MF (x) is the solitonic solution of the NLSE (1), which has the average momentum equal to 2π h L K/N, The normalization factor N from Eq. (B1) is evaluated from the normalization condition: does not have well-defined momentum; i.e., it is a wave packet of eigenstates with different momenta.In contrast, the ansatz (B1) is an eigenstate of the momentum operator.To prove that we begin by showing that translation of all particles by an arbitrary shift x is equivalent to multiplication by a global phase factor: In the second to 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 ψ ansatz is indeed an eigenstate of P with corresponding eigenvalue 2π h L K: Pψ ansatz (x 1 , . . ., x N ) = Pψ ansatz (x 1 + x, . . ., Note that the ansatz is an eigenstate of the total momentum operator, irrespective of the choice of the orbital φ(x).As discussed in the main text, only with an appropriate choice of orbital does the ansatz become good approximation of the yrast state.

Comparisons between the mean-field product state and ansatz in the limit g → 0
As stated before, the MF product state (B4) 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 insight into how these two are related, it is convenient to discuss their properties in the limit g → 0.Here we discuss the case with mean total momentum 2π h L N 2 , as for this one the analytical formulas are the simplest.The MF orbital is which we denote symbolically as |φ g→0 ).The corresponding N-particle product state is a superposition of the yrast states with coefficients given by square roots of the binomial distribution coefficients: where |n 0 : N − k, n 2π h/L : k denotes the Fock state with N − k atoms in orbital with momentum 0 and k atoms in orbital with momentum 2π h/L.Similar analysis may also be done for any gray soliton [24].The expectation value of kinetic energy is, as expected, the same as for the exact yrast state |n 0 : N/2, n 2π h/L : N/2 : However, for the mean value of interaction energy situation turns out to be slightly more complicated.We get where in the first step we have used the fact than interaction energy operator Êint preserves the total momentum of the system and therefore ∀ k =k n 0 : N − k, n 2π h/L : k| Êint |n 0 : N − k , n 2π h/L : k = 0. On the other hand, interaction energy of the yrast state We see that the expectation value of the mean-field 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 : ] takes its maximum in K = N/2, and 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.

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 ˆ (x) on the ansatz.To some extent, it can be evaluated analytically.where we used the fact that ˆ 1) , and introduced the proportionality symbol ∝ because the state ˆ (x 1 )|ψ Ansatz 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 . . ., x m : 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: . (B15)

Interaction and kinetic energy
Here we discuss our procedure of computing the interaction energy: where the square of the norm is in fact a double integral: where the last term under integral is an overlap between product states of orbitals φ MF occupied by (N − 2) atoms, shifted by (y − y ).In the limit N → ∞, the term quickly decay at points y = y .If this term was approximated by a δ function, precisely Lδ(y − y ), then the interaction energy (B16) 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 where we used indistinguishability of bosons.

Overlap between the state after m measurements and yrast state
In Fig. 5, we present the projection of the state | ψm ansatz on the subspace of yrast states.This projection is evaluated as , (B20) where we introduce another notation for solitonic solution of NLSE (1), φ MF;K , to indicate the momentum 2π hK/(NL).The parameters x j are the random positions, at which particle were measured, drawn from probability density function (B14), and N m , N K are the normalization factors of | ψm ansatz and |K , respectively.

Numerical methods
We discuss our numerical methods on an example of the interaction energy Eq.(B16). 2 and 3 were evaluated on a numerical grid 1000 points, i.e., with 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 as sensitive to the numerical grid (as we were not interested in small differences of energies evaluated in different approaches); therefore, we used a numerical grid with 80 points only, i.e., with 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 quasimomenta, satisfying (see Eq. 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 Excitations of the type-I (Bogoliubov states) and the type-II (yrast states), as well as the first excited state with total momentum 2π h L K (for K 2), may be generated by increasing the value of appropriate I j by K: (1) Bogoliubov states: One may check that in the limit g → 0 the formula for the first excited states correctly reproduces the Fock state with N − K + 1 particles in orbital with momentum 0, K − 2 particles with momentum 2π h/L, and one with 2 × 2π h/L, i.e., |n 0 : N − K + 1, n 2π h/L : K − 2, n 2•2π h/L : 1 .

Numerical evaluation
Choosing L, h2 /mL 2 , and h/L as the units of length, energy, and momentum respectively, we get which may be straightforwardly solved numerically even for the number of particles N on the order of 500.Total energy of such a state is given as

FIG. 4 .
FIG. 4. Left:The mth-order correlation functions ρ m(11).Right: Phase arg{g 1 } of the correlation function(12).Both are plotted with respect to the position of mth particle x.Different solid lines correspond to different orders m from the shallowest to the deepest dip in the left panels m = 1 (violet), 25 (blue), 50 (green), 100 (yellow), and 150 (red).The thick gray line corresponds to the mean-field (MF) solitonic solution of Eq. (1).Parameters: ng = 25, N = 500, and K = 50, 125, and 250 for the top, middle and bottom rows, respectively, as indicated in the bottom-right corner of each graph.

FIG. 5 .
FIG. 5. Top: Sum of weights of the yrast subspace K |a K | 2 as a function of the number of atoms lost m.Each value for a given m corresponds to a different stochastic sequence of particle positions measured {x i } obtained for N = 100 (triangles), 250 (circles), 500 (rhombus), and 1000 (crosses), with K set to N/10 in the left panel and N/4 in the right one.Bottom: Five representative distributions of weights a K as functions of the total momentum for m = 10 (left) and m = 40 (right).Each symbol corresponds to a different set {x i }.Parameters N and K are set to 1000 and N/10, respectively.The mean-field parameter ng is equal to 25 for every graph.

FIG. 6 .
FIG.6.Evolution of the density |φ(x, t )| 2 governed by nonlinear Schrödinger equation (1).The mean-field black soliton φ(x, t = 0) with density dip at x = 0, used as an initial state, was found numerically according to steps (A8)-(A11).The black line is the reference line-a trajectory of a point moving with the speed v black = hπ/mL.

1 = ψ ansatz |ψ ansatz = N 2 L0 1 .
2π L K (y−y ) [ φ MF | ⊗N e −i P(y−y )/ h|φ MF ⊗N ] MF (x) φ MF (x − y + y ) dx φ * MF (x) φ MF (x − y + y ) and integrals over y and y are evaluated numerically.The Ansatz as an eigenstate of the total momentum operator Let us start with a comment that MF product state N j=1 φ MF (x j ) ( B 4 ) If the particle has been measured at random position x 1 , then the conditional wave function | ψ1 ansatz is given by ψ1 ansatz ∝ ˆ (x 1 )|ψ ansatz = N L 0 dy e i 2π L Ky ˆ (x 1 )e −i Py/ h|φ MF ⊗N = N L 0 dy e i 2π L Ky ˆ (x 1 ) |φ MF (x − y) ⊗N = N L 0 dy e i 2π L Ky √ Nφ MF (x 1 − y)|φ MF (x − y) ⊗(N−1) , (B12) wherea K := K | ψm ansatz (B19)is the overlap between a state after m measurement and an yrast state |K with N − m atoms.We compute the values of a K under an assumption that the ansatz (B1) is a fair approximation of the yrast state.Then a K readsa K = K ψm ansatz = N m N (z − y )φ MF;K (z − y) (N−m) The computation requires evaluation of overlap integrals A[w]:= dz φ * MF (z) φ MF (z − w) and B[w] := dz [φ * MF (z) φ MF (z − w)]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 byE int ≈ g N (N − 1) N 2 2 y y B[y − y ] (D[y − y ]) N−2 e i 2π L K (y−y ) ( y) 2 , (B21) where integrals dy were discretized to y with discrete values of y separated by y.The function D[y − y ] is equal to A[y − y ] for y > y and it is equal to (A[y − y]) * otherwise.The results presented in Figs.

2 j
(2.15)  in Ref.[2])(−) N−1 e −i L h p j = exp −2i N l=1 arctan h(p j − p l ) mg , ∀ j =l p j = p l , (C1)where total momentum and total energy of such state may be expressed respectively as N j=1 p j and N j=1 p 2m .After taking the logarithm of both sides of (C1) and multiplying by i, we getL h p j = 2π I j − 2 N l=1 arctan h(p j − p l ) mg ; ∀ l = j I j = I l , . The ground state corresponds to {I j } N j=1 satisfying I j+1 − I j = 1, I 1 = −I N , i.e., {I j } N j=1 = − ) yrast states:I N+1−K = I N+1−K + K, (3) the first excited state:I N+2−K = I N+2−K + K.

2 N l=1 arctan p j − p l g 2 ,
numerical convenience, instead of solving the above system of N equations, we translate it to the problem of minimizing the function of N variables: min p 1 ,...,p N N j=1 p j − 2π I j + (C4)