Breakdown of the mean field for dark solitons of dipolar bosons in a one-dimensional harmonic trap

We directly compare the mean-field and the many-body approach in a one-dimensional Bose system in a harmonic trap. Both contact and dipolar interactions are considered. We propose a multi-atom version of the phase imprinting method to generate dark solitons in the system. We begin with a general analysis of system dynamics and observe the emergence of a dark soliton and a shock wave. Center of mass and soliton motion become decoupled because the shock wave oscillates with the trap frequency and soliton does not. A detailed investigation of frequencies reveals significant differences between results obtained in the mean-field and the many-body pictures.


I. INTRODUCTION
Solitons -solutions of non-linear integrable differential equations which propagate without dispersion -appear across many areas of science that range from physics to biology and medicine [1]. Among a number of known equations supporting solitonic solutions an important example constitutes of the non-linear Schrodinger equation called also the Gross-Pitaevski equation (GPE). Although mainly used in the context of weakly interacting ultracold bosons [2], it also describes the properties of the electric field of light in the nonlinear media [3].
If atoms in a Bose-Einstein condensate repel each other by point-like interactions, a resulting solitonic solution of the GPE takes on form of a dark soliton -a density dip with a phase jump across its density minimum. An analytical expression describing dark solitons in a uniform gas characterized by the GPE were already found in the 70s [4,5]. Shortly after achieving the first condensation in 1995, dark solitons were successfully generated by the phase imprinting method in the experimental setups with ultracold 87 Rb [6] and 23 Na [7] atoms trapped in a harmonic confinement.
One of the most classic results about dark solitons in a BEC interacting only by short-range forces concerns its dynamics in a harmonic trap. In the so-called Thomas-Fermi regime for ultracold bosons, a characteristic frequency of soliton's oscillations is expressed by the trapping frequency ω and equals ω/ √ 2 [8] that was also observed in the experiment [9]. However, this remarkably robust results for shortrange interactions, does not hold for the dipole-dipole atomic interactions [10]. In this case, the oscillation frequency of solitonic structures in a repulsive dipolar BEC -predicted only in 2015 [11,12] -highly depends on the strength of the atomic interactions and the interplay between local and non-local contributions to the total energy. With the recent progress on quantum gases consisting of atoms with considerable magnetic moments like 164 Dy [13,14] and Er [15], dipolar systems are now within experimental reach and call for deeper analysis.
The non-linear mean-field (MF) theory of ultracold bosons provides an approximate description of interacting cold atoms. The underlying many-body (MB) model is linear and the state of the system is given by the many-body wave-function depending on positions of all particles. The discussion of correspondence between dark solitons present in the MF and many-body solutions of the full Hamiltonian has a long history. The best known example refers to the link between dark solitons moving on the circumference of a ring and type-II excitations from the seminal Lieb-Liniger model [16,17] in the context of contact interacting particles, see [18][19][20][21][22][23][24][25][26][27][28] and references therein. However, little is known about many-body states possesing features of dark solitons in a one-dimensional harmonic trap both for contact and purely dipolar interactions. Here, we aim to partially bridge this gap by studying many-body dark solitons in weakly interacting trapped systems of only few atoms far beyond the Thomas-Fermi regime. In particular, we investigate the oscillations of the many-body solitons and compare them to the dark solitons described by the GPE. Our results not only establish a link between the MF approach and the full many-body theory, but also can be verified in modern experiments with a precise control over only a few atoms in optical lattices or single traps (see for instance [29][30][31][32][33]).
The work is organized as follows. In Sec. II, we introduce our many-body model for atoms interacting repulsively by contact or purely dipolar forces in a quasi-1D harmonic trap. We also remind the corresponding GPE. Then, in Sec. III, we analyse the many-body eigenstates of the system. Surprisingly, there is no good candidate for a many-body solitonic state among them, in a stark contrast to the ring geometry. Therefore, we introduce the many-body phase imprinting method of creating dark solitons in a harmonic trap. In Sec. IV, we finally present our results. We investigate dynamics of the many-body solitons. We calculate the oscillations frequencies for different coupling strengths for both types of interactions and compare our findings with the GPE. The most important result is a significant disagreement between frequencies obtained in the MB and the MF approaches.

