Poincar\'e index formula and analogy with the Kosterlitz-Thouless transition in a non-rotated cold atom Bose-Einstein condensate

A dilute gas of Bose-Einstein condensed atoms in a non-rotated and axially symmetric harmonic trap is modelled by the time dependent Gross-Pitaevskii equation. When the angular momentum carried by the condensate does not vanish, the minimum energy state describes vortices (or antivortices) that propagate around the trap center. The number of (anti)vortices increases with the angular momentum, and they repel each other to form Abrikosov lattices. Besides vortices and antivortices there are also stagnation points where the superflow vanishes; to our knowledge the stagnation points have not been analyzed previously, in the context of the Gross-Pitaevskii equation. The Poincar\'e index formula states that the difference in the number of vortices and stagnation points can never change. When the number of stagnation points is small, they tend to aggregate into degenerate propagating structures. But when the number becomes sufficiently large, the stagnation points tend to pair up with the vortex cores, to propagate around the trap center in regular lattice arrangements. There is an analogy with the geometry of the Kosterlitz-Thouless transition, with the angular momentum of the condensate as the external control parameter instead of the temperature.


Introduction
A Bose-Einstein condensate is a coherent macroscopic quantum state with unique features that facilitate a high level of experimental control [1][2][3]. In particular condensates that form in a trapped, rotated gas of ultra-cold alkali atoms are studied extensively, both in earthbound and in earth-orbiting laboratory experiments [4]. Among the emerging applications is the development of ultra-sensitive sensors and detectors [5]. The properties of cold atom condensates are also under active investigation, as a potential platform for quantum computation and simulation [6]. At the level of the mean field theory a condensate can be modelled by a macroscopic wave function that solves the time-dependent Gross-Pitaevskii equation [7,8]. This is a nonlinear Schrödinger equation with a quartic nonlinearity that accounts for the interactions between atoms [9][10][11][12][13]. Quantum vortices are the principal topological excitations of a trapped cold atom Bose-Einstein condensate [14][15][16][17][18]. They have been studied extensively in terms of solutions of a stationary Gross-Pitaevskii equation that models a uniformly rotating condensate [12,[19][20][21][22], for detailed reviews, see e.g., [23,24]).
The time-dependent Gross-Pitaevskii equation appears in a broad range of physical systems, beyond condensed matter physics. For example, in combination with a gravitational trapping potential it can model Dark Matter as a Bose-Einstein condensate [25,26]. Using the Madelung transformation, the equation can be sent into a system of equations that describe a barotropic-type fluid, which can model codimension two membranes in higher dimensions [27]. Furthermore, in two and three space dimensions and in the limit of constant and uniform density, one arrives at the fluid dynamical Euler equation that supports pointlike vortex filaments, with dynamics that is described by the localized induction approximation [28]. Finally, in two dimensions these point vortices also appear in an effective field theory that evaluates a conformal field theory operator spectrum at large global charge [29].
Topological methods can be highly effective in the characterization of a physical system, when the goal is to identify and describe robust phenomena. The numerous applications range from fundamental interactions to condensed matter, fluid dynamics and beyond [30].
Here we extend the repertoire, in the case of a cold atom Bose-Einstein condensate, by the Poincaré index formula [31][32][33]. It is a refinement of the Poincaré-Hopf index theorem [34] that is used widely to analyze topological aspects of vectors fields. We employ the Poincaré index formula to identify and analyze previously unreported properties of quantum vortices in the context of cold atom Bose-Einstein condensates. In particular, we identify the presence of a changeover transition that displays a geometry that is quite similar to the geometry of vortices and antivortices in the Kosterlitz-Thouless transition that occurs e.g. in the two dimensional XY-model and in the plane rotor model [35][36][37]. Instead of the temperature, now the external control parameter is the internal angular momentum that is supported by the condensate; in the case of a non-rotating trap its value can be regulated externally, possibly by optical photon beams.
A vortex in the two dimensional Gross-Pitaevskii model has a pointlike core where the modulus of the wave function vanishes. The energy of a trapped vortex is finite, albeit logarithmically increasing in the length scale that characterizes the size of the trap. The vortex is a topological construct with an integer valued circulation of the superflow velocity v(x, t) around its core; in the case of a single vortex core the circulation has the normalized value ±1 where the sign depends on the orientation of the circulation, and the total circulation of a multivortex configuration is the sum of its individual vortex circulations. Commonly the circulation around the vortex core is the only topological invariant that is assigned to a vortex system, in the case of a cold atom Bose-Einstein condensate.
Here we consider the free energy minima of the two dimensional Gross-Pitaevskii model, in the case of a nonrotating harmonic trap. We observe that besides the circulation, there are also other topological structures that can be assigned to a multivortex configuration. In particular, besides the vortex cores that are topologically centers of the vector field v(x, t), there can also be saddles of v(x, t). These are points where the flow velocities around distinct vortices exactly cancel each other, and thus where the superflow vector field stagnates. Unlike in the case of a vortex core, the density of the condensate does not vanish at a stagnation point. Nevertheless, the stagnation points are also topological structures that can not be removed by a generic local deformation of the wave function. Moreover, unlike in the case of a vortex core the circulation of v(x, t) can not directly detect the presence of a stagnation point. But all the fixed point structures of v(x, t) can be assigned an integer valued winding number that is a topological invariant with a content that can be summarized by the Poincaré index formula [31][32][33].
In the common presentation of the Kosterlitz-Thouless transition in the case of the XY-model and the plane rotor model [37], an antivortex has the topology of a saddle; see Appendix for a detailed description. But to our knowledge the properties of the saddles where the superflow stagnates have not been previously analyzed, in the case of cold atoms condensates. Accordingly we apply the Poincaré index formula to analyze both the vortex cores and the superflow stagnation points. To illustrate this, we use a particular set of solutions of the Gross-Pitaevskii equation, those with a minimal energy for fixed values of the particle number and angular momentum, but in a trap that does not rotate. Our methods and results should help to clarify the different role of the circulation and the winding number as complementary integer valued topological invariants, to characterize a multivortex structure in the case of a Bose-Einstein condensate. For this we start and introduce the Poincaré index formula that is relevant for the description of the topology of the velocity superflow in two spatial dimensions. We apply this formula to analyze a set of multivortex solutions of the Gross-Pitaevskii equation that describes a dilute, non-rotating condensate of cold atoms in an harmonic trap. In particular, we observe and describe in detail, using the Poincaré index formula, how the geometric pattern of vortex cores and stagnation points undergo a transition that is analogous to the Kosterlitz-Thouless transition of vortices and antivortices e.g. in the common presentation of the XY-model, but with the value of angular momentum as the external control parameter instead of temperature; in the Appendix we compare the topological structures, as they commonly appear in the two cases.

