Vortex precession and exchange in a Bose-Einstein condensate

Vortices in a Bose-Einstein condensate are modelled as spontaneously symmetry breaking minimum energy solutions of the time dependent Gross-Pitaevskii equation, using the method of constrained optimization. In a non-rotating axially symmetric trap, the core of a single vortex precesses around the trap center and, at the same time, the phase of its wave function shifts at a constant rate. The precession velocity, the speed of phase shift, and the distance between the vortex core and the trap center, depend continuously on the value of the conserved angular momentum that is carried by the entire condensate. In the case of a symmetric pair of identical vortices, the precession engages an emergent gauge field in their relative coordinate, with a flux that is equal to the ratio between the precession and shift velocities.


Introduction
Bose-Einstein condensate is a widely investigated realization of a coherent macroscopic quantum state. In particular, dilute condensates of trapped ultra-cold atoms have unique quantum features that facilitate a high level of experimental control [1][2][3]. Various realizations are studied vigorously, both in earth-bound and in earth-orbiting laboratories [4]. Among the goals is the development of ultra-sensitive sensors and detectors [5], and the properties of cold atom condensates are also employed as a platform for quantum computation and simulation [6].
Quantum vortices are the principal topological excitations in a cold atom Bose-Einstein condensate [7][8][9][10]. A vortex is characterized by an integer valued circulation of the macroscopic wave function ψ(x) 1 2π d ·∇ arg[ψ] ∈ Z around its core. The time-dependent Gross-Pitaevskii equation [11,12] governs the dynamics of the macroscopic wave function, at the level of the mean field theory. This is a nonlinear Schrödinger equation with a quartic nonlinearity that accounts for interactions between the atoms [13][14][15][16][17]. Vortices that appear in rotating cold atom Bose-Einstein condensates have been extensively studied, both experimentally and theoretically [16,18,19]. At the level of the Gross-Pitaevskii equation these vortices are commonly modelled as stationary solutions in a co-rotating frame [13][14][15][16][17][20][21][22]. The rotating condensate accommodates growing angular velocity by forming vortices. This goes with a discontinuous increase in angular momentum: A rotating condensate does not support an arbitrary value of the angular momentum [20].
Here we are interested in solutions of the Gross-Pitaevskii equation in a trapped condensate that does not rotate. In that case, the equation has no time independent solutions besides the trivial one ψ(x) ≡ 0. However, we show that whenever the non-rotating condensate supports a non-vanishing arbitrary valued angular momentum, there is a stable minimum energy solution of the time dependent Gross-Pitaevskii equation. These solutions describe eccentric vortices that precess around the center of the non-rotating axially symmetric trap, as illustrated in Fig. 1. In particular, the single vortex solutions are very similar to the experimentally observed precessing vortices reported in [23].
It is notable that even though the solutions that we present are minima of the ensuing free energy, they are not critical points of the free energy. Thus we employ methods of constrained optimization [24,25], to solve the time dependent Gross-Pitaevskii equation.
We note that the microscopic, individual atom angular momentum is certainly quantized. But as pointed out in [26,27] (see also [20]), in the limit where the number of atoms is very large, the condensate can accommodate arbitrary values of the macroscopic angular momentum that appears as a Noether change in the Gross-Pitaevskii equation. It was pointed out in [26,27] and confirmed by the variational analysis in [27], in the case of a single vortex, that a solution with arbitrary angular momentum must in general be time dependent. Our results confirm this proposal.
We emphasize the following difference between our approach that considers a nonrotating trap and builds on the ideas of [26,27], and the more common approach that considers vortices in a rotating trap, in terms of a co-rotating frame [13][14][15][16][17][20][21][22]: In the rotating trap problem, the minimal energy state is controlled by the value of the angular velocity of the trap that appears as a free parameter in the stationary, time independent co-rotating Gross-Pitaevskii energy functional. The vortices appear as critical points of this energy functional, as solutions of the time independent Gross-Pitaevskii equation they minimize an unconstrained variational problem. In contrast, in our case we use the fact that the angular momentum is conserved, thus we fix its value. The vortices are then solutions to a constrained optimization problem and solve the time dependent Gross-Pitaevskii equation. The difference between these two problems is detailed in Appendix B.
Since the vortex configurations that we find are simultaneously both constrained minima of the Gross-Pitaevskii free energy, and time dependent solutions of the ensuing Hamilton's equation, they bear a certain resemblance to the concept of a classical Hamiltonian time crystal [28,29]. Our approach is modelled on the general theory developed in [30]; see also [31].
We shall further observe that in the case of a symmetric pair of two identical vortices, the evolution of the minimal energy vortex states, at specified angular momentum, brings about an exchange of the vortices. The exchange of identical bodies in two dimensions have been studied extensively, as it can be used to identify anyons. Anyons are two-dimensional quasiparticles that are neither fermions nor bosons [32,33]. The anyon statistics implies that the exchange of the position of two identical particles results in a change of their quantum mechanical phase which differs from 0 (the case of bosons) or π (the case of fermions); for reviews on anyons see e.g. [34,35]). The anyons, that can play an important role for quantum computing [36], have recently been observed by electronic interferometry in quantum Hall effect experiments [37].

