# Vortex precession and exchange in a Bose-Einstein condensate

Julien Garaud, Jin Dai, Antti J. Niemi

► To cite this version:

Julien Garaud, Jin Dai, Antti J. Niemi. Vortex precession and exchange in a Bose-Einstein condensate. Journal of High Energy Physics, 2021, 2021 (7), pp.157. 10.1007/JHEP07(2021)157. hal-03887740

HAL Id: hal-03887740

https://hal.science/hal-03887740v1

Submitted on 22 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Vortex precession and exchange in a Bose-Einstein condensate

Julien Garaud<sup>a</sup> Jin Dai<sup>b</sup> Antti J. Niemi<sup>b,c,d</sup>

<sup>a</sup>Institut Denis Poisson CNRS-UMR 7013, Université de Tours, 37200 France

<sup>b</sup>Nordita, Stockholm University and Uppsala University,
Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden

<sup>c</sup>Pacific Quantum Center, Far Eastern Federal University,
690950 Sukhanova 8, Vladivostok, Russia

<sup>d</sup>Department of Physics, Beijing Institute of Technology, Haidian District, Beijing 100081, China

E-mail: garaud.phys@gmail.com, jin.dai@fysik.su.se, Antti.Niemi@su.se

**ABSTRACT:** 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.

**KEYWORDS:** Vortex dynamics, Bose-Einstein condensate
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Theoretical developments</td><td>3</td></tr><tr><td>3 Numerical simulations</td><td>5</td></tr><tr><td>4 Exchange of vortices</td><td>7</td></tr><tr><td>5 Summary and outlook</td><td>8</td></tr><tr><td>A Nondimensionalization of the Gross-Pitaevskii equation</td><td>9</td></tr><tr><td>B Constrained minimization with fixed angular momentum versus minimization at constant angular velocity</td><td>10</td></tr><tr><td>C Details of the numerical methods</td><td>13</td></tr><tr><td>C.1 Finite-element formulation</td><td>13</td></tr><tr><td>C.2 Constrained minimization: Augmented Lagrangian Method</td><td>14</td></tr><tr><td>C.3 Time evolution: forward extrapolated Crank-Nicolson algorithm</td><td>15</td></tr></table>

## 1 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–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–10]. A vortex is characterized by an integer valued circulation of the macroscopic wave function $\psi(\mathbf{x})$

$$ \frac{1}{2\pi} \oint d\ell \cdot \nabla \arg[\psi] \in \mathbb{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–17].
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Examples of asymmetric vortex configurations of a Bose-Einstein condensate in an axially symmetric harmonic trap. The elevation of the semitransparent surface stems for the condensate density $|\psi|^2$, while the coloring projected onto the *xy*-plane indicates the value of the phase arg[$\psi$] of the macroscopic wave function. The examples shown are minima of the two dimensional Gross-Pitaevskii free energy (2.2), with non-integer values of the angular momentum (2.3) $L_z = 0.80$ (left panel) and $L_z = 1.48$ (right panel). The vortices precess around the trap center in accordance with the time dependent Gross-Pitaevskii equation (2.1).</figcaption></figure>

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–17, 20–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 $\psi(\mathbf{x}) \equiv 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 non-rotating 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–17, 20–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 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 $\pi$ (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].

## 2 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 $\psi(\mathbf{x}, t)$ whose dynamics obeys the (dimensionless) time dependent Gross-Pitaevskii equation [14, 38]:

$$i\partial_t\psi = -\frac{1}{2}\nabla^2\psi + \frac{|\mathbf{x}|^2}{2}\psi + g|\psi|^2\psi \equiv \frac{\delta F}{\delta\psi^*}. \qquad (2.1)$$

Here $g$ is a dimensionless coupling that accounts for the interactions between the atoms. In a typical Bose-Einstein condensate there are $N_a \sim 10^4 - 10^6$ ultracold alkali atoms and the corresponding coupling is typically $g \sim 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

$$F = \int d^2x \left\{ \frac{1}{2} |\nabla \psi|^2 + \frac{|\mathbf{x}|^2}{2} |\psi|^2 + \frac{g}{2} |\psi|^4 \right\}. \qquad (2.2)$$

We observe that F is a strictly convex functional, it can not have any critical points besides the global minimum $\psi \equiv 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.<fnref n="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

$$L_z \equiv \int d^2x \, \psi^* \mathbf{e}_z \cdot (-i\hbar \mathbf{x} \wedge \nabla) \psi \qquad (2.3)$$

and the norm of the macroscopic wave function a.k.a. the number of particles

$$N \equiv \int d^2x \, \psi^* \psi. \qquad (2.4)$$

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

$$F_{\lambda} = F + \lambda_{z}(L_{z} - l_{z}) + \lambda_{N}(N - 1). \qquad (2.5)$$

Here $\lambda_z$, $\lambda_N$ are the Lagrange multipliers that enforce the values $L_z = l_z$ and $N = 1$, respectively. The critical points ($\psi_{cr}$, $\lambda_z^{cr}$, $\lambda_N^{cr}$) of $F_\lambda$ (2.5) obey

$$-\frac{1}{2} \nabla^2 \psi + \frac{|\mathbf{x}|^2}{2} \psi + g |\psi|^2 \psi = -\lambda_N \psi + i \lambda_z \mathbf{e}_z \cdot \mathbf{x} \wedge \nabla \psi \\ \text{together with } L_z = l_z \text{ and } N = 1. \qquad (2.6)$$

