Holding and transferring matter-wave solitons against gravity by spin-orbit-coupling tweezers

We consider possibilities to grasp and drag one-dimensional solitons in two-component Bose- Einstein condensates (BECs), under the action of gravity, by tweezers induced by spatially confined spin-orbit (SO) coupling applied to the BEC, with the help of focused laser illumination. Solitons of two types are considered, semi-dipoles and mixed modes. We find critical values of the gravity force, up to which the solitons may be held or transferred by the tweezers. The dependence of the critical force on the magnitude and spatial extension of the localized SO interaction, as well as on the solitons norm and speed (in the transfer regime), are systematically studied by means of numerical methods, and analytically with the help of a quasi-particle approximation for the soliton. In particular, a noteworthy finding is that the critical gravity force increases with the increase of the transfer speed (i.e., moving solitons are more robust than quiescent ones). Nonstationary regimes are addressed too, by considering abrupt application of gravity to solitons created in the weightless setting. In that case, solitons feature damped shuttle motion, provided that the gravity force does not exceed a dynamical critical value, which is smaller than its static counterpart. The results may help to design gravimeters based on ultracold atoms.


I. INTRODUCTION
Optical tweezers, which were first constructed in the classic works of Ashkin et al. [1,2] for selecting and driving microscopic dielectric particles by means of laser beams, are nowadays used in various experiments aimed at trapping and transfer of atomic Bose-Einstein condensates (BECs) [3]- [6]. Among these settings, is holding BEC under the action of gravity. The interaction of ultracold condensates with the gravitational field is a problem with applications to experiments in gravitational physics and development of precise measurement techniques. In particular, atomic Raman interferometry [7]- [11] and Bloch oscillations [12]- [14] of ultracold atoms trapped in a vertical optical lattice were used for precise measurement of the free-fall acceleration, as well as its gradient, aimed at geological applications. A combination of these techniques was reported to remarkably improve the accuracy of atomic interferometry [15] and fundamental measurements of gravity [16]. Recently, a fountain gravimeter, which achieved a sub-µG degree of accuracy, was fabricated on an atom chip [17], and a compact quantum gravimeter was realized in polarizationsynthesized spin-dependent optical lattices [18]- [20]. Dynamics of BEC under the action of microgravity has also been a subject of many experiments [21]- [22].
The self-attractive intrinsic nonlinearity of BEC, induced by inter-atomic collisions, which is accurately modeled by the Gross-Pitaevskii equation (GPE) [23], gives rise to self-trapped states, in the form of solitons [24] and breathers [25]. As concerns the above-mentioned topics, relevant aspects are manipulations of matter-wave solitons with the help of tweezers [5], motion of solitons [26]- [31] and breathers [32,33] driven by gravity, and Bloch oscillations of gap solitons in optical lattices under the action of a constant (gravity) force [34]. In this context, it is worthy to mention great diversity in the design of matter-wave interferometers offered by matter-wave solitons [35]- [45].
However, the applications of optical tools to BEC are hampered by effects of heating the condensate by laser beams which are used for the creation of tweezers [46]- [48] and optical lattices [49]- [55]. In this work, we aim to elaborate a tweezers scheme for matter-wave solitons based on spatially localized spin-orbit (SO) coupling applied to BEC, by means of properly focused laser illumination. One of incentives for the introduction of the scheme is the fact that the SO interaction may be induced in the condensate without heating it [56]- [59].
The configuration is schematically shown in Fig. 1, which assumes that the soliton is placed in a narrow vertical cigar-shaped potential trap and illuminated by a horizontal Raman laser beam inducing the SO coupling and focused on a segment of the vertical "cigar". In fact, the focusing does not need to be very tight. Indeed, for a typical longitudinal extension of solitons, which is, normally, several microns [60], the Raman-beam illumination should be confined on a comparable segment, which is, thus, much larger than the diffraction limit for the focusing, ≃ λ/2, as determined by the characteristic beam's wavelength, λ ≃ 1 µm [61,62].
FIG. 1: (Color online) The sketch of the setting: a quasi-1D matter-wave soliton, designated by the dark red color, is created in the vertical cigar-shaped potential trap, while horizontal beams, which induce the confined SO coupling, hold the soliton at a particular vertical position. G represents the gravity force.
The SO coupling in BECs, which emulates the eponymous effect in physics of semiconductors, has been an active research subject in recent years [61]- [67]. In the combination with the usual collision-induced nonlinearity, the SO coupling was theoretically shown to support various species of one-dimensional (1D) solitons [68]- [77]. A surprising prediction is that the interplay of the linear SO interaction and cubic self-attraction leads to stabilization of multidimensional solitons in free space, which are completely unstable in the absence of the SO terms [78]- [87].
As concerns the application of the SO coupling to a spatially confined region, such as in the configuration sketched in Fig. 1, it was predicted that 1D solitons and their bound states may be maintained by means of localized SO interaction [88,89]. Further, it was recently demonstrated that spatially confined SO coupling can be used to maintain stable 2D solitons [90]. A related possibility is the use of effectively 1D and 2D (low-dimensional) SO coupling for the stabilization of, respectively, 2D and 3D solitons (in the full dimension), including a vortex component [91].
The subject of the present work is the stability and dynamics of quasi-1D solitons held or transported by localized SO coupling against the action of gravity, as suggested by Fig. 1. The consideration is based on a combination of systematic numerical simulations and the use of an analytical approximation. The results may help to develop a new type of cold-atom gravimeter and, more generally, new techniques for precise measurements.
The remainder of the paper is structured as follows. The setting is introduced in Section II, which is followed by presentation of numerical and analytical results for stationary solitons, obtained in both quiescent and moving reference frames (the latter one pertaining to the transport of solitons), in Section III. Motion (shuttle oscillations) of solitons is numerically studied in Section IV, and the paper is concluded by Section V.

