The nonlinear Dirac equation in Bose-Einstein condensates: II. Relativistic soliton stability analysis

The nonlinear Dirac equation for Bose-Einstein condensates in honeycomb optical lattices gives rise to relativistic multi-component bright and dark soliton solutions. Using the relativistic linear stability equations, the relativistic generalization of the Boguliubov-de Gennes equations, we compute soliton lifetimes against quantum fluctuations and classify the different excitation types. For a Bose-Einstein condensate of $^{87}\mathrm{Rb}$ atoms, we find that our soliton solutions are stable on time scales relevant to experiments. Excitations in the bulk region far from the core of a soliton and bound states in the core are classified as either spin waves or as a Nambu-Goldstone mode. Thus, solitons are topologically distinct pseudospin-$1/2$ domain walls between polarized regions of $S_z = \pm 1/2$. Numerical analysis in the presence of a harmonic trap potential reveals a discrete spectrum reflecting the number of bright soliton peaks or dark soliton notches in the condensate background. For each quantized mode the chemical potential versus nonlinearity exhibits two distinct power law regimes corresponding to the free-particle (weakly nonlinear) and soliton (strongly nonlinear) limits.


Introduction
Vacuum states with broken symmetry play an important role in the study of quantum many-body physics, since they provide clues to the principles that govern the full symmetric theory [1,2,3]. Solitons are finite energy solutions of classical equations of motion and have been studied as nonuniform ground states, i.e., bound states or defects in the fundamental degrees of freedom that provide a launching point for perturbative expansions. Broken translational, rotational, or inversion symmetry, ubiquitous to discrete as well as continuous systems, can usually be cast in terms of a topological framework [4]. When attractive interactions are present non-topological solitons model globally regular bound states of the system [5,6]. Such states owe their existence to an unbroken symmetry of the Lagrangian and thus have a conserved Noether charge. In contrast, topological solitons are defects typically associated with spontaneous symmetry breaking. In this case the defect breaks a discrete symmetry and appears as a boundary separating two degenerate asymptotically flat solutions while retaining a topological charge degree of freedom as a relic of the broken symmetry. Examples of solitons in extant physical systems include domain walls in BCS superconductors [7], superfluid vortices [8,9,10], and quantum Hall states in topological insulators [11,12]. In one spatial dimension dark [13,14,15] and bright [16,17,18,19] solitons in repulsive or attractive Bose-Einstein condensates (BEC) with spontaneously broken U(1) symmetry are examples of broken spatial symmetry. Beyond familiar condensed matter systems solitons emerge in low-energy sectors of the standard model of particle physics as extended particles [20,21,22], and in M-theory as subcritical dimensional D-brane embeddings [23].
In all of these cases, one is typically interested in the properties of the low-energy spectrum since this characterizes the system near equilibrium. The presence of a defect, or soliton, partitions the domain into a core region which spans the size of the defect, and a bulk region far from the core. Excitations in the bulk describe the system's response to the presence of the soliton, whereas fluctuations in the core describe undulations and translations of the soliton itself. In superfluid systems, soliton core bound states may be metastable, possessing a finite lifetime against dissipation through lower energy scattering states, or truly stable if the soliton lies at an energy minimum of the system.
In this article we focus on elementary excitations and stability of a topological defect near the Dirac point of a BEC. At very low temperatures interactions between condensate and non-condensate atoms is minimal, allowing for existence of longlived metastable states. Thus, quantum fluctuations of a kink-like soliton in the nonlinear Dirac equation (NLDE) presents an analog of a domain wall in a gas of Dirac fermions interacting through a local quartic term [24]. Solution profiles for the soliton backgrounds were explored analytically and numerically in a companion paper [25]. We note that similar solitons appear in nonlinear optics [26,27], in graphene [28,29,30,31,32,33], and in various other fields of physics [34,35,36,37,38,39,40,41,42,43,44]. Figure 1 provides a schematic of our setup depicting solitons and their fluctuations in the two distinct quasi-one-dimensional (quasi-1D) reductions of the honeycomb lattice. The solution of the relativistic linear stability equations (RLSE) gives us the linear spectrum from the presence of small quantum fluctuations in the BEC [45,46]. For both dark and bright soliton solutions, far from soliton core the BEC occupies only one of the two sublattices, switching from the A sublattice to the B sublattice when translating through the core. Thus these solution types present a 1D analog of a skyrmion localized to the soliton width. However, the skyrmion analogy does not hold near the core since in our case the total density ρ(x) is a non-constant function of the longitudinal coordinate x. The total density here is defined as the sum of the squared spinor amplitudes, which in the case of the reduced two-spinor formulation is ρ(x) ≡ |ψ A (x)| 2 + |ψ B (x)| 2 , with ψ A (x) and ψ B (x) the wavefunctions corresponding to A and B sublattices of the honeycomb lattice. We will show that quasi-particle excitations far from the soliton exist as scattering states which respect this asymmetry. Because of this feature, it is convenient to think of the switching point from the A to B sublattice as a defect analogous to a domain wall. It is instructive to view the NLDE from a mathematically elegant perspective by recasting it in terms of the covariant pseudospin formalism. As we will see, this approach allows for a domain wall interpretation which connects to other areas of physics. For example, in magnetic systems domain walls appear as topologically stable solitons separating two distinct regions of different magnetic polarization [47]. Another context is the case of two interpenetrating BECs comprised of atoms in different hyperfine states, wherein one finds regions across which the relative phase of the two condensates changes by 2π [48]. In spin-1 BECs, domain walls have been studied extensively as boundaries between regions of pseudospin polarization S z = ±1 [49,50], in addition to investigations into the quasi-particle transmission and reflection properties of such boundaries [24].
Domain walls also play an important role in high energy physics, for example as extended supersymmetric objects which isolate different vacua [51,52]. It is thus not surprising that solitonic objects play an important role in both condensed matter and particle physics settings. A particular example which highlights this fact is the recent simulation of tachyon condensation using two-component BECs [53]. In such analogs one finds that spontaneous symmetry breaking occurs in a two-dimensional subspace of the full system, i.e., a domain wall in the larger space. In each of the examples mentioned here the domain wall is identified with a continuous deformation of the order parameter between two degenerate asymptotically flat states of the system. The key feature of the deformation is that it is localized; it occurs over a finite region in at least one of the spatial dimensions.
This article is organized as follows. Section 2 establishes the full symmetry of the quasi-1D NLDE order parameter manifold. This describes the set of possible order parameters determined by a series of symmetry breaking reductions from the full (3+1)-dimensional Poincaré group. In Sec. 3, we solve the RLSE numerically to determine soliton lifetimes. In Sec. 4, we solve the RLSE analytically through a method of decoupling and derive the phase and density fluctuations in the soliton core region, which include the Nambu-Goldstone mode responsible for U(1) symmetry breaking, i.e., Bose condensation. In Sec. 5, we solve for the continuous spectrum far from the soliton core where we find the Nambu-Goldstone mode and a spin wave, the later corresponding to nonzero density fluctuations. The asymptotic spectrum naturally leads into Sec. 6 where we formulate relativistic solitons in the language of spin-1/2 domain walls. In Sec. 7, we analyze quantum fluctuations in light of the domain wall interpretation. In Sec. 8, we treat the spectral theory of a BEC in a weak harmonic trap. Finally, in Sec. 9 we conclude.

