Dynamics of vortex quadrupoles in nonrotating trapped Bose-Einstein condensates

Dynamics of vortex clusters is essential for understanding diverse superfluid phenomena. In this paper, we examine the dynamics of vortex quadrupoles in a trapped two-dimensional (2D) Bose-Einstein condensate. We find that the movement of these vortex-clusters fall into three distinct regimes which are fully described by the radial positions of the vortices in a 2D isotropic harmonic trap, or by the major radius (minor radius) of the elliptical equipotential lines decided by the vortex positions in a 2D anisotropic harmonic trap. In the “recombination” and “exchange” regimes the quadrupole structure maintains, while the vortices annihilate each other permanently in the “annihilation” regime. We find that the mechanism of the charge flipping in the “exchange” regime and the disappearance of the quadrupole structure in the “annihilation” regime are both through an intermediate state where two vortex dipoles connected through a soliton ring. We give the parameter ranges for these three regimes in coordinate space for a specific initial configuration and phase diagram of the vortex positions with respect to the Thomas-Fermi radius of the condensate. We show that the results are also applicable to systems with quantum fluctuations for the short-time evolution.

experimentally only in the past few decades. With the development of atom chip technology 37 and recent experimental achievements in the monitoring of the dynamics of vorticity states in BEC systems in real time 38 , there has been fast-growing interest in vortex dynamics, and the new physical effects associated with the trap geometry.
A vortex quadrupole is a vortex cluster with zero total vorticity constructed as the product of four singly charged vortices. Previous study showed that this cluster configuration is both energetically and dynamically unstable 49 . However, for some special configurations, stationary vortex quadrupoles were found to exist 53,49 . It has shown that a vortex quadrupole is always structurally unstable against perturbations in some external trap parameters 63 . But a complete picture of quadrupole dynamics has not been given, which provides critical information for understanding how the vortices influence superfluid behaviors, such as evolution to turbulence and decay processes.
In this paper we study, by numerical simulations, the detailed dynamics of two-dimensional (2D) vortex quadrupoles with different initial configurations in trapped BECs, and for the first time we identify the parameter ranges which give the distinct regimes for the revival and disappearance of a quadrupole structure. The paper is organized as follows. In theoretical model section we introduce the model under investigation and define the general initial state of a vortex quadrupole. The method to introduce quantum fluctuations into the systems is also presented. In numerical results and analysis section the trajectory of each vortex in a quadrupole is investigated in detail for isotropic, anisotropic and fluctuating condensates. The mechanism of dipole formation and annihilation is discussed. The ranges of parameters of the initial configurations of quadrupoles corresponding to different dynamical regimes are calculated. A brief conclusion is also given.

