Weakly bound solitons and two-soliton molecules in dipolar Bose-Einstein condensates

TL;DR

This paper derives analytic expressions for dipolar bright soliton interactions and explores weakly bound soliton molecules in dipolar Bose-Einstein condensates, combining variational methods with numerical simulations.

Contribution

It introduces a variational approach to model the evolution of dipolar soliton molecules and compares predictions with numerical simulations, demonstrating good agreement.

Findings

Analytic interaction potential for dipolar bright solitons.

Description of weakly bound two-soliton molecules.

Qualitative agreement between model predictions and simulations.

Abstract

Analytic expressions have been derived for the interaction potential between dipolar bright solitons and the binding energy of a two-soliton molecule. The properties of these localized structures are explored with a focus on their behavior in the weakly bound regime, with a small binding energy. Using the variational approach a coupled system of ordinary differential equations for the parameters of a soliton molecule is obtained for the description of their evolution. Predictions of the model are compared with numerical simulations of the governing nonlocal Gross-Pitaevskii equation and good qualitative agreement between them is demonstrated.

Equations

$$i{\psi }_t+\frac{1}{2}{\psi }_{xx}+q|\psi {|}^2\psi +g\psi (x,t){\int }_{-\infty }^{+\infty }R(|x-x^{\prime }|)|\psi (x^{\prime },t){|}^{2}d{x}^{\prime }=0,$$

$$R(x)=(1+2x^2)\exp (x^2)\mathrm{erfc}(|x|)-2{\pi }^{-1/2}|x|.$$

$$R(x)=({\sigma }^2{x}^2+1{)}^{-3/2}.$$

$$\mathcal{L}$$

$$\psi (x,t)=A\exp \left[-\frac{(x-\xi {)}^2}{2a^2}+ib(x-\xi {)}^2+iv(x-\xi )+i\varphi \right],$$

$$\frac{L}{N}=\frac{1}{2}{a}^2{b}_t-\frac{1}{2}{\xi }_t^2+{\varphi }_t+\frac{1}{4a^2}+a^2{b}^2-\frac{qN}{2\sqrt{2\pi }a}-\frac{gN}{2\sqrt{2}\sigma a}\mathcal{U}\left(\frac{1}{2},0,\frac{1}{2{\sigma }^2{a}^2}\right),$$

$$a_{tt}=\frac{1}{a^3}-\frac{qN}{\sqrt{2\pi }{a}^2}+\frac{gN}{\sqrt{2}\sigma }\frac{\partial F_s(a,\sigma )}{\partial a},$$

$$F_s(a,\sigma )=\frac{1}{a}\mathcal{U}\left(\frac{1}{2},0,\frac{1}{2{\sigma }^2{a}^2}\right).$$

$$a_{tt}=-\frac{\partial U}{\partial a},U(a)=\frac{1}{2a^2}-\frac{qN}{\sqrt{2\pi }a}-\frac{gN}{\sqrt{2}\sigma }{F}_s(a,\sigma ).$$

$$\psi (x,t)=A(x-\xi )\exp \left[-\frac{(x-\xi {)}^2}{2a^2}+ib(x-\xi {)}^2+iv(x-\xi )+i\varphi \right],$$

$$\frac{L}{N}=\frac{3}{2}{a}^2{b}_t-\frac{1}{2}{\xi }_t^2+{\varphi }_t+\frac{3}{4a^2}+3a^2{b}^2-\frac{3qN}{8\sqrt{2\pi }a}-\frac{gN}{16\sqrt{2\pi }}{F}_m(a,\sigma ).$$

$$a_{tt}=\frac{1}{a^3}-\frac{qN}{4\sqrt{2\pi }{a}^2}+\frac{gN}{24\sqrt{2\pi }}\frac{\partial F_m(a,\sigma )}{\partial a},$$

$$F_m(a,\sigma )$$

$$\frac{d{}^2{x}_c}{d{t}^2}$$

$$\Delta n_R$$