Let $\psi_{min}(\mathbf{x})$ be a solution of (2.6) that is also a minimum of (2.2), and let $\lambda_N^{min}$ and $\lambda_z^{min}$ denote the corresponding solutions for the Lagrange multipliers. Whenever the solutions for the Lagrange multipliers do not vanish, $\psi_{min}(\mathbf{x})$ cannot be a critical point of the free

<footnote n="1">¹In the case of a rotating condensate [13–17, 20–22] the angular velocity makes the energy into a non-convex functional with non-trivial critical points that describe vortices with a discrete valued angular momentum.</footnote>
energy F in (2.2). Instead, it generates a time dependent solution $\psi(\mathbf{x}, t)$ of (2.1) which, from (2.6), with the initial condition $\psi(\mathbf{x}, t = 0) = \psi_{min}(\mathbf{x})$ obeys the linear time evolution

$$i\partial_t\psi = -\lambda_N^{min}\psi + i\lambda_z^{min}e_z \cdot \mathbf{x} \wedge \nabla\psi, \qquad (2.7)$$

where both $\lambda_N^{min}$ and $\lambda_z^{min}$ 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 $\psi_{min}(\mathbf{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 \neq 0$ the minimum energy solution always describes a precessing vortex configuration.

The minimum energy wave function $\psi_{min}(\mathbf{x})$ spontaneously breaks both Noether symmetries. The time evolution equation (2.7) describes a symmetry transformation of $\psi_{min}(\mathbf{x})$, that is generated by a definite linear combination of the two conserved Noether charges:

$$i\partial_t\psi = -\lambda_N^{min}\frac{\delta N}{\delta\psi^*} - \lambda_z^{min}\frac{\delta L_z}{\delta\psi^*} \equiv -\lambda_N^{min}\{N, \psi\} - \lambda_z^{min}\{L_z, \psi\}, \qquad (2.8)$$

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].

## 3 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 \in (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 $\psi_{min}(\mathbf{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 \approx 1.8$ with the exact value depending on $g$, additional vortices start entering the condensate; this is the regime E in
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Minimum energy states for angular momenta $l_z \in (0, 2]$. The bottom row shows the dependence of the multipliers $\lambda_N^{min}$ and $\lambda_z^{min}$ on the angular momentum $l_z$ for $g = 5, 100$ and $400$. The panels on the top row display the density of the condensate $|\psi|^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 $|\psi|^2$ vanishes and the phase circulation is $2\pi$. The five regimes A – E are detailed in the text.</figcaption></figure>

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 $\psi_{min}(\mathbf{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 $\psi_{min}(\mathbf{x})$ is not an eigenstate of the angular momentum.

We have simulated the time evolution (2.1) of the minimum energy configuration $\psi_{min}(\mathbf{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 \neq 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.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The evolution of a minimum energy configuration $\psi_{min}(\mathbf{x})$, in the case of a symmetric vortex pair (with $g = 5$) in region D of figure 2. The bottom panels show the evolution of $\theta$, the rotation angle of the pair, and of the value of the phase $\arg[\psi]$ measured at the center of the trap. The three panels A–C on the top row are snapshots displaying the phase $\arg[\psi]$ after rotations $\theta = \theta_0$, $\theta = \theta_0 + \pi$ and $\theta = \theta_0 + 2\pi$ respectively. Since the Lagrange multipliers $\lambda_N^{min}$ and $\lambda_z^{min}$ are not commensurate, $\arg[\psi](t)$ and $\theta(t)$ feature different periodicities. As a result after $\pi$ rotations of the pair the phase profile is not simply an overall phase multiplication.</figcaption></figure>

## 4 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 $\lambda_z^{min} \equiv -\Omega$. In addition the phase of the wave function evolves with a rate that depends on $\lambda_N^{min}$. This coincides with the rate of change in the phase $\arg[\psi](t)$, measured at the center of the trap. Similarly, the rate of change of $\theta(t)$ that measures the rotation of the pair, equals $\lambda_z^{min}$. 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 $\pi$.

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:

$$i\partial_t \psi = \Omega \mathbf{L}_z^{\text{cov}} \psi \equiv \Omega(e_z \cdot \mathbf{L} + \phi e_z \cdot \mathbf{x} \wedge \mathbf{A})\psi. \quad (4.1)$$
Here $\mathbf{L} = -i\mathbf{x} \wedge \nabla$ is the canonical angular momentum, the parameter $\phi = \lambda_N^{min}/\lambda_z^{min}$, and the flux of $\mathbf{A} = \nabla \tan^{-1} x/y$ equals $2\pi$ around any counter-clockwise path that encircles the $z$-axis once. Notably, in the case of two identical vortices $\mathbf{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 $\mathbf{L}_z^{cov}$. We note that when $\phi \neq 0$, the rotations around $z$-axis, i.e. the changes in the polar angle $\theta$, are generated by the covariant angular momentum.

To find the solutions of (4.1) we first remove $\mathbf{A}$ by phase multiplication, and then proceed with separation of variables in terms of polar coordinates $(r, \theta)$.

$$ \psi(\mathbf{x}, t) = e^{-iEt} \mathcal{X}(r, \theta) e^{-\frac{i\phi}{\mathbf{x}_0} \int_{\mathbf{x}_0}^{\mathbf{x}} \mathbf{A} \cdot d\mathbf{l}}, \quad (4.2) $$

The single valuedness of the wave function $\psi(\mathbf{x}, t)$ implies that $\mathcal{X}(r, \theta)$ is generally not single valued as it obeys

$$ \mathbf{L}_z \mathcal{X}_n = (n + \phi) \mathcal{X}_n. \quad (4.3) $$

The solution of the time evolution (4.1) is thus

$$ \psi(\mathbf{x}, t) = e^{-\frac{i\phi}{\mathbf{x}_0} \int_{\mathbf{x}_0}^{\mathbf{x}} \mathbf{A} \cdot d\mathbf{l}} \sum_n f_n(r) e^{-i(n+\phi)(\Omega t - \theta)}. \quad (4.4) $$

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 $\psi_{min}(\mathbf{x})$ of (2.6), (2.1) at $t = 0$.

# 5 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–46]. Additional experimental set-ups can be built on thermal quenching in combination with the Kibble-Zurek mechanism [47–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$Rb atoms [53, 54].

## Acknowledgments

AJN thanks C. Pethick and K. Sacha for discussions. We thank Rima El Kosseifi for discussions. The work by DJ and AJN is supported by the Carl Trygger Foundation Grant CTS 18:276 and by the Swedish Research Council under Contract No. 2018-04411. DJ and AJN also acknowledge collaboration under COST Action CA17139. The research by AJN was also partially supported by Grant No. 0657-2020-0015 of the Ministry of Science and Higher Education of Russia. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at National Supercomputer Center at Linköping, Sweden.

## 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 $\Psi$, whose dynamics obey

$$i\hbar\partial_t\Psi = \left[-\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{x}) + \frac{N_a}{N}g_0|\Psi|^2\right]\Psi, \qquad (A.1a)$$

$$\text{where } N = \int |\Psi|^2, \text{ and } g_0 = \frac{4\pi\hbar^2 a_s}{m}. \qquad (A.1b)$$

Here, $\nabla^2 \equiv \nabla \cdot \nabla$ is the Laplace operator, $m$ is the mass of the atoms, $\hbar$ 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

$$V(\mathbf{x}) = \frac{m}{2}(\omega_x^2 x^2 + \omega_y^2 y^2 + \omega_z^2 z^2), \qquad (A.2)$$

where $\omega_x \le \omega_y \le \omega_z$ are the trapping frequencies in the different spatial directions. The dimensionless units are obtained by the following rescaling

$$t = t_s \tilde{t}, \quad \mathbf{x} = x_s \tilde{x}, \quad \Psi = x_s^{-3/2} \tilde{\psi}, \quad \text{where } t_s = \frac{1}{\omega_x}, \quad \text{and} \quad x_s = \sqrt{\frac{\hbar}{m\omega_x}}. \qquad (A.3)$$
Hence, after dropping the ~, the dimensionless Gross-Pitaevskii equation is

$$i \partial_t \psi = \left[ -\frac{1}{2} \nabla^2 + V(\mathbf{x}) + g |\psi|^2 \right] \psi, \qquad (A.4a)$$

$$\text{where } V(\mathbf{x}) = \frac{1}{2} \sum_{\alpha=x,y,z} \left(\frac{\omega_{\alpha}}{\omega_x}\right)^2 x_{\alpha}^2, \quad N = \int |\psi|^2, \text{ and } g = \frac{4\pi a_s N_a}{x_s N}. \qquad (A.4b)$$

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 $\omega_x \approx \omega_y$ and $\omega_z \gg \omega_x$. This results in a further rescaling of the coupling constant $g_{2d} = g_{3d}\sqrt{\frac{\omega_z}{2\pi\omega_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$.

<table><caption>**Table 1.** Typical physical values for various atomic Bose-Einstein condensates, and the resulting scale and dimensionless coupling. The reduced Planck constant is $\hbar = 1.05 \times 10^{-34}$ J.s.</caption><thead><tr><th>Atoms</th><th>m (kg)</th><th>a<sub>s</sub> (m)</th><th>&omega;<sub>x</sub> (rad.s<sup>-1</sup>)</th><th>N<sub>a</sub></th></tr></thead><tbody><tr><td><sup>87</sup>Rb [1]</td><td>1.44 &times; 10<sup>-25</sup></td><td>5.1 &times; 10<sup>-9</sup></td><td>20</td><td>10<sup>2</sup> &sim; 10<sup>7</sup></td></tr><tr><td><sup>23</sup>Na [2]</td><td>3.817 &times; 10<sup>-26</sup></td><td>2.75 &times; 10<sup>-9</sup></td><td>2&pi; &times; 200</td><td>&times;10<sup>4</sup> &sim; 5 &times; 10<sup>5</sup></td></tr></tbody></table><table><thead><tr><th>Atoms</th><th>x<sub>s</sub> (m)</th><th><sup>4&pi;a<sub>s</sub></sup>&frasl;<sub>x<sub>s</sub></sub></th><th>g (&times;N = 1)</th></tr></thead><tbody><tr><td><sup>87</sup>Rb [1]</td><td>6.0 &times; 10<sup>-6</sup></td><td>1.1 &times; 10<sup>-2</sup></td><td>1 &sim; 10<sup>5</sup></td></tr><tr><td><sup>23</sup>Na [2]</td><td>1.5 &times; 10<sup>-6</sup></td><td>2.3 &times; 10<sup>-2</sup></td><td>10<sup>2</sup> &sim; 10<sup>4</sup></td></tr></tbody></table>

## 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 co-rotating frame [13–17, 20–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 $\Omega$. The free energy $F^{\Omega}$ 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 $\Omega$. 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:

$$F^{\Omega} = F - \Omega L_z, \quad \text{where } L_z \equiv \int d^2x \, \psi^* e_z \cdot (-i\hbar \mathbf{x} \wedge \nabla) \psi. \qquad (B.1)$$

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^{\Omega}$ 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

$$F_{\lambda} \equiv F + \lambda_{z}(L_{z} - l_{z}) + \lambda_{N}(N - 1). \qquad (B.2)$$

Here $\lambda_z$, $\lambda_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^{\Omega}$ is a critical point of the following free energy

$$F_{\lambda}^{\Omega} \equiv F^{\Omega} + \lambda_{N}(N - 1) = F + \Omega L_{z} + \lambda_{N}(N - 1), \qquad (B.3)$$

where $\lambda_N$ is the (only) Lagrange multiplier that enforces $N = 1$. Note that the present Lagrange multiplier $\lambda_N$ is commonly treated as a chemical potential that is an externally controlled parameter, in par with $\Omega$. This is not the case here, as $\lambda_N$ has to be adjusted to satisfy the constraint. The critical points of $F_{\lambda}$ and $F_{\lambda}^{\Omega}$ obey respectively

$$F_{\lambda}: -\frac{1}{2}\nabla^2\psi + \frac{|\mathbf{x}|^2}{2}\psi + g|\psi|^2\psi = -\lambda_N\psi + i\lambda_z e_z \cdot \mathbf{x} \wedge \nabla\psi, \qquad (B.4a) \\ \text{and the conditions } L_z = l_z \text{ and } N = 1,$$

$$F_{\lambda}^{\Omega}: -\frac{1}{2}\nabla^2\psi + \frac{|\mathbf{x}|^2}{2}\psi + g|\psi|^2\psi - i\Omega e_z \cdot \mathbf{x} \wedge \nabla\psi = -\lambda_N\psi, \qquad (B.4b) \\ \text{and the condition } N = 1.$$

Note that in a numerical simulation $\Omega$ is an input parameter that is fixed during the simulation, while $\lambda_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.

Comparison of the minimal energy solutions of the two equations (B.4), can be seen in figure 4. The results are indeed quite different. In particular, the free energy (evaluated in
<table><caption>**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.</caption><thead><tr><th>Critical pts. of</th><th>Obey the equation<br>Obey the equation</th><th>Satisfy</th><th>Control params.</th><th>Mult.</th><th>Output</th></tr></thead><tbody><tr><td>F<sub>λ</sub></td><td>-½∇<sup>2</sup>ψ + |x|<sup>2</sup>ψ + g|ψ|<sup>2</sup>ψ =<br>-λ<sub>N</sub>ψ + iλ<sub>z</sub> e<sub>z</sub> · x ∧ ∇ψ</td><td>L<sub>z</sub> = l<sub>z</sub><br>and N = 1</td><td>g, l<sub>z</sub></td><td>λ<sub>N</sub>,<br>λ<sub>z</sub></td><td>λ<sub>N</sub><sup>min</sup>, λ<sub>z</sub><sup>min</sup><br>ψ<sub>min</sub></td></tr><tr><td>F<sub>λ</sub><sup>Ω</sup></td><td>-½∇<sup>2</sup>ψ + |x|<sup>2</sup>ψ + g|ψ|<sup>2</sup>ψ<br>-iΩ e<sub>z</sub> · x ∧ ∇ψ = -λ<sub>N</sub>ψ</td><td>N = 1</td><td>g, Ω</td><td>λ<sub>N</sub></td><td>λ<sub>N</sub><sup>min</sup>,<br>ψ<sub>min</sub></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Detailed comparison of the solutions to the two equations (B.4) for $g = 100$. The panels on the top row display the energy, in a non-rotating frame, as a function of angular velocity $\Omega$ (left) and value of angular momentum $L_z$ (right). The empty (red) circles denote the results from the constrained minimization at given value of the angular momentum $L_z$. The filled (green) circles show the results from the minimization of a rotated condensate at given angular velocity $\Omega$.</figcaption></figure>

a non-rotating laboratory frame) that is displayed on the panels on the top row of figure 4, show that the minimization with fixed angular frequency $\Omega$ cannot support solutions with arbitrary values of angular momentum. Indeed, there are excluded ranges

$$L_z \in [0, 1[ \cup ]1, 1.457[ \cup ]1.75, 2.68[ \dots, \qquad (B.5)$$

where there are no rotation-induced stable vortex states (see the green dots on the top right panel). This is in agreement with the known results for rotating Bose-Einstein condensates;
compare for example the green dots in the top left panel of figure 4, e.g. with the results shown in figure 2 of [20]. In contrast, in the present case, there cannot be any excluded values of the angular momentum in the equation (B.4a) (see the red dots on the top right panel). Those vortex configurations with angular momentum in the excluded ranges (B.5) describe eccentric asymmetric vortex configurations that precess around the trap center.

Observe that in the case of equation (B.4a), there is multi-valuedness in the free energy (red dots), as a function of the Lagrange multiplier $\lambda_z$; see for example the top left panel at $\lambda_z^{min} \approx 0.55$ which displays 3 different states at the same value of $\lambda_z$. These are a single eccentric vortex, an asymmetric vortex pair, and a symmetric vortex pair for increasing free energy, respectively. These configurations have both a different angular momentum value, and a different value of the Lagrange multiplier $\lambda_N^{min}$ (see the bottom panels) which plays the role of a chemical potential. The present comparisons, summarized in figure 4, show that the two equations (B.4a) and (B.4b) do indeed describe two quite different physical scenarios.

## 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 $\Omega$ which a bounded open subset of $\mathbb{R}^2$ and denote $\partial\Omega$ its boundary. $H(\Omega)$ stands for the Hilbert space, such that a function belonging to $H(\Omega)$, and its weak derivatives have a finite $L^2$-norm. Furthermore, $\mathcal{H}(\Omega) = \{u+iv | u, v \in H(\Omega)\}$ denotes the Hilbert spaces of complex-valued functions. The Hilbert spaces of real- and complex-valued functions are equipped with the inner products $\langle \cdot, \cdot \rangle$, defined as:

$$ \langle u, v \rangle = \int_{\Omega} uv, \quad \text{for } u, v \in H(\Omega), \quad \text{and} \quad \langle u, v \rangle = \int_{\Omega} u^* v, \quad \text{for } u, v \in \mathcal{H}(\Omega). \quad (C.1) $$

The spatial domain $\Omega$ is discretized as a mesh of triangles using for the Delaunay-Voronoi algorithm, and the regular partition $\mathcal{T}_h$ of $\Omega$ 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_h^{(2)}$ as the 2-nd order Lagrange finite-element subspace of $H(\Omega)$, and correspondingly $\mathcal{P}_h^{(2)}$ for $\mathcal{H}(\Omega)$. 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 $\psi \mapsto \psi^{(h)} \in \mathcal{P}_h^{(2)}$. This describes a linear vector
space of finite dimension, for which a basis can be found. The canonical basis consists of the shape functions $\phi_k(\mathbf{x})$, and thus

$$V_h(\mathcal{T}_h, \mathbf{P}^{(2)}) = \left\{ w(\mathbf{x}) = \sum_{k=1}^{M} w_k \phi_k(\mathbf{x}), \phi_k(\mathbf{x}) \in \mathbf{P}_h^{(2)} \right\}. \quad (C.2)$$

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(\mathbf{x}) = \sum_{k=1}^{M} w_k \phi_k(\mathbf{x})$, on a given basis of shape functions $\phi_k(\mathbf{x})$ of the polynomial functions $\mathbf{P}^{(2)}$ for the triangle $T_{i_k}$. The finite element space $V_h(\mathcal{T}_h, \mathbf{P}^{(2)})$ hence denotes the space of continuous, piecewise quadratic functions of $x, y$ on each triangle of $\mathcal{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^{\text{aug}}$ is defined as

$$F^{\text{aug}} := F + \frac{\mu}{2} \left( \sum_j C_j^2 \right) + \sum_j \lambda_j C_j. \quad (C.3)$$

Here $\mu$ is a penalty parameter and $\lambda_j$ are the Lagrange multipliers associated with the conditions $C_j$. In the Augmented Lagrangian Method, the augmented Lagrangian is minimized, and the penalty $\mu$ is iteratively increased while the multipliers $\lambda_j \leftarrow \lambda_j + \mu 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–61].

In this work, the minimized functional $F$ is the dimensionless Gross-Pitaevskii free energy [17]

$$F = \int d^2x \left\{ \frac{1}{2} |\nabla \psi|^2 + V(\mathbf{x}) |\psi|^2 + \frac{g}{2} |\psi|^4 \right\}, \quad (C.4)$$

where $g$ is the dimensionless coupling that accounts for the interatomic interactions, and $V(\mathbf{x})$ is the trapping potential. In the main body we considered an harmonic oscillator trapping potential $V(\mathbf{x}) = |\mathbf{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

$$C_N = N-1 \equiv \int d^2x |\psi|^2 - 1 \quad \text{and} \quad C_z = L_z - l_z \equiv \int d^2x \psi^* [e_z \cdot (-i\mathbf{x} \wedge \nabla)] \psi - l_z. \quad (C.5)$$
Hence, the variations of the augmented Lagrangian with respect to $\psi^*$ and $\lambda_i$ give the set of equations

$$-\frac{1}{2}\nabla^2\psi + [V(\mathbf{x}) + g|\psi|^2]\psi + \tilde{\lambda}_N\psi + \tilde{\lambda}_z(e_z \cdot \mathbf{L})\psi = 0, \quad \text{where } \tilde{\lambda}_i = (\mu C_i + \lambda_i) \quad (\text{C.6a})$$

$$C_i = 0, \quad \text{and} \quad i = N, z, \quad (\text{C.6b})$$

where the angular momentum operator is $\mathbf{L} = -i\mathbf{x} \wedge \nabla$. The associated weak form, is obtained by multiplying the equation (C.6a) by test functions $\psi_w \in H(\Omega)$, integrating over $\Omega$ 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

$$[\psi_w] \cdot F'(\psi) = \frac{1}{2} \langle \nabla \psi_w, \nabla \psi \rangle + \langle \psi_w, [V(\mathbf{x}) + g|\psi|^2] \psi \rangle + \tilde{\lambda}_N \langle \psi_w, \psi \rangle + \tilde{\lambda}_z \langle \psi_w, e_z \cdot \mathbf{L} \psi \rangle. \quad (\text{C.7})$$

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 derived by multiplying (C.6a) by $\psi^*$ and integrating over the domain $\Omega$ (and again integrating by part the Laplace operator). The error is alternatively obtained by replacing the test functions $\psi_w$ by $\psi$ in (C.7):

$$\int \frac{1}{2} |\nabla \psi|^2 + [V(\mathbf{x}) + g|\psi|^2] |\psi|^2 + \tilde{\lambda}_q \int |\psi|^2 + \tilde{\lambda}_z \int \psi^*(e_z \cdot \mathbf{L}) \psi = 0 \quad (\text{C.8a})$$

$$F + \frac{g}{2} \int |\psi|^4 + \tilde{\lambda}_N(C_N + 1) + \tilde{\lambda}_z(C_z + l_z) = 0 \quad (\text{C.8b})$$

$$1 + \frac{\sum_i \tilde{\lambda}_i (C_i - l_i)}{F + \frac{g}{2} \int |\psi|^4} = \text{err.} \quad (\text{C.8c})$$

In practice 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:

$$i \partial_t \psi = -\frac{1}{2} \nabla^2 \psi + [V(\mathbf{x}) + g |\psi|^2] \psi. \quad (\text{C.9})$$
<figure><img src="image_6.png" /><figcaption>**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).</figcaption></figure>

The weak from, obtained by multiplying by test functions $\psi_w \in \mathcal{H}(\Omega)$ and integrating by parts reads as, in terms of the inner products (C.1),

$$ \langle \psi_w, i \partial_t \psi \rangle = \frac{1}{2} \langle \nabla \psi_w, \nabla \psi \rangle + \langle \psi_w, [V(\mathbf{x}) + g|\psi|^2] \psi \rangle . \quad (C.10) $$

The discretization of the time turns the continuous evolution into a recursive series over the uniform partition $\{t\}_{n=0}^N$ of the time variable. The time variable is discretized as $t = n\Delta t$ and the wave function at the step $n$ is $\psi_n := \psi(n\Delta 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

$$ \partial_t \psi := \delta_t \psi_n = \frac{\psi_{n+1} - \psi_n}{\Delta t}, \quad \text{and} \quad \bar{\psi}_n = \frac{\psi_{n+1} + \psi_n}{2}, \quad (C.11) $$

the Crank-Nicolson scheme for the time-dependent Gross-Pitaevskii equation (C.10) reads as:

$$ \langle \psi_w, i \delta_t \psi_n \rangle = \frac{1}{2} \langle \nabla \psi_w, \nabla \bar{\psi}_n \rangle + \langle \psi_w, V(\mathbf{x}) \bar{\psi}_n \rangle + \langle \psi_w, g | \bar{\psi}_n |^2 \bar{\psi}_n \rangle . \quad (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 $\bar{\psi}_n \approx (3\psi_n - \psi_{n-1})/2$. Next, defining the time-discretized
operators:

$$ \mathcal{O}_1 \psi = \left\langle \psi_w, \frac{i\psi}{\Delta t} \right\rangle - \left\langle \nabla \psi_w, \frac{1}{4} \nabla \psi \right\rangle - \left\langle \psi_w, \frac{1}{2} V(\mathbf{x}) \psi \right\rangle, \quad (C.13a) $$

$$ \mathcal{O}_2 \psi = \left\langle \psi_w, \frac{i\psi}{\Delta t} \right\rangle + \left\langle \nabla \psi_w, \frac{1}{4} \nabla \psi \right\rangle + \left\langle \psi_w, \frac{1}{2} V(\mathbf{x}) \psi \right\rangle, \quad (C.13b) $$

$$ \mathcal{U}_n \psi = \left\langle \psi_w, \frac{g}{8} |3\psi_n - \psi_{n-1}|^2 \psi \right\rangle, \quad (C.13c) $$

allows to rewrite Eq. (C.12) in the compact form

$$ \mathcal{O}_1 \psi_{n+1} = \mathcal{O}_2 \psi_n + \mathcal{U}_n (3\psi_n - \psi_{n-1}). \quad (C.14) $$

Hence the time-evolution is formally given by the recursion

$$ \psi_{n+1} = \mathcal{O}_1^{-1} [\mathcal{O}_2 \psi_n + \mathcal{U}_n (3\psi_n - \psi_{n-1})]. \quad (C.15) $$

The spatial discretization is achieved by replacing the wave function $\psi$ by its finite-element space representation $\psi^{(h)} \in V_h(\mathcal{T}_h, \mathcal{P}^{(2)})$ in the time-dependent Gross-Pitaevskii equation (C.15). The test functions $\psi_w$ now take values in the same discrete space as $\psi^{(h)}$. Denoting the matrix representation of the time-discretized evolution operators (C.13), as:

$$ \mathcal{O}_1 \mapsto \mathbf{M}_\psi, \quad \mathcal{O}_2 \mapsto \mathbf{N}_\psi, \quad \text{and} \quad \mathcal{U}_n (3\psi_n - \psi_{n-1}) \mapsto \mathbf{R}_\psi, \quad (C.16) $$

the recursion Eq. (C.15) reduce to a linear system that read as:

$$ [\mathbf{M}_\psi] [\psi_{n+1}^{(h)}] - [\mathbf{N}_\psi] [\psi_n^{(h)}] = [\mathbf{R}_\psi]. \quad (C.17) $$

The vector $\mathbf{R}_\psi$ which is a function of $\psi_n^{(h)}$ and $\psi_{n-1}^{(h)}$, has to be recalculated at each step. The matrices $\mathbf{M}_\psi$ and $\mathbf{N}_\psi$, on the other hand, are constant matrices to be allocated just once and are in principle easily preconditioned. Finally the recursion is thus given by

$$ [\psi_{n+1}^{(h)}] = [\mathbf{M}_\psi]^{-1} ([\mathbf{N}_\psi] [\psi_n^{(h)}] + [\mathbf{R}_\psi]), \quad (C.18) $$

which requires recalculating the vectors $\mathbf{R}_\psi$ 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.

## References

[1] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman and E. A. Cornell, *Observation of Bose-Einstein Condensation in a Dilute Atomic Vapor*, Science 269 (1995) 198.
[2] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn et al., *Bose-Einstein Condensation in a Gas of Sodium Atoms*, Physical Review Letters **75** (1995) 3969.

[3] C. C. Bradley, C. A. Sackett and R. G. Hulet, *Bose-Einstein Condensation of Lithium: Observation of Limited Condensate Number*, Physical Review Letters **78** (1997) 985.

[4] D. C. Aveline, J. R. Williams, E. R. Elliott, C. Dutenhoffer, J. R. Kellogg, J. M. Kohel et al., *Observation of Bose-Einstein condensates in an Earth-orbiting research lab*, Nature **582** (2020) 193.

[5] E. R. Elliott, M. C. Krutzik, J. R. Williams, R. J. Thompson and D. C. Aveline, *NASA's Cold Atom Lab (CAL): system development and ground test status*, npj Microgravity **4** (2018) 16.

[6] I. Bloch, J. Dalibard and S. Nascimbène, *Quantum simulations with ultracold quantum gases*, Nature Physics **8** (2012) 267.

[7] M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S. Hall, C. E. Wieman and E. A. Cornell, *Vortices in a Bose-Einstein Condensate*, Physical Review Letters **83** (1999) 2498.

[8] F. Chevy, K. W. Madison and J. Dalibard, *Measurement of the Angular Momentum of a Rotating Bose-Einstein Condensate*, Physical Review Letters **85** (2000) 2223.

[9] C. Raman, J. R. Abo-Shaeer, J. M. Vogels, K. Xu and W. Ketterle, *Vortex Nucleation in a Stirred Bose-Einstein Condensate*, Physical Review Letters **87** (2001) 210402.

[10] J. R. Abo-Shaeer, C. Raman, J. M. Vogels and W. Ketterle, *Observation of Vortex Lattices in Bose-Einstein Condensates*, Science **292** (2001) 476.

[11] E. P. Gross, *Structure of a quantized vortex in boson systems*, Il Nuovo Cimento **20** (1961) 454.

[12] L. P. Pitaevskii, *Vortex Lines in an Imperfect Bose Gas*, Soviet Journal of Experimental and Theoretical Physics **13** (1961) 451.

[13] L. P. Pitaevskii and S. Stringari, *Bose-Einstein condensation*, no. 116 in International series of monographs on physics. Clarendon Press, Oxford New York, 2003.

[14] E. H. Lieb and R. Seiringer, *Derivation of the Gross-Pitaevskii Equation for Rotating Bose Gases*, Communications in Mathematical Physics **264** (2006) 505.

[15] C. J. Pethick and H. Smith, *Bose-Einstein Condensation in Dilute Gases*. Cambridge University Press, ii ed., 2008, 10.1017/CBO9780511802850.

[16] A. L. Fetter, *Rotating trapped Bose-Einstein condensates*, Reviews of Modern Physics **81** (2009) 647.

[17] W. Bao and Y. Cai, *Mathematical theory and numerical methods for Bose-Einstein condensation*, Kinetic & Related Models **6** (2013) 1.

[18] P. C. Haljan, I. Coddington, P. Engels and E. A. Cornell, *Driving Bose-Einstein-Condensate Vorticity with a Rotating Normal Cloud*, Physical Review Letters **87** (2001) 210403.

[19] A. L. Fetter and A. A. Svidzinsky, *Vortices in a trapped dilute Bose-Einstein condensate*, Journal of Physics: Condensed Matter **13** (2001) R135.

[20] D. A. Butts and D. S. Rokhsar, *Predicted signatures of rotating Bose-Einstein condensates*, Nature **397** (1999) 327.
[21] A. Aftalion and Q. Du, *Vortices in a rotating Bose-Einstein condensate: Critical angular velocities and energy diagrams in the Thomas-Fermi regime*, *Physical Review A* **64** (2001) 063603.

[22] R. Seiringer, *Gross-Pitaevskii Theory of the Rotating Bose Gas*, *Communications in Mathematical Physics* **229** (2002) 491.

[23] B. P. Anderson, P. C. Haljan, C. E. Wieman and E. A. Cornell, *Vortex Precession in Bose-Einstein Condensates: Observations with Filled and Empty Cores*, *Physical Review Letters* **85** (2000) 2857.

[24] R. Fletcher, *Practical Methods of Optimization*. Wiley, Chichester New York, 1987, 10.1002/9781118723203.

[25] J. Nocedal and S. J. Wright, *Numerical Optimization*, Springer Series in Operations Research. Springer, 1999, 10.1007/978-0-387-40065-5.

[26] B. Mottelson, *Yrast Spectra of Weakly Interacting Bose-Einstein Condensates*, *Physical Review Letters* **83** (1999) 2695.

[27] G. M. Kavoulakis, B. Mottelson and C. J. Pethick, *Weakly interacting Bose-Einstein condensates under rotation*, *Physical Review A* **62** (2000) 063605.

[28] A. Shapere and F. Wilczek, *Classical Time Crystals*, *Physical Review Letters* **109** (2012) 160402.

[29] F. Wilczek, *Quantum Time Crystals*, *Physical Review Letters* **109** (2012) 160401.

[30] A. Alekseev, J. Dai and A. J. Niemi, *Provenance of classical Hamiltonian time crystals*, *Journal of High Energy Physics* **2020** (2020) 35.

[31] J. Dai, A. J. Niemi and X. Peng, *Classical Hamiltonian time crystals – general theory and simple examples*, *New Journal of Physics* **22** (2020) 085006.

[32] J. M. Leinaas and J. Myrheim, *On the Theory of Identical Particles*, *Il Nuovo Cimento B* **37** (1977) 1.

[33] F. Wilczek, *Quantum Mechanics of Fractional-Spin Particles*, *Physical Review Letters* **49** (1982) 957.

[34] R. Iengo and K. Lechne, *Anyon quantum mechanics and Chern-Simons theory*, *Physics Reports* **213** (1992) 179.

[35] A. Lerda, *Anyons: quantum mechanics of particles with fractional statistics*, vol. 14 of *Lecture Notes in Physics*. Lecture Notes in Physics, Berlin New York, 1992, 10.1007/978-3-540-47466-1.

[36] A. Y. Kitaev, *Fault-tolerant quantum computation by anyons*, *Annals of Physics* **303** (2003) 2.

[37] J. Nakamura, S. Liang, G. C. Gardner and M. J. Manfra, *Direct observation of anyonic braiding statistics*, *Nature Physics* **16** (2020) 931.

[38] E. H. Lieb and R. Seiringer, *Proof of Bose-Einstein Condensation for Dilute Trapped Gases*, *Physical Review Letters* **88** (2002) 170409.

[39] J. E. Marsden and T. S. Ratiu, *Introduction to Mechanics and Symmetry*. Springer New York, New York, 1999, 10.1007/978-0-387-21792-5.

[40] F. Hecht, *New development in freefem++*, *Journal of Numerical Mathematics* **20** (2012) 251.
[41] see Supplementary Materials for animations.

[42] F. Wilczek, *Magnetic Flux, Angular Momentum, and Statistics*, *Physical Review Letters* **48** (1982) 1144.

[43] K.-P. Marzlin, W. Zhang and E. M. Wright, *Vortex Coupler for Atomic Bose-Einstein Condensates*, *Physical Review Letters* **79** (1997) 4728.

[44] E. L. Bolda and D. F. Walls, *Creation of vortices in a Bose-Einstein condensate by a Raman technique*, *Physics Letters A* **246** (1998) 32.

[45] R. Dum, J. I. Cirac, M. Lewenstein and P. Zoller, *Creation of Dark Solitons and Vortices in Bose-Einstein Condensates*, *Physical Review Letters* **80** (1998) 2972.

[46] A. E. Leanhardt, A. Görlitz, A. P. Chikkatur, D. Kielpinski, Y. Shin, D. E. Pritchard et al., *Imprinting Vortices in a Bose-Einstein Condensate using Topological Phases*, *Physical Review Letters* **89** (2002).

[47] T. W. B. Kibble, *Topology of Cosmic Domains and Strings*, *Journal of Physics A: Mathematical and General* **9** (1976) 1387.

[48] W. H. Zurek, *Cosmological experiments in superfluid helium?*, *Nature* **317** (1985) 505.

[49] W. H. Zurek, *Cosmological experiments in condensed matter systems*, *Physics Reports* **276** (1996) 177.

[50] J. R. Anglin and W. H. Zurek, *Vortices in the Wake of Rapid Bose-Einstein Condensation*, *Physical Review Letters* **83** (1999) 1707.

[51] V. M. H. Ruutu, V. B. Eltsov, A. J. Gill, T. W. B. Kibble, M. Krusius, Y. G. Makhlin et al., *Vortex formation in neutron-irradiated superfluid $^3$He as an analogue of cosmological defect formation*, *Nature* **382** (1996) 334.

[52] V. B. Eltsov, T. W. B. Kibble, M. Krusius, V. M. H. Ruutu and G. E. Volovik, *Composite Defect Extends Analogy between Cosmology and $^3$He*, *Physical Review Letters* **85** (2000) 4739.

[53] C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis and B. P. Anderson, *Spontaneous vortices in the formation of Bose-Einstein condensates*, *Nature* **455** (2008) 948.

[54] D. V. Freilich, D. M. Bianchi, A. M. Kaufman, T. K. Langin and D. S. Hall, *Real-Time Dynamics of Single Vortex Lines and Vortex Dipoles in a Bose-Einstein Condensate*, *Science* **329** (2010) 1182.

[55] D. V. Hutton, *Fundamentals of Finite Element Analysis*, Engineering Series. McGraw-Hill, 2003, 10.1036/0072922362.

[56] J. Reddy, *An Introduction to the Finite Element Method*. McGraw-Hill Education, 2005.

[57] P. E. Gill, W. Murray and M. H. Wright, *Practical optimization*. Academic Press, 1981.

[58] E. G. Birgin and J. M. Martínez, *Improving ultimate convergence of an augmented Lagrangian method*, *Optimization Methods and Software* **23** (2008) 177.

[59] M. R. Hestenes and E. Stiefel, *Methods of conjugate gradients for solving linear systems*, *Journal of Research of the National Bureau of Standards* **49** (1952) 409.

[60] R. Fletcher and C. M. Reeves, *Function minimization by conjugate gradients*, *The Computer Journal* **7** (1964) 149.

[61] E. Polak and G. Ribière, *Note sur la convergence de directions conjuguées*, *Revue française d'informatique et de recherche opérationnelle. Série rouge* **3** (1969) 35.
[62] J. Crank and P. Nicolson, *A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type*, *Advances in Computational Mathematics* **6** (1996) 207.