A method for the dynamics of vortices in a Bose-Einstein condensate: analytical equations of the trajectories of phase singularities

We present a method to study the dynamics of a quasi-two dimensional Bose-Einstein condensate which contains initially many vortices at arbitrary locations. We present first the analytical solution of the dynamics in a homogeneous medium and in a parabolic trap for the ideal non-interacting case. For the homogeneous case this was introduced in the context of photonics. Here we discuss this case in the context of Bose-Einstein condensates and extend the analytical solution to the trapped case, for the first time. This linear case allows one to obtain the trajectories of the position of phase singularities present in the initial condensate along with time. Also, it allows one to predict some quantities of interest, such as the time at which a vortex and an antivortex contained in the initial condensate will merge. Secondly, the method is complemented with numerical simulations of the non-linear case. We use a numerical split-step simulation of the non-linear Gross-Pitaevskii equation to determine how these trajectories and quantities of interest are changed by the presence of interactions. We illustrate the method with several simple cases of interest both in the homogeneous and parabolically trapped systems.

Mathematically, a vortex that appears in complex waves is always associated to a phase singularity. A phase singularity occurs in those positions where the intensity of the complex field is zero and the phase is undetermined. Following circuits around these dark spots which are infinitely close to it, the phase increases or diminishes in integer multiples of 2π (when it increases (diminishes) one says the vortex has a positive (negative) charge). These phase singularities appear in many fields, such as plasma physics [33], fluid physics [34], atmospheric studies [35] or photonics [36] (with an independent branch of optics called singular optics [37][38][39]).
In this paper, we will use a method initially intro- * garciamarch@mat.upv.es duced in the context of photonics [40,41] as a benchmark to study the propagation of initial states containing many vortices with a Gross-Pitaevskii equation. This method has two sides: (1) it allows one to solve a linear Schrödinger equation (a paraxial scalar wave equation in the context of photonics or an ideal Gross-Pitaevskii equation (GPE) with a vanishing coupling constant) when the initial condition contains many phase singularities; (2) it offers naturally equations for the dynamical evolution of the location of each phase singularity, which we call the trajectory of the phase singularity or vortex line. From these equations one can obtain figures of merit of interest, such as merging points for singularities of opposite charge. While in the homogeneous case, the extension from photonics to BECs is merely a pedagogical analogy, in the inhomogeneous (parabolically trapped) case, we generalize here the technique to include potentials with powers of the spatial coordinates. We note that the parabollic case is not common in photonics (except for graded-index optical waveguides [42,43]). Once the technique for the ideal non-interacting case is established, we use the results, both in the homogeneous and trapped cases, to obtain some information insightful to analyze the non-linear cases. To illustrate this, we offer several examples of initial conditions (with one, two and three initial singularities), numerically solving the GPE with a split-step method in the non-linear case, and comparing the results with the ones obtained in the linear case.
The study of the motion of vortices in inhomogeneous (trapped) BECs has an extensive and thorough literature. For a single vortex in a two-dimensional condensate, early numerical studies illustrated that a single (off-axis) vortex precesses around the center of the trap, and the motion was described with effective models from which one can derive a Magnus force, which altogether allowed one to estimate the frequency of the precession [44].
With a variational approach, the effective potential experienced by the vortex was found, proving that indeed the vortex should precess around the cortex core [45], within the regime in which the (inhomogeneous) condensate is large compared with the size of the vortex core (also see Refs. [46,47] for a general study in two and three dimensions and anisotropic harmonic traps, where precession is also discussed, additionally finding, for the case of rotating trap, the angular frequency which stabilizes the vortex). Several works followed, researching diverse aspects of inhomogeneous BECs, generally with a Thomas-Fermi profile, and extending the limits of validity of previous studies or looking to other configurations or systems, e.g., two-component BECs (a non-comprenhesive list includes [48][49][50][51][52][53][54][55][56][57][58][59][60]-see also experimental results in Refs. [61][62][63][64][65]). The research effort also focused in the study of the dynamics of two vortices in a BEC, named also vortex dipoles, or more vortices [66][67][68][69][70][71], researching diverse aspects such as turbulence, chaos, or vortex solitonic structures [72][73][74][75] (see also review in [27]). These two or more vortices configurations were studied experimentally in [64,65,76,77]. We highlight here the case of Jones-Roberts solitons, which show elongated elliptical shape, are immune to the snaking instability, and can sustain imprinting of configurations of vortices [78,79]. Also, we note here that a Magnus force also appeared in some of the co-authors previous work, where it was considered a highly charged vortex located on axis and broken by a sudden turn on of an optical lattice (with a squared symmetry)-see [80].
We emphasize here that the present work diverts generally in one fundamental aspect from this list of works. The initial condition is not a vortex embedded in a BEC with Thomas-Fermi profile or a Jones-Roberts soliton with a phase pattern imprinted. In the examples used to illustrate the technique introduced here we consider one, two or three vortices inside a disc of density (see Figure  1). Generally, the technique is valid for an arbitrary number of vortices at arbitrary positions, see Equation (4). In all initial conditions, the vortices are located within the ring of density, and the dynamics of the singularities occurs mostly there. These initial conditions are more natural for interactions which are attractive, but here we also consider its dynamics in the repulsive case. As discussed in conclusions and outlook this paper has to be considered as a first step of a research program which will extend the method to other initial conditions and compare with the effective models developed in the literature. The paper is organized as follows: In Section II, we introduce the method both for the homogeneous and trapped (linear) cases. In Section III, we exemplify the method in the homogeneous case with some linear and non-linear examples (attractive and repulsive interactions, two and three initial singularities), and compare both. In Section IV, we undertake a similar analysis for the trapped case (here with one, two and three initial singularities). We end the paper with some conclusions and outlook in Section V.  1. (a-c) Amplitude (square root of density) of the initial condition for N = 5, 20 and 30 respectively, when there are initially two vortices of opposite charge (q = ±1), located at (±0.5, 0). (d) corresponding phase. The phase is the same in the three cases because, due to the form of Equation (4), one decides the position of the singularities with the coordinates of the initial vortices, and therefore the phase profile.
(e,f ) amplitude and phase, respectively, for three vortices, one negatively charged (q = −1) in the origin and two, positively charged (q = 1) located at (±1, 0). In all plots, the dashed black square is the window represented in other figures in the paper.. Also notice that the computational box is much larger than the one represented in all figures.