$$\frac{d^2{x}_c}{d{t}^2}=-\frac{4qN}{a^3\sqrt{\pi }}{x}_c{e}^{-4x_c^2/a^2}-\frac{3gN{\sigma }^2}{a\sqrt{\pi }}\underset{-\infty }{\overset{\infty }{\int }}\frac{t{e}^{-(t-2x_c{)}^2/a^2}}{({\sigma }^2{t}^2+1{)}^{5/2}}dt.$$

$$U(x)=\frac{qN}{2\sqrt{\pi }a}{e}^{-4x^2/a^2}+\frac{3gN{\sigma }^2}{4}\underset{-\infty }{\overset{\infty }{\int }}\frac{\tau \mathrm{erfc}[(\tau -2x)/a]}{({\sigma }^2{\tau }^2+1{)}^{5/2}}d\tau ,$$

$$E=\underset{-\infty }{\overset{\infty }{\int }}\left[\frac{1}{2}|{\psi }_x{|}^2+\frac{q}{2}|\psi {|}^4+\frac{g}{2}|\psi {|}^2\underset{-\infty }{\overset{\infty }{\int }}R(|x-x^{\prime }|)|\psi (x^{\prime }){|}^{2}d{x}^{\prime }\right]dx,$$

$$E_s=\frac{N_s}{4a_s^2}-\frac{q{N}_s^2}{2\sqrt{2\pi }{a}_s}-\frac{g{N}_s^2}{2\sqrt{2}\sigma }{F}_s(a_s,\sigma ),$$

$$E_m=\frac{3N_m}{4a_m^2}-\frac{3q{N}_m^2}{8\sqrt{2\pi }{a}_m}-\frac{g{N}_m^2}{16\sqrt{2\pi }}F(a_m,\sigma ),$$

$$\Delta E=E_m-2E_s.$$

$$\psi (x)=A_s\left[e^{-(x-x_0{)}^2/2a_s^2}-e^{-(x+x_0{)}^2/2a_s^2}\right],$$

$$A_s=\frac{2A_m{a}_m}{\sqrt{\pi }}{e}^{-2/\pi },a_s=\frac{\pi a_m}{16}{e}^{4/\pi },x_0=\frac{2a_m}{\sqrt{\pi }}.$$

$$\Delta E$$

Full text

††journal: J. Phys. B

1 Introduction

The realization of Bose-Einstein condensation (BEC) of atomic species with large natural magnetic moments, such as chromium [1, 2], dysprosium [3, 4], erbium [5] and heteronuclear dipolar quantum mixtures of erbium and dysprosium atoms [6] has initiated new directions of research in the field of ultra-cold quantum gases. The physical properties of dipolar condensates are essentially different [7] as compared to other members of the BEC family with contact atomic interactions [8, 9]. The differences originate from the long-range and anisotropic dipole-dipole interaction potential between atoms in the gas, which decay with distance as $\sim 1/x^{3}$. In contrast to this, the atoms in non-dipolar BEC interact only when they are at the same spatial location, the so called contact interactions. The effects of the anisotropic and long range dipolar atomic interactions were observed in diverse forms, such as in intriguing collapse mechanism [10], emergence of a roton-like branch in the excitation spectrum [11, 12], modification of superfluid properties [13, 14], peculiar density fluctuations in the ground state [15], Faraday patterns sensitive to the roton mode [16] and other novel many-body effects in dipolar quantum gases [17, 18]. The properties of dipolar BEC in quasi-2D traps were analyzed in [19] and a number of analytic approximations for the condensate excitation spectrum were introduced.