Symmetries of the order parameter manifold
The order parameter that we study is analogous to metastable vacua in high energy systems with quasi-particles and thermal excitations playing the role of virtual and real particles, respectively. Clarifying the underlying symmetries of the order parameter manifold is key towards identifying the various excitations associated with continuous symmetry breaking. In the quasi-1D NLDE [25], the order parameter manifold comes from a series of symmetry breaking steps. To see this, we begin by noting that noninteracting bosons at the Dirac point of a quasi-2D honeycomb lattice occupy singleparticle states in one-to-one correspondence with massless Dirac states. The 2 × 2 unit and Pauli matrices 1, σ x , σ y are the group generators in 2D consistent with the spin and momentum vector coupled Dirac Hamiltonian H p = c l σ · p. One may think of the absence of the third Pauli matrix σ z a consequence of projecting the full SU(2) group onto the coordinate plane thereby removing one degree of freedom through the reduction SU(2) → U(1) ⊗ Spin (2). Here the factor of U(1) accounts for an overall phase and the spin group Spin(2) is isomorphic to a double covering of U(1), i.e., expressed in terms of the fundamental group π 1 (Spin(2)) ∼ = 2Z ∼ = 2π 1 (U(1)). This can be summarized in a short exact sequence by recalling the isomorphisms U(1) ∼ = SO(2), 2Z ∼ = Z 2 , then from which we write Spin(2)/Z 2 ∼ = SO(2), or equivalently Spin(2) ∼ = Z 2 ⊗ SO(2). The quasi-1D theory then demands a second coordinate reduction which breaks the 2D rotation group into its reflection subgroup along either of two orthogonal directions SO(2) → Z 2 ⊕ Z 2 , where the two copies of Z 2 are the reflection subgroups associated with two orthogonal directions in the plane connected with armchair and zigzag forms of the NLDE [25]. From this we see that the full symmetry of the quasi-1D NLDE order parameter manifold is where K(1+1) is the set of relativistic boosts and space-time translations in (1+1)D. To make this discussion more concrete, we can write the representation of this symmetry reduction from 2D to 1D in terms of the order parameter manifold as where θ(p) ≡ tan −1 (p y /p x ), and in terms of the Hilbert space the reduction in Eq. (3) acts according to H 2D → H zigzag H armchair . Note that on the left side of Eq. (3) the vector p is two-dimensional, whereas the right hand side applies to one spatial dimension. In the reduced space, the direction of p is completely determined by a sign, i.e., p = ±|p| ≡ ± p. We adhere to this convention throughout our work. The first order parameter manifold in Eq. (3) has the full U(1) ⊗ Spin (2), where φ and θ are the U(1) and Spin(2) parameters. To the right of the arrow in Eq. (3) the order parameter takes on the reduced symmetry where φ is the U(1) parameter with the positive/negative eigenvalues and parity inverse factors associated with the Z 2 ⊗ Z 2 products in Eq. (2).
The presence of a soliton in the reduced quasi-1D problem breaks translational symmetry, in which case one would expect to find one zero-energy mode in addition to one massless excitation for each broken continuous symmetry. These include two Goldstone modes, one from condensation in the overall phase and one from the internal phase; and two zero modes, one from breaking rotational symmetry when going from 2D to 1D, and one from the broken translational symmetry due to the soliton. Only two out of the four are in fact present. The Goldstone and zero modes from breaking Spin(2) symmetry are suppressed, since they fluctuate along the direction of the quasi-1D confining potential. We expect therefore to find one Goldstone mode as an overall phase fluctuation and a zero mode from the soliton. It must be noted that in the literature the Goldstone mode is sometimes identified as a zero mode. Technically, the Goldstone mode corresponds to the gapless energetic branch associated with local twists in the phase. When the condensate background is spatially uniform, the Goldstone branch is continuous and connects to a spatially uniform zero mode. In the presence of a defect, however, translational symmetry is broken and the Goldstone branch is discrete with a nonzero momentum lower bound, p ≥ p min . In this case the Goldstone branch connects to a spatially nontrivial zero mode in the limit p → p min .
The nonlinearity in the NLDE allows for asymptotically flat solutions |ψ A |, |ψ B | ∼ 0, µ/U , for |x| much larger than the soliton core size. The Z 2 ⊗ Z 2 symmetry in the zigzag and armchair states of Eq. (3) leads to four distinct asymptotic states but only two are independent because of an overall phase redundancy. These are for the zigzag solution, and for the armchair case, an overall complex constant prefactor is omitted for clarity. As we showed in [25], NLDE solitons interpolate between two asymptotic states that are linear combinations of for the zigzag case, and for the armchair case. We will see in this article that the presence of a soliton partially breaks the inversion symmetry implicit in Eqs. (6)-(7), splitting the spectrum into massless modes with linear dispersion, which retain the full symmetry, and massive modes with quartic dispersion, which break parity inversion symmetry. The central focus of this article is to understand the nature of these quantum fluctuations, both asymptotically and in the transition region inside the soliton core.

