Dynamics of Ultracold Bosons in Artificial Gauge Fields—Angular Momentum, Fragmentation, and the Variance of Entropy

We consider the dynamics of two-dimensional interacting ultracold bosons triggered by suddenly switching on an artificial gauge field. The system is initialized in the ground state of a harmonic trapping potential. As a function of the strength of the applied artificial gauge field, we analyze the emergent dynamics by monitoring the angular momentum, the fragmentation as well as the entropy and variance of the entropy of absorption or single-shot images. We solve the underlying time-dependent many-boson Schrödinger equation using the multiconfigurational time-dependent Hartree method for indistinguishable particles (MCTDH-X). We find that the artificial gauge field implants angular momentum in the system. Fragmentation—multiple macroscopic eigenvalues of the reduced one-body density matrix—emerges in sync with the dynamics of angular momentum: the bosons in the many-body state develop non-trivial correlations. Fragmentation and angular momentum are experimentally difficult to assess; here, we demonstrate that they can be probed by statistically analyzing the variance of the image entropy of single-shot images that are the standard projective measurement of the state of ultracold atomic systems.


Introduction
Since the first realization of Bose-Einstein condensates in 1995 [1][2][3], ultracold atoms have become a standard probe for analog quantum simulations-due to the tunability and flexibility of these quantum states of matter, they can be manipulated to behave like other systems, for instance, condensed matter systems which are not as flexible or easy to observe. Popular examples include the realization of the quantum simulation of the superfluid-to-Mott-insulator transition [4,5], quantized conductance [6,7], the Dicke model [8,9], and magnetism realized via artificial gauge fields for ultracold atoms [10].
Such artificial gauge fields can make the neutral ultracold atoms behave as if they were charged particles experiencing a magnetic field and were investigated experimentally and theoretically with an external lattice potential [11][12][13] or without one [14][15][16].
Our paper is structured as follows-in Section 2 we introduce the Hamiltonian and the MCTDH-X method we use, in Section 3 we discuss the observables that we are using in Section 4 to investigate the dynamics of ultracold atoms in an AMF; Section 5 summarizes our conclusions and provides an outlook.

Hamiltonian and Methods
We consider a system of bosonic particles with two-body interactions in two spatial dimensions. The state of the bosons is initialized in the ground state of a parabolic trap without an AMF present. Subsequently, the system is quenched by turning on suddenly an artificial gauge field corresponding to a homogeneous AMF perpendicular to the plane in which the bosons are trapped.
For the sake of clarity of presentation, we will omit the dependence of quantities on time t throughout this work, where it is obvious.

Setup
To setup the time-dependent many-body Schrödinger equation (TDSE), we use the Hamiltonian Here, we work in atomic units (h = m = 1), the potential V(x) [with x = (x, y)] is chosen to be harmonic, V(x) = 1 2 x 2 , and we consider contact interactions W(x, x ) = g 0 δ(x − x ). Formally, one cannot use a contact interaction in two spatial dimensions with a complete basis set since the outcome would be that of the noninteracting bosons [75,76]. In simple terms, this is because the integral measure of the support of the Dirac-δ is zero for two and more spatial dimensions. For a proof, the interested reader is deferred to Ref. [76] and for an example for finding non-zero ranged Gaussian interaction potentials with similar physical behavior for ultracold bosons, see Ref. [75].
In the present work, we employ a finite truncation of the many-body basis (chiefly M = 4 orbitals) and aim to demonstrate that beyond-mean-field phenomena do emerge. The kinetic energy is augmented with an artificial gauge field A(x; t): For simplicity, we consider the case of unit charge g = 1 and a homogeneous magnetic field B in z-direction of strength B(x; t): Here,ê z denotes the unit vector in z-direction. In the following, we work in Landau gauge, and consider a quench scenario in the following, that is, Here, Θ(t) denotes the Heaviside step function, that is, the magnetic field is suddenly turned on at t > 0 after the system has been initialized, see Figure 1 for a sketch. is switched on suddenly. This quench triggers many-body dynamics of the state (label "t > 0").
In our investigation, we analyze the many-body dynamics by monitoring observables as a function of the effective magnetic field strength B at t > 0 after the quench.