II. THE MODEL
We investigate a system of N repulsive bosons trapped in a harmonic potential We assume that the transverse confinement is tight, so the wave function stays in the lowest energy level in the Y-and Z-directions meaning our system is quasi-1D. The aim of our study is to compare contact and dipolar solitons in many-body and mean-field approaches. The system of N repulsive bosons is often approximated by the using the quasi-1D GPE where U (r) is the trapping potential and V (r) is an interaction potential. Throughout the paper we use oscillatory units with ω, mω , √ mω, 1 ω as units of energy, length, momentum and time respectively. In the context of ultracold atoms, this equation is also called the MF description of weakly interacting bosons. This approach provided correct predictions of many properties of the Bose-Einstein condensate, including its shape, energy, normal modes of excitations and many other nonlinear phenomena. Moreover, the GPE supports dark solitons in a form of solutions with a density notch and a quickly changing phase, which have also been experimentally produced with the phase imprinting method. However, the MF model assures a simplified description of system of N repulsive bosons based on a naive assumption that every atom is in the same state. It is only an approximation of the more fundamental many-body approach. As the nonlinear GPE supports solitons, it is important to look for solitons in the linear many-body approach and compare them. It has been done in the Lieb-Liniger model [16,17] but, to the best of our knowledge, this is the first paper where multi-atomic solitons are considered for the harmonically trapped system.
In order to describe a system of N bosons following the many-body approach, one needs to derive the wave function depending on positions of all particles. One of possible ways of deriving the many-body wave function is to diagonalize a Hamiltonian matrix. The Hamiltonian of the system under investigation can be written in a form: with T i being the single-particle kinetic energy operator and V ij the two-body interaction operator. The Hamiltonian (3) can be rewritten as a sum of the Hamiltonian of the noninteracting quantum harmonic oscillator and the interaction term. Therefore, in the second-quantization the Hamiltonian (3) reads: withΨ (x) being a bosonic field operator and H osc = We study systems where V(r) is either a short range V sr (r) or a long range dipolar V dd (r) potential. The short range potential is V sr (r) = g sr δ(r) with parameter g sr = V sr (r)dr defining strength of the interaction. We study repulsive systems with g > 0. In order to obtain the explicit formula describing V dd , we follow the procedure described in paper [11]. We introduce an aspect ratio of the trap σ = ω ⊥ ω and a dipolar coupling strength A dd yielding This effective quasi-1D potential comes from the integration of the full 3D dipolar interaction over both transverse variables. The area of our interest is a long range part of interaction, so we assume that the contact term of the effective dipolar potential is exactly cancelled possibly with the help of Feshbach resonances. As we want to investigate similarities and differences between systems interacting via contact and dipolar forces, we define dipolar strength parameter g dd From now on, we will keep σ = 0.1 and compare systems described by the same strength parameters g dd = g sr = g.
It is important to mention that in the case of the MF approach only a gas parameter (N −1)g defines the system. However, in the multi-atom approach both N and g are separately relevant. From the numerical point of view, the most optimal basis to describe the system are quantum harmonic oscillator eigenstates. We use the second quantization formalism and define a basis of Fock states. We take into account all states with a given number of particles, which energies are smaller than a cut-off energy. Defining the cut-off in the energy space, rather than in the momentum space is more efficient method of obtaining numerical convergence [34].
In order to obtain eigenstates and energies of the system, we diagonalize the Hamiltonian matrix where |i , |j are states belonging to the Fock space.