Theoretical model
In the frame of zero temperature mean-field theory, a trapped BEC system is well parameterized by ψ(x, y, z, t) which satisfies the Gross-Pitaevskii (GP) equation , the condensate can be regarded as a 2D disk-shaped cloud, where the motion in the z-direction is squeezed. The 2D order parameter ψ(x, y, t) is governed by the corresponding 2D-GP equation, where the coupling constant g 3D in Eq. (1) is substituted by In our simulations of a 87 Rb BEC containing N = 1.5 × 10 5 atoms, the bulk s-wave scattering length is a s = 5.4 nm. For the trap frequencies, we used a 2D disk-shaped configuration, where ω x = ω y = 2π × 5 rad · s −1 and ω z = 2π × 100 rad · s −1 . It is convenient to turn all quantities dimensionless by introducing a time scale , which can provide better precision in simulations. For these parameters, the 2D Thomas-Fermi (TF) radius of the condensate is 0 . The order parameter of the system can be well described by a single particle wave function where ρ(r, t) and θ(r, t) are density and phase distributions of atoms respectively, with r = (x, y). Vortices can be seeded into the system by employing the phase imprinting techniques 24,64 experimentally. Numerically, the quadrupole states of interest in our simulations are obtained by using the method in refs 41 and 65. Firstly the vortex-free ground state of the condensate, ψ g , is obtained by the imaginary time evolution. The initial phase factor of the order parameter of the quadrupole state is j j j j where s j = ± 1 is the topological charge of the j th vortex, imprinting the velocity field of a vortex located at (x j , y j ). The radial position of each vortex is then = , the quadrupole vorticity state of a BEC system can be obtained by evolving the GP equation again in imaginary time while forcing the phase distribution during the calculation to be θ(r, t = 0) until the solution is numerically well converged, i.e. the relative change in norm of number of atoms and in energy of the system is smaller than 10 −6 . The converged quasi-equilibrium state is then the desired quadrupole configuration ready for a real-time evolution. In this paper we focus on a subset of possible configurations, namely, the initial vortex quadrupole structures are axisymmetric. An example density profile and the corresponding phase diagram are shown in Fig. 1e, where the vortices with opposite charges are nested in the condensate alternately forming dipole structures along each axis. The thermal and quantum fluctuations are experimentally only important for temperatures which are larger or comparable to the mean-field interactions. Since we are studying systems in the frame of zero temperature mean-field theory, it is reasonable to think that these fluctuations will not affect the dynamics significantly, especially during the short-time evolution. However, taking the symmetry of systems into account, quantum fluctuations might have obvious effects on the long-time evolution of the systems. To introduce quantum fluctuations into our systems, the coordinate space order parameter equivalent to our choice of stochastic initial state is easily shown to be where ψ 0 (r) is chosen to be the initial state containing the quadrupole structure used in the situations without fluctuations, i.e. wave function of the condensed atoms. ξ(r) is a complex Gaussian random field, which represents the uncorrelated quantum fluctuations 66,67 In practice, one can implement this initial condition by utilizing a complete set of plane-wave basis ik r 1 j , such that the random field can be written as Here ξ j are complex gaussian random variables By using the truncated Wigner method 68,69 , and integrating the fully coordinate space stochastic differential equation of a 2D system T 2 which includes both the condensed atoms and the virtual particles introduced into the field. It is obvious that even though Eq. (11) has a form similar to that of the GP equation (1), the important difference is that quantum fluctuations (virtual particles) are included in this stochastic differential equation while the GP equation is only valid for the condensed atoms.
In our calculations, only the 'low-energy' modes in momentum space contribute to the fluctuations, which means the noise is just added into the modes with amplitude ξ ∈ j L , j , where L = [− 25, 25] is a subspace containing 51 × 51 modes, only about 1% of the 501 × 501 momentum space used for our simulations. All other ξ ∉ j L , j are set to be 0.

