Dipolar dark solitons
Abstract
We numerically generate, and then study the basic properties of dark solitonlike 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.
pacs:
03.75.Hh , 03.75.Lm , 05.45.YvI Introduction
The main motivation to our work comes from recent progress concerning the role of dark solitons in the physics of the gas with contact interaction only. First, let us remind, that solitons are localized waves propagating through nonlinear 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 nonlinear one dimensional Schrödinger equation:
(1) 
in which one can find bright (the case with attracting forces, ) and dark (the repulsive case ) solitons.
The Eq. (1) applies, among others, to the dynamics of bosons interacting via a like potential, i.e. the interaction potential between two bosons is equal to . In this context, the Eq. (1) is only an approximation valid in the variously defined weak interaction limit Lieb et al. (2005). Such a gas of bosons, can be described by more fundamental manybody model, the LiebLiniger model, which has been used to find rigorously the ground state Lieb and Liniger (1963) as well as the elementary excitations Lieb (1963). 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, within research performed by the next generations of physicists turned out to be solitonic branch Kulish et al. (1976); Kolomeisky et al. (2000); Jackson and Kavoulakis (2002); Kanamoto et al. (2010); Karpiuk et al. (2015); Syrwid and Sacha (2015). It is now clear that the dark solitons, as the elementary excitations, appear spontaneously in a gas at finite temperature Karpiuk et al. (2012).
The equation (1), is not capturing many important and common situations, even in the weakly interacting limit. For example, it is neglecting possible external traps in which the gas is confined, and it is assuming contact interaction only. In particular, it misses the gas with a dipolar interaction, which is recently studied in more and more laboratories Griesmaier et al. (2005); Aikawa et al. (2012); Lu et al. (2011); Ni et al. (2008). Among many motivations to study the dipolar gas is definitely its unusual excitation spectrum, which in certain situations, close to instabilities, can have local minimum for nonzero quasi momenta, so called roton spectrum Santos et al. (2003). It seems natural to ask if here is a place for the second branch and from what excitations it can be build. Up 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 reach (see parts of the review Lahaye et al. (2009)), 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 Adhikari (2014); Cuevas et al. (2009). Surprisingly, the life time of the dark solitons can be increased, with some tricks even in the higher dimensional space Nath et al. (2008). The dark solitons in nonlocal nonlinear media were the subject of research in optics. It has been found that two optical dark solitons interact with each other via nontrivial potential, they can even form a bound state Kong et al. (2010). In this context, although the effective forces between photons can be nonlocal, it is always assumed that they are not longrange.
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 Sec. II. After presentation of the basic features of single soliton (shapes, dynamics) in Sec III, we discuss a motion of two dipolar solitons in Sec. IV. Sec. V is devoted to the analysis of the thermal samples.
Ii The Model
We assume that the dipolar gas is confined in the potential
(2) 
with a tight harmonic trap in the and the directions. The space in the direction is assumed to be finite, with the length and with periodic boundary conditions. These conditions impose a quantization of the momenta in this direction, is a multiple of .
Throughout the paper we use the box units, where , , and are the units of length, time, energy and temperature, respectively. We assume that the dipole moments of the constituent atoms are polarized along the axis , namely when the dipolar potential in the momentum space has the form ^{1}^{1}1We give the dipolar interaction potentials in the momentum space because this form is more useful in numerical treatment.
(3) 
where is a dipolar coupling strength. Thus, the atoms move on the circumference of a circle, having the dipole moments perpendicular to the circleplane. In this configuration dipoles mostly repel each other, what favors the dark solitons. The geometry of our problem is sketched in Fig. 1.
We will investigate the system within the mean field model. The timedependent mean field wavefunction obeys the three dimensional GrossPitaevskii equation (GPE)
We assume that the transverse confinement is so tight, that the wavefunction stays in the lowest energy level in the and directions, namely when the following singlemode approximation holds:
(4) 
where is the Gaussian wave function:
(5) 
where is the transverse oscillator width. The necessary condition underlying this approximation is that both, the thermal energy and the interaction energy, are much below the energy of the first excited state in the transverse direction.
(6) 
Assuming such a situation we multiply the three dimensional GPE by , integrate the result with respect to and and finally we reach the one dimensional effective GPE
(7) 
The Fourier transform of the effective one dimensional dipolar potential reads Nath et al. (2008); Cai et al. (2010).
(8) 
where the constant term with mimics a contact interaction, we introduced one dimensional dipolar coupling strength and is equal to the aspect ratio of the trap. The function which appears in Eq. (8) is equal to , where Ei is the exponential integral. In this paper we aim at studying the role of pure dipolar interactions only, hence from now on, we will assume ^{2}^{2}2This can be done via Feshbach resonances, setting .
The spatial representation of this potential is visualized in Fig. 2.
It is worth stressing that Eq. (7) describes not a movement of single dipoles, but rather the movement of slices of the dipolar gas as marked in Figs. 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 . The situation changes, when the distance between them is of the order of the transverse width . 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.
Iii 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 the 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:
(9) 
such that the phases at the edges of the box are equal. Due to this additional term the gas is moving with the speed . When showing the result for the dynamics, we would unrotate this motion.
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). Densities computed with this method for three different aspect ratios , and are shown in Fig. 3. We compare these profiles with the density of a dark soliton from purely contact interacting case Zakharov and Shabat (1973) marked in Fig 3 with the thick gray line. As in both cases, a contact interacting gas and the considered here quasiD dipolar gas, the lowest lying excitations are phonons, then comparing their speed of sound we translate the dipolar coupling coefficient into an effective coupling coefficient . We find . The speed of sound in the dipolar media is equal to . The width of the ShabatZakharov soliton is proportional to the healing length, , where is the effective scattering length defined by the relation . The GrossPitaevskii equation is supposed to be valid in the weakly interacting limit, namely when , what implies .
One can see that the more elongated is the trap, the closer is the result to the soliton in contact interacting gas, marked in Fig. 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 , equal in the box units to . Hence, in the limit of a very elongated trap, the range of the effective potential is going to , and the evolution is dominated by its values relatively close to . The derivative of the effective potential Eq. (8) at vanishes, what 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 . Thus if the trap is very elongated and the coupling strength is small, we expect that the dipolar soliton will be close to the ShabatZakharov 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 ShabatZakharov solitons. The crucial feature of solitons is their ability to propagate without changing their shape. To verify if this is the case also here, we were setting 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 Fig. 4, where we unrotated the constant speed due to the linear dependance of the phase. The evolution time in this figure is about times longer than a time 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 noticable distortion.
Iv Two solitons
The real solitons obey an equivalent of the superposition principle: although the system is nonlinear one can construct many solitonic solution based on a single one. As shown in the paper of Shabat and Zakharov Zakharov and Shabat (1973) 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 and , and without linear background
(10) 
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 Fig 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 though, these excitations resemble the shape for a long evolution times, thus we will use the name solitonlike or, for brevity, even the name soliton.
To understand the motion of such a pair we monitor the position of one of the solitons versus time: . Then we numerically compute both the velocity and the acceleration . Combining these quantities we obtain acceleration versus intersolitonic distance , where . We repeat this procedure for many initial states. One can see that results of many simulations are lying on a common curve, see Fig. 6. It suggests us to use the following classical mechanics picture, in analogy to many cases for dark solitons in the contact interacting gas Theocharis et al. (2005); Busch and Anglin (2000) we treat each soliton as a pointlike (negative) mass obeying the Newton law. Assuming , we can numerically integrate to extract intersolitonic potential divided by a mass .
We show the results for two cases: in Fig. 7 for very elongated trap with and weak dipolar interaction and in Fig. 8 for and . In both cases the intersolitonic potential has a repulsive core and the tail scaling with the distance as . We tested this classical picture starting again from different initial states and watching their evolution. Samples of this investigation are given at the bottom panels in Figs. 7 and 8. Two solitons placed close to each other strongly repel each other. As the intersolitonic 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, like in the panel marked with the letter B, they hardly see each other, but after longer time they finally start oscillations in the common potential.
In Figs. 7 and 8 we compare the resulting potential just with the rescaled density of a single soliton , where is the density far from the soliton, approximately equal to ^{3}^{3}3We compute numerically, it is not fitting parameter.. The agreement between solitonic shape and intersolitonic potential is surprisingly good. This indicates the origin of the intersolitonic 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 could notice some peculiarity in Fig. 9 bottom right panel. There, the solitons seem to slightly accelerate. This counterintuitive effect is not a numerical artefact.
To illustrate clearly what is happening we present the evolution for the extreme conditions, , in Fig. 9. Then one can see that the collision is not elastic: colliding solitons radiate phonons. In consequence the energy of 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 acceleration of solitons, as clearly shown in Fig. 9. Actually this effect, although not directly described, seems to exist for the optical solitons with nonlocal nonlinearity studied in Kong et al. (2010) (Fig. 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 . In the case of a gas with contact interactions only, this parameters can be translated into the ShabatZakharov soliton speed, equal to , and its energy, equal to Jackson and Kavoulakis (2002). The evolution is shown in Fig. 10. As we started with two black solitons, initially . Then, at short times,the grayness increases, first quadratically and then exponentially. At a very long time the quantity 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 Bach et al. (2002). 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. If the environment would be infinitely large then eventually we would loose 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 phonons is equal to the energy absorbed, as in the qubit interacting with a thermal reservoir. Indeed we see in the inset of Fig. 10, that even after a long evolution time the two solitons survive. The large fluctuations of are due to overlaps between phonons and solitons, but also due to statistical fluctuation of the flow of the energy between solitons and phonos.
V Thermal solitons
As mentioned in the introduction, the ultracold gas with contact interaction possesses not only phonons but also dark solitons Karpiuk et al. (2012). These solitons are the elementary excitations forming the second branch of excitations derived in 1963 by Elliot Lieb Lieb (1963). We would like to check if nonBogoliubov excitations exist also in a quasi 1D dipolar gas. To fully verify it, one should solve a manybody problem, an equivalent of the LiebLiniger 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 the first glance, the presence of such solitons seem unlikely: the system is not integrable and we have already shown that the dipolar solitonlike 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 Karpiuk et al. (2012), we use the classical field approximation (CFA), in which one replaces the quantum field operator with the classical field . This method is commonly used at zerotemperature, where 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 nonzero temperature relies on the fact that for the degenerate gas there are many substantially occupied singleparticle orbitals, which sum up to a classical field .
To compute the values of target variables averaged over a statistical ensemble one performs averaging over the appropriately chosen set of the classical fields or/and timeaveraging over the fields propagated with the help of the GPE equation or its generalization.
There is variety of different classical field methods devised to describe ultracold gases publication (2013). Here we use a simple version described in Brewczyk et al. (2007). 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 wavefunction reads
(11) 
where the wave vector is with an integer. The sum over wavevectors is limited by the with , where is the temperature and is the chemical potential. The cutoff is chosen such that with the classical field method one would recover the quantum results at least for the ideal gas Witkowska et al. (2009, 2010). The expansion coefficients 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 of drawing initial set of are given in the Appendix. As the procedure gives 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 shorttime windows of the evolution during the last stages are illustrated in Fig 11. We see there three examples of the evolution of the density of the gas, with the fraction of the momentum mode from (top panel), through (middle panel), (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 momentum mode , the fluctuation of this fraction and the number of the deepest dips, with the density below . The last quantity, called the number of dark solitons, is presented in Fig. 12. We estimated the temperature of each sample using formulas for average number of motionless particles, i.e. , for the ideal gas Wilkens and Weiss (1997). 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 . The agreement was within a few percent, hence we concluded that: 1) the interactions in the gas are still so small that the statistics seems 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 only Jackson and Kavoulakis (2002) as done in the previous sections. This gives the estimation for the energy of the dipolar dark soliton (in box units). Finally, having the estimated energy of the dark soliton, and estimated temperature of each sample we computed the estimated number of dark solitons , assuming their equipartition as it should be for quasiparticles. This expected number of solitons is given by the blue dashed line in Fig. 12. We also repeat all the steps for the gas with the contact interaction only.
Vi Conclusions
We investigated the 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 solitonlike dipolar excitations mutually interact, and are able to form twosoliton bound states. Moreover, the collisions between such pairs are not elastic: during collision they emit phonons, hence they loose energy and because of this they accelerate. After thousands of subsequent collisions they can reach dynamical equilibrium with the already emited 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 nonzero temperatures is an evidence that they are the nonBogoliubov 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 needs also benchmarks from other methods used for the finitetemperature simulations, as they are quite developed for the dipolar gases Ticknor (2012); Bisset et al. (2012); Blakie et al. (2009); Bisset et al. (2011). Probably the most important task would be a verification if these excitations are really elementary and if they are present in the manybody model.
Acknowledgements.
We are grateful to T. Karpiuk, T. Sowiński, P. Deuar and M. Gajda for valuable discussions. The work was supported by the (Polish) National Science Center Grant No. DEC2012/04/A/ST2/00090.Vii Appendix: Generating thermal states within classical field approximation
In the Bogoliubov method one approximates the grand canonical Hamiltonian with , where is the Bogoliubov spectrum and the operators and are the annihilation and creation bosonic ladder operators for quasiparticles, respectively. The operator , annihilating a quasiparticle with quasimomenta is related to the bosonic operators annihilating a particle with a given momenta via relation:
(12) 
where and . In the classical version, we replace , , and with the complex numbers denoted as , , and , respectively. We introduce the classical probability distribution with the partition function defined such that . Within this approximation the number of quasiparticles with momenta is given by the exponential distribution with the parameter . To generate the initial state foe a given temperature and for a given number of atoms we follow the steps:

For we draw the number of quasiparticles using the exponential distribution.

For we draw a phase with the uniform distribution on the interval .

We construct the amplitude .

We compute the amplitudes .

We compute the amplitude of the momentum mode, .

We construct the classical field according to Eq. (12).
References
 Lieb et al. (2005) E. H. Lieb, R. Seiringer, J. P. Solovej, and J. Yngvason, The Mathematics of the Bose Gas and its Condensation, vol. 34 (Birkhäuser Basel, 2005).
 Lieb and Liniger (1963) E. H. Lieb and W. Liniger, Physical Review 130, 1605 (1963).
 Lieb (1963) E. H. Lieb, Phys. Rev. 130, 1616 (1963), URL http://link.aps.org/doi/10.1103/PhysRev.130.1616.
 Kulish et al. (1976) P. Kulish, S. Manakov, and L. Faddeev, Theoretical and Mathematical Physics 28, 615 (1976), ISSN 00405779, URL http://dx.doi.org/10.1007/BF01028912.
 Kolomeisky et al. (2000) E. B. Kolomeisky, T. J. Newman, J. P. Straley, and X. Qi, Phys. Rev. Lett. 85, 1146 (2000), URL http://link.aps.org/doi/10.1103/PhysRevLett.85.1146.
 Jackson and Kavoulakis (2002) A. D. Jackson and G. M. Kavoulakis, Phys. Rev. Lett. 89, 070403 (2002), URL http://link.aps.org/doi/10.1103/PhysRevLett.89.070403.
 Kanamoto et al. (2010) R. Kanamoto, L. D. Carr, and M. Ueda, Phys. Rev. A 81, 023625 (2010), URL http://link.aps.org/doi/10.1103/PhysRevA.81.023625.
 Karpiuk et al. (2015) T. Karpiuk, T. Sowiński, M. Gajda, K. Rzążewski, and M. Brewczyk, Phys. Rev. A 91, 013621 (2015), URL http://link.aps.org/doi/10.1103/PhysRevA.91.013621.
 Syrwid and Sacha (2015) A. Syrwid and K. Sacha, Liebliniger model: emergence of dark solitons in the course of measurements of particle positions (2015), eprint arXiV:1505.06586.
 Karpiuk et al. (2012) T. Karpiuk, P. Deuar, P. Bienias, E. Witkowska, K. Pawłowski, M. Gajda, K. Rzążewski, and M. Brewczyk, Phys. Rev. Lett. 109, 205302 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.109.205302.
 Griesmaier et al. (2005) A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau, Phys. Rev. Lett. 94, 160401 (2005), URL http://link.aps.org/doi/10.1103/PhysRevLett.94.160401.
 Aikawa et al. (2012) K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.108.210401.
 Lu et al. (2011) M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, Phys. Rev. Lett. 107, 190401 (2011), URL http://link.aps.org/doi/10.1103/PhysRevLett.107.190401.
 Ni et al. (2008) K.K. Ni, S. Ospelkaus, M. H. G. de Miranda, a. Pe’er, B. Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Julienne, D. S. Jin, and J. Ye, Science (New York, N.Y.) 322, 231 (2008), ISSN 10959203, URL http://www.ncbi.nlm.nih.gov/pubmed/18801969.
 Santos et al. (2003) L. Santos, G. V. Shlyapnikov, and M. Lewenstein, Phys. Rev. Lett. 90, 250403 (2003), URL http://link.aps.org/doi/10.1103/PhysRevLett.90.250403.
 Lahaye et al. (2009) T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Reports on Progress in Physics 72, 126401 (2009).
 Adhikari (2014) S. K. Adhikari, Phys. Rev. A 89, 043615 (2014), URL http://link.aps.org/doi/10.1103/PhysRevA.89.043615.
 Cuevas et al. (2009) J. Cuevas, B. A. Malomed, P. G. Kevrekidis, and D. J. Frantzeskakis, Phys. Rev. A 79, 053608 (2009), URL http://link.aps.org/doi/10.1103/PhysRevA.79.053608.
 Nath et al. (2008) R. Nath, P. Pedri, and L. Santos, Phys. Rev. Lett. 101, 210402 (2008), URL http://link.aps.org/doi/10.1103/PhysRevLett.101.210402.
 Kong et al. (2010) Q. Kong, Q. Wang, O. Bang, and W. Krolikowski, Phys. Rev. A 82, 013826 (2010), URL http://link.aps.org/doi/10.1103/PhysRevA.82.013826.
 Cai et al. (2010) Y. Cai, M. Rosenkranz, Z. Lei, and W. Bao, Phys. Rev. A 82, 043623 (2010), URL http://link.aps.org/doi/10.1103/PhysRevA.82.043623.
 Zakharov and Shabat (1973) V. Zakharov and A. Shabat, Zh. Eksp. Teor. Fiz. 64, 1627 (1973).
 Theocharis et al. (2005) G. Theocharis, P. Schmelcher, M. K. Oberthaler, P. G. Kevrekidis, and D. J. Frantzeskakis, Phys. Rev. A 72, 023609 (2005), URL http://link.aps.org/doi/10.1103/PhysRevA.72.023609.
 Busch and Anglin (2000) T. Busch and J. R. Anglin, Phys. Rev. Lett. 84, 2298 (2000), URL http://link.aps.org/doi/10.1103/PhysRevLett.84.2298.
 Bach et al. (2002) R. Bach, M. Trippenbach, and K. Rza¸żewski, Phys. Rev. A 65, 063605 (2002), URL http://link.aps.org/doi/10.1103/PhysRevA.65.063605.
 publication (2013) J. publication, Stochastic approaches for degenerate Bose gases at finite temperature in the canonical ensemble (Imperial College Press, 2013).
 Brewczyk et al. (2007) M. Brewczyk, M. Gajda, and K. Rzążewski, Journal of Physics B: Atomic, Molecular and Optical Physics 40, R1 (2007).
 Witkowska et al. (2009) E. Witkowska, M. Gajda, and K. Rzążewski, Phys. Rev. A 79, 033631 (2009), URL http://link.aps.org/doi/10.1103/PhysRevA.79.033631.
 Witkowska et al. (2010) E. Witkowska, M. Gajda, and K. Rzążewski, Optics Communications 283, 671 (2010), ISSN 00304018, URL http://www.sciencedirect.com/science/article/pii/S0030401809010724.
 Wilkens and Weiss (1997) M. Wilkens and C. Weiss, Journal Of Modern Optics 44, 1801 (1997).
 Ticknor (2012) C. Ticknor, Phys. Rev. A 85, 033629 (2012), URL http://link.aps.org/doi/10.1103/PhysRevA.85.033629.
 Bisset et al. (2012) R. N. Bisset, D. Baillie, and P. B. Blakie, Phys. Rev. A 86, 033609 (2012), URL http://link.aps.org/doi/10.1103/PhysRevA.86.033609.
 Blakie et al. (2009) P. B. Blakie, C. Ticknor, A. S. Bradley, A. M. Martin, M. J. Davis, and Y. Kawaguchi, Phys. Rev. E 80, 016703 (2009), URL http://link.aps.org/doi/10.1103/PhysRevE.80.016703.
 Bisset et al. (2011) R. N. Bisset, D. Baillie, and P. B. Blakie, Phys. Rev. A 83, 061602 (2011), URL http://link.aps.org/doi/10.1103/PhysRevA.83.061602.