Analysis of a Trapped Bose–Einstein Condensate in Terms of Position, Momentum, and Angular-Momentum Variance

: We analyze, analytically and numerically, the position, momentum, and in particular the angular-momentum variance of a Bose–Einstein condensate (BEC) trapped in a two-dimensional anisotropic trap for static and dynamic scenarios. Explicitly, we study the ground state of the anisotropic harmonic-interaction model in two spatial dimensions analytically and the out-of-equilibrium dynamics of repulsive bosons in tilted two-dimensional annuli numerically accurately by using the multiconﬁgurational time-dependent Hartree for bosons method. The differences between the variances at the mean-ﬁeld level, which are attributed to the shape of the BEC, and the variances at the many-body level, which incorporate depletion, are used to characterize position, momentum, and angular-momentum correlations in the BEC for ﬁnite systems and at the limit of an inﬁnite number of particles where the bosons are 100% condensed. Finally, we also explore inter-connections between the variances.


Introduction
Bose-Einstein condensates (BECs) made of ultra-cold atoms offer a wide platform to study many-body physics [1][2][3][4][5]. Here, there is a growing interest in the so-called particle limit [6][7][8][9][10][11][12][13][14][15][16], in which the interaction parameter (i.e., the product of the interaction strength times the number of particles) is kept fixed while the number of particles is increased to infinity. At the particle limit, the energy per particle, density per particle, and reduced density matrices [17] per particle computed at the many-body level of theory boil down to those obtained in mean-field theory [7][8][9][10]14,16], despite the fact that the respective many-boson wavefunctions are (much) different [13,15]. It turns out that variances of many-particle operators are a useful tool to characterize correlations (namely, differences between respective many-body and mean-field quantities) that exist even when the interacting bosons are 100% condensed [11,12].
The variance of a many-particle operator of a trapped BEC generally depends on the trap shape, strength and sign of the interaction and, in out-of-equilibrium problems, on time. Consequently, the difference between variances computed at the many-body and mean-field levels of theory also depends on these variables and, of course, on the observable under examination. The first examples [11,12] concentrated on one-dimensional problems and the position and momentum variances, and investigated conditions and mechanisms for the differences between the respective many-body and mean-field variances at the particle limit. In two spatial dimensions, further types of trap topologies come into play, and respective many-body and mean-field variances can exhibit additional phenomena, such as opposite anisotropy [18] and distinct (effective) dimensionality [19]. The many-body variance of a trap BEC has been applied to extract excitations [20], analyze the range of inter-particle interaction [21], examine the effects of asymmetry of a double-well potential [22], and to assess numerical convergence [23,24].
So far, only the position and momentum variances were studied for BECs in rather general traps. In [25,26], the angular-momentum variance is studied for BECs in two-dimensional isotropic traps, and scenarios were the mean-field angular-momentum variance has less [25] or more [26] symmetry (in terms of its conservation) than the many-body angular-momentum variance are identified. Going beyond these works, in the present work we study, analytically and numerically, the angular-momentum variance of a trapped BEC in a two-dimensional anisotropic trap for static and dynamic scenarios, and analyze the difference between the many-body and mean-field variances for finite systems and at the limit of an infinite number of particles. Furthermore, we also study the respective position and momentum variances, and thereby offer a comprehensive characterization of the BEC in terms of its variances. This would allow us to put forward inter-connections between the variances.
Let us elaborate on the strategy of exposition chosen in the paper. We first study the ground state of a many-particle model which is exactly solvable, i.e., integrable, both at the many-body and mean-field levels of theory. A couple of symmetries are also used in the analysis. These would allow us to obtain exact and transparent results for any number of particles and particularly to analyze the variances at the particle limit. The merit of analytical closed-form results and, in the context of interacting bosons, their explicit evaluation at the limit of an infinite number of particles is obvious. Then, as is the usual case in many realistic systems, we continue to explore a set-up which is not integrable, and more so, examine its out-of-equilibrium dynamics which is rather complicated already at the mean-field level of theory, let alone at the many-body level of theory. The later necessitates the state-of-the art numerical tools for the accurate integration of the Schrödinger equation and a careful interpolation of properties to the particle limit. All in all, we show below that the combination of analytics and numerics, i.e., of completely opposite methodologies, provides substantial and complementary novel knowledge on the position, momentum, and angular-momentum variances of anisotropic trapped BECs in two spatial dimensions.
The structure of the paper is as follows. In Section 2 we study the position, momentum, and angular-momentum variances of the ground state within an exactly solvable model, the anisotropic harmonic-interaction model. In Section 3 we study numerically the time-dependent variances of an out-of-equilibrium BEC sloshing in a tilted annulus. Summary and outlook are given in Section 4. Finally, Appendix A discusses translations of variances and inter-connections of the latter.
Transforming from Cartesian to Jacobi coordinates, the many-body solution for the ground state is given by where are the interaction-dressed frequencies of the relative-motion degrees-of-freedom, and are parameters arising in the transformation from Jacoby coordinates back to Cartesian coordinates.
Equation (4) prescribes the range of interactions for which the system is trapped, λ 0 > − ω 2 x 2N , i.e., from moderate repulsion to any attraction. Clearly, the many-body solution (3) in two spatial dimensions factorizes to a product of respective one-dimensional many-body solutions.
All properties of the ground state can in principle be obtained from Ψ, such as the energy, densities, and reduced density matrices, see [29]. Here, as mentioned above, we concentrate on variances and their inter-connections. The many-particle positionX = ∑ N j=1 x j ,Ŷ = ∑ N j=1 y j variance per particle is given by Due to the symmetry of center-of-mass separation in the Hamiltonian (1), the many-particle position variance per particle is independent both of the interaction strength and the number of bosons in the system. Similarly, the many-particle momentumP reflecting the minimal uncertainty product 1 4 of the interacting system in the anisotropic harmonic trap.
The many-particle angular-momentumL variance per particle is, at least for bosons, a less familiar and more intricate quantity. After some lengthy but otherwise straightforward algebra it is given by where we have made use of the bosonic permutational symmetry, the structure of Ψ, and the inverse coordinate transformations to evaluate the various integral terms contributing to (8).
The angular-momentum variance per particle of the ground state (3) depends on the dressed frequencies, Ω x and Ω y , and the number of particles N. Namely, unlike the respective position and momentum variances it depends explicitly on the interaction strength and the number of particles. , the value for a single particle in the anisotropic trap 1 2 ω 2 x x 2 + 1 2 ω 2 y y 2 , which only depends on the trap anisotropy. Opposite to the non-vanishing of the angular-momentum variance, we note that the expectation value of the angular-momentum operator, 1 N Ψ|L Z |Ψ , vanishes for any anisotropy ω y ω y , interaction strength λ 0 , and number of particles N. This is straightforward to see since Ψ is even under reflection of all coordinates X → −X and separately of Y → −Y, whereaŝ L Z is odd under reflection.
The anisotropic harmonic-interaction model (1) can be solved analytically at the mean-field level of theory as well, like in [29], also see [41]. Starting from the ansatz where each and every boson resides in one and the same orbital, the mean-field solution is given by where Λ = (N − 1)λ 0 is the interaction parameter and Λ > − ω 2 x 2 the condition for a trapped solution. Like the many-body solution (3), the mean-field solution (11) in two spatial dimensions factorizes to a product of respective one-dimensional mean-field solutions.
The many-particle position variance computed at the mean-field level is given by and seen to be dressed by the interaction. Similarly, the many-particle momentum variance computed at the mean-field level is dressed by the interaction and given by Interestingly, because the mean-field solution (11) is made of Gaussian functions, it satisfies the minimal uncertainty product 1 as well. The many-particle angular-momentum variance computed at the mean-field level is given by where we have made use of the structure and symmetries of Φ GP , to arrive at the final expression. The relation between the mean-field and many-body variances deserves a discussion. Their difference is used to define position, momentum, and angular-momentum correlations in the system. For the position and momentum variances, the following ratios hold, obviously for any number of particles N. These ratios simply imply that, since repulsion (Λ < 0) broadens the position density, the many-body position variance is smaller than the corresponding mean-field one for repulsive interaction, and vise verse for attraction (Λ > 0). Inversely, since repulsion narrows the momentum density, the many-body momentum variance is larger than the corresponding mean-field one for repulsive interaction, and vise versa for attraction. Furthermore, both the position and momentum variances per particle exhibit the same anisotropies as the respective densities for any interaction parameter Λ, namely, if 1 We shall return to these relations and the anisotropy of the variance in the numerical example below.
We now extend the above discussion to the particle limit, in which the energy per particle, densities per particle, and reduced densities per particle at the mean-field and many-body levels of theory coincide, see in the context of the harmonic-interaction model [16]. Particularly, the system of bosons becomes 100% condensed. The results (16) for the position and momentum variances hold at the particle limit as well, owing to the center-of-mass separability for any number of particles, namely, For the angular-momentum variance the limit has to be taken explicitly for each of the terms in (8). First are the frequencies (4), for which we have at the limit of an infinite number of bosons when Λ is held fixed lim Then, the angular-momentum variance takes on the appealing form Comparing (18) to the mean-field expression (14), it is instrumental to prescribe their ratio at the limit of an infinite number of particles (where, as mentioned above, the density per particle and other properties coincide), which is always smaller than 1 for interacting bosons in the anisotropic trap. Furthermore, we see that for attractive interaction the many-body variance can become much larger than the mean-field quantity in the anisotropic trap, signifying the growing necessity of the many-body treatment, even when the system is 100% condensed. This concludes our investigation of a solvable anisotropic many-boson model in which the variances of the momentum, position, and angular-momentum many-particle operators can be computed and investigated analytically and their values at the many-body and mean-field levels of theory compared and contrasted.