II. THE MODEL
The setting sketched in Fig. 1, with the spatially confined SO coupling of the Rashba type [92], which can be created by the laser beam [93], is represented by the following system of coupled GPEs for a pseudo-spinor (two-component) mean-field wave function, (ψ + , ψ − ), written in the scaled form: where the strength of the attractive self-interaction in each component is normalized to be 1, γ is the relative strength of the cross-interaction between the components, and G is the gravity force. In the experiment, values of G may be controlled by means of the angle between the direction of the quasi-1D trap holding the condensate and the vertical direction, as well as by running the experiment in microgravity settings [21,22]. Because the natural transverse structure of laser beams is Gaussian, we adopt the corresponding shape of the spatial localization of the SO-coupling coefficient, cf. Refs. [88] and [90]: where L defines the confinement size, while amplitude Λ 0 may be fixed to 1 by means of additional rescaling. Nevertheless, it is convenient to keep Λ 0 as a free parameter, to explore dependence of the predicted effects on the strength of the SO interaction. Terms ∼ dΛ/dx in Eq. (1) appear when the GPE system is derived from the respective Hamiltonian, Solutions are characterized by their total norm, which is a dynamical invariant of the system, Moving tweezers, which can transfer trapped solitons, are defined as profiles (2) moving at velocity c, In particular, we aim to find a critical value of velocity, up to which the trapped solitons may be stably transferred. This is a nontrivial issue because Eq. (1) is not Galilean invariant. To address it, we follow Ref. [78] and rewrite Eqs. (1) in the moving reference frame, with coordinate x ′ = x − ct and transformed wave function, Thus, the underlying equations indeed change in the moving frame. Note that another transformation, ψ ′ ± → ψ ′ ± * , t → −t, casts Eq. The linear-mixing terms ∼ icΛ, which are generated by the transformation in Eq. (6), strongly affect the structure of solitons, but, generally, do not make them unstable, as shown below.

III. STATIONARY SOLUTIONS
A. Stationary states in the quiescent reference frame