Theoretical developments
We consider a two dimensional single-component Bose-Einstein condensate on the xy-plane with an axially symmetric, non-rotating, harmonic trap. This approximates for example an anisotropic three dimensional trap, resulting in an oblate condensate. Following [14,38] we assume that the condensate describes a dilute limit of a large number of atoms, with a two-body repulsive interaction potential between the atoms. The condensate is described by a macroscopic wave function ψ(x, t) whose dynamics obeys the (dimensionless) time dependent Gross-Pitaevskii equation [14,38]: Here g is a dimensionless coupling that accounts for the interactions between the atoms. In a typical Bose-Einstein condensate there are N a ∼ 10 4 − 10 6 ultracold alkali atoms and the corresponding coupling is typically g ∼ 1-10 4 ; for discussion of the nondimensionalization, see, e.g, [17] and the Appendix A. The free energy F in eq. (2.1) is We observe that F is a strictly convex functional, it can not have any critical points besides the global minimum ψ ≡ 0. Thus there are no nontrivial time independent solutions of the Gross-Pitaevskii equation (2.1), only time dependent nontrivial solutions of (2.1) can exist. 1 Provided appropriate conditions can be introduced, these solutions can also minimize the free energy F in the framework of constrained optimization. Thus, we search for time dependent solutions of (2.1) that are also minima of the free energy (2.2), where the latter is subject to an appropriate condition.
Indeed, besides the free energy (2.2) the time evolution equation (2.1) supports two additional conserved quantities, as Noether charges: The macroscopic angular momentum along the z-axis and the norm of the macroscopic wave function a.k.a. the number of particles Accordingly, we search for solutions of the time dependent Gross-Pitaevskii equation (2.1), with fixed values of both the angular momentum L z and the number of particles; in the following, with no loss of generality, we normalize the wave function so that N = 1 (see Appendix A for details). On the other hand, the angular momentum (2.3), can assume arbitrary values L z = l z , and even in the case of a single vortex the values l z of the angular momentum are not expected to be integer valued [26,27]. The Lagrange multiplier theorem [39] states that with fixed angular momentum (2.3) and fixed number of particles (2.4), the minimum of F is also a critical point of Here λ z , λ N are the Lagrange multipliers that enforce the values L z = l z and N = 1, respectively. The critical points (ψ cr , λ cr z , λ cr N ) of F λ (2.5) obey Let ψ min (x) be a solution of (2.6) that is also a minimum of (2.2), and let λ min N and λ min z denote the corresponding solutions for the Lagrange multipliers. Whenever the solutions for the Lagrange multipliers do not vanish, ψ min (x) cannot be a critical point of the free energy F in (2.2). Instead, it generates a time dependent solution ψ(x, t) of (2.1) which, from (2.6), with the initial condition ψ(x, t = 0) = ψ min (x) obeys the linear time evolution where both λ min N and λ min z are time independent [30]. Whenever either of the Lagrange multipliers is non-vanishing, the minimum energy solution is time dependent. But the time evolution amounts to a simple phase multiplication only when ψ min (x) is an eigenstate of the angular momentum operator; below we show that this can only occur for l z = 1, for all other values of l z = 0 the minimum energy solution always describes a precessing vortex configuration.
The minimum energy wave function ψ min (x) spontaneously breaks both Noether symmetries. The time evolution equation (2.7) describes a symmetry transformation of ψ min (x), that is generated by a definite linear combination of the two conserved Noether charges: where { , } is the Poisson bracket. This is an example of spontaneous symmetry breaking, but now the symmetry breaking minimum energy configuration is time dependent [30].