Among a rich variety of nonlinear wave phenomena in dipolar BEC a particularly interesting one is the existence of stable isotropic and anisotropic 2D solitons [20, 21, 22], given the fact that multi-dimensional solitons in conventional BEC with contact atomic interactions are unstable against collapse or decay. In quasi-1D BEC with competing dipolar and local interactions bright soliton families were found in free space and in the presence of optical lattices [23]. In this work and subsequent papers [24, 25, 26] by analyzing different regimes of soliton collisions the regions in the parameter space where colliding solitons either merge or pass through each other were investigated. Dynamics and head-on collisions of solitons in dipolar BEC confined to 2D and 3D traps were investigated by variational approach and numerical simulations [27]. The transmission, reflection and splitting regimes for a dipolar soliton moving in a parabolic trap with a barrier potential at its center were also studied [28]. The existence, stability and collisions of dipolar bright solitons, showing a number of novel features, were analyzed [29]. The dipolar BEC with repulsive nonlinearity supports dark solitons, appearing as a traveling local reduction of the field amplitude. Distinctive features of dark solitons in dipolar BEC, such as interaction sensitive oscillations in external potentials, were reported [30].

During the last few years the existence, stability and dynamics of soliton bound states in dipolar BEC, so called soliton molecules, have attracted considerable interest. The possibility of molecular type of interaction between two solitons, reflected through the fact that they attract each-other at long distances and repel at short distances was theoretically studied [31, 32]. In these models individual solitons, forming the bound state, reside in separate layers. Numerical analysis of soliton bound states in quasi-2D dipolar BEC was presented [33]. The vibration spectrum of a two-soliton molecule in dipolar BEC confined to a quasi-1D trap was studied in [34], while the conditions when two solitons can form a stable bound state has been analyzed using the variational approach and numerical simulations [35]. Recently similar bound states of dark solitons in dipolar BEC were reported in [36], where the shape of the interaction potential between two dark solitons was numerically generated. The emission of phonons by colliding dark solitons and their simultaneous acceleration was another important finding in that work. The existence, stability and collision dynamics of dark solitons in quasi-1D dipolar BEC along with formation of bound states of dark solitons were investigated [37].

Nowadays research on dipolar quantum gases in general, and dipolar solitons in particular, is actively pursued due to novel physical phenomena observed in dipolar BEC [38, 39, 40] and many promising practical applications [41, 42]. For instance, the authors of a recent work [43] have demonstrated that oscillations of a dipolar condensate near the superconducting surface (atom chip) can induce eddy currents in it, which act back on the BEC. The interaction of a dipolar BEC with a superconductor via eddy currents provides a new physical mechanism for the coupling between these two macroscopic quantum systems. The eddy currents can change the oscillation frequency of the condensate moving near the superconducting surface and excite collective modes in BEC [44].

Despite the good progress achieved in understanding the physics of dipolar condensates and propagation of nonlinear waves in these media, there are many open issues awaiting their clarification. In particular, when the short-range contact interactions and long-range dipolar interactions between atoms nearly balance each-other an additional term in the equation, the so called beyond mean field term, becomes important [45, 46]. This term originating from quantum fluctuations is responsible for the existence of quantum droplets, a novel ultra dilute liquid state of matter [38, 39, 40]. The existence regimes, stability and dynamics of quantum droplets in dipolar BEC are among the actively explored topics in contemporary physics of ultracold gases.

Our objective in this work is to derive the analytic form of the interaction potential between two bright solitons in dipolar BEC and estimate the binding energy of a two-soliton molecule. Special attention will be given to investigation of weakly bound states of these localized structures with a small binding energy. Next we explore the dynamics of bright solitons and two-soliton molecules under slowly varying strength of dipolar interactions. We address the problem analytically using a variational approach (VA), and numerically by solving the governing nonlocal Gross-Pitaevskii equation (GPE). The resulting equations will allow to reveal the features of dipolar solitons and two-soliton molecules in the weakly bound regime, when the competing short-range contact interactions ($q<0$) and long-range dipolar interactions ($g>0$) nearly balance each-other. However, despite the smallness of the net attractive interactions, the above mentioned beyond mean field effects are supposed to be still negligible, and omitted from the governing equation.

2 Model equations and variational approach

We consider the dynamics of single solitons and two-soliton molecules in a dipolar BEC confined to a quasi-1D trap. The governing equation of our model is the one-dimensional nonlocal GPE for the mean-field wave function of the condensate $\psi(x,t)$

[TABLE]