Numerical results
In the 2D system, two types of fundamental solitons exist under the action of the SO coupling. One is built as a composite state including a fundamental soliton in component ψ + (one with zero topological charge, S + = 0), and a solitary vortex (with charge S − = 1) in component ψ − . States of this type are classified as semivortex (SV) solitons [78] (similar composite modes, found in a model with repulsive nonlinearity, were called "half vortices" [94]). Due to the symmetry of the SO coupling, SVs with the opposite chirality, defined by the vorticity set (S + = −1, S − = 0), exist too. The other type of 2D fundamental solitons is the mixed-mode (MM) soliton, which is built as a superpositions of states with vorticities (0, 1) and (−1, 0) in the two components.
A 1D counterpart of the SV soliton is a semidipole (SD) composite state, which is composed of spatially even and odd modes in its two components, and can be found numerically starting from the following input: with width W and real amplitudes A 1,2 . A 1D version of MM solitons is built as a superposition of even and odd modes in both components [95]. Accordingly, the numerical search for such solutions may be initiated by input The form of these is suggested by their 2D counterparts used in Ref. [78] for constructing solitons of the SV and MM types.
In the numerical analysis, we chiefly focus on the Manakov's nonlinearity, with γ = 1 in Eq. (1) [96] (although other values of γ are considered too). This case is relevant for experimental realization of the SO coupling in binary mixtures of two different hyperfine atomic states [62]- [67], for which scattering lengths of collisions between atoms belonging to the same or different atomic states are nearly equal. The interplay of the SO coupling and nonlinearity of the Manakov's type in the 2D system features special properties. In particular, the SD and MM solitons are included, in this case, in a broader family of 2D solitons with equal energies and equal chemical potentials for a fixed soliton's norm (4) [78]. A systematic numerical analysis demonstrates below that states maintained by the SO trap (2) are especially robust against the effect of gravity in the system with the Manakov's nonlinearity.
Stationary solutions of Eq. (1) for solitons can be obtained as follows.
(i) Use the imaginary-time method (ITM) [97,98] to produce a stationary solution at G = 0. Because the SD and MM solitons represent mutually degenerate ground states in the Manakov's case, either of them can be generated by the ITM with inputs chosen as per Eqs. (7) and (8), respectively (at γ < 1 and γ > 1, the ground state is represented by solitons of the SD and MM types, respectively [78]). Due to the SO coupling between the component of the soliton, they share a common value of the chemical potential µ, taking the form of The chemical potential of the ITM-produced solutions can be accurately identified by the substitution of the solutions in Eq. (1) and looking at the central point (x = 0). (ii) Use the wave function and chemical potential of the solution with G = 0 as an initial guess, and then apply the squared-operator method (SOM) [99], to produce a stationary solution with G (1) = δG, where δG is a small value.
(iii) Using SOM, continue to build solutions, increasing the gravity force by small steps, In the case of G > 0, the formal ground state (the one with the minimal energy) corresponds to the free fall of the wave packet towards x → ∞. For this reason, ITM, which seeks for the system's ground state, is irrelevant for constructing solitons at G > 0, and it was necessary to use SOM at stages (ii) and (iii). The results show that the shape of MM solitons quickly tends to become close to the SD type with the increase of G, therefore the consideration is focused below, chiefly, on the SD modes.
Results of the numerical analysis are summarized in Fig. 3. First, panel (a) displays a stability area for the stationary soliton solutions in the (N, G) plane. They exist and are stable in the yellow area, between G = 0 and the critical (largest) value G = G cr . Horizontal lines connecting points G = 0 and G = G cr show that the variation of G at fixed values of the chemical potential, µ, practically does not affect the respective value of the soliton's norm, N . At G > G cr , solutions with N < 3 become delocalized, developing nonvanishing tails which extend to the boundary of the integration domain, see Fig. 4(a). This result implies that, under the action of strong gravity, the condensate leaks from the finite region at which the SO coupling is applied. Further, in the interval of 3 < N < 5, at G > G cr the numerical procedure yields only zero solution, while unstable solitons are obtained for N > 5 and G − G cr positive but not too large. As shown in Fig. 4(b), unstable solitons eventually escape from the trap and start the free fall towards x → +∞. Only zero solution is produced by the numerical method at any value of N for G − G cr large enough. In the case of γ = 1 (the Manakov's nonlinearity), values of G cr for the SD and MM solitons are naturally found to be identical. Typical examples of such stable solitons at the boundary of their stability area are shown in Figs. 2(a3,b3). Figure 3(a) clearly demonstrates that the stability area expands with the increase of N , hence stronger nonlinearity makes the solitons more stable against the action of the gravity. Additional numerical results reveal that G cr is smaller at both γ < 1 and γ > 1 than in the Manakov's case. In particular, Fig. 3(a) shows the G cr (N ) curve for γ = 0, with values G cr essentially smaller than for γ = 1. Thus, as mentioned above, the Manakov's nonlinearity is optimal for the stabilization of the solitons against gravity. Indeed, it is natural to expect that effects of the SO coupling are strongest when the nonlinearity is equal for the nonlinear self-and cross-interactions of components of the pseudo-spinor wave function.
Further, Fig. 3(b) shows the dependence of G cr on size L of the confinement region in which the SO coupling is applied, as per Eq. (2), for the SD soliton with N = 5. The dependence features a maximum at L ≈ 0.6, which is, quite naturally, comparable to the width of the corresponding soliton, W ≃ ∼ 0.3. The numerical results demonstrate that  G cr decays as L −1 with the increase of L. This finding is explained below by means of the analytical approximation, see Eq. (15).
Finally, the growth of G cr with the increase of the strength of the localized SO-coupling term, i.e., Λ 0 in Eq. (2), is displayed in Fig. 3(c). Roughly linear initial increase of G cr with Λ, as well as its increase with N , observed in the figure, also complies with Eq. (15) derived below. The increase of Λ 0 leads not only to the growth of critical gravity force, G cr , but also to growth of the amplitude of the nonzero tail of the wave function extending towards the right edge of the solution domain at G > G cr (not shown here in detail).

