The nonlinear Dirac equation in Bose-Einstein condensates: Vortex solutions and spectra in a weak harmonic trap

We analyze the vortex solution space of the $(2 +1)$-dimensional nonlinear Dirac equation for bosons in a honeycomb optical lattice at length scales much larger than the lattice spacing. Dirac point relativistic covariance combined with s-wave scattering for bosons leads to a large number of vortex solutions characterized by different functional forms for the internal spin and overall phase of the order parameter. We present a detailed derivation of these solutions which include skyrmions, half-quantum vortices, Mermin-Ho and Anderson-Toulouse vortices for vortex winding $\ell = 1$. For $\ell \ge 2$ we obtain topological as well as non-topological solutions defined by the asymptotic radial dependence. For arbitrary values of $\ell$ the non-topological solutions are bright ring-vortices which explicitly demonstrate the confining effects of the Dirac operator. We arrive at solutions through an asymptotic Bessel series, algebraic closed-forms, and using standard numerical shooting methods. By including a harmonic potential to simulate a finite trap we compute the discrete spectra associated with radially quantized modes. We demonstrate the continuous spectral mapping between the vortex and free particle limits for all of our solutions.


Abstract
We analyze the vortex solution space of the 2 1 + ( )-dimensional nonlinear Dirac equation for bosons in a honeycomb optical lattice at length scales much larger than the lattice spacing. Dirac point relativistic covariance combined with s-wave scattering for bosons leads to a large number of vortex solutions characterized by different functional forms for the internal spin and overall phase of the order parameter. We present a detailed derivation of these solutions which include skyrmions, halfquantum vortices, Mermin-Ho and Anderson-Toulouse vortices for vortex winding we obtain topological as well as non-topological solutions defined by the asymptotic radial dependence. For arbitrary values of ℓ the non-topological solutions include bright ring-vortices which explicitly demonstrate the confining effects of the Dirac operator. We arrive at solutions through an asymptotic Bessel series, algebraic closed-forms, and using standard numerical shooting methods. By including a harmonic potential to simulate a finite trap we compute the discrete spectra associated with radially quantized modes. We demonstrate the continuous spectral mapping between the vortex and free particle limits for all of our solutions.

