The nonlinear Dirac equation in Bose-Einstein condensates: Superfluid fluctuations and emergent theories from relativistic linear stability equations

We present the theoretical and mathematical foundations of stability analysis for a Bose-Einstein condensate (BEC) at Dirac points of a honeycomb optical lattice. The combination of s-wave scattering for bosons and lattice interaction places constraints on the mean-field description, and hence on vortex configurations in the Bloch-envelope function near the Dirac point. A full derivation of the relativistic linear stability equations (RLSE) is presented by two independent methods to ensure veracity of our results. Solutions of the RLSE are used to compute fluctuations and lifetimes of vortex solutions of the nonlinear Dirac equation, which include Mermin-Ho and Anderson-Toulouse skyrmions, with lifetime $\approx 4$ seconds. Beyond vortex stabilities the RLSE provide insight into the character of collective superfluid excitations, which we find to encode several established theories of physics. In particular, the RLSE reduce to the Andreev equations, in the nonrelativistic and semiclassical limits, the Majorana equation, inside vortex cores, and the Dirac-Bogoliubov-de Gennes equations, when nearest-neighbor interactions are included. Furthermore, by tuning a mass gap, relative strengths of various spinor couplings, for the small and large quasiparticle momentum regimes, we obtain weak-strong Bardeen-Cooper-Schrieffer superconductivity, as well as fundamental wave equations such as Schr\"odinger, Dirac, Klein-Gordon, and Bogoliubov-de Gennes equations. Our results apply equally to a strongly spin-orbit coupled BEC in which the Laplacian contribution can be neglected.


Introduction
Two contemporary themes in the study of cold atomic gases are the creation of new exotic forms of quantum matter, and quantum simulations of systems already present in nature [1,2,3]. By tuning the parameters for a collection of atoms and lasers one may address problems in quantum many-body systems or in high-energy physics [4]. In the first case degeneracy, quantum correlation, and entanglement are essential ingredients, whereas the latter case usually focuses on low-energy fluctuations of systems where a macroscopic fraction of particles reside in a single quantum state, often amenable to Landau descriptions. The versatility of Bose-Einstein condensates (BECs) allows the freedom to specify the geometry and topology of the order parameter to suit a particular purpose. For example, spinor BECs provide one way to realize order parameters with large symmetry groups, and hence exotic topologies [5,6,7,8,9,10,11,12,13,14,15]. It follows that the inclusion of spin-orbit coupling in such systems increases their complexity and introduces topological order [16,17,18,19,20,21], a distinct classification for the order parameter. However, in order to access interesting physics and to simulate new regimes it may be necessary to extend beyond the usual notion of stability to metastable non-ground-state or non-equilibrium BECs.
In this article, we develop some of the fundamentals underlying non-ground-state BECs in quasi-two-dimensional (quasi-2D) honeycomb lattices and the associated longwavelength emergent theories [22,23,24,25,26,27,28,29]. We focus in particular on superfluid fluctuations in the presence of Dirac points from a semiclassical perspective and by including lowest-order quantum effects. Quantum fluctuations are determined by solving the partial differential equations which describe dynamics of the low-energy modes for an arbitrary condensate profile. These equations are Lorentz invariant and comprise a relativistic generalization of the Bogoliubov-de Gennes equations (BdGE); thus we call them relativistic linear stability equations (RLSE). The RLSE provide a means of calculating vortex stabilities, yet their versatility extends beyond stability calculations to simulating a large number of established theories in addition to some exotic ones. This is because quasiparticles in BECs with inherent relativistic structure (e.g., linear dispersion, CPT invariance, multicomponent order parameter, etc.) can be tuned to have linear or quadratic dispersion with a zero or finite gap coupled to a condensate reservoir with a large number of possible internal symmetries. † Moreover, the "no-node" theorem originally proposed by Feynman [34], which constrains conventional BECs, is circumvented for non-ground state (metastable) systems and in the case of son-orbit coupling, as the order parameter in these systems is generally not positive-definite [35,36]. This property is a fundamental feature of quasirelativistic condensed matter systems. In particular, lifting the "no-node" theorem restriction leads to time-reversal symmetry breaking, which allows for exotic bosonic systems such as p-wave superfluids [37], chiral Bose liquids [38], complex unconventional BECs in high orbital bands of optical lattices [39], and BECs with repulsive interactions that support bright solitons and vortices as well as skyrmions [40,41,17]. We point out that our system is identical to a quasi-2D BEC with spin-orbit coupling in either the long-wavelength limit or the strong tunable spin-orbit coupled limit, provided the interactions are also chosen to retain only the intra-component terms. To map to the strong spin-orbit coupled limit, however, the strength of the spin-orbit coupling term must be much larger than the quadratic term but still below the quantum critical point separating the spin-balanced and spin-polarized ground states [42].
Our results focus on three main topics. First, the physical parameters and necessary constraints to construct a non-ground-state condensate at Dirac points are explained in detail. The BEC is tightly confined in one direction and loosely confined in the other two directions. More precisely stated, magnetic trapping along the z-direction is such that excitations along this direction have much higher energy, by at least an order of magnitude, compared to the lowest excitations in the x and y-directions. Thus, an important step is to calculate the precise renormalization of all relevant physical parameters when transitioning from the standard 3D BEC to a quasi-2D system. In addition to this step we also account for renormalization due to the presence of the optical lattice potential which introduces an additional length scale from the lattice constant. We point out that microscopically the BEC obeys the threedimensional nonlinear Schrödinger equation and we consider temperatures well below the BKT transition energy associated with two-dimensional systems. Nevertheless, throughout our work we often use "2D" for brevity, keeping in mind the quasi-2D picture. Condensation at Dirac points of the honeycomb lattice requires additional techniques beyond ordinary condensation, which we have detailed in our previous work [43]. In addition to the fields needed to construct the lattice one requires a resonant field which provides the time-dependent potential to "walk" atoms from the ground state (zero crystal momentum) to the Dirac point. The result is a transient configuration since a macroscopically occupied nonzero Bloch mode is not in thermodynamic equilibrium.
Care must be taken when transferring atoms from the ground state to a Dirac point in order to minimize depletion out of the condensate. In general, one might expect some dissipation to occur due to secondary interactions within the condensate, and between condensed atoms and the lattice, quantum fluctuations, and thermal excitations, the latter two comprising the normal fluid. However, at the mean-field level repulsive atomic interactions within the condensate itself produce a single Hartree term which just shifts the total energy upward without causing additional depletion. Lattice effects are accounted for completely through the band dispersion, since we are not considering the presence of disorder or artificial impurities. Moreover, we consider only the zerotemperature case. There is certainly finite leakage into energetic modes lower as well as higher than the condensate energy. However, such losses can only occur in the presence of higher-order dissipative terms in the Hamiltonian. In this article we restrict our analysis to the effects of first-order quantum corrections and apply our results to the special case of vortex background.
The second major topic in this article addresses linear stability of vortices near a Dirac point. We first provide a detailed derivation of the RLSE then solve them for vortex solutions of the nonlinear Dirac equation [41]. The resulting eigenvalues determine the characteristic lifetimes of each vortex type. Solutions of the RLSE are inherently massless Dirac spinors with components that couple only through the Dirac kinetic terms. For a vortex background, RLSE solutions describe the quantum density and phase fluctuations near the vortex core. Physically, these are local undulations in the density profile, rigid translations of the vortex itself, and fluctuations in the speed of rotation. Although the latter is topologically protected, at the mean-field level quantum effects introduce small admixtures of different winding numbers into the vortex. These admixtures, which take the form of phase fluctuations, comprise the Nambu-Goldstone modes of the system. Near the vortex core they appear as bound states, the lowest of which are zero-energy modes (zero modes): static modes with zero energy associated with spatial translations of the center of the vortex. From a symmetry perspective, zero modes account for the fact that a vortex breaks the translational and rotational symmetry of an otherwise uniform system. We will address the various modes in generality when we discuss the associated reductions of the RLSE to other well known equations. Our work culminates with the connection to several other important areas of physics including relativistic Bardeen-Cooper-Schrieffer (BCS) theory. In addition to continuous space-time dependence, quasiparticle solutions of the RLSE are labeled by two indices associated with the lattice pseudospin valley and the particle-hole structure analogous to Nambu space from BCS theory. In order to avoid confusion we will refer to these as valley and Nambu indices, and reserve the particle-hole terminology to distinguish between the two states in either the valley or Nambu space. The RLSE are formulated to describe excitations in a repulsive Bose gas but can be reinterpreted as excitations in a theory comprised of attractive particles upon pseudospin valley particle-hole exchange. This symmetry is a consequence of the combined symmetry of charge conjugation (C), parity transformation (P), and time reversal (T), which is fundamentally related to the structure of the Dirac operator. Retaining a mass gap, an intermediate step in Dirac-point condensation [43], and adding atomic interactions between nearest-neighbor lattice sites [44] extends the RLSE to the Dirac-Bogoliubovde Gennes equations (DBdGE) [45,46,47,48], provided valley particles and holes are interchanged. This connection is significant as DBdGE are required for a broader description of superconductivity beyond the standard BCS formalism, particularly for superconductors with a high Fermi velocity. Indeed, a relativistic formulation of BCS becomes important for elements with large atomic number (Z ≥ 40), in neutron stars where superfluidity is expected to play a major role in "glitches" [49,50,51], and in color superconductivity where the strong nuclear force provides the attraction between fermions [52]. In the nonrelativistic and semiclassical limits the RLSE reduce to the Andreev equations. These equations were originally formulated to address physics of nonuniform superconductors, for instance a type-I superconductor near a normalsuperfluid interface or a vortex in a type-II superconductor [53,54]. Interestingly, we find that the RLSE reduce to the Majorana equation inside the core of NLDE vortices. From a fundamental standpoint the Majorana equation describes relativistic fermions that are their own anti-particle [55]. In condensed matter systems, and in particular our problem, finite-energy phase fluctuations inside the core of an NLDE vortex connect smoothly to Majorana zero modes. This is significant as Majorana zero modes are presently of great interest in such fields as topological quantum computation [56], topological insulators [57], and more generally in the study of non-Abelian anyons and fractional statistics [58,59,60,61]. Figure 1 provides a schematic overview of some of the theories and physical regimes encapsulated in the RLSE. This article is organized as follows. In Sec. 2, we discuss physical parameters, constraints, and regimes. In Sec. 3, we analyze superfluid excitations for a Bose gas in the honeycomb lattice from a semiclassical perspective. Section 4 contains two derivations of the RLSE according to paths dictated by two possible orderings of the tight-binding and continuum limits of lattice Bloch functions. In Sec. 5, stability analysis is performed for vortex solutions of the nonlinear Dirac equation by solving the RLSE for quasiparticle functions and eigenvalues. In Sec. 6, we examine several reductions of the RLSE to other well known equations. We map the RLSE to the equations for relativistic BCS theory and demonstrate the non-relativistic limit to standard BCS theory. In Sec. 7, we conclude.