Numerical simulations
To search for minimum energy solutions of the time dependent Gross-Pitaevskii equation, we numerically minimize the free energy (2.2) subject to the conditions L z = l z and N = 1. The problem is discretized within a finite-element framework [40], and the constrained optimization problem is then solved using the Augmented Lagrangian Method (see details of the numerical methods in Appendix C). Minimal energy states for angular momentum values l z ∈ (0, 2] with N = 1, and for three representative values of the interaction coupling g = 5, g = 100 and g = 400 are displayed in figure 2; see also additional animations in Supplementary materials [41]. For non-vanishing values of the angular momentum 0 < l z < 1 the minimum energy configuration ψ min (x) is an eccentric vortex that precesses around the trap center (see regime A). As l z increases the vortex core approaches the trap center. This kind of eccentric precessing vortices have been previously proposed theoretically in [27]. They appear very similar to the precessing vortices that have been observed experimentally [23].
As shown in regime B, when l z = 1, the core position coincides with the center of the trap, and the Lagrange multipliers feature a discontinuity. At this value of l z , our solution coincides with the single vortex solution that has been described extensively in the literature [13-17, 21, 22]. When l z becomes larger than one, a second vortex appears, and there are two qualitatively different two-vortex configurations. As l z increases, these are sequentially asymmetric (regime C) and then symmetric (regime D) with respect to the trap center. Because it depends on g, there is no universal value of l z that separates these two regimes. At higher values of angular momentum, l z ≈ 1.8 with the exact value depending on g, additional vortices start entering the condensate; this is the regime E in The panels on the top row display the density of the condensate |ψ| 2 for qualitatively different solutions with g = 400, obtained for different values of l z (zoomed to relevant data while the actual numerical domain is larger). At the vortex core the density |ψ| 2 vanishes and the phase circulation is 2π. The five regimes A -E are detailed in the text. figure 2. Remark that for g = 100 and g = 400 a third vortex moves towards the trap center as l z increases, while for g = 5 a pair of vortices enters. Further increase of the angular momentum introduces additional vortices in the condensate (not shown).
The case l z = 1, labelled by B in figure 2, is a special case. There, the core of the vortex coincides with the center of the trap. The corresponding ψ min (x) is an angular momentum eigenstate with eigenvalue l z = 1, and the solution of the time dependent Gross-Pitaevskii equation (2.7) is merely an overall phase rotation with no change in the core position. However, for generic l z the minimum energy configuration ψ min (x) is not an eigenstate of the angular momentum.
We have simulated the time evolution (2.1) of the minimum energy configuration ψ min (x) using a Crank-Nicolson algorithm. Simulation details together with animations of the time evolution can be found in the Appendix C. For all values of l z = 1 the vortex configuration precesses uniformly around the trap center. This uniform rotation is a consequence of the equation (2.7) which determines the time evolution of the minimal energy configurations, in terms of the two conserved charges.