Introduction
Vortices appear in physical settings which span a wide range of energy scales and disciplines ranging from hurricane phenomena [1] and turbulent flow in classical fluids [2], to superfluidity in He 4 [3] and general quantum fluids [4,5], down to the smallest length scales of Bogomol'nyi-Prasad-Sommerfield states in supersymmetric field theories [6]. Vortices are relevant from a technological standpoint in quantum computing for example [7], as well as in more theoretical areas of research such as galactic halos in Bose-Einstein condensate (BEC) theories of dark matter [8]. In BECs vortices as stable rotating solutions of the nonlinear Schrödinger equation (NLSE) are ubiquitous [9,10]. Indeed, the presence of a persistent quantized rotation may be considered a defining property of superfluidity. Although a variety of types of vortices are possible in atomic condenstates [11], the Laplacian inherent to the Schrödinger Hamiltonian constrains the number of observable vortex structures in the BEC. For instance, a BEC made up of spin-F bosons gives rise to a F 2 1 + ( )-component order parameter which nevertheless solves a multi-component NLSE. An alternative method of constructing a multicomponent BEC is to add the internal degrees of freedom by placing a condensate in a two-dimensional honeycomb optical lattice [12,13]. Atoms are condensed into the lowest energy Bloch state then translated to a corner of the Brillouin zone using laser assisted Bragg scattering [14]. At wavelengths large compared to the lattice constant, the microscopic details of the lattice manifest as an additional SL C 2, ( ) spin group symmetry, i.e., a Dirac point structure emerges [12].
In a previous paper we derived a BEC version of the nonlinear Dirac equation (NLDE) for the case of weak interparticle interactions [12]. This has attracted attention from diverse fields of research . Solutions of the NLDE are effectively long wavelength limit lattice envelopes, and so can cover a large number of sites. The healing length in our effectively quasi-2D system is typically about 10 times the lattice constant, and we require the healing length to be small compared to the size of the cylindrical container. In typical experiments, the number of lattice sites is on the order of 100 in a linear direction, thus vortex solutions of the NLDE may be described accurately. In this article we study relativistic quantum vortices in the superfluid phase of a Bose gas at the Dirac point of a honeycomb optical lattice focusing on explicit vortex solutions of the NLDE. Our vortices are solutions of equations reminiscent of relativistic systems, but pertain here to ultracold atomic gases, where we work in the mean-field limit throughout. Solution types differ by asymptotic conditions on the amplitude for radius r much greater than the healing length, far from a vortex core, or r much less than the healing length deep inside a vortex core. A general schematic for vortices in our problem is shown in figure 1. Solutions may also differ by use of the internal degree of freedom in the spinor structure and with respect to quantization of rotation. However, we find that spinor components within a particular solution always differ by one unit of rotation.
Quantized rotation in a 2 1 + ( )-dimensional spin-1/2 order parameter is best characterized in terms of the topology of the circulation or phase winding around a singular core. Consider an overall macroscopic phase which wraps around the circular boundary consisting of radial lengths much larger than the healing length. When the internal spin degrees of freedom remain topologically unconstrained we simply referred to this as a vortex. In contrast, when the wrapping around the circular boundary includes the internal spin degrees of freedom in a nontrivial way the order parameter is usually called a texture or skyrmion [11,42]. Throughout our work we adhere to the case where the Berry phase is spliced into the overall U(1) symmetry upon condensation of the order parameter. Thus, we do not consider solutions with a macroscopic π-phase disinclination. We categorized our solutions according to their topological properties based on symmetries of the order parameter manifold. These symmetries must take into account the internal spin degrees of freedom as well as the overall U (1) phase. In particular, localized solutions may asymptotically approach either a local minima, maxima, or saddle point of the effective potential in the mean-field energy functional. The first leads to topological categorization, whereas the latter two cases are not topologically protected configurations.
The presence of the Dirac operator in the NLDE adds new features over conventional spinor BECs. In some cases the derivative terms act to defocus the solution leading to conventional 'dark' vortices with density minima in an otherwise uniform background. In other cases the cross-mixing of spinor components in the Dirac term has a confining effect resulting in a 'bright' vortex made of spinning localized density rings over a zero density background. These localized solutions are similar to ones found in ordinary BECs with attractive interactions [43,44], although we emphasize that in the present work we consider only repulsive contact interactions. It is significant that our solutions apply to certain spin-orbit coupled BECs [45]. In the low-energy limit of these theories the linear (Dirac) kinetic term is dominant and one may safely neglect the second order (Schrödinger) contribution. Tuning the interactions to eliminate cross terms in the spinor components leads to our NLDE. The various types of NLDE vortices are summarized in figure 2.
To put relativistic vortices in a more general context, the large symmetry groups encountered in high energy physics and cosmology admit vortex solutions as fundamental ingredients. Examples include cosmic superstrings [46][47][48], mirror symmetry of Calabi-Yau manifolds [49,50], and brane world scenarios [51]. Vortices are generic to gauge theories where spontaneous symmetry breaking occurs. This includes both the Abelian case [48,52], as with Nielsen-Olesen vortices in the Abelian Higgs model, in addition to the non-Abelian case [53], non-Abelian Yang-Mills for example. Supersymmetric field theories which exhibit weak- Figure 1. Relativistic vortex in the honeycomb lattice. The vortex current is depicted by red arrows with the core towards the interior of the red circle. The relativistic current is quantized around the core and incorporates the overall U(1) phase into the internal SU(2) rotation. strong duality rely fundamentally on the presence of solitons and vortices, as these provide the necessary degrees of freedom to make the dualities possible [48,54].
Our results are organized as follows. In section 2, we express the full time-dependent NLDE in (2+1)D in radial form appropriate for stationary vortices. This reduces the problem to the solution of two coupled first order nonlinear ordinary differential equations. In section 3, symmetries of the order parameter are determined in order to classify the various types of solutions. Section 4 develops the Lagrangian formulation of the nonlinear Dirac equation and connects to the standard theory of spontaneous symmetry breaking. In section 5, we use the mean-field energy functional to examine in detail the interplay between the effective potential and Dirac kinetic terms to gain deeper insight into quasi relativistic vortices. In section 6, we solve the NLDE analytically using asymptotic Bessel and algebraic methods. Density and phase plots depicting each solution type are provided. Section 7 connects our solutions to vortices in the NLSE with correction terms, through a semiclassical reduction. In section 8 we present numerical methods and solutions for vortices. This approach allows us to explore the full vortex landscape by tuning the various physical parameters using the method of numerical shooting. Section 9 presents vortex solutions and spectra in a weak harmonic trapping potential building on results from the previous sections. Finally, in section 10 we conclude.

