Quantum dark solitons in Bose gas confined in a hard wall box

TL;DR

This paper investigates quantum dark solitons in a one-dimensional Bose gas confined within a hard wall box, analyzing eigenstates, their relation to mean-field solitons, and the quantum nature of these excitations across interaction regimes.

Contribution

It extends the understanding of quantum dark solitons by analyzing Gaudin's eigenstates in a hard wall potential and exploring their properties in weak and strong interaction limits.

Findings

Dark solitons are present in the weak interaction regime.

Quantum nature of solitons becomes prominent under strong interactions.

Eigenstates are parametrized by their relation to single and multiple dark solitons.

Abstract

Schr\"odinger equation for Bose gas with repulsive contact interactions in one-dimensional space may be solved analytically with the help of the Bethe ansatz if we impose periodic boundary conditions. It was shown that in such a system there exist many-body eigenstates directly corresponding to dark soliton solutions of the mean-field equation. The system is still integrable if one switches from the periodic boundary conditions to an infinite square well potential. The corresponding eigenstates were constructed by M. Gaudin. We analyze weak interaction limit of Gaudin's solutions and identify parametrization of eigenstates strictly connected with single and multiple dark solitons. Numerical simulations of detection of particle's positions reveal dark solitons in the weak interaction regime and their quantum nature in the presence of strong interactions.

Equations

$$H=\underset{0}{\overset{L}{\int }}\mathrm{d}x\left[{\partial }_x{\hat{\psi }}^{†}{\partial }_x\hat{\psi }+c{\hat{\psi }}^{†}{\hat{\psi }}^{†}\hat{\psi }\hat{\psi }\right],$$

$$\gamma =\frac{c}{n},$$

$$\Psi (x_1=0,x_2,\dots ,x_N)=\Psi (x_1,x_2,\dots ,x_N=L)=0.$$

$$\times \prod _{i<j}\left(1-\frac{ic}{q_i+q_j}\right)\left(1+\frac{ic}{q_{\sigma (i)}-q_{\sigma (j)}}\right),$$

$$L{k}_i=\pi n_i+\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N\left[\arctan \frac{c}{k_i-k_j}+\arctan \frac{c}{k_i+k_j}\right],$$

$$E_{\{k\}}=\sum _{j=1}^N{k}_j^2.$$

$$\Psi (\{x\},\{k\})\stackrel{c\to 0^{+}}{\propto }\sum _{\sigma \in {\mathcal{S}}_N}\prod _{s=1}^n\sin \left(k_{\sigma (s)}{x}_s\right),$$

$${\Phi }_{sol}^{j-1}(\{x{\}}_N)\propto \prod _{s=1}^N\sin \left(\frac{\pi j{x}_s}{L}\right),$$

$$\mathcal{T}(j):n_1=n_2=\dots =n_N=j.$$