Renormalized Parameters and Physical Constraints
To obtain the correct renormalized parameters for the NLDE we proceed by two steps. First, we follow the transformation of the 3D NLSE parameters as we reduce to the 2D NLSE. Second, we take the long-wavelength limit of the 2D theory at the Dirac point to get the NLDE, which induces a second renormalization of the parameters.

Transition from 3D to 2D nonlinear Schrödinger equation
A BEC comprised of N atoms of mass M is described by a wavefunction ψ(r, t) which solves the time-dependent nonlinear Schrödinger equation. The single-particle density is defined as |ψ(r, t)| 2 , the BEC density ρ(r, t) 2 ≡ N |ψ(r, t)| 2 , and the phase is φ ≡ arg[ψ(r, t)], with the superfluid velocity given by v s ≡ ∇φ. The two-particle interaction strength is g = 4π 2 a s /M and the healing length is ξ = 1/ √ 8πna s , where a s is the s-wave scattering length for binary collisions between atoms. We take a s > 0 so that g > 0, i.e., we consider only repulsive interactions, leaving attractive interactions for future studies. Throughout our work, we treat the case of an axisymmetric system associated with a harmonic trapping potential with two large dimensions described by a radius R = x 2 + y 2 , and a small dimension transverse to the plane described by the length L z . The average density which appears in ξ is thenn ≡ N/(πR 2 L z ). Note that ψ(r, t) has dimensions of length −3/2 so that g has dimensions of energy×length 3 . Another important quantity is the speed of sound in the condensate, which is defined as c s = gn/M .
Transforming to the 2D regime requires that a s L z ξ [62,63], which ensures that the condensate remains in the ground state in the transverse direction, and L z R, which ensures that excitations along the plane have much lower energy than those in the transverse direction. The wavefunction can then be separated into longitudinal and transverse modes, following similar arguments as in Ref. [22] ψ(r, t) = ( where f (x, y) and h(z) are the dimensionless spatial functions that describe the longitudinal and transverse normal modes, respectively, and µ is the chemical potential. Projecting onto the ground state of the transverse dimension h gs (z), gives us an effectively 2D wave equation. In the case where L z ∼ ξ, h gs (z) is just the ground state of the one-dimensional particle-in-a-box solution [22], we then have h gs (z) = √ 2 sin(πz/L z ). This reduces the 3D nonlinear Schrödinger equation to the 2D form. It may be convenient to express L z and R in terms of the trap frequencies ω x , ω y , and ω z , in which case we may write L z = ( /M ω z ) 1/2 , R = M −1 (1/ω x + 1/ω y ). The transformation is then completed by defining the renormalized 2D chemical potential and interaction as The 2D renormalized average density can be related to the 3D average density using the transverse oscillator length or frequencȳ Using this definition and the 2D single-particle wavefunction, ψ(x, y) = A −1/2 f (x, y), we can write the 2D condensate density as ρ 2D (x, y) = N |ψ(x, y)| 2 . The 2D renormalized healing length can also be constructed which we find acquires only an extra numerical factor Similarly, we find the 2D speed of sound to be c s2D = g 2Dn2D /M = (3/2) 1/2 c s . It is important to keep track of the effect of the reduced dimensionality on the dimensions of the constants: ψ(x, y) now has dimensions of length −1 , g 2D has dimensions energy×length 2 , andn 2D has dimensions length −2 .

Derivation of nonlinear Dirac equation from 2D nonlinear Schrödinger equation.
The derivation of the nonlinear Dirac equation begins with the second quantized Hamiltonian for a 2D system with the bosonic field operatorsψ ≡ψ(r, t) =ψ(x, y, t) obeying bosonic commutation relations in the Heisenberg picture. We then expand in terms of Bloch states belonging to A or B sites of the honeycomb lattice which breaks up the bosonic field operator into a sum over the two sublattices. The spatial dependence in this expansion is encapsulated in the exponential Bloch wave and the Wannier functions w(x, y) which are then integrated out leaving only number-operator terms in the form of a Dirac-Hubbard Hamiltonian, Eq. (10) in Ref. [27]. Finally, the operator terms are reduced to c-numbers by averaging over on-site coherent states and the long-wavelength limit is taken. We again recover a continuum theory but with a Weyl spinor wavefunction Ψ = (ψ A , ψ B ).
The key point in discerning the correct normalization (and thus other related quantities) is the contraction of the many-body bosonic operators between localized coherent states. The parameter |c i,j | 2 which labels the coherent state at site (i, j), emerges as the number of atoms at each site, so that c i,j itself becomes the continuous amplitude ψ A (r, t) and ψ B (r, t) in the long-wavelength limit. Note that the complex moduli of these amplitudes are pure dimensionless particle numbers, not densities, since they result from taking the spatial integral over the lattice. With the area per lattice site given by A l = √ 3a 2 /4, the local time-dependent sublattice densities can be reconstructed as: ρ A(B) (r, t) = |ψ A(B) (r, t)| 2 /A l . Then, the dimensionally correct sublattice mean-field wavefunctions must be given by where a is the usual lattice spacing. The correct normalization procedure can now be deduced by writing down the total number of particles in the system where the upper limit of the radial integral is taken large enough so that the integrand is negligible. The total number of atoms of the system, N , appears on the left-hand side.
The 3D to quasi-2D reduction and continuum regime result in an effective atomic interaction U , a renormalized version of the usual interaction g. We arrive at the explicit form for U by first approximating the lowest band on-site Wannier functions by the ground state of the harmonic oscillator potential. Integrating over the area of one site, we obtain a new local interaction strength where is the oscillator length of a lattice potential well. It is often more practical to express the area of one site in terms of the lattice constant π 2 = √ 3a 2 /4, and all other parameters in terms of the corresponding 3D parameters. Using Eqs. (2)-(3), the interaction takes the form Note that U has dimensions of energy.
We can now identify the main parameters which appear in the NLDE. The dimensionful coefficient which multiplies the Dirac kinetic term is the effective speed of light c l ≈ 5.31×10 −2 cm/s (compare to the analogous coefficient for relativistic electrons c ≈ 3.00 × 10 8 m/s). In terms of fundamental constants we find c l ≡ t h a √ 3/2 , where a is the lattice constant and t h is the hopping energy. The natural length scale of the NLDE is the Dirac healing length ξ Dirac ≡ c l /U = t h a √ 3/2U , which characterizes the distance over which a disturbance of the condensate will return to its uniform value. We see that ξ Dirac has the correct dimension of length. To simplify the notation, for the remainder of our paper we will omit the 2D subscript on all parameters with typical values as can be achieved in present experiments [43]. Finally, the quantity U which appears in the NLDE determines the strength of the nonlinearity. We have provided a full list of relevant parameters associated with the NLDE in Table 1.

Physical constraints
The realization of the NLDE in a condensate of 87 Rb atoms requires that several constraints are satisfied which we now list and discuss:  (i) Landau Criterion. In order to avoid the instabilities associated with propagation faster than the sound speed in the condensate, we require that the effective speed of light is less than the 2D renormalized speed of sound.
(ii) Long-wavelength Limit. The NLDE describes propagation of the long-wavelength Bloch envelope of a BEC near the Dirac point. Thus, a necessary condition for realizing the NLDE in the laboratory is that the Dirac healing length must be much larger than the lattice constant.
(iii) Relative Lengths for 2D Theory. In order to obtain an effectively 2D system, the vertical oscillator length must be much smaller than the trap size along the direction of the plane of the condensate.
(iv) Relative Energies for 2D Theory. Analogous to the previous restriction, this condition relates to the 2D structure but pertains to the energies of the system. The key point is that we must avoid excitations vertical to the plane of the condensate while enabling them along the plane: the chemical potential and temperature must be less than the lowest transverse excitation energy.
(v) Weakly Interacting Regime. The NLDE and RLSE are derived for a weakly interacting Bose gas. This ensures both the stability of the condensate as well as the effective nonlinear Dirac mean-field description. We then require the interaction energy to be significantly less than the total energy of the system. Yes, we are always in the weakly interacting regime. I added a statement in the overview paragraph in Section 6. In the abstract and introduction I was more specific in descriptions of the limits involved.
(vi) Dirac Cone Approximation. For a condensate in the regime where the NLDE description is valid, we require that the linear approximation to the exact dispersion remain valid. As in the case of graphene, large deviations from the Dirac point induce second order curvature corrections to the dispersion. Thus, we must quantify the parameter restrictions which allow for a quasi-relativistic interpretation. † (vii) Lowest Band Approximation. We derive the NLDE and RLSE assuming that the lowest band is the main contribution to the dispersion.
Having stated each constraint, we can now address each one in detail and explore the conditions under which each is satisfied. In the following, we consider a BEC comprised of 87 Rb atoms where all numbers used are listed in Table 1 and are experimentally realistic [66]. First, the Landau criterion pertains to the effective velocities in the BEC. Stated mathematically, the Landau criterion requires that c l /c s2D < 1. Using the definitions for the effective speed of light and the sound speed found in the first part of this section, we compute c l /c s2D = 0.90, which satisfies the inequality.
The length constraints are as follows. The long-wavelength limit is defined by ξ Dirac /a 1, for which we find that ξ Dirac /a = 6.91. For an effectively 2D system, the required length constraint implies the condition L z R. Taking R ≈ 100 a (a typical condensate size), and using a realistic value for the vertical oscillator length (Table 1), we obtain L z = 2.73 a, which satisfies the constraint. Moreover, we require a healing length close to or less than the transverse oscillator length. With ξ = 1.10 µm and L z = 1.50 µm, we find that this condition holds.
The energy constraints may be stated as µ, k B T ω z . We can solve the NLDE for the lowest excitation to obtain an expression for the chemical potential µ = c l k + U |Ψ| 2 [67]. Next, we evaluate this expression using the lowest excitation in a planar condensate of radius R ≈ 100a, which has wavenumber k ≈ π/2R = 2.86 × 10 4 m −1 . The interaction U is computed using Eq. (7) for the binary interaction g and mass M pertaining to a condensate of 87 Rb atoms. Finally, for a uniform condensate we take |Ψ| 2 ≈ 4/ √ 3 (Eq. (5)) and the constraint on the chemical potential becomes µ = 2.59 nK < 22.17 nK, which is satisfied. For the temperature, we require T ω z /k B . Using the data in Table 1 for the vertical oscillator frequency, we obtain the upper bound for the temperature T 22.17 nK. This is a reasonable requirement given that BEC occurs for T in tens or hundreds of nanoKelvins or as low as picoKelvins.
Next, we examine constraints on the particle interaction. To check that we are in the weakly interacting regime, i.e., that U/µ 1, we use the value for the chemical potential µ which we have just computed and compare this to the interaction energy U , whereby we find that U/µ = 0.41. An essential feature of the NLDE is that † We note that the Dirac cone approximation is not necessarily adhered to in analogous honeycomb photonic lattice systems. See for example Refs. [64,65]. characteristic fluctuations are close enough to the Dirac point so that the linear Dirac cone approximation remains valid. Expanding the exact dispersion near the Dirac point, where k is a small deviation away from the Dirac point. The first term gives the linear Dirac dispersion. Higher order corrections describe curvature of the band structure away from the Dirac point. From the second order term we see that the NLDE description is valid for ak/ √ 8 1. This determines a lower bound on the wavelength for fluctuations of the condensate: λ min (2π/ √ 8)a. Linear dispersion places an additional constraint on the chemical potential: |µ| U + 6t h 101.9 nK. From the value of the chemical potential already obtained, we find µ = 2.59 nK 101.9 nK. Finally, weak short range interactions at very low temperatures justifies a lowest-band approximation to describe the physics of the NLDE.

Superfluid excitations near a Dirac point
The mean-field physics of single-particle states for a collection of fermions with Fermi energy near a Dirac point of a honeycomb lattice has been studied exhaustively and is discussed in various comprehensive articles [68,69,70,71]. For systems of bosons, however, one must carefully consider the meaning of condensation in the presence of Dirac points. To discuss BECs and Dirac points together one must address the compatibility of single-valuedness for phase functions required for stable vortex formation in a proper superfluid description, with the half-angle phase winding when circumnavigating a single Dirac point, i.e., the geometric or Berry phase [72].

Geometric and dynamical phase structure
To address these issues, we first review some relevant information treated in most review articles on graphene, as this information is true for cold bosonic atoms as well [27]. The single-particle spectrum of the honeycomb lattice exhibits zero-points, or Dirac points, in the reciprocal lattice associated with crystal momentum K = (0, ±4π/3a) rotated by 0, 2π/3, 4π/3, where a is the lattice constant shown in Fig. 2. Dirac points occur when the crystal momentum is tuned to the natural periodicity of the lattice with standing waves established due to Bragg scattering of the wave function. Reflection at the Brillouin zone edge is shown in Fig. 2, where one adds up projections of the vectors n 1 and n 3 along the direction of the crystal momentum vector K connecting points on the A sublattice of equal phase, A 1 and A 2 , to a third point A 3 . In particular, at the Dirac point this sum results in a net 2π accumulated phase angle at A 3 . In Fig. 2, the A and B sublattice wavefront density peaks are shown as red and blue dashed lines, respectively. In the tight-binding limit, the full lattice Hamiltonian reduces to two operators which couple the degenerate triangular A and B sublattices. The single-particle dispersion is computed by solving the 2 × 2 eigenvalue problem in momentum space determined by the Hamiltonian where the matrix elements come from computing the sublattice hopping energies Specifically, one finds The eigenfunctions ofĤ in Eq. (8) are with eigenvalues ±|E(k)|. Physically, the parameter ∈ R comes from an extra U(1) phase degeneracy and reflects the gapless symmetry of the system under spatial translations of the atomic density at the Dirac point. The matrix in Eq. (8) describes the amplitude and phase associated with real-particle tunneling between neighboring lattice sites. In particular, the phase of the wavefunction gets multiplied alternately by factors of e ±iφ(k) , so that no net phase is accrued when circumnavigating a closed path in the lattice. In contrast, long wavelength modes propagating in the lattice are described by linearizing the phase angle φ(k) so that the local lattice scale variations in the phase structure are neglected, in which case one should expect a net phase accumulation.
We are particularly interested in this net geometric or Berry phase since we must factor it into the phase winding for vortex solutions of the NLDE. Although most treatments of the subject use a momentum space argument, here we use instead a more direct analysis in real space. We expand the Hamiltonian and eigenstates near the Dirac point by taking k = K + δk, with K = (4π/3a)ŷ and δk the small expansion parameter, i.e., we consider small deviations from the sublattice Brillouin zone corner. In real space, this amounts to a derivative expansion of Eq. (9) in terms of the directional derivatives n 1 · ∇ and n 2 · ∇. The first-order term gives the massless Dirac Hamiltonian and Dirac equation [27], while higher products of derivatives provide corrections that probe the finer details of Bragg scattering around the Dirac point. To isolate the geometric phase, we consider adiabatic transport around a closed loop. Adiabaticity ensures that we do not accumulate a dynamical contribution to the phase and restricts the path to energy eigenstates states nearest the Dirac point. A direct way to accomplish this is by linearizing Eq. (8) in real space, solving for the eigenstates in plane-polar coordinates r and θ, and restricting to paths with large radii R ≡ r a. The large radius limit allows us to access only the longest wavelength modes that vary mainly tangentially with minimal radial contribution. Equation (8) then reduces tô from which we find the eigenstates and energies ± ω/2, with ω ≡ c l /R and c l the effective speed of light. Note that in this limit the degeneracy in Eq. (11) is lifted and eigenstates are forced into the form Eq. (13), which acquires a net phase of π under a full 2π rotation. Thus, linearization of Eq. (8) leads to a double wrapping of the phase angle φ around the polar angle θ.
In the case of a vortex, we can compensate for the Berry phase by requiring half-winding in the overall dynamical phase multiplying the spinor order parameter: exp[i( + 1/2)θ]. As a result, the geometric phase becomes identified with the relative phase between the two sublattices. Hence, stable vortices are required to have halfinteger internal geometric winding plus an overall half-integer dynamical winding such that the superfluid velocity is the sum of the gradients of both phases, effectively splicing together the internal and external phase.

Superfluid regime and dimensional analysis
To address superfluidity in the honeycomb lattice using a semiclassical approach, consider a thermal excitation with crystal momentum p (measured from the Dirac point) interacting with an atomic gas at the Dirac point (p Dirac ≡ 0) producing an excitation in the gas with momentum p . It follows from energy conservation, ∆E = 0, that In Eq. (14), for generality we assume that the incoming thermal mode may be in a Bloch state far enough removed from the Dirac point such that second-order corrections are important. Thus, m * is the effective mass related to the dispersion curvature, E(p ) is the energy of the quasiparticle excitation in the gas, and the upper and lower signs refer to negative and positive dispersion branches, respectively. We first examine the linear regime for which p ≈ p << c l m * , in which case we can neglect quadratic terms in Eq. (14). Keeping only linear terms and using |p − p | = p 2 − 2pp cosθ + p 2 , with p = |p| and θ the angle between p and p , Eq. (14) forces the constraint and four conditions determined by the different sign combinations for incoming and scattered modes. When the signs of the incoming thermal and scattered condensate modes are the same we find that θ = 0. On the other hand, if the energies of incoming and scattered modes have opposite sign, we obtain θ = π, thus scattering in the reverse direction occurs between Dirac particles and anti-particles as one should expect. Notice that in this linear regime conservation of energy places no additional constraints on p and p , so that in our mean-field analysis an equilibrium between condensate and non condensate atoms is maintained: the incoming mode transfers all of its energy and momentum to an excitation of the condensate leaving a single outgoing excitation with the same energy and momentum. A second regime of Eq. (14) corresponds to the condition p << p < c l m * , which yields the constraint leading to the same conditions as in the previous case for the scattering angle θ, but now with an additional constraint for an upper bound critical momentum p c below which no excitation can be created in the condensate Equation (17) recovers Landau's criterion for superfluidity but here in terms of the absolute value of the quasiparticle energy to account for scattering into negative energy states. We point out that the absolute value in Eq. (17) is a strictly a consequence of energy conservation and the presence of quadratic terms in Eq. (14); the various sign combinations in Eq. (16) are taken into account through the scattering angle θ.
With E(p ) = ±c l p, the upper critical bound is just p c = m * c l . Below this value (and for p >> p ) thermal modes cannot interact with the condensate, thus superfluidity is preserved. For p > p , however, as expected we see a breakdown of superfluity. In our analysis we have nowhere included details of the interaction; only a knowledge of states near the Dirac point was needed. Once we consider quantum effects and details of the interaction our results will change significantly. Consider again a non-condensate excitation with initial momentum p interacting with the condensate by transferring all of its momentum and energy to the condensate. A secondary excitation is then emitted with exactly the same momentum and energy. Since the initial and final excitations are indistinguishable, we can view this process as a single excitation interacting weakly with the condensate and continuing on its way with only an average self-energy correction. At the quantum level and to first order in the interaction U , a single interaction point must be averaged over the volume (area in quasi-2D) of the condensate. Since we are dealing with very long wavelengths the result is a nonlocal collective excitation formed as a composite of the initial incoming mode dressed by the condensate background. This effectively pairs orthogonal degrees of freedom in k-space. In contrast, at shorter wavelengths (higher energies) the incoming excitation couples with a quasiparticle locally, so that the available states for thermal and condensate modes remain distinct. Elaborating further, the number of accessible states less than k for the undressed excitation plus condensate is Ω(k) = a r k r , and the dispersion is E = ±c l k s . Here for our argument we leave the constants r, s, and a r general. Thus, Ω(E) = 2a r E r/s /( c l ) r/s where the extra factor of 2 accounts for both positive and negative eigenvalues. This yields the density of states In order to maintain D(E) constant when transitioning from short to long wavelengths, c l p > U → c l p < U , imposing the dimensional reduction Ω(k) = a r k r → a r k r/2 requires that we also take s → s/2. The renormalized energy is then E ∝ ±k 1/2 (for s = 1). The proportionality constant must involve the quasi-2D interaction U which we determine through dimensional analysis to be Note that Eq. (19) leaves out the possible form E(p) = ± √ −U c l p, which displays a low-momentum dynamical instability. However, this is regularized by accounting for a finite system size which imposes a lower momentum cutoff |p| min = 2π /R, where R is the radial size of the condensate. For the usual harmonic trap with frequency Ω we have R = ( /M Ω) 1/2 . By dimensional analysis, in terms of the quasi-2D renormalized average particle densityn and interaction U , we obtain the stability requirement for the oscillator length R ≤ c l π /(nU ). From a practical standpoint the lower bound |p| min removes the longest wavelength modes which opens an insulating buffer between the positive and negative parts of the spectrum in addition to regulating the dynamical instability.

Relativistic linear stability equations
Bogoliubov's method was originally introduced in his 1947 paper [73] (see also [74,75] for thorough contemporary treatments), and the concept later generalized by Fetter [76] to accommodate nonuniform condensate profiles. The latter formulation gives a convenient method for computing quasiparticle states and the associated eigenvalues by substituting the spatial functions for a particular background condensate into a pair of coupled differential equations, and then solving the resulting eigenvalue problem. Fetter's extended method was designed with a vortex profile in mind, and has proven successful for computing stability of vortices in trapped condensates, but also for gaining a deeper understanding of general vortex dynamics [77,78,79]. The set of equations that we derive in this section form the counterpart to Fetter's equations, but for trapped condensates that exhibit Dirac points in their dispersion [67]. We call them relativistic linear stability equations because of the quasi-relativistic context here and the similarity to equations that appear in relativistic fluid dynamics. It is noteworthy that our result is not limited to the honeycomb optical lattice but applies generically to any system where the linear dispersion and bipartite structure are present, and where the contact interaction between constituent bosons is weak. Mathematically, the essence of our derivation is contained in two steps: (1) transformation of a spatially continuous second quantized Hamiltonian into a spatially discrete one through an operation F; and (2) diagonalization of the Hamiltonian with an appropriate unitary transformation G. The effect of F is to take the system from the continuum to the tight-binding limit on the lattice, and G is equivalent to a Bogoliubov rotation. We will see that the final result is independent of the order of these operations so that the full procedure can be summarized abstractly in Fig. 3.
Our derivation of the RLSE relies fundamentally on Bogoliubov's method [73] as the underlying principle, and refers to Fetter's work [76] for technical considerations regarding nonuniform condensates. First, we recall the second-quantized many-body Hamiltonian for weakly interacting bosonŝ where Here, V (r) is the lattice potential and g is the strength of the contact interaction. The first step is to decompose the wavefunction as the sumψ(r) = ζ(r)â 0 +φ(r), where we have split the wavefunction into a part that describes the condensate (first term) and satisfies the bosonic commutation relation [â 0 ,â † 0 ] = 1, and a second part that describes small quasiparticle fluctuations. The operator in the first term destroys a particle in the mean-field ζ, which, by itself, is a good approximation tô ψ. The second term destroys a particle in a number of single particle basis states of the noninteracting system, and describes the part ofψ that deviates from the mean field. Taking the Bogoliubov limit requiresâ 0 → N 1/2 0 , where N 0 is the total number of condensed atoms, but we choose to compute the commutator before taking this limit in order to retain the effect of the presence of a macroscopic condensate field. We can obtain the commutation relations forφ andφ † by knowing thatψ, ψ † , andâ 0 andâ † 0 obey bosonic commutation relations. We obtain the quasiparticle commutation relations: φ (r) ,φ † (r ) = δ(r − r ) − ζ(r) ζ * (r ), φ (r) ,φ(r ) = 0, φ † (r) ,φ † (r ) = 0. In the Bogoliubov limit the condensate wavefunction has no operator part, in which caseψ may be written asψ(r) = Ψ(r) +φ(r). The condensate wavefunction has well defined phase and particle density and so may be expressed as: Ψ(r) = N 0 /A e iS(r) ρ(r), where A is the area covered by the planar condensate. Note that the radial part is normalized as A −1 d 2 r ρ(r) = 1. With these definitions, the usual bosonic commutation relations become: φ (r),φ † (r ) = e iS(r) e −iS(r )δ (r, r ), wherē δ(r, r ) = δ(r − r ) − A −1 ρ(r) ρ(r ).
Next, we transform to the new Hamiltonian defined byK =Ĥ − µN = H − µ d 2 rψ †ψ , then expand through second order in the operator part eliminating the linear terms by forcing the condensate wavefunction to satisfy the constraint (H 0 − µ + g |Ψ| 2 )Ψ = 0. We arrive at the Bogoliubov HamiltonianK =K 0 +K 2 , wherein zero-order and second-order operator terms are grouped intoK 0 andK 2 respectively. These are defined aŝ Note that in addition to the kinetic operator we also have an arbitrary external potential in the first two terms, which in our case will be the periodic potential of the optical lattice. Equation (23) is quadratic in the field operators and so may be diagonalized with the appropriate field redefinition. To diagonalize Eq.(23) we first apply the linear transformationφ(r) = e iS(r) where the prime notation on the summation sign indicates that we are omitting the condensate from the sum. Theα j 's andα † j 's inherit standard bosonic commutation relations fromφ andφ † , and the spatially dependent transformation coefficients u j (r) and v j (r) obey the completeness relations So far, our discussion has taken place in two continuous spatial dimensions constrained only at the boundary by a trapping potential. We now want to translate to a formalism that fits a two-dimensional periodic optical lattice potential with honeycomb geometry. This is done by assuming a tight-binding limit at each lattice site. Formally, this corresponds to expanding the wavefunction in terms of a Wannier basis: functions which are localized and centered on each lattice site. The nearestneighbor approximation then allows for a decomposition of the condensate and operator parts in terms of individual sublattices labeled A and B. In this new basis, the spatial dependence of the condensate and quasiparticle functions follows

First method: Tight-binding limit followed by diagonalization of quasiparticle Hamiltonian
We substitute Eqs. (27)-(28) into the Hamiltonian, Eq. (23), then take the longwavelength limit while translating the exponential (crystal) momentum factors to coincide with the Dirac point. The continuum limit effectively converts the sublattice sums into integrals. By performing one of the integrations, over the A sublattice, say, while adhering to nearest neighbor overlaps, we obtain the affective Hamiltonian for the condensate and quasiparticlesĤ =K 0 +K 2 wherê Here we have defined the condensate two-spinor in terms of the A and B sublattice components Ψ(r) ≡ [ψ A (r), ψ B (r)] T , and the Dirac operator is defined as D ≡ ∂ x + i∂ y . Next, we isolate the first six terms (terms with the daggered operator to the right) and write them as a matrix contraction of two pure operator valued vectors where The eigenvalues are then obtained by and the corresponding eigenvectors follow The unitary matrix that diagonalizes Eq. (31) is The first six terms in Eq. (30) may be expressed in the new basis as where we have included the j, k subscripts on the eigenvalues to be fully descriptive.
The new quasiparticle operators can be written in terms of the old ones aŝ Note that the right hand side is k-dependent which is implied on the left. The substance of the transformation is contained in the momentum and space-dependent eigenvalues The next step is to constrain the quasiparticle amplitudes in Eq. (45) (the u's and v's) in order to diagonalize the Hamiltonian with respect to the momentum indices j and k. First, we let and then substitute these into Eq. (45), which reduces the two eigenvalues to , and where we have reinserted the chemical potential terms. It is important that Eq. (45) depend only on one index so that quasiparticle amplitudes for different eigeneneregies are not coupled. Dividing Eq. (45) through by v j,A and v j,B , respectively, cancels all j-index terms except for ones that appear as u j,A /v j,A and u j,B /v j,B . To completely decouple the j-k modes, we must ensure that u j,A /v j,A = u j,B /v j,B = η(r ), i.e., the amplitudes for any given quasiparticle mode have the same relative spatial form. We can then rewrite λ +{jk} as . Finally, we impose the constraints Multiplying Eqs. (48)-(49) by v j,A and v * k,A , respectively, and using the property that , we may separate out 1/2 of each derivative term in Eq. (48), which reduces the non-derivative terms in the first line of Eq. (48) to We may reduce the second line using the other half of each derivative term, thereby condensing the eigenvalues down to The next six terms in Eq.(30) may be diagonalized in a similar way yielding the eigenvalues and Following our previous steps, we obtain Combining Eqs. (51) and (54), and inserting the quasiparticle operators, reduces the first twelve terms in Eq. (30) to the expression For the special case where j = k, we may further combine the terms at the cost of an extra c-number term to arrive at Applying the completeness relations Diagonalizing the rest of Eq. (30) (terms with no daggered operators and ones with only daggered operators) by capitalizing on the j-k symmetry of terms such as d 2 r u k,A v j,A , and anti-symmetry of the (E j − E k ) factor, we obtain the final form of the interacting Hamiltonian with the resulting constraints on quasiparticle amplitudes given by

Second method: Diagonalize quasiparticle Hamiltonian then impose tight-binding
Although the first method is cumbersome it is the more rigorous approach and instils confidence in the final constraint equations. A shorter approach is to first obtain the usual Bogoliubov equations for a condensate not confined in a lattice, and then apply the tight-binding limit directly. The Bogoliubov Hamiltonian iŝ with the constraint equations (BdGE) given by In Eqs. (65)-(66), L is a differential operator containing terms that couple the quasiparticle and condensate velocities. An additional implicit constraint is that Ψ satisfies the nonlinear Schrödinger equation. To pass to the tight-binding limit we express all spatial functions in Eqs. (64)- (66) in terms of Wannier functions for the individual sublattices, and evaluate the Bloch plane wave factors at the Dirac point momemtum. Adhering to nearest-neighbor overlap for on-site Wannier functions, we integrate out spatial degrees of freedom (which splits the honeycomb lattice into A and B sublattices), regroup terms into finite differences, and then take the continuum limit. Equation (64) where φ is the condensate phase. After going through the steps that culminate in the tight-binding continuum limit, these terms transform to where the coefficients encapsulate the spatial integrals as follows: These extra terms depend on the condensate phase φ A(B) , and so couple the superfluid velocity to the quasiparticle excitations. In particular, the term with coefficient τ 1 depends on the direction of quasiparticle emission relative to the motion of the condensate. The relativistic linear stability equations, Eqs. (60)-(63), may be expressed in compact notation as where and 1 2 is the 2 × 2 unit matrix.

Stability of vortex solutions
Two independent derivations of the RLSE in Sec. 4 and their reduction to the BdGE, which we discuss in Sec. 6, establishes Eqs. (69)-(74) as the correct method for computing the low-energy structure (quasiparticle states and eigenenergies) for arbitrary vortex solutions of the NLDE [43]. The most immediate and pragmatic concern is the combined effect of the honeycomb lattice geometry and the inter-particle interaction on the lifetime of a vortex. It should be emphasized that the presence of an infinite tower of negative energy states below the Dirac point seems to imply that a condensate residing there will eventually decay provided there is a mechanism for energy dissipation into noncondensate modes (i.e., secondary interactions with thermal atoms). † Generically, negative energy states are present for moving condensates for which excitations subtended by a backward cone have negative frequencies [76]. Moreover, when a vortex is present small displacements of the core from the symmetry axis of the trap results in a precession of the core, which, when combined with dissipation, causes the vortex to spiral to the edge of the condensate. In the absence of a periodic lattice potential this dynamical process is known to be driven by the anomalous modes in the linear spectrum, i.e., modes with negative energy and positive norm [79] also called Goldstone modes. The time for a vortex to spiral to the edge of the trap would then define its lifetime. In the absence of the lattice this precessional motion is canceled by introducing rotation to the trap [79,77], a result which we suspect to be true in the lattice case as well.
To undertake a full treatment of the lifetime would mean computing this spiraling time and then comparing it with the lifetime that we compute here due to the dynamical instability from the complex frequencies. The lifetime of the vortex would then be the smaller of the two values. Nevertheless, in cases where dissipation is weak and the vortex is centered on the symmetry axis of the trap, the dominant source of instability arises from the complex eigenvalues associated with RLSE modes. We will limit our analysis to the effect of the latter, and regard the negative real part of the eigenvalues from a standpoint of metastability. Physically, the complex eigenvalue gives rise to fluctuations in the angular rotation of the vortex spinor components [67]. In the case of the NLDE this is a result of internal "friction", i.e., energy exchange, between the two spinor components displayed in the complex derivative terms of the Dirac kinetic energy. This drag force between the two vortex components (or between vortex and soliton) eventually causes substantial degradation of the vortex. This is the measure that we will use to compute vortex lifetimes.

Numerical solution of the relativistic linear stability equations and vortex lifetimes
The stability of a particular condensate density and phase profile such as a vortex is arrived at by expanding Eq. (69) and expressing differential operators in terms of suitable coordinates, for example polar coordinates for a vortex, then using separation of variables for the quasiparticle amplitudes with the appropriate form of ψ A(B) . This yields a set of first-order coupled ODE's in the radial coordinate to be solved consistently for the functions u A(B) (r), v A(B) (r) and the eigenvalues E k . We discretize the derivatives and functions using a forward-backward average finite-difference scheme, then solve the resulting discrete matrix eigenvalue problem using MATLAB function Eig. In Fig. 4 we have plotted the real and imaginary parts of the first 20 eigenvalues, labeled by the † We remind the reader that this infinite tower of negative energy states is only in the Dirac cone approximation.
quantized quasiparticle rotation number n ∈ Z, for the vortex/soliton solution which we discuss in previous work [43]. The lowest modes are anomalous with negative real parts and positive, nonzero but small, imaginary parts.  Convergence of RLSE eigenvalues for the l = 1 vortex/soliton background as a function of the grid size N used in the 4N × 4N matrix problem is displayed in Fig. 5, where we have plotted the real and imaginary parts of the eigenvalue for the lowest excitation mode.  The lifetime of a particular vortex solution can be computed by examining the lowest quasiparticle rotation mode n = 1, since at very low temperatures this mode dominates the spectrum. The lifetime is then characterized by the reciprocal of the imaginary part of the associated eigenvalue, i.e., lifetime ≡ /Im (E −1 ). Here, the −1 subscript refers to quasiparticle rotation relative to the vortex rotation. Eigenvalues for the lowest quasiparticle rotational mode and the associated lifetimes for all of our solutions are listed in Table 2.  Table 2. Stability properties of NLDE vortices. Lifetimes are computed using the value of the interaction U in Table 1 and the formula lifetime = /Im (E −1 ).
To understand the character of the quasiparticle modes we must consider the spatial functions associated with each eigenvalue. The radial quasiparticle functions have the forms previously determined in Fig. 3 of Ref. [67], which shows in particular the lowest excitation mode (n = 1) of the vortex/soliton ( = 1) solution. They are bound states near the core of the vortex localized specifically at the point where the soliton and vortex components of the background are equal (for radial plots of NLDE vortex solutions see Ref. [41]). Physically, the imaginary part of the eigenvalues imply a transfer of energy between the vortex and soliton components through quantum fluctuations. In particular, each component acquires quantum admixtures from different rotational modes as well as local shifts in amplitude from phase and density fluctuations, respectively. Mathematically, the full quasiparticle operator with time and spatial dependence for this mode isφ(r, t) φ A,−1 (r, t),φ B,−1 (r, t) T , where the quasiparticle spinor operators arê As discussed previously, relative to the vortex background the quasiparticle has rotation = −n = −1, which has the effect of reducing the rotation of the vortex. Note that the expression for the operatorφ(r, t) is approximate since we have truncated the sum over quasiparticle modes after the lowest mode. We recall that the spatial functions have the properties u A,−1 (r), u B,−1 (r) ∼ 10 −2 and v A,−1 (r), v B,−1 (r) ∼ 10 −5 (see Ref. [67]), where all are peaked in the "notch" region ξ Dirac < r < 2ξ Dirac , and where the absolute values of the slopes of the soliton and vortex are maximum. In this region, the normalization integrals (one for each sublattice) are given by This combination of positive norm and negative Re(E −1 ) signals the presence of the anomalous mode. In Sec. 6, we will see that these bound quasiparticle modes solve the Majorana equation, which predicts an additional zero energy mode localized at the same distance from the center of the vortex.

Connection to other theories
In this section we examine several reductions of the RLSE to other equations familiar to BECs, superconductivity, graphene, and high energy physics. Our results demonstrate the variety of substructures contained within the RLSE framework. Note that we adhere to the weakly interacting regime through all of our derivations as explained in Sec. 2.

Reductions of the relativistic linear stability equations
To begin, we look for a non-relativistic reduction of Eq. (69) by working first from the lattice form of the NLDE since the hopping terms are the same as those for the RLSE. We recall that the standard massive Dirac equation has a well defined non-relativistic limit to the Schrödinger equation. The proof uses the fact that in the low-energy limit the mass term (proportional to mc 2 ) is the largest contribution to the energy. The two-dimensional formulation reduces to two coupled equations. The mass term is isolated in one equation then substituted into the second equation. The substitution converts the first-order spatial gradient to a second-order Schrödinger kinetic term, and pushes the mass dependence into smaller correction terms. A similar procedure may be implemented in our case but we must first introduce an offset between the sublattice potential well depths (a mass gap) so that we obtain the desired curvature in the spectrum near the Dirac points, effectively opening up a non-relativistic regime. Starting from the discrete NLDE for a single Dirac point and following similar steps as in our previous work [27], we obtain where t h , t 0 , U , and k are the hopping, same site, and on-site interaction energies and crystal momentum, respectively. The δ's, n's, and 2D vector indices j indicate the same lattice vectors described in our original derivation of the NLDE [27]. In Eqs. (78)-(79), t 0 is the sublattice offset equivalent to a spectral gap 2|t 0 |. For weak interactions, the on-site energy can be made much larger than the contact interaction strength by tuning the lattice potential so that |µ ± t 0 | >> U . After inserting the correct values for the lattice vectors and solving Eq. (79) for ψ B j , to zeroth-order in U/|µ − t 0 | we obtain From Eq. (80) we may write analogous expressions for neighboring sites by shifting the indices using the lattice vectors n j , Substituting Eqs. (80)-(82) into Eq. (78), expanding complex factors and regrouping the terms to form finite differences, we arrive at the expression Equation (83) is a discrete nonlinear Schrödinger equation for the honeycomb lattice in the sense that it has as its continuum limit the usual nonlinear Schrödinger equation with cubic nonlinearity. Substituting the correct continuum forms for the finite differences and then expressing the result in rectangular coordinates, we obtain with imaginary mass in Eq. (89), and an ordinary Klein-Gordon mode with real mass in Eq. (90). In contrast, if we tune the lattice potential offset so that t 0 ∼ µ the mode described by Eq. (85) has a very small effective mass and large energy, whereas the mode in Eq. (86) will have a very large mass and small energy. In this case the mode in Eq. (86) gets "frozen out" and we are left with only one propagating mode in Eq. (85). Here multiplication by the total energy µ + t 0 does not cancel the effective mass µ − t 0 in the denominator of the gradient term. Reintroducing the time dependence by µ + t 0 → i ∂ t and the effective mass m = (µ − t 0 )/2c 2 l , Eq. (85) reduces to the nonlinear Schrödinger Thus, tuning t 0 interpolates between a Dirac and a Schrödinger structure with Klein-Gordon bridging the two. Applying the same steps to the RLSE, Eq. (69) yields the BdGE with ∆ m = −µ + 2U |ψ| 2 , ∆ p = U |ψ| 2 , where ψ is the condensate wavefunction for either of the decoupled sublattices and the effective mass is the same as in Eq. (91), m = (µ − t 0 )/2c 2 l . Note that we have suppressed explicit space-time dependence in Eq. (92) for clarity. In the particle regime for large characteristic momentum, c l |p| U , the particle and hole amplitudes satisfy u v and Eq. (92) reduces to the standard Schrödinger equation for a particle moving in the potential V ≡ −∆ m .
Next, we look at the case of single-mode approximation for the pseudospin degrees of freedom in Eq. (69), i.e., where the sublattice backgrounds are equal ψ A ≡ ψ B which also implies that u A = u B ≡ u(r, t) and v A = v B ≡ v(r, t). One then finds that the system Eq. (69) reduces to the Andreev equation The unit vectorp in Eq. (93) points in the direction of quasiparticle propagation.
Here we have chosen the case of zero background flow ∇φ A,B = 0 as we will do for the remainder of this section except for the vortex background. Equation (93) is the Andreev equation for propagation through a medium comprised of both normal and superconducting regions [53]. The spatially dependent pairing and mass terms are ∆ p (r) = U |ψ(r)| 2 and ∆ m (r) = 2U |ψ(r)| 2 − µ. In this analogy the condensate wavefunction ψ(r) stands in for the order parameter in a superconducting medium. Equation (93) describes slowly varying particle and hole functions u(r) and v(r) split off from an overall rapidly oscillating plane wave portion which moves in the direction † This step can be justified formally from the Heisenberg equation of motion for the wavefunction starting from the operator formalism, but such justification is well known from theory of NLSE.
p. Thus, we should expect similar exotic scattering such as specular and retroreflection [80]. Next, we look at the particle regime where the particle component is dominant, (B) . In this regime Eq. (69) reduces to the Dirac equation where a potential term appears ∆ A(B) (r) = 2U |ψ A(B) (r)| 2 − µ. Equation (94) further reduces to the massless Dirac equation in the case of a constant background |ψ A(B) | 2 ≡ µ/(2U ). Interestingly, zero mode solutions (E = 0) of the RLSE occur as well and we find that these solve the Majorana equation which is implicit in the RLSE for certain background configurations. To see this we set E k = 0 in Eq. (69), which decouples the system into two sets of equations in the extreme long-wavelength regime characterized by |u A(B) | = |v A(B) |. In this regime Eq. (69) gives two copies of the form where the potential terms are ∆ A(B) = U |ψ A(B) | 2 − µ. For a uniform condensate, i.e., far from any vortex cores, the asymptotic choices are U |ψ A(B) | 2 → µ, 0. In both cases Eq. (95) offers no solution. However, for the vortex/soliton there is a "notch" in the order parameter near the core, where |ψ A | 2 = |ψ B | 2 < µ/U ⇒ ∆ A(B) < 0, in which case Eq. (95) reduces to the Majorana equation , m ≡ |∆ A(B) |, and ψ : R 2 → R 2 . Equation (96) supports real solutions with linear dispersion and has been studied extensively in its original mathematical form [55] and more recently in condensed matter physics intimately associated with topological insulators [57]. In their present incarnation these Majorana zero modes also occur in the core of nonlinear Dirac vortices with higher winding ( > 1 in Ref. [43]) where both spinor components vanish |ψ A(B) (0, θ)| 2 = 0. In this case the mass term in Eq. (96) reduces to the condensate chemical potential m = µ. In the superfluid context the meaning of the Majorana zero mode is of a zero-energy pure spatial density fluctuation associated with rigid translations of the vortex core. Here phase fluctuations only appear as finite-energy fluctuations in the vortex rotational and translational motion. For the vortex/soliton the zero mode is a circular ring reflecting the symmetry under both rigid rotations as well as translations of the vortex. In Fig. 6 we summarize the various types of reductions of the RLSE indicating the conditions or limits for each equation type.

Mapping to relativistic Bardeen-Cooper-Schrieffer theory
In this section we discuss the modifications needed to connect the RLSE to relativistic BCS theory. Here we capitalize on an important property of the NLDE and RLSE.  This is that repulsive interactions for bosons in the honeycomb lattice break the valley particle-hole exchange symmetry at the Dirac point in a significant way such that an additional sign change of the interaction restores the symmetry. More properly stated, the noninteracting theory is invariant independently under charge conjugation (C), parity inversion (P), and time reversal (T). Repulsive interactions break T and C, but the symmetry-breaking cancels in such a way as to preserve the full CPT symmetry [27]. Consequently, a parity inverted positive energy solution (valley particle) can be interpreted as a negative energy solution (valley hole) in a theory with attractive interactions but without parity inversion. Stated differently, a theory of particles with repulsive interactions is equivalent to a theory of holes with attractive interactions.
To complete the mapping to BCS theory we introduce a mass term and nearestneighbor interactions at the lattice scale to couple the different spinor components. The mass term is obtained through an asymmetry in the honeycomb sublattice potential depths, an intermediate step in populating Dirac points, as we have explained in [43]. The various types of relativistically invariant interactions may be constructed using nearest-neighbor interactions as follows. Specifically, the symmetry of the nonlinearity in the NLDE determines the symmetry of the superconducting order parameter and pair potential in the corresponding BCS analog equations [45]. The vector-vector interaction can be obtained by including repulsive nearest-neighbor interactions. A scalar-scalar type coupling can be realized similarly, but by using attractive (instead of repulsive) nearest-neighbor interactions in addition to the repulsive on-site interactions. The spin and pseudo-spin symmetric terms are characterized by an alternating sign for the coupling between the two spinor components. This type of coupling may be realized in a lattice setting via Feshbach resonances using a beam with the proper spatial modulation to produce interactions whose sign alternates from between neighboring lattice site. Pseudo-scalar forms can be realized by eliminating on-site interactions while retaining repulsive nearest-neighbor interactions.
The case of scalar-scalar coupling in the NLDE with equal on-site and nearestneighbor interactions U = U nn and mass term m s c 2 l (see ref. [43]) elevate the RLSE to the form i c l σ · ∇ + m s c 2 l · 1 2 + qσ µ A µ −i∆ p σ y i∆ p σ y −i c l σ · ∇ + m s c 2 l · 1 2 + qσ µ A µ u k v k where ∆ p (r) ≡ U [|ψ A (r)| 2 + |ψ B (r)| 2 ] is the scalar pairing function, the effective polarized 4-vector potential in (2+1) dimensions (so reduced to 3 components) is A µ (r) ≡ (U/q) [|ψ A (r)| 2 − |ψ B (r)| 2 − µ/U, |ψ A (r)| 2 − |ψ B (r)| 2 , 0], and q is an effective charge. As before 1 2 is the two-dimensional unit matrix. Equations (97) comprise the relativistic Bogoliubov-de Gennes equations also known as the Dirac-Bogoliubov-de Gennes equations [45,46]. In the special case of a uniform condensate which solves the nonlinear Dirac equation we have |ψ A | 2 = |ψ B | 2 = µ/U and Eq. (97) yields the eigenvalues E k = ± ( c l k) 2 + (m s c 2 l ) 2 ± (m s c 2 l + µ) where the magnitude of the quasiparticle momentum k = |k| labels the eigenstates. The signs outside of the radical relate to pseudospin valley states and those inside the radical to the particle-hole Nambu states. The spectrum Eq. (98) is plotted in Fig. 7.

Conclusion
In this article we have delineated the various constraints required for stabilizing a BEC at Dirac points of a honeycomb optical lattice. Energetically, we find that the Bose gas must be weakly interacting with excitations in the transverse direction suppressed relative to longitudinal ones. The latter condition can be implemented by using a relatively small vertical trap size. Additionally, Bloch states for the Bose gas must remain near enough to the Dirac point crystal momentum so that secondorder band distortions are negligible. This condition is equivalent to the requirement that quasiparticle momenta remain much less than the Dirac point crystal momentum. Length constraints include a large quasi-2D effective healing length relative to the lattice spacing so that a continuum theory is physically sensible. Atomic and lattice parameters are related primarily by imposing the usual Landau criterion for dynamical stability, which relates the effective speed of light (lattice parameters) to the quasi-2D renormalized speed of sound (atomic parameters). We performed a detailed analysis of lifetimes for nonlinear Dirac vortices, elucidating the low-energy landscape for each solution type.
Vortex lifetimes were computed based on dynamical instabilities induced by quantum fluctuations: complex eigenvalues appear in the linear spectrum for all vortex types. These include a complex topological vortex, topological vortices with generic winding, ringvortex, ring-vortex/soliton, vortex/soliton, Mermin-Ho, Anderson-Toulouse, and halfquantum vortices. The longest lived vortices are the ring-vortex, ring-vortex/soliton, vortex/soliton, and Anderson-Toulouse vortex with lifetimes 0.5295 s, 4.043 s, 3.841 s, and 4.041 s, respectively. A significant part of our work was devoted to the derivation and analysis of the relativistic linear stability equations (RLSE). We demonstrated that the RLSE reduce to several well known equations. The presence of a mass gap through an offset in the sublattice potential depths allows for an interpolation between the RLSE and the BdGE. By tuning the ratio of the gap to the chemical potential between small and large values, the governing equations for quasiparticles vary continuously between RLSE and Bogoliubov-de Gennes equations (BdGE) passing through a Klein-Gordon type structure associated with fluctuations of the nonlinear Klein-Gordon equation. In the particle regime where momenta are large compared to the interaction strength, the three types of stability equations reduce to the standard Dirac and Schrödinger equations with the Klein-Gordon equation interpolating between these. In the singlemode approximation, where the pseudospin valley spatial functions are equal, the RLSE reduce to the Andreev equations for electrons in inhomogeneous superconductors. For zero-energy modes residing at the core of a defect such as a vortex, the RLSE reduce to the Majorana equation with the Majorana mass determined by the local density of the condensate at the "notch" in the case of the vortex/soliton, and equal to the chemical potential in the general case of higher winding vortices ( > 1).
By including nearest-neighbor interactions and a mass gap we have shown that the RLSE transform to the Dirac-Bogoliubov-de-Gennes equations, which describe Cooper pairing of relativistic fermions. The additional Nambu space elevates the two-spinors in two spatial dimensions to a four component object consistent with our RLSE. The nonrelativistic limit is defined for quasiparticle momenta much smaller than the momentum scale set by the mass gap, in which case we recover standard BCS theory. In the analog picture the BCS pairing function is mapped to the total local condensate density, that is, the sum of squared moduli of the sublattice amplitudes. Superconductivity is strong or weak depending on the magnitude of the pairing function relative to the mass gap energy. We have shown that when the pairing function transforms as a scalar under the Lorentz group the absence of internal structure for the scalar term leaves an extra degree of freedom in the form of a vector potential. The difference in sublattice densities acts as an additional polarized vector potential acting on the pseudospin-Nambu spinor.
Interesting research directions that extend the work presented in this article could include elevating the boson-honeycomb lattice problem to a relativistic field theory. The lowest-band approximation would still be viable provided the theory is regularized by imposing an upper momentum cutoff at the lattice scale. The various classes of Lorentz quartic interactions may be constructed by including nearest-neighbor interactions in the lattice, as we have outlined in Sec. 6 of this article. It has been demonstrated that quartic interactions are fundamentally constrained by the conformal structure of all the terms of a particular relativistic Lagrangian [81]. Thus, by tuning the sign and strength of nearest-neighbor interactions it may be possible to observe quantum phase transitions in the superfluid phase between different conformal theories associated with various relativistic field theories.