Vortices in the NLDE
In this section we express the NLDE in the appropriate coordinates for finding vortex solutions and analyze the solution space at the energy functional level. We show that in different regions of solution space the effective potentials switch between focusing and defocusing forms, providing conditions for either bright (localized ring) or dark (conventional) vortices. The NLDE describes the dynamics of a four-spinor in 2 1 ) with the upper (+) and lower (−) two-spinors relating to opposite K and K ¢ points of the honeycomb lattice and : .
The full NLDE can be expressed as Figure 2. Schematic of relativistic nonlinear Dirac vortices. At the Dirac point (horizontal direction) changing the interaction, lattice constant, particle hopping, and density tunes the order parameter between bright vortices (red dots over yellow) with different quantized rotation .
The peaks and valleys denote topological transitions between elements of the fundamental group S 1 1 p ( ) around the vortex core. Changing the chemical potential away from the Dirac point (vertical direction) transitions the order parameter into skyrmion solutions with a single quantum of rotation 1. = ℓ Higher energy then leads to multiple phase winding states 2  ℓ (red dots over blue), the energy increasing with winding . ℓ =w here a and t h are the lattice spacing and hopping energy, and L z , g, n are the vertical oscillator length, two-body interaction, and average three-dimensional atomic density, respectively. The associated long-wavelength healing length associated with the NLDE is then found to be t a U 3 2 .
h Dirac x = For a complete list of the relevant renormalized physical parameters with experimental values see [14].
We seek planar stationary solutions to equation (3), and thus we require the upper two-spinor to be expressed in factorized form t t r r r , exp i , , where μ is the chemical potential of the system. The spinor components r A B y ( ) ( ) here are complex functions on the plane, i.e., : .

For stationary solutions in rectangular coordinates, equation (3) further reduces to
with the full solution expressed as a linear combination of solutions from Dirac points K and K . ¢ To obtain cylindrically symmetric solutions with arbitrary integer phase winding, i.e. vortices, we transform to plane-polar coordinates. The time-independent spinor wavefunctions are then written as where ℓ is the angular momentum quantum number (integral phase winding) and the radial functions f A B ( ) are real functions of r. The asymmetric dependence on the winding number ℓ between spinor components in equation (7) is required. In polar coordinates, the derivatives in equations (5), (6) acquire factors of i exp q -( ) and exp i , q ( ) respectively. An offset factor of exp iq -( ) in A y is then required in order to separate the rotational and radial dependence. Applying the method of separation of variables forces the particular form of angular dependence in equation (7). Equations (5), (6) then reduce to the two coupled radial equations Equations (8), (9) comprise the key system of coupled nonlinear equations of this article. In the sections that follow we will solve these using combinations of analytical and numerical techniques for various asymptotic and core boundary conditions. Localized solutions of equations (8), (9) may be categorized by their asymptotic forms, i.e., whether the amplitudes f A B ( ) have zero or nonzero limits far from the core, a relevant distinction in our problem. In particular, we note that in conventional BECs with repulsive interactions the nonlinearity is always defocusing which leads to a vortex with the familiar constant ambient profile that vanishes only at the core. This contrasts radically to the physics in our problem, which we will see allows for (bipartite) vortices with nonzero asymptotic profiles in addition to ones with asymptotically vanishing profiles. The latter structures are the bright vortices or ring-vortices previously mentioned.

Symmetries of the order parameter
To gain a foundational understanding of vortex solutions of the NLDE, we study first the symmetries of the order parameter. Condensed bosons near the Dirac point of a honeycomb lattice are characterized by U(1) symmetry breaking in the full many body theory, and transfer from the ground state to the metastable state at the edge of the Brillouin zone edge [14]. The order parameter is comprised of 2 complex amplitudes , A B T y y y = ( ) associated with the bipartite lattice, forming a representation of the two-dimensional special unitary group The Spin group in 2D forms a double cover of the rotation group SO 2 U 1 . @ ( ) ( ) Spin components A y and B y are coupled through the kinetic part of the single particle Hamiltonian but are independent with respect to onsite interactions. Thus, the full symmetry group G for a Dirac-point BEC that leaves the mean-field energy functional invariant is where the subscripts f and S denote the gauge and spin degrees of freedom, respectively. The 2D representation of G is generated by the 2 × 2 unit and Pauli matrices , , , x y s s  whose Lie group is given by a Î θ is the polar angle, and f is the gauge degree of freedom. In matrix form, G can be parameterized as where again f is the gauge angle, θ is the polar angle, and where the discrete cyclic (off-diagonal matrix) part of G is isolated by setting 0. f q = = One sees that an arbitrary order parameter ψ, obtained by applying G to a representative spinor , 0 y has an additional discrete symmetry H. This is evident by taking f f p  + and 2 q q p  + in equation (11) described by the discrete symmetry , For the gapless system considered here, we generally require both sublattices to be occupied, so that a representative standard spinor is given by 1,1 .
Requiring the order parameter to be single valued under spatial rotations places a constraint on the gauge angle.
In particular, we require f to contain a part that depends on the polar angle by defining 2 f f q ¢ = + and taking f ¢ as the new gauge parameter. Note that we could have factored equation (13) another way and obtained e e 1 e e 1 , 1 4 The discrete choice of redefinition 2 f f q ¢ =  splits the order parameter space into right and left chiral spinor manifolds displayed in equation (13) and (14), respectively. Imposing the condition of single valuedeness for the order parameter has effectively eliminated the extra discrete symmetry .