Circulation
In a theoretical approach, a trapped Bose-Einstein condensate of cold atoms is commonly modelled by a macroscopic, complex-valued and normalizable wave function ψ(x, t). The modulus |ψ(x, t)| describes the density of the condensate and its phase determines the velocity vector field of the superflow v(x, t) = ∇ arg[ψ](x, t) . (2.1) The principal topological excitations in these condensates are vortices [18]. In a two dimensional model, the core of a vortex is a center of the superflow velocity v(x, t); at the vortex core the modulus of the wave function vanishes. As a topological invariant, a vortex is commonly characterized by the integer valued circulation of the superflow velocity around its core. For a simple closed, counterclockwise curve Γ a single vortex core contributes ±1 to the circulation (2.2), with sign depending on the direction of v(x, t) along the curve; the right-hand panel of Figure 1 shows an example of a vector field v(x) with n v (p; Γ) = +1. In the case of several vortex cores in the interior of the curve, the circulation (2.2) counts their net number. In addition of vortex cores which are centers of the vector field v(x, t), there are also other topologically invariant fixed points that can be assigned to a two dimensional vector field including sinks, sources and saddles. From the knowledge of all its fixed points, a phase portrait of the vector field v(x, t) can then be constructed. It characterizes the entire vector field v(x, t), in a manner that is invariant under continuous local deformations. A phase portrait commonly describes adequately the solutions of the underlying dynamical equation, often without any need to actually solve the equation.

Winding number
In addition of vortex cores that are mathematically centers of the vector field v(x, t), stagnation points of v(x, t) can also be present in cold atom Bose-Einstein condensates. Stagnation points are mathematically saddles of v(x, t). At a stagnation point the superflow vanishes as the velocity (2.1) of distinct vortices exactly cancel each other while the density of the condensate does not need to vanish. But unlike centers, saddles can not be captured by the circulation (2.2). For example, the middle panel of Figure 1 shows a vector field with a saddle. The index (2.2) vanishes when evaluated over any simple counterclockwise curve Γ that encircles the saddle.
Besides the circulation (2.2), as a differentiable two dimensional vector field, the superflow velocity (2.1) also supports the integer valued winding number Here Γ is a simple closed curve on the (x, y)-plane that does not pass through any fixed point of v(x). For a curve that encircles a single isolated center p once in the counterclockwise direction, the winding number is i v (p; Γ) = +1 independently of the sign of (2.2). For a source and a sink, the winding number is also one. A saddle is the only fixed point with a negative winding number. For a single isolated saddle the winding number is i v (p; Γ) = −1; see Appendix.

