Dipolar dark solitons

We numerically generate, and then study the basic properties of dark soliton-like excitations in a dipolar gas confined in a quasi one dimensional trap. These excitations, although very similar to dark solitons in a gas with contact interaction, interact with each other and can form bound states. During collisions these dipolar solitons emit phonons, loosing energy, but accelerating. Even after thousands of subsequent collisions they survive as gray solitons and finally reach dynamical equilibrium with background quasiparticles. Finally, in the frame of classical field approximation, we verified, that these solitons appear spontaneously in thermal samples, analogously to the type II excitations in a gas of atoms with contact interaction.


Introduction
The main motivation for our work comes from recent progress concerning the role of dark solitons in the physics of a gas with contact interaction only. First, let us recall that solitons are localized waves propagating through non-linear media without any spreading and also moving through each other without changing the shape. Among equations of broad interest for physicists, there are a few which exhibit solitons. One of them is the non-linear one-dimensional Schrödinger equation: in which one can find bright (the case with attracting forces, g 0 < ) and dark (the repulsive case g 0 > ) solitons.
Equation (1) applies, among others, to the dynamics of bosons interacting via a δ-like potential, i.e. the interaction potential between two bosons V z z -¢ ( )is equal to g z z . d -¢ ( ) In this context, equation (1) is only an approximation valid in the variously defined weak interaction limit [1]. Such a gas of bosons can be described by a more fundamental many-body model, the Lieb-Liniger model, which has been used rigorously to find the ground state [2] as well as the elementary excitations [3]. It turns out that there are two branches of elementary excitations, staying on the same footing. The first branch has been identified as the Bogoliubov excitations. The second one, consisting of the so-called type II excitations, in the research performed by the next generations of physicists turned out to be the solitonic branch [4][5][6][7][8][9][10]. It is now clear that the dark solitons, as the elementary excitations, appear spontaneously in a gas at finite temperature [11].
Equation (1) does not capture many important and common situations, even in the weakly interacting limit. For example, it neglects possible external traps in which the gas is confined, and it assumes contact interaction only. In particular, it misses the gas with a dipolar interaction, which has been studied recently in more and more laboratories [12][13][14][15]. Among many motivations to study a dipolar gas is definitely its unusual excitation spectrum, which in certain situations, i.e. close to instabilities, can have a local minimum for non-zero quasi momenta, the so-called roton spectrum [16]. It seems natural to ask if here is a place for the second branch and from what excitations it can be built. To our knowledge, there is no study of the elementary excitation spectrum beyond the Bogoliubov approximation, neither by solving the underlying many-body problem, nor on the mean field level.
Moreover, although the literature about the bright solitons is very rich (see parts of the review [17]), relatively little is known about dark solitons. For the time being the shapes of dipolar solitons in a gas with both contact and dipolar interaction has been presented [18,19]. Surprisingly, the life time of the dark solitons can be increased with some tricks even in the higher dimensional space [20]. The dark solitons in non-local nonlinear media have been the subject of research in optics. It has been found that two optical dark solitons interact with each other via non-trivial potential, and that they can even form a bound state [21]. In this context, although the effective forces between photons can be non-local, it is always assumed that they are not long-range.
Within this work, limited to the mean field description, we look closer at the basic properties of the dark solitons in the dipolar media. The aim of this paper is to study dark dipolar solitons and contrast them with the true solitons existing in the gas with contact interaction only. The analysis is performed within the mean field description, but we include also the finite temperature case. We start with the detailed description of our model in section 2. After presentation of the basic features of single soliton (shapes, dynamics) in section 3, we discuss a motion of two dipolar solitons in section 4. Setion. 5 is devoted to the analysis of the thermal samples.

