Three-dimensional splitting dynamics of giant vortices in Bose-Einstein condensates

We study the splitting dynamics of giant vortices in dilute Bose-Einstein condensates by numerically integrating the three-dimensional Gross-Pitaevskii equation in time. By taking advantage of tetrahedral tiling in the spatial discretization, we decrease the error and increase the reliability of the numerical method. An extensive survey of vortex splitting symmetries is presented for different aspect ratios of the harmonic trapping potential. The symmetries of the splitting patterns observed in the simulated dynamics are found to be in good agreement with predictions obtained by solving the dominant dynamical instabilities from the corresponding Bogoliubov equations. Furthermore, we observe intertwining of the split vortices in prolate condensates and a split-and-revival phenomenon in a spherical condensate.


I. INTRODUCTION
Quantized vortices are archetypal topological objects that play important roles in various branches of physics, ranging from superconductors [1] and helium superfluids [2] to cosmology [3] and optics [4].Quantized vortices exist in matter fields described by a smooth complex-valued scalar field.The essential idea is that, while the complex field itself is single valued, its phase is defined only modulo 2π.Hence, the contour integral of the phase around a closed loop need not vanish, but may in fact be any integer multiple κ of 2π.A nonzero κ implies the presence of a quantized vortex within the loop and is referred to as the winding number of the vortex.
Previous theoretical studies of vortex splitting have been limited to relatively small winding numbers |κ| ≤ 5 [16][17][18][19][20][21][22][23][24]28] or to quasi-two-dimensional models pertaining to highly oblate BECs [25][26][27]31].In Ref. [34], vortex splitting was studied in three dimensions up to κ = 45, but only for small BECs in isotropic harmonic traps.Splitting patterns exhibiting up to tenfold rotational symmetry were observed in the numerical simulations.In this work, we carry out a more comprehensive investigation of giant-vortex splitting in threedimensional BECs.Considering all three different types of cylindrically symmetric harmonic traps (oblate, spherical, and prolate) and a wide range of repulsive interaction strengths, we simulate the temporal evolution of axisymmetric giant vortex states subjected to small random perturbations.In general, we find good agreement between the splitting patterns observed in the evolution and those predicted by linear stability analysis.Vortex splitting in prolate BECs is found to result in branched intertwining of the vortices, and spherical BECs are observed to exhibit a split-and-revival effect.
Importantly, we also find that the splitting patterns appearing in the simulated time evolution can be prone to numerical artifacts stemming from the symmetry of the underlying spatial grid.As a result, particular care should be taken when discretizing the time-dependent Gross-Pitaevskii equation (GPE) for the condensate.Specifically, the Cartesian grids used in the previous investigations tend to favor the fourfold splitting pattern, which may explain why, in Ref. [26], the higher-symmetry splitting patterns predicted by the linear stability analysis were not observed to arise from random perturbations.We solve this problem by basing our time integration scheme on discrete exterior calculus [35][36][37] with tetrahedral tiling.
The remainder of this article is organized as follows: In Sec.II, we present the time-dependent GPE, derive the Bogoliubov equations used for the linear stability analysis, and outline our numerical integration method.Section III begins with an analysis of the integration method and presents our numerical results.Finally, we conclude the paper in Sec.IV.

A. Mean-field model
The complex-valued order parameter Ψ of a dilute BEC at low temperatures satisfies the GPE where i is the imaginary unit, is the reduced Planck constant, m is the atom mass, and g is the effective interaction strength.The order parameter is normalized such that |Ψ(r, t)| 2 d 3 r = N is the number of condensed atoms.We employ a cylindrically symmetric harmonic trapping potential V(r) = m ω 2 r r 2 + ω 2 z z 2 /2, where ω r and ω z are the radial and axial trapping frequencies, respectively.
To have generally applicable results, we employ dimensionless units and measure position in the units of the radial harmonic oscillator length a r = √ /mω r , time in units of 1/ω r , the order parameter in units of N/a 3 r , and the effective interaction strength in units of a 3 r ω r /N.Thus, the conversion into the dimensionless units (denoted by a bar) is given by r Consequently, the dimensionless order parameter is normalized as | Ψ(r, t)| 2 d 3 r = 1, and it satisfies the dimensionless GPE i∂t Ψ(r, t) = − 1 2 ∇2 + V(r) + ḡ| Ψ(r, t)| 2 Ψ(r, t). ( The dimensionless potential is given by V(r) = r2 + λ 2 z2 /2, where λ = ω z /ω r is referred to as the aspect ratio.In cylindrical coordinates, the Laplacian is given by ∇2 1) has stationary vortex solutions Ψλ,ḡ,κ , which depend on λ, ḡ, and the integer winding number κ.These stationary states can be written as where f is a real-valued function and µ is the chemical potential.The stationary vortex states satisfy the time-independent equation 1 2 which can be solved using a relaxation method [38].

