Dark solitons in the unitary Bose gas

We study the dilute and ultracold unitary Bose gas, which is characterized by a universal equation of state due to the diverging s-wave scattering length, under a transverse harmonic confinement. From the hydrodynamic equations of superfluids we derive an effective one-dimensional nonpolynomial Schr\"odinger equation (1D NPSE) for the axial dynamics which, however, takes also into account the transverse dynamics. Finally, by solving the 1D NPSE we obtain meaningful analytical formulas for the dark (gray and black) solitons of the bosonic system.


Introduction
The possibility of tuning the interaction between atoms near Feshbach resonances [1] almost at will, by precision control of external magnetic fields, is allowing the experimenters to explore the physics of ultracold gases in a wide range of values of the coupling constant characterized by the s-wave scattering length. The possibility to control the scattering length via Feshbach resonances makes ultracold atomic gases an excellent setting for studies of strongly correlated behavior. In particular it made possible to observe first a Fermi gas [2][3][4][5][6] and then also a Bose gas at unitarity [7][8][9][10], that is the limit in which the s-wave scattering length is infinitely large. In this so-called unitary regime, the interactions are as strong as allowed by quantum mechanics, and the physics cannot explicitly depend on the scattering length, leading to the possibility of new types of universal behavior. In the case of a Bose gas, the presence of three body recombination processes [12] complicates the experimental realization of the unitary limit, which can be held for a limited amount of time. Nevertheless, the peculiar properties of a unitary gas, which behaves in a universal manner and exhibits scale invariance, encourage us to study the physics of such systems more in depth.
In this paper we investigate solitonic configurations of the unitary Bose gas. More specifically, we analyze axial and transverse density profiles of dark solitons in bosonic alkali-metal atoms at unitarity under the action of a transverse harmonic confinement. A dark soliton is a self-bound solitary wave, whose existence is made possible by the interplay between the repulsive interaction and the presence of a phase gradient. A dark soliton is characterized by a local decrease in density with respect to a uniform background. The depth of the dark soliton crucially depends on the phase angle of the complex field which describes the Bose gas.
The realizability of black and gray solitons in ultracold gases is well documented in the weak-coupling limit [13,14] and in a unitary Fermi gas [15]. Moreover, the properties of black solitons have been theoretically investigated for superfluid fermions also in the BCS-BEC crossover [16][17][18][19]. It is then a natural question to ask, whether black and gray solitons can be realized also in the case of Bose gases at unitarity, and what are their specific properties. It is important to stress that dark (black and gray) solitons of the unitary Bose gas have not yet been observed experimentally nor theoretically analyzed. Our theoretical paper can help experimental groups because we are giving analytical and semi-analytical formulas for density profiles, widths, and velocities of these dark solitons. All these quantities are universal because, contrary to the dark solitons of weakly-interacting bosonic gases, they do not depend on the s-wave scattering length.
We first derive the equations which govern the dynamics of these systems. The resulting equation is a nonlinear Schrödinger equation for a complex order parameter ψ(r, t) = |ψ(r, t)|e iα(r,t) [27], where n(r, t) = |ψ(r, t)| 2 is the local number density of the bosonic system and v(r, t) = ( /m)∇α(r, t) is the local velocity with α(r, t) the phase angle of the complex field. Then, we find the solutions which are in the form of dark solitons, specifically as stationary objects in a moving frame.
Mimicking the experimental practice of confining ultracold gases in the minimum of a potential, we study the case in which the Bose gas is placed in a harmonic potential with cylindrical symmetry, which is an approximation to a more realistic cigar shaped potential [13,15]. A Gaussian variational approach, whose effectiveness in similar cases has already been proved [20][21][22][23][24][25], is then deployed to obtain an effective onedimensional nonpolynomial Schrödinger equation (1D NPSE) for the axial wavefunction. The relevance of using the 1D NPSE in the study of solitons is due to the fact that the obtained solutions are a reliable generalization of familiar strictly one-dimensional results. These generalized solutions take properly into account that the transverse width of the cigar-shaped bosonic cloud depends on the axial coordinate.
In particular, we focus on the case of weak transverse confinement, in which more naive attempts to reduce the dimensionality of the problem fail. We integrate analytically the 1D NPSE thanks to the presence of a constant of motion, leaving behind a solution in terms of an integral, which we finally evaluate by means of numerical techniques. We are then able to find the axial density profile, transverse width and phase of the bosonic system. Moreover, we find the relation between density at the minimum and velocity of the soliton. While doing so, we also find that the theory actually depends on just one free parameter. Finally, we comment on the limits of the weak-coupling approximation.

