{'fraction_non_alphanumeric': 0.038606042018095386, 'fraction_numerical': 0.012421407759546081, 'mean_word_length': 4.392805457928468, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Relativistic linear stability equations for the nonlinear Dirac equation in Bose-Einstein condensates

We present relativistic linear stability equations (RLSE) for quasi-relativistic cold atoms in a honeycomb optical lattice. These equations are derived from first principles and provide a method for computing stabilities of arbitrary localized solutions of the nonlinear Dirac equation (NLDE), a relativistic generalization of the nonlinear Schr\"odinger equation. We present a variety of such localized solutions: skyrmions, solitons, vortices, and half-quantum vortices, and study their stabilities via the RLSE. When applied to a uniform background, our calculations reveal an experimentally observable effect in the form of Cherenkov radiation. Remarkably, the Berry phase from the bipartite structure of the honeycomb lattice induces a boson-fermion transmutation in the quasi-particle operator statistics.

Progress in condensed matter and particle physics has been periodically marked by significant mutual exchanges between the two disciplines, many proposals for which are realized in model systems of ultracold quantum gases in optical lattices [1][2][3]. Recent active areas of research include holographic dualities such as AdS/CFT [4], theoretical constructions of superstrings in ultracold quantum gases [5], chiral confinement in quasi-relativistic Bose-Einstein condensates (BECs) [6], and our own derivation of the nonlinear Dirac equation (NLDE) describing ultracold bosons in a honeycomb optical lattice [7]. Our investigation into relativistic effects in BECs is motivated by this spirit of cross fertilization with the aim of tying in theory to experiment.
In this Letter, we develop the relativistic linear stability equations (RLSE) for the NLDE. Moreover, we find emergent nonlinear localized solutions [8] to the NLDE, including solitons, vortices, skyrmions, and half-quantum vortices, the latter so-far unobserved in BECs. Although most of these objects have been studied in multicomponent BECs, such models lie within the usual Schrödinger many-body paradigm. In contrast to this paradigm, our investigations reside within a relativistic framework in which the elementary excitations are governed by a Dirac-like equation. This provides a fundamentally different context distinguished by the presence of a Berry phase, so that exchange of two vortices leads to integer or halfinteger exchange statistics. It is not surprising that the elementary excitations in our theory exhibit a rich structure: a Dirac-like dispersion which obeys either bosonic or fermionic statistics depending on the strength of the contact interaction. Consequently, in order to determine the quasi-particle states and energies we cannot rely on the Bogoliubov-de Gennes equations (BdGE) since these are based on nonrelativistic quantum mechanics. Instead, we derive, from first principles, the RLSE which give the correct low energy dynamics for an arbitrary background condensate. The RLSE are reducible to the BdGE in certain limits, and so may naturally be considered relativistic generalizations of the latter. Based on the RLSE we predict Cherenkov radiation that can be measured in experiments: the combination of lattice and particle interactions results in a rich spatial distribution that is not seen in the BdGE for the uniform case [9].
In the laboratory, the NLDE can be obtained by cooling bosons into the lowest Bloch band of a honeycomb optical lattice [10]; the lattice is constructed by establishing three phase-locked interfering laser beams in a plane while freezing out excitations in the vertical direction as in Fig. 1. To obtain the desired Dirac structure, particles are first condensed into the lowest energy state (zero crystal momentum) of the lattice and then adiabatically trans- lated to the Dirac point at the band edge (see Fig. 1) by adiabatically tuning the relative phases between the laser beams. We emphasize that the Dirac point, which is key to the NLDE and our predictions, is maintained even in the presence of the shallow harmonic trap endemic to atomic BECs [11]. Notably, the NLDE may also be obtained by the alternative method of using a square optical lattice with a staggered gauge field induced by a time-dependent optical potential [12]. A similar arrangement using only bosons is expected to show the same low energy structure as in our model. Nonlinear phenomena in BECs have been studied extensively over the past decade [8], from single-component vortices in rotating, trapped BECs [13] to complex multi-component order parameters [14,15] resulting from interactions between the different components and the possibility of nontrivial topological windings of the internal symmetry space around a singular vortex core. Some form of BdGE analysis plays a central role in such constructions as a means of probing stability as well as for gaining a deeper understanding of the low-energy fluctuations.
Derivation of the RLSE. -Since the system we describe in this letter is a BEC confined strictly to two spatial dimensions, it is appropriate to recall the justification for such a construction before presenting the RLSE. For uniform 2D systems the Mermin-Wagner theorem forbids the formation of a true condensate defined by an infinite phase coherence length. This comes from the fact that the density of states diverges in the 2D case for finite T . Instead, one sees the formation of a quasi-condensate characterized by local phase coherence restricted to finite size regions. The size of these regions greatly exceeds the healing length so that all of our solutions are realizable in this picture. However, the inclusion of a harmonic confining potential allows the formation of a true 2D condensate. The potential places a lower bound on the energy for fluctuations and, since it is these long wavelength fluctua-tions that are responsible for destroying long range order, the trap provides a means of expanding the spatial range of validity of the mean field description.
To obtain the low energy excitations of solutions of the NLDE, we must find the correct set of equations that describe quasi-particle states analogous to the BdGE equations for the general case. These are obtained from the Hamiltonian for a weakly interacting Bose gas,Ĥ = drψ † H 0ψ + U 2 drψ †ψ †ψψ , H 0 ≡ 2 ∇ 2 /2M + V (r), and working through four steps [16]. (1) Takeψ = Ψ c (r) + δψ q (r) (condensate + quasi-particles), with δψ q small. (2) Impose a constraint on Ψ c to eliminate linear terms in δψ q (r), keep only quadratic terms in δψ q (r), and expand as a sum of particle and hole creation operators. (3) Invoke Bloch-state expansions for Ψ c (r) and δψ q (r) and take the lowest band. (4) Take the long-wavelength limit while taking momentum with respect to the Dirac point K, such that k ≪ q ≪ K, where q is the momentum of the condensate relative to the Dirac point K and k is the quasiparticle momentum measured relative to q, and finally diagonalize the quasi-particle part of the Hamiltonian. One finds the RLSE: where the matrix coefficients are defined as Here, D = (∂ x + i∂ y ) is the single particle Dirac operator. Also, Ψ = (Ψ A , Ψ B ) is the BEC order parameter at the K Dirac point, with normalization on sublattice components dr (|Ψ A | 2 + |Ψ B | 2 ) = 1. Analogous equations hold for the inequivalent Dirac point at −K. Cast in this highly compact form, Eqs. (1)-(2) are reminiscent of the BdGE and may be solved for the spinor quasi-particle amplitudes , v k,B (r)] T and the quasi-particle energy E k . The components of these 2-spinors represent quantum fluctuations of the sublattice condensate order parameters Ψ A and Ψ B which in general are nonuniform Cfunctions on the plane. The presence of the local phase of the condensate φ A(B) (r) indicates the complex interaction between the local superfluid velocity of the condensate v s,A(B) (r) ≡ ∇φ A(B) (r) and the spinor quasiparticles u k (r) and v k (r). We have taken = c l = 1 for simplicity, where c l is the effective speed of light in the NLDE. Note also that we have included an effective mass m eff (anisotropic lattice) that competes with the chemical potential µ.
The RLSE exhibit a positive-negative energy inversion symmetry which is found by complex conjugating Eqs. (1) and (2) followed by a spatial parity inversion. This is similar to the BdGE where negative energy modes are allowed by symmetry but are ignored since the underlying Hamiltonian is positive-definite. It is important to note that for a moving condensate, the negative-energy modes cannot be removed and are crucial indicators of Cherenkov radiation. However, in our case, the Dirac Hamiltonian is not positive-definite since our theory is defined at zero lattice energy, not the lowest energy Bloch state, so we must respect the presence of both energy raising and lowering modes. Another important feature is that the RLSE are reducible to the BdGE when the local lattice potential energy is the main contributor to the condensate chemical potential and the condensate is slowly varying (quasiuniform background), i.e., |µ| >> U , ) the Wannier functions.
Quasi-particle Operator Statistics. -Since we use only bosons in our construction, and we have shown that these collectively produce the spinor structure of Dirac theory, it is natural to ask what are the quantum statistics for collective excitations. When quantum fluctuations are weak, quasi-particle states are superpositions of direct products of bosonic operators and spinor singleparticle states: operators retain their bosonic structure. In contrast, when the depletion of the condensate becomes significant, quasi-particle operators at each lattice site become sums of superpositions of terms each with a different phase, i.e., the single-particle bosonic operatorsb andb † become entangled with the on-site phase. We can illustrate this more rigorously. Consider a BEC in the thermodynamic limit: 1/N → 0 and N/M = constant, where N is the total number of particles and M is the number of lattice sites. Also, we focus on the regime where U/t h ≪ 1 at T = 0. The ground state is then described by the superfluid density order parameter and the quasi-particles are coherent long wavelength fluctuations in the phase. A Gutzwiller ansatz in terms of on-site number states provides an adequate formulation of the wavefunction. If we consider increasing the particle interaction U towards the critical value U c that separates the superfluid from the Mott-insulating phase, we would expect to see an increase in depletion of the condensate. The part of the wavefunction that describes particles outside the condensate has complete phase decoherence so that the phase fluctuations at the i th site, for N ′ i particles outside the condensate, can be described using the eigenstates of the phase operator in the number state basis [17]: where |0 i is the ground state wavefunction for the i th site. With this definition, the i th site phase operator isθ and θ 0 is a reference angle which must be averaged over when computing observables. In general, the depletion number N ′ i is a function of the total number of particles N i and the interaction U . When N ′ i is large, this dependence on N and U can be incorporated into a coefficient f i (N, U ) that multiplies each term of the sum while maintaining the relative sizes N i << N ′ i << 1: where α i = f i e iθ0 and we have used the symmetric form for the phase coherent stateD(α) ≡ e αib † i −α * ibi . The normalized quasi-particle operators can now be defined having the correct Bogoliubov limit: With this definition, the average number of bosons forming a quasi-particle is b †b = 1 + |α| 2 , where |α| encodes the degree to which the bare quasi-particle operatorsb andb † are dressed or fused to the background phase. When working at the lattice scale it is clear that the reference angle θ 0 must be single-valued under one full rotation so that the dependence on the polar angle is e i(θ0+θ) . However, when translating to the continuum limit, the two triangular sublattices of the honeycomb lattice become identified with the same spatial point so that a 2 − 1 mapping of parameters is needed consistent with the double covering map p : SU (2) → SO(3). The continuum quasi-particle statistics can then be determined by computing the Berry phase (holonomy) by adiabatically transporting a quasi-particle through a suitable closed path and correctly accounting for the SU(2) structure of the single-particle states. The path must remain within a degenerate subspace of the Hamiltonian so that we may isolate the geometric phase from the dynamical phase.
Physical Parameters and Regimes. -We list first the fundamental dimensionful parameters that we use. They are as follows: the average particle density n 0 , the chemical potential µ, the lattice spacing a, the s-wave scattering length a s , the mass of the constituent bosons M , and the lattice well depth V 0 . Several relevant composite quantities may be constructed from these. These are the effective speed of light c l = t h a √ 3/2 , the sound speed c s = U n 0 /M , the interaction strength U = 4π a s /M , the healing length ξ = t h a √ 3/2 n 0 U , and the hopping energy t h = d 2 r w * Ĥ 0 w, where t h depends on a and V 0 , respectively, through the overlap of Wannier functions and the lattice potential insideĤ 0 . Two fundamentally important constraints regarding these quantities should be stated. First, in order to avoid reaching the Landau velocity at the band edge and creating unwanted excitations we require c l < c s , where c s is the sound speed. Thus we require t h a √ 3/2 < U n 0 /M or (t h a √ 3/2 )(4π 2 a sn /M ) −1/2 < 1. For 87 Rb with t h = × 10 3 Hz , a = 0.5 × 10 −7 m , a s = 5 × 10 −9 m ,n = 2 × 10 12 cm −3 , we get c l /c s 0.17. Second, in order for our long-wavelength approximation to be correct, we require the NLDE healing length ξ ≡ t h a √ 3/2 n 0 U ≫ a; using the same values for the physical parameters, we find ξ ≈ 9.43 a.
Next, we discuss the physical regimes for our theory. First, consider the extreme weakly interacting regime, i.e., n 0 U/t h ≪ 1. Then excitations of the ground state obey bosonic statistics. Furthermore, at length scales much larger than the healing length, ξk ≪ 1 where k is a characteristic quasi-particle momentum, excitations are comprised of correlated particle-hole pairs that propagate with a dispersion given by E ∝ k 1/2 . This is a Bose gas of composite particles in the sense that excitations of opposite spin are paired up (albeit non-locally) to form bosons. In contrast, for ξk ≫ 1, excitations are particle-like which corresponds to the case where spin eigenstates are excited independently. These states reflect the bipartite structure of the lattice (multi-component) but are local objects and so also reflect the bosonic nature of the fundamental constituent particles. They exhibit a mixture of fermionic and bosonic properties having a Dirac-like dispersion ∝ k but with quantum operators that obey bosonic statistics. This is a hybrid Dirac-Bose gas. On the other hand for moderate interactions, n 0 U/t h = 0 and close to the critical point on the superfluid side, higher order processes become significant and quasi-particles become highly nonlocal objects. Here quasi-particles are heavily "dressed" in the presence of the background condensate and so acquire a geometric phase resulting in anticommutation relations for their associated creation and annihilation operators (see previous section). In this regime, the Bogoliubov approximation fails to provide a good model but we expect that the Dirac spinor structure will remain robust against many-body effects with quasi-particles undergoing a renormalization due to many-body interactions similar to the case of graphene [18]. This describes an interacting Dirac gas.
Uniformly Moving Condensate. -Now we return to the RLSE and solve them for the simplest case of a uniform background Ψ(r) ≡ √ n 0 e iq·r (1, C 0 ) T , where C 0 ∈ C contains a relative phase, n 0 is the average particle density, and q is the condensate momentum measured with respect to the Dirac point. In order to obtain the coherence factors and quasi-particle dispersion, we must then solve a 4 × 4 eigenvalue problem; the RLSE yield In keeping with the usual Bogoliubov notation found in the literature, we may write E k = (c ′ l /c l ) q · ǫ 0 k ± E 0 k where ǫ 0 k ≡ c l k is the single quasi-particle energy for zero interaction and E 0 k = (ǫ 0 k ) 2 + n 0 U ǫ 0 k is the quasiparticle energy for a static background. The associated coherence factors can then be written as |u k,A(B) | =