Poincaré index formula
All the fixed points of a vector field v(x) can be detected by the Poincaré index formula. It states that for any simple closed curve Γ that encloses finitely many fixed points p 1 , p 2 , · · · , p k of v(x) in a counterclockwise direction, the sum of their winding numbers equals the following index [31][32][33] Here X Γ is the Euler characteristic of the area that is bounded by the curve Γ, and I Γ is the number of concave curve segments with tangencies that are internal to the area that is bounded by Γ, and E Γ is the number of convex curve segments where the tangencies are external to Γ; see Figure 1. The Poincaré index formula is a generalization of the Poincaré-Hopf index theorem [34]. The latter assumes, that the vector of v(x) always points towards the outward (or inward) normal direction along the boundary curve of a region. We remind that for any disk the Euler characteristic is X Γ = 1. , and the cyan lines represent its streamlines. The first two panels, respectively, show a center that corresponds to a single vortex core with n v = 1 and i v = 1, and a saddle that corresponds to a stagnation point with n v = 0 and i v = −1. The rightmost panel illustrates what are internal and external tangencies. Since the curve segment between i 1 and i 2 is concave, it ads one to I Γ . Likewise, the curve segment between e 1 and e 2 is convex and it ads one to E Γ . It is notable that both centers and saddles appear in the common description of Kosterlitz-Thouless theory of phase transitions in the 2D XY model. In that model vortex cores appear topologically as centers and antivortex cores appear as saddles, and both have an energy that diverges logarithmically when the short distance cutoff goes to zero [37] (See Appendix.) Besides the superflow velocity v(x, t) we introduce another vector field, that depends only on the modulus of the complex wave function ψ(x, t): Unlike the phase, the modulus of a normalizable wave function is a smooth, single valued and strictly non-negative function that vanishes both at x → ∞ and at vortex cores. Consider any simple closed trajectory around the origin, with a sufficiently large radius to encircle all the critical points of |ψ(x, t)| but so that the vector field (2.5) on the trajectory does not vanish. Along this trajectory the vector field (2.5) always points towards its interior, so that the winding number (2.3) has the value i w = +1, independently of the number of vortex cores that are encircled (note that the circulation (2.3) of this vector field vanishes). A single vortex core is an isolated minimum value critical point of |ψ(x, t)|, thus it is a source of the vector field (2.5) and as such it contributes i w = +1 to the index (2.4). Since a stagnation point is the only fixed point of a vector field with a negative winding number i w = −1, each vortex core must always be balanced by an accompanying stagnation point of |ψ(x, t)|. Note that for certain vortex configurations there can be degeneracies in the fixed point structure of (2.5). This is the case e.g. when a single vortex is located at the center of the trap, and the vortex core is encircled by a nodal line where the modulus |ψ(x, t)| has a maximum value. The vector field (2.5) then has a degeneracy circle and instead of (2.3) the appropriate index theoretic analysis is based on the more general concept of a degeneracy index introduced in [38].

