Mean-field model of interaction between bright vortex solitons in Bose–Einstein condensates

Using the explicit numerical solution of the axially symmetric Gross–Pitaevskii equation we study the dynamics of interaction among vortex solitons in a rotating matter-wave bright soliton train in a radially trapped and axially free Bose–Einstein condensate to understand certain features of the experiment by Strecker et al (2002 Nature 417 150). In a soliton train, solitons of opposite phase (phase δ = π) repel and stay apart without changing shape; solitons with δ = 0 attract, interact and coalesce, but eventually come out; solitons with a general δ usually repel but interact inelastically by exchanging matter. We study this and suggest future experiments with vortex solitons.


Introduction
Solitons are solutions of wave equation where localization is obtained due to a nonlinear interaction. Solitons have been noted in optics [1], high-energy physics and water waves [2], and more recently in Bose-Einstein condensates (BEC) [3]- [5]. The bright solitons of BEC represent local maxima [4]- [6], whereas dark solitons represent local minima [3,7]. In addition to the observation of an isolated bright soliton in an expulsive potential [5], a number of bright solitons constituting a soliton train were also observed in an experiment by Strecker et al [4], 137.2 where they suddenly turned a repulsive BEC of 7 Li atoms attractive. Consequently, the BEC collapsed, exploded and generated a soliton train. Similar collapse and explosion in 85 Rb BEC have been studied before from a different point of view [8]. These experiments were performed by manipulating the background magnetic field near a Feshbach resonance [9]. It was found [4] that solitons in such a train usually stay apart. Also, a soliton was often found to be missing from a train [4]. Hasegawa [10] considered the generation of a train of optical solitons by an induced modulational instability. There have been theoretical attempts [11]- [15] to simulate essentials of these experiments [4,5] on bright solitons in BEC. Al Khawaja et al [11] used a variational approach to describe the experiment by Strecker et al [4] whereas Salasnich et al [12] and Leung et al [13] used an effective one-dimensional model for the same purpose. In a more recent work Salasnich et al [14] considered a three-dimensional mean-field model to study the production of and interaction between bright solitons to account for different aspects of the experiment by Strecker et al [4]. Carr and Castin [15] used an approximate approach to describe the experiment of Khaykovich et al [5]. Elyutin et al and Shchesnovich et al [16] studied solitons in one dimension.
We use the explicit numerical solution of the axially symmetric mean-field Gross-Pitaevskii (GP) equation [17] to study the dynamics of bright solitons in a soliton train [4]. It is found that in an axially symmetric configuration with no axial trap, the GP equation permits solution in the form of a radially confined soliton train. We also give an explanation of missing solitons in the experiment of Strecker et al [4].
The present approach is also extended to the study of vortex solitons with an angular momentum ofh per atom in the axial direction. Vortex solitons are rotating solitons of an attractive condensate. Due to experimental observation [18] of a vortex state in a rotating BEC, the experimental generation of vortex solitons seems possible. Attractive BECs may not form vortices in a thermodynamically stable state. However, such vortices may be created via a Feshbach resonance [9]. Due to the conservation of angular momentum, a vortex soliton train could be generated by suddenly changing the interatomic interaction in an axially symmetric rotating vortex condensate from repulsive to attractive near a Feshbach resonance in the same fashion as in the experiment by Strecker et al [4] for a non-rotating BEC. Alternatively, a single vortex soliton could be prepared and studied in the laboratory by forming a vortex in a small repulsive condensate and then making the interaction attractive via a Feshbach resonance and subsequently reducing the axial trap slowly. Already, there have been theoretical considerations about these vortex states [19,20]. It would be worthwhile exploring these possibilities experimentally.
In particular we study in some detail the interaction between two bright vortex solitons for different relative phases between them. The interaction is usually found to be repulsive and inelastic with exchange of particles. The interaction turns attractive for small values of relative phase. In one dimension the interaction is usually elastic without exchange of particles [2].
In section 2 we present the mean-field model of solitons that we use in the numerical analysis. In section 3 we present the results of our numerical study. Finally, in section 4 we present the conclusions.