The Model
We assume that the dipolar gas is confined in the potential with a tight harmonic trap in the X-and the Y-directions. The space in the Z-direction is assumed to be finite, with the length L and with periodic boundary conditions. These conditions impose a quantization of the momenta k z in this direction, k z is a multiple of .  are the units of length, time, energy and temperature, respectively. We assume that the dipole moments of the constituent atoms are polarized along the X-axis, namely when the dipolar potential in the momentum space has the form 2 where g dd 3D is a dipolar coupling strength. Thus, the atoms move on the circumference of a circle, having the dipole moments perpendicular to the circle-plane. In this configuration dipoles mostly repel each other, which favors the dark solitons. The geometry of our problem is sketched in figure 1.
We will investigate the system within the mean field model. The time-dependent mean field wave-function obeys the three-dimensional Gross-Pitaevskii equation (gpe) 3 We assume that the transverse confinement is so tight, that the wave function stays in the lowest energy level in the X-and Y-directions, namely when the following single-mode approximation holds: x y z t x y z t , , ; e , , 4 where 0 f is the Gaussian wave function: ) is the transverse oscillator width. The necessary condition underlying this approximation is that both, the thermal energy k T B and the interaction energy, are much below the energy of the first excited state in the transverse direction.
Assuming such a situation we multiply the three-dimensional GPE by x y , 0 0 f f ( ) ( ) integrate the result with respect to x and y and finally we reach the one-dimensional effective GPE The Fourier transform of the effective one-dimensional dipolar potential V k eff ( ) reads [20,22]. = p^, and l L l º^is equal to the aspect ratio of the trap. The function f which appears in equation (8) is equal to f a a a e Ei , a = ( ) ( ) where Ei is the exponential integral. The interaction potential, written in the spatial representation, has two parts: the smooth one which starts at finite value for zero distance and the one proportional to the Dirac delta function, precisely g z z . con d -¢ ( ) In this paper we aim to study the effects due entirely to the long-distance nature of the dipolar interaction, hence from now on, we will assume that there is no delta-like contribution to the potential, so g 0. con = Hence we assume that the van der Waals part of the interaction (possibly with the help of Feshbach resonance) exactly cancels the dipolar contact contribution, namely g g dd which implies g 0. con = The spatial representation of the potential V eff is visualized in figure 2. It is worth stressing that equation (7) describes the movement of slices of a dipolar gas as marked in figures 1 and 2. If such slices are far from each other, then the details of the distribution of dipoles within such slices plays no role, and we can treat them as single large dipoles. In such situations the interaction between them scales as z 1 . 3 The situation changes when the distance between them is of the order of the transverse width l ⊥ . This gives the characteristic range of the effective 1D potential (8). Finally in the limit of slices going to the same place, the effective dipolar potential converges to a constant.

Single solitons
In the gas with contact interaction only, the phase of a solitonic wave function is changing by π around the region with minimal density. This motivates us to look for the solitons in a dipolar gas within the wave function which minimizes the energy but under constraint of a jump of the phase equal to π. If we only force the jump of the phase equal to π, then the periodic boundary condition would not be satisfied. Hence, except the π jump required by the soliton, we add to the constraint on the phase a linearly changing function of the position: such that the phases at the edges of the box are equal. Due to this additional term the gas is moving with the speed L. p When showing the result for the dynamics, we would unrotate this motion, i.e. at each time t we present the density shifted back by Lt p and accounting for the periodic boundary condition. To find the state with the lowest energy but with the phase equal to the constraint (9) we use the method of the imaginary time evolution, forcing at each step the phase (9) 4 . Precisely, the evolution from the time t to t t + D reads: • Evolution in the imaginary time: We repeat the time step until the chemical potential will stop to change on the first 10 digits. The algorithm without the step with forcing the phase is a well known and mathematically understood procedure, which leads to the condensate wave-function. In case with forcing the phase, and when the phase is appropriate, procedure leads to states which also do not change in the evolution. We tested the method for a gas with contact interaction only, generating single and a few Shabat-Zakharov dark and gray solitons.
Densities computed with this method for three different aspect ratios 10 , 3 l = -10 −2 and 10 −1 are shown in figure 3. We compare these profiles with the density of a dark soliton from the purely contact-interacting case [23] marked in figure 3 with the thick gray line. As in both cases, a contact interacting gas and the quasi1D dipolar gas considered here, the lowest lying excitations are phonons, then comparing their speed of sound we translate the dipolar coupling coefficient g dd into an effective coupling coefficient g . eff We find g g 3 . The Gross-Pitaevskii equation is supposed to be valid in the weakly interacting limit, namely when na 1, eff  which implies l . x One can see that the more elongated the trap, the closer the result is to the soliton in the contact interacting gas, marked in figure 3 with the thick gray line. To understand this observation mathematically we look at the system in the momentum space. The range of the effective dipolar potential in the momentum space is equal to l 1 , equal in the box units to 1 . l Hence, in the limit of a very elongated trap, 0 l  the range of the effective potential is going to , ¥ and the evolution is dominated by its values relatively close to k = 0. The derivative of the effective potential equation (8) at k = 0 vanishes, which means that the potential is flat at the origin. Thus locally, where this 'local' region spans even to high momenta, the interaction potential looks like a constant function, as it would be in the contact interacting case. The ratio between the typical width of the soliton and the range of the potential, written in the momentum space and in the box units, is equal to ng .
Thus if the trap is very elongated and the coupling strength is small, we expect that the dipolar soliton will be close to the Shabat-Zakharov solution and moreover that our system will evolve similarly to the gas with contact interactions only. The presented solutions are only candidates for solitons, having similar jump of the phase as the Shabat-Zakharov solitons. The crucial feature of solitons is their ability to propagate without changing their shape. To verify if this is also the case here, we have set as initial states the results of the imaginary time evolution, and trace their dynamics obtained via numerical solution of GPE. An example is given in figure 4, where we unrotated the constant speed due to the linear dependence of the phase. The evolution time in this figure is about 30 times longer than a time t phon at which a phonon would travel across the whole box. When testing the system in a wide parameter range and for a longer evolution time we always observe the dips which propagate without noticeable distortion.