The analytical approximation
Characteristics of the stationary state under the combined action of gravity and SO-coupling-induced trap may be analyzed by means of the approximation which treats the soliton as a quasi-particle. To illustrate the application of the approximation in a simple form, in Appendix we formulate it for the simplest setting based on the single GPE, with the localized trap represented by an attractive delta-functional potential. Here we apply similar analysis to the system under the consideration.
If the soliton's width is L, then the effective soliton's potential, induced by the localized SO coupling as per Eq. (3), is approximated as cf. Eq. (A5), where ξ is the coordinate of the soliton's center, U 0 is a constant, and use is made of the fact that the soliton's width and norm scale with its amplitude A as W ∼ A −1 and N ∼ A 2 W . This potential gives rise to the trapping force, The maximum value of the force is attained at the value of ξ determined by condition and accordingly, the largest absolute value of the force corresponds to The competing gravity force acting on the soliton is estimated as where F 0 is another constant. Finally, the largest value of G, up to which the equilibrium condition, F SOC (ξ)+F G = 0, may hold, is cf. Eqs. (A6) and (A7) in appendix. Equation (15) explains the above-mentioned L −1 dependence for relatively large L. On the other hand, if L is too small, the SO-coupling trap is obviously weak, therefore G cr vanishes at L → 0.