where the subscripts denote partial derivatives. The dimensionless quantities, entering this equation are scaled using the frequency of the radial confinement $\omega_{\bot}$, atomic mass $m$ and radial harmonic oscillator length $\l_{\bot}=\sqrt{\hbar/m\omega_{\bot}}$ as follows: time $t\rightarrow t\omega_{\bot}$, space $x\rightarrow x/l_{\bot}$, wave function $\psi\rightarrow\sqrt{2a_{s0}}\psi$, the coefficient of contact interactions $q\rightarrow a_{s}/|a_{s0}|$, the coefficient of long-range dipolar interactions $g\rightarrow a_{d}/|a_{s0}|$. The background values of the atomic scattering length $a_{s0}$ and the length scale of dipolar interactions $a_{d}=md\,^{2}/\hbar^{2}$ are supposed to be known characteristics of the condensate. Although, in experiments both the atomic scattering length $a_{s}$ and dipole moment of atoms $d$ can be varied through external magnetic, electric or optical fields.

The response function $R(x)$ in the integral term of Eq. (1) characterizes the degree of non-locality of the medium. It shows how strongly the properties of the medium at a given location depend on the properties of its neighborhood. An analytic expression for the response function of a quasi-1D dipolar BEC was obtained in a single mode approximation [47].

[TABLE]

While relevant analytic calculations are complicated it is worth mentioning that the response function is characterized by a cusp at $x=0$. For practical reasons, however, another function was proposed in [23], which behaves smoothly at the origin and is more convenient for analytical treatment

[TABLE]

The value of the cutoff parameter $\sigma=\pi^{1/2}$ was found from the condition of equal areas beneath the curves of the single mode kernel and the one given by Eq. (3). A very good correspondence between the two response functions was noted (see Fig. 1 in Ref. [23]). Hence, our variational calculations will be based on the response function given by Eq. (3).

For an arbitrary external potential and nonlocal term Eq. (1) does not possess analytic solution. The variational approach [48, 49] is frequently used to study the dynamics of solitons governed by this type of equations. To develop the VA for our model we note that Eq. (1) can be generated from the following Lagrangian density

[TABLE]

The success of the VA strongly depends on the right choice of the trial function. On one hand the selected function should closely match the shape of the localized state (soliton or molecule), and on the other hand it should lead to analytically tractable integrals associated with this Lagrangian. Below we develop the suitable models for a single soliton and two-soliton molecule using the Gauss-Hermite trial functions.

2.1 Variational approach for single solitons

To derive variational equations for a single soliton moving in the external potential we use the Gaussian trial function

[TABLE]

where $A(t),a(t),b(t),\xi(t),v(t),\phi(t)$ are time-dependent variational parameters, meaning the amplitude, width, chirp parameter, position of the center-of-mass, velocity and phase of the soliton, respectively. The velocity is defined as a time derivative of the soliton’s center-of-mass position $v=\xi_{t}$. The norm $N=A^{2}a\sqrt{\pi}$, which is a conserved quantity of the governing equation, is proportional to the number of atoms in the condensate.

Substitution of Eq. (3) and Eq. (5) in the Lagrangian density (4), followed by spacial integration, $L=\int_{-\infty}^{\infty}{\cal L}dx$, leads to the effective Lagrangian

[TABLE]

where ${\cal U}({\rm a},{\rm b},{\rm c})=\frac{1}{\Gamma({\rm a})}\int_{0}^{\infty}e^{-{\rm c}t}t^{{\rm a}-1}(t+1)^{{\rm b}-{\rm a}-1}dt$ is the confluent hypergeometric function [50], with $\Gamma({\rm a})$ being the gamma function. Now using the associated Euler-Lagrange equations $d/dt(\partial L/\partial\nu_{t})-\partial L/\partial\nu=0$ for the variational parameters $\nu~{}\rightarrow~{}a,~{}b,~{}\xi,~{}\phi$ one can derive the following equation for the width

[TABLE]

where the non-locality function $F_{s}(a,\sigma)$ is given by

[TABLE]