From Euler to Schrödinger
Euler equations describe the dynamics of a non-viscous and irrotational fluid, such as a superfluid Bose gas at zero temperature [32]. In presence of an external potential U(r), these equations are given by ∂n ∂t The equation of state of the unitary Bose gas at zero temperature is assumed to be where n is the density of particles and A = u 2 m , with u a universal and adimensional (positive) coefficient. In fact, we expect this to be the case in the unitary limit a s → ∞ because of dimensional analysis, taking into account that the only important length scale is the mean distance between atoms ≈ n − 1 3 . The value of u = ξ(6π 2 ) 2 3 5 6 has been derived theoretically in many ways [26,[28][29][30][31] and, depending on the procedure, values in the range ξ ≈ 0.4 ÷ 1.75 have been reported.
We shall now derive a nonlinear Schrödinger equation (NLSE) with a 4/3 nonlinearity exponent [27], starting from equations (1)-(2) and adding a term, which will account for quantum effects. This is done via the mapping where α(r, t) is a scalar field. Equation (5) is indeed fully justified by the fact the fluid is irrotational, i.e. ∇ ∧ v = 0. If we substitute (4)-(5) into (1)-(2) we obtain Let us integrate equation (7) and multiply both sides by |ψ| The integration constant C appeared, but we shall see that it is not going to affect the dynamics. In order to get the Schrödinger equation we now add the term − 2 2m ∇ 2 |ψ(r, t)| on the right hand side of equation (8). This is chosen to obtain the same kinetic term which appears in the weak-coupling limit P ∝ n. In that case we would obtain the Gross-Pitaevskii equation, whose validity for weak coupling is well established [33,34].
We can express the system of real equations as a single complex equation by adding equation (6) multiplied by i to equation (8). Then we multiply every term by e iα , to Finally, by defining last equation can be written as a NLSE C results in a potential energy offset. It can be removed by substituting ψ → e − i Ct ψ.
It is important to stress that equation (10) implies that the scalar field α(r, t) is an angle and the circulation of the velocity field v(r, t) = ( /m)∇α(r, t) is quantized.

Nonpolynomial Schrödinger equation
Equation (11) can also be derived minimizing the action functional with the constraint Let us suppose that the confinement potential is harmonic, with cylindrical symmetry along the z axis, and consider the following variational ansatz for ψ where the variational parameters are the axial wavefunction φ(z, t) and the width σ(z, t) of the transverse Gaussian wavefunction. The choice of a Gaussian function in the (x, y) plane is justified by the presence of a harmonic potential in the (x, y) plane. Only in the strictly one-dimensional case one has σ = a ⊥ with a ⊥ = /(mω ⊥ ) the characteristic length of the transverse harmonic confinement.
By plugging equation (15) into the action in equation (12) we obtain under the approximation of neglecting the space derivative of σ(z, t). In the case of weakly-interacting bosons, which are very well described in three dimensions by the Gross-Pitaevskii equation (i.e. the 3D Schrödinger equation with cubic nonlinearity) this approximation has been found to be quite good also when σ(z, t) depends strongly but monotonically on z [20,25].
Let us rewrite the equations in terms of adimensional quantities. We can do that with the substitutions Equation (17) then becomes while equation (18) becomes Remember that in this equation σ is adimensional: it is in units of the characteristic length a ⊥ of the transverse harmonic confinement. In equation (20) one can identify two regimes: a regime of strong transverse confinement where σ ≃ 1 and consequently the system is truly one-dimensional, but also a regime of weak transverse confinement where σ ≫ 1. Thus, the regime of weak transverse confinement is obtained by neglecting the first addend on the right hand side of equation (20), supposing that AN which plugged back in equation (19) In the same approximation we can neglect the term ∝ (AN where we defined Finally, let us express also σ in terms of γ It will prove useful to study also the limit of strong confinement, that is instead given by σ = 1:

Dark solitons in the 1D NPSE
A soliton is a solution to equation (23) of the form where f (ζ) ≥ 0 and α(z, t) are real functions. Once we substitute (28) in (23) we obtain the equations for f and α Upon inspection of equation (30), we find that the right hand side should depend only on ζ, so ∂α ∂z is only function of ζ. Let us assume the following α(z, t) = θ(z − vt) + β(t).
If we plug this back into equation (29)-(30) we get

Phase
Let us solve equation (33) in order to obtain θ ′ in terms of f . We can rewrite it as for f (ζ) = 0. After integration with respect to ζ we have where D is an integration constant. If we assume that and we get D = vf 2 ∞ . The expression for θ ′ we are looking for is then valid for f (ζ) = 0.

Let us plug equation (38) back into equation (32)
We can choose β in order to get simpler equations. Let us suppose so that equation (39) becomes Equation (41) can be written in the form of a Newton equation for a conservative force field The total energy K is then a conserved quantity for all ζ and energy balance gives us the following first order equation which is satisfied by the two branches This equation can be put in the integral form with the separation of variables, giving where, since the integral is positive, we have chosen the + sign when ζ − ζ 0 ≥ 0 and the − sign otherwise. Therefore, f (ζ) = f (|ζ − ζ 0 |) has even parity with respect to reflections around ζ 0 , and we just need to study the case ζ − ζ 0 ≥ 0.
We have now to find the values of the parameters µ, K, f ∞ . This can be done by choosing boundary conditions for f .

Black solitons
Let us find an odd parity solution (which features a node in the origin for the density n) with a positive horizontal asymptote at +∞. We set ζ 0 = 0, and since f is even, the phase θ will account for the change of sign. Then, the following boundary conditions are needed This is possible only if v = 0 (ζ = z) in order to get rid of the term in equation (45) which diverges as f → 0. The condition f ′′ (+∞) = 0 applied to equation (41) gives .
It is possible to get rid of f ∞ with the rescaling and the parameter Plots of the numerical integration of equation (54) are shown in figure 1, for three values of δ. Accordingly with those integrations, also σ(z) is computed and shown in the plot.
We may definez = √ δz andσ = σ √ δ , so that the theory has no free parameters. This way there is only one plot to make, but we want to explicitly show in the plots the effects of changing the interaction, so for the moment we are keeping z and σ.

Gray solitons
Let us find an even parity solution with a positive horizontal asymptote. We set ζ 0 = 0 once again, and require the following boundary conditions The condition f ′′ (+∞) = 0 applied to equation (41) gives f ′ (+∞) = 0 applied to equation (45) gives while by equating K inside equation (45) evaluated at 0 and at ∞, applying f ′ (0) = 0 we get a relation between v, γ and f 0 (58) Equation (47) then becomes Also in this case we can get rid of f ∞ in the same way as before, obtaining (60) where of course h 0 = f 0 f∞ and v s = 2 5 δ is the speed of sound, i.e. the propagation velocity of a infinitesimal perturbation in the fluid density.
Now the theory has two free parameters, namely h 0 and δ, while v is given in terms of them. In order to be more concise, we now defineζ = √ δζ,σ = σ so that we are left with a theory that depends only on h 0 (orṽ). The scaled velocityṽ ranges from 0 (black soliton) to 1 (sound wave). Plots of the numerical integration of equation (61) are shown in figure 2, for three values ofṽ. Accordingly with those integrations, alsoσ(ζ) is computed and shown in the plot.

Phase
In the case of the black soliton, equation (38) becomes θ ′ = 0. Since we want the solution to be odd, we choose In the general case, we can now compute h and substitute inside equation (38) to obtain θ ′ . Let us express also θ in terms ofζ: it follows that Notice that δ = γ|f ∞ | 4/5 , with γ given by equation (24) and f ∞ the bulk value of the soliton wavefunction.ṽ = v/v s is the velocity rescaled by the sound velocity v s .
Last equation upon integration yieldsθ. A plot of the numerical integration of equation (66) is shown in figure 3 for four values ofṽ, includingṽ = 0, which is the black soliton case.

Weak vs strong transverse confinement
In the previous sections we have mainly investigated the regime of weak transverse confinement where σ ≫ 1. See the discussion of equation (20). Let us now briefly analyze some properties of the regime of strong transverse confinement where σ ≃ 1. In this case the system is truly one dimensional and the equations are much simpler. By performing computations which are very close to the ones we have already shown in the weak confinement case, one can derive from (26) that in the case of a black soliton the following relation holds where δ SC = γ SC f where v SC s = 2 3 δ SC . In dimensional units, taking into account equation (18), denoting ρ 0 = ρ(0) the minimal axial density of the dark soliton and ρ ∞ = ρ(±∞) the bulk axial density of the dark soliton, one finds that under the condition ρ ∞ a ⊥ ≪ 1 the dark soliton is surely in the strictly 1D regime of strong transverse confinement at any point of the axial coordinate. Instead, under the condition ρ 0 a ⊥ ≫ 1 the dark soliton is surely in the 3D regime of weak transverse confinement at any point of the axial coordinate. Clearly, in the case of the black soliton, where ρ 0 = 0, near the minimum the bosonic cloud cannot be in the regime of weak transverse confinement.
It is important to stress that the comparison of our theoretical results with future experiments is constrained not only by the short lifetime of unitary bosonic gas due to fast three-body recombinations but also by the lifetime due to the snake instability [16][17][18][19]. A dark soliton, that is not strictly 1D, has a snaking transverse oscillation which breaks the axial symmetry and eventually destroys the soliton. Since the NPSE preserves axial symmetry, it is not suitable for investigating snaking oscillations and the dynamics of the dissolving dark soliton. However, it is possible to write a triaxial NPSE by using two transverse width σ x (z, t) and σ y (z, t) in the ansatz of equation (15); see [35] for details in the case of weakly-interacting bosons.

Conclusion
We have derived the dynamical equations that govern the motion of a unitary Bose gas at zero temperature. Using these equations, we have studied dark solitons in a unitary Bose gas, confined in a cylindrically symmetric potential. We have reduced the dimensionality of the 3D nonlinear Schrödinger problem by means of the 1D nonpolynomial Schrödinger equation approach, which keeps dynamically into account also the transverse width. We have solved the equations employing for the most analytical techniques, in order to find axial density, transverse width, phase and velocity of both black and gray solitons, in the weak and strong confinement regimes. We have found that the weak confinement approximation breaks down in the case of a black soliton.
Our theoretical predictions could be tested experimentally employing Bose gases of alkali-metal atoms at ultra-low temperatures, whose scattering length can be tuned by means of Feshbach resonance techniques to reach the unitary limit. Despite the difficulties which presently limit the time during which such systems can be observed due to three-body losses but also to the snake instability, we are confident that the recent developments in the experimental techniques are going to allow in the future to put our predictions at test.