The full interacting Hamiltonian is given byĤ
where A is the area of the plane. The first three terms are the mean-field and quantum corrections to the condensate energy and the last accounts for the number of quasiparticles present in the system. The constant c ′ l is defined in terms of the overlap integral between Wannier states at neighboring lattice sites by c ′ l = √ 3aτ /2 , where τ A,B ≡ − dr w * A ∇w B and τ = |τ A,B | 1 . The low energy behavior of the uniform condensate has a rich structure. The q = 0 case corresponds to a condensate with zero crystal momentum measured from the Dirac point but with momentum K relative to the lowest Bloch state of the crystal. The idea of a condensate in motion relative to its background has been treated in both free-space as well as the case of a moving background lattice [19,20]. Physically, the lattice potential is moving relative to the stationary condensate (laboratory frame). Two-body collisions reduce the momentum of some particles relative to the lattice (slowing down) and increase the momentum of others (speeding up) corresponding to a finite depletion of the condensate. In the laboratory frame, a two-particle collision appears as one particle gaining a component of momentum to the left and the other a component to the right. This is consistent with the well known particle-hole symmetry of the Dirac Hamiltonian: negative energy states can be interpreted as positive energy states that propagate in the opposite direction. In our theory these are quasi-particles with momentum K − k (for the K-Dirac point) relative to the lowest Bloch state.
For q = 0 then, we get E The two energy regimes evident here are separated by the condition c l k/n 0 U ≡ ξk ≈ 1. At short wavelength, kξ ≫ 1 so that E (±) k ≈ ±(c l k+n 0 U/2), where the dominant first term reflects only the presence of the honeycomb lattice, while the second term is a small mean-field Hartree shift due to the interaction with the background. When kξ ≪ 1, we find E (±) k ≈ ± k/ξ. These are collective excitations induced by the particle interactions just above the condensate energy. The presence of negative energy modes means that the condensate may lower its energy through spontaneous emission of radiation. This process can be suppressed by introducing an anisotropy in the lattice by breaking the A-B sublattice degeneracy with a deeper optical lattice in one direction [21]. This results in an additional term in the dispersion opening up a mass gap 2m eff at the Dirac point. For the negative energy modes we then have E citations require a minimum momentum determined by c l k min = 4m 2 eff + n 2 0 U 2 . Alternatively, we can consider the effect of the confining potential: this sets a lower bound for quasi-particle energy given by |E (−) k(min) | ∼ (c l 2π/R ⊥ ) 2 + n 0 U c l 2π/R ⊥ , where R ⊥ is the characteristic trap radius in the 2D plane.
Cherenkov Radiation. -This usually refers to the anisotropic emission of electromagnetic radiation from a source whose speed exceeds the local speed of light in some medium [22]. This concept generalizes to any source moving through a medium at a speed that exceeds the phase velocity of the elementary excitations of the medium. For example, a BEC moving in the laboratory frame, or with respect to a background, will "radiate" (emit particles) when its speed exceeds the sound speed. Moreover, the radiation will be emitted in a cone subtended by a specific angle in the direction opposite the motion of the BEC. The RLSE can be used to demonstrate this effect in the present context of a BEC in a honeycomb optical lattice.
For a BEC with momentum q > 0 measured from the Dirac point, examination of the angular dependence of E k reveals an intriguing structure for the emission of Cherenkov radiation. We observe the following properties for E k . (1) When v < c l , where v = c ′ l q is the condensate speed, all excitations have positive energy regardless of the angle of emission. (2) When v > c l , quasi-particle energies are positive only for emission angles (measured relative to q) for which θ < θ c ≡ cos −1 (−c l /v) while all other modes have negative energy corresponding to the emission of radiation in a backwards cone bounded by θ c . When v = c l , θ c = π marks the onset of radiation, in which case radiation is only emitted in the direction opposite q. This unique directional property of the radiation suggests an obvious detectable signature in the laboratory: a time-offlight analysis of a BEC prepared with precise values of the parameters should show a predictable shift in the momentum distributions between the forward and backward directions.
Nonlinear Localized Modes. -To obtain solutions of the NLDE that are localized in x, y for U > 0, we substitute the plane-polar ansatz sities. For the l = 1 case, we also obtain a skyrmion solution for which the pseudospin S =Ψ(r)σΨ(r) (with Pauli vector σ) exhibits an integral number of flips near the core and approaches a constant value far from the core. This feature is encoded in a topologically conserved charge (1/8π) Ω dr ǫ ij S · ∂ i S × ∂ j S which one recognizes as the Pontryagin index that classifies the mapping S 1 spin → S 1