Quantities of Interest
Here, we define the observables that we use to quantify the dynamics of N-boson systems: the one-body density, the eigenvalues of the reduced one-body density matrix (1BDM), the angular momentum, and the image entropy and its variance evaluated from simulated single-shot images.

Density, One-body density matrix, and Natural Occupations:
The 1BDM is a hermitian matrix defined as Here, we used the matrix elements ρ kq = Ψ|b † kb q |Ψ to represent the 1BDM using the orbitals corresponding to the creation and annihilation operatorsb † k andb q , respectively. The diagonal of the one-body density matrix is referred to as the density ρ(x): The density ρ(x) is a real quantity and has no phase, because it is the diagonal of a hermitian matrix, ρ (1) (x, x ). The eigenvalues of the 1BDM, Equation (9), can be obtained via a diagonalization that corresponds to a unitary transformation of the orbitals φ j (x; t) to the so-called natural orbitals φ The eigenvalues λ j are normalized, ∑ M j=1 λ j = 1 and, without loss of generality, we consider them to be sorted by size, λ 1 ≥ λ 2 ≥ . . . , throughout this work. The eigenvalues λ j (natural occupations) determine the degree of condensation and fragmentation of the system. Bosons with a 1BDM with only a single contributing eigenvalue λ 1 are condensed [82] and bosons with a 1BDM with multiple macroscopic eigenvalues contributing λ 1 ∼ O(N); λ 2 ∼ O(N); . . . are fragmented [60,61].
The eigenvalues of the 1BDM can be used as a precursor of correlations in the state |Ψ . This can be seen using, for instance, the Glauber first-order correlation function [83], If the system is in a condensed state, the 1BDM has only a single eigenvalue and is therefore a product of a single (complex-valued) orbital, ρ (1) -for particles in a condensed state correlations are absent. It is instrumental to note here, that this condensed, single-orbital case with absent correlations is presupposed in mean-field approaches like the time-dependent Gross-Pitaevskii theory [17,18]. Similarly, when ρ (1) (x, x ) is a sum of two or more orbitals as in Equation (9) then it can no longer be represented simply using its diagonal, the density ρ (1) . Furthermore, in this fragmented case, the denominator of Equation (12) is a product of two weighted sums of orbitals, ρ(x), and ρ(x ), respectively. This product involves non-trivial cross terms and the correlation function attains a value |g (1) (x, x )| ≤ 1 for all (x, x )-for particles in a fragmented state, correlations are present and the single-orbital picture of mean-field approaches like the time-dependent Gross-Pitaevskii theory [17,18] cannot be applied.