Lagrangian formulation and spontaneous symmetry breaking
So far we have focused primarily on the symmetries of the order parameter associated with the Dirac point structure. We now examine the effect of the nonlinearity associated with contact atomic interactions. Before we obtain exact vortex solutions to the nonlinear Dirac equation it is helpful to study the theory at the Lagrangian level in order to gain insight into the interplay of kinetic and interaction terms and to acquire a qualitative sense of the kinds of vortex solutions one should expect. From a field theory viewpoint, equation ( Note that a factor of the two-dimensional average density n 2D appears in equation (15), which gives  the correct units of energy (see [56]). The positive sign on the second term is correct, since we are taking derivatives with respect to the conjugate of the field and its derivatives. The more common approach is to differentiate the Lagrange density with respect to the field and its derivatives and take the conjugate of the resulting equations afterwards. This leads to the same result. Equation (15) describes the dynamics of two self-interacting, scalar fields coupled through the spatial derivative terms, with interaction strength U. The corresponding Hamiltonian is derived using the formula where A p and B p are the canonical momenta associated with A y and B y respectively, defined in the usual way by The Lagrangian density in plane polar coordinates reads n n c r r U Since we consider only the case of repulsive interaction where U 0, > and given the linear derivatives in the Lagrangian density, the condensate can lower its energy by finding favorable combinations of fields and gradients that render the derivative terms negative. This leads to the formation of a relativistic vortex on an otherwise uniform positive energy background. Spontaneous symmetry breaking is exhibited in equation (19) by expressing the field energy as the sum of a uniform contribution equal to the chemical potential of the system μ, plus an additional piece that accounts for the presence of a vortex. Substituting the ansatz r t v r t , e , (19), where the factor of 2 accounts for an even sublattice asymptotic particle distribution, we obtain Reading off the effective potential energy density in equation (20), Decomposing the sublattice order parameters in terms of amplitude and phase v e , ( ) we see that the minima of the effective potential energy are located at four points in spin space with infinite U(1) degeneracy at each point: Thus, we see that topological solutions are constrained beyond the symmetry classifications presented so far. At ultra low temperatures the nonlinearity arising from local onsite atomic interactions locks the order parameter into a one of four possible spin configurations characterized by different relative internal phases. However, taking onto account the phases , f ( ) the symmetry group of the NLDE nonlinearity reduces to where the discrete transformation simply exchanges v A and v B in equation (21). So far, our analysis shows that topological solutions of the NLDE are possible and should involve mapping of one or both U 1 A B ( ) ( ) phases to the circle S 1 at spatial infinity. When both phases are involved, the vortex core describes ) We note that non topological solutions are also possible. These are solutions that extremize equation (15), but asymptotically approach saddle points or local maxima of equation (21). Such solutions are possible even in the presence of strictly repulsive interactions. This is because of the rich structure contained in the spin-orbit type (kinetic) coupling in equation (15), the effects of which we shall examine in further detail.

Energy functional analysis for vortices
To study the interplay between the effective potential and kinetic terms of the NLDE, we would like to map out the energy landscape for a vortex configuration in more detail. Here, we work with the mean-field energy functional for a generic vortex by eliminating the angular dependence and incorporating centripetal terms into the effective potential energy. Assuming the vortex form in equation (7) and taking 0 = ℓ in our analysis, the total energy is given by where, and we regulate the energy by introducing an upper cutoff R. Note that radial derivatives are indicated by prime notation. Evidently the f f In the f f , A B -plane, the centripetal term in V eff dominates for small r and forms a saddle point at the origin (of the f f , A B -plane) becoming singular when r = 0. Because of the saddle-point, points on the f A and f B axes have zero potential energy but are unstable. If we want to study solutions that begin at the saddle-point, i.e., f f 0 0 0, or on the f B -axis, say, near the saddle-point, we must include the full contribution from the kinetic terms. In fact, we see that a path along the f B -axis defined by p t t e : ,

) This can be made arbitrarily large and negative by
adjusting the value of f A ¢ at r = 0. The V eff -landscape flattens rapidly as r increases from zero, so that it may be possible to attempt a solution for which f 0 0, If we consider finite energy solutions at large r, f A ¢ and f B ¢ must go to zero and the centripetal term may be neglected. In this limit, we find that V eff has nine local extrema, four of which are minima for the total energy E and thus associated with topological solutions. We obtain the following asymptotic (large r) critical points for )-local maximum. For our chosen case study with 0, = ℓ i.e., with vorticity in only of spinor component, we then expect topological solutions with properties f 0 0, ) This behavior describes a Mermin-Ho vortex. Note that more general topological solutions exist for arbitrary winding with rotation in both spinor components which must both vanish at r = 0, in order to maintain the finite energy requirement. Moreover, we find that several nontopological localized solutions exist as well, with asymptotic limits coinciding with local extrema that do not minimize the potential energy. Such solutions may satisfy the conditions f 0 0, ) which describes a ring-vortex/soliton. Another possible nontopological solutions occurs for the same initial conditions but with asymptotic behavior ) which describes the vortex/soliton solution. The appearance of non-topological solutions that are nevertheless localized is a consequence of the various possible combinations of configurations for the spinor components at the core of a vortex. Some of these result in attractive contribution to the total energy leading to self-trapping of the order parameter, which is particularly dramatic in the case of ring-vortex solutions. A look at the energy functional for equations (8), (9) will help us gain insight into this effect. The NLDE energy functional associated with a vortex with arbitrary circulation ℓ is Taking a different viewpoint, we can define the component-specific effective potentials as where we have introduced condensed notation for the angular momentum and atom-atom interactions U U , , a a ℓ respectively. The energy functional equation (25) reduces to the form The effective potentials U A B eff ( ) encapsulate the angular momentum, binary interactions, and derivativeamplitude bilinear terms. These effective potentials exhibit transitions between focusing and defocusing forms for different regions of spinor solution space. A convenient way to characterize such transitions is through the eigenvalue of the 2D massless radial Dirac operator defined as ) we can further divide the spectrum into the negative, null, and positive subspaces , 0 , The nature of the effective potentials U A eff and U B eff (and thus vortex solutions) is linked to the character of the subspace of D ℓ that dominates a particular region of the NLDE solution space. Table 1 displays the various regions based on analysis of the effective potentials in equations (26), (27).
The conditions in table 1 are not exhaustive but provide some evidence for the simultaneous presence of dark and bright vortices in the same physical system, particularly in the presence of repulsive atomic interactions. Upon obtaining explicit solutions, we will see that condition 1 and 2 describe the upswing and downswing behavior near the apex of the bright vortex and excitations of the asymptotically flat vortex. This reversal in trend for the radial spinor profiles characterizes the self-trappin effect of focusing potentials. Condition 3 applies to all texture solutions whose spinor degrees of freedom are topologically constrained. Condition 4 gives rise to conventional dark vortices with arbitrary phase windings and non-vanishing tails but lacking the topological constraint of the skyrmion.