II. MODEL AND SYSTEM
The dynamics of a system of ultracold Bose-Einstein condensed bosons is governed by the Gross-Pitaevskii with r = (x, y), ∇ 2 the Laplacian, m the mass of the bosons, V (x, y) the external potential, and g the coupling constant. We assume a 2D system characterized by the complex wave function φ(r, t) of the condensed bosons, with dynamics frozen in the third direction z due to tight confinement, giving rise then to a quasi-two-dimensional Bose-Einstein condensate. In the following, we will use complex position operators,ŵ andŵ, which we define aŝ w =x + iŷ andŵ =x − iŷ. Their associated momenta arep = −i∂/∂ŵ andp = −i∂/∂ŵ. The commutations relations are From here on we omit the hats in the operators. These definitions allow one to write the GPE as Now the Laplacian is, explicitly, ∇ 2 w = ∂ 2 ∂w 2 + ∂ 2 ∂w 2 . We will consider initial conditions of the form with which is a Gaussian of amplitude determined by A and width σ (see some examples of initial conditions described by this equation in Figure 1). We note that in the case of BECs the amplitude to the square is the density of the condensate. We normalize the initial condition (4) to the total number of atoms, N . We keep N in the initial condition normalization to show explicitly how it appears in a quantity of interest like the merging point in the noninteracting case and because it determines the distance of the ring of density to the positions of the singularities (because we fix the position of singularities-see the examples for two singularities in Figure 1). With this form, the initial condition (4) is a combination of singularities embedded in a Gaussian wave function. Here, there are N A (N B ) singularities of charge +1 (-1) located at positions a i = (w i ,w i ) (b j = (w j ,w j )). The topological charge q of one singularity is defined as the circulation of the gradient of the phase of the complex field φ(w,w) = |φ(w,w)|e iθ(w,w) in circuits around the positions of the phase singularities a i (b j ), and very close to them, divided by 2π, i.e., In plane words, the phase grows (decreases) 2π along a circuit around a i (b j ) assuming this circuit does not encircle any other singularity. For simplicity of notation we consider in this point initial conditions which include only individual singularities of charge q = ±1 (as in the initial condition (4)). To consider an initial condition embedding individual singularities with larger (in modulus) charge is possible: by including powers of factors, i.e., (w − a i ) k or (w − b j ) s , with k and s integers, we would consider singularities of charge k or −s. We note that the total winding number of the initial condition is the sum of all topological charges, as it has to be calculated with the same formula but in a circle that surrounds all singularities inside. We refer to [81] for a pedagogical discussion on symmetry, angular momentum, topological charge and winding number, with the same definitions that we use in this manuscript.
Let us now discuss a method to easily find the analytical solutions of the GPE (1) for multisingular initial conditions, Equation (4), in the linear case, g = 0. This method was introduced in the context of the paraxial wave equation in optics [40,41] (see also [82] for an application to a different system). In the following sections we will compare the solutions analytically obtained for the linear case, with some interesting examples calculated numerically for the full non-linear equation.
The first step is to expand the initial condition (4) as powers of w andw as follows φ(w,w, 0) = {n,n} t n,n w nwn φ 00 (w,w), where {n,n} is the set which includes all powers of the form w nwn after the expansion and t n,n are complex coefficients. Let us define the functions φ n,n (w,w) = w nwn φ 00 (w,w), which are defined by means of two quantum numbers, n andn.
To give a meaning of these quantum numbers, let us see that they relate to the angular momentum and radial nodes quantum numbers. We define = n −n and p = min(n,n).
Let us first consider that n ≥n. Then, one can write Equation (8) as which, in polar coordinates, is This is a scattering mode (as termed in [41]) labelled by two quantum numbers: the angular momentum quantum number = n −n > 0 and the radial quantum number p =n. For the case n <n one obtains that p = n and = n −n < 0, and one can write Equation (8) as which, in polar coordinates, is which is again a scattering mode. Then, we can expand the initial condition in terms of the scattering modes (which form a basis) as that is, as a linear combination of scattering modes with a vortex of charge at the origin. This function (which contains the same information as Equations (7) and (4)) we normalize to the number of atoms N . As a second step, let us show how to find the analytical expression for the evolution of the functions in this basis. First we notice that, in the complex variables, we can write the effective Hamiltonian leading to Equation (1) with g = 0 as H = pp + V (w,w). We will consider two cases: i) the homogeneous case V (w,w) = 0; and ii) that the atoms are trapped in a parabolic potential, V (w,w) = 1 2 mω 2 |w| 2 (in cartesian coordinates, ). The first case corresponds to the optical case discussed in [40,41], where the technique used in this paper was introduced for photonic systems.
Here we generalize this technique to the trapped case and compare with the non-linear case, which is the most relevant for ultracold atom systems.
The evolution operator is U (t) = exp [itH/ ]. Using the commutation relations (2), one has that [w, pp] = ip and [w, pp] = ip. Also, we use that [w, ww] = 0 and [w, ww] = 0. Then, we obtain (We notice the The last step is to use this commutators to perform evolution. Consider an initial condition which is a scattering mode φ ,p (w,w, 0). To evolve it to time t we apply the evolution operator, φ ,p (w,w, t) = U (t)φ ,p (w,w, 0). We now use that the scattering mode can be written as (10) and (12)). Then, using the commutation relations (15) we obtain: if ≥ 0 and if < 0 In Equations (16) and (17) one can obtain analytically the evolved scattering mode φ ,p (w,w, t) if one can evolve the function φ 0,0 (w,w, 0). That is if one can obtain the function We call this function the generating function of the dynamics. For the homogeneous case, V (w,w) = 0, with an initial condition as (4) with ψ 0,0 in (5), the generating function is conventionally found in simple systems (like free particle in Schrödinger equation). For the inhomogeneous case, V (w,w) = 0, with the same initial condition, and for the particular case where the potential is parabolic, V (x, y) = 1 2 mω 2 (x 2 + y 2 ) that is V (w,w) = 1 2 mω 2 |w| 2 , one can also perform the evolution analytically, provided g = 0. We will use the two-dimensional harmonic oscillator propagator (see e.g., [80]) which we have written for time adimensionalized by the frequency ω of the harmonic potential and for space adimensionalized with the harmonic oscillator length a ho = /mω. Then, we can calculate the explicit expression for the generating function G(x, y, t) in cartesian coordinates where s(t) = σ 2 sin(t) + 2iπ cos(t). We note that this last expression is 2π-periodic. Now, to obtain any evolved scattering mode with Equations (16) and (17) one has to apply several operations to the generating function. To write this in a succinct manner let us define the raising and lowering operators: As shown in [41], these operators actually increase/decrease the angular momentum operator of a general scattering mode φ ,p by one unit. Also let us define the operator ∆ = + − which was also shown in [41] to increase its radial quantum number by one unit. Then, the evolution of a scattering mode can be written as follows: Let us summarize the method to solve the GPE, Equation (1) when g = 0 and the initial condition is a multisingular field of the form Equation (4). The method is as follows: 1. Expand the initial condition Equation (4), as a linear combination of terms w nwn of the form (7).
2. Use definition (9) to write this expansion as an expansion in terms of the scattering modes, Equation (14).
3. Evolve each scattering mode. To this end, find the generating function (which is to evolve the fundamental mode Φ 0,0 in this case).Then use Equation (22) to find each evolved scattering mode.
4. Use the evolved scattering mode to find the solution, via Equation (14), that is, Step 3 involves several repetitive operations on the generating function. This can be done analytically in the homogeneous case, V (w,w) = 0. The evolved scattering modes can be obtained systematically in this case, using the F -polynomials introduced in [41]. These give closed expressions for the evolved scattering modes. We remark here that the F -polynomials, which are tabulated in [41], obey recurrence relations which facilitate its computation. We also note that, in the inhomogeneous case, if the potential involves any combination of w andw and its powers, obtaining the generating function cannot be, in general, performed analytically. However, one can still use the method evolving φ 0,0 numerically and then, to evolve the scattering modes, perform the operations included in Equation (22) again numerically. Finally, one builds the solution φ(r, θ, t) as in step 4 (Equation (23)).
We finally mention that the second side of this method, both for the homogeneous and inhomogeneus cases in the non-interacting g = 0 case, is that it permits one to obtain algebraic equations for the singularity trajectories. These in turn allow one to determine several properties related to the position of the phase singularities of the evolved solution via solving these algebraic equations. To illustrate the method, we offer several examples in the next section. Also, we compare our results with the non-linear case. To this end, we use the split-step method to solve the Equation (1) in the non-linear (g = 0) homogeneous and inhomogeneous cases.