Exchange of vortices
Finally, an interesting phenomenon occurs for the time evolution of a symmetric pair of identical vortices (corresponding to the regime D of figure 2). The time evolution of such a vortex pair, with g = 5, is displayed in figure 3. The vortex pair rigidly rotates around the trap center at a constant distance, with an angular velocity specified by λ min z ≡ −Ω.
In addition the phase of the wave function evolves with a rate that depends on λ min N . This coincides with the rate of change in the phase arg[ψ](t), measured at the center of the trap. Similarly, the rate of change of θ(t) that measures the rotation of the pair, equals λ min z .
Since the two Lagrange multipliers are in general not commensurate, it follows that when the time evolution exchanges the position of two identical vortices, the change in the phase of the wave function is not by an integer multiple of π.
The quantum mechanical exchange of two identical bodies, with ensuing phase change, is familiar from description of anyons [32,33,42]. The similarity to that scenario can be elaborated further, by following [42] to rewrite the time evolution (2.7) as follows: Here L = −ix ∧ ∇ is the canonical angular momentum, the parameter φ = λ min N /λ min z , and the flux of A = ∇ tan −1 x/y equals 2π around any counter-clockwise path that encircles the z-axis once. Notably, in the case of two identical vortices x coincides with their relative coordinate in a center-of-mass frame. Furthermore, the equation (4.1) is akin a time dependent linear Schrödinger equation that describes the evolution of a charged point particle, moving on the (x, y) plane under the influence of a unit strength magnetic vortex filament along the z-axis. The corresponding Hamiltonian comprises the z-component of the covariant angular momentum L cov z . We note that when φ = 0, the rotations around z-axis, i.e. the changes in the polar angle θ, are generated by the covariant angular momentum.
To find the solutions of (4.1) we first remove A by phase multiplication, and then proceed with separation of variables in terms of polar coordinates (r, θ). 2) The single valuedness of the wave function ψ(x, t) implies that X (r, θ) is generally not single valued as it obeys L z X n = (n + φ)X n . (4. 3) The solution of the time evolution (4.1) is thus Here the f n (r) form an orthogonal set of normalizable, real valued functions, they are selected so that (4.4) describes the energy minimizer ψ min (x) of (2.6), (2.1) at t = 0.

Summary and outlook
In summary, we have investigated the time dependent two-dimensional Gross-Pitaevskii equation, that models dynamics in an axially symmetric harmonic trap. The equation supports two Noether charges that correspond to the number of atoms and to the axial component of the macroscopic angular momentum, respectively. We have used the methodology of constrained optimization to construct solutions that minimize the Gross-Pitaevskii free energy, with a fixed number of atoms and as function of the macroscopic angular momentum to which we have assigned arbitrary numerical values. For generic angular momentum values, a minimum energy configuration always exists and it describes eccentric vortices that precess uniformly around the center of the trap. In particular, no time independent solution can exist in a non-rotating cylindrically symmetric trap. Moreover, whenever the solution describes two identical vortices, the time evolution exchanges their positions but the phase of their common wave function evolves in a nontrivial fashion.
We foresee that the vortex precession in a non-rotating condensate that we have described, can be realized and observed in a variety of laboratory experiments; the precessing vortices already reported in [23] appear to be very similar to the vortex precession that we have simulated. There are several conceivable experimental setups, where vortices can form in a non-rotating condensate. These include vortex formation by optical photon beams that carry a non-vanishing angular momentum, with appropriate topologies [43][44][45][46]. Additional experimental set-ups can be built on thermal quenching in combination with the Kibble-Zurek mechanism [47][48][49][50]. This mechanism describes how topological defects form during rapid symmetry breaking transitions. There are already several laboratory examples that are relevant to the scenario we have studied, including neutron irradiation in the case of superfluid 3 He [51,52], evaporative cooling in the case of ultracold 8 7Rb atoms [53,54].

Appendices A Nondimensionalization of the Gross-Pitaevskii equation
In the main body of the paper, the dynamics of a Bose-Einstein ultracold atoms is described by the dimensionless Gross-Pitaevskii equation. Here, we derive this dimensionless form from the original equation (see e.g. [17]). In the ultracold dilute regime, the Bose-Einstein condensate state of N a identical atoms is described by the macroscopic wave function Ψ, whose dynamics obey where N = |Ψ| 2 , and g 0 = 4π 2 a s m . (A.1b) Here, ∇ 2 ≡ ∇·∇ is the Laplace operator, m is the mass of the atoms, is the reduced Planck's constant, and a s is s-wave scattering length of the atoms. Note that for practical purposes, we introduce N to be the norm of the condensate, while N a is the actual number of atoms in the condensate. The harmonic oscillator trap potential reads as where ω x ≤ ω y ≤ ω z are the trapping frequencies in the different spatial directions. The dimensionless units are obtained by the following rescaling Hence, after dropping the˜, the dimensionless Gross-Pitaevskii equation is The dimensionless time dependent Gross-Pitaevskii equation (A.4), is obtained for a three dimensional harmonic trap. In the main body of the paper, we consider the case of a two Bose-Einstein condensate. The two dimensional problem fairly approximates a disk-shaped condensate with small height in z-direction, for the trap frequencies ω x ≈ ω y and ω z ω x . This results in a further rescaling of the coupling constant g 2d = g 3d ωz 2πωx . For detailed derivation see, e.g., [17].
The dimensionless coupling g is thus the parameter of the accounts for the interactions between atoms. In the numerical investigations we choose values g that are representative of a Bose-Einstein condensate of ultracold alkali atoms. The Table 1 shows typical experimental values of trapped ultracold atoms, and the corresponding estimates for the length scales and the dimensionless coupling g.