III. SOLUTIONS
Having access to both eigenenergies and eigenstates we are ready to look for solitons in the system. Following recent papers investigating many-body solitons in the Lieb-Liniger model [24,35], one could think that also in this case dark solitons could be identified among eigenstates of the Hamiltonian (4). Studying the repulsive case, we are interested in properties of dark solitons. In this situation density forms a single notch and in the area of the notch phase exhibits a jump of π. Keeping this in mind, one can ask a question if any single particle state fulfil these conditions.
We start our analysis with the ideal gas. In this case, the excited state of quantum harmonic oscillator (QHO), N |0, 1, 0, ... = |0, N, 0, ... , seems like a reasonable candidate because it has both the density dip and the phase jump, so one can try to find the eigenstate of the interacting system with the highest contribution of the aforementioned state among all the Fock states.
This turnes out to be non-trivial even for weak interactions. We expected the maximum |0, N, 0, ... occupation tend to one as the interactions become weaker. Instead, it was approaching different values depending on a number of particles considered. This is a peculiar property of the harmonic trap caused by evenly spaced energies of H osc . Once the interactions become weaker, the energy of the eigenstate with the highest |0, N, 0, ... contribution approaches 3N 2 but there are multiple other states with the same energy leading to degeneracy. For example in the case of N=2, |0, 2, 0 has the same energy as |1, 0, 1 namely 3. Hence in contrast with the Lieb-Liniger model, these eigenstates remain the combination of several other states with the energy equal to 3N 2 even for vanishing interactions.
Even if the excited eigenstate would form the dark soliton, still it would be hard to realize this state in the experiment. Therefore we decided to follow a different approach and try to replicate the experimental procedure of phase imprinting. This method creates a dark soliton in a BEC via a pulse of a far detuned laser applied on one half of the condensate and so creates a phase difference between the left and the right side. The length of this pulse is tuned to create a phase difference of π and hence causes emergence of the dark soliton. There are not many papers discussing the phase-imprinting method in the many-body approach [36]. To the best of our knowledge it was only applied together with density engineering, which is not the case in the experimental realisation. As in real life, our implementation of phase imprinting modifies only the phase of the wave function and is equivalent to multiplying the ground state wave function by an arbitrary phase factor Ψ(x 1 , x 2 , ..., x N ) = Φ(x 1 , x 2 , ..., x N )e iφ(x1,x2,...xN ) , (9) where Ψ is the many-body wave function of a solitonic state, Φ is the ground state wave function and φ is an arbitrary phase factor. For the numerical convenience we choose where α is the parameter changing sharpness of phase jump which can also be controlled in experiments. The optimal sharpness of the phase jump has to be tuned to fit healing length of the soliton. It means that the stronger the interaction the narrower the phase jump has to be. As for now, the solitonic wave function Ψ is merely an initial condition. We derive the time evolution of the system by expressing Ψ in the basis of eigenstates of the system as follows: where ψ i (x 1 , ..., x N ) and E i are the eigenstates and the eigenvalues of the Hamiltonian (4). In order to visualise a soliton, we derive a one-particle density ρ(x 1 , t) = dx 2 ...dx N Ψ * (x 1 , ..., x N , t)Ψ(x 1 , x 2 , ..., x N , t). (12) As we aim to obtain the MF dark solitons from Eq. (2) and compare them with MB solutions, we employ an analogous scheme. Firstly, we find a ground state ψ GPE (x) for given parameters by using the well-known imaginary time evolution (ITE) technique. At this point, we can compare ground states obtained in MF and MB approaches. Both density profiles and ground-state energies (up to 2 % difference for the highest g) are in a very good agreement for both dipolar and contact interactions and for all coupling strenghts considered in this work. Then, we imprint the same phase as in the MB calculations, namely Ψ GPE (x) = ψ GPE (x)e iφMF (x) , with φ MF = tan −1 (αx). Finally, we evolve Eq. (2) in a standard real-time evolution with Ψ GPE (x) as an initial condition.
Note that we can calculate the quantum depletion for any many-body state, in particular for a many-body ground state before and after phase imprinting, by diagonalazing a single particle density matrix constructed from the many-body wave function. It provides a tool for comparing mean-field and many-body results.