Nonlinear mean-field model for solitons
In the following one-dimensional nonlinear free Schrödinger equation in dimensionless units solitons are bound states due to the attractive nonlinear interaction with wavefunction at time t and position y: (y, t) = √ 2| | exp(−i t)sech(y √ | |), with the energy [21]. Equation (2.1) can sustain any number of such solitons at different positions. Although, the one-soliton solution of (2.1) may remain stationary at a fixed position, the many solitons of a multiple-soliton solution of this equation generally move around because of the interaction among them [22]. The time-dependent BEC wavefunction (r; τ ) at position r and time τ is described by the following mean-field nonlinear GP equation [17] −ih where m is the mass and N the number of atoms in the condensate, g = 4πh 2 a/m the strength of interatomic interaction, with a the atomic scattering length. For an axial trap V (r) = 1 2 mω 2 (r 2 +λ 2 z 2 ) where ω is the angular frequency in the radial direction r and λω that in the axial direction z, with λ the aspect ratio. The normalization condition is dr | (r; τ )| 2 = 1.
In a quantized vortex state [20], with each atom having angular momentum Lh along the z axis, (r, τ ) = ψ(r, z, τ ) exp(iLθ) where θ is the azimuthal angle. Now transforming to dimensionless variables where nonlinearity n = N a/l. For solitonic states n is negative. In terms of the one-dimensional probability P(y, t) defined by the normalization of the wavefunction is given by ∞ −∞ dy P(y, t) = 1 We solve the GP equation (2.3) numerically using a variation of the split-step time-iteration method using the Crank-Nicholson discretization scheme described recently [23]. Typical space and time steps for discretization are 0.1 and 0.001. The variation of the standard approach is required for λ = 0. The time iteration is started with the known harmonic oscillator solution for a small λ ≡ λ 0 ≈ 0.05 and nonlinearity n = 0: [20]. The desired value of nonlinearity n is then slowly switched on and λ slowly switched off from λ = λ 0 to 0 in the course of time iteration. The solution then corresponds to the trapped BEC for λ = 0. Then, without changing any parameter, this solution is iterated several thousand times so that a solution for λ = 0 is obtained independent of the initial input λ = λ 0 .