B Constrained minimization with fixed angular momentum versus minimization at constant angular velocity
In the main body of the paper, we emphasize that there is a difference between the common set-up that describes vortices in a rotating trap in terms of a free energy expressed in the corotating frame [13][14][15][16][17][20][21][22], and the present set-up where vortices occur in a non-rotating trap; the present set-up was first addressed in [26,27]. Here we clarify the difference, by a direct comparison of the two cases: In the case of a rotating trap, vortices are created by rotating the trap, with a constant angular velocity Ω. The free energy F Ω is expressed in the co-rotating frame, and the vortices are the critical points of the free energy. The vortex structures can then be investigated as a function of the externally applied angular velocity Ω. In contrast, in our case the trap does not rotate. Instead, we use the fact that the angular momentum L z of the condensate is a conserved quantity, and the condensate carries a corresponding angular momentum value L z = l z . The ensuing minimal free energy configuration is a solution of a constrained optimization problem that minimizes the free energy F , subject to the condition L z = l z . The two free energies are related as follows: Due to the constraint, the free energy minimizer is in general not a critical point of F . But it can be located using the Lagrange multiplier theorem [39], in the framework of constrained optimization [24,25]. Whenever the minimizer of the free energy F is not a critical point, it describes a solution of a time dependent Gross-Pitaevskii equation.
In addition, in both cases the norm of the condensate wave function (i.e. the particle number) needs to be accounted for. Here the goal is to clarify the difference between the rotating and non-rotating cases. Therefore, in our comparison we treat the constraint N = 1 the same way, both when searching for the minimum of F Ω and F , by using the Lagrange multiplier theorem. For the problem we consider, the Lagrange multiplier theorem [39] states that with fixed angular momentum and fixed number of particles N , the minimum of F is also a critical point of Here λ z , λ N are the Lagrange multipliers that enforce the values L z = l z and N = 1 respectively. On the other hand, for a rotating condensate, the minimum of F Ω is a critical point of the following free energy where λ N is the (only) Lagrange multiplier that enforces N = 1. Note that the present Lagrange multiplier λ N is commonly treated as a chemical potential that is an externally controlled parameter, in par with Ω. This is not the case here, as λ N has to be adjusted to satisfy the constraint. The critical points of F λ and F Ω λ obey respectively and the conditions L z = l z and N = 1 , and the condition N = 1 .
Note that in a numerical simulation Ω is an input parameter that is fixed during the simulation, while λ z is adjusted sequentially, to satisfy the constraint (according to the Augmented Lagrangian Method, see details in Sec. C.2). The differences between the two approaches are summarized at the end of this section in Table 2.  Table 2. Summary of the differences: The red font shows the parameters that are given as input to the corresponding minimization problem. The purple font highlights the Lagrange multipliers. In a numerical simulation, these are sequentially adjusted according the Augmented Lagrangian Algorithm, until the constraints are satisfied (and the minima of F λ are found). The olive color fonts show the output quantities obtained in the numerical simulations.

C Details of the numerical methods
In the numerical investigations in the main body of the paper, we use Finite-Element Methods (FEM) (see e.g. [55,56]) both for the minimization and for the time-evolution. In practice we use the finite-element framework provided by the FreeFEM library [40]. Within the finite-element framework, the constrained minimization is addressed using an Augmented Lagrangian Method together within a non-linear conjugate gradient algorithm. The time-dependent Gross-Pitaevskii equation, is integrated using a Crank-Nicolson algorithm with a forward extrapolation of the nonlinear term.