Analytical vortex solutions
In this section we derive explicit vortex solutions of the NLDE using a combination of analytical and numerical techniques. Some of the solutions have been presented in our previous work [14] but without the mathematical detail in the present article. Different methods are applicable depending on whether the chemical potential is zero or finite-valued and if the vortex has a single unit of winding or arbitrarily large winding .
ℓ For arbitrary ¢ ( ) i.e., the largest contribution to a vortex solution may come from the positive (pos.) or negative (neg.) part of the spectrum or may be an equal admixture of the two (null). Note that conditions 1-5 do not exhaust all possibilities but are listed to demonstrate the variability needed to observe bright and dark vortices.
winding and nonzero chemical potential we use both an approximate asymptotic method and numerical shooting. In the case where 0 = ℓ or 1 exact analytical solutions are possible, and in particular when 0 m = concise algebraic forms occur. For all of the solutions in this section we consider only the case of a spatially infinite condensate. In experiments, this is a good approximation when the healing length, and thus the core of the vortex, is small compared to the oscillator length associated with a trapping potential which sets the size of the BEC.
6.1. Asymptotic bessel solutions for large ℓ Returning to the cylindrically symmetric form of the NLDE, equations (8), (9), we can see that for winding values 2  ℓ both f A and f B must vanish at r = 0 due to the presence of the centrifugal terms. Thus, we will treat the special cases 0, 1 = ℓ in a separate section. To obtain vortex solutions we require spatial derivatives to vanish at infinity; this means that for r ,  ¥ equations (8), (9) yield the asymptotic behavior The important point to note here is whether the f A (B) tends to zero or a non-zero constant as r .  ¥ To determine asymptotic forms for the radial functions we first look for exponential decay to zero, or growth or decay to the constants U , m  in which case we find no consistent asymptotic solution to equations (8), (9). Algebraic decay to zero may be discerned by substituting f r r For the case of non-vanishing boundary conditions and 2  ℓ the asymmetry in the angular momentum terms complicates the task of finding solutions. We recall that the factors r r ¶ + ℓ and r 1 r ¶ +ℓ ( ) act as index raising and lowering operators for the Bessel functions J n . One can verify that Bessel functions are nearly exact solutions in the case of a weak nonlinearity given by U n c 1.
where A and B are normalization constants, J n (r) is the Bessel function of the first kind of order n, a 0 and a n are the expansion coefficients. Anticipating an asymptotic expansion, we include the real function Q(r) in the complex prefactor in order to cancel oscillations in the Bessel series at large argument. In addition, we have included the radial function F r , ( ) which we require to cancel the r 1 asymptotic decay of the Bessel functions.
The series that we have chosen to use runs over the Bessel index rather than the usual form where the summation runs over the zeros of a single Bessel function with a fixed index. This choice of expansion is valid but does not offer the convenience of using the standard orthonormal relations for Bessel functions when computing the coefficients a n . Substituting the ansatz equations (32), (33) into (8), (9), we then consolidate the angular momentum and derivative terms in the resulting series expansion by using the recurrence relations for Bessel functions: C a a n n a n n 1 2 C a a n n a n n 1 2 In equations (34), (35) we have absorbed a fraction C 1 1 -< from the chemical potential terms into the recursion relations. Solving the recursion relations leads to where after integration the oscillating part of equation (40) The overall constant A is determined by normalizing the wavefunction to the particle number. We emphasize here that equation (43) is an approximate solution which applies only very near to or far from the vortex core, and for large winding number. Density and phase plots for the complex Bessel solution are shown in figure 3 on the left with radial profile plots shown in figure 4(a).