Stability of soliton solutions
The combination of the honeycomb lattice geometry and the atom-atom interactions results in a characteristic signature effect on soliton stabilities. In particular, the presence of negative energy states below the Dirac point means that a BEC will eventually decay by radiating into the continuum of negative energy scattering states. However, this requires a mechanism for energy dissipation into non-condensate modes which must come about from secondary interactions with thermal atoms. Thus, as long as the system is at very low temperatures our main concern for depletion of the BEC comes from potential imaginary eigenvalues in the linear spectrum. The situation is analogous to dark solitons in quasi-1D BECs described by the nonlinear Schrödinger equation: in practice such excited states can easily have a lifetime longer than that of the BEC [54]. In this section, we compute the linear spectrum for soliton solutions of the zigzag and armchair NLDE. Before proceeding it is useful to elaborate on units and dimensions of some of the physical quantitates key to our discussion. The main composite parameters relevant to the NLDE, and hence the RLSE, are the effective speed of light c l = t h a √ 3/2 and the quasi-1D renormalized atom-atom binary interaction strength U 1D = U 2D /(π 1/2 L y ), expressed in terms of its quasi-2D counterpart U 2D = L z gn 2 3 √ 3a 2 /8. The presence here of the trap oscillator lengths, L y and L z , reflect the fact that U 1D and U 2D come from integrating over the degrees of freedom transverse to the single large dimension in our problem. For instance, U 1D is obtained by integrating over the ground state in the y-direction in the quasi-2D NLDE [46] where the oscillator length is related to the frequency ω y and atomic mass M by L y = /M ω y . The parameters that comprise U 2D and c l are the vertical oscillator length L z (in the quasi-2D problem), the average particle densityn = N/V , the interaction g = 4π 2 a s /M , the lattice constant a, and the hopping energy t h . Throughout our work we take the atomic mass M and scattering length a s = 5.77 nm to be those of 87 Rb. A complete discussion of NLDE parameters and constraints can be found in [46]. With these parameter definitions one finds that the spinor order parameter Ψ = (ψ A , ψ B ) is dimensionless and the quasi-1D interaction strength U 1D has dimensions of energy. To simplify the notation, from here on we will omit the subscript on U 1D and write U for the quasi-1D interaction strength.
To compute soliton lifetimes we must solve the relativistic linear stability equations (RLSE) modified for our quasi-one-dimensional problem [55]. This allows us to account for quantum mechanical perturbations to the mean-field result by using the corrected order parameter with the condensate spinor wavefunction and quantum correction given by whereα † andβ † (α andβ) are the creation (destruction) quasi-particle operators and u A(B) and v A(B) are the associated spatial functions, respectively. Linear stability of a particular soliton solution is determined by substituting the spatial function for that solution (i.e., the dark or bright soliton) into the RLSE as a background for the quasiparticle functions. This substitution gives a set of first-order coupled ODEs in one independent variable to be solved consistently for the quasi-particle energies E k and amplitudes u k and v k , where the subscript denotes the mode with momentum p = |k|, defined in terms of the magnitude of the wavevector k. We remind the reader that we are working in one spatial dimension, thus there is at most a sign difference between the bold vector notation and the corresponding norm: k = ±|k|. Since we are perturbing from a spin- have vector form describing quasi-particle and quasi-hole excitations of the A and B sublattice, as indicated by the A(B) sublattice subscripts. We discretize the derivatives and spatial functions in the RLSE using a forward-backward average finite-difference scheme, then solve the resulting discrete matrix eigenvalue problem using the Matlab function eig.
Solutions of the RLSE are perturbations of the NLDE four-spinor components and respect the same decoupling to two-spinor form. Thus, focusing on equations for the upper two-spinor, the 1D RLSE is Equations (12)-(15) inherit the linear derivative structure on the sublattice particle and hole functions u A(B) and v A(B) . The constant chemical potential µ and particle interaction U appear as coefficients in addition to the spatially dependent condensate profiles ψ A(B) (x) and eigenvalues E k . The parameters in Eqs. (12)-(15) are already renormalized due to dimensional reduction from 2D to quasi-1D as described in Sec. 3 of Ref. [25], specializing to the zigzag or armchair geometries. We point out that Eqs. (12)- (15) pertain to the zigzag NLDE. In the armchair version the momentum terms have identical complex coefficients, −i c l , which comes from rotating the Dirac operator by 90 degrees. This transformation between zigzag and armchair forms is equivalent to the two-spinor Pauli transformation discussed in Sec. 2 in Ref. [25], and the four equations of the RLSE inherit this feature: choosing to work in one form leads to no loss of generality. Alternatively, one may argue that since the RLSE are linear in the amplitudes u A(B) and v A(B) , absorbing a factor of i into either pair of the sublattice amplitudes, i.e., u A and v A or u B and v B , simply converts between the zigzag and armchair forms. Thus, for a given condensate spatial profile the RLSE for the zigzag and armchair geometries have the same linear eigenvalues, and the stability properties of solitons in both cases are the same. We find the lowest excitation energies for the two types of solitons E DS  eigenvalues correspond to modes which decrease the energy of the solitons into states below the Dirac point. What is significant is the absence of imaginary modes; thus our solitons are dynamically stable in both zigzag and armchair geometries. This means that at very low temperatures we expect solitons to remain viable over the lifetime of the BEC. To obtain the next order correction due to finite temperature effects would require a modified version for the RLSE analogous to the Hartree-Fock-Bogoliubov treatment which takes into account interactions between condensate and non-condensate atoms [56].