Numerical results
A classic soliton in three-dimensional BEC can be realized for an attractive nonlinear potential (n < 0) by setting λ = 0 [6] in (2.3). For a fixed L, the BEC is then governed by the single parameter n. The absence of a trap in the axial y direction will allow free movement of the solitons in this direction. Consequently, the study of soliton interaction will be trap independent. A localized BEC so created should be the three-dimensional analogue of the one-dimensional soliton. However, for calculational or experimental convenience, in the recent studies some weak potential was applied in the axial y direction. In the classic experiment of Strecker et al [4] an optical trap was maintained in the y direction. In the recent theoretical investigation by Al Khawaja et al [11] on this experiment a weak harmonic trap was applied in this direction, whereas Salasnich et al [12] employed infinite walls in an effective one-dimensional model. In the experiment by Khaykovich et al [5] and in the related theoretical study by Carr and Castin [15] an expulsive potential (λ 2 < 0) was applied in the y direction. We note that for λ 2 < 0, only a meta-stable and no stable soliton of (2.3) is possible as the potentials in this equation including the attractive nonlinear term do not provide confinement [15]. For the same reason no bright soliton can be generated for n ≥ 0 in (2.3) (repulsive condensate).
Although, under the conditions n < 0 and λ = 0, the potentials of (2.3) lead to confinement, a soliton-type BEC state can be generated only for n greater than a critical value (n cr ): n cr < n < 0. For n < n cr , the system becomes too attractive and collapses and no stable soliton could be generated. The collapse was first confirmed in the pioneering experiment by Gerton et al [24] for a trapped BEC of 7 Li. The actual value of n cr is a function of the trap parameter λ. For the spherically symmetric case λ = 1, and n cr = −0.575 [17,20]. The value of n cr should slightly change for λ = 0.
We solve the GP equation for solitons with λ = 0 and L = 0 and 1 [20]. We find numerically that the critical n for collapse of a single soliton is n cr = −0.67 for L = 0 and n cr = −2.10 for L = 1. The result for L = 0 is in good agreement with (a) the numerical result −0.676 obtained by Gammal et al [25] using the Crank-Nicholson method, (b) the numerical result −0.676 obtained by Pérez-García et al [6] using the steepest-descent method to minimize the mean-field Hamiltonian (these authors quote Q = −8π N a/l = 17 instead of n = N a/l) and (c) the result −2/3 obtained by Salasnich et al [12] using an approximate analytic onedimensional model.
By solving the three-dimensional GP equation Salasnich et al [12] found their result to be very accurate. However, there is some discrepancy between these results for L = 0 and 137.5 the value n cr = −0.6268 ± 0.0035 obtained by Carr and Castin [15] using the imaginary time relaxation method. The variational estimate n cr = −0.78 of Salasnich et al [12] seems not to be entirely consistent with other results. Further independent studies are necessary to resolve the discrepancy. For ω = 2π × 800 Hz and final scattering length −3a 0 as in the experiment of Strecker et al [4], n cr = −0.67 corresponds to about 6000 7 Li atoms. One can have proportionately about three times more atoms in the L = 1 state. An L = 0 soliton with n = −0.2 is illustrated in figure 1(a) where we plot the threedimensional wavefunction |ϕ(x, y)/x| versus x and y. Next we consider L = 1. In this case we calculated the soliton for n = −1 and plot the three-dimensional wavefunction |ϕ(x, y)/x| in figure 1(b). Because of the radial trap the soliton remains confined in the radial direction x, although free to move in the axial y direction. In either case the single soliton remains stable for more than 400 000 time iterations.
After having demonstrated the formation of a single soliton we next consider the dynamics of two solitons in a soliton train. Two solitons are prepared at positions y 1 and y 2 and then superposed with a phase difference δ. Specifically, we consider the following superposition of two normalized solitonic wavesφ at y = ±y 0 with phase difference δ between them with 2y 0 the initial separation between the solitons in the axial direction. To conserve the total number of atoms we normalize the superposed wavefunction (3.1) of the two solitons to 2 as it contains twice as many particles as a single soliton. In this fashion one can also construct the superposition of several solitons with a specific phase difference between them. The time evolution of the soliton train so formed upon superposition is found using the iterative solution of (2.3). We present results for a soliton train with two vortex solitons of the same nonlinearity with a phase difference δ of π, 3π/4, π/2, π/4, π/8 and 0 between them. Consideration of two equal solitons does not lead to a specialization and we shall see that for a general δ two equal initial solitons generally lead to two unequal solitons due to exchange of atoms.
In one dimension the solitons attract for δ = 0 and repel for δ = π [11,22]. In the present three-dimensional case for a general δ a more complicated motion emerges. In the repulsive δ = π case, two solitons stay away from each other. In the attractive δ = 0 case, they come close, coalesce, interact and come out. In a one-dimensional model of two trapped (λ = 0) 137.6 three-dimensional solitons, it has been shown that these solitons repel for δ = π [11,12]. In the following we study the rich dynamics of soliton interaction in a bright soliton train in three dimensions. The solitons can easily be set into motion by applying a perturbation or a constant (gravitational) force in y direction. However, for calculational convenience we chose not to do that and thus we studied the interesting relative motion between solitons suppressing their centre-of-mass motion. In all the cases studied the initial velocity of the solitons is zero. The numerical simulation was performed on a lattice 7 > x ≥ 0 and 250 > y > −250.
We studied the dynamics of interaction of two solitons in view of the experiment by Strecker et al [4]. For simplicity we first prepared two solitons of equal mass corresponding to the same nonlinearity n = −0.2 centred at points y = ±15 with angular momentum L = 0 and introduced them as the input to the GP equation with the initial phase difference δ = π . The time evolution of the train of two such solitons is exhibited in figure 2(a) where we plot the one-dimensional probability P(y, t) of (2.4) versus y and t. Because of mutual repulsion for δ = π , the two solitons stay apart and move away from each other. In an interval of time 200, the two solitons moved from positions y = ±15 to ±25, respectively.
We studied soliton trains with both L = 0 and 1 solitons and they lead to physically similar results. In the rest of the study we consider only L = 1 vortex solitons as they have never been studied before. Also, an effective one-dimensional model as usually employed in other studies does not seem to be applicable to the study of a vortex soliton. The rotational degree of freedom responsible for the generation of vortex does not exist in a one-dimensional model. We consider four L = 1 vortex solitons each with n = −0.4 at positions ±45 and ±15 with phase difference δ = π between two neighbouring and, hence, repelling, solitons. Obviously, one can accommodate as many repelling solitons as one prefers in a train. Because of the repulsion, the solitons stay apart and move forward without interaction and maintaining shape. Such a soliton train is illustrated in figure 2(b). In an interval of time 200, the four solitons have moved from positions y = ±15 and ±45 to ±17 and ±52 respectively. Now we consider the interaction between two vortex solitons in some detail. For that we reduce the phase δ between two L = 1 vortex solitons at y = ±15 each with n = −0.4 from π to 0. In figures 3(a)-(f) we plot one-dimensional probability P(y, t) of (2.4) for two interacting vortex solitons versus y and t for δ = π , 3π/4, π/2, π/4, π/8 and 0 respectively. As δ is reduced from π to π/2 the two solitons continue to repel and stay apart without much change in shape. However, while δ is reduced from π the repulsive interaction between the solitons gradually becomes 'inelastic' and the exchange of atoms between the solitons increases, resulting in a change of shape. Because of exchange of atoms the motion of solitons in figures 3 is asymmetric in general except for δ = π and 0.
In figure 3(a) for δ = π in the interval of time 400, the solitons moved from y = ±15 to ±38 respectively, corresponding to a final separation between the two solitons of 76. In figure 3(b) for δ = 3π/4 in the interval of time 400, the solitons moved from y = ±15 to 35 and −41 respectively, again leading to a final separation of 76. Although the repulsion is same in both cases, an asymmetry has appeared in figure 3(b) due to inelastic collision with exchange of atoms. In figure 3(c) for δ = π/2 in the interval of time 400, the solitons moved from y = ±15 to 38 and −26 respectively, corresponding to a final separation of 64. As δ is reduced further a similar trend is maintained: reduction in repulsion due to inelasticity (exchange of atoms between the two solitons). In figure 3(d) we see that for δ = π/4 the two solitons maintain their identity, exchange atoms and the overall interaction continues to be repulsive. corresponding to a final separation of 49. For a further reduction of δ from π/4 the inelasticity and asymmetry of interaction reduces. The mutual overall interaction becomes less and less repulsive and finally becomes attractive as δ → 0. In figure 3(e) we show the relative motion of two solitons for δ = π/8. The two solitons have a tendency to combine to form a single soliton before they lead to two separate solitons. In time 400, solitons moved from y = ±15 to 22 and −11, respectively, corresponding to a separation of 33. The relative movement between the solitons has reduced monotonically as δ is changed from π to π/8. From figure 3(f) we see that for δ = 0 there is strong overall attraction and the two solitons interact and coalesce to form a single soliton which eventually gives birth to two solitons of the same shape as the initial ones. The final separation between the two solitons in figures 3(a)-(e) for δ = π, 3π/4, π/2, π/4 and π/8 respectively, are 76, 76, 64, 49 and 33 showing a gradual reduction in repulsion as δ is reduced from π to π/8. In all cases the final separation is larger than the initial separation 2y 0 (= 30) between solitons, reflecting an overall repulsion. The asymmetry in the position of the two final solitons is zero for δ = 0 and π and largest for δ = π/4. Two features of figures 3 deserve some comment. First, the value of the phase δ for which the interaction between two solitons is attractive. In the one-dimensional model study by Salasnich et al [12] it is attractive for all cos δ > 0 in agreement with a previous analysis [26] whereas in the present three-dimensional study with an additional radial trap there is no sharp transition from attraction to repulsion for a fixed δ. The presence of the radial trap in our study makes a slow transition from repulsion to attraction with the change of δ. However, there is clear attraction for δ ≈ 0 and repulsion for cos δ < 0 in both models. Further studies are needed for a complete understanding of the present three-dimensional solitons under transverse confinement. This transverse confinement is absent in the one-dimensional model of Salasnich et al [12]. It is true that in the one-dimensional model of Salasnich et al the radial trap frequency enters as a parameter (see their equation (2.1)), it does not include real three-dimensional nonlinear dynamics. The present three-dimensional solitons cannot be considered to be classic one-dimensional integrable text-book solitons of the type considered in [12,26]. Hence the text-book wisdom is not directly applicable to the present analysis.
The second feature is the asymmetry in figures 3 for an arbitrary δ because of inelastic exchange of atoms between the two solitons which is absent in the one-dimensional model [12]. Although we report the results for angular momentum L = 1 in this study we have verified that similar asymmetry is present for L = 0. The asymmetry is also present in the analysis of [14] of two interacting solitons. One of the two equal solitons increases and the other decreases in size due to inelastic exchange of atoms. In figures 3 the left soliton grows in size and moves more slowly, the right soliton becomes smaller in size and moves faster. The direction of particle exchange is determined by the way the phase is introduced between the two solitons in (3.1). We verified that as δ is changed to −δ the role of asymmetry is reversed. For a positive δ due to particle exchange in figures 3 the left soliton grows in size. For a negative δ the right soliton grows in size in a symmetric fashion so that P(y, t, −δ) = P (−y, t, δ). The asymmetry 137.9 of interaction is absent in figures 2 and 3(a), where the interaction is purely repulsive with no exchange of atoms.
To study systematically the asymmetry of the interaction between two solitons we plot in figure 4 the ratio N R /N as a function of time for the cases presented in figures 3. Here N R is the number of atoms in the right soliton and N the total number in two solitons. The variation of N R /N is quite similar to that found in [14]. For δ = π/4 in both studies N R /N first attains a minimum and then oscillates, but remains less than 0.5. A similar behaviour is noted in the present study for δ = π/8. For δ > π/4, N R /N attains a minimum and then increases a little to attain a constant value smaller than 0.5. However, a quantitative comparison between the two studies is not to the point as the initial conditions of the two studies are different, e.g. the initial size and angular momentum L and the initial separation between the solitons. The details of interaction dynamics should depend on these initial conditions. In the present simulation we find that, for the change δ → −δ, is the number of atoms in the left soliton.
Next we exhibit the profile of two interacting solitons undergoing particle exchange. From figure 4 we find that δ = π/4 leads to a large exchange of particles. For this purpose we consider the situation depicted in figure 3(d) corresponding to two L = 1 vortex solitons of n = −0.4 each, which are set at y = ±15 at time t = 0 with a phase difference of δ = π/4 according to (3.1). The profiles of the system at times t = 0, 100, 200, 300 and 400 are shown in figures 5(a)-(e), respectively. Due to exchange of particles the left soliton gradually becomes larger and larger consistent with figure 4. At t = 100 the solitons come closer, interact strongly by exchanging atoms and become wider and deformed. Then they gradually move apart and become narrower. Finally, at t = 400 they re-acquire Gaussian shapes. However, at t = 400 the right soliton is wider and shorter and accommodates a smaller number of atoms consistent with figures 3(d) and 4. At that stage the left soliton is narrower and taller and contains a larger fraction of the atoms.
The interaction between more than two solitons with a general phase δ among different pairs is much too complicated to be studied exhaustively here. If the phase difference δ between two neighbouring solitons is not close to zero, they experience overall repulsion and stay apart. However, for δ close to zero they interact attractively and a soliton could often be lost, as observed in the experiment of Strecker et al [4]. This is illustrated in figure 6 where we consider three solitons all with the same phase (δ = 0). Due to attraction the three solitons come closer, interact and form two solitons which after some time form a single soliton. This single soliton next decays to two, recombines to one, and decays to two again and the three original solitons are never recovered. It is also possible that for some values of nonlinearity n and initial separation 2y 0 two solitons with relative phase δ = 0 may coalesce to form a single soliton without ever decaying to two solitons again. These could explain the missing soliton observed in the experiment by Strecker et al [4].
Throughout this investigation, in the interaction of two equal solitons we assumed that the nonlinearity |n| for each is less than |n cr |/2, so that a stable solitonic condensate with total |n| < |n cr | exists when the two coalesce. However, if two solitons each with |n| > |n cr |/2 encounter for δ = 0, the system is expected to coalesce, collapse and emit atoms via three-body recombination [8]. It is possible that in this case only a smaller single soliton survives. This might also explain some missing soliton(s) in the experiments.
In the investigation by Al Khawaja et al [11], to justify the repulsion between solitons in a train [4] a phase difference of δ = π was suggested between neighbours. They made a model for the evolution of δ along the axial direction which they used for explaining the formation of a soliton train, with neighbouring solitons of phase difference π , after changing the scattering length from repulsive to attractive in a BEC as in the experiment [4]. Such an order in phase seems not to be necessary for an overall repulsion between solitons. Almost any δ, except δ near 0, is found to lead to overall repulsion. In [12], using an approximate analysis, it has been found that repulsive interactions require cos δ < 0. In the present three-dimensional analysis, we find that attraction starts at δ = π/2 and increases as δ is reduced. However, the overall interaction remains repulsive until a small value of δ close to zero is attained, when the solitons attract, lose their separate identity and coalesce.

