Interplay between binary and three-body interactions and enhancement of stability in trapless dipolar Bose-Einstein condensates

We investigate the nonlocal Gross-Pitaevskii (GP) equation with long-range dipole-dipole and contact interactions (including binary and three-body collisions). We address the impact of the three-body interaction on stabilizing trapless dipolar Bose-Einstein condensates (BECs). It is found that the dipolar BECs exhibit stability not only for the usual combination of attractive binary and repulsive three-body interactions, but also for the case when these terms have opposite signs. The trapless stability of the dipolar BECs may be further enhanced by time-periodic modulation of the three-body interaction imposed by means of Feshbach resonance. The results are produced analytically using the variational approach and confirmed by numerical simulations.

The theoretical description of a dilute weakly interacting DBEC is based on the Gross-Pitaevskii (GP) equation with the nonlocal DD-interaction term [1,2,23,24,25]. In particular, the combination of local and nonlocal terms in the GP equation can support various species of bright and dark matter-wave solitons. In the alkali BECs, bright solitons exist when the negative (attractive) contact interaction exactly balances the dispersion (kinetic-energy) term [26,27]. In DBEC, the nonlocal DD interaction term may reinforce local ones originating from the s-wave contact interaction. DD interactions have caught the limit with the atomic condensates of lanthanide series elements such as 168 Er and 164 Dy. Long-range interactions also play a crucially important role in optical of nonlocal media, where they induce modulational instability, solitons, and vortices [28,29]. The DD interactions in BECs can support bright solitons even for the positive (repulsive) contact interaction [23].
Following the scheme of the stabilization of the inverted (Kapitza) pendulum [30], scenarios for stabilization of two-dimensional (2D) optical [31] and matter-waves [32,33,34] by means of the "nonlinearity management" [38], i.e., the cubic nonlinearity periodically switching between self-attraction and repulsion, have been elaborated. This concept has been subsequently applied to 3D vortex solitons [32,34] and extended to the model containing the three-body interaction [35] and quantum fluctuations [36]. The stabilization of self-repulsive BECs with periodically varying higher-order interactions has been explored recently [37]. The investigation of spatial and temporally varying nonlinearities has also drawn considerable interest in optics [39,31,40,41] and other BEC settings [42,43,44,45,46,47,48,49,50] The above investigations addressed systems with local interactions. Recently, the possibility of the stabilization of DBECs by temporal modulation of the contact interaction was demonstrated too [51]. In the latter case, a potential minima necessary for self-trapping is found to be absent without the time-periodic modulation of the two-and three-body interaction. The current work primarily addresses this situation, attempting to stabilize DBECs by adding three-body terms to the binary interaction, cf. Refs. [52,53,54,55,56,57]. As a result, the present work offers a simple protocol to stabilize multidimensional trapless DBECs.
The organization of the paper is as follows. In Section 2, we set the mean-field model based on the nonlocal GP equation. Then, we introduce the variational method in Section 3 and demonstrate the stabilization of a trapless DBEC without time-periodic nodulation of the contact interaction. Numerical results produced with the help of Runge-Kutta (RK4) method are also presented in Section 3. In Section 4, we report numerical results for the 3D time-dependent GP equation simulated by means of the Crank-Nicholson method. Concluding remarks are presented in Section 5.