Bound state fluctuations of the soliton core
We would like to solve for the quasi-particle structure of the NLDE using analytical methods. Towards this end, in this section we reduce the RLSE down to four decoupled second order equations. We begin by changing variables using symmetric and antisymmetric functions defined as where we have suppressed the mode index k in order to simplify the notation. For the zigzag geometry and using the transformation defined by Eqs.
which we obtain by adding and subtracting the transformed versions of Eqs. (12)- (13) and Eqs. (14)- (15). Note that we use the prime notation to denote derivatives with respect to x. These equations lead to a further decoupling into two pairs of equations.
Transforming to the dimensionless form using the new variable and constantsμ = µ/U , where we have defined new variables Note that now Eqs. (23)-(24) are not coupled to Eqs. (25)- (26). Next, we make the substitution Substituting Eqs. (27)-(28) into Eqs. (23)-(24), we obtaiñ Back substitution yields the decoupled second order equations These equations may be further simplified by transforming to standard form using the substitutionũ Eq. (23)-(24) are reduced to where Following a similar path yields the reduced forms for Eqs. (25)-(26) where Equations (35)-(36) and (39)-(40) comprise the final reduced form of the RLSE in terms of four decoupled second order equations. The quasi-particle amplitudes can be obtained from the functions w, z, q and s by working backwards through each transformation step by which we obtain It is also interesting to note that Eqs. (35)- (36) and (39)-(40) have a Schrödinger-like form for a particle with zero total energy subject to the potential functions Q, R, S and T . For soliton solutions the potential functions develop minima which support stable localized bound states. It is also important to note that Eqs. and where we have used the fact that the square modulus of a quasi-particle amplitude is much smaller than that of the condensate wavefunction. Next, we solve the decoupled equations (35)-(36), (39)-(40) using approximate forms for the dark soliton spinor components The approximate forms in Eqs. gives the quasi-particle amplitudes The envelope functions in these expressions decay exponentially for x → ±∞, with the specific details of the fluctuations determined by solving for the functions w, z, q, and s subject to the particular forms of the potentials Q, R, S, T . Computing these potentials analytically, we find three regimes: (i) all potentials are negative on the interval −∞ < x < +∞ when |Ẽ| < 2 − √ 3, which lead to exponentially growing and decaying solutions; (ii) two of the potentials, Q and R, are identically zero and two are negative on the interval −∞ < x < +∞ when |Ẽ| = 2 − √ 3; and (iii) potentials Q and R become positive when 2− √ 3 < |Ẽ| <μ, for which complex oscillating solutions exist. The upper boundμ is the energy of the gapped branch of the continuous spectrum at zero quasi-particle momentum (p → 0) which we will discuss in Sec. 5.
We can interpret the regimes (i) and (ii) in light of the density and phase fluctuations, Eqs. (47)- (48). In regime (i), the fluctuations reduce to and A bound state zero mode (E = 0) is supported associated with the breaking of translational symmetry by the soliton [57,58]. In this case, we find that Q(x) = R(x) and S(x) = T (x) over the infinite interval −∞ <x < +∞ so that w = z and s = q. Thus, in the standard normalization the zero mode has a zero norm, [59]. Moreover, very near the core where |x|/ξ 1, where ξ = c l /U is the healing length which gives the approximate width of the core, the potential functions satisfy the relations Q ≈ T + O(x/ξ), R ≈ S + O(x/ξ), so that w ≈ s, z ≈ q. In this case, the total combined density and phase fluctuations from both sublattices are Thus, near the soliton core the quasi-particle mode contributes to an overall phase fluctuation but does not contribute to density fluctuations for which the A and B sublattice contributions cancel exactly. This is the Nambu-Goldstone mode associated with simultaneous U(1) symmetry breaking in both A and B sublattices. In regime (ii) the situation is slightly different since here the functions w, z, q, s oscillate, but the general results of our analysis hold. Physically, fluctuations of the soliton core correspond to an additional quantum uncertainty in a single measurement of value of the phase at the core and to a nonzero average phase over large time scales which imparts a net translational motion to the soliton.

