Phantom vortices: hidden angular momentum in ultracold dilute Bose-Einstein condensates

Vortices are essential to angular momentum in quantum systems such as ultracold atomic gases. The existence of quantized vorticity in bosonic systems stimulated the development of the Gross-Pitaevskii mean-field approximation. However, the true dynamics of angular momentum in finite, interacting many-body systems like trapped Bose-Einstein condensates is enriched by the emergence of quantum correlations whose description demands more elaborate methods. Herein we theoretically investigate the full many-body dynamics of the acquisition of angular momentum by a gas of ultracold bosons in two dimensions using a standard rotation procedure. We demonstrate the existence of a novel mode of quantized vorticity, which we term the phantom vortex. Contrary to the conventional mean-field vortex, can be detected as a topological defect of spatial coherence, but not of the density. We describe previously unknown many-body mechanisms of vortex nucleation and show that angular momentum is hidden in phantom vortices modes which so far seem to have evaded experimental detection. This phenomenon is likely important in the formation of the Abrikosov lattice and the onset of turbulence in superfluids.

In this supplementary document, we define the MCTDHB method we use to solve the time-dependent many-body Schrödinger equation and specify the numerical details of our computations. Furthermore, we give a detailed description of the quantities of interest analyzed in the main text, and an analytical discussion of the existence of fragmented vortices.
We then compare our results to a GP mean-field description and assess the generality of our findings for 10 to 10 4 particles. Finally, we provide captions for the enclosed videos and compare our results for short-ranged interactions to computations with a zero ranged Dirac-delta pseudopotential.

MCTDHB Method and Numerical Details
We solve the time-dependent Schrödinger equation with the Hamiltonian described in the main text using the MCTDHB method [1]. MCTDHB uses the ansatz, where |n; t⟩ = |n 1 , ..., n M ; t⟩ is a permanent of M orthonormal orbitals which explicitly depend on time, i.e. a fully symmetrized product of at most M distinct single particle states.
Using a time-dependent variational basis ensures the highest possible accuracy for a given M at all times [2]. In fact, each coefficient and each orbital are treated as independent variational parameters. The basis of the considered Hilbert space is hence a set of (M +N −1)! N !(M −1)! permanents. In the limit M → ∞, the configurations |n; t⟩ span the complete N −particle Hilbert space and the expansion in Eq. 1 becomes formally exact. In practice, a small number of time-dependent single-particle states is enough to obtain numerically exact results [3,4]. In the case of M = 1, Eq. 1 is the GP ansatz and the coupled system of nonlinear integro-differential equations of motion of MCTDHB also simplify to the time-dependent GP equation [1].
We use the recursive implementation of the MCTDHB theory [1] in the MCTDH-X package [5] for the simulations presented in this work. Computations using MCTDH-X employ a discrete variable representation (DVR) to represent the time-dependent orbitals and the one-body Hamiltonian operator as an expansion of a time-independent basis set, see for instance [6]. For the present work, plane waves were used as the primitive basis.

Quantities of Interest
Here we provide the quantities which were analyzed throughout the main text, namely the density, one-body reduced density, natural occupations, natural orbitals, and 1 st -order correlation function. The density, ρ(r; t), is given as the diagonal of the one-body reduced density matrix: The one-body reduced density matrix, ρ (1) , for an N -body system is defined as the partial trace of the N -body density [9]: where the Ψ is assumed to be normalized to 1. The one-body density matrix can be expanded in a complete basis {ξ j (r; t)}, where we have diagonalized the matrix elements ρ kj in the last line to obtain the eigen-  (1) , is defined in terms of the one-body density matrix [10]: The off-diagonal, r ̸ = r ′ , terms in g (1) and the orbitals ϕ i are complex-valued. Thus it is instructive to plot their magnitudes and phases separately. We hence define the phase of and the orbital phases to be For the sake of completeness, we likewise define the phase of the many-body wavefunction, Ψ, as We choose to plot all phases on the interval [−π, π], but any interval of the same length is equally valid. Furthermore, the global phase transformation, leaves the many-body wavefunction (Eq. 1) invariant. Phase shifting an orbital ϕ i by θ i is equivalent to plotting it in the interval [−π + θ i , π + θ i ] with all structures unchanged.
) depends only on phase differences within the same orbital, and is thus invariant to global phase shifts. Therefore, the topological defects we see in the orbitals (Figs. 2 and 3 in the main text, Videos S1 and S2) and g (1) (Fig. 4 and Video S3) are unaffected by global phase transformations.