C.1 Finite-element formulation
We consider the domain Ω which a bounded open subset of R 2 and denote ∂Ω its boundary. H(Ω) stands for the Hilbert space, such that a function belonging to H(Ω), and its weak derivatives have a finite L 2 -norm. Furthermore, H(Ω) = {u + iv | u, v ∈ H(Ω)} denotes the Hilbert spaces of complex-valued functions. The Hilbert spaces of real-and complex-valued functions are equipped with the inner products ·, · , defined as: The spatial domain Ω is discretized as a mesh of triangles using for the Delaunay-Voronoi algorithm, and the regular partition T h of Ω refers to the family of the triangles that compose the mesh. Given a spatial discretization, the functions are approximated to belong to a finite-element space whose properties correspond to the details of the Hilbert spaces to which the functions belong. We define P (2) h as the 2-nd order Lagrange finite-element subspace of H(Ω), and correspondingly P (2) h for H(Ω). Now, the physical degrees of freedom can be discretized in their finite element subspaces. And we define the finite-element description of the degrees of freedom as ψ → ψ (h) ∈ P (2) h . This describes a linear vector space of finite dimension, for which a basis can be found. The canonical basis consists of the shape functions φ k (x), and thus Here M is the dimension of V h (the number of vertices), the w k are called the degrees of freedom of w and M the number of the degrees of freedom. To summarize, a given function is approximated as its decomposition: w(x) = M k=1 w k φ k (x), on a given basis of shape functions φ k (x) of the polynomial functions P (2) for the triangle T i k . The finite element space V h (T h , P (2) ) hence denotes the space of continuous, piecewise quadratic functions of x, y on each triangle of T h .

C.2 Constrained minimization: Augmented Lagrangian Method
In the main body of the paper, we aim to minimize the free energy while enforcing a set of two conditions. Such problems, referred to as constrained optimization are studied in great details, see for example textbooks [24,25,57,58]. Here we describe the numerical algorithms that were used in the main body of the paper, to solve the constrained optimization problem. The Augmented Lagrangian Method (ALM) used to solve the constrained optimization is based on the following. In terms of the original energy functional to be minimized F , and the set of conditions C i , the augmented Lagrangian F aug is defined as Here µ is a penalty parameter and λ j are the Lagrange multipliers associated with the conditions C j . In the Augmented Lagrangian Method, the augmented Lagrangian is minimized, and the penalty µ is iteratively increased while the multipliers λ j ← λ j + µC j until all conditions are satisfied with a specified accuracy. The ALM algorithm has the property to converge in a finite number of iterations. Our choice for the minimization algorithm within each ALM iteration is a non-linear conjugate gradient algorithm [59][60][61].
In this work, the minimized functional F is the dimensionless Gross-Pitaevskii free energy [17] where g is the dimensionless coupling that accounts for the interatomic interactions, and V (x) is the trapping potential. In the main body we considered an harmonic oscillator trapping potential V (x) = |x| 2 /2. Since the derivations here do not depend on the specific form of the potential, we keep in general in the following. The two conditions specifying the particle number N = 1 and the value l z of the angular momentum are Hence, the variations of the augmented Lagrangian with respect to ψ and λ i give the set of equations where the angular momentum operator is L = −ix ∧ ∇. The associated weak form, is obtained by multiplying the equation (C.6a) by test functions ψ w ∈ H(Ω), integrating over Ω and integrating by parts the Laplace operator. Alternatively the weak form is calculated as the Fréchet derivative of the augmented Lagrangian as (C.3). In terms of the inner products (C.1), we find This formulation, hence corresponds to the variations of the Lagrangian with respect to all degrees of freedom. It can be seen as the gradient of a function to be minimized. To do so we choose a nonlinear conjugate gradient algorithm.