Bosons in an Annulus Subject to a Tilt
In most scenarios of interest, the position, momentum, and angular-momentum variance cannot be computed analytically. This in many cases is the situation when symmetries are lifted. Moreover, even when the variances can be computed for the ground state, like in the previous Section 2, their values for an out-of-equilibrium scenario are rarely within analytical reach. This would be the situation of the present investigation.
Bosons in rings, annuli, and shells have attracted considerable attention . Here we consider weakly interacting bosons initially prepared in the ground state of a two-dimensional annulus. The annulus is then suddenly slightly tilted, leading to an out-of-equilibrium dynamics in an anisotropic setup. We build on and extend the study of bosons' dynamics in an annulus within an isotropic setup [19] (for which, e.g., the angular-momentum variance is 0). We analyze the BEC dynamics in terms of its time-dependent variances and other quantities of relevance, see Figures 1-7 below.
The bosons are initially prepared in the ground state of the annulus, see Figure 1 in [19]. The trap potential is given byV(r) = 0.05r 4 + V 0 e − r 2 2 , with a barrier of heights V 0 = 5 and 10 throughout this work. The interaction between the bosons is repulsive and taken to be where the interaction strengths are λ 0 = 0.02 and 0.04 throughout this work. The form and extant of the interaction potential do not have a qualitative influence on the physics to be described below.
At time t = 0 a linear term is added such that V(r) = 0.05r 4 + V 0 e − r 2 2 + 0.01x. The physical meaning of the added potential is that a constant force pointing to the left is suddenly acting on the interacting bosons. Geometrically, the annulus can be considered to be slightly tilted to the left. Symmetry-wise, the isotropy of the potential is lifted and anisotropy sets in. All in all, the interacting bosons are not in their ground state any more and out-of-equilibrium dynamics emerges.
To compute the time-dependent many-boson wavefunction we use the multiconfigurational time-dependent Hartree for bosons (MCTDHB) method [67][68][69]. MCTDHB represents the wavefunction as a variationally optimal ansatz which is a linear-combination of all time-dependent permanents generated by distributing the N bosons over M time-adaptive orbitals. The quality of the wavefunction increases with M and convergence of quantities of interest is attained. The theory, applications, benchmarks, and extensions of MCTDHB are extensively discussed in the literature, see, e.g., Refs. . Here we employ the numerical implementation in [96,97] both for preparing the ground state [98] (using imaginary-time propagation) and real-time dynamics. Finally, we mention that MCTDHB is the bosonic version of the nearly three-decades-established distinguishable-particle multiconfigurational time-dependent Hartree method frequently used (alongside its extensions) in molecular physics [99][100][101][102][103][104][105].
From the time-dependent wavefunction Ψ(r 1 , . . . , r N ; t), here normalized to 1, we compute properties of interest. The reduced one-particle density matrix is defined as The number of particles residing outside the condensed mode φ 1 (r; t), i.e., the total number of depleted particles, is given by ∑ j>1 n j (t) = N − n 1 (t).
Analogously, the reduced two-particle density matrix is given by ρ(r 1 , r 2 , r 1 , , from which the variance of a many-particle operatorÂ = ∑ jâ (r) is computed, To compute the various terms for the position, momentum, and angular-momentum variance numerically we work in coordinate representation and operate on orbitals first with coordinate derivatives and then with coordinate multiplications. Thus, for the position operatorâ(r) =x and a 2 (r) =x 2 and likewise forâ(r) =ŷ, for the momentum operatorâ(r) = 1 For the numerical solution we use a grid of 64 2 points in a box of size [−8, 8) 8) with periodic boundary conditions. Convergence of the results with respect to the number of grid points has been checked using a grid of 128 2 points.
We begin with the dynamics of N = 10 bosons in the annulus. Following the sudden tilt of the potential, the bosons start to flow to the left. To quantify their sloshing dynamics, Figure 1 shows the time-dependent center-of-mass, 1 N Ψ|X|Ψ (t), for the two barrier heights, V 0 = 5 and V 0 = 10, and the two interaction strengths, λ 0 = 0.02 and λ 0 = 0.04 [we mention that 1 N Ψ|Ŷ|Ψ (t) = 0 due to the Y → −Y reflection symmetry]. The dynamics of 1 N Ψ|X|Ψ (t) appears to be almost periodic and rather simple. We examine the amplitude and frequency of oscillations. It is useful to compare the amplitude of the center-of-mass motion with the radius of the (un-tilted) annulus. The radius of the density at its maximal value, R, is determined numerically using a computation with a resolution of 256 2 grid points as R = 1.75(0) for V 0 = 5, λ 0 = 0.02, and R = 2.06(2) for V 0 = 10, λ 0 = 0.02 [19]. From Figure 1 we see that the amplitude is about 13-25% of the radius, implying a mild sloshing of the density along the tilted annulus. The amplitude increases with the radius of the annulus and decreases with the interaction strength, where the latter implies that it is more difficult to compress the BEC for a stronger interaction. The decrease of the frequency of oscillations with R (V 0 ) and increase with λ 0 are compatible with angular excitations, also see [19]. Last but not least, convergence with M is clearly seen. In fact, here already M = 1 orbitals accurately describe the center-of-mass dynamics for short and intermediate times, and M = 3 orbitals for all times.  Figure 2 depicts the total number of depleted particles, N − n 1 (t), out of N = 10 bosons in the tilted annulus. During the dynamics, the depletion is rather small, ranging from less than 0.012 of a particle out of N = 10 particles (0.12%) for V 0 = 5, λ 0 = 0.02 to less than 0.065 of a particle out of N = 10 particles (0.65%) for V 0 = 10, λ 0 = 0.04. Generally, the depletion increases with the annulus radius and interacting strength, implying angular excitations, see [19]. Finally, convergence with M is clearly seen. Now, M = 3 orbitals nicely follow and M = 5 orbitals accurately describe the depletion dynamics, see  Figure 3 plots the time-dependent many-particle position variance per particle, 1 N ∆ 2X (t) and 1 N ∆ 2Ŷ (t), for the two barrier heights and two interaction strengths. There are several features that immediately are seen. First, since rotational symmetry is lifted, the dynamics of respective quantities along the x-axis and y-axis are different [note that at t = 0 the variances 1 N ∆ 2X = 1 N ∆ 2Ŷ because the initial condition is the ground state of the un-tilted, isotropic annulus]. The mean-field (M = 1) and many-body (M ≥ 3) values are clearly separated from each other, and the former lie about 10-25% above the latter depending on the repulsion strength and barrier height, also see [11,19]. This is despite the small amount of depletion, see Figure 2. Furthermore, the many-body and mean-field variances do not cross each other, see Figure 3, indicating that the dynamics is mild and sufficiently close to the ground state and low-lying manifold of excited states (compare to [18] with interaction-quench dynamics in a single trap).
The mean-field position variance accounts for the geometry of the annulus and shape of the density and weakly depends on the interaction strength. The many-body position variance incorporates the (small amount of) depletion and hence strongly depends on the interaction strength. Both the mean-field and many-body variances oscillate with a relatively small amplitude, albeit with a different frequencies' content, see in this respect [20]. This amplitude slightly decreases with the repulsion strength, which correlates with the dependence of the center-of-mass dynamics on the interaction strength, see Figure 1. Moreover, the amplitude of oscillations of the y-axis variances is smaller than that of the x-axis variances, since the sloshing dynamics is primarily along the x direction. Last but not least is the so-called opposite anisotropy of the (position) variance [18]. During the dynamics, there can occur instances where 1 N ∆ 2X > 1 N ∆ 2Ŷ at the many-body level (M ≥ 3) whereas 1 N ∆ 2X < 1 N ∆ 2Ŷ at the mean-field level (M = 1) [or, in principle, vice versa, i.e., 1 N ∆ 2X < 1 N ∆ 2Ŷ at the many-body level whereas 1 N ∆ 2X > 1 N ∆ 2Ŷ at the mean-field level]. Examples for the former can be readily found for V 0 = 10, λ 0 = 0.02, see Figure 3c,g around t = 70, and for V 0 = 10, λ 0 = 0.04, see Figure 3d,h around t = 100, signifying among others that correlations 'win' over shape. Finally, we see that already M = 3 orbitals accurately describe the dynamics of the position variance.
We move to the momentum variance and also make contact with the results of the position variance. Figure 4 displays the many-particle momentum variance per particle, 1 N ∆ 2P X (t) and 1 N ∆ 2P Y (t), for V 0 = 5, V 0 = 10 and λ 0 = 0.02, λ 0 = 0.04. Just like the results of the position variance, since rotational symmetry is lifted the dynamics of respective quantities along the x-axis and y-axis are different [the initial conditions imply 1 The mean-field (M = 1) and many-body (M ≥ 3) values are, again, separated from each other, but now the former lie below the later, and there is only about 1-4% of a difference depending on the repulsion strength and barrier height, also see [12,19]. Thus, the momentum variance rather weakly depends on the (small amount of) depletion. This is because the matrix elements in (20) are typically smaller with the momentum operator than with the position operator. Yet, despite their small difference, the many-body and mean-field momentum variances do not cross each other, see Figure 4 (contrast with the interaction-quench dynamics in a single trap in [18]).
It is instructive to analyze the momentum-variance dynamics at short times. Whereas ∆ 2P This matches the geometry of the sloshing dynamics in the tilted annulus, in which bosons from the 'north' and 'south' poles (on the y-axis) start to move to the left and accumulate in the 'west' pole (on the x-axis), and that the cross section of the rim of an annulus is enlarged when moving away from the center of the annulus. In other words, the dynamics of the momentum variances at short times when moving to the left reflects the relative localization of the bosons in the x direction and the effective broadening of the wavepacket along the y direction. Both the mean-field and many-body variances oscillate with a very small amplitude, note the scale on the y-axis in Figure 4. The high-frequency oscillations mark high-energy radial excitations across the (tight) annulus rim [19]. Like for the position variance, the amplitude of oscillations of the y-axis momentum variances is smaller than that of the x-axis momentum variances. Finally, we see that already M = 3 orbitals accurately describe the dynamics of the momentum variance; the difference to the M > 3 results is lower than 1%.
We now move to the angular-momentum variance and an interesting inter-connection with the momentum variance. Figure 5 presents the many-particle angular-momentum variance per particle, (t) appears to be almost periodic and rather regular, more than that for the respective position and momentum variances, compare to Figures 3 and 4. On the other end, focusing on the dynamics of the center-of-mass in Figure 1, one can clearly observe correlation between the two quantities; Whenever 1 N Ψ|X|Ψ (t) has a minimum, i.e., the bosons are maximally localized to the left, 1 N ∆ 2 L Z (t) has a maximum, and whenever 1 N Ψ|X|Ψ (t) has a maximum (which value is about 0), i.e., the bosons are momentarily, approximately equally distributed along the annulus, 1 N ∆ 2 L Z (t) has a minimum (which value, as mentioned above, is close to 0). Furthermore, the frequencies of the two quantities as well as their relative amplitudes as a function of the barrier height and interaction strength are alike. These observations call for a dedicated analysis.     To shed light on the above dynamics of the angular-momentum variance, see Figure 5, we analyze the translational properties of variances in Appendix A. Whereas the position variances and, trivially, the momentum variances, are translationally invariant, this invariance does not hold for the angular-momentum variance. If a wavepacket prepared in the origin has angular-momentum variance 1 N ∆ 2 L Z , then several terms are added when the wavepacket is translated to the point (a, b) in plane, and angular-momentum variance is thereafter computed, see Equation (A3). Now, if this wavepacket is rotationally symmetric, i.e., 1 N ∆ 2 L Z = 0, then several of the terms in (A3) vanish due to spatial symmetry and we are left with the appealing relation, 1 , connecting the angular-momentum variance of Ψ(a, b) localized at (a, b) and of Ψ at the origin. The meaning of this relation is that the momentum variances, 1 N ∆ 2P X and 1 N ∆ 2P Y , together with the spatial translations along the y-axis and x-axis, respectively, determine the angular-momentum variance of a translated wavepacket (rotationally symmetric at the origin).
Returning to and combining Figure 5 for the angular-momentum variance, Figure 1 for the center-of-mass dynamics, and The momentum variance helps us understand the deviations between the many-body and mean-field results in Figure 5. We see that the maxima of the many-body 1 N ∆ 2 L Z (t) (M > 3) are larger than the maxima of the mean-field 1 . The difference is about 7-25% (compare to the low depletion, Figure 2), depending on V 0 and λ 0 , and follows the respective trend of the many-body and mean-field results for 1 N ∆ 2P Y , see Figure 4e-h. We note that, although the wavepacket describing the bosons dynamics in the tilted annulus is not a translated, rotationally invariant wavepacket, and the values of deviations (in percents) between the many-body and mean-field results are actually larger for ∆ 2 L Z (at the maxima) than for ∆ 2P Y , we find the above analytically based analysis to well explain the numerical findings and trends. Last but not least, a close inspection of the many-body and mean-field curves of the angular-momentum variance in Figure 5 shows that there are instances when they cross each other, i.e., one is smaller or larger than the other. This is in contrast with the non-crossing of the many-body and mean-field position and momentum variances, see Figures 3 and 4, respectively. Finally, we find that already M = 3 time-adaptive orbitals accurately describe the dynamics of the angular-momentum variance.
Our investigations are nearing their end, what is left to explore is the behavior of the position, momentum, and angular-momentum variances at the particle limit. Which of the above-described detailed findings, plotted in Figures 1-5 for a rather small (N = 10 bosons) yet weakly depleted BEC, survive this limit? To answer the question, we concentrate on the system with the higher barrier, V 0 = 10, and stronger interaction (for N = 10 bosons), λ 0 = 0.04. We hence fix the interaction parameter Λ = λ 0 (N − 1) = 0.36, and compute and compare the dynamics for N = 10, N = 100, and N = 1000 bosons using M = 3 time-adaptive orbitals. We have seen for N = 10 bosons that M = 3 time-adaptive orbitals accurately describe the variances. This implies that, keeping the interaction parameter Λ fixed while increasing the number of particles N, using M = 3 time-adaptive orbitals for calculating the variances will be (at least) as accurate as for N = 10 particles, see in this respect [24]. Before we proceed, a methodological remark. Examining the convergence of properties with the number of particles for N = 10, N = 100, and N = 1000 bosons is (still) far away from infinity, see in this respect [15]. We hence use, interchangeably, the term en route to the particle limit. We shall see below that, in effect, the particle limit is practically well achieved for the variances already for N = 1000 bosons. Figure 6 prints the total number of depleted particles, N − n 1 (t), for N = 10, N = 100, and N = 1000 bosons for Λ = 0.36 and V 0 = 10 using M = 3 time-adaptive orbitals. Convergence of the number of depleted particles with N is nicely seen. Since N − n 1 (t) converges to a finite (and small) value with N, the bosons are becoming 100% condensed in the limit of an infinite number of particles, i.e., n 1 (t) N → 1 as N → ∞, at least up to the maximal time of the computation, t = 100. Figure 7 exhibits the position variances per particle, 1 N ∆ 2X (t) and 1 N ∆ 2Ŷ (t), momentum variances per particle, 1 N ∆ 2P , angular-momentum variance per particle, 1 N ∆ 2 L Z (t), and the expectation value of the center-of-mass, 1 N Ψ|X|Ψ (t), for N = 10, N = 100, and N = 1000 bosons and for Λ = 0.36 and V 0 = 10 using M = 3 time-adaptive orbitals. Once again, convergence of each of the quantities with N is clearly seen. Yet, whereas the center-of-mass dynamics converges to the mean-field dynamics when the number of particles is increased, the variances exhibit many-body dynamics which converges nicely with N, but not to the respective mean-field dynamics. Beyond that, all the above results, for the frequencies, amplitudes, anisotropies, inter-connections, and particularly the differences between the many-body and mean-field position, momentum, and angular-momentum variances persist at the limit of infinite number of particles, despite the bosons becoming 100% condensed. This brings the present analysis to an end.  Figure 6. Depletion dynamics following a potential quench en route to the particle limit. The time-dependent total number of depleted particles, N − n 1 (t), of N = 10, N = 100, and N = 1000 bosons with interaction parameter Λ = λ 0 (N − 1) = 0.36 for an annulus with barrier height V 0 = 10 following a sudden potential tilt by 0.01x. The number of time-adaptive orbitals is M = 3. The respective position, momentum, and angular-momentum variances along with the expectation value of the center-of-mass are plotted in Figure 7. See the text for more details. The quantities shown are dimensionless.   Figure 7. Position, momentum, and angular-momentum variance dynamics following a potential quench en route to the particle limit. The mean-field (M = 1 time-adaptive orbitals) and many-body (using M = 3 time-adaptive orbitals) time-dependent position variances per particle, (a) 1 N ∆ 2X (t) and , and angular-momentum variance per particle, (f) 1 N ∆ 2 L Z (t), of N = 10, N = 100, and N = 1000 bosons with interaction parameter Λ = λ 0 (N − 1) = 0.36 for an annulus with barrier height V 0 = 10 following a sudden potential tilt by 0.01x. (e) The time-dependent expectation value of the center-of-mass, 1 N Ψ|X|Ψ (t). The respective depletions are plotted in Figure 6. See the text for more details. The quantities shown are dimensionless.

Summary and Outlook
In the present work we studied, analytically and numerically, the position, momentum, and especially the angular-momentum variance of interacting bosons trapped in a two-dimensional anisotropic trap for static and dynamic scenarios. Explicitly, we investigated the ground state of the anisotropic harmonic-interaction model in two spatial dimensions analytically and researched the out-of-equilibrium dynamics of repulsive bosons in tilted two-dimensional annuli numerically accurately by using the MCTDHB method. The differences between the variances at the mean-field level, which are attributed to the shape of the density per particle, and the respective variances at the many-body level, which incorporate a small amount of depletion outside the condensed mode, were used to characterize sometimes large position, momentum, and angular-momentum correlations in the BEC for finite systems and at the limit of an infinite number of particles where the bosons are 100% condensed. Finally, we also explored and utilized inter-connections between the variances, particularly between the angular-momentum and momentum variances, through the analysis of their translational properties.
There are many intriguing directions to follow out of which we list three below. First, variances of BECs in the rotating frame of reference in which high-lying excitations become low-energy excitations and even the ground state. Second, angular-momentum variance of a BEC flowing past an obstacle in which the mean angular-momentum variance vanishes. And third, variances in three-dimensional geometries lacking lower-dimensional analogs, such as a Möbius strip. In all these cases, whether considering a few interacting bosons or a BEC in the particle limit, interesting and exciting results are expected.

Appendix A. Variances and Translations
Consider the many-particle translation operator in two spatial dimensions e −i(P X a+P Y b) , wherê P X = ∑ N j=1p x,j andP Y = ∑ N j=1p y,j . Its operation on a multi-particle wavefunction Ψ is given by e −i(P X a+P Y b) Ψ(x 1 , y 1 , . . . , x N , y N ) = Ψ(x 1 − a, y 1 − b, . . . , x N − a, y N − b) ≡ Ψ(a, b). What are the implications on the variances when computed with respect to the translated wavefunction Ψ(a, b)?
From Ψ(a, b)|L Z |Ψ(a, b) Equation (A3) deserves a discussion. In turn, even for the ground state of an interacting many-boson system in a rotationally symmetric [for which 1 N ∆ 2 L Z = 0 holds] but otherwise translated trap, the angular-momentum variance differs at the many-body level and mean-field level of theory, i.e., when a, b = 0 and λ 0 = 0. This is, as can be seen in (A4), because of the respective many-body and mean-field momentum variances, The analytical result (A4) is employed to analyze the numerical findings for the time-dependent angular-momentum variance in the main text. Generally in the absence of spatial symmetries, see Equation (A3), more terms contribute to the translated angular-momentum variance.