Numerical results and analysis
To get insight into the dynamics of vortex quadrupoles in a disk-shaped BEC, we first review the motion of a single vortex and a vortex dipole in such a system. A 2D disk-shaped condensate in an axisymmetric harmonic trap is a nonuniform system. In the TF limit, for a single vortex generated in this nonrotating system the precession velocity is a function of the vortex displacement r from the centre of the condensate as is the metastable frequency for the appearance of a central vortex in a disk-shaped condensate 70 , with R ⊥ being the mean condensate radius and R c being the vortex core radius. It is obvious that a central vortex does not undergo any evolution over time. At zero temperature, an off-centre vortex with positive (negative) charge will perform a circular motion around the centre of the condensate in a anti-clockwise (clockwise) direction in a constant radius with a constant energy 60,71 without considering the interaction between the vortex and the excited phonon mode. This pattern of behaviour repeats as the density wave oscillations induced by the vortex motion are all confined within the condensate, and no net phonon radiation is expected. However, due to the interaction of the vortex with the collective phonon mode, the vortex trajectory is not a perfect circle 70,72 . At finite temperatures it spiral out of the condensate due to its interaction with the thermal component 65 .
For a dipole structure of vortex-antivortex pair in a trap, the dynamics is the result of the competition between the density inhomogeneity driven and the dipole interaction driven (mutually driven of each vortex of the pair). In the case where the mutually driven and inhomogeneity driven balance with each other, a stationary dipole structure generates. For configurations, where a dipole placed symmetrically about the trap centre, vortices make close orbits about fixed points 46 . Nearly the same as the motion in a single vortex system, each of the vortices begins its orbit along the equipotential line as if they move independently, but once the vortices begin to approach each other they begin to interact. Instead of annihilating each other, both vortices alter their paths and begin to travel down the condensate from one side to the other side, parallel to each other. This continues until the vortices begin to reach the boundary of the condensate, then both of them move back on to their original orbits.
The situation gets more complicated when a quadrupole structure is taken into account. By setting different initial configurations, i.e. varying the radial position of the vortices, we find that the dynamics of the vortices shows three distinct regimes, which are named as "recombination", "annihilation" and "exchange".
In the recombination regime as shown in Fig. 1a, the vortices take quarter-circular motion and near-linear motion repeatedly, where the solid dots are the initial positions of the vortices, and the arrows show the direction of movement of each vortex. We can see clearly that, when they approach each other, two vortices in the left half-plane move as a dipole towards the trap centre along the x-axis, while the other two in the right half-plane make a dipole movement as well facing to the left pair. However, after the four vortices reach their minimal separation, two of them in the upper half-plane combine into a new dipole structure and move away the trap centre towards the edge of the condensate. So do the remains in the lower-half plane. The trajectories of the upper and lower dipoles are mirror symmetric to the x-axis. During the movement of vortices, the density wave oscillations arise across the condensate, and are severe around the edge of the condensate where the density is low. Therefore, the oscillation of the quarter-circular parts of the vortex trajectories in the low density region is obvious as seen in Fig. 1a-d. For different initial positions the time evolution of the x− and y-coordinates of the down-left vortex is plotted in Fig. 1b. Typical density and phase distributions in this regime are plotted in Fig. 1e-k in time series. We note that for given parameters in Fig. 1a-c and e-k four vortices never contact. After the first period the vortices will move away from their initial position a little bit due to the interactions with the density wave oscillations arising in the system, unlike a single vortex or a vortex dipole movement which generates close orbits. In Fig. 1c we plot time series for the radial positions of the s = − 1 vortices (blue dots). If we increase the initial radial displacement of vortices, the dynamics of the quadrupole can still stay within the recombination regime for the first period. But the dipole annihilation and revival can occur, which makes the deviation of the vortex positions from the initial positions sizable as seen in Fig. 1d ((x, y) = (± 4.75, ± 3.35)). During the following time evolution, this may induce a transition of the dynamical regime from recombination to exchange (annihilation). Figure 2 shows the trajectory of the down-left vortex (s = − 1) for the first period in the exchange regime. From the symmetry of the system we know that the s = 1 vortices will give the same trajectories but with different starting position and in opposite direction (anticlockwise). Therefore the left dipole and the right dipole move towards each other along the x-axis to the trap centre at the beginning, which is the same as the situation in the recombination regime. As they approach each other the distance between the two vortices of the left (right) dipole decreases gradually in the y-direction as shown in Fig. 2f. However, instead of separate in the y-direction after getting their minimal separation in the recombination regime, the two dipoles join together, and finally The vortices continue to move in a quarter-circular orbit until they approach each other again. But at this time two vortices in the upper half-plane and the other two in the lower half-plane combine into two dipoles respectively. They move in the y-direction and repeat what happened in the x-direction. After four times charge flipping, the system will restore to its original configuration.
The dynamics in the annihilation regime is the most simplest, where the vortices vanish permanently as long as they collide in the trap centre for the first time as shown in Fig. 3a. In this process there is also a ring structure formed. However, different from the situation in the exchange regime, the ring expands across the condensate quickly and smear at the edge of the condensate, which induces severe density wave oscillations as shown in Fig. 3k.
As is well known, the energy and angular momentum associated with a vortex increase as the vortex is moved to higher density region 70 , which coincides with the fact that the energy cost to create a vortex is high when the local density is high 41 . So we note that in all the three regimes, there are two timescales in the motion: the slow quarter-circular precession around the trap centre and the fast motion nearly parallel to the axis. We note that the density wave oscillations are greater after the annihilation of vortex dipoles and the decay of the soliton ring, which indicates that the decay of vortex dipoles or soliton rings can make much severe excitations in dynamical systems. The revival of vortex dipoles can depress the excitations in the system as stated in refs 34,35. In Fig. 4 we identify the coordinate parameters which separate regions of "recombination" (circles), "annihilation" (dots) and "exchange" (crosses) dynamics. The boundaries between these regions are both nearly quarter-circles which are denoted by the dashed-lines. For a relatively small radial position r j < 5.91a 0 the dynamics of the quadrupole is characterized by the combination of the vortices into dipoles in x− and y-directions alternately, while for a large radial position r j > 6.02a 0 exchanging of positions of the vortices and charge flipping occur during the dynamical process. In between these two regions both the quadropole and dipole structures vanish due to the annihilation of the vortices with opposite charges. For r j > 6.6a 0 (black solid curve in Fig. 4 the density of the condensate is too low to identify a complete vortex structure. In Fig. 5, we give the radial position indicating the boundaries of three dynamical regimes with respect to varying TF-radius, which is applicable at any time step during the dynamical process. We find that the fitting curves of the boundaries (the blue and red lines) are both proportional to R TF .
Through the analysis of isotropic systems, we know that the quarter-circular precessions are nearly along the equipotential lines, which are concentric ellipses in anisotropic traps (with different trapping frequencies in x− and y− axis, ω x /ω y = 1.25). We find that the dynamics of vortices still depends on the configuration of the quadrupole structure in anisotropic systems. The parameter ranges can also be determined according to Fig. 5. However, r which is chosen to identify the dynamical regimes is not the radial position of the vortices as in the isotropic systems, but the major radius or minor radius of the ellipse decided by the vortex positions, since the TF-radius of the condensate in each direction is different. By releasing the trapping potential in one direction (both directions) suddenly or adiabatically (the trapping frequencies are decreased linearly within 10t 0 ) to a smaller value (form 2π × 5 Hz to 2π × 3 Hz or 2π × 4 Hz in our calculations), we confirm that the dynamics of quadrupole structures still obeys the same rule, which supports the conclusion in ref. 75, where they showed that the sudden reduction in the number of atoms following expansion is expected to excite a small-amplitude  To examine the robustness of the quadrupole dynamics, we introduce quantum fluctuations into the systems we studied above as stated in Sec. 2. We find that in the recombination regime with fluctuations the trajectories of vortices are nearly the same as those in the unperturbed systems during the first few periods. However, after that the quadrupole structure is eventually destroyed and the vortices take chaotic trajectories as shown in Fig. 6a. In the exchange regime with fluctuations, the dynamics of the vortices still obeys the rule of the unperturbed systems but with an additional rotation of the quadrupole structure as shown in Fig. 6b. This rotation also occurs if a dipole is placed asymmetrically about the origin 46 . Moreover the vortices will not maintain a rectangular configuration anymore, leading to the destroy of the quadrupole structure and chaotic trajectories of vortices eventually. In the annihilation regime the quantum fluctuations introduced will not change the dynamics of the system since the vortices annihilate permanently as soon as they meet at the trap centre during the first period.
In the frame of zero temperature mean-field theory we expect that quantum fluctuations will not change the dynamics of the system dramatically especially for the short-time evolution. However, symmetries of the system are affected by quantum fluctuations, which leads to chaotic motion of vortices for the long-time evolution. As we focused on the first few periods of quarupole dynamics in our discussions, it shows that our results are robust, i.e. quantum fluctuations do not change the dynamical regimes of the system we studied at least for the short-time evolution. This also demonstrates that vortices have very stable topological structures.
We note that even though our calculations are limited at zero temperature, the vortices of the quadrupole cluster do deviate from their original locations during the periodical movement, which is similar to the situation in finite temperatures. The dynamics of the next period is fully decided by the fact that which region in Fig. 5 the  radial positions (major radius or minor radius) of the vortices will fall into. Transition from one regime to another of the system dynamics may occur during the time evolution if the vortices do not annihilate permanently, which means that the initial configurations whose dynamics is ruled by "recombination" in the first period may transit to the other two regimes in the following period. And this may finally result in expulsion of the vortices from the condensate.

Conclusion
In this paper, we have shown that there are three characteristic dynamical regimes in a trapped BEC according to the varying radial positions of vortices in a quadrupole structure. The closer the vortices are located to the periphery of the condensate, the higher the kinetic energy they will obtain when approaching the centre of the condensate, which will compete with the interatomic interactions and the vortex interactions and leading to distinct dynamical regimes. In the recombination regime the kinetic energy is too small to overcome the interatomic interactions, which prohibits the formation of soliton rings, and the vortices recombine into dipoles. In the annihilation regime the ring structure arises and expands fast out of the condensate before it decays into vortex dipoles, while in the exchange regime the emerged soliton ring can break into dipoles with charge flipping due to higher kinetic energy possessing by each vortex. These are also valid for anisotropic systems and systems with quantum fluctuations during the short-time evolution, where the symmetry of the system is altered. For the long-time evolution, the symmetry of the system is important for maintaining the rectangular configuration of a quadrupole structure.
By providing the relation between TF-radius of the condensate and the configuration of the vortex quadrupoles, our results disclose detailed dynamics of vortex quadrupole in a trapped BEC and a road map to distinguish them, and can be extended to other special structures, which helps to study more complex vortex settings. The numerical simulations also demand new theoretical models about the generation and expanding speed of soliton rings, by the combination of individual vortices, which are crucial to understand the dynamics of vortex clusters.