The stationary width $a_{0}$ and amplitude $A_{0}=\left[N/(a_{0}\sqrt{\pi})\right]^{1/2}$ of the soliton with a given norm $N$ can be found from the fixed point of Eq. (7).

The Eq. (7) has a formal analogy with the equation of motion for a unit mass particle moving in the anharmonic potential

[TABLE]

The presence of a local minimum in this potential for a given set of parameters indicates the existence of a self-bound localized state, which is also called soliton, despite the non-integrability of the governing Eq. (1). In Fig. 1 we illustrate the numerical experiment with Eq. (1) when the VA predicted wave profile Eq. (5) is introduced as initial condition, and the strength of the dipolar interactions is linearly decreased with time. Since the variational solution is approximate one, at the beginning the amplitude while decreasing oscillates at a high frequency, which is an evidence for a strongly bound soliton. At some critical value of the dipolar strength $g=g_{cr}$ the soliton disintegrates, and the local minimum in the potential Eq. (9) disappears. Below this critical value of the strength of dipolar interaction the net nonlinearity in Eq. (1) becomes repulsive, and the system does not support self-bound bright localized states. Throughout the paper we consider the repulsive contact interactions with $q=-1$.

2.2 Variational approach for two-soliton molecules

As a trial function for two-soliton molecules we select the Gauss-Hermite function proposed in [51, 52]

[TABLE]

where the meaning of the variational parameters are similar to those of a single soliton considered above. It should be noted that Eq. (10) adequately describes the properties of two-soliton molecules when the solitons move near their equilibrium positions. At large separations of solitons this waveform deviates from the superposition of two anti-phase solitons, leading to a lower accuracy of the method.

Substitution of Eq. (10) into Eq. (4) and integration over space variable gives the effective Lagrangian

[TABLE]

The equation for the width can be obtained from the Euler-Lagrange equations

[TABLE]

where

[TABLE]

$K_{0}(x),K_{1}(x)$ - are the modified Bessel functions of the second kind. The fixed point of Eq. (12) gives the stationary shape of the two-soliton molecule. The waveforms of a two-soliton molecule according to these equations are presented in Fig. 2 for two different strengths of dipolar interactions. The center-of-mass positions of the solitons are located at $x_{0}=\pm\,2a_{0}/\sqrt{\pi}$, where $a_{0}$ is the fixed point of Eq. (12).

The properties of two-soliton molecules in the weakly bound regime can be probed by slowly reducing the strength of the dipolar atomic interactions. In the experiments tuning the dipolar interactions can be achieved by fast rotation of the orientation of the dipoles in a magnetic field [53, 54]. When the competing long-range dipolar and short-range contact interactions nearly balance each-other, the bound state of two solitons becomes unstable against small perturbations. A slight imbalance may lead to disintegration of the molecule, as shown on the right panel of Fig. 2. In numerical simulations we induce mild vibrations of the molecule by chirping the initial waveform. The oscillatory component in the dynamics helps to detect precisely the time instance when the dissociation of the molecule occurs. Just before the disintegration the soliton molecule has a very small binding energy. Hence, without the oscillatory component in the dynamics, it would be difficult to detect the time instance, when the molecule dissociates into freely moving solitons.

3 Potential of interaction between solitons

From numerical simulations of the GPE (1) it was noticed, that the potential of interaction between solitons at large separations ($a~{}\gg~{}a_{0}$) notably deviate from the prediction of variational approach (see the right panel of Fig. 3). Obviously, this is due to the fact that the ansatz (10) does not well represent two separate solitons when the distance between them is large.

The potential of interaction between two solitons, shown on the right panel of Fig. 3 by blue symbols, has been constructed as follows: Two sufficiently far separated solitons (but still “feeling” each - other) are prepared in anti-phase configuration as initial condition in Eq. (1). Due to the long - range dipolar attractive interaction they start to move toward each - other. When the inter-soliton distance becomes small, they repel due to the anti-phase configuration. The maximum and minimum of the center of mass positions correspond to turning points $x_{p_{1,2}}$ in the potential well. Each time when the soliton passes its equilibrium position ($x=\pm\,2a_{0}/\sqrt{\pi}$), the potential energy fully transforms into kinetic energy. By measuring the maximal velocity $v_{m}$ one can define the potential energy at the turning points $U(x_{p_{1,2}})=v_{m}^{2}/2$, where the mass of the soliton is set to one. Repeating this procedure with different initial separations, one can construct the shape of the potential.