Error estimates
To estimate the quality of the solution, a global error can be derive by multiplying (C.6a) by ψ and integrating over the domain Ω (and again integrating by part the Laplace operator). The error is alternatively obtained by replacing the test functions ψ w by ψ in (C.7): In pratice during our constrained minimization simulations, we obtain the typical values for the relative error to be around 10 −5 . Numerical minimization of the weak formulation (C.7) of the augmented Lagrangian (C.3) gives the minimal energy states under the specified values of the constraints (C.5). In all generality we consider unit particle number N = 1 for various values of the angular momentum l z . The corresponding values of the Lagrange multipliers are displayed in figure 5.

C.3 Time evolution: forward extrapolated Crank-Nicolson algorithm
To address the question of the time-dependent Gross-Pitaevskii equation, the strategy is to use a Crank-Nicolson algorithm [62] to iterate the time series. More precisely, for efficient calculations, we write a semi-implicit scheme where the nonlinear part is linearized using a forward Richardson extrapolation. The time-dependent Gross-Pitaevskii reads as: (C.9) Figure 5.
Details of the results of simulations where the particle number is N = 1. The dependence of the Lagrange multipliers on the constrained angular momentum l z are displayed on the leftmost panel. The right panels illustrate that the time-evolution algorithm (C.15) used here indeed preserves the energy (bottom), the angular momentum (top) and the particle number (not shown).
The weak from, obtained by multiplying by test functions ψ w ∈ H(Ω) and integrating by parts reads as, in terms of the inner products (C.1), ψ w , i∂ t ψ = 1 2 ∇ψ w , ∇ψ + ψ w , V (x) + g|ψ| 2 ψ . (C.10) The discretization of the time turns the continuous evolution into a recursive series over the uniform partition {t} N n=0 of the time variable. The time variable is discretized as t = n∆t and the wave function at the step n is ψ n := ψ(n∆t). The Crank-Nicolson scheme, uses (Forward-)Euler definition of the time derivative, while the r.h.s of Eq. (C.10) is evaluated at the averaged times. Introducing the notations ∂ t ψ := δ t ψ n = ψ n+1 − ψ n ∆t , andψ n = ψ n+1 + ψ n 2 , (C.11) the Crank-Nicolson scheme for the time-dependent Gross-Pitaevskii equation (C.10) reads as: ψ w , iδ t ψ n = 1 2 ∇ψ w , ∇ψ n + ψ w , V (x)ψ n + ψ w , g|ψ n | 2ψ n . (C.12) This fully implicit scheme results in a nonlinear algebraic system which is very demanding to solve. The alternative to solving the nonlinear algebraic problem is to approximate the nonlinear part in terms of values at previous time steps. Very schematically, the idea of the modified algorithm is to approximate the fields in the nonlinear term by using an extrapolation of the previous time steps, that retains the same order of truncation error as the rest of time series. Thus, using the forward extrapolation the averaged wave function in the non-linear term becomesψ n ≈ (3ψ n − ψ n−1 )/2. Next, defining the time-discretized operators: U n ψ = ψ w , g 8 |3ψ n − ψ n−1 | 2 ψ , (C.13c) allows to rewrite Eq. (C.12) in the compact form O 1 ψ n+1 = O 2 ψ n + U n (3ψ n − ψ n−1 ) . (C.14) Hence the time-evolution is formally given by the recursion The spatial discretization is achieved by replacing the wave function ψ by its finite-element space representation ψ (h) ∈ V h (T h , P (2) ) in the time-dependent Gross-Pitaevskii equation (C.15). The test functions ψ w now take values in the same discrete space as ψ (h) . Denoting the matrix representation of the time-discretized evolution operators (C.13), as: the recursion Eq. (C.15) reduce to a linear system that read as: The vector R ψ which is a function of ψ which requires recalculating the vectors R ψ at each step and then multiplying by the inverse matrices. The numerical simulations of the time-evolution algorithm (C.18) representing the discretized evolution scheme (C.15) accurately reproduces the intrinsic physical properties of the time-dependent Gross-Pitaevskii equation. Namely it preserves the conserved quantities like the energy, the angular momentum, and the particle number. A consistency check, as displayed on figure 5, shows that these are indeed (exactly) conserved during the time-evolution.