B. Stationary states in the moving reference frame
Stationary solutions in the moving reference frame were obtained by solving Eq. (6) with velocity c = 0. The procedure is similar to that outlined above for c = 0: (i) For a given value of c, ITM is used to produce a stationary solution with G = 0, as the ground state.
(ii) Stable solitons with G = 0 are obtained by means of the continuation of the latter solution by small steps to G > 0, with the help of SOM because, as well as in the case of c = 0, the stationary soliton cannot be a ground state at G > 0.
An example of a stable soliton with c = 0 is displayed in Figs. 5(a1,a2). The stability is corroborated by direct simulations of the perturbed evolution, as shown in Fig. 5(b). Comparison of the field profiles displayed in Fig. 5(a) with their counterparts, shown in Fig. 2(b2) for the same values N = 4 and G = 0.12, but zero velocity, demonstrates that the moving soliton, remaining a stable solution, develops a complex structure, in comparison with the essentially real one at c = 0. On the other hand, Fig. 5(c) suggests that the total-density profile of the moving soliton, shown in terms of n(x ′ ) = |φ + (x ′ )| 2 + |φ − (x ′ )| 2 , varies quite slowly with the increase of the velocity at a fixed value of the total norm. The latter figure demonstrates that the overall shape of the soliton gradually becomes sharper, which enhances the strength of trapping by the localized SO coupling. Stability of the moving solitons of the SD type is summarized in the parameter chart in the (c, G) plane, which is displayed in Fig. 5(d). The symmetry of the stability area with respect to c ←→ −c corresponds to the abovementioned invariance of Eq. (6) for the replacement of c by −c. The figure demonstrates that the increase of |c| leads to additional stabilization of the trapped solitons (growth of G cr ). This effect may be explained by the fact that the linear-mixing terms ∼ icΛ (x ′ ) in Eq. (6), being proportional to the spatial-modulation coefficient, Λ(x), induce an additional effective potential which enhances the trapping action of the spatially confined SO coupling. In particular, considering stationary solutions to Eq. (6), ψ ′ ± = exp (−iµt) φ ± (x ′ ), with |φ | ≪ |φ + |, and focusing on the effect of terms ∼ icΛ (x ′ ) (neglecting, for the time being, terms ∼ dΛ(x ′ )/dx ′ ), the use of the second equation in system (6) makes it possible to eliminate φ − in favor of φ + , as φ − ≈ −icµ −1 Λ(x)φ + . Then, the substitution of this result in the first equation of system (6) produces an effective potential, Because soliton states exist at µ < 0, the potential given by Eq. (16) is indeed a trapping one, with the negative sign, thus helping to additionally stabilize the bound state. The action of the extra potential also explains the shift of the overall-density profile back towards x = 0 and sharpening of the profile with the increase of c in Fig. 5(c).
FIG . 6: (a,b) Profiles of the two components of the moving soliton with the same values (N, L, γ) = (4, 2, 1) as in Fig. 5(a1,a2,c), but much larger gravity force, G = 1.36 (instead of G = 0.12 in Fig. 5). This soliton is located close to the stability boundary, in terms of Fig. Fig. 5(d). (c) The comparison of the total density of the same soliton (the red short dotted curves) and its counterpart pertaining to G = 0.12 (the solid black line).
Lastly, while the moving solitons shown in Fig. 5(a1,a2,c) are located deep inside the stability area, in terms of Fig. 5(d), Fig. 6 displays an example of a soliton with the same values of the norm, N = 4, and speed, c = 4, but taken at a much stronger gravity force (G = 1.36 instead of G = 0.12), which places the soliton close to the stability boundary. Comparison of the latter profile with its counterpart from Fig. 5(a1,a2), see panel (c) in Fig. 6, provides another proof of the robustness of the solitons: the increase of G by a factor ≈ 11.3 leads to a minor change in the shape, viz., reduction of the overall peak density by a factor ≈ 0.81, and a shift of the peak towards x > 0. The dynamics of solitons in the present system may be tested by abruptly applying gravity to a mode created in a weightless environment, i.e., generated by Eq. (1) with G = 0 (in the experiment, this may be realized by creating a soliton in a horizontally oriented cigar-shaped trap, which is then quickly turned in the vertical plane). Figure 7 exhibits a typical example of the subsequent evolution of a soliton of the MM type, under the action of the suddenly applied gravity with different values of G. It is seen that the soliton develops shuttle motion with an amplitude and period depending on G (oscillatory motion of nonlinear SO-coupled wave packets in a trapping potential was considered in Ref. [100]). The resulting period of the shuttle motion for the SD and MM solitons is shown in Fig. 8 as a function of G. The figure demonstrates that solitons of both types develop identical dynamics in the case of the Manakov's nonlinearity, γ = 1. The divergence of the period at the critical value, G = G d cr , implies that, at G > G d cr , the soliton can no longer be held in the shuttle state by the spatially confined SO-coupling trap, see an example in Fig. 7(d)]. The dynamical critical value of the gravity strength, G d cr , is presented, as a function of N and L, in Figs.  3(a,b). Shapes of dependences G d cr (N ) and G d cr (L) are similar to their counterparts for the stationary solutions, while the magnitude of G d cr is, naturally, smaller, as maintaining stable solitons in the dynamical regime requires to add some margin to the critical value defined for the static state. Long-time dynamics in the shuttle regime may be characterized by the soliton's center-of-mass coordinate, defined as N and G. It is observed that the oscillations are, generally, damped, which may be explained by weak radiation losses of the soliton moving with acceleration [101]. The damping rate increases with the increase of G and decrease of N . The former peculiarity is explained by the fact that rapid oscillations at smaller values of G (see Fig. 8) produce little radiation [101], while, on the other hand, heavier solitons with a larger norm are less amenable to the action of the weak recoil force produced by the emission of radiation.