Two solitons
The real solitons obey an equivalent of the superposition principle: although the system is non-linear one can construct many solitonic solutions based on a single one. As shown in the paper of Shabat and Zakharov [23] in the contact interacting gas two solitons, when colliding, do not change the shape. To test if the dipolar solitons can behave similarly, we generate and propagate pairs of them. We follow the method described for a single soliton, but this time we enforce two jumps of the phase, at positions z 1 and z 2 , and without linear background Having such constraint on the phase we reach with the help of the imaginary time evolution two solitons and then we propagate them using GPE. An example of such dynamics is shown in figure 5. The immediate conclusion is that the dipolar solitons interact; in the example shown they attract each other. Hence, the excitations we study are not solitons in the strict mathematical sense. Even so, these excitations resemble the shape for a long evolution times, thus we will use the name soliton-like or, for brevity, soliton.
To understand the motion of such a pair we monitor the position of one of the solitons versus time, z t , sol ( ) defined as the position of the local minimum in the density. Then we numerically compute both the velocity and the acceleration a t . sol sol D = | | We repeat this procedure for many initial states. One can see that results of many simulations are lying on a common curve; see figure 6. This suggests that we should use the following classical mechanics picture, so in analogy to many cases for dark solitons in the contact interacting gas [24,25]    In both cases the inter-solitonic potential has a repulsive core and tail scaling with the distance as z 1 . 3 We tested this classical picture starting again from different initial states and watching their evolution. Samples of this investigation are given in the bottom panels in figures 7 and 8. Two solitons placed close to each other strongly repel each other. As the inter-solitonic potential has a local minimum, thus there is a possibility for bound states as shown in the panels marked with the letter A. Two far away lying solitons, such as those in panel B, hardly see each other, but after a longer time they finally start oscillations in the common potential.
In figures 7 and 8 we compare the resulting potential V m sol just with the rescaled density of a gas with a single soliton (namely the corresponding density plotted in figure 3), z , where 0 r is the density far from the soliton, approximately equal to l 1 5 . The agreement between solitonic shape and intersolitonic potential is surprisingly good. This indicates the origin of the inter-solitonic potential: each soliton is simply moving in the mean field potential created by its partner. As the dark solitons have negative masses, they are climbing towards the local maxima of the density.
The attentive reader might notice some peculiarity in the bottom right panel of figure 9 . There, the solitons seem to slightly accelerate. This counter-intuitive effect is not a numerical artefact. To illustrate clearly what is happening we present the evolution for the extreme conditions, ng 0.1 1000, dd l = = in figure 9. Then one can see that the collision is not elastic: colliding solitons radiate phonons. In consequence the energy of the solitons decays, and the dark solitons are transformed into the gray solitons. As in the case with contact interacting gas, the gray solitons have some characteristic speed. This extra velocity is thus gained by the solitons due to the inelastic collision.
Hence, the dissipation during the collision leads in fact to the acceleration of solitons, as clearly shown in figure 9. Actually this effect, although not directly described, seems to exist for the optical solitons with nonlocal   5 We compute 0 r numerically, it is not a fitting parameter. nonlinearity studied in [21] (figure 3 therein). Until collision the movement of a pair of solitons can be captured by a classical model. Then, collision goes clearly beyond the classical picture: it is an event when the part of the energy is transferred from solitons into initially unoccupied excitations, i.e. phonons. What is the fate of these solitons?
To find the answer we let the system evolve for a long time. During the evolution we monitored minimal density in the system, called here a grayness H.
) [6]. The evolution is shown in figure 10. As we started with two black solitons, initially H = 0. Then, at short times, the grayness H increases, first quadratically and then exponentially. At a very long time the quantity H starts to saturate. We understand the dynamics in terms of interaction between the subsystem (solitons) with an environment (all phonons and other quasiparticles), taking the idea from the paper [26]. Initially the environment is empty: this is the regime of spontaneous emission, where a coupling between subsystem and the unoccupied modes leads to a slow dissipation. Then, we enter the regime of stimulated emission, where the modes after the stage of spontaneous emission are already occupied, and hence the rate of phonon emission rapidly grows.
To model the first stages of the evolution we use the analogy of a two-level atom interacting with the electromagnetic field. In that case, the energy remaining in the atom decays initially as E t t 1 , 0 0 2 -( ( ) ) but for longer times the decay is exponential. We model it with E t E A Bt sinh , 0 2 = -( ) ( ) and use it as an ansatz for the energy remaining in the solitons, keeping E 0 , A and B as fitting parameters. The comparison is shown in figure 10.
If the environment would be infinitely large then eventually we would lose the energy from the subsystem completely. However, here we deal with a small environment, in which the phonons always overlap with solitons. This allows for the dynamical equilibration, in which the amount of the energy radiated to the phonons is equal to the energy absorbed, as in the qubit interacting with a thermal reservoir. Indeed we see in the inset of figure 10 that even after a long evolution time the two solitons survive. At this stage, around time t = 40, there are phonons and quasiparticles visible as the shallow dips moving with or faster than the speed of sound (compare with the superimposed red line in the inset of figure 10), the two gray but deep solitons being reminiscent of the initial two dark solitons, and a few shallow gray solitons which appear next to phonons. The large fluctuations of H are due to overlaps between phonons and solitons, but also due to statistical fluctuation of the flow of the energy between solitons and phonons.