The model
At ultra-low temperatures, a DBEC with two-body, three-body and nonlocal DD interactions can be described by the following time-dependent GP equation [1,2,23,24,25,58]: where φ(r, t) refers to the condensate wave function, N the number of particles with magnetic dipole moment µ, m the mass, ∇ 2 the Laplacian operator, the reduced Planck constant, V(r) = 1 2 m ω 2 x x 2 + ω 2 y y 2 + ω 2 z z 2 is the trapping potential with frequencies ω x , ω y and ω z , and a(t) is the time-modulated atomic scattering length. The potential of the dipolar interaction for magnetic dipoles is U dd (R) = µ 0 µ 2 4π 1−3 cos 2 θ |R| 3 , where R = r − r ′ determines the relative position of dipoles and θ is the angle between R and the vertical (z) polarization direction, µ 0 is the free-space permeability [23]. The parameter η(t) is the strength of the threebody interaction which depends on both the s-wave scattering length a(t) and the universal constants d 1 , d 2 , a 0 and s 0 [58]. This parameter reads η(t) = 12π 2 a(t) 2 where the numerical values of the universal constants are given in Refs. [58,59]. However, the strength of the three-body interaction is very small when compared with strength of the two-body interaction [60]. In some other regimes, the three-body interaction may become a dominant factor under special conditions as in the cases of higher densities and Effimov resonance [58]. The normalization is dr|φ(r, t)| 2 = 1. ( To compare the dipolar and contact interactions, it is often useful to introduce the length scale a dd ≡ µ 0 µ 2 m/(12π 2 ), whose experimental values for 52 Cr, 164 Er and 168 Dy condensates are 16a 0 , 66a 0 and 131a 0 , respectively, a 0 being the Bohr radius [1,2]. In the present investigation, we use the parameters of 52 Cr. It is convenient to rescale Eq. (1) into a dimensionless form. For this purpose, we apply the transformationr = r/l,R = R/l,ā(t) = a(t)/l,ā dd = a dd /l,t = tω,x = x/l,ȳ = y/l,z = z/l,φ = l 3/2 φ, where the harmonic-oscillator length is l = √ /(mω)and rewrite Eq. (1) (dropping the overbar) as where For an axially-symmetric (with ν = γ in Eq. (4)) disk-shaped DBEC with a strong axial trap (λ > ν, γ), we presume that the structure of DBEC in the axial direction amounts to the axial ground state, where ρ ≡ (x, y), ψ(ρ, t) is the effective 2D wave function and d z is the axial harmonic-oscillator length. To obtain the effective 2D equation for the disk-shaped DBEC, we substitute expression (5) in Eq. (3), multiply it by φ axial (z) and integrate the result over z to arrive at the 2D equation [23,24,25] i where we have introduced the time-dependent factor d(t), which may provide temporal modulation of the effective 2D trapping potential defined as Next, the external trap may be removed adiabatically. For that purpose, the parameter d(t) is slowly ramped down from 1 to 0. The dipolar term is written in the Fourier space after calculating the convolution of the corresponding variables [23] In Eq. (6), the length is measured in units of the characteristic harmonic oscillator length l = √ /mω, time t in units of ω −1 and energy in units of ω.
The quasi-2D and full 3D nonlocal GP equations are then simulated using the split-step Fourier method. In our dimensionless units, we use the space and time step sizes 0.05 and 0.001.

The variational method
To gain an analytical insight into the condensate dynamics, we use the variational approximation (VA) with the Gaussian ansatz as a trial wave function for the solution of Eq.(6) in the absence of an external trap [23]: Here the variables R(t) and β(t) are real-valued and stand for the radius and radial chirp of the self-trapped state while the amplitude is determined by the normalization condition (2). The Lagrangian density generating Eq. (6) in the absence of the trap (d(t) = 0) is given by where ψ * denotes the complex conjugate of the wave function ψ. The trial wave function (8) is substituted into the Lagrangian density (9) and the effective Lagrangian is calculated by integrating the Lagrangian density. Then, the following Euler-Lagrangian equations for the variational parameters can be derived from the effective Lagrangian as By associating Eqs. (10) and (11), we obtain the following equation for the evolution of the width R(t) as We plan to consider the effects of periodic modulation of the contact interaction in the form of on the stability of DBEC where ǫ 0 , χ 0 are constant parts of the two-body and three-body interactions respectively and ǫ 1 , χ 1 are the amplitudes of their temporal modulation. The Kapitza averaging method can be employed to treat the oscillatory parts [30]. With the explicitly included oscillating nonlinearity, we get the following equation for the evolution of the radial width (cf. Ref. [23]): Now, R(t) may be separated into a slowly varying part R 0 (t) and a rapidly varying one R 1 (t), Keeping the terms of the order of Ω −2 in R 1 (t), one may obtain the following equations of motion for R 0 (t) and R 1 (t), where the < R 1 sin(Ωt) > indicates the time average of the rapid oscillation. From Eq.(15), we obtain We now substitute R 1 (t) into Eq.(16) to obtain the following equation of motion for R 0 (t) The VA suggests that the effect of the DD interaction is to minimize the impact of the constant part of the contact interaction for a dd > 0. Thus, one can conclude that the system may becomes 5 effectively attractive depending on the choice of a dd and ǫ 0 . Hence one can explore the possibility of creation of bright solitons even for positive (repulsive) scattering lengths provided a dd > ǫ 0 . Thus, it is pretty obvious that DD interaction minimizes the effect of contact interaction as dipolar BECs really need non-zero dipolar strength to compete with contact interaction. In addition, DD interaction (a d d > 0) always counteracts with two-body contact interaction. We do agree under the present circumstances that the dipoles assume a pancake that are oriented in the z-direction and dipolar BEC expands along the xy plane. To arrest the expansion, here we use the attractive two-body contact interaction. The effective potential U(R 0 ) corresponding to the above equation of motion can be written as where δ(ξ 0 ) = [1 + 2ξ 2 0 − 3ξ 2 0 d(ξ 0 )]/(1 − ξ 2 0 ) and ξ 0 = R 0 /d z . Now, one can scrutinize the nature of the effective potential for different system parameters, namely the strengths of the two-body, three-body and DD interactions. In the present analysis, we explore the possibility of stabilizing the trapless dipolar condensates by manipulating the constant and oscillatory part of the two-and three-body interactions.
In Fig. 1, we display the effective potential (19) of the 52 Cr condensate in the absence of the temporal modulation of both nonlinearities. One may infer from Fig. 1(a) that the potential energy curves do not show any minimum for both repulsive (dotted and solid curves) and attractive (the dashed curvraction (ǫ 0 ) strengths. In the repulsive case, this indicates that the condensate leads to expansion due to the repulsive two-body interaction and the kinetic pressure. On the other hand, no stable counter balance point (i.e., potential minimum) exists for the attractive case either, which indicates at the attractive constant two-body interaction overcomes the kinetic pressure, leading to the collapse of the condensates. These conclusions are commonly known for the cubic GP/nonlinear-Schrödinger equation in 2D [61] Next, in Fig. 1(b), we show how the interplay between the two-and three-body collisions 6 along with the DD interaction can lead to a stable condensate. We first get back to one of the potential curves in Fig. 1(a) for the attractive case corresponding to ǫ 0 = −25a 0 , ǫ 1 = 0, χ 0 = 0, χ 1 = 0. In this case, the absence of a repulsive force to balance the self-attraction force arising due to the constant part of the two-body and DD interactions initiates the collapse of the condensates. Now, the combined impact of the repulsive force arising due to the repulsive three-body interaction safely counterbalances the self-attraction as is evidenced by the presence of the potential minima in the respective curve in Fig. 1(b) while the nonlinearity management is not included. Further, in Fig 2, we show the stability domains in the plane of the constant parts of the attractive two-body interaction in competition with repulsive three-body interactions. The diagrams show stable regimes between the collapse and expansion regions. From figs. 1-2, we conclude that with help of the repulsive three-body interaction, one can stabilize the trapless DBEC with attractive two-body interaction. This has been done without the contribution of the oscillatory part of the contact interactions. Next, we plan to stabilize the trapless DBECs with repulsive two-body contact interaction. We first think about one of the potential energy curves from Fig. 1(a) for the repulsive case corresponding to ǫ 0 = +25a 0 , ǫ 1 = 0, χ 0 = 0, χ 1 = 0. Under these circumstances, the additional repulsive force caused by the constant part the three-body interaction along with the constant part the two-body interaction and the DD interaction makes the condensate decay through the expansion as shown in Fig. 3(a) by the dotted red curve. Now to stabilize the system in this case, we introduce an attractive force arising from the constant part of the three-body interaction, χ 0 < 0. We observe the presence of an inverted dashed blue-colored potential curve in Fig. 3(a) implying that we have not yet reached a stable regime. Moreover, a stable region has been obtained neither with the help of an attractive nor with that of a repulsive three-body interaction. However, the introduction of the time-dependent part of the two-body interaction along with the attractive three-body interaction recovers potential energy minima represented by the black curve in Fig. 3(a). Here, the two-body repulsion is balanced by the interplay between the attractive three-body interaction and the time-modulated part of the two-body interaction, ǫ 1 . Figure 3(b) depicts the time evolution of the width of the condensate corresponding to Fig. 3(a). Figure 3(b) clearly shows the dynamical regimes corresponding to the expansion, collapse and stability of DBEC. In particular, the stable case with χ 0 = −0.025 and ǫ 1 = 4ǫ 0 displays robust small-amplitude oscillations of the width. On the other hand, the expansion and collapse are observed, respectively, at χ 0 = 0.025, ǫ 1 = 0, and χ 0 = −0.025, ǫ 1 = 0. Further, in Figs. 4(a) and (b), we show the effective potential for different values of the constant part of the threebody interaction and the corresponding VA-predicted evolution of the width of the condensate. Once again, robust oscillations of the widths confirm the stability of the trapless DBECs which 8 corresponds to the presence of the minimum in the potential curve. We further show that the introduction of the time-dependent part of the three-body interaction, χ 1 may further enhance the stability of the trapless dipolar repulsive BECs. Figures 5  (a) and (b) display the effective potential for various values of χ 1 and the evolution of the width of the condensate respectively. Figure 5(a) confirms that one can make the potential well deeper to stabilize the condensate. In Fig.4(a), the potential minimum starts to appear for χ 0 = −0.020ǫ 0 which is plotted as dash-dotted blue curve. However, no potential minimum exists for χ 0 = −0.010ǫ 0 which is plotted as dotted black curve. From the above, we conclude that the potential depth starts to appear for χ 0 ≥ | − 0.020ǫ 0 -. Also, if we increase χ 0 to −0.025ǫ 0 , and −0.030ǫ 0 , the potential depth increases further. From this, it is obvious that if the attraction becomes too strong, the system becomes highly unstable leading to its collapse. But, the inclusion of the χ 1 along with the case discussed in Fig.4(a) reduces the strong attraction. If the attraction is weak for its size, the system becomes weakly attractive in the final stage and it expands to infinity. On the other hand, if we include a suitable χ 1 , then the attraction due to χ 0 is balanced by the oscillation due to the effect of χ 1 and the system becomes stable as shown in Fig.5. In Table 1

Three-dimensional Numerical Results
In this section, we consider the full 3D model for direct numerical simulations of Eq. (20). First, we find the ground-state solution by solving Eq. (20) by means of the imaginary-time (t → −it) propagation in the presence of the external trap and two-body contact plus dipolar interactions in the absence of three-body interactions. Then, the evolution of the input provided by the ground state is simulated in the real time. In the course of the real-time evolution, the external trap is adiabatically removed and the three-body interaction is ramped up which brings the system to the same form for which the VA was elaborated above. We now consider the following 3D GP equation and numerically solve by gradually decreasing d(t) and increasing g(t) , χ(t). The DD interaction was evaluated by means of the fast Fourier transform [8]. Typical space and time steps for the numerical grid were taken as 0.05 and 0.005, respectively. In the course of time evolution, we increase the nonlinearity coefficients from 0 at each time step as with f (t) = t/τ for 0 ≤ t ≤ τ and f (t) = 1 for t > τ. At the same time, the trap is removed by changing d(t) from 1 to 0 by d(t) = 1 − f (t). Thus, the trap is eventually removed and the final values of the coefficients are attained at t = τ and the periodically oscillating nonlinearities become g(t) = g f [a 1 − b 1 sin (Ωt)] and χ(t) = χ f [a 2 − b 2 sin (Ωt)] is introduced for t > τ, cf. Refs. [35,34,33]. After the complete removal of external trap, the condensates will collapse for 11 strong attraction due to the final nonlinearity. If the attraction due to the final value of nonlinearity after switching off the external trap is too weak, the system becomes weakly attractive and the condensates start to expand. In Fig. 6, we present the results for the condensate of 52 Cr atoms which have a moderate dipole moment corresponding to a dd = 16a 0 [1,2]. The figure demonstrates how one can stabilize the condensate with ǫ 0 = −25a 0 by varying the strengths of the repulsive three-body interactions (without the contribution of the oscillatory part of the interactions). In panel (1,1) the condensates are found to be stable up to 2ms for χ 0 = -0.013ǫ 0 . Similarly, in panels (1,2), (2,1) and (2,2), the condensates are stable up to 2.5ms, 3ms and 5ms for χ 0 = -0.012ǫ 0 , -0.011ǫ 0 and -0.01085ǫ 0 , respectively. Further evolution of the condensates may perhaps lead to their collapse or expansion.

Conclusion
We have considered the dynamics of BEC with parameters corresponding to 52 Cr with the motivation of studying the impact of the three-body interaction on the stabilization of trapless dipolar BECs with two-body contact interactions. By means of the variational approximation, we have produced phase diagrams showing stable states of the trapless dipolar BECs which indicates the enhancement of the stabilization of the condensates by the reinforcement of binary interaction with three-body collisions. The stability is indicated by the presence of the minimum of the effective potential and stable oscillation of the widths in the course of the evolution. Further, we have performed full 3D simulations to confirm the stability. In particular, the stability may be enhanced by the time-periodic modulation of the three-body interaction added to the attractive binary interaction.
It may be relavant to extend the present analysis to 3D self-trapped condensates with embedded vorticity. In that case, stability of the solitary vortices against spontaneous splitting is a crucially important issue [62].