Exploring vortex formation in rotating Bose-Einstein condensates beyond mean-field regime

The formation of quantized vortices is a natural consequence of angular momentum in rotating fluids. These quantized vortices can be produced and observed in Bose-Einstein condensate (BEC) setups [1, 2, 3, 4]. Various techniques are used to generate quantized vortices in BECs, such as laser stirring [6], rapid rotation [5, 7], direct phase-defect imprinting [8].

Vortex signatures in BEC manifest themselves in density nodes and phase discontinuity, observed in the density or momentum distribution. A variety of fascinating phenomena has been theoretically predicted and experimentally observed, such as vortex bending [9] Tkachenko oscillations, quantum Hall states, vortex lattice formation [10], vortex amalgamation[11] and melting [12, 2].

A rotating BEC is generally analyzed from a mean-field perspective by solving the Gross–Pitaevskii equation. The mean-field framework, which assumes that all the particles are coherent and condensed in a single state, can aptly predict a variety of phenomena, including the formation of a vortex lattice [13, 14].

The mean-field theory necessitates a system with a huge number of particles having negligible interactions and becomes inexact when the quantum correlations become significant. In particular, both for extremely strong interactions and rapid rotations, fragmentation in BECs occurs, i.e., multiple orbitals become macroscopically occupied [15]. Hence, the fundamental assumptions of the mean-field method are no longer fulfilled, and the corresponding solutions become inaccurate. Therefore, it is necessary to employ a general quantum many-body framework to account for fragmentation and quantum correlations appropriately in such situations.

The multi-configuration time-dependent Hartree method for bosons (MCTDHB)  [16] is a suitable general quantum many-body method capable of addressing the strong interaction regime [17, 18]. The MCTDHB method has been applied to compute the static and dynamics properties in a variety of strongly interacting and fragmented system [19, 20, 21, 22, 23, 24, 25, 26, 27, 28].

This work explores the vortex formation in a strongly interacting, rapidly rotating BEC from a general ab initio quantum many-body perspective by numerically solving the Schrödinger equation using the MCTDHB method. Primarily through the analysis of the one-body density, we identify the characteristic vortex structure signatures that emerge as the interaction coupling and angular velocity are varied. Our work identifies the distinct impact of angular momentum and interaction energy on vortex formation and its structure.

The paper is structured as follows. In Section II we discuss the model Hamiltonian and the setup of our system.

In Section III, we describe our numerical method. In Section IV, we display and discuss the results. Section V concludes the paper.

We consider a two-dimensional rotating BEC comprising $N$ atoms, trapped in a quasi-two-dimensional circular disk potential. The many-body Hamiltonian in the co-rotating frame is given by

	$\displaystyle H={}-\sum_{i=1}^{N}\frac{\hbar^{2}}{2M}\nabla_{i}^{2}+\sum_{i=1}^{N}V_{pot}(\bm{r}_{i})+\sum_{i=1}^{N}V_{rot}(x_{i},y_{i})$		(1)
	$\displaystyle+\sum_{i<j}V_{int}(\bm{r}_{i}-\bm{r}_{j}).$

Here, $N$ is the number of atoms each of mass $M$. The first term represents the kinetic energy of the setup. The trapping potential $V_{pot}$ is a hard-wall circular disk potential in the $X-Y$ plane modeled as

$$V_{pot}(r)=\left\{\begin{array}[]{ll}0&r\leq a\\ 1000&r>a\end{array}\right.,$$

(2)

In the ultracold temperature regime, the energies of the atoms are low and consequently only the s-wave scattering is relevant. Therefore, the exact (Van-der wall-like) potential between the atoms can be modeled as a contact delta-like interaction between two atoms at $\bm{r}_{i}$ and $\bm{r}_{j}$ as $V_{int}(\bm{r}_{i}-\bm{r}_{j})=g\delta(\bm{r}_{i}-\bm{r}_{j})$. Here, $g$ is the interaction coupling constant and regulates the strength of the interaction. However, owing to the computational convergence problems pertaining to the usage of actual delta potentials in two-dimension [29], we model the delta function as a narrow normalized Gaussian function as

$$\delta(\bm{r=\bm{r}_{i}-\bm{r}_{j}})=\frac{\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)}{2\pi\sigma^{2}},$$