Conclusion
To conclude, employing a numerical solution of the GP equation with axial symmetry, we have performed a realistic mean-field study of the interaction among two vortex solitons in a train and find the overall interaction to be repulsive except for phase δ between neighbours close to 0. For phase δ = π between neighbouring solitons in a train, they are found to repel and move away without exchanging atoms for both normal (L = 0) and vortex (L = 1) solitons. For phase δ = 0 in a two-soliton train, the solitons attract, collide, coalesce and eventually come out. For other δ, overall repulsion prevails. However, there is an inelastic exchange of atoms between two solitons resulting in a change of size and shape. These unequal solitons travel in general with different speeds: the smaller soliton travels faster and the larger soliton travels slower. The unequal speed leads to an asymmetry in final positions of the solitons. By changing the sign of the phase between the two solitons the asymmetry in the final size and shape of the two solitons can be reversed. Except in the δ ≈ 0 case, the solitons in a train stay apart and never cross each other as observed in the experiment by Strecker et al [4]. For δ ≈ 0 a single soliton can often disappear as a result of the attractive interaction among solitons, as observed experimentally by Strecker et al [4]. These features of soliton interaction are present for both normal (L = 0) [14] vortex (L = 1) solitons. Although the present study is performed in the absence of an axial trap as in one dimension [2], these general conclusions should remain valid for a weak axial trap as in the experiments [4,5]. The L = 1 vortex solitons can accommodate a larger number of atoms and the present study may motivate future experiments with them.