Thermal solitons
As mentioned in the introduction, the ultracold gas with contact interaction possesses not only phonons but also dark solitons [11]. These solitons are the elementary excitations forming the second branch of excitations derived in 1963 by Elliot Lieb [3]. We would like to check if non-Bogoliubov excitations exist also in a quasi 1D dipolar gas. To fully verify it, one should solve a many-body problem, an equivalent of the Lieb-Liniger model for a dipolar gas. This task is definitely beyond this paper. Instead, still relying on the mean field model, we can at least check if the solitons, as candidates for type II excitations in a dipolar gas, appear spontaneously at finite temperatures and if their number obeys the thermal statistics. At first glance, the presence of such solitons seems unlikely: the system is not integrable and we have already shown that the dipolar soliton-like excitations are not robust, as they collide inelastically. On the other hand, we also demonstrated that two solitons immersed in the gas full of phonons can stay in the dynamical equilibrium. The latter observation encourages us to look also at the thermal equilibrium.
As in [11], we use the classical field approximation (CFA), in which one replaces the quantum field operator z Ŷ ( ) with the classical field z .
Y( ) This method is commonly used at zero-temperature, where z Y( ) is interpreted as the macroscopically occupied (by all particles) single-particle orbital, called the GPE wave function. The idea of using the same replacement to describe a gas at non-zero temperature relies on the fact that for the degenerate gas there are many substantially occupied single-particle orbitals, which sum up to a classical field z . There is variety of different classical field methods devised to describe ultracold gases [27]. Here we use a simple version described in [28]. First, we generate many initial states using a method based on the Bogoliubov description. As our method of generating thermal samples is not sufficiently good, we evolve the field with the help of GPE, we trace the system for sufficiently long time and compute time averages.
More precisely, in each realization, the initial wave-function reads )is the chemical potential. The cut-off k max is chosen such that with the classical field method one would recover the quantum results at least for the ideal gas [29,30]. The expansion coefficients k a should be drawn from the thermal classical distribution. Even on the classical level this drawing can be a cumbersome task. The details of our procedure for drawing the initial set of k a are given in the appendix. As the procedure gives an unsatisfactory effect (especially for higher temperatures), we have to propagate each guessed classical field over a long time to allow for thermalization. Finally, after this auxiliary stage of thermalization, in the next stage we further propagate the state to perform averaging over time. The short-time windows of the evolution during the last stages are illustrated in figure 11. We see there three examples of the evolution of the density of the gas, with the fraction of the 0-momentum mode from 0.68 (top panel), through 0.17 (middle panel), to 0.08 (bottom panel).
The speed of sound is depicted with the solid black line. In the middle and the bottom panels one can see dips in the density, propagating much slower than the speed of sound. We checked that these dips were stable even at long timescales. Thus we expect that they are dark solitons.
To convince ourselves we prepare the following test: for each thermal sample we performed a long evolution, during which we average over time the fraction of the least energetic 0-momentum mode N , ) and the number of the deepest dips, with the density below L 0.05 . The last quantity, called the number of dark solitons, is presented in figure 12. We estimated the temperature of each sample using formulas for the average number of motionless particles, i.e.
N , for the ideal gas [31]. As the crosscheck, we compute what would be fluctuation of the latter average in the ideal gas and compare it with the previously numerically found N .
The agreement was within a few percent, hence we concluded that 1) the interactions in the gas are still so small that the statistics seem close to the ideal gas case and 2) our estimation for temperature is reasonable within a few percent. We also estimated the energy of dark solitons, relying on the formulas for the gas with contact interactions  figure 12. We also repeated all the steps for the gas with the contact interaction only.