(3)

The width of the Gaussian is taken as $\sigma=0.25$, chosen to correspond with the physics of the true delta potential. The BEC rotates in the $z$ direction with an angular velocity $\Omega$. The corresponding Hamiltonian term is given as $V_{rot}=\Omega\bm{L_{z}}$, where $\bm{L_{z}}=-i\hbar(x\partial_{y}-y\partial_{x})$. The ground-state of the co-rotating Hamiltonian 1 provides the steady-state solution of the system. For our numerical computation, we rescale the Hamiltonian 1 by $\frac{\hbar^{2}}{2ML^{2}}$ where $L$ is a typical length-scale of the system. This makes all quantities in the rescaled Hamiltonian dimensionless and defines the units of computed quantities.

We use the numerically exact multi-configuration time-dependent Hartree method for bosons (MCTDHB) implemented in the MCTDH-X package [18] to solve the many-body Schrödinger equation. In this method, the many-body wavefunction is expanded in terms permanents

$$|\Psi\rangle=\sum_{\vec{n}}C_{\vec{n}}(t)|\vec{n};t\rangle.$$

(4)

The permanents are constructed as a linear combination of time-dependant, symmetrized, variationally optimal bosonic basis states as:

$$|\vec{n};t\rangle=\frac{1}{\sqrt{\prod_{\alpha=1}^{M}n_{\alpha}!}}\prod_{\alpha=1}^{M}\left(b_{\alpha}^{\dagger}(t)\right)^{n_{\alpha}}|\textrm{vac}\rangle,$$

(5)

Here, $N$ bosons are distributed in $M$ single-particle states. Subsequently, a time-dependent variational method is used [30] to generate two sets of equations of motion  [16] which are solved to obtain the many-body solutions. Fundamentally a time-dependent method, the ground-state is obtained using propagation in imaginary time.

We focus on the strong interaction and rapid rotation regime wherein the particles are no longer entirely condensed, and hence the mean-field approximation is no longer accurate. We consider $N=100$ in the disk potential. To explore the formation of quantized vortices, we analyze the one-body density as a function of the interaction strength $g$ and the angular velocity $\Omega$.

We aim to examine the impact of angular velocity and interaction strength on vortex formation. To that end, we first study the non-rotating $\Omega=0$ and non-interaction case $g=0$ as a reference. Figure 1 (top panel) depicts the one-body densities for no rotation $\Omega=0$. Absent both interaction and rotation, the density shows a layered profile – the density is largest at the center and reduces along the radius being lowest at the rim. This can be attributed to the dominating influence of the kinetic energy and the hard-wall potential forcing the particles towards the center.

As interaction increases ($g=0.1$), the repulsive interaction results in the particles being pushed out, countering the kinetic energy influence. As the interaction becomes significant $g=0.5--1$, the interaction energy dominates, resulting in a uniform disk density.

Figure 1 (bottom panel) depicts the one-body densities for no interaction $g=0$. For significant angular velocity $\Omega\geq 0.5$, an extensive region of very negligible density is formed at the center. This can be attributed to the centrifugal effect of the rotation pushing the particles radially outward. The size of this region increases with increasing angular velocity and, at $\Omega=0.9$, covers most of the disk.

We now introduce both rotation and interaction in our system. First, we investigate the physics with low angular velocity ($\Omega=0.2,0.5$).

The introduction of interaction and rotation leads to their interplay, impacting the vortex production.

Fig 2(top panel) depicts the one-body density for $\Omega=0.2$. For interaction strength $g=0.1$,the one-body density is roughly similar to the non-interacting case – a layered profile. For $g=0.5$, the rotation results in a vortex formation exactly at the center of the disk. This is revealed as a density node or a ‘hole’ in the one-body density. The density assumes a uniform pattern in the other regions. The vortex disappears for stronger interaction $g>1$, and we acquire a roughly uniform density profile.

Fig 2(bottom panel) depicts the one-body density for moderate angular velocity $\Omega=0.5$. The higher angular velocity enhances the vortex production in the BEC.