Continuous spectrum far from the soliton core
Equations (35) which yields the spectrum The long wavelength limit of Eq. (58) defined by c l p/µ 1, gives the four branches of the continuous spectrum where the two modes E ± g have a gap equal to the condensate chemical potential µ, and the modes E ± 0 are gapless linear Dirac-like excitations for p → 0. The continuous spectrum Eq. (58) is plotted in Fig. 3. We point out that the results in this section apply to the dark and bright solitons alike since both share the same asymptotic form.

Relativistic solitons as domain walls
In the case of NLDE solitons the asymptotic values of the spin components correspond to two different asymptotically flat vacua connected through the core region where a spin twist occurs. In this section we study the internal spin space rotations that occur inside the core regions of NLDE solitons, and thus firmly establish such solutions as domain walls. We work in the chiral representation since this is the most natural approach as particle interactions do not mix upper and lower Dirac two-spinors. For the two-spinor order parameter ψ(x, t) = [ψ A (x, t), ψ B (x, t)] T , the quasi-1D reduction of the NLDE is expressed concisely in terms of spin-orbit coupling and contact interaction terms where the interactions are contained in the 2 × 2 diagonal matrix N defined as and the commonly used "spin-orbit" terminology refers to the derivative terms expressed as a contraction between the Pauli matrices σ µ and the space-time derivatives ∂ µ . The contraction is over the Greek index µ = 0, 1 following the usual Einstein convention σ µ ∂ µ = g µν σ µ ∂ ν , with the (1 + 1)-dimensional Minkowski metric and Pauli matrices The additional Pauli matrix σ 1 is chosen for the specific case of the zigzag or armchair geometry, A thorough description of this formalism may be found in Ref. [55].
To make the connection between superfluidity and relativistic current apparent in the pseudospin formalism, we use the Gordon decomposition approach [60] which lets us separate the phase gradient from the magnitude gradient in Eq. (61). The former describes conventional superfluidity and the latter contains the pseudospin dependent contribution to the overall current. In our previous work [55] we found that the NLDE has the associated current whereψ ≡ ψ † σ 0 . Solving Eq. (61) for ψ and the conjugate of Eq. (61) forψ allows us to express Eq. (65) as where N −1 is the inverse matrix of Eq. (62). Working out the properties of the matrix products σ µ N −1 σ ν , we find the decomposition where D αν and O µν are 2 × 2 diagonal and off-diagonal matrices with respect to the spin indices, respectively, and encapsulate the interaction terms. The current in Eq. (66) can now be expressed as where we use the ∂ ± ν notation to mean the antisymmetric differentiationψ M µν ∂ ± ν ψ ≡ ∂ νψ M µν ψ −ψM µν (∂ µ ψ). The first term in Eq. (68) is the orbital contribution to the current which does not mix pseuodspin indices and describes the familiar superfluid velocity as the gradient of an overall phase, whereas the second term is the pseudospin portion which mixes spinor components and describes the phase-independent density current.
For solutions of the zigzag and armchair NLDE, Eq. (68) gives the following results: (i) For the zigzag NLDE with the spinor solution Ψ( where f A and f B are real functions of x and solve the time-independent armchair NLDE, we have (ii) For the armchair NLDE the relative phase between spinor components is zero and Thus, for the ordinary dark soliton and bright solitons found in Ref. [25], the currents are peaked at the soliton cores with the current j x pspin for the zigzag solitons pointing in the negative x direction and in the positive direction for armchair solitons. To understand results (i) and (ii) we emphasize that NLDE solitons have a constant (spatially uniform) overall phase which reflects the zero superfluid current j µ orb = 0, whereas the relative phase between spinor components varies through the soliton core leading to nonzero values for j t pspin and j x pspin . The structure of NLDE solitons is reminiscent of domain walls as in the interface between 3 He-A and 3 He-B superconductors [61] but more closely resembling the boundary between pseudospin domains in two-component BECs. The analogy extends only loosely to spin-1 BECs [24] in the case where the asymptotic order parameter (far from the domain wall, i.e., x → ±∞) is in oppositely polarized states. However, we must point out that two-dimensional BECs in honeycomb lattices are Dirac spin-1/2 analogs which couple spinor components through the kinetic term; in contrast, multi-component BECs are coupled through the interaction term.
Our system possesses a pseudospin structure described by the vector operator S = (S x , S y , S z ), where the component spin operators are expressed in terms of the Pauli matrices S x,y,z = ( /2) σ x,y,z , and the spin quantization axis is in the z-direction. The appearance of the Pauli matrices here underscores the marked contrast to ordinary twocomponent BECs: the latter are not associated with a Clifford algebra, i.e., fermionic anticommutativity. Although strictly speaking there is no z-component of spin in our problem and we must choose between σ x or σ y for the 1D zigzag or armchair NLDE, the transverse spin operator σ z is still relevant in characterizing our solutions.

Spin waves and Nambu-Goldstone modes
In this section we will extend the discussion from Sec. 6 to the asymptotic linear spectrum far from the soliton core. As discussed in Eqs. (9)-(11), in general the quantum perturbed wave function can be written as Ψ(x) +φ( T , from which the fluctuation in the spin density δF = (δF x , δF y , δF z ) T is computed to be Note that we have used the quasiparticle averages φ A(B) = φ A(B) and condensed the notation by not writing out the space-time dependence of the component fluctuations.
In the case of a spin-1 BEC there are three types of spin fluctuations: transverse, quadrupolar, and Nambu-Goldstone. The first mixes the S z = ±1 into the S z = 0 hyperfine channel, the second mixes the S z = +1 and −1 hyperfine states, and the third is the phase fluctuation associated with U (1) symmetry breaking (see Ref. [24]).  In our spin-1/2 problem, the physical picture is slightly different. Here there are only two modes: a spin wave mode which mixes spin-up and spin-down states and the Nambu-Goldstone mode. In the zigzag and armchair geometries, the spin wave modes are, respectively, δF SW = δF x and δF SW = δF y ; in both cases the Nambu-Goldstone mode is δF NG = δF z . Near the soliton core, the spin wave mode vanishes as we showed in Sec. 4, leaving only the Nambu-Goldstone mode. Far from the core the situation is different. For instance, in the zigzag geometry where x/ξ → −∞, we have F z = + /2, In contrast, in the other limit . We may verify that our definitions for δF SW and δF NG make sense physically by computing their momentum dependence. Expanding the quasi-particle functions gives , and solving the RLSE in the long wavelength limit (c l p/µ 1), shows that the gapless excitation, with energy The formulation of fluctuations in terms of spin wave modes sets the ground work for detailed analysis of scattering from the soliton core, an interesting topic in itself which may provide insight into the question of integrability of the quasi-1D reduction of the NLDE. The pseudospin domain structure which we have discussed in this section is summarized in Fig. 4.

Discrete spectra for solitons in a harmonic trap
In experiments, solitons reside in a BEC within a harmonic magnetic trap. Consequently, in the case of extended solitons the trap boundary affects the soliton in a nontrivial way. The result is quantization of spatial modes, which has a significant effect at large healing length (comparable to the trap size) or equivalently for weak nonlinearity. In this section we study the behavior of our dark and bright solitons in the presence of a harmonic trap by computing the chemical potential spectra for the single, double, and triple soliton states. As throughout this Article, these results apply to both the zigzag and armchair geometries.
For the case of a highly oblate harmonic confining potential which defines a 2D system, the oscillator frequencies satisfy ω z ω ≡ ω x , ω y . If in addition to this condition we also take ω y ω x , as we discussed in [25] with the soliton in the xdirection, we encounter only the effects resulting from quantization of spatial modes along the soliton direction in either the zigzag or armchair geometries shown in Fig. 1. We then take the longitudinal trapping potential to be V (x) = 1 2 M ω 2 x 2 . We proceed numerically by incorporating this potential into the NLDE and then transforming to dimensionless equations by defining thereby obtaining the dimensionless form of the NLDE Here the two rescaled physical parameters in the NLDE are The analogous rescaling for the case of the NLSE in an oblate harmonic trap differs from our problem in a fundamental way. For the NLSE, energies are scaled to the trap energy ω and lengths to the oscillator length = /M ω (see Ref. [62]). In contrast, our problem retains the same scaling to the trap energy but lengths are scaled to the ratio c l /ω = t h a √ 3/2 ω as can be seen in Eq. (81), where the natural scales of the lattice appear in the hopping energy t h and lattice constant a. This particular choice of scaling is forced on us because of the single spatial derivative in the NLDE; we cannot completely scale away the lattice information. The mass energy factor M c 2 l in Eq. (84) is thus a direct result of the relativistic linear dispersion of the NLDE. Length scales are defined in Table 1 along with associated momentum and energy scales for quasi-1D Physical scale Length Momentum Energy Table 1. Physical Scales. Length, momentum, and energy scales for the quasi-1D NLDE in a harmonic trap. Scales are determined by the 3D healing length ξ; the transverse oscillator length ⊥ ; the large-momentum healing length µ ; the scale latt associated with the lattice constant and hopping energy; the low-momentum quasi-1D healing length ξ 1D ; and the harmonic trap length trap . All other fundamental parameters were defined in Sec. 3. Momentum and energy scales are related to their associated length scales by: momentum ∼ /length and energy ∼ 2 /(M × length 2 ). We have included the transverse oscillator frequency ω ⊥ which defines the transverse size of the condensate.
NLDE solitons in a harmonic trap. Note that the last three scales in Table 1 are key in our calculations since they contain information about the first two scales. Realization of the NLDE solitons in a harmonic trap requires the particular ordering of length scales due to the long-wavelength approximation used to obtain the NLDE. † The lengths in the hierarchy Eq. (85) are: the lattice scale latt , which contains the lattice constant and hopping energy; the quasi-1D effective healing length ξ 1D , which incorporates the atom-atom interaction and the transverse length; and the harmonic trap length trap , which defines the overall size of the BEC. For a typical scenario for a 87 Rb BEC [25] with a trap oscillator frequency ω = 2π × 0.039 Hz, we obtain latt ≈ 2.3 µm, ξ 1D ≈ 10 µm, and trap ≈ 55 µm.
To connect with dark and bright soliton solutions of the NLDE [25] we consider the limit for zero trap energy. This amounts to taking the trap size to infinity, i.e, letting trap to be the largest length scale in Table 1. Expressed in terms of the lengths in Table 1 one finds that the derivative, interaction, and chemical potential terms in Eqs. (82)-(83) scale as trap / latt , 2 trap / latt ξ 1D , and 2 trap / µ latt , respectively, while the harmonic term does not scale with the trap size. Thus in the large trap limit the harmonic term can be neglected and we regain the continuum theory as expected.
We use a numerical shooting method to solve the NLDE in the presence of the harmonic trap [25]. This is done by first expanding the spinor wavefunction Ψ(x) = [ψ A (x), ψ B (x)] T in a power series about the center of the trap at x = 0. The leading coefficient a 0 in the expansion for ψ A is then tuned to obtain a stable solution. The second free parameter b 0 , from ψ B , is held fixed between 0 and 1, for the dark soliton, or between 1 and the value at the peak, for the bright soliton. For b 0 = 0,    the single soliton in order to push excitations out to larger values of x. See also [63] for a study of precision issues in shooting methods related to BEC in harmonic traps. Next, we solve Eqs. (82)-(83) with a nonzero oscillator length. In the presence of the trap potential solutions are spatially quantized and labeled by a discrete index. In particular, for Q = 10 3 , corresponding to a longitudinal oscillator frequency ω = 2π × 0.0305 Hz, we find the free parameter a 0 for the single dark soliton at a 0 = 0.94640402384±10 −9 and for the two and three soliton states a 0 = 0.89882708125±10 −9 and a 0 = 0.8523151 ± 10 −13 , respectively. The lowest multiple dark solitons in a trap are plotted in Fig. 7 Near the origin, the trap potential is weak and the chemical potential term dominates so that we have η B < 0 and η A > 0. However, as we move away from the origin and into the strong potential region the quadratic terms in χ grow eventually overwhelming (e) the other terms. In this asymptotic region η A and η B solve the limiting equations whose solutions are These functions oscillate with a spatially increasing frequency k ≡ Q χ 2 , so it is clear that the tail oscillations are coming from the unbounded potential barrier. Physically, the barrier potential forces a positive energy particle into the continuum of negative energy states below the Dirac point. In contrast, this effect does not arise for an ordinary Schrödinger-like particle in a quasi-1D harmonic potential. There the particle is described by a single component wavefunction which must decay exponentially inside the potential barrier. Nevertheless, in terms of the density the oscillations average to zero. This phenomenon is known as Zitterbewegung and is associated with relativistic fermions [64,65]. To obtain the functional relation between the chemical potential µ and the interaction U for a particular excitation inside the harmonic trap, we first derive an expression for the normalization of the wavefunction for the new rescaled NLDE in Eqs. (82)-(83), which is found to be where the right hand side is given by where N is the number of atoms in the system and we have formulated the expression after the second equality in terms of the lattice constant a and effective speed of light c l . To compute the chemical potential spectra, we fix Q (which is the same as fixing the relative effects of the lattice geometry and the trap) and varyμ, calculating the norm N for each value ofμ. We thus obtain paired values (N ,μ). These values for the single, double, and triple soliton states are shown in Fig. 9. The plots in Fig. 9(a)-(b) show two regimes: weakly nonlinear at small N versus strongly nonlinear for large N . Note that N depends on both the total number of atoms and the interaction U , as one would expect. For dark solitons in particular shown in Fig. 9(a), at small N (∼ 10 −3 ) solutions are weakly nonlinear and correspond closely to the single-particle bound states of a massless Dirac spinor trapped inside a harmonic potential. Here the quantization of spatial modes can be seen by noting that the three lowest multiple dark soliton states intersect the vertical axis atμ = 2.83, 3.80, 4.88, or in terms of the oscillator frequency ω: µ = 2.83 ω, 3.80 ω, 4.88 ω. We see that these quantized modes display approximate integer multiples n of the energy ω: µ ≈ (2.8+n) ω. For large N (∼ 1), solutions are strongly nonlinear bound dark solitons  with spectra characterized by a power law:μ ∝ N α . A similar analysis applies to the bright soliton case exhibited in Fig. 9(b).
To obtain values for α in the power law fitμ = β N α , where β and α are real constants, we use the Matlab function polyfit to obtain a linear fit of the data values (N , log 10μ ) which returns a two-component vector p. With x the vector of N values and y the vector ofμ values, p = polyfit[x, log 10 (y), 1] has vector components p(1) = log 10 (α) and p(2) = log 10 (β), from which we extract the exponent α = exp[p(1)]. We find α = 0.86, 0.88, 0.94, for the three lowest multiple dark soliton states, and α = 1.4, 2.1, 8.6, for the lowest bright soliton states. Soliton density profiles and energies are important for visibility in experiments and we list these in Table 2. The density peaks and notch depths for both soliton types were computed for the solitons in Figs. 7-8. In Fig. 10 we have plotted the densities for the dark and bright multi-solitons in the zero trap limit. Both solitons extend in a series along the horizontal direction with tight confinement in the vertical direction. Different vertical and horizontal scales are used for ease of viewing. Note also that we use different density color scales in each panel, since dark solitons dip below the asymptotic value of the density, here set to 1 in our units, while bright solitons rise above it.

Conclusion
We have presented soliton stability properties in the zigzag and armchair nonlinear Dirac equation and characterized the various excitations in the core and in the bulk. Solitons in both the zigzag and armchair configurations are stable with positive or negative real eigenvalues. At finite temperatures, when non-condensate modes are appreciably populated, the negative eigenvalues allow for dissipation into lower energy Bloch states. It is important to note that suppression of these modes at low temperatures is consistent with our interpretation in terms of a metastable background condensate; one sees the same kinds of effects in analogous experiments on non-relativistic dark solitons in BECs described by the nonlinear Schrödinger equation [54]. Finite temperature corrections may be modeled by incorporating a stochastic term or more simply by including a temperature-dependent Bose distribution function to account for finite occupation of higher energy modes. However, we reserve such questions for future investigations.
We have analyzed the quasi-particle spectrum for modes localized in the soliton core and found an anomalous mode, i.e., a massless Nambu-Goldstone mode associated with phase fluctuations of the core from U(1) symmetry breaking (condensation). Moreover, inside the core we find one zero-energy mode (zero mode) associated with translational symmetry breaking by the soliton. Far from the core the spectrum consists of exotic massive excitations with quartic dispersion in addition to massless Dirac-like excitations. Hence, at low energies and near zero momentum the integrity of the Dirac point is preserved. Moreover, casting our problem in terms of pseudospin degrees of freedom places our results in the context of other domain wall theories. We found that in our case the continuous spectrum far from the core lies in the same universality class as excitations in theories which contain Fermi points such as 3 He and the Standard Model of particle physics [66].
We have computed the discrete chemical potential spectra for dark and bright solitons bound in a weak harmonic trap. Our results show two clearly distinct asymptotic regions: one for weak nonlinearity where the chemical potential for multiple soliton states differ by a constant multiple of the oscillator energy; and the other limit for strong nonlinearity where the chemical potential obeys a power law. Our numerical solutions confined in the harmonic trap yield ratio values for the notch to bulk contrast in total particle density of 0.66 − 0.77, for the dark soliton, and 2.63 − 2.86 for the peak to bulk contrast of the bright soliton. These values were computed for single, double, and triple soliton solutions. In addition, we calculated the range of the total energy for the three lowest multi-soliton states and found these to be 53.61 − 71.95 nK and 100.49 − 131.66 nK for dark and bright solitons, respectively, for a reasonable experimental parameter set for 87 Rb [25]. Density contrasts and energies offer vital comparative experimental predictions.