∂Ω
where the two circles S 1 spin and S 1 ∂Ω parameterize the rotations between the densities ρ A(B) and the polar angle rotation on the 2D boundary ∂Ω at spatial infinity. In general, similar types of solutions exist for (ii) above. Analytical and numerical solutions are plotted in Fig. 2 for which different values of µ/U and l allow us to obtain the different asymptotic forms.
Besides vortices with integer phase winding, we also find solutions with fractional phase winding, called halfquantum vortices (HQVs). Ordinarily, analyticity (singlevaluedness) of the order parameter forbids the rotation of the phase of Ψ around the core to take on fractional values. In the NLDE, Ψ can acquire a coherent internal Berry phase in addition to an external phase whose angles are identified with the polar angle θ [23,24]. Such states may have half-integer winding in both the internal and external phase angles while remaining singlevalued overall. We obtain HQVs with asymptotic form lim r→∞ Ψ HQV (r) = 2i n 0 /2 e −iθ/2 [cos(θ/2), i sin(θ/2)] T ; the complete solution is shown in Fig. 2 (c).
We also obtain one-dimensional kink-soliton, skyrmion, and line-soliton solutions. The kink and skyrmion solutions are obtained by a straight-forward substitution of the ansatz Ψ(x) = η [cos(ϕ) , sin(ϕ)] T into the NLDE and then considering the distinct cases where ϕ = constant p-5 (kink) or η = constant (skyrmion). The line-soliton solution is obtained when both η and ϕ are functions of x with the additional condition that, at the origin, η remains below a certain value. This ensures that µ 2 < U/8 < (U + 1)/2U ⇒ U < 3.365 and µ 2 < 0.649, which allows the wavefunction to collapse away from the y-axis while the nonzero wavefunction near and along the y-axis has a Lorentzian form in the x-direction due to the attractive effect from the kinetic terms. Localized Mode Stability. -To check the stability of our localized solutions we substitute the corresponding solution into the RLSE and obtain the low energy spectrum. For the case of the vortex/soliton (l = 0), we solve the eigenvalue problem numerically and obtain E 0 = −3.9274+0.0020 i in units of n 0 U . Similar to the case of a trapped rotating BEC [25], we find that the lowest mode in the spectrum is anomalous with negative energy and positive norm. The coherence factors shown in Fig. 3, near the core of the vortex, appear as |u k,A(B) | 2 ∼ 10 −2 , |v k,A(B) | 2 ∼ 10 −5 so that u k,A(B) >> v k,A(B) . This mode is dynamically unstable due to the presence of a small imaginary component which gives a decay rate relative to characteristic oscillation time of 0.0020/3.9274 = 0.0005. The lowest quasi-particle energies for the other localized solutions are: −3.9276 + 0.0019 i; 2.634 × 10 2 + 9.96 × 10 4 i; −3.9274 + 0.0019 i; 7.8409 × 10 −3 − 9.9993 × 10 2 i; 7.9349 × 10 −3 − 9.9993 × 10 2 i; for the ring-vortex/soliton, halfquantum vortex, planar skyrmion, line-skyrmion, and linesoliton respectively, where all quantities are given in units of n 0 U .
In conclusion, we have shown that an effective relativistic fermionic system may be designed using ordinary cold bosonic atoms as the underlying degrees of freedom. We solved the resulting NLDE for different classes of nonlinear modes including half-quantum vortices. We derived and solved relativistic linear stability equations and gave explicit criteria for experimental observation of Cherenkov radiation, as well as predicting an anomalous mode for the vortex/soliton solution. Density profiles may be observed by time-of-flight techniques to detect both massive and massless Dirac fermions in the laboratory [21,26]; nonlinear modes involving phase winding can be created by techniques analogous to those used at JILA [27]; and we anticipate that Bragg scattering can be used to populate the Dirac cones at both K and K ′ points, leading to arbitrary superpositions over our localized solution types between the two cones, and thereby populating all four components of the Dirac spinor.
We thank Ken O'Hara, Jeff Steinhauer, and Michael Wall for useful discussions. This work was supported by the NSF and the Aspen Center for Physics.