Algebraic solutions
Here we treat solutions for 1 > ℓ at zero chemical potential 0. m = For zero chemical potential, algebraic closed forms exist for arbitrary values of the winding number. In contrast, for general values of μ closed forms exist only for 0, 1, = ℓ which we discuss in the last paragraph of this section. By setting 0 m = in equations (8), (9) we obtained exact non-topological vortex solutions. From a technical standpoint eliminating the chemical potential terms in the NLDE simplifies the problem considerably. This is because solutions to the homogeneous noninteracting equations are algebraic forms, not Bessel functions. Thus, vortex solutions do not connect to Bessel functions in the zero-interaction limit. We note here that the precise relationship between the chemical potential μ and interaction strength U is determined by computing the spectrum of each solution confined within a harmonic trap. This is the subject of section 9 of the present article. To obtain algebraic solutions we start from an ansatz f Ar r 1    Equations (47), (48) describe a vortex whose density peaks in the shape of a ring with a bright soliton (no rotation) located at its center. This solution corresponds to the ring-vortex/soliton from our previous work [13], but obtained here as a special case of the more general result in equations (44), (45). Note that Setting 0 = ℓ in equations (44), (45) effectively interchanges the forms for f A and f B but leads to the same solution type. The ring-vortex/soliton radial profiles are plotted in figure 4(d) with the density and phase shown on the right of figure 5.
Addressing now the case 0 m ¹ and 1, = ℓ an algebraic solution is obtained using an ansatz similar to that used to obtain equations (44), (45), which gives

Skyrmion solutions
To obtain Skyrmion solutions we choose an ansatz of the form f cos , where the parameters η and j are functions of the radial coordinate. For background on skyrmions in two-dimensions see [42,57]. Substituting these forms into equations The first point to note is that the centripetal terms place a restriction on the behavior of j for r 0.  for the constant envelope factor. The radial profiles for both the Mermin-Ho and Anderson-Toulouse vortices are obtained by solving equations (54), (55) using a straightforward shooting method [43]. We have plotted the radial solutions for both types of skyrmions in figures 4(d) and (g), with plots of the density and phase shown in figure 6 (right side) and figure 7 (left side).

Half-quantum vortices
The NLDE supports half-quantum vortices, i.e., vortices with half-integer winding 1 2. = ℓ Such solutions are obtained by forming superpositions of the components f A and f B which solve the Mermin-Ho condition [61,62]. The necessary requirement is that both spinor components approach the same nonzero value at spatial infinity far from the vortex core (see figure 4(g)). The appropriate spinor combinations read From equation (58) we compute the geometric phase by circling the vortex core through the Berry phase prescription exp d , 59