B. Bogoliubov equations and stability
To study the local stability properties of a given stationary vortex solution Ψλ,ḡ,κ , we decompose the order parameter as where χ is a function describing a small perturbation such that |χ(r, t)| 2 d 3 r ≪ 1.By substituting Eq. ( 3) into Eq.( 1), neglecting the second-and third-order terms in χ, and seeking oscillatory solutions of the form χ(r, t) = q∈N l∈Z u q,l (r, z)e ilφ−iω q,l t + v * q,l (r, z)e iω * q,l t−ilφ , (4) we obtain the Bogoliubov equations where the linear differential operator is defined as The integer l specifies the angular momentum of the excitation with respect to the condensate, and q ∈ N is an index for the different eigenmodes with a given l.Equation ( 5) can be used to determine the stability characteristics of the stationary vortex state in question.If the excitation spectrum {ω q,l } contains at least one eigenfrequency with a positive imaginary part Im(ω q,l ) > 0, the state is dynamically unstable; otherwise, the state is dynamically stable.If the spectrum contains an excitation for which Re(ω q,l ) < 0 and |u q,l | 2 − |v q,l | 2 r dr dz ≥ 0, the state is energetically unstable; if no such excitations exist, the stationary state is (locally) energetically stable.We emphasize that energetic stability is a stronger condition than dynamical stability, since the former implies the latter.
As can be observed from Eq. ( 4), the occupations of excitation modes with Im(ω q,l ) > 0 are predicted to increase exponentially over time, and, consequently, small perturbations of a dynamically unstable stationary state typically lead to large changes in its structure.For dynamically unstable multiquantum vortices, in particular, the complex-frequency modes usually induce instability against splitting of the multiply quantized vortex into singly quantized ones.In fact, the quantity max l max q [Im(ω q,l )]/2π and the maximizing winding number l can be used to predict, respectively, the inverse lifetime of a vortex and the symmetry of its typical splitting pattern [18].Note, however, that the dynamically unstable modes quickly drive the system beyond the linear regime of the Bogoliubov analysis.As a result, the long-time dynamics of dynamically unstable states must be described with the time-dependent GPE, Eq. (1), instead.
Typically, these methods rely on Cartesian spatial discretization, even though there are strong reasons to prefer simplicial grids [46,47].
This work, on the contrary, utilizes a time integration method based on discrete exterior calculus (DEC) [35][36][37], which naturally segregates the differentiable and metric structures [48,49].This approach can be regarded as a generalized finite-difference technique that closely resembles the finite integration technique [50] or the finite-difference timedomain method [51,52].The DEC method is applicable to unstructured grids, while being stable and conserving the particle number.
The discretization is based on a pair of interlocked threedimensional meshes: a primal (Delaunay) mesh and its dual (Voronoi) mesh.We assign each dual node with a floating point number to obtain a column vector ψ k that represents the discrete order parameter at a time instance k∆t/2, where k is an integer and ∆t is the length of the time step.With the notation of Ref. [53], the discrete Laplacian is denoted as where ⋆ p is a diagonal matrix called the discrete Hodge and d 2 is a sparse matrix called the discrete exterior derivative.The time integration of Eq. ( 1) is carried out using the central-difference method Here r j and z j denote the radial and axial coordinates of the jth dual node.The method is numerically stable if ∆t < M −1 , where M is the maximal diagonal element of the matrix

A. Evaluation of time integration
First, we test our numerical solver by numerically integrating a stationary vortex state forward in time and investigating its stability during the simulation.We consider the normalized GPE (1) with parameters λ = 1, ḡ = 300, and κ = 10.The time integrator is initialized at time instances −∆t/2 and 0 by letting ψ k j = Ψλ,ḡ,κ (r j , k∆t), where k = −1, 0 and r j is the jth dual node position of the mesh.Let us vary the spatial mesh and consider its effects on the solution.We employ three qualitatively different grids, which correspond to Delaunay meshes generated by the node positions illustrated in Fig. 1.The simplest and most commonly used grid is the one with the cubic tiling.Its popularity is mainly based on its ease of implementation.Second, we employ body-centered cubic (BCC) tiling [54,55], which is preferred by certain numerical studies [46,56].The third option is the C15 structure, which is one of the tetrahedrally closepacked tilings [57][58][59][60].The C15 structure has been found to be a high-quality grid for the solution of the Maxwell equations [47,53].For each of these three grid types, we employ three discretization levels, where tasks are scaled to involve 10 9 , 10 10 , or 10 11 floating point multiplications for integration over a unit time interval.