$$\psi(x)=\left\{\begin{array}[]{ll}\!\!-\sqrt{n}\tanh\left(\frac{x}{\xi}\right)&\textrm{for $0\leq x\ll L/2-\xi$}\\ \!\!\sqrt{n}\tanh\left(\frac{x-L/2}{\xi}\right)&\textrm{for $\xi\ll x\ll L-\xi$}\\ \!\!\sqrt{n}\tanh\left(\frac{L-x}{\xi}\right)&\textrm{for $L/2+\xi\ll x\leq L$.}\end{array}\!\!\!\right.$$

$$\psi (x)\propto \Psi (x_1,\dots ,x_{i-1},x,x_{i+1}\dots ,x_N).$$

$$\mathrm{F}\equiv \rho (L/2),$$

$$L\delta k=2\arctan \frac{c}{\delta k},$$

$$\frac{L}{2}(\delta k{)}^2\approx c⟹\underset{c\to 0^{+}}{\lim }\frac{c}{\delta k}=0.$$

$$\underset{c\to 0^{+}}{\lim }\Psi (\{x\},\{q\})=\sum _{\sigma \in {\mathcal{S}}_N,\{ϵ{\}}_N}{ϵ}_1\cdot \cdot \cdot {ϵ}_N\exp \left[i\sum _{s=1}^N{q}_{\sigma (s)}{x}_s\right].$$

$$\sum _{\sigma \in {\mathcal{S}}_N}\sum _{\{ϵ{\}}_N}{ϵ}_1\cdot \cdot \cdot {ϵ}_N\exp \left[i\sum _{s=1}^N{ϵ}_{\sigma (s)}{k}_{\sigma (s)}{x}_s\right]=\sum _{\sigma \in {\mathcal{S}}_N}\prod _{s=1}^N\left({\mathrm{e}}^{i{k}_{\sigma (s)}{x}_s}-{\mathrm{e}}^{-i{k}_{\sigma (s)}{x}_s}\right).$$

$$\sum _{\{ϵ{\}}_N}{ϵ}_1\cdot \cdot \cdot {ϵ}_N\exp \left[i\sum _{s=1}^N{ϵ}_{\sigma (s)}{k}_{\sigma (s)}{x}_s\right]=\prod _{s=1}^N\left({\mathrm{e}}^{i{k}_{\sigma (s)}{x}_s}-{\mathrm{e}}^{-i{k}_{\sigma (s)}{x}_s}\right).$$

$$\sum _{\{ϵ{\}}_N}{ϵ}_1\cdot \cdot \cdot {ϵ}_N\exp \left[i\sum _{s=1}^N{ϵ}_{\sigma (s)}{k}_{\sigma (s)}{x}_s\right]=\left(\sum _{{ϵ}_{\sigma (1)}}{ϵ}_{\sigma (1)}{\mathrm{e}}^{i{ϵ}_{\sigma (1)}{k}_{\sigma (1)}{x}_1}\right)\left(\sum _{{ϵ}_{\sigma (2)}}{ϵ}_{\sigma (2)}{\mathrm{e}}^{i{ϵ}_{\sigma (2)}{k}_{\sigma (2)}{x}_2}\right)\cdot \cdot \cdot \left(\sum _{{ϵ}_{\sigma (N)}}{ϵ}_{\sigma (N)}{\mathrm{e}}^{i{ϵ}_{\sigma (N)}{k}_{\sigma (N)}{x}_N}\right),$$

$$\sum _{{ϵ}_{\sigma (r)}}{ϵ}_{\sigma (r)}{\mathrm{e}}^{i{ϵ}_{\sigma (r)}{k}_{\sigma (r)}{x}_r}$$

$$=$$

$$\underset{c\to 0^{+}}{\lim }\Psi (\{x\},\{k\})={\left(2i\right)}^N\sum _{\sigma \in {\mathcal{S}}_N}\prod _{s=1}^n\sin \left(k_{\sigma (s)}{x}_s\right).$$

Full text

I Introduction

Self-reinforcing solitary solutions of non-linear wave equations maintaining their shape during time evolution are called solitons and appear in various fascinating phenomena KivsharOpticalSol . Solitons are particularly investigated in non-linear optics and ultra-cold atomic gases. In the latter case bosons may form a Bose-Einstein condensate (BEC) where all atoms occupy the same single-particle state and the many-body wave function factorizes into the product of identical single-particle states pethicksmith . In the presence of inter-particle interactions the product state describes atoms living in an averaged potential coming from the milieu of other identical particles (mean-field description). The behaviour of the single-particle state is determined by the Gross-Pitaevskii equation (GPE) which has analytical bright and dark soliton solutions in one-dimensional (1D) space for attractively and repulsively interacting system, respectively pethicksmith . GPE gives very accurate description of the solitons realized experimentally so far KivsharOpticalSol ; burger1999 ; denschlag2000 ; strecker2002 ; khaykovich2002 ; becker ; Stellmer2008 ; Weller2008 ; Theocharis2010 .

Theoretical investigations of quantum nature of dark solitons, i.e. properties which go beyond the standard Bogoliubov corrections dziarmaga2004 ; Mishmash2009_1 ; Mishmash2009_2 ; Dziarmaga2010 ; Mishmash2010 ; delande2014 ; kronke15 ; Hans2015 , employ the phase imprinting method burger1999 ; denschlag2000 ; becker ; Stellmer2008 ; Yefsah — starting with the system in the ground state one can carve a soliton by shining a short laser pulse on a half of the atomic cloud. Sufficiently long time evolution reveals quantum character of the dark soliton when the soliton position fluctuates on a length scale which can be much greater than the soliton width. Hence, the soliton location has to be considered as a quantum degree of freedom described by a probability density whose standard deviation increases with time dziarmaga2004 ; delande2014 . It is not easy to observe the many-body effects experimentally because it requires relatively small number of particles in the system in order to reduce atomic losses that are able to kill the quantum character of the solitons roberts2000 ; plodzien2012 . In the field of ultra-cold atomic gases the experimental techniques evolve very rapidly giving an opportunity to investigate systems for which the many-body effects play a key role Theocharis2010 ; dziarmaga2004 ; Lai89 ; Lai89a ; castinleshouches ; Weiss09 ; delande2013 ; corney97 ; corney01 ; martin2010b ; delande2014 ; kronke15 ; Mishmash2009_1 ; Mishmash2009_2 ; Hans2015 ; Boisse2017 .

The Schrödinger equation for some quantum many-body systems possesses analytical solutions. Generally, it happens for systems in lower dimensions when one can use a brilliant concept of the Bethe ansatz Bethe31 . It turns out that the Bethe ansatz approach is crucial in the analysis of bosonic and fermionic 1D systems interacting via point-like $\delta$-potentials Korepin93 ; Gaudin ; Oelkers2006 . Non-relativistic ultra-cold bosons in 1D space with contact interactions are described by renowned Lieb-Liniger model Lieb63 ; Lieb63a . One can expect dark soliton solutions in the presence of repulsive inter-particle interactions. In contrast to attractive interactions for which bright soliton solution corresponds to the lowest energy state, dark solitons have to be a reflection of excitations. Furthermore, if the system satisfies periodic boundary conditions all energy eigenstates are also eigenstates of the unitary operator that translates all particles by the same distance. Hence, they satisfy translation symmetry which is broken in the case of the mean-field soliton solutions. For those reasons the identification of dark soliton-like many-body eigenstates was particularly difficult. The conjecture and various evidences that the eigenstates belonging to the so-called type II excitation spectrum of the Lieb-Liniger model are strongly connected with the mean-field soliton solutions can be found in the literature kulish76 ; ishikawa80 ; komineas02 ; jackson02 ; kanamoto08 ; kanamoto10 ; karpiuk12 ; karpiuk15 ; sato12 ; sato12a ; sato16 ; Gawryluk2017 . Especially, it has been shown that dark soliton signatures are visible in the reduced single-particle density if the system is prepared in a proper superposition of the eigenstates coming from the second branch of the excitation spectrum sato12 ; sato12a ; sato16 . Ultimately, it was demonstrated that dark soliton signatures emerge in the course of measurement of particle positions if the system is prepared initially in a type II eigenstate Syrwid2015 ; Syrwid2016 .

In the presence of hard-wall boundary conditions the system does not possess the translation symmetry but the question if there are many-body eigenstates directly corresponding to mean-field dark soliton solutions in this case is still open. It is clear that in the limit of infinitely weak repulsion, one can identify many-body eigenstates that reveal density notches present in plots of the single-particle probability densities. Genuine solitons may appear if the inter-particle interactions are turned on and solutions of the GPE reveal phase flips and density notches whose widths are much smaller than the size of the hard-wall box. The goal of the present work is to show that there are many-body eigenstates that correspond directly to the solitonic structures. We will show that the phase flips and signatures of the density notches can be observed for arbitrary strength of the repulsive particle interaction. Therefore, in order to simplify the nomenclature in the present manuscript this class of states we will be sometimes dubbed dark solitons even in the nearly non-interacting and strongly interacting cases.

II Lieb-Liniger model

Ultra-cold system consisting of $N$ bosons with repulsive contact interactions in the 1D space may be described by the following Lieb-Liniger Hamiltonian Lieb63

[TABLE]

where $\hat{\psi}$ is the canonical Bose field operator and $L$ is the system size. The units have been chosen so that $2m=\hbar=1$ where $m$ is the particle mass. The dimensionless parameter

[TABLE]

reflects the strength of the interactions. The coupling constant $c>0$ and $n=\frac{N}{L}$ denotes the average density of particles. One deals with the weak interaction limit when $\gamma\ll 1$ and strongly interacting impenetrable bosons for $\gamma\gg 1$ Lieb63 ; Lieb63a .

II.1 Gaudin’s solution in the presence of hard walls

The Hamiltonian (1) with periodic boundary conditions has analytical solutions that can be found with the help of the Bethe ansatz Korepin93 . The Lieb-Liniger eigenstates are parametrized by a collection $\{k\}$ of $N$ real (if $c>0$) numbers called quasimomenta Gaudin . It is clear that the integrability of the model can be easily broken by switching on a trapping potential. However, the model is still integrable in the presence of hard-wall boundary conditions. The solutions were constructed by M. Gaudin Gaudin71 ; Gaudin ; Batchelor05 ; Tomchenko17 . Assuming that $0\leq x_{1}\leq x_{2}\leq\ldots\leq x_{N}\leq L$ we have to fulfil two conditions

[TABLE]

Firstly, using McGuire’s optical analogy McGuire , one constructs solutions in the semi-infinite 1D space $x_{i}\geq 0$ vanishing at $x_{1}=0$. The solutions turn out to be superpositions of elementary waves and take the following form

$\displaystyle{\Psi(\{x\},\{q\})=\sum_{\sigma\in\mathcal{S}_{N}}\sum_{\{\epsilon\}}\epsilon_{1}\epsilon_{2}\cdot\cdot\cdot\epsilon_{N}\exp\left[i\sum_{s=1}^{N}q_{\sigma(s)}x_{s}\right]}$

[TABLE]

where $\epsilon_{i}=\pm 1$, $q_{i}=\epsilon_{i}|k_{i}|$ (with $0<|k_{1}|<\ldots<|k_{N}|$). The first sum has to be taken over all possible $N$-element combinations of $\pm 1$. The second sum refers to all permutations $\sigma$ from the permutation group $\mathcal{S}_{N}$. Therefore, calculation of a single value of (4) requires summation over $2^{N}N!$ elements Gaudin71 ; Gaudin ; Batchelor05 ; Tomchenko17 .

Secondly, imposing the vanishing of the wave function at $x_{N}=L$ one obtains elegant system of coupled equations

[TABLE]

called Gaudin’s equations where integer numbers $n_{i}$ are parameters that define an eigenstate. In the case of repulsive inter-particle interactions the system of Gaudin’s equations (5) has unique real solutions $k_{i}$ Tomchenko17 . In order to deal with admissible physically different solutions one considers only $k_{i=1,\ldots,N}>0$ and the set of integers $\{n\}$ where $1\leq n_{1}\leq n_{2}\leq\ldots\leq n_{N}$ Gaudin71 ; Gaudin ; Batchelor05 ; Tomchenko17 . By construction the set $\{k\}$ determines the energy of the eigenstate (4) Gaudin ; Gaudin71 ; Batchelor05

[TABLE]

We would like to stress that even if two or more $n_{i}$ parameters are equal, the resulting solutions $k_{i}$ of (5) are always distinct.

II.2 Collective soliton-like excitations

For the non-interacting system ($c\rightarrow 0^{+}$) the solutions (4) reduce to superposition of products of sine functions (see Appendix)

[TABLE]

with $k_{j}\stackrel{{\scriptstyle c\rightarrow 0^{+}}}{{\longrightarrow}}\pi n_{j}/L$ for $n_{j}=1,2,3,\ldots$ (negative values are physically equivalent to non-negative ones). Obviously, they coincide with the well known solutions of the problem of non-interacting particles in an infinite square well potential. If we choose identical integer parameters $n_{i}=j$ in the Gaudin’s equation (5), all elements of the set $\{k\}$ approach $\pi j/L$ in the limit $c\rightarrow 0^{+}$ but they are always slightly different if $c$ is not strictly equal to 0.

In the non-interacting case of Bose gas confined in the 1D box of length $L$ the solutions resembling solitons correspond to product states where all particles occupy the same excited eigenstate of the single-particle problem,

[TABLE]

where $j=2,3,\ldots$. It is clear that the wave function (8) reveals $j-1$ density notches and phase flips and thus resembles signatures of $j-1$ dark solitons. When $j=1$ it reproduces the ground state. In the $c\rightarrow 0^{+}$ limit the eigenstates (4) reduce to (8) for the following collections of the integer numbers in (5)

[TABLE]

We believe that in the presence of particle interactions there is a range of the system parameters where the eigenstates (4), with quasimomenta $\{k\}$ parametrized by (9), not only resemble but reveal dark solitons unambiguously.

Looking at the expression (9) we notice that $\mathcal{T}(1+s)$ corresponds directly to $s$-fold collective excitation of the aforementioned non-interacting system. For convenience, we will use this terminology also in the presence of inter-particle repulsion. The excitation given by $\mathcal{T}(j)$ will be dubbed as the first collective excitation for $j=2$ and the higher ($j-1$)th collective excitation for $j>2$.

II.3 Numerical method

Properties of the many-body eigenstates $\Psi(\{x\},\{q\})$ are hidden in the structure of terms in (4) whose number dramatically grows with increasing number of particles $N$. Hence, despite the fact that one has the analytical solutions it is not easy to understand their features. In order to study the eigenstates of the many-body system we will simulate measurements of particle positions.

In order to simulate results of the measurements one has to choose randomly positions of $N$ particles according to the $N$-dimensional probability density. In practice it is often done sequentially: one by one particle at each step calculating conditional probability density for a choice of a next particle javanainen96 ; dziarmaga03 ; dziarmaga06 ; Syrwid2015 ; Dagnino09 ; Kasevich20016 . This procedure is very difficult in the present case. Instead, we follow an equivalent idea of usage of the Monte Carlo algorithm of Metropolis et al. Metropolis1953 . As in Refs. Gajda_PauliCrystal ; Syrwid2016 we perform Markovian walk in the configuration space. That is, by sampling $N$-dimensional probability distribution $\left|\Psi(x_{1},\ldots,x_{N},\{q\})\right|^{2}$ we generate collections of sets $\mathcal{X}=$ $\{x_{1},\dots,x_{N}\}$ of particles’ positions. The procedure is based on step by step acceptance of sets $\mathcal{X}^{\prime}$ with probability $p={\rm min}\left(1,|\Psi(\mathcal{X}^{\prime})|^{2}/|\Psi(\mathcal{X})|^{2}\right)$, where $\mathcal{X}$ is the previously accepted set of particles’ positions. Although the method is very efficient, only a few body systems are numerically attainable. In the studies of the Lieb-Liniger model with periodic boundary conditions we were able to investigate 8-particle system Syrwid2015 ; Syrwid2016 . There, number of terms in the Bethe ansatz solutions increases like $N!$ with an increase of the total number of particles $N$. In the present case we need to use the Gaudin’s solution (4) where the number of terms proliferates $2^{N}$ times faster than in the previously used Bethe ansatz solutions. Hence, in order to collect reliable statistics, we perform simulations for the system consisting of $N=6$ and 7 bosons only.

III The analysis of collectively excited many-body eigenstates

The ground and collectively excited states of a Bose system are usually described with the help of the mean-field approximation if we deal with the weak interaction regime ($\gamma\ll 1$). The assumption that all particles occupy the same single-particle state leads to GPE pethicksmith . Analytical solutions of GPE in the presence of the hard wall box potential exist and are given by Jacobi elliptic functions Carr_HWSoliton . If $L\gg\xi$, where $\xi=1/\sqrt{cn}$ is the so-called healing length pethicksmith , the analytical solution related to a dark soliton can be approximated by the following function dziarmaga2004

[TABLE]

The healing length $\xi$ describes a typical distance over which the condensate wave function forgets about, e.g., boundaries of the system. Therefore, in a finite system with hard walls positioned at the boundaries we expect disturbance of the particle density close to the boundaries on the same length scale as the width of the dark soliton notch.

When quantum many-body effects are taken into account, one may expect that the position of the dark soliton starts fluctuating in different realizations of the particle positions measurements. The fluctuations are the more significant the stronger interactions are present Syrwid2015 ; Syrwid2016 . Hence, we expect that the averaged particle density $\rho(x)$, i.e. averaged over many realizations of the measurement process, should reveal the shallower dark soliton notches the stronger repulsion is.

III.1 First collective excitation

The signatures of single dark soliton are expected to be observed in the weak interaction limit when one prepares the system in the eigenstate (4) choosing the parameters in (5) so that $\forall_{i=1,\ldots N}:n_{i}=2$. We have performed numerical simulations of the particles positions measurement for 6- and 7-body systems confined in a box of length $L=1$. Single realization of the detection process produces a sequence of $N$ positions only. Therefore, in order to investigate averaged particle density, we have repeated the measurement procedure starting with the same many-body eigenstate many times. Collecting all results of the simulations in a histogram allows us to look at the average particle density.

In panels (a) and (c) of Fig. 1 we present histograms of the obtained particle positions in many realizations of the detection process for $N=6$ and 7 and for a wide range of the interaction strength ($\gamma=0.01,1,100$). Smearing of the density notch with increasing $\gamma$ is clearly visible at the center of the box. It is caused by the fact that the soliton position varies randomly from realization to realization and the range of the fluctuations is the larger the stronger interactions are Syrwid2015 ; Syrwid2016 . In a single realization of the measurement process one obtains a set of particle positions $\{x_{1},\dots,x_{N}\}$. Let us choose from this set arbitrary $N-1$ positions and consider the eigenstate (4), parametrized by the set of $n_{i}$ that satisfy (9) with $j=2$, as a single-particle wave function of the last remaining particle $x_{i}=x$,

[TABLE]

It turns out that the wave function $\psi(x)$ reveals density notch and a phase flip at the center of the notch. We identify the position of the phase flip with the position of the dark soliton. The distributions of the phase flip positions in the case of different interaction strengths are depicted in Fig. 1 for $N=6$ (b) and $N=7$ (d). Additionally, in Fig. 1(e), we present the quantity

[TABLE]

which measures the degree of filling of the notch observed in the averaged particle densities $\rho(x)$ versus $\gamma$ parameter in the case of $N=7$. The degree of filling $\mathrm{F}$ saturates already for $\gamma\approx 5$. For strong interactions ($\gamma=100$) one can also observe $N+1$ oscillations in the histograms, cf. panels (a) and (c) of Fig. 1. In the ground state case and for a small particle number one can expect $N$ local maxima in the average particle density that correspond to average locations of $N$ bosons with strong repulsive interactions — the density should be identical to the density of non-interacting fermions when $\gamma\rightarrow\infty$ (Tonks-Girardeau limit) girardeau1960 ; Paredes2004 . In Fig. 1 we see one local maximum more which seems to be related to the fact that if there is a phase flip between neighbouring particles, then their relative distance is modified and on average it results in an additional oscillation in the density profile. We expect that the oscillations can be visible for small $N$ only and for large particle number they will become negligible, i.e. the averaged particle density profile will be almost flat (except the edges of the system).

The distributions of positions of the phase flip in the case of strong repulsion ($\gamma=100$) also reveal oscillations [dotted histograms Fig. 1(b) and (d)]. We note that the positions of local maxima of such distributions roughly coincide with positions of local minima of corresponding averaged particle densities (dotted histograms in Fig. 1(a) and (b), respectively). Moreover, dealing with even (odd) number of particles we observe that in the presence of strong inter-particle interactions the distribution of the phase flip positions has the minimum (maximum) in the center of the box $x=L/2$. The presence of the oscillations and the correlations between positions of the maxima in the averaged particle densities and positions of the minima in the distributions of phase flip positions results from the fact that for large $\gamma$ at a space point where there is greater probability to observe a particle, there must be smaller probability to find the phase flip.

The significant quantum many-body effects appear for $\gamma>0.1$ but they do not destroy the signatures of solitons like the density notch and phase flip in single realizations of the detection process. Especially the phase flip can be clearly observed even in the very strong interaction regime ($\gamma\gg 1$). In Fig. 2 we show probability density and phase of the wave function (11) for the last 7th particle (provided 6 particles have been already detected) in the 7-particle system obtained in a single realization of the measurement for $\gamma=1$ and $100$. One easily notices that every measurement leaves its mark on the wave function. In the Fig. 2(a) we observe 6 slight incisions in the profile of the probability density for $\gamma=1$. The stronger inter-particle repulsion is, the deeper incisions are observed. In the case of strong repulsion ($\gamma=100$) we note that the probability density is essentially nonzero only in regions away from the measured positions of particles. The reason is the strong repulsion does not allow for detection of two particles close to each other. Therefore, it is very difficult to establish where is the phase flip of the wave function (11) by looking at the profile of the probability density only if $\gamma\gg 1$ — compare Fig. 2(a) and Fig. 2(b).

In the non-interacting case the many-body eigenstate that we analyse here reduces to a simple product state (8) with $j=2$. Then, the parity symmetry, which is fulfilled by the Hamiltonian (1), becomes apparent because such a state is also an eigenstate of the operator that transforms each $x_{i}$ in the following way $(x_{i}-L/2)\rightarrow-(x_{i}-L/2)$. The resulting average particle density vanishes at the centre of the box, cf. Fig. 1. Average particle densities related to $\gamma=1$ and $\gamma=100$ do not possess this property but it is not in contradiction with the parity symmetry. Indeed, the interactions between particles introduce coupling between states where all particles occupy anti-symmetric modes with states where some even number of particles occupy symmetric modes — superposition of such states is also an eigenstate of the parity operator.

III.2 Higher collective excitations

Increasing the number $j$ in the eigenstate parametrization (9) one expects to observe increasing number of notches that appear in averaged particle densities for weak and intermediate interaction strength. We performed numerical simulations of particles detection for $N=7$ in the system prepared initially in the eigenstates parametrized by $n_{i}$ numbers that fulfil (9) with $j=3$ and $j=4$. The results confirm that there are $j-1$ notches in the average probability densities for $\gamma=0.01$ and 1, see Fig. 3. Analysis of single measurement realizations reveals $j-1$ phase flips in the wave function (11) for the last particle in the system. The phase flips are key signatures of dark solitons for the intermediate interactions and indicates where the solitons are localized in single realizations. The phase flips are also present in the strong interaction regime.

Similarly to the case of the first collective excitation, considered in the previous subsection, in the presence of intermediate and strong inter-particle repulsion the density notches are blurred because of fluctuations of positions of the phase flips. We also observe oscillations in the averaged particle densities for $\gamma=100$, see Fig. 3. One notices that the number of the oscillations is equal to $N+j-1$, i.e. the number of particles plus the number of phase flips. In the strong interaction case and for small $N$ the average particle density related to the ground state would reveal $N$ local maxima. The presence of a phase flip between a pair of particles modifies its relative distance. It turns out that the presence of $j-1$ phase flips between different neighbouring particles results in $j-1$ additional local maxima in the average particle density as compared to the ground state case.

In the case of the eigenstate parametrized by $n_{i}$ numbers that fulfill (9) with $j=3$, the distribution of distances between positions of two phase flips is obviously concentrated at $L/3\approx 0.33$ in the non-interacting case, cf. Fig. 4. For weak and intermediate interactions it is still localized around $L/3$ but its width increases with $\gamma$. This can be explained by the observation that the wave function (11) must drops at the edges on a length scale similar to half of the width of the density notches. If we now imagine that we merge the edges of the box, the shape of the modulus squared of the wave function (11) resembles not 2 but 3 solitons on a ring. Thus the expected mean separation between solitons should be approximately equal to $L/3$ but due to fluctuations of the soliton positions the width of the distribution increases with $\gamma$.

For $\gamma\rightarrow\infty$ we approach the Tonks-Girardeau regime where impenetrable bosons tend to localize at $L/N$ distances one from each other. Then, the phase flips are expected to be located half the way between two neighboring particles. Thus, the distributions of distances between two phase flips may reveal peaks at integer multiple of $L/N$. In Fig. 4 we show the distribution for $\gamma=100$ where the peaks structure emerges. Similar structures were also observed in the case of periodic boundary conditions Syrwid2016 .

IV Conclusions

We have considered gas of bosons interacting via point-like $\delta$-potential confined in a one-dimensional hard wall box. Eigenstates of the system were analytically constructed by M. Gaudin and they were parametrized by a set of positive integers $\{n\}$. We show that Gaudin’s solutions, in the limit of infinitely weak interactions, factorize to simple eigenstates of non-interacting particles in the box. In this case each number $n_{i}$ corresponds simply to a number of an excited eigenstate of a single-particle in the square well potential. If all numbers $n_{i}$ are equal to the same integer $j>1$, the average particle density for infinitely weak interactions resemble $j-1$ dark soliton-like notches.

Genuine dark solitons may appear only for non-vanishing inter-particle interactions. In this case the numbers $\{n\}$ do not have clear interpretation and should be treated just as parameters defining a many-body eigenstate uniquely. Nevertheless, the eigenstates parametrized by $\forall_{i}$ $n_{i}=j$ do correspond to $j-1$ dark solitons. In order to show it we have performed numerical simulations of measurements of particle positions. It turns out that the wave function before the measurement of the last particle always possesses dark soliton signatures. That is, there are $j-1$ phase flips of the wave function and $j-1$ probability density notches for the weak and intermediate interactions. The small numbers of particles ($N=6,7$) we consider do not allow us to fulfil the conditions, $\gamma\ll 1$ and $\xi\ll L$. Hence, shape of the density notches do not reproduce a hyperbolic tangent function like in Eq. (10). For strong interactions, the positions of the phase flips fluctuate strongly from one realization of the measurement process to another one that results in smearing of the notches in the average particle density.

We have also investigated the relative distance between neighboring phase flips in the case of $j=3$. The results show that depending on interaction strength there are specific distances the phase flips localize more eagerly.

Acknowledgments

The authors express their sincere gratitude to D. Delande for fruitful discussions. Support of the National Science Centre, Poland via projects: No.2016/21/B/ST2/01095 (A.S.) and No.2015/19/B/ST2/01028 (K.S.) is acknowledged. A.S. acknowledges support in the form of a special scholarship from the Marian Smoluchowski Scientific Consortium "Matter-Energy-Future", from the KNOW funding. This work was performed with the support of EU via Horizon2020 FET project QUIC (nr. 641122)

Appendix

Let us consider the Gaudin’s equations for 2-particle system with identical numbers $n_{1}=n_{2}$ in (5). Subtracting the particular equations one obtains

[TABLE]

where $\delta k=k_{2}-k_{1}$. Taking the tangent of both sides and using the series expansion, $\tan(x)=x+{\cal O}(x^{3})$, one notices that in the limit $c\rightarrow 0^{+}$

[TABLE]

Last result holds true also in the case of $N$-particle system parametrized by $s\leq N$ identical numbers $n_{j}$. It obviously coincides with the statement that physically relevant solutions may always be ordered such that $0<k_{1}<k_{2}<\ldots<k_{N}$ Tomchenko17 .

Using the fact that $\delta k$ vanishes slower than $c$ in the limit $c\rightarrow 0^{+}$ the equation (4) reduces to

[TABLE]

Now we want to show the following identity

[TABLE]

It is clear that it is enough to show the identity for a single arbitrary permutation $\sigma$,

[TABLE]

In general, we can change positions of epsilons (or equivalently permute indices) in the following way $\epsilon_{1}\epsilon_{2}\cdot\!\cdot\!\cdot\epsilon_{N}=\epsilon_{3}\epsilon_{1}\cdot\!\cdot\!\cdot\epsilon_{N}\cdot\!\cdot\!\cdot\epsilon_{5}=\epsilon_{\sigma(1)}\epsilon_{\sigma(2)}\cdot\!\cdot\!\cdot\epsilon_{\sigma(N)}$. Now, the equality (17) can be shown using the fact that

[TABLE]

which ends the proof because every single term of the product,

[TABLE]

Hence, in the considered limit, the equation (4) can be rewritten in terms of sine functions

[TABLE]