Vortex Patterns of atomic Bose-Einstein condensates in a density-dependent gauge potential

We theoretically examine the vortex states of a gas of trapped quasi-two-dimensional ultracold bosons subject to a density-dependent gauge potential, realizing an effective nonlinear rotation of the atomic condensate. The nonlinear rotation has a two-fold effect; as well as distorting the shape of the condensate it also leads to an inhomogeneous vorticity resulting in novel morphological and topological states, including ring vortex arrangements that do not follow the standard Abrikosov result. The dynamics of trapped vortices are also explored, which differs from the case of linear rotation due to the absence of a global laboratory reference frame.


Introduction.
Quantum vortices represent the fundamental excitations of superfluid systems, appearing in response to the rotation of atomic condensates. Unlike their classical counterparts a quantum vortex's rotational properties are more restricted, as their velocity field is quantized. The first generation of experiments used a laser to induce rotation by stirring the atomic cloud [1][2][3]. Vortices and in particular their interaction dynamics play a central role in condensed matter systems. A single vortex constitutes a core region where the density vanishes around which there is a circulation of the quantum fluid, consequentially there is a defect in the phase of the order parameter. Dimensionality plays a central role in the physics of vortices -for two-dimensional quantum fluids experiments have demonstrated different topological phase transitions, prominent examples being the quantum spin Hall effect and the Berezinskii Kosterlitz Thouless transition, both of which have been realized with ultracold atoms [4,5]. In the three-dimensional context, more elaborate topological configurations are available, such as knots [6][7][8], skyrmions [9,10] and the related problem of engineering analogies of the magnetic monopole [11,12].
Ultracold atomic gases constitute exemplar physical systems for examining quantum mechanical phenomena, since experiments afford high controllability. Complementary to this, quantum gases represent dilute systems, which facilitates their accurate theoretical modelling via the celebrated Gross-Pitaevskii model and its numerous extensions. Vortices central role in superfluidity continues to attract theoretical and experimental interest in the macroscopic dynamics of these excitations. Early work focussed on studying the fundamental properties of the rotating system [13][14][15], while more recent work has focussed on understanding the effect of anisotropic trapping [16], vortex lattice [17,18], and chaotic [19] dynamics. Focus has also been on the structure and dynamics of vortices in condensates at finite temperature including non-equilibrium effects [20,21], multi-component systems which have been shown to possess a rich vortex physics [23][24][25][26][27][28] and the on-going quest to understand quantum turbulence [29,30]. For a detailed description of the basic properties of rotating gases, the interested reader is directed to the review of Fetter [31].
The last decade has seen the experimental realization of artificial electromagnetism in quantum gases [32,33]. This provides an important new tool for accessing some of the paradigmatic effects associated with condensed matter physics, including simulating orbital magnetism [34][35][36][37]. Ultracold atomic gases provide a platform for implementing not only these traditional manifestations of magnetism; but also for realizing forms of magnetism that are not naturally occurring, such as spin-orbit [38,39] and spin-angular momentum coupling [40,41] with bosons, synthetic dimensions [42], 'knotted' gauge theories [43] and density-dependent magnetism [44]. The realization of spin-orbit coupled bosons in-particular yields a new route to investigate the interplay of synthetic electromagnetism and rotation [45][46][47]. Induced electromagnetism in atomic gases typically produces a static gauge potential, here there is no feedback between the light and the matter-wave. To address this, several proposals have been put forward [48][49][50] to instead simulate dynamical gauge potentials. Very recently the first generation of experiments to realize density-dependent magnetism have appeared, for bosons [51] and fermions [52] trapped in two-dimensional optical lattices, as well in a system of Rydberg atoms [53].
Theoretical model. We consider a system comprising N two-level atoms coupled via a coherent light-matter interaction, forming a Bose-Einstein condensate. The Hamiltonian describing our setup can be written aŝ where the light-matter interaction is defined cos θ e −iφ sin θ e iφ sin θ − cos θ .
The other quantities appearing in Eq. (1) are the external trapping potential V ext (r), and the mean-field interac- , where ∆ j = g jj n j +g jk n k and g jk = 4π 2 a jk /m defines the scattering parameter between atoms in the internal states j and k, while n j = |ψ j (r)| 2 defines the population density of state j.
Meanwhile the harmonic trapping potential is given by V ext (r) = m(ω 2 x x 2 + ω 2 y y 2 + ω 2 z z 2 )/2, where ω j defines the strength of the trapping potential in each coordinate direction. To build an interacting gauge theory, we construct interacting dressed states using perturbation theory which is valid when the mean-field interactions are weak compared to the strength of the light-matter coupling. Denoting the (unperturbed) eigenstates ofÛ MF (Eq. (2)) as |± , the interacting dressed state basis can be written as The dressed mean-field detuning is ∆ d = sin θ 2 cos θ 2 (∆ 1 − ∆ 2 )/2. To build the interacting gauge theory, we construct a state vector from the two perturbed dressed states defined by Eq. (3). The qualitative features do not depend on this choice, here we project the atoms motion into the + state. Then, the effective Hamiltonian becomeŝ where the two geometric phases that arise due to the adiabatic motion of the atoms are given by A + = i ψ + |∇ψ + and W + = 2 | ψ + |∇ψ − | 2 /2m respectively, while the dressed mean-field interactions are defined by ∆ + = (∆ 1 cos 2 θ 2 + ∆ 2 sin 2 θ 2 )/2. These terms can be shown to be defined as To obtain a mean-field equation of motion for the atoms, we extremize the energy functional E = ψ + |Ĥ + |ψ + , leading to the following generalized Gross-Pitaevskii equation [54][55][56] Where the coupling to the gauge field is defined as a 1 = (∇φ∆ d sin θ/n + Ω r ). Then, the density-dependent dressed basis, Eq. (3) gives rise to additional terms in the generalized Gross-Pitaevskii equation for atoms in the ψ + state, including the current nonlinearity j which appears as The current nonlinearity Eq. (7) gives rise to novel topological states; in-particular proposals to simulate exotic spacetime geometries [57] in the three-dimensional context. In the one-dimensional limit the theory violates Kohn's theorem [58,59], and exact chiral soliton solutions [60] can also be constructed in this limit, which have recently been shown to constitute quantum time crystals [61]. We identify two small parameters: θ, the ratio of the Rabi frequency to the detuning, and ε = nγ 12 / ∆ that encompasses the collisional and coherent interactions. After expanding Eqs. (5a)-(5b) to linear order in ε and θ we obtain simplified expressions for A + and W + .
To build the interacting gauge theory, we choose a laser beam carrying units of angular momentum such that the phase φ = ϕ, where ϕ is the polar angle [62]. The spatial profile meanwhile satisfies θ = θ 0 r, where r is the radial distance. Then, using Eqs. (6) and (7) along with the simplified expressions for A + and W + , the equation  of motion for the condensate becomes here the angular momentum operator is defined asL z = −i ∂/∂ϕ. Then we define the density-dependent rotation frequency as Ω(r, t) = Ω + Cn(r) where the strength of the nonlinear term is C = θ 2 0 (g 11 − g 12 )/(2mΩ r ), while the effective scattering parameter is given by g eff = g 11 + θ 2 0 (g 11 − g 12 )(3 − 2/ 2 )/(mΩ r ). As we are interested in the quasi-two-dimensional situation, we assume a pancake geometry for the atomic cloud (ω z ω x,y ), which allows us to factorize the atomic wave function as ψ(r, t) = ψ(x, y, t)e −z 2 /2σ 2 z /( 4 πσ 2 z ), which has the effect of rescaling the two nonlinear terms g eff and C by 1/ √ 2πσ z . Ground states and vorticity. To gain an understanding of the basic physics associated with a trapped condensate under nonlinear rotation, we begin with the energy functional associated with Eq. (8) which can be written in the hydrodynamic prescription using the Madelung transformation ψ = √ ne iφ where n ≡ n(ρ, t) is the quasi-two-dimensional density while φ ≡ φ(ρ, t) is the corresponding phase, then after dropping the quantum pressure term (valid for ma x g eff / √ 2πσ z 2 1) we obtain here v = ( ∇φ − A)/m defines the kinetic velocity with A = mΩ × ρ, Ω =ê z (Ω + Cn/2), and V 2D (x, y) = m(ω 2 x x 2 + ω 2 y y 2 )/2. The energy defined by Eq. (9) can be minimized to obtain the superfluid velocity v sf = Ω × ρ which couples to both the linear rotation Ω and the density n(ρ) of the gas. Using Eq. (9) and the expressions for the superfluid velocity, we can obtain a generalized Thomas-Fermi distribution in the rotating frame as where ρ 2 = x 2 + y 2 and µ is the chemical potential in the rotating frame, and the radius of the cloud is defined as R 2 x,y = 2µ /m(ω 2 x,y − Ω 2 ), while n(ρ) is subject to the normalization d 2 ρ n(ρ) = N . Figure 1 (a) and (b) compute the chemical potential µ and density n(ρ). Panel (a) shows that the solutions for µ exist on finite regions, between a minimum and a maximum C. The maximum value of C in each case corresponds to the point where the gas locally exceeds the trapping frequency, while the minimum C corresponds instead to the point when the approximations leading to Eq. (10) breakdown. The solid blue curve in panel (b) shows an example for large positive C, where the density profile exhibits a large plateau region, reminiscent of a quantum droplet [63]. For large negative C, the tails of the distribution appear to decay as they approach the edge of the cloud, with a small central plateau region. The radial density n(ρ) is plotted as a function of C in panel (c), for Ω/ω x = 0.6. The final panel (d) of Fig. 1 computes the boundaries that the solutions occupy in the (Ω, C) parameter space. Accompanying this is the critical rotation frequency Ω c at which it is energetically favorable for a vortex to enter the cloud [56,65], calculated from Ω c = ( E v − C nL z /2)/ L z , (green dashed line) using the non-rotating (Ω(r, t) = 0) solution. We define the average vorticity as ω v = ∇ × v, which can be calculated from v = Ω(r, t)ê z × ρ giving and the synthetic magnetic field is related to Eq. (11) since B(ρ) = mω v (ρ). Numerical simulations. Figure 2 presents calculations of the ground states of Eq. (8) as a function of the strength of the nonlinear rotation C. Throughout we take the van der Waals strength as ma x g eff N/( √ 2πσ z 2 ) = 300, while the linear rotation strength Ω/ω x = 0.6, and trap anisotropy is ω y /ω x = 1.01. Panel (a) shows the number of vortices as the strength of C is changed. The labels correspond to the panels (c) to (j) showing the atomic density n(ρ) and phase φ(ρ) of the ground state in each case. For large positive C, an unusual cylindrical density is observed (panel f), coroborating the analytical prediction, Fig. 1(b) which in the presence of vortices resembles a Renkon lotus root. For negative C, the triangular lattice is no-longer observed. Instead the vortices arrange into rings with increasing radius as C is decreased (for larger Ω, concentric vortex ring arrangements are observed whenC 0). The localization of the vortices can be interpreted from Eq. (11) due to the vorticity being maximum in the center of the cloud whenC > 0, while forC < 0 it is maximal at the edges of the cloud, leading to ring arrangements. Panel (b) shows the energy E (gray triangles) and angular momentum L z (solid blue circles) along with a linear fit for the angular momentum (see also Ref. [64] concerning vortex patterns in a related anyonic model).
Next we explore the effect of a general elliptical trapping potential (ω y ω x ) on the ground state in fig. 3. The parameters are chosen such that ω y /ω x = 1.5 and Ω/ω x = 0.8. Panels (c) and (d) is broken, instead the remaining vortices are located at either edges of the cloud. The nonlinear rotation gives rise to unusual vortex dynamics. For the case of linear rotation, it is always possible to define and move between a lab (LB) and co-moving (CM) frame via the transformation |ψ CM =Û (t)|ψ LB , yielding the HamiltonianĤ CM =Û (t)Ĥ LBÛ (t) † + i ∂ tÛ (t)Û (t) † whereÛ (t) = exp(itΩ · L/ ). Since the density-dependent gauge theory facilitates a timedependent rotation frequency, this simple transformation does not hold for C = 0. This in turn means that this nonlinear system does not possess a global laboratory frame of reference. A simple illustration of the consequence of this is presented in Fig. 4. Here the dynamics of two vortices are prepared such that initially they are in a nonsymmetric arrangement such that (x v1 , y v1 ) = (0, a x ) and (x v2 , y v2 ) = (0, −2a x ) (panels (a,d)). Then one can see that the time evolution (panels (b,c,e,f)) causes the dynamics with Ω = 0 to differ depending on the value of C. Panel (g) presents the space-time trajectories of both situations.