III. SOME EXAMPLES IN THE HOMOGENOUS SYSTEM
A. Two Initial Singularities, One Positive and One Negative First, we will discuss the case with a field with a positive singularity in a = (a, 0), and a negative one in We note that the same method can be used for singularities at arbitrary initial conditions, but the results obtained are long and we omit them here for the sake of simplicity. We remark here that in all calculations below, we use the adimensionalization mentioned above for the trapped case (see paragraph after Equation (19)), that is time is scaled by a trapping frequency ω and positions by the related harmonic oscillator length a ho ≡ /mω, with a generic mass and generic frequency. Then, by choosing mass and frequency, one can recover units.
We first expand the initial condition (24) in powers of w andw Following step 2 we produce And now, step 3 is accomplished using Equation (22), obtaining where q(t) = t + iσ 2 and γ(t) = πσ 2 tq(t) . Now the solution can be built from Equation (25), realizing step 4, Equation (28) is the solution of the evolution for all t, given g = 0. We compared this solution with those obtained with a numerical simulations of the GPE (1) for g = 0 obtaining complete agreement. The numerical simulations were performed with a split-step method, which is a spectral method used conventionally when solving non-linear Schrödinger equations, appearing in various fields, like non-linear optics (see e.g., [83], where code can be found in appendices). This method is valid for g = 0, for g > 0 (repulsive case) and g < 0 (attractive case). For all simulations we use a computational domain which is a box where x, y ∈ [−10, 10] in adimensional units. We use a discretization with M = 1024 points in each side (and then ∆x = ∆y ≈ 0.02 a.u. For the time step we use also ∆t = 0.02 a.u. We calculate for very long times which depend on the simulations. That is, the simulation time is limited for the repulsive and attractive case when the density reaches the boundaries, and starts to interfere with the central dynamics. This maximum computational time is almost in all (homogeneous) cases t max = 15 a.u. In the attractive case, in some instances the computational time is shorter than that because simulation stops when instability occurs. For the trapped case, we can perform much longer simulations. For every example shown in this paper, we also calculated with the split-step method in the non-interacting case to check that the results coincide with the analytical solution. We checked in all calculations that the energy is conserved during the whole simulation.
Once the linear solution is established, it is interesting to compare with the non-linear cases. We show a typical non-linear evolution in Figure 2, where we plot the amplitude and phase after some time evolution for the repulsive case, with N = 30 and g = 0.4. We present two exemplary times, at t = 1.5 a.u. and t = 2.5 a.u., which are a bit before and a bit after the merging time (the time at which the two singularities have annihilated each other). After merging time there is no singularity present in the phase profile. In Figure 3, we plot the same but for larger interactions strength, g = 0.6. Now for the same time as in previous example, t = 2.5 a.u., merging of the two singularities has not occurred yet.
For the attractive case, we plot in Figure 4 the amplitude and phase after some time evolution of different initial conditions, that is, for t = 2.2 a.u., g = −0.5 and N = 5, N = 20 and N = 30. First we see that the dynamics of the amplitude is very different in this case when compared to the repulsive one, as expected. In the repulsive one the amplitude spreads and widens, while in the attractive case it focuses more and more, leading eventually to instability. The case N = 30 shows an example just a bit before and after instability occurs (see panels (e) to (h)). Also, in view of Figure 4 we see that for this g, merging time is larger than t = 2.2 a.u. for N = 5 and shorter for N = 20 (there is no singularity in the profile for N = 20 yet there are two singularities for N = 5). For N = 30 there are two singularities because N has the effect to make the merging to take place at larger |g|, for large enough N (that is for this g merging time is larger for N = 30). In summary, the behavior is very different for different N : whilst for N = 5 and N = 30 we see the singularities has not merged yet, for N = 20 merging has occurred. Also for N = 30 we see an example just before instability takes place.
We see here that the trajectories followed by the singularities, alternatively referred to as vortex lines, can show intricate behavior. It is helpful to determine these trajectories for the linear case and the compare with the non-linear cases. In these non-linear cases, as we will see,  ((a,b)) and after t = 2.5 a.u. ((c,d)). Notice we plot a box which is smaller than the computational box. For t = 1.5 a.u. merging has not yet occurred but for t = 2.5, a.u. the two singularities have merged leaving a phase profile without singularities. The dashed black squares mark the plotting box in other figures below.
Here, merging has not occurred at t = 1.5 a.u. nor at t = 2.5 a.u., contrarily as in the case with g = 0.4. The dashed black squares mark the plotting box in other figures below. the trajectories change due to the effect of the density in the dynamics. It would be very interesting to study, in the context of the literature mentioned in the last paragraph of the introduction, the behavior in the non-linear case, and whether this method can shed some light in this direction. However, here we focus on introducing the method to find equations for the vortex lines in the non-interacting case, and use it to compare with the interacting case.
To obtain equations of the vortex lines we note that Equation (28) gives more information. Equating to zero We also plot in panels (g,h) the amplitude and phase for N = 30 and t = 4.0 a.u., a time in which instability has occurred and the simulation is not valid anymore.
this solution we can find in which points the solution is zero, for all t. This will provide equations for the trajectories of the singularities. For this particular case, the equation is which we solve together with its complex conjugate. From here we can find some figures of merit. In the case of two singularities, we expect that at certain merging time t m , the two singularities converge and annihilate each other. This merging time can be found analytically from Equation (29). For b = −a it is Notice that this expression depends on σ and a. Given a, for different numbers of atoms, changing σ one can implement the normalization of the initial condition, that is, the number of atoms N . In Figure 5 we plot how the merging time changes with the number of atoms in the initial condition, N . As shown, it is a decreasing function of N . Notice that, due to the way we wrote the initial condition, Equation (4), N governs the distance between the singularities and the ring of density surrounding them, and the shape of the density. The dependence of the initial conditions density shape with N is illustrated in Figure 1. In Figure 6a we plot the position of the two singularities as a function of time (the trajectories or vortex lines), for a = 0.5, b = −0.5, and N = 30 atoms (which gives σ 2 = 10), from the analytical solution, Equation (23). We notice that to find the position of the singularities requires a dedicated numerical method. For a large enough number of singularities it may require sophisticated methods, see, e.g., Ref. [84,85] but in our case we used a simple method. The analytically calculated merging time is t m = 1.9 a.u. In Figure 6b-d we plot the position of the two singularities as a function of time for repulsive interactions, when g = 0.4, 0.54 and g = 1, with N = 30. The vortex lines shown in Figure 6b, for g = 0.4 illustrate how merging time is increased with g. The case plotted in Figure 6c, calculated for g = 0.54, shows a case where the singularities do not merge and seem to stay parallel to the time axis for the whole computational time (here only plotted up to t = 4a.u. to facilitate location of merging time). For larger g, the singularities seem to travel outwards (e.g., as in Figure 6d, for g = 1). However, when calculating for longer times we observe that they bend inwards again. Nevertheless we cannot study longer dynamics, because of the presence of the computational boundary, so we cannot be conclusive about the long-time behavior of the singularities.
In Figure 6e,f we plot the position of the two singularities as a function of time for the non-interacting case and for attractive interactions, when g = −0.75, with N = 5. In the attractive case, for N = 30 instability occurs for very small interactions before merging time, so we calcu- lated for various cases with N < 20 (see also Figure 7). For N = 5 the merging time in the non-interacting case is t m = 3.4 a.u. (see Figure 5). As seen in Figure 6f for attractive interactions and N = 5 the merging time is decreased. Nevertheless, this is not always the case, as we discuss below.
For better visualization, we plot in Figure 7a the xprojection of the trajectories of the singularities for the same cases as in Figure 6a-d, that is, for N = 30 and for g = 0, 0.4, 0.54 and 1 (blue, magenta, red, green curves respectively). Here, we also see the behavior of the cases at g = 0.54 and g = 1, where the singularities seem to travel parallel or outwards, respectively. Nevertheless, as said before, we cannot conclude that they stay traveling parallel or outwards due to computational limitations. We plot in Figure 7b the merging time as a function of g for N = 5, N = 20 and N = 30 for the repulsive case. We see that for N = 30 the merging time increases only or g > 0.38 approximately while for N = 5 and N = 20 it increases for every g > 0.
In Figure 7c, we plot the x-axis coordinate of the singularities for the attractive case, with N = 5 and for g = 0, −0.25, and −0.75. As seen, the merging time decreases as |g| is increased. Nevertheless, this behavior changes as N is increased. In Figure 7d we plot the merging time as a function of |g| for N = 5, N = 10 and N = 20 for the attractive case. We see that for N = 5 (blue curve) the merging time is indeed reduced as the strength of interactions is increased. For N = 10 (green curve) instead it decreases slightly. For N = 20 (orange curve) the merging time increases with the strength of interactions. It is not the goal of this paper to study this effect in depth. We conjecture this is an effect of a nontrivial interplay which involves the dynamical behavior of the density.
Now we use Equation (22), to obtain the one wave function we do not have from previous example (34) We obtain the solution from Equation (23) Equation (35) represents the solution for all time in the linear case. We also solved numerically for interacting cases. In Figure 8 we present some illustrative cases, for N = 150. In panels (a) to (d) we present amplitude and phase for g = 0.3 at t = 1 a.u. and at t = 4 a.u. ((c) and (d)). As seen one positive singularity has merged with the central negative singularity, leaving in the origin a positive singularity. We plot in (e) to (h), amplitude and phase for g = 0.6 at the same times. Very interestingly, now we observe that at t = 4 a.u. merging has not yet occurred. Very remarkably, two pairs of positive/negative have appeared far from origin. We discuss this effect succinctly after obtaining the vortex lines.
Then, to look in some more depth the behavior of the phase singularities, let us study the trajectories in the linear case. From Equation (35) we find the equations for the trajectories of the singularities and the merging point. For this particular case the equation is To find the three zeros we need also to solve the complex conjugate of this equation. In the case that a 2 = −a 1 the evolution is such that at certain merging time t m one positive charged singularity merges with the central negatively charged one. The merging time is In Figure 9 we plot the position of the three singularities as a function of time (the vortex lines), for a 1 = 1, a 2 = −1, and N = 150 atoms (giving σ 2 ≈ 10.5) for (a) the linear g = 0 case and two exemplary interacting cases, (b), (c) and (d) with g = 0.3, g = 0.45 and g = 0.6, respectively. Again, we use the analytical method for the linear case and numerical split-step simulations both for the linear (to check it coincides with the analytical) and interacting cases. The analytically calculated merging time is now t m = 1.9 a.u., and again it corresponds with the numerically calculated one. As in previous example we observe that for increasing interactions, the merging time gets larger and larger. We observe a second effect which is that the singularities emerge again after certain time (see panel (b) for g = 0.3. We have calculated for a collection of values between g = 0.3 and g = 0.6. We observed that this re-emergence of the singularities occurs closer and closer to the merging time and eventually the singularities do not merge again in very long calculations (we calculated for values pf g up to 1). Also, we observed a third effect of the interactions, really remarkable: at certain points, vortex/antivortex pairs appear. They occur in areas of low density, but outside of the central ring. From Figure 8c, one see that for g = 0.3 the low density area occurs away from origing. This is the area where the singularities appear in Figure 8g,h for a larger interactions, g = 0.6. This effect is then not at all trivial, but it falls out of the scope of this paper. It will be the subject of future research. Finally, we mention that in this example we have fixed the positions of the vortices in the initial condition in the x axis. To consider arbitrary positions for the two positive singularities is possible, but it provides lengthly expressions that we do not reproduce here.  ((a,b)) and after t = 4 a.u. ((c,d)). In (d) we see that merging of two of the singularities has occurred, leaving only one on-axis singularity. (e) to (h), same for g = 0.6. Now, at t = 4 a.u. merging has not yet occurred. Also two pairs of singularities have appeared far from origin. The dashed black squares mark the plotting box in Figure 9.

IV. SOME EXAMPLES IN THE PARABOLICALLY TRAPPED SYSTEM
In this section we illustrate the dynamical evolution of initial conditions containing vortices when the external potential is a parabolic trap. The first example is a single singularity located outside the center of the potential trap, i.e., φ(w,w, 0) = (w − a)φ 00 (w,w). (38) In the second example the initial condition contains two singularities, both away from the center of the potential trap, one positively charged and one negatively charged, In the third example the initial condition contains three singularities, one negatively charged located in the minimum of the potential, and two away from its center, both positively charged, φ(w,w, 0) = (w − a)(w − a )wφ 00 (w,w). (40) To solve analytically in the linear g = 0 case, we use the method described in Section II, for the trapped case. The main difference with the homogeneous case is that one propagates φ 00 using propagator (19) to obtain Equation (20). In Figure 10 we present the dynamical evolution of the singularities for the three examples (we emphasize that in all cases we checked that energy is conserved). In Figure 10a we present the case of a single phase singularity, when a = (1, 0) and N = 12 (σ 2 = 10). The trajectory shows the influence of the trap, as it spins around approaching more the center of the trap. For the second example, shown in Figure 10b, the trajectories of the two phase singularities with opposite charge bend around and eventually merge and annihilate each other. From that time on, no singularity persists in the field. Here, a = (1, 0), b = (−1, 0) and N = 30 (σ 2 = 10). In the third example (Figure 10c), the positively charged singularities perform a spiral movement before reaching the center of the trap, where one of them annihilate with the central negatively-charged one. From that time on, there is only one phase singularity in the field which stays in the potential minimum of the trap. Here, a = (1.5, 0), a = (−1.5, 0) and N = 150 (σ 2 = 10). These three examples show the rich variety of dynamics one can explore with this method in a parabolic trap. In the panels above these figures we present these evolutions as three dimensional plots, for a complementary visualization.
To consider interactions shows similar effects as in the homogeneous case, that is, to contribute to avoiding annihilation of singularities or to create pairs, for instance. To illustrate this, we show in Figure 11 the evolution with g = 0 and 0.3 when a = (0.5, 0) and b = (−0.5, 0). This is similar to the second example (shown in Figure 10b), but we put the singularities initially closer to the minimum of the potential to make the merging time shorter. As shown in Figure 11a, the singularities tend to center of the trap as in Figure 10b, and annihilate each other at a merging time of t m = 1.5a.u., leaving a field without singularities from that time. In the interacting case, shown in Figure 11b, they do not annihilate at this merging time. Instead, the singularities repel, and perform complicated trajectories around the center of the trap from that time on. We also observe creation of pairs which eventually annihilate. The results for larger g (not shown) are qualitatively similar to the homogeneous case, that is, merging in the center of the trap is avoided and pairs are created. Here, the trajectories become very intricate and for that reason we do not include them. There is also a re-appearance of the singularities (see panels (b,c)). For large enough singularities, these singularities do not merge again for very long simulations (see panel (d)). There is also one more effect visible in panel (d): the generation of vortex-antivortex pairs. Green (red) tubes represent positively (negatively) charged singularities.

V. CONCLUSIONS AND OUTLOOK
We presented a method which allows one to study different aspects of the dynamics of a quasi-two dimensional Bose-Einstein condensate with a number of vortices located at arbitrary points. The method has an analytical part, which is only valid when one considers a vanishing coupling constant. We considered two possible realizations: a homogeneous system and a Bose-Einstein condensate in a parabolic trap. For the homogeneous system, the method is similar to that introduced in the context of photonics [41]. Here, we adapted it to the case of atoms. In this paper, we extended the method to the trapped case, which is more conventional in the context of Bose-Einstein condensates. This requires us to determine a generating function, which is very different than the one in the homogeneous case.
Once this first part is summarized in a recipe with four steps, we showed that the analytical solution can be used to determine the trajectories of the singularities (vortex lines), via a set of algebraic equations. These equations can be used to determine some figures of merit. Then, for the non-linear case, it is necessary to find the numerical solution of the Gross-Pitaevskii equation via a split-step simulation. The role of the analytical solution found in the first part is to determine a non-interacting scenario to compare with. The numerical simulation of the nonlinear case allows one to observe how interactions mod-ify the dynamics obtained in the non-interacting case. For small interactions and depending on the number of atoms, the dynamics does not change appreciably in some cases.
We illustrated the method with several examples. In the homogeneous system, we discussed the cases of two singularities with opposite charge and three singularities, one, located at the origin, with opposite charge than the other two. We found the equations for the trajectories of the singularities in the two cases. With this we found the merging time at which two singularities of opposite charge will collide and annihilate. We then turned on interactions, both attractive and repulsive, and found, numerically, how the trajectories change and how the merging time changes as a function of the number of particles and the coupling constant. For two singularities, we exemplified the kind of study one can perform numerically with this method. For repulsive interactions, we observed that the number of particles reduces the merging time for the non-interacting case. The merging time increases with the coupling constant, but for larger N the merging time is less sensitive to small values of the coupling constant. For N = 30 atoms we observed that the merging time does not change until g = 0.38. On the contrary, for larger g the behavior seems similar irrespectively of the number of atoms. For attractive interactions we see that the behavior is very affected by the presence of the instability. We show that there is a complex interplay with the density, which for small N makes the merging time decrease with the strength of the interactions and for larger N makes the merging time increase with the strength of the interactions. For three singularities, we observed qualitatively the effect of interactions (only repulsive) in merging time and an additional effect of the interactions, which is the spontaneous creation of vortex/antivortex pairs at certain times.
For the inhomogeneous, trapped system, we included examples for one, two oppositely charged singularities, and three singularities, the one in the origin with opposite charge than the other two. We showed, within the non-interacting scenario, that the trajectories are influenced by the trap, making the trajectories spiral-like in the trap. We illustrated the introduction of interactions in one case, only to show that similar effects as in the homogeneous case are observed, that is, modification of merging time and creation of vortex/antivortex pairs.
The results presented in this paper are only intended to illustrate the utility of the method. We envisage several future research directions. For example, to study systematically the dependence of different quantities of interest, not only merging time, with number of atoms and interactions, for different cases, such as those included here, both in the homogeneous and trapped cases. We also consider it interesting to study the phenomena of creation of vortex/antivortex pairs and its dependence with number of atoms, coupling constant, and trap frequency.
Another research problem which can be approached with this method is to study initial conditions with many singularities at random positions; we note that this last case will require us to determine its locations numerically (which as commented requires a dedicated numerical method as [84,85]). In addition, it will require us to solve a large set of equations even in the non-interacting case. The perturbative study and comparison with the analytical solution may give some hints as well. Last but not least, we see that one can pursue adapting the method to initial conditions which are more conventional in the context of BECs. That is, with a Thomas-Fermi profile and vortices separated by their vortex cores or Jones-Roberts solitons sustaining vortices as in [78,79]. The comparison with the effective models established in the literature may give interesting results. For the initial conditions studied here, it is still of interest to study the interplay of the singularities with the density, probably with similar models as those used for initial conditions with a Thomas-Fermi profile. All these possible directions show that this method can be of utility in the study of vortices in Bose-Einstein condensates.