Analytical Discussion
In the main text we consider the possibility that a fragmented vortex exists in the total density. The simplest way this can happen is through the presence of one phantom vortex in each orbital with all cores coinciding, say at the origin, at some fixed time. This can be formally expressed in polar coordinates as where A j (r) ∈ R is the radial amplitude and q j ∈ Z is the vortex charge of the j th orbital.
Furthermore, the constraint of a single vortex in each orbital requires that A j (0) = 0, Orthonormality of the orbitals then requires The radial integral is strictly nonzero since the radial amplitudes cannot change sign. Thus, the angular integral must vanish. Using the orthonormality of Fourier modes, we conclude that q k ̸ = q j , i.e. the different phantom vortices must have different charges. This case is partially exemplified in ϕ 1 and ϕ 4 (Fig. 3b,e in the main text), which contain a charge 1 and 3 phantom vortex at the origin, respectively. No density vortex is observed because there are no phantom vortices at the origin in ϕ 2 or ϕ 3 . There are many ghost vortices at the edge of the orbitals which do not affect the present argument because they occur in areas where A j is very small. If all the particles were to occupy only the fragments, ϕ 1 and ϕ 4 , there would be a density vortex at the origin.
The other possibility to observe a fragmented vortex in the density is to have coincident phantom vortices of the same charge as well as noncoincident phantom vortices in regions where A j is large enough to ensure orthonormality. For example, ϕ 2 and ϕ 3 (Fig. 2c,d in the main text) each have a charge 1 phantom vortex at the origin and several noncoincident charge 1 phantom vortices near the origin. If all the particles were to occupy only these two orbitals, there would be a density vortex at the origin.
In either case, the existence of a fragmented vortex visible in the total density requires more than unit angular momentum per particle. This directly contradicts the mean-field prediction. Our findings indicate that the many-body system tends to absorb more angular momentum from a perturbation than the mean-field system (see next paragraph).

Different Particle Numbers, Trap Anisotropies, Rotation Frequencies and Interparticle Interactions and Comparison to Mean-Field
Phantom vortices were found in various systems where the particle number N , trap anisotropy η max , rotation frequency ω and interparticle interaction strength λ 0 were varied as shown in Table 1.  In order to estimate the generality of the emergence of fragmentation, we would like to discuss a scan through different particle numbers in detail. We keep all other dynamical parameters fixed, with ω = 0.78 and η max = 0.10, and vary the particle number from N = 10 to N = 10 000 with M = 2 orbitals. The interaction strength g = λ 0 (N − 1) = 17.1 is kept fixed, as it is the single relevant parameter in GP dynamics. We emphasize that GP theory would predict exactly the same dynamics as N is varied with fixed g. Two orbitals are used because this is the minimum required to capture fragmentation. Should fragmentation not occur, the dynamics would be identical for all M ≥ 2, and exactly mimic GP theory.

Parameter Range
We studied systems with up to N = 10 4 and found that fragmentation occurs on relevant time scales (Fig. 1). Importantly, after the trap becomes symmetric, all these calculations show a phantom vortex at the origin in the first orbital, and ghost vortices in both orbitals.
Only ghost vortices are seen in the density. At the end of the rotating procedure, the stable angular momentum per particle ranged from 0.88 to 1.55. This further corroborates the claim that the time-dependent GP approximation systematically underestimates the quantity of angular momentum absorbed from the perturbation. It is important to stress that the phenomenon of phantom vortices exists over a large range of particle numbers also relevant to current experiments.

Comparison to Gross-Pitaevskii
We ran a calculation with M = 1 orbital, with the same one-body and two-body potential, to compare our results directly to the time-dependent GP theory. This check is necessary because the literature does not display the results of subcritical rotations in our parameter regime in detail [11,12]. Furthermore, these treatments used a contact (delta) interaction while we use a Gaussian. We found that the GP system does not nucleate vortices in the density, but does show ghost vortices along the edge of the cloud. This agrees with our manybody treatment. However, the GP treatment results in a much smaller angular momentum L z /N = 0.46 as compared to our computations which yielded L z /N = 1.25 for N = 100 and M = 4. The many-body system absorbs significantly more angular momentum from an external perturbation.

Comparison to contact interaction
In the main manuscript we presented calculations with a short-range interaction, modeled by a Gaussian function. It is natural to ask whether the correlations and the appearance of phantom vortices have been triggered by the finite-range interparticle interaction. The answer is negative. A comparison to a simulation of an identical system but with a con-tact interaction (Dirac-delta pseudopotential rather than a Gaussian) shows that the main characteristics of the system are unchanged and phantom vortices do form in a very similar configuration (Fig. 3). The initial total energies are close to identical in the two cases, whereas the contact-potential system accumulates slightly more energy and angular momentum at the end of the process (Fig. 2). The largest natural occupation number is around 65%, which is higher than its Gaussian counterpart. This shows that the contact interaction system exhibits a smaller fragmentation. The orbital densities are found to have a very similar structure to those in the main text with a Gaussian interaction: the first orbital shows a charge-1 vortex, the second two charge-1 vortices, the third is vortex-free and the last orbital features a large charge-3 vortex. However, since the first occupation is higher, the first phantom vortex is visible as a dip in the total density. c. Video S3: http://youtu.be/gG7dprvRWGg Scan of reference point in g (1) at t = 450: Movie corresponding to Fig. 4 of the main text. At a fixed time, t = 450, we plot g (1) (r|r ′ ) as r ′ is swept from −6 to 6 along the x-axis (panels a,c) and y-axis (panels b,d). We plot the magnitude in panels a,b and the phase, S g , in panels c,d. There are topological defects that are associated with phantom vortices at each value of r ′ . Densities at t = 450.0 (towards the end of the rotation routine) for the system interacting via a contact pseudopotential (i.e. delta function) (cf. Fig. 3 of the main text -all other parameters besides the interaction type are identical). Qualitatively the system behaves in the same fashion as the one originally presented; the density appears vortex-free (first panel left) while phantom vortices have been nucleated on the first, second and fourth natural orbitals of the same charges as in the case of short-range interactions. Notably, in the case of the contact interaction the first natural occupation is higher than its counterpart in the Gaussian system, resulting in a more pronounced dip in the density in the trap center. All quantities shown are dimensionless.