IV. RESULTS
Having a model of the experimental method of phase imprinting and being able to calculate the evolution of dark solitons, we can focus on properties of contact and dipolar manybody solitons and compare them with the MF results. Firstly, we would like to focus our attention on general aspects of the evolution of many-body solitons in the harmonic trap. In order to study the evolution of the system we plot the one-particle density as a function of time and space in Fig. 1 for N = 6 dipolar bosons and g = 0.3. It reveals that the phase imprinting method causes not only the dark soliton to emerge but also a shock wave in the form of a density peak initially moving in the opposite direction to the soliton. Plots of density profiles in consecutive time steps shown in Fig. 2 reveal more details of soliton evolution. The local density minimum moves from the center of the trap towards the left side as long as the density of the notch is greater than zero. When the soliton becomes black, namely when the one-particle density in the dip reaches zero, its velocity also equals zero, both indicating a turning point. The soliton begins to move right and becomes shallow in the center of the trap. The relation between the depth of the dark soliton and its velocity is one of the fundamental properties of solitons and have been studied in a number of papers [4,5,37].
As we pointed before the phase imprinting method creates a soliton but also a shock wave. This effect has been already observed in experiments implementing phase imprinting [6,7]. Both the shock wave and the soliton oscillate in the trap harmonically but the shock wave oscillates with the     trap frequency. It is then worth to study the frequency of the soliton movement as it was one of the factors differentiating contact and dipolar solitons in the mean-field approach [10]. One of properties of contact solitons revealed by number of studies [8,37,38] is that the frequency of oscillation does not depend on the strength of interactions and equals ω TF = 1 √ 2 ω. However, this result is obtained in the Thomas-Fermi (TF) limit assuming that the background density varies slowly on the scale of the soliton. The kinetic energy of particles in the TF limit can be neglected compared to the potential energy. However, satisfying this condition in the case of small systems would demand very strong interactions causing atoms to deplete the ground state and thus making our and the mean-field result incomparable. Our studies focused on the case of small systems in the many-body approach far from the TF limit and with the quantum depletion of the ground state before phase imprinting not exceeding 5%.
In order to analyse the frequency of oscillation, we trace the position of a local minimum of the one-particle density (in the many-body approach) or condensate wave function (in the mean-field picture) in consecutive time steps. We trace the minimum until it crosses the center of the trap which gives us half of the period.
Recent paper investigating dipolar solitons in the meanfield approach revealed significant differences between contact and dipolar solitons in the TF regime [10], one of them exhibited by the frequency of oscillations. In contrast to previously described mean-field contact solitons, the frequency of dipolar solitons depends on the interaction strength and, in general, the dipolar soliton frequency is smaller than the one obtained for contact solitons. It is then worth asking if manybody solitons exhibit similar behaviour already for weak interactions.
To answer this question we directly compare many-body and mean-field solitons for contact and dipolar interactions. We plot frequencies obtained for system of N=6 particles for the increasing coupling strength g in Fig. 3. Firstly, we note that the frequencies of MF and MB solitons differ significantly. In the MF picture, contact and dipolar solitons are almost indistinguishable. In opposition, MB dipolar solitons oscillate much slower than contact ones which agrees with the recent MF analysis within the TF approximation in [10].
It is then important to ask why the MF model far away from the TF regime is inconsistent when applied to the studied system. We can indicate two factors that may play a significant part. On the technical level, it seems like rapidly changing excited QHO states contribute to the difference in the energy between short-range and dipolar interactions in the many-body picture. On the other hand, the MF wave-function does not vary at the scale given by the range of dipolar interaction. Hence, for systems far from the TF regime, the dipolar interacting scenario almost does not differ from the short-range interaction case. The other factor contributing to the difference between MB and MF solutions is the phase imprinting method. Before imprinting, the systems are comparable as the depletion of the MB ground state does not exceed 5%. After the procedure, it rises by 10-20% depending on the sharpness of the phase jump. It means that the excited fraction cannot be neglected anymore in the MB calculations while it is not present at all in the MF case. This is a very important observation as the depletion of the ground state is fundamentally bound with phase-imprinting method and need to be taken into consideration when studying small systems both theoretically and experimentally.
Having discussed the differences between MF and MB results, we focus on the properties of MB contact and dipolar solitons. The system proves to be interaction-sensitive as the frequency varies significantly with both the strength and the range of the interaction. In both cases, the frequency decreases with the increasing coupling strength g. The dipolar solitons always oscillate slower than their short-range counterparts. As we are far from the TF limit, the frequency of contact solitons does not converge to 1 FIG. 3: (color online) Comparison of frequency as a function of a coupling strength g for contact and dipolar interactions in the many-body (circles) and mean-field (squares) approaches for N = 6 particles. Mean-field solitons are almost indistinguishable with dipolar ones being only slightly slower. On the other hand, many-body solitons differ significantly. Dipolar solitons always oscillate slower as for the MF dark solitons in the TF regime studied in [10].
cillation but also their lifetimes. While it is hard to define a sharp condition for the soliton to be indistinguishable from the background, we noted that contact solitons live significantly longer than dipolar counterparts, with the lifetime strongly dependent on the coupling strength g.

V. CONCLUSIONS
The goal of this paper is to compare dark solitons in the mean-field and many-body approaches for contact and dipolar interactions. We have began our many-body analysis with calculating eigenstates and energies of the system via the numerical diagonalization of the Hamiltonian matrix. We have found that one cannot identify the dark solitons with a specific eigenstate of the system, in stark contrast to the well-known situation of atoms in a ring trap. This follows directly from quantum degeneracy of multi-particle eigenstates of the noninteracting gas in the harmonic trap because the single particle eigenenergies are spaced evenly.
In order to study many-body solitons in the harmonic trap we have introduced the multi-atom version of the phase imprinting method. Just as in the classic experiment it causes not only the soliton but also the shock wave to appear. Those waves oscillate with different frequencies and thus movements of the soliton and the center of mass are decoupled. We investigate the frequency of oscillation of dark solitons to reveal similarities and differences between contact and dipolar solitons and compare our many-body analyses with mean-field results.
Although calculated for small and weakly interacting systems, our studies comparing many-body contact and dipolar solitons uncover similar features as previously discussed mean-field results within the Thomas-Fermi regime [10]. The frequency of oscillations for dipolar solitons strongly depends on the coupling strength and is lower compared to contact solitons for the corresponding interaction strength. For comparison, we also analyzed contact and dipolar solitons in our small system induced by the phase imprinting method at the MF level. The MF approach fails in the case of our system as the dipolar and contact solitons are almost identical and their properties differ significantly from the MB solitons. We can define two factors that cause significant differences between our MB and MF solitons. Firstly, the quickly oscillating excited QHO states play an important role in the case of dipolar interaction. Secondly, the phase imprinting method enlarges a depletion of the ground state and the excited fraction is no longer negligible.