Here, for $g=0.1$, we observe the formation of a large density hole at the center covering about one-third the radius of the disk. The formation of this density hole can be attributed to the dual effect of rotation and interaction pushing the particles to the rim.

As interaction increases to $g=0.5$, we observe the formation of vortex structures. The large central hole transforms into a low-density triangular region with a vortex at the vertices. This triangular triple-vortex pattern transforms into a double-vortex structure for $g=1$. For very strong interaction $g=10$, we observe a solitary vortex exactly at the center. This single vortex persists for stronger interactions. We observe that the higher angular velocity enhances the vortex production, especially for the intermediate interaction regime. However, very strong interaction reduces the number of vortices produced.

We now explore the larger angular velocity regime. Fig 3(top panel) depicts the one-body density for $\Omega=0.7$. Because of the larger angular velocity, the interplay between the angular momentum and interaction occurs at a much lower interaction coupling. For small interactions, we have a similar central region hole as before (not shown here).

Increasing the interaction coupled with the faster rotation leads to an enhancement of vortex production and the richness of various structures formed. At $g=0.7$, the central circular hole transforms into a pentagonal shape of low density. The vortexes arrange themselves into the vertices of this pentagonal geometry. The pentagonal vortex geometry transforms into a diamond-shaped configuration for $g=1$ with vortexes at the four vertices. Further interaction $g=2$ reduces the structure into a triangular geometry comprising 3 vortexes. Finally, at strong interaction $g=8$, we obtain a two-vortex structure that persists for stronger interaction, demonstrating the suppression of vortex formation with very high interactions.

Fig 3(bottom panel) depicts the one-body density for $\Omega=0.9$. Here, the rapid rotation results in a large hole region for $g=1$ – the rapid rotation effectively destroys the previously obtained pentagonal and diamond vortex lattices. At $g=2$ instead, we obtain a square structure comprising 8 vortexes. Interestingly, the eight vortexes do not arrange as an octagon but a square – this specific geometrical configuration minimizes the energy. Further interaction $g=3$ reduces the 8-fold vortex structure into a 4-vortex diamond. The consistent reduction of vortex formation for higher interactions is also seen here. At $g=8$, the structure is reduced into two vortexes, which persist for larger interactions.

These simulations display the contrasting influence of the interaction strength compared to the rotational velocity. While the increasing angular velocity $\Omega$ leads to an increased formation of vortexes having a rich plethora of structures, the influence of interaction strength $g$ is more involved. Whereas increasing from small to intermediate interaction enhances the vortex production, very strong interactions destroy the vortex structures, leading to a reduction of the vortexes.

We explored the vortex formation in a rotating Bose-Einstein condensate through an ab initio numerical simulation of the many-body Schrödinger equation. The vortex structure and formation depend on the interplay between the angular velocity and the interaction strength.

Only a single vortex is generated for slow rotation and intermediate interaction strength. The vortex production is enhanced for faster angular velocity, and various vortex structures form. Here, the impact of interaction becomes critical – increasing the interaction to intermediate strength assists in the vortex structure formation. However, very strong interaction strength results in suppressing the vortex formation and the melting of the corresponding structures.

For intermediate rotation ($\Omega=0.5$), we observe a transformation of a large central hole into a triangular vortex geometry when the interaction increases. For faster rotations ($\Omega=0.7$), richer structures emerge for intermediate interaction strength ($g=0.7-3$), that includes a pentagon, diamond, triangle geometry finally melting into a double vortex. For high-speed rotations ($\Omega=0.9$), we observe an octet square structure transforming into a triangle and a double vortex as the interaction increases. The simulation demonstrates the vortex formation and destruction tendency of interactions. The number and structures of vortices are most significant for intermediate interactions and decrease as the interaction increases.

In general, we observe that while both interaction strength and angular velocity facilitate the formation of the quantized vortex of various structures, a powerful interaction leads to the breakup of the vortex formation. The higher rotational rate favors the formation of a greater number of vortices having various structures, which is especially evident in the intermediate interaction regime.

We thank A.U.J. Lode and the other contributors for the development and support of the MCTDHX package.