Gross-Pitaevskii equation
We proceed to apply the Poincaré index formula (2.4) to characterize the minimal energy configurations of a trapped Bose-Einstein condensate of N alkali atoms, with specified values of the particle number and of the macroscopic angular momentum. The trap is axially symmetric and non-rotating, and without loss of generality, we choose the trapping potential to be harmonic V (x) = |x| 2 /2. The set-up can model an anisotropic three dimensional trap, resulting in an oblate, essentially two dimensional spheroid condensate D. The atoms interact with each other via a repulsive short range pair potential, so that in the limit where N becomes large, the condensate can be modelled by a solution of the time-dependent two dimensional Gross-Pitaevskii equation. This is a nonlinear Schrödinger equation for the complex valued wave function ψ(x, t) that describes the condensate as a coherent macroscopic quantum state.
With an appropriate choice of various scales and parameters, the relevant two dimensional time-dependent Gross-Pitaevskii equation reads as follows [13]: The dimensionless coupling g specifies the strength of the pairwise interatomic interactions.
In a typical experiment with 10 4 -10 6 atoms its numerical values are g ∼ 10 1 -10 3 . For a detailed discussion of the conventions used here, see [39]. We are interested in the ground state solution of (2.6) that minimizes the Gross-Pitaveskii free energy Since F is a strictly convex functional, its only critical point is the absolute minimum ψ(x) ≡ 0. But in this case there is no cold atom condensate in the trap. Thus, in the presence of any cold atom condensate, the free energy can not have any critical point. As a consequence it follows that, as detailed in [40], the ground state solution of the Gross-Pitaevskii equation can be considered "time-crystalline" [41] provided there are conserved quantities: Besides the free energy F the time evolution (2.6) conserves two additional quantities, as Noether charges. One of these is is also conserved. The numerical value l z of the canonical angular momentum of the condensateL with θ the polar angle around the z-axis, is a free parameter that can take an arbitrary value; see [39]. The r.h.s. of (2.6) can not vanish unless ψ ≡ 0, so that for non-vanishing values of the Noether charges, the minimum of F can not be a critical point of F . Thus, to construct the ground state wave function of (2.6), we use methods of constrained optimization and minimize the free energy (2.7) at the fixed values of the Noether charges (2.8) and (2.9). The Lagrange multiplier theorem [42] states that the minimum of (2.7) can be found as a critical point of where the Lagrange multipliers λ N , λ z respectively enforce the values N = 1 and L z = l z of the Noether charges. The critical points of F λ obey together with the two conditions (2.8) and (2.9). From these three equations we solve for the critical point wave function ψ cr (x) and for the ensuing Lagrange multiplier values λ cr N , λ cr z . Note that both Lagrange multipliers are time independent [40], and that they cannot vanish simultaneously.
We focus on the critical points of (2.11) that are also minima of the free energy (2.7). Let {ψ min (x), λ min N , λ min z } denote such a configuration. If ψ min (x) is an initial value of the Gross-Pitaevskii equation (2.6), then it obeys the linear time evolution We recall (2.10), denote σ = λ min N /λ min z , and define This is a vector field on the plane with a center at the origin (x, y)=(0,0). The winding number (2.3) is i A = +1 for any trajectory that encircles the origin once in counterclockwise direction; we note that A θ is akin the azimuthal component of the vector potential of a line vortex along the z-axis. The time evolution (2.13) can further be written as (2.14) In the presence of A θ the rotations around the z-axis i.e. changes in the polar angle θ are generated by the covariant angular momentum operatorL In the following we are interested in the consequences of the Poincaré index formula (2.3), in the case of a vector field v(x, t) that describes the velocity of the superflow (2.1) of the minimal energy solution of the Gross-Pitaevskii equation (2.6)) with initial condition ψ(x, t = 0) = ψ min (x). We focus on the topological and geometrical properties of its fixed point structure, as the value of the angular momentum l z is changed. For this we analyze the profile of the superflow vector field, by employing the Poincaré index formula in combination with various different trajectories Γ.

Numerical analysis
We have numerically constructed the critical points of (2.11) that minimize the Gross-Pitaevskii free energy (2.7). For clarity we choose positive l z > 0 and with g = 5, 100 and 400. Since the topological properties can not change with continuous changes of g, the results are displayed only for g = 400. (For detailed results with other values of g, see [39]). The problem is discretized within a finite-element framework [43], and the constrained optimization problem is solved using the Augmented Lagrangian Method; numerical methods are discussed in details in [39]. Whenever l z = 0 the minimum free energy solution is a configuration with vortices and superflow stagnation points; the former are centers of the vector field v(x), and the latter are saddles of v(x). The ensuing timeevolution (2.6), (2.13) is timecrystalline, as it always depends on the time variable t in a nontrivial fashion [39].
The number of vortices increases with l z . Each new vortex core enters the condensate at the boundary of the disk D, and moves towards the trap center as l z increases. Since the Euler characteristic of a disk is X Γ = 1, the Poincaré index formula ensures that, on the entire trapping disk D, there is always one more vortex core i.e. center than there are stagnation points i.e. saddles. That is, when the closed curve Γ coincides with the perimeter of the entire disk D in counterclockwise direction, the winding number (2.3) always has the value i v = 1.

