See the published version in J. High Energ. Phys.

Bose-Einstein condensates, in two spatial dimensions are described by the (dimensionless) Gross-Pitaevskii free energy $F=\int {\cal F} d^2{\bf r}$, whose density reads as: \begin{equation*} {\cal F}= \frac{1}{2}|\Grad\psi|^2+V(\x)|\psi|^2+\frac{g}{2}|\psi|^4 \,, ~~\text{where}~~V(\x)=\frac{|\x|^2}{2} \,. \end{equation*} Here, $\psi(\x)$ is the macroscopic wave function, the potential $V(\x)$ is an harmonic oscillator trapping potential, and the dimensionles coupling constant $g$ parametrizes the strength of the interatomic interactions. The theory features two conserved quantities. First is the particle number \begin{equation*} N=\int |\psi|^2d^2{\bf r}\,, \end{equation*} which in dimensionless units is normalized: $N=1$. The second conserved quantity is the axial compenent of the angular momentum: \begin{equation*} L_z=\int \psi^*\left[{\bf e}_z\cdot(-i\x\times\Grad)\right]\psi\, d^2{\bf r}\,. \end{equation*} $L_z$ is the macroscopic angular momentum. The investigated values of the angular momentum below are $L_z=l_z\in[0,2)$.

The numerical algorithms that is used solve the constrained optimization problem is the Augmented Lagrangian Method (ALM), which combines penalty methods with the method of Lagrange multipliers. In terms of the original free energy $F$, and the set of conditions $C_j$, the augmented Lagrangian $F^\mathrm{aug}$ defines the unconstrained objective function as: \begin{equation*} F^\mathrm{aug}:= F + \frac{\mu}{2} \left( \sum_j C_j^2 \right) +\sum_j \lambda_jC_j \,, \end{equation*} where $\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 Lagragian is minimized, and the penalty $\mu$ is iteratively increased while $\lambda_j\leftarrow\lambda_j+\mu C_j$ until all constraints are satisfied with a specified accuracy. The ALM algorithm converges in a finite number of iterations, and the minimization algorithm within each ALM iteration is a non-linear conjugate gradient algorithm.

The displayed quantities are : the minimal energy configurations that simultaneously satisfy $N=\int|\psi|^2=1$ and $L_z=l_z\in[0,2)$. Here, the interaction coupling is $g=400$. The elevation of the semitransparent surface stems for the condensate density $|\psi|^2$, while the coloring indicates the value of the phase $\arg[\psi]$ of the macroscopic wave-function. The coloring displays both the surface that is associated with the density, and is also projected onto the $xy$-plane.

The time evolution of the Bose-Einstein condensate is determined by the Time-Dependent Gross--Pitaevskii equation \begin{equation*} i\partial_t\psi=-\frac{1}{2}\nabla^2\psi +\left[V(\x)+g|\psi|^2\right]\psi \,. \end{equation*} The strategy is to use a Crank-Nicolson algorithm 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 extrapolation. \begin{equation*} \partial_t\psi :=\delta_t\psi_n = \frac{\psi_{n+1} - \psi_{n}}{\Delta t} \,,~~\text{and}~~ \bar{\psi}_n = \frac{\psi_{n+1} + \psi_{n}}{2}\,. \end{equation*} For faster computations, the fields in the nonlinear terms are expressed by using an extrapolation of the previous time steps that retains the same order of truncation error as the rest of time series. The averaged wave function in the non-linear term thus becomes $\bar{\psi}_n\approx(3\psi_{n}-\psi_{n-1})/2$.

The displayed quantities are : The left panels show the density $|\psi|^2$, while the right panels display the phase $\arg[\psi]$ of the condensate. The starting configuration is the minimum energy configuration $\psi_{min}(\x)$ for $N=1$, and an angular momentum $L_z=l_z$.

1. J. Garaud, J. Dai and A. Niemi,
   Timecrystalline vorticity and anyonic exchange in a cold atom Bose-Einstein condensate.
   J. High Energ. Phys. 2021, 157 (2021) Link , arXiv .