Angular momentum in interacting many-body systems \\ hides in phantom vortices
Storm E. Weiner$^{*
Marios C. Tsatsos
Lorenz S. Cederbaum
Axel U. J. Lode
Vortices are essential to understand angular momentum in quantum systems such as superfluid Helium, ultracold atomic gases, and type-II superconductors. 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 interacting many-body systems is enriched by the emergence of quantum correlations whose description demands a more elaborate method. Herein we theoretically investigate the full many-body dynamics of the acquisition of angular momentum in 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 that, contrary to the conventional 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 vortex modes which have so far evaded experimental detection. This phenomenon is likely important in the formation of the Abrikosov lattice and the onset of turbulence in superfluids.
Weiner_PhantomVortices_2014.tex
\documentclass[preprint]{revtex4}
\usepackage{graphicx}  % needed for figures
\usepackage{dcolumn}   % needed for some tables
\usepackage{bm}        % for math
\usepackage{amssymb}   % for math
\usepackage{sidecap}   % for side caption
\usepackage{doi}%<----------
\usepackage{hyperref}
\usepackage{changes}
\usepackage{amsmath}


\begin{document}



\title{Angular momentum in interacting many-body systems \\ hides in phantom vortices}

\author{Storm E. Weiner$^{*}$\footnote[0]{$^*$
Corresponding author: StormWeiner@Berkeley.edu}}
\affiliation{Department of Chemistry, University of California at Berkeley, CA, USA}
\author{Marios C. Tsatsos}
\affiliation{Instituto de F\'isica de S\~ao Carlos, Universidade de S\~ao Paulo, S\~ao Carlos, S\~ao Paulo, Brazil}                          
\author{Lorenz S. Cederbaum}
\affiliation{Theoretische Chemie, Physikalisch--Chemisches Institut, Universit\"at Heidelberg,
Im Neuenheimer Feld 229, D-69120 Heidelberg, Germany}
\author{Axel U. J. Lode}
\affiliation{Department of Physics, University of Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland}

\date{\today}

\begin{abstract}
 Vortices are essential to understand angular momentum in quantum systems such as superfluid Helium, ultracold atomic gases, and type-II superconductors. 
 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 interacting many-body systems is enriched by the emergence of quantum correlations whose description demands a more elaborate method.
 Herein we theoretically investigate the full many-body dynamics of the acquisition of angular momentum in 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 \textit{phantom vortex} that, contrary to the conventional vortex, can be detected as a topological 
 defect of spatial coherence, but \textit{not} of the density. 
 We describe previously unknown many-body mechanisms of vortex nucleation and show that angular momentum is hidden in phantom vortex modes which have so far evaded experimental detection. This phenomenon is likely important in the formation of the Abrikosov lattice and the onset of turbulence in superfluids.
 \end{abstract}



\pacs{ 03.75.Hh, 05.30.Jp, 03.65.w, 03.75.Nt}
\maketitle


Quantized vortices are perhaps the most interesting way angular momentum is known to manifest in quantum many-body systems \cite{Donnelly,Pethick,Ketterle2001}. Quantized vortices appear in a variety of systems including superfluid Helium \cite{Donnelly}, type-II superconductors \cite{Abrikosov1957}, atomic Bose-Einstein condensates (BECs) \cite{Pethick},  and exciton-polariton condensates \cite{Baas2008} to name but a few. Typically, a quantum vortex is characterized by a density node and a $2\pi$ phase discontinuity. 
Yet, coreless vortices have been recently created in spinor \cite{Lovegrove2014} and multicomponent \cite{Liu2013} condensates, where the core of a vortex in one (spin-)component is filled with the other (spin-)component. However, coreless vortices have not been seen in single-component condensates.
Other known signatures of angular momentum, such as center-of-mass rotation, surface waves, and quadrupole modes for instance, are also features of the density profile~\cite{Dalibard2000,Pethick}. However, it is well known that there exist aspects of many-body dynamics, like fragmentation \cite{James1982,Spekkens1999}, that may not be visible in the density \cite{glauber,Kaspar2008}.

Traditionally, the dynamics of BECs are aptly treated using the time-dependent Gross-Pitaevskii (GP) equation \cite{Dalfovo1999} which assumes that the many-body wavefunction  is coherent and condensed for all time. Naturally, this assumption neglects the existence of quantum correlations and precludes the possibility of fragmentation, i.e.~macroscopic occupation of more than one single particle state \cite{Onsager1956}.  More recently, there has been much effort devoted to exploring the role of fragmentation in stationary states of bosonic systems \cite{James1982,Spekkens1999,Cederbaum2003,Streltsov2006,Mueller2006,Lewenstein2009,Fischer2009,Fischer2014} including those with spin degrees of freedom \cite{Qi2013,Song2014} using methods such as best mean-field theory or general variational principles. In particular, it has been shown that there exist energy eigenstates with definite angular momentum that are fragmented in two dimensions \cite{Marios2014}. Furthermore, in the pioneering study of adiabatic vortex nucleation in few-body systems \cite{Lewenstein2009}, it was shown that fragmentation is unavoidable. Additionally, nontrivial dynamics involve many eigenstates, each of which is generally fragmented. Hence, it stands to reason that a true many-body method is necessary to describe the dynamics of a two-dimensional condensate as it acquires angular momentum.

Such a many-body generalization is the multi-configurational time-dependent Hartree for bosons (MCTDHB) \cite{MCTDHB} method, which has been shown to accurately capture the transition from coherence to fragmentation with great success in a variety of systems \cite{Alexej2007,Axel2012,streltsova2014}.

Here we investigate the dynamics of interacting bosons in two dimensions by numerically solving the time-dependent many-body Schr\"odinger equation for a standard rotation scenario using the recursive implementation of MCTDHB (R-MCTDHB) \cite{RMCTDHB}. To the best of our knowledge, this is the first work which explores the many-body dynamics of this type of system beyond mean-field theory.
The system that we describe here is similar to that simulated in Ref.~\cite{Suominen2003} using the time-dependent GP equation, and analyzed experimentally in Ref.~\cite{Madison2000}. There, regions in the parameter space were established for which vortices are observed in the density. Even though the GP treatment for the herein chosen parameter regime predicts no vortex nucleation (Supplementary Information), the many-body analysis exhibits rich vortex dynamics that are invisible in the density. Our results demonstrate that angular momentum in fragmented condensates manifests itself in a new type of coreless vortex, which we name the \textit{phantom vortex} due to its elusive nature. An analysis of the energy eigenstates of a six-boson system in the corotating frame has shown signatures of phantom vortices (see Fig.~3 in Ref. \cite{Lewenstein2009}). Yet, their importance in the dynamical acquisition of angular momentum by many-body systems, even in parameter regime where no vortices nucleate in the density, has been overlooked. 



\section{System Description}
We begin our analysis, \textit{in silico}, by first computing the many-body ground state of interacting bosons in an isotropic harmonic trap.  Subsequently, we transfer angular momentum gently into the system by elliptically deforming the harmonic trap while rotating the axis of anisotropy. The level of anisotropy is linearly increased to a maximum value, held constant over a fixed time interval, and then linearly decreased back to zero.  The trap is then kept isotropic for the remainder of the simulation. See Fig.~\ref{fig:Lz}a and Supplementary Information for details.

The dynamics of an interacting gas of $N$ bosons are governed by the time-dependent many-body Schr\"odinger equation
\begin{equation}
 i \partial_t \vert \Psi\rangle = \mathcal{H} \vert \Psi\rangle .
 \label{eq:shro}
\end{equation}
MCTDHB solves Eq.~\ref{eq:shro} allowing the bosons to dynamically populate all configurations of $N$ particles in $M$ time-dependent fully variationally optimized single-particle states (Supplementary Information). 
The many-body Hamiltonian reads in dimensionless units \cite{dimensions}:
\begin{equation*}
 \hat{\mathcal{H}}(\mathbf{r}_1,...,\mathbf{r}_N;t)=\sum_{i=1}^{N} \hat{h}_i (\mathbf{r}_i;t)+ \lambda_0 \sum_{i<j}^N\hat{W}(\mathbf{r}_i-\mathbf{r}_j).
\end{equation*}
We use as a two-body potential $\hat{W}$ a normalized Gaussian of width $\sigma=0.25$ and interaction strength $\lambda_0 (N-1)=17.1$. The one-body Hamiltonian $\hat{h}_i(\mathbf{r}_i;t)$ is the sum of the kinetic energy $\hat{T}_i=-\frac{1}{2}\partial_{\mathbf{r}_i}^2$ and the trapping potential,
\begin{equation*}
\hat{V}(\mathbf{r},t)= \frac{1}{2} \left(x(t)^2+y(t)^2\right)+\frac{1}{2} \eta(t) \left(x(t)^2-y(t)^2\right),
\end{equation*}
where 
\begin{equation*}
 \left(
\begin{array}{c}
x(t)\\
y(t)\\
\end{array}
\right) =
\left(
\begin{array}{c c}
cos(\omega t) & sin(\omega t)\\
-sin(\omega t) & cos(\omega t)\\
\end{array}
\right)
\left(
\begin{array}{c}
x(0)\\
y(0)\\
\end{array}
\right)
\end{equation*} 
and $\omega=0.78$. The anisotropy parameter varies in time so that it is linearly ramped up from $0$ to $\eta_{max}=0.1$ over $t_r=80$, held fixed at $\eta_{max}$ for $t_f=220$, linearly ramped down to zero over $t_r$, and then held fixed at $0$ for the rest of simulation as shows in Fig.~\ref{fig:Lz}a. 

The primary object of our analysis is the one-body reduced density matrix  (RDM) \cite{Kaspar2008,glauber}, which is written in its eigenbasis as
\begin{equation*}
 \rho^{(1)}(\textbf{r}\vert \mathbf{r}';t)=\sum_{i=1}^M \rho^{(NO)}_i(t) \phi_i^{*}(\mathbf{r}';t) \phi_i(\textbf{r};t),
\end{equation*}
and its diagonal $\rho(\mathbf{r};t)=\rho^{(1)}(\textbf{r}\vert \mathbf{r};t)$, which is called the density.
The eigenvalues, $\rho^{(NO)}_i(t)$, are ordered in decreasing magnitude and called \textit{natural occupations} and the eigenfunctions, $\phi_i(\mathbf{r};t)$, are called the \textit{natural orbitals}. Note that the density is equal to the sum of the amplitude squared of the natural orbitals, weighted with their natural occupations.
The first order correlation function, $g^{(1)}(\mathbf{r}\vert\mathbf{r}';t)$, is defined in terms of the RDM as
\begin{equation*}
 g^{(1)}(\mathbf{r}\vert\mathbf{r}';t)=\frac{\rho^{(1)}(\mathbf{r}\vert\mathbf{r}';t)}{\sqrt{\rho^{(1)}(\mathbf{r}\vert\mathbf{r};t) \rho^{(1)}(\mathbf{r}'\vert\mathbf{r}';t)}} .
\end{equation*}
A coherent state at time $\tau$ must satisfy $\vert g^{(1)}(\mathbf{r}\vert\mathbf{r}';\tau)\vert=1$, while any other value indicates fragmentation~\cite{glauber}.

\section{Simulation Results}
We have found that there may exist many vortices in the natural orbitals despite there being no visible vortices in the density (Figs.~\ref{fig:115} and \ref{fig:450}) even with angular momentum per particle $>1\hbar$ (Fig.~\ref{fig:Lz}). Hence we name them \textit{phantom vortices}. 
These phantom vortices persist for long time scales compared to the period of the harmonic trap. Phantom vortices that exist for long times near the center of the trap either nucleate on existing topological defects present in the initial natural orbitals or nucleate via a transfer of vorticity between natural orbitals (Video S1).
It is interesting to note that the appearance of vortices in different fragments necessitates the definition of an order parameter for \textit{each} fragment. 
Since the simulation starts with a condensed state and both mechanisms of phantom vortex nucleation rely on the dynamics and occupation of several fragments, they are inherently many-body phenomena that cannot be described by mean-field methods.

We now give the chronological description of the evolution of $N=100$ particles with $M=4$ orbitals starting in the ground state at $t=0$ to $t=500$ (see also Video S1 in Supplementary Information). 
At $t=0$, $\phi_1,...,\phi_4$ resembled 1s, 2p$_{x}$, 2p$_{y}$, and 2s orbitals, respectively.  The system began with $\rho_{1}>99\%$ occupation in $\phi_1$ and thus was almost entirely condensed (Fig.~\ref{fig:Lz}b).  As the anisotropy was ramped up, $\phi_1 ,..., \phi_4$ deformed smoothly and rotated with the potential.  
By $t=29.6$, due to a reordering of the natural occupations, $\phi_3$ and $\phi_4$ switched labels. At this time, $\phi_1$, $\phi_3$, and $\phi_4$ already showed faint hints of vortices at the edge of their density similar to those seen at later times in Figs.~\ref{fig:115} and \ref{fig:450}. These vortices are the fragmented counterparts of \textit{ghost vortices} reported in~\cite{Tsubota2002}.  
By $t=45$, many ghost vortices were established in $\phi_1$ and $\phi_4$. Since the state was still more than $99\%$ condensed, the ghost vortices in $\phi_1$ manifested in the density and would thus have been detectable in high-fidelity absorption imaging.  A one-dimensional cut along a core of a ghost vortex in $\rho(\mathbf{r})$ showed that the outer density maximum was less than $1\%$ of the density maximum at the center of the cloud. By this time, the two initial lobes in $\phi_2$ had spread out and closed off the angular node into an elliptic shape. This was the first sign of a true phantom vortex in the bulk of the cloud (Fig.~\ref{fig:115}c). We term this first mechanism of phantom vortex nucleation \textit{node mutation}, since a node in the fragment deforms and mutates into a phantom vortex.
Around $t=70$, ghost vortices in $\phi_4$ fused with a distorted angular node to split the two lobes into four with an ``I'' shaped node. This node mutated to nucleate three persistent phantom vortices which merged into the charge 3 vortex shown in Fig.~\ref{fig:450}e,j. 

By $t=90$, the system was significantly fragmented,  $\rho^{(NO)}_1 < 90\%$ (Fig.~\ref{fig:Lz}b). Fragmentation obscured the vortex-induced density nodes of $\phi_1$ from being visible in  $\rho(\textbf{r})$ because density from other orbitals filled the vortex cores present in $\phi_1$. We comment that although larger particle numbers delay the onset of fragmentation, we have observed phantom vortices in simulations with up to $N=10^4$ (Supplementary Information).
From $t=90$ to $t=140$, each orbital had a complicated vortex structure while the density $\rho(\textbf{r})$ was largely featureless (Figs.~\ref{fig:115}a and \ref{fig:Lz}a). This marks an important aspect of the dynamics: there were many phantom vortices in each orbital, but remarkably no vortices were detectable by directly observing the density (Figs.~\ref{fig:115} and \ref{fig:450}). In this time interval, the orbital angular momenta, $(L_z)_{ii}\equiv \langle\phi_i \vert L_z\vert\phi_i \rangle$, attained their maximal values and started to decay to their equilibrium values despite the maximal anisotropy.  

After $t=150$, the trap was still maximally anisotropic, yet the system energy and total angular momentum became saturated (Fig.~\ref{fig:Lz}d).
By $t=144$, the vortex structures of $\phi_2$ and $\phi_4$ reached their steady state. $\phi_2$ had a single phantom vortex near the center of the cloud which was nucleated from the single node of its initial shape. This node-mutated phantom vortex persisted for the length of the simulation. At $t=150$, $\phi_4$ had three phantom vortices in a linear arrangement about the center of the cloud, which persisted for the duration of the anisotropy. This triplet fused into a triply charged phantom vortex at the center of the cloud at $t=380$ when the trap became symmetric again (Fig.~\ref{fig:450}e,j). Near $t=220$, $\phi_2$ and $\phi_1$ swapped labels due to occupation reordering (see $\rho^{(NO)}_1$ and $\rho^{(NO)}_2$ in Fig.~\ref{fig:Lz}b), therefore the steady state of the phantom vortex structure can be seen in Fig.~\ref{fig:450}b,e,g,j. From $t=220$ on, the labeling of the orbitals remained fixed because the occupation numbers did not change order anymore (cf. Fig.~\ref{fig:Lz}b).  
 

At $t=220$, there were two corotating vortices near the center of the cloud in $\phi_2$ and none present in $\phi_3$. This phantom vortex pair nucleated at the edge of the cloud and then moved towards the center, resembling mean-field vortex nucleation.  However, in our treatment, this pair was transient. The phantom vortices then transferred from $\phi_2$ to $\phi_3$, marking a second mechanism unique to phantom vortex nucleation: slow orbital-orbital vorticity transfer, which we now describe in detail.
By $t=275$, there were two prominent ghost vortices in $\phi_3$ and the phantom vortex pair in $\phi_2$ had returned from the center to the edge of this fragment.  Gradually, the vortices in $\phi_2$ disappearred by exiting the edge of the orbital while the prominent phantom vortex pair in $\phi_3$ entered the orbital bulk. By $t=340$, the phantom vortex nucleation in $\phi_3$ was complete, marked by the intervortex separation being as small as previously in $\phi_2$ at $t=220$, before the transfer (Video S1).  This direct interaction between $\phi_2$ and $\phi_3$ appears as a strong correlation in the orbital angular momentum, $(L_z)_{22}$ and $(L_z)_{33}$, as seen from $t=250$ to $t=500$ in Fig.~\ref{fig:Lz}c.  The transfer of vorticity between fragments occurred on a time scale much slower ($\tau \approx 50-100$) than all other observed dynamics ($\tau\approx5-10$). 

As the anisotropy was ramped down, the phantom vortex pair in $\phi_3$ approached the trap center and reached a minimal separation of $\approx\hspace{-0.5mm}1.5$ by $t=380$ when the trap became symmetric.  This was the only instance where we observed a phantom vortex nucleate at the edge of a fragment and persist in its bulk (Fig.~\ref{fig:450}d,i). 

Although phantom vortices are not detectable directly from the density, they are strikingly pronounced in the one-body correlation function, $\vert g^{(1)}(\mathbf{r}\vert \mathbf{r}';t)\vert$ (Fig.~\ref{fig:g1}).  We have two key observations about the coherence of phantom vortices. First, we fix $\mathbf{r}'$ at the core of a phantom vortex and see that $\vert g^{(1)}(\mathbf{r}\vert \mathbf{r}';t)\vert  \approx 1 $ for $\mathbf{r}\approx \mathbf{r}'$. Remarkably, $\vert g^{(1)}(\mathbf{r}\vert \mathbf{r}';t )\vert$ is close to $0$ for all $\mathbf{r}$ near and inside the cores of other phantom vortices. Second, we fix $\mathbf{r}'$ away from the cores of all phantom vortices and again observe that $\vert g^{(1)}(\mathbf{r}\vert \mathbf{r}' ;t)\vert$ is close to $0$ for all $\mathbf{r}$ near and inside the cores of phantom vortices (Fig.~\ref{fig:g1}, Videos S2 and S3). Even phantom vortices in the same fragment are found to be incoherent.  We conclude, therefore, that phantom vortices are incoherent both with respect to each other and the remaining bulk density. This implies that they are distinct objects in the fragmented condensate. Furthermore, phantom vortices are observable since $g^{(1)}(\mathbf{r}\vert \mathbf{r}';t)$ is measurable via interference experiments, e.g. Ref.~\cite{Hofferberth2007}. 


\section{Discussion}

Since we analyzed a subcritical rotation, we did not detect any bulk vortices in $\rho(\mathbf{r};t)$. However, such vortices are routinely produced, see e.g. Refs.~\cite{Anderson2008,Fetter2009,Hall2010} and references therein. We found fragmentation well before there was sufficient angular momentum to nucleate a single vortex in $\rho(\mathbf{r};t)$, so it is a natural question to ask whether or not the experimentally observed vortices are coherent or fragmented  objects. If they are fragmented, phantom vortex cores across all relevant orbitals must coincide. Our calculations show that if the coincident phantom vortices are of the same charge there must be additional phantom vortices because the fragments must maintain orthonormality, see for instance $\phi_2$ and $\phi_4$ in Fig.~\ref{fig:115}c,e and Supplementary Information for detailed discussion.  Furthermore, if the only phantom vortices are coincident, they must be of different charge, see for instance $\phi_1$ and $\phi_4$ in Fig.~\ref{fig:450}b,e (see Supplementary Information for discussion). In either case, the angular momentum per particle must be \textit{greater than unity} in order to nucleate a fragmented vortex in $\rho(\mathbf{r};t)$.

In summary, we have observed rich vortex dynamics within individual orbitals of a fragmented single component BEC. Since these vortex dynamics cannot be observed in the density, we termed the vortices in the fragments \textit{phantom vortices}. We have identified two mechanisms of phantom vortex nucleation, namely \textit{node mutation} and \textit{slow orbital-orbital vortex transfer}, which have no mean-field analogue. In node mutation, phantom vortices nucleate on preexisting topological defects in a fragment, whereas in slow orbital-orbital transfer, vortices are transferred between fragments. Phantom vortices are clearly visible in the correlation function. A detailed analysis of the correlation function shows that phantom vortices are completely incoherent \textit{both} with each other \textit{and} the bulk density between other phantom vortices. Phantom vortices are thus distinct quantum objects that are experimentally observable via interference experiments. We suspect that phantom vortices play an important role in the formation of the Abrikosov lattice of type-II superconductors, as well as the onset of turbulence in superfluid Helium and BECs. In particular, our results suggest that there would be phantom vortex dynamics well before flux-vortices are visible in the superconductor. In turn, the existence of phantom vortices in type-II superconductors would prove the existence and importance of fragmentation in these systems.


\input{Weiner_PhantomVortices_2014.bbl}
\clearpage

\acknowledgments{We thank Alexej I. Streltsov, Ofir E. Alon, and Tomos Wells for helpful discussions. Computation time on the Cray XE6 cluster Hermit at the High Performance Computing Center Stuttgart (HLRS), and financial support by the DAAD-RISE program, the Swiss SNF, the NCCR Quantum Science and Technology, DFG, and FAPESP are gratefully acknowledged. We also acknowledge the hospitality of Vanderlei Bagnato and the IFSC-USP, where part of the work was completed.}
\\[5mm]
\indent \textbf{Author Contributions} \\[4mm]
\indent All authors contributed equally to this work.

\begin{figure}[t]
 \includegraphics[width=1\textwidth]{Lz.png}
 \caption{ \textbf{Emergence of fragmentation and accumulation of angular momentum in the dynamics of the rotating two-dimensional Bose gas}:
 (a) Sketch of the anisotropy ramping procedure with plots of $\rho(\mathbf{r})$ at representative times in each part
 (b) The onset of fragmentation occurs around t=80, which corresponds to the time of maximum anisotropy. By the end of the time of maximal anisotropy, the state is completely fragmented.
 (c) Orbital angular momenta reach their maximum values early in the period of maximum anisotropy ($t=80-150$), but evolve to their equilibrium values before the anisotropy is removed.  Discontinuities in $(L_z)_{ii}$ arise from occupation reordering. Fluctuations in $(L_z)_{22}$ and $(L_z)_{33}$ from t=250 to t=500 are strongly correlated due to vorticity transfer between fragments $\phi_2$ and $\phi_3$.
 (d)  Energy and angular momentum are saturated by t=200, while the anisotropy is still maximal. The angular momentum stabilizes at $L_z\approx1.25$ per particle. The strong correlation between energy and $L_z$ indicates that the perturbation strictly excites angular momentum modes.
 To guide the eye, $L_z=1$ on the bottom plot is marked with a horizontal dotted line. All quantities shown are dimensionless.}
 \label{fig:Lz}
\end{figure}



\begin{figure}[thp!]
 \includegraphics[width=1\textwidth]{115.png}
 \caption{\textbf{Phantom vortices manifest in condensate fragments but are obscured in the density.}
 (a) The density $\rho(\mathbf{r})$ shows some density minima due to phantom vortices in $\phi_1$, but no true density node.
 (b)-(e) The natural orbital densities, $\vert \phi_i (\mathbf{r})\vert ^2$ ($i=1,2,3,4$), are plotted.  There are many phantom vortices present in each orbital.
 (f) The phase, $S_g(\mathbf{r}\vert0)$, of $g^{(1)}(\mathbf{r}\vert0)$ is plotted. Note, $S_{g}$ is not the many-body phase, $S_{\text{MB}}(\mathbf{r}_1,...,\mathbf{r}_N;t)$, which is too complicated to visualize. Since $\phi_1$ carries most of the particles at this time, $S_g$ bears strong resemblance to $S_1$.
 (g)-(j) The phases, $S_i$, of the natural orbitals $\phi_i$ ($i=1,2,3,4$) are plotted. Each $2\pi$ phase discontinuity marks a phantom vortex core.
 The central phantom vortex in $\phi_2$ (c),(h) was mutated from its initial angular node, and persisted for the duration of simulation.  The two central phantom vortices in $\phi_3$ (d),(i) were mutated from its initial angular nodes and are transient. The three centermost phantom vortices in $\phi_4$ (e),(j) were mutated from an ``I'' shaped node and persisted for the length of the simulation.  All other phantom vortices nucleated at the edge of their orbital density and are transient. 
 All panels are plotted at t=115.0, when $\rho^{(NO)}_1=82.0\%$, $\rho^{(NO)}_2=11.8\%$, $\rho^{(NO)}_3=4.1\%$, and $\rho^{(NO)}_4=2.1\%$. See complementary Video S1 in Supplementary Information. All quantities shown are dimensionless.}
 \label{fig:115}
\end{figure}

 \begin{figure}[!]
  \includegraphics[width=1\textwidth]{450.png}
  \caption{\textbf{Persistence of phantom vortices in the stationary trap.} 
  (a) The density $\rho(\mathbf{r})$ shows a density minimum at the origin which is $75\%$ the maximum density.
  (b)-(e) The natural orbital densities, $\vert \phi_i (\mathbf{r})\vert ^2$ ($i=1,2,3,4$), are plotted. 
  (f) The phase, $S_g(\mathbf{r}\vert0)$, of $g^{(1)}(\mathbf{r}\vert0)$ is plotted. Note, $S_{g}$ is not the many-body phase, $S_{\text{MB}}(\mathbf{r}_1,...,\mathbf{r}_N;t)$, which is too complicated to visualize. There are singularities in $S_g$  that are associated with phantom vortices which do not correspond to features in $\rho$.
  (g)-(j) The phases, $S_i$, of the natural orbitals $\phi_i$ ($i=1,2,3,4$) are plotted. Each $2\pi$ phase discontinuity marks a phantom vortex core. 
  The depicted phantom vortex configuration is stable for more than 100 trap periods. The phantom vortices in $\phi_1$ (b),(g) and $\phi_4$ (e),(j) nucleated along preexisting topological defects through node mutation whereas the pair of phantom vortices in $\phi_3$ (d),(i) were nucleated by slow orbital-orbital vortex transfer as ghost vortices which then moved to the center. The nucleation of the phantom vortex pair was coupled to the destruction of a pair in $\phi_2$ (c),(h)  See the coupled oscillation of $(L_z)_{22}$ and $(L_z)_{33}$ (Fig.~\ref{fig:Lz}c).  
  All panels are plotted at t=450.0, when $\rho^{(NO)}_1=40.8\%$, $\rho^{(NO)}_2=25.0\%$, $\rho^{(NO)}_3=20.8\%$, and $\rho^{(NO)}_4=13.5\%$.
  See complementary Video S1 in Supplementary Information. All quantities shown are dimensionless.}
 \label{fig:450}
 \end{figure}



\begin{figure}
 \includegraphics[width=1\textwidth]{g1.png}
 \caption{\textbf{Phantom vortices are exposed in the correlation function.} 
 To visualize the four-dimensional single-particle correlation function, $g^{(1)}(\textbf{r}\vert\textbf{r}')$, we fix a reference point at $\textbf{r}'=(0,0)$ (a),(c) and $\textbf{r}'=(1.25,0)$ (b),(d). The function is complex, so we plot the magnitudes (a),(b) and phases (c),(d) separately.  In all panels, we mark the reference point, $\mathbf{r}'$, with a white circle.
 In panels (a) and (c) [(b) and (d)], the $\mathbf{r}'$ is colocated with a phantom vortex core in $\phi_2$ [$\phi_1$]. 
 In both cases, phantom vortex cores in all orbitals appear as spots of complete incoherence, $|g^{(1)}|\approx 0$, while the core colocated with the reference point has a $|g^{(1)}|\approx 1$ (full coherence). All panels are plotted at t=219.4, when $\rho^{(NO)}_1=35.3\%$, $\rho^{(NO)}_2=35.2\%$, $\rho^{(NO)}_3=17.5\%$, and $\rho^{(NO)}_4=12.0\%$. All quantities shown are dimensionless.}
 \label{fig:g1}
\end{figure}


\clearpage
\renewcommand\thefigure{S\arabic{figure}}
\renewcommand\theequation{S\arabic{equation}}

\section{Supplementary Information}
\setcounter{figure}{0}
\setcounter{equation}{0}
In this supplementary document, we elaborate on the model we used to study the accumulation of angular momentum in a bosonic system. We define the MCTDHB method we use to solve the time-dependent many-body Schr\"odinger 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.
\subsection{Model}
We solve the time-dependent Schr\"odinger equation,
\begin{equation}
 i\partial_t \vert \Psi\rangle = \hat{\mathcal{H}} \vert\Psi\rangle,
 \label{eq:ShroSUP}
\end{equation}
for $N$ interacting bosons using the dimensionless~\cite{dimensions} Hamiltonian
\begin{equation}
 \hat{\mathcal{H}}(\mathbf{r}_1,...,\mathbf{r}_N;t)=\sum_{i=1}^N \hat{h}_i(\mathbf{r}_i;t)+\lambda_0\sum_{i<j}^N \hat{W}(\mathbf{r}_i-\mathbf{r}_j),
 \label{eq:hamiltonian}
\end{equation}
where $\hat{h}_i(\mathbf{r};t)$ is the one-particle Hamiltonian
\begin{equation*}
 \hat{h}_i(\mathbf{r}_i;t)=-\frac{1}{2}\partial^2_{\mathbf{r}_i}+V(\mathbf{r}_i;t).
\end{equation*}
The trapping potential $V(\mathbf{r};t)$ we use in the present work,
\begin{equation*}
 V(x,y;t)=\frac{1}{2} (x^2+y^2)+\frac{\eta(t)}{2}(x^2-y^2),
\end{equation*}
depends on time through the anisotropy parameter $\eta(t)$ and the rotating coordinates
\begin{equation*}
 \left(
\begin{array}{c}
x(t)\\
y(t)\\
\end{array}
\right) =
\left(
\begin{array}{c c}
cos(\omega t) & sin(\omega t)\\
-sin(\omega t) & cos(\omega t)\\
\end{array}
\right)
\left(
\begin{array}{c}
x(0)\\
y(0)\\
\end{array}
\right).
\end{equation*} 
In our simulations, we set $\omega=0.78$ and $\eta_{max}=0.1$.  We control the anisotropy, $\eta(t)$, according to the following scheme:
 \begin{equation}
   \eta(t) = \left\{
     \begin{array}{lcr}
       \eta_{max}\frac{t}{t_r} & :& t < t_{r} \\
       \eta_{max} & :& t_{r} < t < t_{r} + t_{f} \\
       \eta_{max}-\eta_{max}\frac{t-t_{r}-t_{f}}{t_r}  & :& t_{r}+t_{f} < t < 2t_{r}+t_{f} \\
       0  & :& 2t_{r}+t_{f} < t ,\\
     \end{array} 
   \right.
   \label{eq:eta}
\end{equation}
where $t_r=80$ and $t_f=220$. In this scheme, Eq.~\ref{eq:eta}, we linearly ramp up $\eta(t)$ from $0$ to $\eta_{max}$ over a ramping time $t_r$, and hold it constant at $\eta_{max}$ for a time $t_f$. Subsequently, we linearly ramp $\eta(t)$ back down to $0$ in the same time $t_r$, and then hold it constant at $0$. The initial and final trapping potentials are thus isotropic (see also Fig.~\ref{fig:Lz}a).

For the two-body interaction, we use a normalized Gaussian,
\begin{equation*}
 \hat{W}(\mathbf{r})=\frac{1}{\pi\sigma^2} e^{-\frac{r^2}{2\sigma^2}},
\end{equation*}
with characteristic width $\sigma=0.25$. We use a Gaussian rather than the typical contact interaction because the contact  interaction has convergence issues in more than one dimension \cite{rosti2013}.

\subsection{MCTDHB Method and Numerical Details}
We solve the time-dependent Schr\"odinger equation (Eq.~\ref{eq:ShroSUP}) with the Hamiltonian, Eq.~\ref{eq:hamiltonian}, using the MCTDHB method \cite{MCTDHB}. MCTDHB uses the ansatz,

\begin{equation}
 \vert \Psi \rangle = \sum_{\{\mathbf{n}\}} C_\mathbf{n} (t) \vert \mathbf{n};t \rangle
 \label{eq:anzatz}
\end{equation}
where $\vert \mathbf{n};t \rangle=\vert n_1,...,n_M;t\rangle$ is a permanent of $M$  orthonormal orbitals  which \textit{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 ~\cite{Kramer1981}. 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 $\frac{(M+N-1)!}{N!(M-1)!}$ permanents.
In the limit $M\rightarrow \infty$, the configurations $\vert \mathbf{n};t\rangle$ span the complete $N-$particle Hilbert space and the expansion in Eq.~\ref{eq:anzatz} becomes formally exact. In practice, a small number of time-dependent single-particle states is enough to obtain numerically exact results~\cite{Lode2012,AxelThesis}.  In the case of $M=1$,  Eq.~\ref{eq:anzatz} 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~\cite{MCTDHB}. 

We use the recently developed recursive implementation of the MCTDHB
theory \cite{MCTDHB}, the R-MCTDHB package \cite{RMCTDHB} for the
simulations presented in this work. Computations using R-MCTDHB 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 \cite{MCTDH-rev}. For the present
work, plane waves were used as primitive basis.

To check convergence of our computations with respect to the DVR we
computed the energy of the ground state of $N=100$ bosons in a parabolic trap
when varying the number of
DVR points used to represent the simulated domain
of extension $[-8,8] \times [-8,8]$ at different widths of the short-range gaussian potential, $\sigma \in [0.1,0.35]$.
We varied the number of DVR points from $64\times 64$ to $512 \times 512$.
We found that the energy was converged up to the $16$ significant digits with respect to the DVR
using $\sigma=0.25$ and $128 \times 128$ grid points. This
choice of width is well in the range of widths found to match the
physics of contact interaction potentials in Ref.~\cite{rosti2013}.

The bulk of our simulations considers $N=100$ and $M=4$, which yields
$N_C=176~851$ configurations. The integration of the
MCTDHB equations of motion comprises the solution of two coupled sets
of generally non-linear integro-differential equations of motion, see
\cite{MCTDHB,AxelThesis} for details. The set of coefficient
equations is linear and of dimension $N_C$ and the set of orbital
equations is non-linear, integro-differential and of dimension
$M\times128\times128$. It is important to note that the number of
configurations in the mean-field case, $M=1$, is $N_C=1$. Since the
configurations are mutually orthonormal functions, the number of
coefficients directly measures the size of the Hilbert space in
comparison to the mean-field case.
In dimension-full units, the simulation time $t=500.0$ corresponds to
$1.08$ seconds and the domain $[-8,8]\times [-8,8]$ corresponds to
$8\times 8$ micrometers when considering \cite{dimensions}.

\subsection{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, $\rho(\mathbf{r};t)$, is given as the diagonal of the one-body reduced density matrix:
\begin{equation*}
 \rho(\mathbf{r};t)=\rho^{(1)}(\mathbf{r}\vert \mathbf{r};t)
\end{equation*}

The one-body reduced density matrix, $\rho^{(1)}$, for an $N$-body system is defined as the partial trace of the $N$-body density~\cite{Kaspar2008}:
\begin{equation*}
\begin{split}
 \rho^{(1)}(\mathbf{r}_1\vert\mathbf{r}'_1;t) & = N \int \Psi(\mathbf{r}_1,\mathbf{r}_{2},...,\mathbf{r}_N;t) \\
 \times & \Psi^*(\mathbf{r}'_1,\mathbf{r}_2,...,\mathbf{r}_N;t)d\mathbf{r}_{2}...\mathbf{r}_N,
\end{split}
\end{equation*}
where the $\Psi$ is assumed to be normalized to 1.  
The one-body density matrix can be expanded in a complete basis $\{\xi_j(\textbf{r};t)\}$,
\begin{align*}
 \rho^{(1)}(\mathbf{r}\vert\mathbf{r}';t) & =\sum_{j,k} \rho_{jk}(t)\xi_j^*(\mathbf{r};t) \xi_k(\mathbf{r}';t) \\
                                          & =\sum_{j} \rho^{(NO)}_j(t) \phi_j(\mathbf{r};t)\phi_j(\mathbf{r}';t),
\end{align*}
where we have diagonalized the matrix elements $\rho_{kj}$ in the last line to obtain the eigenvalues $\rho^{(NO)}_j(t)$ and eigenfunctions $\{\phi_j(\mathbf{r};t)\}$, which are called the natural occupations and natural orbitals, respectively.
The $1^{st}$-order correlation function, $g^{(1)}$, is defined in terms of the one-body density matrix~\cite{glauber}:
\begin{equation*}
 g^{(1)}(\mathbf{r},\vert\mathbf{r}';t)=\frac{ \rho^{(1)}(\mathbf{r}\vert\mathbf{r}';t)}{\sqrt{\rho^{(1)}(\mathbf{r}\vert\mathbf{r};t) \rho^{(1)}(\mathbf{r}'\vert\mathbf{r}';t)  }}.
\end{equation*}
The off-diagonal, $\mathbf{r} \neq \mathbf{r}'$, terms in $g^{(1)}$ and the orbitals $\phi_i$ are complex-valued.  Thus it is instructive to plot their magnitudes and phases separately. We hence define the phase of $g^{(1)}(\mathbf{r}\vert\mathbf{r}';t)$ to be 
\begin{equation*}
S_{g}(\mathbf{r}\vert\mathbf{r}';t)\equiv \text{arg}\left[g^{(1)}(\mathbf{r}\vert\mathbf{r}';t)\right],
\end{equation*}
and the orbital phases to be
\begin{equation*}
 S_i(\mathbf{r};t)=\text{arg}\left[\phi_i(\mathbf{r};t) \right].
\end{equation*}
For the sake of completeness, we likewise define the phase of the many-body wavefunction, $\Psi$, as
\begin{equation*}
 S_{MB}(\mathbf{r}_1,...,\mathbf{r}_N;t)=\text{arg}\left[\Psi(\mathbf{r}_1,...,\mathbf{r}_N;t)\right].
\end{equation*}
We choose to plot all phases on the interval $[-\pi,\pi]$, but any interval of the same length is equally valid. Furthermore, the global phase transformation, 
\begin{align*}
\phi_i(\mathbf{r};t) & \rightarrow e^{i \theta_i(t)} \phi_i(\mathbf{r};t) \\
C_{\mathbf{n}} & \rightarrow e^{-i(\gamma+  \sum_{i=1}^M \theta_i(t)n_i)} C_{\mathbf{n}} & \forall \gamma,\theta_i\in\mathbb{R}
\end{align*}
leaves the many-body wavefunction (Eq.~\ref{eq:anzatz}) invariant.  Phase shifting an orbital $\phi_i$ by $\theta_i$ is equivalent to plotting it in the interval $[-\pi+\theta_i,\pi+\theta_i]$ with all structures unchanged.  Whereas, $g^{(1)}(\mathbf{r}\vert\mathbf{r}')\propto \rho^{(1)}(\mathbf{r}\vert \mathbf{r}')=\sum_{i=1}^M \rho^{(NO)}_i \vert\phi_i(\mathbf{r}) \vert \vert\phi_i(\mathbf{r'}) \vert e^{i(S_i(\mathbf{r})-S_i(\mathbf{r}'))}$ 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.\ref{fig:115} and \ref{fig:450}, Videos S1 and S2) and $g^{(1)}$ (Fig.\ref{fig:g1} and Video S3) are unaffected by global phase transformations.



\subsection{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
\begin{equation*}
 \phi_j(\mathbf{r})=A_j (r)e^{i q_j\theta },
\end{equation*}
where $A_j(r)\in \mathbb{R}$ is the radial amplitude and $q_j \in \mathbb{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$, $ A(r\rightarrow \infty)=0$ and $A(r \neq 0)\neq 0 $. Thus, either $A_j(r) \leq 0 $ or $A_j(r) \geq 0$ for all $r$. Orthonormality of the orbitals then requires 
\begin{equation*}
 \int_{0}^{\infty} A_j(r) A_k(r) dr \int_0^{2\pi} e^{i (q_k-q_j) \theta }d\theta = 0, \hspace{2.5mm}\forall j\neq k.
\end{equation*}
The radial integral is strictly nonzero since the radial amplitudes cannot change sign. Thus, the angular integral must be zero. Using the orthonormality of Fourier modes, we conclude that $q_k\neq q_j$, i.e. the different phantom vortices must have different charges. This case is partially exemplified in $\phi_1$ and $\phi_4$ (Fig.~\ref{fig:450}b,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 $\phi_2$ or $\phi_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, $\phi_1$ and $\phi_4$, there would be a density vortex.

The other possibility to observe a fragmented vortex in the density is to have coincident phantom vortices of the \textit{same} charge as well as noncoincident phantom vortices in regions where $A_j$ is large enough to ensure orthonormality. For example, $\phi_2$ and $\phi_3$ (Fig.~\ref{fig:115}c,d) 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 \textit{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).

\subsection{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 \cite{Suominen2003,Lewenstein2009}. 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 many-body 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.
\subsection{Different Particle Numbers}
To assess the generality of our findings about the formation of phantom vortices during the acquisition of angular momentum in a rotating gas of bosons, we vary the particle number while keeping all other dynamical parameters fixed. We vary the particle number from $N=10$ to $N=10~000$ with $M=2$ orbitals. The interaction strength $g=\lambda_0(N-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\geq 2$, and exactly mimic GP theory. We found that fragmentation does occur on relevant time scales for up to $N=10^4$ (Fig.~\ref{fig:frag}). Surprisingly, when $N=2\times10^n$, $L_z/N$ is consistently smaller than for other $N$ investigated.  Curiously, these $N$ values also fragment at later times (Fig.~\ref{fig:frag}). 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.


\subsection{Movies}
\paragraph{Video S1:}
\url{http://youtu.be/ezbdLWvSbBI}

\textbf{Time evolution of the density, natural orbitals, and $S_g$:}
(a) shows the density $\rho(\mathbf{r};t)$, (b)-(e) and (g)-(j) show the squared magnitudes and phases of the natural orbitals, respectively, and (f) shows the phase of the one-body correlation function, $S_g(\mathbf{r}\vert0;t)$.  After the onset of fragmentation at $t=80$, the individual orbitals show a complicated vortex structure while the density is largely featureless. $S_g(\mathbf{r}\vert0;t)$ shows phase discontinuities characteristic of vortices even when no vortices are present in $\rho$.

\paragraph{Video S2:} \url{http://youtu.be/whRL8haF4RA}

\textbf{Scan of reference point in $g^{(1)}$ at $t=115$:}
At a fixed time, $t=115$, we plot $g^{(1)}(\textbf{r}\vert\textbf{r}')$ as $\textbf{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 $\textbf{r}'$.

\paragraph{Video S3:} \url{http://youtu.be/gG7dprvRWGg}

\textbf{Scan of reference point in $g^{(1)}$ at $t=450$:}
At a fixed time, $t=115$, we plot $g^{(1)}(\textbf{r}\vert\textbf{r}')$ as $\textbf{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 $\textbf{r}'$.



\begin{figure}
 \includegraphics[width=\textwidth]{FragTime.png}
 \caption{\textbf{Onset of fragmentation at different particle number:} We repeated our calculations varying the particle number, but keeping the mean-field coupling strength $g=\lambda_0(N-1)$ constant with $M=2$ orbitals. Generally, higher particle numbers fragment at later times. Surprisingly, particle numbers that are not multiples of 10 (we tested $2\times10^n$ and $5\times10^n$) are particularly resilient to fragmentation. However, all particle numbers tested up to $N=10^4$, do fragment during the anisotropy ramp-up or during the period of maximum anisotropy.}
 \label{fig:frag}
\end{figure}



\end{document}

Weiner_PhantomVortices_2014.bbl
\begin{thebibliography}{10}
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem{Donnelly}
\bibinfo{author}{Donnelly, R.~J.}
\newblock \textit{\bibinfo{title}{Quantized Vortices in Helium II}}.
\newblock No. \bibinfo{number}{2} in \bibinfo{series}{Cambridge Studies in
  American Literature and Culture} (\bibinfo{publisher}{Cambridge University
  Press}, \bibinfo{year}{1991}).

\bibitem{Pethick}
\bibinfo{author}{Pethick, C.} \& \bibinfo{author}{Smith, H.}
\newblock \textit{\bibinfo{title}{Bose-Einstein Condensation in Dilute Gases}}
  (\bibinfo{publisher}{Cambridge University Press}, \bibinfo{year}{2002}).

\bibitem{Ketterle2001}
\bibinfo{author}{Abo-Shaeer, J.~R.}, \bibinfo{author}{Raman, C.},
  \bibinfo{author}{Vogels, J.~M.} \& \bibinfo{author}{Ketterle, W.}
\newblock \bibinfo{title}{Observation of vortex lattices in Bose-Einstein
  condensates}.
\newblock \textit{\bibinfo{journal}{Science}}
  \textbf{\bibinfo{volume}{\textbf{292}}}, \bibinfo{pages}{476}
  (\bibinfo{year}{2001}).

\bibitem{Abrikosov1957}
\bibinfo{author}{Abrikosov, A.~A.}
\newblock \bibinfo{title}{On the magnetic properties of superconductors of the
  second group}.
\newblock \textit{\bibinfo{journal}{Sov. Phys. JETP}}
  \textbf{\bibinfo{volume}{\textbf{5}}} (\bibinfo{year}{1957}).

\bibitem{Baas2008}
\bibinfo{author}{Lagoudakis, K.~G.} \textit{et~al.}
\newblock \bibinfo{title}{Quantized vortices in an exciton-polariton
  condensate}.
\newblock \textit{\bibinfo{journal}{Nature Phys.}}
  \textbf{\bibinfo{volume}{\textbf{4}}}, \bibinfo{pages}{706}
  (\bibinfo{year}{2008}).
  
 \bibitem{Lovegrove2014}
\bibinfo{author}{Lovegrove, J.}, \bibinfo{author}{Borgh, M.~O.} \&
  \bibinfo{author}{Ruostekoski, J.}
\newblock \bibinfo{title}{Energetic stability of coreless vortices in spin-1
  Bose-Einstein condensates with conserved magnetization}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. Lett.}}
  \textbf{\bibinfo{volume}{\textbf{112}}}, \bibinfo{pages}{075301}
  (\bibinfo{year}{2014}).

\bibitem{Liu2013}
\bibinfo{author}{Liu, C.~F.}, \bibinfo{author}{Fan, H.}, \bibinfo{author}{Gou,
  S.~C.} \& \bibinfo{author}{Liu, W.~M.}
\newblock \bibinfo{title}{Crystallized and amorphous vortices in rotating
  atomic-molecular Bose-Einstein condensates}.
\newblock \textit{\bibinfo{journal}{Sci. Rep.}}
  \textbf{\bibinfo{volume}{\textbf{4}}}, \bibinfo{pages}{4224} (\bibinfo{year}{2013}).

\bibitem{Dalibard2000}
\bibinfo{author}{Chevy, F.}, \bibinfo{author}{Madison, K.~W.} \&
  \bibinfo{author}{Dalibard, J.}
\newblock \bibinfo{title}{Measurement of the angular momentum of a rotating
  Bose-Einstein condensate}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. Lett.}}
  \textbf{\bibinfo{volume}{\textbf{85}}}, \bibinfo{pages}{2223}
  (\bibinfo{year}{2000}).

\bibitem{Spekkens1999}
\bibinfo{author}{Spekkens, R.~W.} \& \bibinfo{author}{Sipe, J.~E.}
\newblock \bibinfo{title}{Spatial fragmentation of a Bose-Einstein condensate
  in a double-well potential}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{\textbf{59}}}, \bibinfo{pages}{3868}
  (\bibinfo{year}{1999}).

\bibitem{James1982}
\bibinfo{author}{Nozi\`eres, P.} \& \bibinfo{author}{James, D.~S.}
\newblock \bibinfo{title}{Particle vs. pair condensation in attractive Bose
  liquids}.
\newblock \textit{\bibinfo{journal}{J. Phys. France}}
  \textbf{\bibinfo{volume}{\textbf{43}}}, \bibinfo{pages}{1133}
  (\bibinfo{year}{1982}).  
  


\bibitem{glauber}
\bibinfo{author}{Glauber, R.~J.}
\newblock \textit{\bibinfo{title}{Quantum Theory of Optical Coherence}}
  (\bibinfo{publisher}{Wiley-VCH Verlag GmbH \& Co. KGaA},
  \bibinfo{year}{2007}).

\bibitem{Kaspar2008}
\bibinfo{author}{Sakmann, K.}, \bibinfo{author}{Streltsov, A.~I.},
  \bibinfo{author}{Alon, O.~E.} \& \bibinfo{author}{Cederbaum, L.~S.}
\newblock \bibinfo{title}{Reduced density matrices and coherence of trapped
  interacting bosons}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{\textbf{78}}}, \bibinfo{pages}{023615}
  (\bibinfo{year}{2008}).

\bibitem{Dalfovo1999}
\bibinfo{author}{Dalfovo, F.}, \bibinfo{author}{Giorgini, S.},
  \bibinfo{author}{Pitaevskii, L.~P.} \& \bibinfo{author}{Stringari, S.}
\newblock \bibinfo{title}{Theory of Bose-Einstein condensation in trapped
  gases}.
\newblock \textit{\bibinfo{journal}{Rev. Mod. Phys.}}
  \textbf{\bibinfo{volume}{\textbf{71}}}, \bibinfo{pages}{463}
  (\bibinfo{year}{1999}).

\bibitem{Onsager1956}
\bibinfo{author}{Penrose, O.} \& \bibinfo{author}{Onsager, L.}
\newblock \bibinfo{title}{Bose-Einstein condensation and liquid helium}.
\newblock \textit{\bibinfo{journal}{Phys. Rev.}}
  \textbf{\bibinfo{volume}{\textbf{104}}}, \bibinfo{pages}{576}
  (\bibinfo{year}{1956}).



\bibitem{Cederbaum2003}
\bibinfo{author}{Cederbaum, L.~S.} \& \bibinfo{author}{Streltsov, A.~I.}
\newblock \bibinfo{title}{Best mean-field for condensates}.
\newblock \textit{\bibinfo{journal}{Phys. Lett. A}}
  \textbf{\bibinfo{volume}{\textbf{318}}}, \bibinfo{pages}{564}
  (\bibinfo{year}{2003}).

\bibitem{Streltsov2006}
\bibinfo{author}{Streltsov, A.~I.,} , \bibinfo{author}{Alon, O.~I.} \& \bibinfo{author}{Cederbaum, L.~S.}
\newblock \bibinfo{title}{General variational many-body theory with complete self-consistency for trapped bosonic systems}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{73}}, \bibinfo{pages}{063626}
  (\bibinfo{year}{2006})
  
  
\bibitem{Mueller2006}
\bibinfo{author}{Mueller, E.~J.}, \bibinfo{author}{Ho, T.-L.},
  \bibinfo{author}{Ueda, M.} \& \bibinfo{author}{Baym, G.}
\newblock \bibinfo{title}{Fragmentation of Bose-Einstein condensates}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{\textbf{74}}}, \bibinfo{pages}{033612}
  (\bibinfo{year}{2006}).

\bibitem{Lewenstein2009}
\bibinfo{author}{Dagnino, D.}, \bibinfo{author}{Barber\'an, N.},
  \bibinfo{author}{Lewenstein, M.} \& \bibinfo{author}{Dalibard, J.}
\newblock \bibinfo{title}{Vortex nucleation as a case study of symmetry
  breaking in quantum systems}.
\newblock \textit{\bibinfo{journal}{Nature Phys.}}
  \textbf{\bibinfo{volume}{\textbf{5}}}, \bibinfo{pages}{431}
  (\bibinfo{year}{2009}).

\bibitem{Fischer2009}
\bibinfo{author}{Bader, P.} \& \bibinfo{author}{Fischer, U.~R.}
\newblock \bibinfo{title}{Fragmented many-body ground states for scalar bosons
  in a single trap}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. Lett.}}
  \textbf{\bibinfo{volume}{\textbf{103}}}, \bibinfo{pages}{060402}
  (\bibinfo{year}{2009}).

\bibitem{Fischer2014}
\bibinfo{author}{Kang, M.-K.} \& \bibinfo{author}{Fischer, U.~R.}
\newblock \bibinfo{title}{Revealing single-trap condensate fragmentation by
  measuring density-density correlations after time-of-flight}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. Lett.}} \textbf{\bibinfo{volume}{113}}, \bibinfo{pages}{140404}
  (\bibinfo{year}{2014}).

\bibitem{Qi2013}
\bibinfo{author}{Zhou, Q.} \& \bibinfo{author}{Cui, X.}
\newblock \bibinfo{title}{Fate of a Bose-Einstein condensate in the presence of
  spin-orbit coupling}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. Lett.}}
  \textbf{\bibinfo{volume}{\textbf{110}}}, \bibinfo{pages}{140407}
  (\bibinfo{year}{2013}).

\bibitem{Song2014}
\bibinfo{author}{Song, S.~W.}, \bibinfo{author}{Zhang, Y.~C.},
  \bibinfo{author}{Zhao, H.}, \bibinfo{author}{Wang, X.} \&
  \bibinfo{author}{Liu, W.~M.}
\newblock \bibinfo{title}{Fragmentation of spin-orbit-coupled spinor
  Bose-Einstein condensates}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{\textbf{89}}}, \bibinfo{pages}{063613}
  (\bibinfo{year}{2014}).
  
\bibitem{Marios2014}
\bibinfo{author}{Tsatsos, M.~C.}
\newblock \bibinfo{title}{Attractive Bose gas in two dimensions: An analytical
  study of its fragmentation and collapse}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{\textbf{89}}}, \bibinfo{pages}{043604}
  (\bibinfo{year}{2014}).  

\bibitem{MCTDHB}
\bibinfo{author}{Alon, O.~E.}, \bibinfo{author}{Streltsov, A.~I.} \&
  \bibinfo{author}{Cederbaum, L.~S.}
\newblock \bibinfo{title}{Multiconfigurational time-dependent Hartree method
  for bosons: Many-body dynamics of bosonic systems}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{\textbf{77}}}, \bibinfo{pages}{033613}
  (\bibinfo{year}{2008}).

\bibitem{Alexej2007}
\bibinfo{author}{Streltsov, A.~I.}, \bibinfo{author}{Alon, O.~E.} \&
  \bibinfo{author}{Cederbaum, L.~S.}
\newblock \bibinfo{title}{Role of excited states in the splitting of a trapped
  interacting Bose-Einstein condensate by a time-dependent barrier}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. Lett.}}
  \textbf{\bibinfo{volume}{\textbf{99}}}, \bibinfo{pages}{030402}
  (\bibinfo{year}{2007}).


\bibitem{Axel2012}
\bibinfo{author}{Lode, A.~U.~J.}, \bibinfo{author}{Streltsov, A.~I.},
  \bibinfo{author}{Sakmann, K.}, \bibinfo{author}{Alon, O.~E.} \&
  \bibinfo{author}{Cederbaum, L.~S.}
\newblock \bibinfo{title}{How an interacting many-body system tunnels through a
  potential barrier to open space}.
\newblock \textit{\bibinfo{journal}{Proc. Natl. Acad. Sci.}}
  \textbf{\bibinfo{volume}{\textbf{109}}}, \bibinfo{pages}{13521}
  (\bibinfo{year}{2012}).

\bibitem{streltsova2014}
\bibinfo{author}{Streltsova, O.~I.}, \bibinfo{author}{Alon, O.~E.}, \bibinfo{author}{Cederbaum, L.~S.}, \& \bibinfo{author}{Streltsov, A.~I.}
\newblock \bibinfo{title}{Generic regimes of quantum many-body dynamics of trapped bosonic systems with strong repulsive interactions}
\newblock \textit{\bibinfo{journal}{Phys. Rev. A (R)}}
  \textbf{\bibinfo{volume}{\textbf{89}}}, \bibinfo{pages}{061602}
  (\bibinfo{year}{2014})
  
  
\bibitem{RMCTDHB}
\bibinfo{author}{Lode, A.~U.~J.} \& \bibinfo{author}{Tsatsos, M.~C.}
\newblock \bibinfo{title}{The recursive multiconfigurational time-dependent
  Hartree for bosons package, version 1.0}.
\newblock \bibinfo{howpublished}{http://ultracold.org; http://r-mctdhb.org; http://schroedinger.org} (\bibinfo{year}{2014}).

\bibitem{Suominen2003}
\bibinfo{author}{Lundh, E.}, \bibinfo{author}{Martikainen, J.-P.} \&
  \bibinfo{author}{Suominen, K.-A.}
\newblock \bibinfo{title}{Vortex nucleation in Bose-Einstein condensates in
  time-dependent traps}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{\textbf{67}}}, \bibinfo{pages}{063604}
  (\bibinfo{year}{2003}).

\bibitem{Madison2000}
\bibinfo{author}{Madison, K.~W.}, \bibinfo{author}{Chevy, F.},
  \bibinfo{author}{Wohlleben, W.} \& \bibinfo{author}{Dalibard, J.}
\newblock \bibinfo{title}{Vortex formation in a stirred Bose-Einstein
  condensate}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. Lett.}}
  \textbf{\bibinfo{volume}{\textbf{84}}}, \bibinfo{pages}{806}
  (\bibinfo{year}{2000}).

\bibitem{dimensions}
  \bibinfo{year}{To convert to dimensionful units, the Hamiltonian is
  multiplied by $\frac{\hbar^2}{m L^2}$. With the mass
  $m=1.44316\times10^{-25}$Kg corresponding to $Rb^{87}$, and length scale
  $L=0.750\mu$m, we get a transverse trapping frequency of $\omega_0=1.30 $kHz and a time scale of $4.84 $ms.
  With this choice of scale, $t=500$ corresponds to $2.42 s$, and the linear extent of the simulation is $12 \mu$m. The interaction parameter $\lambda_0$ is related to the
  scattering length $a_s$ and transverse confinement $l_z=\sqrt{\frac{\hbar}{m
  \omega_z}}$ by $\lambda_0=2\sqrt{2\pi} \frac{a_s}{l_z}$. For $Rb^{87}$,
  $a_s=90.4a_0$. Using this scattering length
  and $\lambda_0(N-1)=17.1$ gives $\omega_z=37.9 $kHz for $N=100$. So the
  transverse aspect ratio is $\frac{\omega_0}{\omega_z}=3.43 \times 10^{-2}$
  and the maximum energy per particle, $E/\hbar N=1.41 $kHz is merely $3.72\times
  10^{-2} \omega_z$. Thus, we can truly treat the dynamics as
  two-dimensional.}

\bibitem{Tsubota2002}
\bibinfo{author}{Tsubota, M.}, \bibinfo{author}{Kasamatsu, K.} \&
  \bibinfo{author}{Ueda, M.}
\newblock \bibinfo{title}{Vortex lattice formation in a rotating Bose-Einstein
  condensate}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{\textbf{65}}}, \bibinfo{pages}{023603}
  (\bibinfo{year}{2002}).

\bibitem{Hofferberth2007}
\bibinfo{author}{Hofferberth, S.}, \bibinfo{author}{Lesanovsky, I.},
  \bibinfo{author}{Fischer, B.}, \bibinfo{author}{Schumm, T.} \&
  \bibinfo{author}{Schmiedmayer, J.}
\newblock \bibinfo{title}{Non-equilibrium coherence dynamics in one-dimensional
  Bose gases}.
\newblock \textit{\bibinfo{journal}{Nature}}
  \textbf{\bibinfo{volume}{\textbf{449}}}, \bibinfo{pages}{324}
  (\bibinfo{year}{2007}).

\bibitem{Anderson2008}
\bibinfo{author}{Weiler, C.~N.} \textit{et~al.}
\newblock \bibinfo{title}{Spontaneous vortices in the formation of
  Bose-Einstein condensates}.
\newblock \textit{\bibinfo{journal}{Nature}}
  \textbf{\bibinfo{volume}{\textbf{455}}}, \bibinfo{pages}{948}
  (\bibinfo{year}{2008}).

\bibitem{Fetter2009}
\bibinfo{author}{Fetter, A.~L.}
\newblock \bibinfo{title}{Rotating trapped Bose-Einstein condensates}.
\newblock \textit{\bibinfo{journal}{Rev. Mod. Phys.}}
  \textbf{\bibinfo{volume}{\textbf{81}}}, \bibinfo{pages}{647}
  (\bibinfo{year}{2009}).

\bibitem{Hall2010}
\bibinfo{author}{Freilich, D.~V.}, \bibinfo{author}{Bianchi, D.~M.},
  \bibinfo{author}{Kaufman, A.~M.}, \bibinfo{author}{Langin, T.~K.} \&
  \bibinfo{author}{Hall, D.~S.}
\newblock \bibinfo{title}{Real-time dynamics of single vortex lines and vortex
  dipoles in a Bose-Einstein condensate}.
\newblock \textit{\bibinfo{journal}{Science}}
  \textbf{\bibinfo{volume}{\textbf{329}}}, \bibinfo{pages}{1182}
  (\bibinfo{year}{2010}).


\bibitem{rosti2013}
\bibinfo{author}{Doganov, R.~A.}, \bibinfo{author}{Klaiman, S.},
  \bibinfo{author}{Alon, O.~E.}, \bibinfo{author}{Streltsov, A.~I.} \&
  \bibinfo{author}{Cederbaum, L.~S.}
\newblock \bibinfo{title}{Two trapped particles interacting by a finite-range
  two-body potential in two spatial dimensions}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{\textbf{87}}}, \bibinfo{pages}{033631}
  (\bibinfo{year}{2013}).

\bibitem{Kramer1981}
\bibinfo{author}{Kramer, P.,} \& \bibinfo{author}{Saraceno, M.}
\newblock \textit{\bibinfo{title}{Geometry of the Time-Dependent Variational
  Principle in Quantum Mechanics}} (\bibinfo{publisher}{Springer},
  \bibinfo{year}{1981}).

\bibitem{Lode2012}
\bibinfo{author}{Lode, A.~U.~J.}, \bibinfo{author}{Sakmann, K.},
  \bibinfo{author}{Alon, O.~E.}, \bibinfo{author}{Cederbaum, L.~S.} \&
  \bibinfo{author}{Streltsov, A.~I.}
\newblock \bibinfo{title}{Numerically exact quantum dynamics of bosons with
  time-dependent interactions of harmonic type}.
\newblock \textit{\bibinfo{journal}{Phys. Rev. A}}
  \textbf{\bibinfo{volume}{\textbf{86858}}}, \bibinfo{pages}{063606}
  (\bibinfo{year}{2012}).

\bibitem{AxelThesis}
\bibinfo{author}{Lode, A.~U.~J.}
\newblock \textit{\bibinfo{title}{Tunneling Dynamics in Open Ultracold Bosonic
  Systems}} (\bibinfo{publisher}{Springer}, \bibinfo{year}{2014}).
  
  
\bibitem{MCTDH-rev}
\bibinfo{author}{Beck, M.~H.}, \bibinfo{author}{Jackle, A.},
  \bibinfo{author}{Worth, G.~A.} \& \bibinfo{author}{Meyer, H.-D.}
\newblock \bibinfo{title}{The multiconfiguration time-dependent Hartree (MCTDH)
  method: A highly efficient algorithm for propagating wavepackets}.
\newblock \textit{\bibinfo{journal}{Phys. Rep.}}
  \textbf{\bibinfo{volume}{\textbf{34}}}, \bibinfo{pages}{1}
  (\bibinfo{year}{2000}).


\end{thebibliography}