Small angular momentum
The Figure 2 shows examples of vortices and stagnation points for different values of the macroscopic angular momentum l z . The corresponding solutions of Eq. (2.6), describe how these structures rotate uniformly around the symmetry axis, with constant angular velocity −λ min z ; see Eq. (2.14). For 0 < l z < 1 the condensate features a unique vortex core that is located off the trap center. This is the regime l z = 0.8 of Fig. 2, which features a single eccentric vortex core as the unique fixed point of v(x). Next, when l z = 1.4, the minimal energy wave function features two additional fixed points of v(x), a vortex core and a stagnation point. Upon increasing l z , the second vortex core together with the stagnation point move toward the trap center, until they form a symmetric pair of vortex cores with a single stagnation point in between (not shown). With a further increase of l z , a third vortex core and a second stagnation point of v now appear, as can be seen when l z = 2.0. Finally, for l z = 2.5 the third vortex core enters closer to the center of the condensate, to form a symmetric triangle. This is accompanied by a merging of the two stagnation points into a higher degree saddle; there is a degenerate saddle with winding number i v = −2 at the center of the disk and it is surrounded by the three vortex cores as centers.

A Kosterlitz-Thouless analogy
The Figure 3 shows the evolution of the minimum energy state, when l z is further increased: We observe a remarkable transition in the geometric pattern of vortex cores and stagnation points; the topology closely resembles the geometry of vortices and antivortices  When the value of l z is further increased, we observe that the pattern identified in Figures 3 are repeated: The difference between the vortex cores and stagnation points is always one. When l z increases the new vortex cores and stagnation points enter the disk together, as tightly bound pairs, and the pairs form structures that are akin co-centric Abrikosov lattices. With the additional structure, that as dictated by the Poincaré index formula there is commonly either a single vortex core, or a triplet of two vortex cores and a single stagnation point near the center of the trap.

Degeneracy circle
Finally, the figure 4 summarizes the main features that can be observed for the vector field w(x) (2.5), that is associated with density gradients, in the case of a single vortex: The configuration l z = 0.6 in the Figure features a single eccentric vortex core as a source, and a sink of w(x). The Poincaré index formula dictates that there is also a saddle. In the case where l z = 1, there is a nodal line encircling the source. This is a degeneracy circle with winding number i w = −1 [38].

Conclusion
We have used the Poincaré index formula in combination with numerical simulations, to study the topology and geometry of solutions to the two dimensional Gross-Pitaevskii equation at different angular momentum values. The minimum energy solution of the equation models the ground state of a cold atom Bose-Einstein condensate in an axially symmetric disk-like, non-rotating harmonic trap. We have varied the angular momentum that is supported by the ground state wave function, and confirmed that the difference in the number of vortex cores and in the number of stagnation points is always equal to one. In particular we have observed how the vortex cores and the stagnation points enter the condensate, always in pairs, when the angular momentum increases. We have observed how the vortex cores repel each other to form co-centric Abrikosov lattices, while for small values of the angular momentum the individual stagnation points tend to aggregate into higher degree stagnation ipoints at the trap center. But when the angular momentum increases, there is a transition and the stagnation points start to pair up with vortex cores. The ensuing structures always rotate around the trap center at a constant angular velocity.
This transition that occurs when the angular momentum increases has a geometric pattern that closely resembles the geometric pattern of the Kosterlitz-Thouless transition in e.g. two dimensional XY-model; see Appendix. But note that in the case of the standard XY-model, the energy of the vortices and antivortices diverges when the cutoff scale for the core size vanishes. In the present case the energy of the individual vortices and stagnation points is always finite but increases logarithmically in the length scale that characterizes the range of the confining trap potential. Furthermore, while in the constant density approximation the core contribution to the energy of vortices diverges logarithmically, in the case of stagnation points their core contribution to energy remains finite even in this limit. Moreover, instead of temperature the external control parameter is the angular momentum l z , such that at large value of l z the geometrical pattern resembles that of small temperature XY-model, and for small values of l z the topological charge structure resembles that of the XY-model high temperature phase. It remains to be clarified whether the present transition is truly a topological phase transition [44]. On the other hand, vortices and antivortices are both centers of the superflow vector field v in (2.1). Note that the phase ϕ is defined up to a constant value, and thus v XY is defined up to global rotations. It follows that a counterclockwise center in v XY as shown in Table II can be continuously rotated to a source, a sink or a clockwise center. The velocity v of the superflow ( In the Table 2 we compare the vector fields v XY (A1) and v (2.1) for a vortex-vortex pair, and for a vortex-antivortex pair. Unlike the XY-vector field (A1) for which there is no saddle-point in the case of a vortex-vortex pair, the superflow vector field v has a saddle-point. Remarkably the phase portrait of v XY , in the case of a vortex-vortex pair, has the same topology as the vector field v (2.1) in the case of a vortex-antivortex pair.
Moreover, the XY phase portrait in the case of a vortex-antivortex pair is similar to the vortex-vortex pair phase portrait of (2.1), but with one of the vortices (the one on the left in the Figure) moved to infinity.