Conclusions
We investigated a dipolar gas confined in the potential having the shape of a long tube, with the dipole moments polarized perpendicular to the long axis. In this geometry we find robust localized excitations, very similar to the Shabat-Zakharov solitons appearing in the gas with contact interaction only. In contrast to the contact interacting case, pairs of these soliton-like dipolar excitations mutually interact, and are able to form two-soliton bound state. After thousands of subsequent collisions they can reach dynamical equilibrium with the already emitted quasiparticles. Finally we show that these excitation appear spontaneously at finite temperatures and their number obeys thermal statistics. The latter observation is closely related to the contact interacting case, where the appearance of the dark solitons at the non-zero temperatures is evidence that they are the non-Bogoliubov elementary excitations predicted by Elliot Lieb in 1963. There are many open questions concerning these excitations. One should check their stability in the real three-dimensional calculation. It is interesting if they exist also in the roton regime. There is the question of the effect of a harmonic trap along the long axis (here the box with periodic boundary was assumed). One also needs benchmarks from other methods used for the finite-temperature simulations, as they are quite developed for dipolar gases [32][33][34][35]. Probably the most important task would be a verification if these excitations are really elementary and if they are present in the many-body model. To generate the initial state for a given temperature and for a given number of atoms N we follow the steps: 1) For k 0 ¹ we draw the number of quasiparticles n k k 2 b = | | using the exponential distribution.
2) For k 0 ¹ we draw a phase k f with the uniform distribution on the interval , .   and ng 0. con = The green points are computed analogously to the dipolar case, but for ng 12000 con = and g 0. dd = Finally the blue dashed line is the expected number of solitons assuming the equipartition of energy and assuming dark solitons as quasi-particles, i.e. the number of solitons is approximated with k T E , B s o l where T is the estimated temperature and E sol is the estimated energy of the dark soliton. The details of the estimation are in the text. In each case we assumed N = 1000 atoms. The results for a dipolar gas and the gas with contact interaction comes from time averaging performed on 1000 realizations of the classical field, each drawn as described in the appendix.