V. CONCLUSION
The objective of this work is to study the dynamics of two-component matter-wave solitons under the action of the gravity field, with the aim to hold the solitons in a fixed position, or move them at a constant speed, by means of an effective trapping potential imposed by the spatially confined SO (spin-orbit) coupling. In comparison with usual laser-beam tweezers, advantage offered by the SO-coupling trap is the absence of detrimental heating effects. First, stationary solutions of the system of coupled GPEs for solitons of the SD (semi-discrete) and MM (mixed-mode) types were produced, in the quiescent and moving reference frames. In both cases, the critical value of the gravity force, G cr , which is the boundary of the stability area for the solitons trapped by the localized SO coupling, was found. By means of systematically collected numerical results and the analytical approximation, which treats the soliton as a quasi-particle, it was demonstrated that the stability area can be expanded (i.e., G cr made larger) by increasing the strength of the SO-coupling term, the norm, and the speed, in the case of the moving system. In the second part of the work, dynamical regimes were addressed, by abruptly applying gravity to solitons created in the weightless setting. If the gravity force is smaller than the respective dynamical critical value (which, in turn, is smaller than its stationary counterpart), the solitons feature damped shuttle motion, with the period and damping rate affected by the gravity strength. The results produced by the consideration of both stationary and dynamical settings demonstrate the ability of the effective trapping potential, induced by the spatially confined SO coupling, to hold solitons against the gravity, and transfer them at a constant speed, the latter regime being actually more stable. These findings may help to design new gravimeters and improve other techniques for precise measurements, based on the use of cold atoms.
The present analysis may be extended in other directions. First, a natural possibility is to consider this problem in the 2D geometry, where conditions for the soliton stability are drastically different from those in the 1D case [78]- [87]. The 2D matter-wave solitons stabilized by the SO coupling include vortex components, which suggests a possibility to consider the interaction of vortex solitons with gravity. Further, one can add the beyond-mean-field Lee-Huang-Yang corrections to the GPE system [102], and thus consider holding and transfer of "quantum droplets" (self-trapped modes with the flat-top shape [24,102]) in the presence of gravity. Note that the interplay between the spatially uniform SO coupling and LHY corrections was considered in Refs. [103][104][105]. A challenging option is to extend the current setting to the full 3D geometry.
If the gravitational and trapping potentials in Eq. (A2) are considered as small perturbations, the respective terms in the soliton's energy are U ε (ξ) + U G (ξ) = − (N/2) 2 ε sech 2 (N ξ/2) − GN ξ, which predicts a stationary state at a local minimum of the potential, defined by (d/dξ) (U G (ξ) + U ε (ξ)) = 0, i.e., at ξ which is a smaller root of the equation As follows from Eq. (A6), the stationary state exists at values of the gravity force Alternatively, one can conclude that, for given G, the trapping potential is able to hold solitons with norms exceeding a certain minimum value,