Vortices in the semiclassical reduction to NLSE with correction terms
We now study the connection of NLDE vortices to their nonlinear Schrödinger counterpart. Working again from the ansatz t v t r r , e , , where we have factored out the dominant energy contribution from the spinor functions, and substituting into the NLDE gives where we have abbreviated the differential operators by using i .
x y  = ¶ + ¶ For clarity of notation, we have set c 1 l  = = and used the abbreviated subscript notation to indicate differentiation. We then make the following approximations Note that equality from this point on must be understood within approximations equations (63)- (65).
is the NLSE with derivative correction terms and attractive constant background potential . m -Note that the effective particle mass is 2.
m The correction terms in equation (70) are smaller than the first three terms on the right hand side due to the extra factor of U m and powers of derivatives of v A . This tells us that localized solutions will differ from standard NLS form only in regions where gradients are steep, i.e., near the core of a vortex.
Next, we look for vortex solutions of equation ( Table 2. Vortex solutions of the NLDE. Solutions are described by their phase winding, closed-form expression, and asymptotic properties. The Mermin-Ho, half-quantum vortex, and general topological solutions have conserved topological charge associated with azimuthal phase winding classified according to elements of the fundamental group S .
In all other cases, the spinor components asymptotically approach either a saddle point or a local maximum of the effective potential. Note that r 0 is the length scale associated with the chemical potential or the interaction strength depending on the particular solution.
Inserting equation (71) gives Canceling some terms and using condensed notation, we have There are several limits of this equation that are interesting: (i) For small r we expect that v 0.  In this limit the nonlinear terms on the right hand side become negligible This is Bessel's equation and the solutions are the well known Bessel functions J n (x) and Y x . n ( ) The ones we are interested in are the Bessel functions of the first kind, J x , n ( ) since these are regular at the origin. Of these, we further restrict solutions to those for which n 0 > since these are zero at the origin.
(ii) We seek solutions with constant square modulus for large r. In this limit, all derivatives may be set to zero in equation (78). This allows us to solve for C in terms of the constants U, μ, and ω. For large ξ (large r), v 1,  and equation (78) reduces to For U 0,  we should retrieve the standard Schrödinger result for the radial part of the wavefunction. This is indeed the case since, in this limit, all terms on the right hand side of equation (78) vanish and we are left with Bessel's equation as expected.

Vortices by the method of numerical shooting
In general, vortex profiles for arbitrary phase winding ℓ and chemical potential μ may be obtained using a numerical shooting method as described in [43] and [63]. These include topological solutions whose tails do not vanish asymptotically, as well as non-topological vortices where one or both spinor functions have asymptotically vanishing profiles. To proceed we first convert the NLDE to dimensionless form by introducing the rescaled radial coordinate and spinor components Here the dependence of the solution on the choice of angular quantum number ℓ is implied. The asymptotic form of these equations shows that moduli for convergent solutions A B h | | ( ) can approach either 0 or 1 for large χ.
To study the analytic structure as well as a practical starting point for application of numerical methods, we expand the solution in a power series in χ  a  a  3  3  3 , 86 Choosing a particular value for ℓ determines the values of both a 0 and b 0 and we find that there is only one independent parameter. Thus, for a given value of ℓ a vortex is found by tuning a 1 ℓ towards a critical value a . The associated radial profiles for these solutions are plotted in figure 8. In figure 9 we illustrate convergence for both spinor components for the 2 = ℓ vortex. Convergence to the vortex profile is seen by overlaying an excitation of the vortex (a) and the free-particle Bessel-like solution (b). Note that there are two types of radially excited solutions characterized by the property that A B    figure 10(c). This is what we expect to see for weak nonlinearity. This illustration highlights the property of the Dirac kinetic terms which act by confining the solution in the presence of strong nonlinearity. So far in this section we have shown how tuning the initial condition above some critical value forces the solution into the strongly nonlinear regime, whereas below this value leads to solutions in the weakly nonlinear regime. In our analysis we held the chemical potential and interaction fixed, U 1. m = = We now study the effect of tuning the chemical potential for 1 0 , m m >  = while maintaining the initial value in the shooting process and the interaction fixed. For arbitrary phase winding we will see a progression from Bessel solutions through a topological vortex, then finally ending in a single bright ring-vortex with vanishing tail. This solution resembles those found in attractive BECs and also in optics. Although we demonstrate such solutions by analytical methods in the next section, the numerical approach allows us to demonstrate the transition from the free-particle limit, U m  (weak nonlinearity), to the vortex limit, U m  (strong nonlinearity). To show this transition we choose a particular rescaling of the NLDE such that the limit 0 m  makes sense, i.e., The resulting dimensionless NLDE is Starting with the case 1, Here k is a parameter tuned to give the desired values of functions at the origin, analogous to a 0 for the forward shooting. The Anderson-Toulouse vortex is obtained similarly but with the boundary condition cos 2 10 , The half-quantum vortex is obtained by forming linear combinations of the numerical Mermin-Ho components, as previously discussed.

Discrete spectra in a harmonic trap
We now extend our numerical studies from section 8 to the treatment of vortices bound within a weak harmonic potential similar to the analysis in our previous work on solitons (a brief summery of these results can be found in [14], presented here in full detail). Physically this accounts for an additional external harmonic trapping potential present in real experiments. We follow a similar procedure as before in preparing the NLDE for numerical analysis by first converting to a dimensionless form. Here we examine the same solutions as in section 8 but with tight confinement assumed in one of three spatial dimensions and with additional weak harmonic confinement in the remaining two directions. The oscillator frequencies thus satisfy The sequence of plots depicts the transition from the weakly-interacting Bessel-like solution in (a) for which U , m > to the strongly nonlinear case for which U m < in (i). Note that 1 m = is the boundary between solutions which oscillate around 0 and solutions which oscillate around 1. The results in these plots show that general ring-vortex solutions may be obtained by starting from excited states such as in (e), then reducing m towards zero. for all of our calculations. The spectrum for a fixed quantized mode tracks the flow of the chemical potential as one tunes the system between the free-particle and strongly nonlinear limits. In the former limit radial profiles for the lowest excitations resemble decaying Bessel functions as in harmonically trapped massless Dirac particles, whereas in the latter limit the nonlinearity dominates and the lowest trapped solutions flatten, consistent with the Thomas-Fermi regime. For higher quantized modes radial derivatives force the solution to maintain the Bessel-like form even as one tunes the parameters into the strongly nonlinear regime. Convergence of our solutions is provided in appendix.

Conclusion
In this article we have solved the NLDE by a variety of different methods for both the idealized case of zero trapping potential and in the harmonic case. The latter context provides the opportunity to study spatial quantization of vortex solutions in the radial dimension. A preliminary asymptotic solution using Bessel function expansions provides insight into the structure of the NLDE itself. Algebraic solutions were then obtained by considering the case of zero chemical potential. For these solutions the derivative and nonlinear terms are perfectly balanced leading to bright vortex rings over a zero-density background, a direct consequence of the Dirac operator. More generally, for nonzero chemical potential, we found that analytical solutions are only possible when one unit of winding is considered. In this case we used a numerical approach to obtain vortices with arbitrarily large winding number. A combination of numerical and analytical techniques yields skyrmion and half-quantum vortices, i.e., textures. Table 3. Numbers for computing spectra of 2 = ℓ topological vortex ground state solution. For a fixed value of the chemical potential , m the free parameter a 0 is tuned until the desired radially quantized state is reached, for which one then computes the normalization  by equation (100). As m increases into the Thomas-Fermi regime the dependence of the solution on a 0 becomes more sensitive, requiring a greater degree of tuned accuracy as shown in the column on the left.  Having obtained our solutions, we computed their discrete spectra in the presence of a weak harmonic potential. This gives us the low-temperature μ versus U landscape for relativistic vortices, where μ and U are the chemical potential and lattice renormalized particle interaction, respectively. For example, for finite μ we found that a series of phase transitions occur as U is tuned from zero upward: we encounter a Mermin-Ho skyrmion transitioning into a half-quantum vortex, followed by the Anderson-Toulouse skyrmion, then finally into a vortex/soliton, i.e., a bright soliton at the core of a singly-wound vortex. Some of our solutions are similar to those obtained in spinor BECs, skyrmions in particular. However, presence of the Dirac operator contrasts heavily with the Laplacian case. This difference is most obvious in our bright ring-vortex solutions. Components of these solutions resemble the bright vortex which occurs in attractive BECs, but in our case the confining (focusing) regime of the Dirac operator is responsible for the effect in spite of repulsive atomic interactions. Furthermore, we should emphasize the significant distinction between our results and similar confinement effects in spinor BECs and in traditional models such as Thirring and Gross-Neveu. In each of these other theories the presence of confinement relies on attractive interactions in addition to a mass term, whereas we consider strictly the repulsive case with zero mass.
The interdisciplinary nature of our work suggests several future research directions. For instance, the form of our solutions implies possible mappings to finite energy solutions to classical gauge field equations. In particular, there is a deep connection between the NLDE and Chern-Simons terms relevant to quantum Hall fluids and more generally to relativistic field theories [64][65][66]. To see this connection consider that equation (3) resembles the low-energy effective theory for massless Dirac fermions interacting through a gauge field, where the latter has been absorbed into the local fermion contact terms. The Dirac term by itself is symmetric under global phase and chiral rotations which are made local through the addition of a gauge field. Quantization results in the well known axial anomaly which one finds to have exactly the Chern-Simons structure [67][68][69]. Thus, quantized theories of interacting fermions generally result in couplings to Chern-Simons terms. In our case we have focused only on the fermion part of the argument which naturally retains imprints of the omitted Chern-Simons terms. Another argument for the NLDE/Chern-Simons connection hinges on the well established duality between the Thirring model in (2+1)-dimensions and Maxwell-Chern-Simons theory. The mapping is arrived at through bosonization, i.e., reformulation in terms of paired particle and antiparticle field operators effective at strong coupling or low energy. These points merit further inquiry into potentially fruitful questions. Finally, in light of the successful analogs to date connecting condensed matter and nonlinear optics we expect that our results should be reproducible within an optics setting. Another obvious direction is to go deeper into the mathematics of nonlinear partial differential equations. The vast amount of work in this area provides an array of solution and classification methods which may be used to fully understand critical bounds for wellposedness as well as general solutions to the full NLDE with time dependence.

Acknowledgments
This material is based in part upon work supported by the National Science Foundation under grant numbers PHY-1067973, PHY-1011156, and the Air Force Office of Scientific Research grant number FA9550-08-1-0069. LDC thanks the Alexander von Humboldt foundation and the Heidelberg Center for Quantum Dynamics for additional support.

Appendix. Convergence of numerical vortex solutions
We demonstrate convergence of solutions in the harmonic trap for three of the 2 = ℓ topological vortices associated with the black curve in figure 14 These values interpolate between the free-particle and strongly nonlinear limits (small to large μ values). These solutions were obtained by finite differencing using a shooting method to tune the precision of the initial value of A y such that