Conclusions.
We have calculated the ground state configurations of a trapped atomic Bose-Einstein condensate in a density-dependent gauge theory. Depending on the sign and magnitude of the corresponding nonlinear rotation, different vortex configurations were obtained. When the strength of the nonlinear rotation is large and positive, almost flat-topped vortex distri-butions were found, whereas for large negative nonlinear rotation strengths the vortices instead arrange into ring structures. For condensates experiencing a highly anisotropic confining potential, the vortices found at large negative rotation frequencies separate into two regions at opposing edges of the trap, breaking the ring configuration. We also examined the dynamics of the vortices in this system, which are unusual due to the lack of a global laboratory reference frame. It would be interesting in the future to study nonlinear rotation in the three-dimensional situation as well as turbulent vortex dynamics.
Acknowledgements. MJE would like to thank Sven Bjarke Gudnason for assistance with the HPC aspects of this work and Salvatore Butera for comments on the manuscript. Numerical calculations were carried out using the TSC-computer of the "Topological Science" project at Keio University. This work was supported by the Ministry of Education, Culture, Sports, Science (MEXT)-Supported Program for the Strategic Research Foundation at Private Universities "Topological Science (Grant No. S1511006). This work of M.N. is also supported in part by JSPS KAKENHI Grant Numbers 16H03984, and 18H01217 and by a Grant-in-Aid for Scientific Research on Innovative Areas "Topological Materials Science (KAKENHI Grant No. 15H05855) from MEXT of Japan.