C15 BCC cubic
During the integration, we monitor the deviation S(t) = 1 − Ψ * λ,ḡ,κ (r, 0) Ψ(r, t) d 3 r from the stationary state and terminate the simulation when S(t) exceeds 0.1.The duration before the termination is referred to as the time span of stability.The evolution of S(t) is illustrated in Fig. 2. The time span of stability appears to be very sensitive to the grid type used.The BCC grid offers the longest time spans, since it is numerically the most isotropic of the three grids [53].With the finest discretization level, BCC leads to threefold splitting, which is the most likely physical solution for the used parameter values (see Sec. III B).In other cases, the fourfold symmetry of the cubic base grid steers the numerical solution into fourfold splitting.This demonstrates the importance of the tiling in obtaining correct physical results.
The BCC grid also offers the smallest early-stage errors before the actual vortex splitting occurs.The early-stage error seems to approximately obey the function h 4 , where h is the dual edge length.With the lowest discretization level (10 9  operations/unit time), the average dual edge lengths are 0.20, 0.16, and 0.17 for the cubic, BCC, and C15 grids, respectively.The edge lengths of the finest (10 11 ) and second finest (10 10 ) discretization levels are about 0.38 and 0.61 times the abovementioned edge lengths, respectively.
Owing to these results, we choose to employ the BCC grid in the remaining numerical simulations presented in this work.

B. Dominant splitting symmetries
Even the smallest random perturbation to a dynamically unstable stationary vortex state triggers the splitting of the vortex.To find the most likely physical splitting symmetries, the stationary vortex states are perturbed slightly by adding low-amplitude random noise in the beginning of the compu- tation.The discrete order parameter is initialized at instances k = −1, 0 by where ρ j is a random variable chosen uniformly from the unit disk in the complex plane.
The spatial discretization employs the BCC grid, whose dual edge lengths are < 5% of the effective wavelength This corresponds to the second finest discretization level of Sec.III A. The computational domain is a rectangle that contains all points r for which | Ψλ,ḡ,κ (r, 0)| is greater than 10 −5 times its maximum.Zero particle density is employed as the boundary condition.
The following procedure is applied to find dominant splitting symmetries.During a time integration, splitting indicators P l (t) = e ilφ | Ψ(r, t)| 2 d 3 r are computed at each time instance t.The number l dom ∈ N + , for which P l dom (t) ≥ P l (t), ∀l ∈ N + , indicates the dominant splitting symmetry.Vortex dynamics is divided into three categories: If P l dom exceeds 0.1 before the time reaches 200, we classify the case as vortex splitting with l dom -fold symmetry (see Fig. 3).Otherwise, if S(t) < 0.1 for the entire integration interval 0 ≤ t ≤ 200, we detect a relatively stable vortex and label this case as no split.Otherwise, we observe an unstable vortex without any obvious dominant splitting symmetry; this case is called unclear.
Three representative trapping ratios λ are employed to simulate oblate (λ = 10), spherical (λ = 1), and prolate (λ = 0.1) condensates.In addition, we vary the effective interaction strength ḡ and the winding number κ to obtain a comprehensive understanding of the splitting process.The observations from the time integrator are not entirely unique, since the results depend slightly on the seed of the random number generator.To reduce variation, we simulate each splitting process twice with different seeds and choose the splitting symmetry that is closer to the prediction of the Bogoliubov stability analysis.The splitting symmetry predicted by the Bogoliubov equation is defined as the one corresponding to the value of |l| for which max q Im(ω q,l ) is largest.Visual inspection of Fig. 4 shows that the results of the time integration mostly coincide with the predictions of the Bogoliubov equation.
The characteristics of the splitting symmetries as functions of ḡ and κ are similar for different aspect ratios.With lower  The symbol indicates the result of the time integration, while the background color corresponds to the prediction of the Bogoliubov equation, namely, the value of |l| for which max q Im(ω q,l ) is largest.aspect ratios, a given splitting symmetry is found at higher interaction strength, which is explained by the increased size of the condensate.The most significant difference is that the unclear splitting symmetries appear only in prolate and spherical condensates.This phenomenon will be studied in more detail in the next section.