Angular momentum:
The angular momentum operator inê z -direction for a two-dimensional system is defined asL Bosonic quantum systems with angular momentum are rich in physics: they feature condensed vortices [17,18,84], phantom vortices [69,85], spatially partitioned many-body vortices [68,86], and fragmentation [68,69,[85][86][87][88]. Since phantom vortices are the most pronounced characteristic feature of angular momentum which we find in our study below, we discuss them in the following. Typically, the term "vortex" refers to a topological defect in the density of quantum system connected with a discontinuity in the phase. A phantom vortex is an analog of such a conventional vortex, but for a natural orbital. A phantom vortex thus represents a topological defect connected with a discontinuity in the phase in one of the field modes of a many-particle states that corresponds to a natural orbital. Phantom vortices were shown to emerge as topological defects in the correlation function in Reference [69]. Moreover, in the common detection scheme for cold atoms, single-shot images (see below) they show as topological defects whose position fluctuates from image to image, see Reference [85].

Single shots, image entropy and its variance:
To assess the observability of the emergent physics in experimental setups with ultracold atoms, we simulate the detection of our numerical model wavefunctions in absorption or single-shot images [9,85,89]. A set of N s single shots, is nothing but N s random samples that are N-variate and distributed according to the N-particle probability given by |Ψ| 2 , To generate images from these single shots, we convolute them with a point spread function. Typical choices include Gaussian (see [9,44,49,85]) or even quantum point spread functions [90]. Here, for simplicity, we consider the idealized case of a δ-shaped point spread function to obtain our single-shot images: We will consider the image entropy ζ of single-shot images of the state |Ψ : In the limit of large N s , the image entropy ζ is equivalent to the density-entropy studied, for instance, in Reference [91]. Fundamentally, the image entropy is a measure for the information content in the particle distributions detected in single-shot images. While Ref. [91] found the entropy to be connected to the presence of correlations in the state, we found it not to be a conclusive pointer in our present work. The variances of observables, however, serve as a precursor of quantum fluctuations and correlations in many-body systems [39][40][41][42][43]; we are thus motivated to also analyze the variance σ ζ of the image entropy ζ:

Results
We now carve out the connection between artificial gauge fields and many-body correlations. For this purpose we focus on the dynamics of a model system of N = 100 two-dimensional ultracold bosonic atoms with an interaction strength of g 0 = 0.05 [cf. Equations (1)-(5) for t ≤ 0]. We now make an example for realizing this interaction g 0 with a real trap configuration with 87 Rb in analogy to the study [69]-the Hamiltonian in Equation (1) is multiplied by¯h mL 2 , where m = 1.44 × 10 −25 kg is the mass of a 87 Rb and L is a length scale that we choose to be L = 0.75 × 10 −6 m. This sets our trapping frequency to ω = (2π)207 Hz and yields a unit of time, mL 2 h of 4.84 ms. The total interval of time we consider in the following t ∈ [0, 200] thus corresponds to 0.97 s. In quasi-twodimensional setups, the interaction parameter g 0 = 2 √ 2π a s l z depends on the transversal oscillator length, l z = h mω z , and the scattering length of 87 Rb, a s = 90.4a 0 ; here, ω z is the transversal trapping frequency and a 0 the Bohr radius. The two-dimensional interaction strength we use, g 0 = 0.05, is obtained for a transversal trapping of ω z ≈ (2π)3.178 kHz. We remark here, that our choice of g 0 corresponds to a weak interaction; the healing length ξ = 1 2g 0 (N−1) ≈ 0.101 is comparable to the oscillator length, that is, 1 in our units. The system is initialized in its ground state and its dynamics (t > 0) are then triggered by suddenly turning on an AMF of strength B [Equation (5)]. In what follows, we aim at an understanding of how the strength of the AMF affects the emergent dynamical behavior.  For a comparatively strong AMF, B = 6.5, in contrast, vortices at the edges of the density (so-called "ghost vortices" [92]) and phantom vortices [69] in the orbitals emerge in To quantify the dynamics of angular momentum triggered by quenches of the AMF a bit better, we plot L z /N for our system as a function of evolution time and as a function of the strength of the AMF in Figure 3.
We find from Figure 3a that a threshold AMF strength of about B 6 is required to generate states with significant angular momentum at long evolution times (here, t = 200). Furthermore, the average angular momentum content increases as the strength of the AMF does. We remark here that the angular momentum features an oscillatory behavior for all AMF strengths that we investigated. This can be understood as a consequence of our quench scenario-the initial state is an eigenstate of the Hamiltonian in Equation (1) for time t < 0 and has vanishing angular momentum L z = 0. At t = 0, the Hamiltonian and its spectrum is changed abruptly due to presence of the AMF, Equation (5). The dynamics of the initial state is thus dependent on these changes in the spectrum of the Hamiltonian by the AMF. At low AMF (B 6) strength, the time-evolution of |L z | quasi-periodically goes back to L z = 0 and is dominated by oscillations at a single frequency [ Figure 3b]; we infer that only very few states with |L z | = 0 are contributing. The dynamics of the state |Ψ of the many-boson system is thus a superposition of a very small number of eigenstates of the Hamiltonian after the quench in this quasi-adiabatic case. The situation changes for B 6, where the time-evolution of |L z | does not return to L z = 0 quasi-periodically. In this case, the time-evolution of |L z | still has its shortest-time oscillations of roughly the same amplitude as for the quasi-adiabatic case B 6 [ Figure 3c]. However, several large-amplitude oscillations with other frequencies contribute. We infer, that the dynamics of the state |Ψ of the many-boson system is thus a superposition of a large number of states of the Hamiltonian after the quench, in this genuinely non-adiabatic evolution for B 6.
We remark here that our results on the oscillatory behavior of the angular momentum in time render it impossible to approach the physics of the many-body state using a co-rotating frame at a certain angular frequency as done in References [62,93].
In References [62,68,69,85,87,88,[93][94][95][96][97], an intricate connection of angular momentum content and the presence of correlations or the fragmentation of many-boson states has been pointed out. This motivates us to analyze the time-evolution of the eigenvalues of the reduced one-body density matrix as a precursor of correlations and the departure of the analyzed state from a mean-field description; we, thus, underpin the limitations of a mean-field description, see Figure 4 for a plot of λ j as a function of time and strength B of the AMF.
Our present findings for a quenched AMF are in line with results obtained with a time-dependent and slow transfer of angular momentum via a rotating asymmetry of the harmonic trapping potential [69,85]: the dynamical departure from a single-eigenvalue 1BDM to a fragmented many-body state with correlations accompanies the dynamical acquisition of a significant amount angular momentum. In particular, the region in time and AMF strength, with a quasi-adiabatic evolution of small |L z | in Figure 3 agrees roughly with the regions with small fragmentation in Figure 4. Conversely, the region in time and AMF strength, with a non-adiabatic evolution of large |L z | in Figure 3 agrees roughly with the regions with substantial fragmentation in Figure 4. The dynamical emergence of multiple significant eigenvalues, i.e., the fragmentation of the many-body state, is quantified here via the time-dependent depletion 1 − λ 1 = ∑ M k=2 λ k and is in sync with the dynamics of angular momentum (cf. Figure 3).
We now turn to the question of the possibility of an experimental detection of the emergent behavior of angular momentum and the eigenvalues of the one-body density matrix. For this purpose we simulated N s = 1000 single-shot images for all the many-body wavefunctions |Ψ(t) for every time in t ≡ k · dt ∈ [0, 200] in steps of dt = 1.0. From this dataset of single-shot images, we computed the image entropy and its variance σ ζ [Equation (18)]; see Figure 5. The variance σ ζ becomes significant for the same values of time and AMF strength where the natural occupations in Figure 4 herald the emergence of fragmentation and correlations. We emphasize here, that fragmentation implies the inapplicability of mean-field descriptions that are restricted to fully condensed and uncorrelated states, as already described in Section 3. By itself it is an interesting result that the variance of image entropy represents an experimentally feasible way to study the emergence of fragmentation in ultracold bosonic systems. Moreover, we observe that the region in time and AMF strength B, where |L z | features quasi-adiabatic (non-adiabatic) dynamics in Figure 3, coincides roughly with the region where fragmentation and image entropy are low (large) in Figures 4  and 5, respectively. We therefore conjecture that an experimental detection of the presence of a non-adiabatic evolution of angular momentum is feasible by measuring the variance of the image entropy.
We remark that we have not shown complementary results on the image entropy ζ [Equation (17)] as a function of AMF strength and time, because it shows only very little overall variation. We noted that the entropy ζ is lowest, where its variance ( Figure 5) is changing the most. The exploration of a possible fundamental connection of the gradient of the image entropy and the variance σ ζ goes beyond the exploratory scope of our present investigation.

Conclusions and Outlook
We analyzed the dynamics of interacting two-dimensional ultracold bosonic particles triggered by a quench of an artificial gauge field. We used the multiconfigurational timedependent Hartree method for indistinguishable particles software (http://ultracold.org) (accessed on π-day, 14 March 2021) to solve the many-body Schrödinger equation from first principles. Our exploratory investigation demonstrates that fragmentation emerges due to the quench if the artificial magnetic field is sufficiently strong. Such an emergence of fragmentation entails the breakdown of conventionally-used mean-field descriptions and, therewith, the occurrence of many-body correlations. We underpin our results by checking their consistency across different quality levels of our MCTDH-X approximation and different numbers of particles.
We have portrayed how correlations show in the expectation value of the angular momentum operator and in the orbitals and their phases as phantom vortices. Using simulations of single-shot images, we demonstrate that the fragmentation and correlations can be detected via the variance of the entropy of the images.
Our work highlights the importance of deploying modern computational and theoretical many-body approaches like the MCTDH-X to systems with artificial gauge fields as well as the necessity to consider not only the wavefunctions of ultracold atoms themselves, but also their detection.
Straightforward continuations of this investigation include the deployment of the developed analysis and computational tools to other many-body systems. As examples of interest, we name here the variance of the image entropy in ultracold dipolar atoms as discussed in Refs. [44,48,49] and an exploration of the competition of long-ranged dipolar interactions and artificial magnetic fields for two-dimensional ultracold atoms. Another direction of physical interest is the dynamics of two-dimensional bosonic Josephson junctions [99] subject to gauge fields and the resulting tunneling of (phantom) vortices or the emergence of quantum turbulence [100] via entropy production [101], which, in turn, as we have shown above, could result from the presence of artificial gauge fields.

Data Availability Statement:
The MCTDH-X software to recompute all the results in this study is available at http://ultracold.org (accessed on π-day, 14 March 2021). Results can be made available upon request; the total amount of data is several TB.
Acknowledgments: Insightful discussions with Saurabh Basu and Ofir E. Alon are gratefully acknowledged.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The Equations of Motion of the Multiconfigurational Time-Dependent Hartree Method for Indistinguishable Particles
We now discuss the equations of motion of the MCTDH-X. For simplicity, we present them using the gauge condition φ i |φ j = 0; for other possible gauges, see for instance [22,102].
i∂ t C n (t) = ∑ n n; t|H| n ; t C n . (A2) Here, we introduced the matrix elements of the reduced one-body and two-body density matrices, respectively. The {ρ} −1 jk notation denotes the jk-th matrix element the inverse of a matrix filled with the entries defined in Equation (A3).
Moreover, the time-dependent local interaction potentials were used In this work, we solve the Equations (A1) and(A2) with the MCTDH-X software [74] hosted at http://ultracold.org (accessed on π-day, 14 March 2021).

Appendix B. Investigation of Convergence
Here, we discuss the convergence of our numerical simulations with respect to the number of variational parameters in the ansatz of the MCTDH-X wavefunction. For this purpose, we consider the same systems as investigated in the main text, but perform the solutions of the Schrödinger equation with MCTDH-X with a different number of orbitals M. Figures A1 and A2 show the result for the first eigenvalue λ 1 of the 1BDM for B = 1 and B = 5, respectively, as a function of time for M ≤ 5.
The time-evolution of λ 1 is not fully converged for M ≤ 5. However, all shown timeevolutions with M > 1 show a significant departure from the mean-field case, i.e., λ 1 < 1 for M > 1. Importantly, the magnitude of the departure from the mean-field case increases drastically as the strength of the AMF is increased (cf. ranges of λ 1 in Figures A1 and A2). We note that the number of time-dependent configurations |n 1 , . . . , n 5 ; t for N = 100 and M = 5 is 4598126 ≈ 4.6 × 10 6 and the total number of optimized and fully time-dependent parameters in these demanding computations is thus 4598126 + 5 × 128 2 × 4 = 4680046 ≈ 4.7 × 10 6 .
To verify this finding about λ 1 in a setting, where a larger number of orbitals can be included in our simulations, we repeated our convergence study in a system with N = 10 particles. To keep the mean-field interaction g 0 (N − 1) and predictions comparable, we investigated a stronger interaction of g 0 = 0.5 in this case. Figures A3 and A4 show the N = 10 results for the first eigenvalue λ 1 of the 1BDM for B = 1 and B = 5, respectively, as a function of time for M ≤ 13.
A complete convergence of the eigenvalue λ 1 (t) is still not achieved; however, the discrepancies for different M are very small, cf. the M = 12 and M = 13 cases in Figure A3 and the M = 11 and M = 12 cases in Figure A4. The tendency to an increased amount of fragmentation and correlations for an increased AMF field strength inferred in the N = 100 case in Figures A1 and A2 is also underpinned by the results for N = 10: B = 5 yields a drastically smaller λ 1 (t) than B = 1 (cf. ranges of λ 1 in Figures A3 and A4). At this stage, the discrepancies may well arise due to the formal problems with the true contact interaction potential [75,76] that we have chosen to use for our present investigation. We note that the number of time-dependent configurations |n 1 , . . . , n 1 3; t for N = 10 and M = 13 is 646646 ≈ 6 × 10 5 and the total number of optimized and fully time-dependent parameters in these demanding computations is thus 646646 + 13 × 128 2 × 4 = 859638 ≈ 9 × 10 5 .
In summary, our convergence study corroborates our main observations about the emergence of fragmentation and correlations in the main text and underpins the accuracy of our MCTDH-X computations.