An alternative approach to evaluation of the interaction potential between two solitons consists in analyzing the equation for the center-of-mass $x_{c}$ of one soliton moving in the potential created by the other soliton [55]. Suppose we have “left” and “right” solitons positioned at $x_{c}$ and $-x_{c}$ respectively. In the case of attractive interactions the left soliton moves towards the right soliton due to its modification of the “refractive index” of the medium

[TABLE]

Using the response function (3) and ansatz (5) we obtain the equation of motion for the center-of-mass coordinate of the left soliton

[TABLE]

On the right hand side of this equation we have the force, acting on the left soliton. The potential is the integral of the force with respect to $x_{c}$ (below we write it by dropping the subscript)

[TABLE]

where ${\rm erfc}(z)$ is the complimentary error function[50]. The first term with $q<0$ is responsible for the short-range repulsive interaction, while the second term with $g>0$ corresponds to the long-range attractive interaction. The potential Eq. (16) is depicted in Fig. 4.

By comparing this figure with the right panel of Fig. 3 one can observe the qualitative agreement between the numerically constructed potential and the prediction of Eq. (16).

4 The binding energy of a two-soliton molecule

The total energy of a single soliton or two-soliton molecule can be evaluated from the GPE energy functional

[TABLE]

where $\psi(x)$ is the waveform of the localized state, while $R(x)$ is the response function.

By substituting in Eq. (17) the waveforms, corresponding to the single soliton Eq. (5) and two-soliton molecule Eq. (10) we get the corresponding energies

[TABLE]

where $N_{s}$, $a_{s}$ - are the norm and width of a single soliton, respectively, while the function $F_{s}(a_{s},\sigma)$ is given by Eq. (8).

Similar expression for the total energy of the molecule has the form

[TABLE]

where $N_{m}$, $a_{m}$ - are the norm and width of a two-soliton molecule, respectively, while the function $F_{m}(a_{m},\sigma)$ is given by Eq. (13).

The binding energy of a two-soliton molecule can be estimated as the difference between the energy of the molecule, and the energies of two solitons well separated from each other, i.e. non-interacting solitons.

[TABLE]

To simplify the resulting expression, we use the approximate relations between the parameters of the molecule and single solitons which form this molecule being placed at the equilibrium distance.

The waveform of a two-soliton molecule (10) can be well approximated by two anti-phase solitons, as illustrated in Fig. 5

[TABLE]

with the following amplitude, width and half separation [35]

[TABLE]

Due to the destructive interference the norm of the molecule will be slightly less than twice of the soliton’s norm. We shall neglect this effect and set $N_{s}=N_{m}/2$. These relations lead to the expression for the binding energy of the molecule

[TABLE]

where $N_{m}$, $a_{m}$ are the norm and equilibrium width of the molecule, given by the fixed point of Eq. (12).

5 Conclusions

The analytic expressions for the potential of interaction between bright solitons and the binding energy of a two soliton molecule in quasi-1D dipolar BEC have been derived. The mathematical model, based on the variational approach, has been corroborated with the results of numerical simulations of the nonlocal Gross-Pitaevskii equation. The developed model allows to explore the properties of solitons and two-soliton molecules in the weakly bound regime, when the short-range contact atomic interactions and long-range dipole-dipole interactions nearly balance each-other. The obtained results can be of interest in future experiments with dipolar BEC, including the conditions close to the formation of quantum droplets.

Acknowledgments

This work has been supported by the KFUPM research group projects RG1503-1 and RG1503-2. B.B.B and U.A.K thank the Physics Department at KFUPM and Saudi Center for Theoretical Physics for the hospitality provided during their visits.