C. Intertwining of vortices
In prolate condensates, we observe vortices to intertwine as they split, as illustrated in Fig. 5. Similar intertwining processes of doubly quantized (κ = 2) vortices have already been discovered in Refs.[18,20,21].Our study demonstrates that intertwining also occurs for large winding numbers.The branched intertwining of a five-quantum vortex (κ = 5) is illustrated in Fig. 5  The intertwining of vortices does not occur in the oblate condensates with the aspect ratio λ = 10, but the phenomenon seems to become observable when λ is close to 1.To investigate this further, we consider the dynamics of three-quantum (κ = 3) vortices for different aspect ratios.To equalize the local peak interaction strengths, the effective interaction strength ḡ is chosen to be inversely proportional to the aspect ratio as ḡ = 1000/λ.
The simulations indicate that the vortices in the oblate condensates of λ ≥ 1.5 are stable.In the prolate condensates with λ ≤ 0.5, the vortices seem to be unstable and exhibit intertwining.In between the oblate and the prolate, no prevalent behavior of the vortices is detected.Nevertheless, in a condensate with λ = 1.0, we discover a cyclic splitting process, where the vortex begins to split but then returns nearly to its initial state.This split-and-revival effect is illustrated in Fig. 6.

D. Computational performance
The time integrations of this paper were executed on central processing units (CPUs), but we have also implemented the solver with graphics processing units (GPUs).The performances of the two implementations are studied here by measuring the simulation times in the case λ = 0.1, ḡ = 5000, and κ = 20.We use up to 96 12-core Intel (Xeon) Haswell (E5-2690v3, 64bits) CPUs and up to four NVIDIA Tesla P100 GPUs.The results in Fig. 7 indicate that the performance of the GPU implementation on one GPU corresponds to the performance of the CPU implementation executed on at least 60 CPU cores.

IV. CONCLUSION
In summary, we have studied the splitting dynamics of giant vortices in dilute BECs with a particular focus on the time integration of the three-dimensional GPE.We showed that a significant reduction of the numerical error is achieved when a tetrahedral spatial tiling is utilized instead of the routine Cartesian grid.Importantly, the careful choice of the numerical method provides us with the physically correct splitting symmetry.
Comprehensive maps of vortex splitting symmetries were presented for oblate, spherical, and prolate BECs.The solutions of the time integrations were found to agree with the linear stability analysis based on the Bogoliubov equation.
The splitting-induced intertwining of vortices in prolate condensates is demonstrated.The aspect ratios for which the intertwining becomes observable are also studied.A splitand-revival phenomenon, where the vortex almost returns to its initial state after splitting temporarily, was observed in the  crossover from a dynamically stable vortex into an unstable one as a function of the aspect ratio.
The performance study presented in Sec.III D indicates nearly optimal scalability of the CPU implementation and promising performance for the GPU implementation.In the future, we will study how the GPU performance scales with a larger number of GPUs.This will allow us to accomplish even more challenging tasks than is currently possible with CPUs.These tasks may include solving the dynamics of a lattice of monopole-antimonopole pairs [61,62].

Figure 2 .
Figure 2. Error S(t) induced by the numerical implementation of the GPE as a function of time for different tilings and discretization levels.The parameters for the stationary state are λ = 1, ḡ = 300, and κ = 10.

Figure 4 .
Figure 4. Observed splitting symmetries in (a) the oblate (λ = 10), (b) spherical (λ = 1), and (c) prolate (λ = 0.1) condensates.The symbol indicates the result of the time integration, while the background color corresponds to the prediction of the Bogoliubov equation, namely, the value of |l| for which max q Im(ω q,l ) is largest. (b).

Figure 6 .
Figure 6.Effect of the aspect ratio λ on the stability of a threequantum vortex with the interaction strength set to ḡ = 1000/λ.(a) The deviation S(t) from the stationary state as a function of time.(b) Particle density isosurfaces visualizing the split-and-revival effect observed for λ = 1.0.

Figure 7 .
Figure 7. Performance of the CPU implementation (blue) and the GPU implementation (red) as a function of the computing resources.One iteration corresponds to the integration over one unit of time.