The nonlinear Dirac equation in Bose-Einstein condensates: I. Relativistic solitons in armchair nanoribbon optical lattices

We present a thorough analysis of soliton solutions to the quasi-one-dimensional nonlinear Dirac equation (NLDE) for a Bose-Einstein condensate in a honeycomb lattice with armchair geometry. Our NLDE corresponds to a quasi-one-dimensional reduction of the honeycomb lattice along the zigzag direction, in direct analogy to graphene nanoribbons. Excitations in the remaining large direction of the lattice correspond to the linear subbands in the armchair nanoribbon spectrum. Analytical as well as numerical soliton Dirac spinor solutions are obtained. We analyze the solution space of the quasi-one-dimensional NLDE by finding fixed points, delineating the various regions in solution space, and through an invariance relation which we obtain as a first integral of the NLDE. We obtain spatially oscillating multi-soliton solutions as well as asymptotically flat single soliton solutions using five different methods: by direct integration; an invariance relation; parametric transformation; a series expansion; and by numerical shooting. By tuning the ratio of the chemical potential to the nonlinearity for a fixed value of the energy-momentum tensor, we can obtain both bright and dark solitons over a nonzero density background.

Abstract. We present a thorough analysis of soliton solutions to the quasi-onedimensional zigzag and armchair nonlinear Dirac equation (NLDE) for a Bose-Einstein condensate in a honeycomb lattice. The two types of NLDEs correspond to quasione-dimensional reductions in two independent directions, in direct analogy to the narrowest of graphene nanoribbons. Analytical as well as numerical soliton Dirac spinor solutions are obtained. We analyze the solution space of the quasi-onedimensional NLDE by finding fixed points, delineating the various regions in solution space, and through an invariance relation which we obtain as a first integral of the NLDE. For both the zigzag and armchair geometries we obtain spatially oscillating multi-soliton as well as asymptotically flat single soliton solutions using five different methods: by direct integration; an invariance relation; parametric transformation; a series expansion; and by numerical shooting. By tuning the ratio of the chemical potential to the nonlinearity for a fixed value of the Dirac spinor kinetic energy, we can obtain both bright and dark solitons over a nonzero background. The density contrast between the dark (bright) soliton notch (peak) and the background is about 2/3 (3/1). For both soliton types, we compute the discrete spectra for several spatially quantized states in a harmonic potential. We interpret our solitons as topologically protected domain walls in a quasi-one-dimensional system which separate distinct regions of pseudospin-1/2 with S z = ±1/2. By solving the relativistic linear stability equations we obtain the low-energy spectrum for excitations in the bulk region far from the soliton core and for bound states in the core and find that excitations occur as spin waves and as a Nambu-Goldstone mode. For a Bose-Einstein condensate of 87 Rb atoms, we find that our soliton solutions are stable on time scales relevant to experiments.

Introduction
The nonlinear Dirac equation (NLDE) appears in a variety of physical settings, typically as classical field equations for relativistic interacting fermions [1,2]. In fact, the (1+1)dimensional NLDE with scalar-scalar or vector-vector interaction is the prototypical effective model for interacting fermions, and has been the subject of much analysis over the past decades [3,4,5,6,7,8]. Dirac-like spin-orbit couplings for interacting cold atoms have also been investigated, simulating some features of quark confinement [9]. Moreover, soliton solutions of the NLDE appear in quasi-one-dimensional (quasi-1D) nonlinear optical structures [10,11], acoustic physics [12], and electron propagation in graphene [13,14,15,16,17,18]. In all of these cases the combination of Dirac kinetic term and nonlinearity leads to a plethora of solitary wave solutions whose properties depend on the particular form of the interaction term [6,19]. Our own recent work has placed the NLDE in the context of a Bose-Einstein condensate (BEC) [20]. Significantly, our particular form of the NLDE has opened up research in other fields of physics [21,22,23,24,25,26,27,28,29,30,31]. For the NLDE in a BEC, the relativistic structure arises naturally as bosons propagate in a shallow periodic honeycomb lattice potential, and yields a rich soliton landscape which we explore in detail in this Article.
The quasi-1D reduction of the NLDE may be obtained by isolating a single direction in the full two-dimensional (2D) honeycomb lattice theory [20]. In practice this is accomplished by starting with the 2D lattice. One then increases the trap potential in one of the planar directions until a desired effective width is obtained so that excitations of the condensate undergo a quasi-one-dimensional (quasi-1D) reduction to a single spatial dimension in the lattice. A schematic of the harmonic magnetic trap, interfering lasers, and BEC required to realize our set up is shown in Fig. 1.
Details of the experimental construction can be found in [32]. Thus, there are two independent forms of the quasi-1D reduction of the NLDE corresponding to the armchair and zigzag patterns found in the narrowest possible graphene nanoribbons [33], but here related to the honeycomb optical lattice potential. From here on we will refer to these two forms as the armchair NLDE and zigzag NLDE, respectively [34]. The NLDE operator is complex in the armchair direction and real along the zigzag direction. Consequently, spinor solutions associated with propagation along these two directions are related by a complex Pauli matrix rotation. Linear stability of the NLDE is studied using the relativistic linear stability equations (RLSE) [20], with the NLDE+RLSE system providing a relativistic generalization of the NLSE+BDGE (nonlinear Schrödinger equation + Boguliubov-de Gennes equations) system [35].
To obtain soliton solutions we first integrate the zigzag NLDE to obtain an invariance relation which describes solutions at fixed values of the Dirac kinetic energy, a quantity which may be positive or negative valued, in relation to zero energy set at the Dirac point. The invariance relation provides a vantage point which offers insight into general solutions of the quasi-1D reduction of the NLDE. In particular, for a fixed value of the conserved kinetic energy, we find soliton solutions residing at the boundary between two oscillating solution regimes. We will refer to this boundary in parameter space as the soliton boundary. The soliton boundary appears for a particular value of the ratio µ/U = (µ/U ) SB , where µ is the chemical potential of the system and U is the quasi-one-dimensional renormalized interaction. Tuning µ/U towards (µ/U ) SB while maintaining the local particle density above some critical value, we encounter there a bright soliton, whereas a dark soliton is obtained for densities less than the critical value. Oscillating solutions away from the soliton boundary correspond to multiple dark or bright soliton and are not necessarily stable. However, the single solitons at the soliton boundary are robust objects.
To better understand the two types of solitons, we consider the two-dimensional solution space that results from fixing the internal and overall phase of a Dirac twospinor. The two types of solitons correspond to paths in solution space that interpolate between two fixed points and pass along either the small amplitude (dark soliton) or the large amplitude (bright soliton) side of a third fixed point. The two paths (and associated solitons) are topologically distinct. Our analysis centers on solutions of the zigzag NLDE but extension of our results to the armchair NLDE via the aforementioned Pauli rotation is straightforward. We point out that a substantial part of the work presented in this article is devoted to finding single and multi-soliton solutions of the NLDE. We have chosen to do this using several methods to emphasize the multiple lines of evidence for dark and bright solitons.
We determine the stability of our solitons as in our previous work [36,32], but with the appropriate dimensional reduction from quasi-2D to quasi-1D. The solution of the RLSE gives us the linear spectrum from the presence of small quantum fluctuations in the BEC. For both dark and bright soliton solutions, far from the center of each soliton the BEC occupies only one of the two sublattices, switching from the A sublattice to the B sublattice when translating through the soliton core. Thus these solution types are, far from the core, 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 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.
The work that we present here is related to a number of parallel studies in condensed matter and cold atomic gases. The dimensional reduction of the quasi-2D honeycomb lattice to a quasi-1D lattice provides a novel way to study BECs. Quasi-particles near a Dirac point acquire features similar to Dirac fermions in graphene connecting naturally to graphene nanoribbons [37,34]. For ultracold quantum degenerate bosons near the Dirac point, quasi-particles are fermionic, a feature which when combined with the typical thermal isolation of a BEC reduces the tendency for particles to dissipate into lower energy states. This partly simulates the effect of a Fermi surface. Another approach which has been proposed for simulating Dirac fermions using cold bosonic atoms relies on laser-induced spin-orbit coupling in a spinor BEC [9]. The hyperfine structure provides the internal degrees of freedom needed to simulate spin while the additional lasers couple spinor states to the spatial degrees of freedom. We note that in our case both of these effects come from the lattice background and are therefore geometric in origin. It is the combination of nonlinearity and Dirac spin structure which allows for self-localization similar to chiral confinement in relativistic models such as the massive Thirring and Gross-Neveu models [9,38,39,40,41]. Our main interest is not in simulating Dirac fermions per se, but exploring a relativistic nonlinear system in the highly tunable and controllable context of BECs, where effective relativistic velocities are 10 orders of magnitude slower than the speed of light [32].
This article is organized as follows. In Sec. 2, we provide an introduction to the NLDE and discuss the key physical parameters. In Sec. 3, we explain how the NLDE zigzag and armchair geometries are realized starting from the 2D honeycomb lattice. This step is essential in order to establish an experimental foundation for the rest of this paper. In Sec. 4, we determine general properties of the NLDE solution space. It is important to note that we treat only stationary solutions in this article; the question of dynamics may be treated in future work. Focusing on the time-independent zigzag NLDE, we find all the fixed points and regions of solution space according to the character of the associated direction fields, i.e., the vectors formed from the spatial first derivatives of the two-spinor components. We also derive the main invariance relation governing the NLDE which leads to explicit soliton solutions. In Sec. 5, we use the insight obtained by our study of fixed points and invariance relations to map out the phase diagram for NLDE solutions. In Sec. 6, we solve the NLDE analytically using a trigonometric ansatz, through detailed analysis of our invariance relation, using a parametric transformation and by a power series expansion. In Sec. 7, we obtain solitons using a numerical shooting method. In Sec. 8, we solve the RLSE numerically to determine soliton lifetimes. In Sec. 9, 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. 10, 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. 11 where we formulate relativistic solitons in the language of spin-1/2 domain walls. In Sec. 12, we analyze quantum fluctuations in light of the domain wall interpretation. In Sec. 13, we treat the spectral theory of a BEC in a weak harmonic trap. Finally, in Sec. 14 we conclude.

The nonlinear Dirac equation
In this section we introduce the NLDE, a nonlinear extension of the massless Dirac equation, and discuss the experimentally relevant physical parameters. The parameters which enter directly into the NLDE and will therefore appear in all of our solutions, are the effective speed of light c l = t h a √ 3/2 and the atom-atom binary interaction strength, U 2D = L z gn 2 3 √ 3a 2 /8, where U 2D is the 2D optical lattice renormalized version of the usual interaction g = 4π 2 a s /M . Appearing in these definitions are the average particle densityn = N/V , the s-wave scattering length a s , the vertical oscillator length L z , the mass M of the constituent atoms in the BEC, the lattice constant a, and the hopping energy t h . For the hopping energy, we use a semiclassical estimate given by [42], where V 0 and E R are the lattice potential depth and recoil energy. Typical practical values of key physical parameters in order to realize NLDE solitons in the quasi-one-dimensional regime (which we will discuss in detail in Sec. 3) are T = 8.0 nK,n = 1.5 × 10 9 cm −3 , L y ≈ L z = 3.0 µm, t h = 4.31 nK = 1.04 kHz, V 0 /E R = 16, U 2D = 0.391 pK = 0.051 Hz, a = 0.55 µm, where T is the temperature and L y is the transverse oscillator length parallel to the plane of the honeycomb lattice. For this typical parameter set we took the atomic mass M to be that of 87 Rb, and the associated scattering length a s = 5.77 nm. A complete discussion of NLDE parameters and constraints can be found in [32].
The NLDE for two inequivalent Dirac points describes the dynamics of a Dirac four-spinor of the form Ψ ≡ (Ψ + , Ψ − ) T , with the upper (+) and lower (−) two-spinors relating to opposite K and K points of the honeycomb lattice (see Refs. [20,36]). We remind the reader that Dirac points are locations in the single particle spectrum where the upper and lower energy bands become degenerate (the energy bands cross) with a linear structure, i.e., E(k) ≈ vk, a consequence of the underlying symmetry of the honeycomb lattice. For graphene the proportionality constant v is just the Fermi velocity, and the quasi-particle group velocity in the case of a BEC. In BECs v is the quasi-particle group velocity, and required to be less than the speed of sound. Note that in both cases the Dirac point is a kinetic single-particle effect, where v is determined by the microscopic physics and plays the role of an effective speed of light. In terms of the A and B sublattice wavefunctions, we have The matrices γ µ are the usual Dirac matrices and the interaction terms are encapsulated in the summation with the matrices M i constructed to give the correct cubic nonlinearities, local to each spinor component [20]. Explicitly, the interaction matrices are The interactions do not couple different spinor components, which are only coupled through the kinetic term. ‡ We focus on the equations for a two-spinor in rectangular coordinates while omitting the Dirac point subscript with the full solution expressed as a linear combination of solutions from each Dirac point. Note the presence of the effective speed of light, c l , and interaction strength, U 2D . Equations (4)-(5) allow for quasi-one-dimensional (quasi-1D) solutions by confining the BEC tightly in one of the planar directions. From here on we will assume tight confinement in the y-direction. The experimental construction is shown in Fig. 1 where the BEC resides within a weak magnetic harmonic trap and an optical potential created by three laser beams offset by relative angles of 120 o in the x, y-plane. With tight confinement in the y-direction, ψ A and ψ B are functions only of x. Then the stationary states of interest are obtained by taking the time-dependence to be the usual exponential ‡ Note that were we to consider nearest-neighbor interaction terms in the quantum Hamiltonian, we would have such a coupling; however, here we make the usual tight-binding, lowest band approximation, as such intersite interaction terms are small -see Ref. [32].
factor with the chemical potential µ as the frequency: The system, Eqs. (6)- (7), is the time-independent quasi-1D NLDE for the armchair direction of the honeycomb lattice. Notice that taking f A → if B and f B → f A converts the armchair NLDE to the zigzag version as can also be obtained by choosing the ydirection in Eqs (4)- (5). This transformation is equivalent to multiplication by a linear combination of Pauli matrices with an overall complex factor [(1 + i)/2](σ x − σ y ). It is natural to retain the x notation when discussing either form of the NLDE; thus we write the zigzag NLDE as In Sec. 3, we will see that the interaction U 1D in Eqs. (6)- (9) is the renormalized version of corresponding 2D interaction U 2D in Eqs. (4)-(5).

Quasi-one-dimensional reduction in the honeycomb lattice
The zigzag and armchair NLDE are physically realized by starting with the 2D honeycomb lattice, then adding tight confinement in the y-direction. The potentials for the honeycomb lattice with harmonic confinement for the armchair and zigzag geometries are given by where α is the polarizability of the atoms, E 0 is the electric field strength, r = (x, y) is the planar coordinate vector, M is the atomic mass of the constituent bosons, and ω y is the frequency of the harmonic potential which adds the additional confinement. The wave vectors k 1 , k 2 , and k 3 in Eqs. (10)-(11) are defined as The different forms of V armchair and V zigzag in Eqs. (10)-(11) reflect a uniform rotation of the beams by 90 o . The wavelength of the laser light that forms the lattice is λ = 2π/k which defines the distance between sites on a hexagonal sublattice as a = 2λ/3. The underlying electric field which produces the potentials in Eqs. (10)-(11) is the superposition of the fields from each beam and is given explicitly by where ω i = c|k i | (c = speed of light ≈ 2.99×10 8 m·s −1 ) here must be distinguished from the trapping frequency, andẑ indicates the polarization of the beams perpendicular to the plane of the lattice. The lattice potential is then obtained by taking the average of the square of the electric field where V here can stand for either potential in Eqs. (10)- (11). Plots of the armchair and zigzag potentials are shown in Fig. 1(c)-(d).
Transforming to the quasi-one-dimensional regime requires that a s L y ξ [43,44], where L y characterizes the width of the condensate along the y-direction. This condition ensures that the condensate remains in the ground state in the transverse direction. In addition, we require that L y L x , which ensure that excitations along the x direction have much lower energy than those in the y direction. Equivalently, these constraints may be expressed in terms of the trap frequencies and the atom-atom interaction: ω x ω y , ω x U . For the moment we neglect the harmonic trap in x in Eqs. (10)-(11), but will come back to its effects in Sec. 13. Finally, we require that the 3D s-wave scattering length a s is far enough removed from potential resonances for the purposes of this Article. A modified renormalized interaction U and careful consideration of phase coherence in quasi-1D BECs should be sufficient to account for extremely tight confinement, but are not necessary here [45,43,46,47].
Dimensionally reducing the honeycomb lattice translates to the spinor component wavefunctions. To see this explicitly, we separate ψ A (r, t) and ψ B (r, t) into longitudinal and transverse modes following similar arguments as in Ref. [48]: where f A(B) (x) contains the x-dependence and h(y) is the dimensionless function that describes the transverse part of the wavefunction. Projecting onto the ground state h gs (y) in the transverse dimension gives us an effectively quasi-1D wave equation. With h gs (y) = π −1/4 exp (−M ω y y 2 /2 ) and L y = /M ω y , substituting Eqs. (17)-(18) into the NLDE obtained for the lattice potential in Eq. (10), multiplying through by h gs / L y , and integrating over y gives Similarly, for the potential in Eq. (11) where the prime derivative notation is clear since we only have x-dependent functions. Note that we have introduced the renormalized interaction U 1D defined as A key result here is that the chemical potential is not renormalized since the Dirac equation is first order in the spatial derivatives and the extra term proportional to dh gs /dy is an antisymmetric function of y which vanishes upon integration. To simplify the notation, for the rest of this article we will use the plain notation U and understand that this refers to the quasi-one-dimensional renormalized interaction.

General properties of NLDE solutions: fixed points and invariance relations
As a first step towards solving the NLDE, we map out the solution landscape by understanding the character of the various solution types. In this section and throughout our work we confine our analysis to stationary solutions, leaving the case of dynamics for future investigations. In particular, we require a clear understanding of the points where solutions are constant (zero spatial derivative), and the flow of solutions near these fixed points. Working from the zigzag NLDE we look for real solutions and write Eqs. (8)-(9) as a derivative field (or direction field) where the dependence on x is implied. Together, Eqs.  Table 1 along with the corresponding signs for the derivatives f A and f B . Nine fixed points exist in the solution space (f A , f B ): (0, 0), (± µ/U , 0), (± µ/U , ± µ/U ), (0, ± µ/U ). A qualitative analysis of each region in Table 1 provides a guide to the types of soliton. Two main types of solutions exist: oscillating solutions which do not flatten out asymptotically, and localized solutions whose derivatives vanish asymptotically. The former turn out to be variants on multi-solitons, while the latter turn out to be varieties of single solitons. In addition, solutions differ qualitatively depending on which regions in Table 1 are involved. At one extreme, strongly oscillating solutions exist for which f A and f B both have amplitudes greater than µ/U , in which case each component experiences three turning points during a half cycle. We give an example of this type of solution in this section, which we obtain analytically. Other oscillating solutions occur around each of the turning points (± µ/U , 0) which may be discerned from Eqs. (25)-(26) using linear perturbation theory by substituting f A (x) = coskx ± µ/U or f A (x) = coskx, with the small amplitude such that / µ/U 1. These solutions cross two region boundaries with each derivative f A and f B changing sign once per half cycle. We also find non-oscillating asymptotically flat solutions having a similar form to the function (1/2)[1 ± tanh(x)], which we study in Sec. 5. This spinor solution has a constant total density everywhere except near a localized region where a dip, or notch, in the density occurs, a form which describes a dark soliton. In addition, we find that a bright soliton occurs in the combined region I+II. Figure 2 gives a qualitative schematic depiction of the various regions, oscillating solutions centered on fixed points, and asymptotically flat solutions which interpolate between fixed points.
To deepen our analysis, a detailed map of the solution space of Eqs. (25)- (26) can be arrived at by uncovering spatially invariant quantities. To obtain the first invariant quantity we multiply Eq. (25) by the right hand side of Eq. (26) and vice versa, then add the resulting equations to obtain which then gives where the prime notation indicates differentiation with respect to x. Equation (28) expresses the equal contribution of the positive and negative kinetic energy states to the total spectrum. Integrating Eq. (28) gives a relation between the functions f A and interpolating between fixed points. The regions described in Table 1 are highlighted in yellow.
where we have simplified the expression by completing the squares in f A and f B , and C is the integration constant. Equation (28) has a clear meaning when we consider that for a Dirac system the total energy density is given in terms of the functions f A and f B by so that the kinetic energy density is Thus we interpret C as the local kinetic energy density which may be positive or negative for any given solution, as is evident from the antisymmetric form of Eq. (31). This can be understood in light of the Dirac point structure: positive (negative) energy means excitations above (below) the Dirac point. Solutions for which Eq. (29) holds will have spatially invariant total kinetic energy E kin = C. A second invariant (or fixed) quantity is arrived at by expressing Eqs. (25)- (26) in the form then adding, integrating, and combining terms to get a total derivative on the left hand side: Equation (35) expresses the gradient of the local total density ρ( In particular, Eq. (36) allows for solutions which asymptotically approach a constant value, i.e., with fixed quantity Thus, such stationary solutions are topologically stable in that local fluctuations in f A and f B will cancel out such that Q remains fixed. Note that Q depends only on the difference between the asymptotic values of the total density ρ.

Transition of solutions across the soliton boundary
We now combine the qualitative information and the invariance relation (Eq. (29)) from Sec. 4 to arrive at a more technical understanding of the relationship between oscillating and asymptotically flat solutions of the NLDE. In particular, by choosing specific values for the chemical potential µ and strength of nonlinearity U , Eqs. (25)-(26) may be solved by numerical or analytical methods. In general though, such solutions may be highly oscillatory and not necessarily stable. In this section we study the evolution of oscillating solutions as µ and U are tuned to obtain particular soliton solutions which are stable. We may better understand solutions of Eqs. (25)-(26) by plotting Eq. (29) for a particular value of C. This allows us to see how solutions evolve as we vary the ratio µ/U . For example, for C = −1 we have plotted the solution space for µ/U = 0.8, 0.9, 0.98, 1.0, 1.02, 1.25, 1.5, 1.75, 2.0, in Figs. 3(a)-(i). In Fig. 3(a)-(c), solutions oscillate about the fixed points (f A , f B ) = (± µ/U , ± µ/U ) with the orbits beginning to coalesce in Fig. 3(c). In Fig. 3(d), µ/U = 1 and solution paths begin and end at the saddle fixed points (0, ± µ/U ), (± µ/U , 0), asymptotically flattening out for large positive and negative x. In this case there are two distinct types of solitons indicated by the paths 1 and 2 that connect the points (1, 0) and (0, 1). We shall see that these solutions correspond to dark and bright solitons as previously mentioned. As µ/U is increased from unity, orbits bifurcate into solutions which oscillate about (0, 0), with small and Figure 3. Evolution of NLDE solution orbits in parameter space. For a chosen value of the fixed quantity, C = −1, orbits evolve as the ratio of chemical potential to interaction, µ/U , is tuned from values less than one to greater than one. For the particular value µ/U = 1, shown in (d), solutions labeled as paths 1 and 2 begin and end on the fixed points (0, 1) and (1, 0) (refer to Table 1). Solutions must approach (or leave) a fixed point asymptotically for x → +∞ (−∞), so that the non-trivial (nonzero spatial derivative) part of the solution is spatially localized. In general, soliton solutions such as those depicted in Fig. 3(d) correspond to the case C = −(µ/U ) 2 requiring C < 0 in Eq. (29). In addition, the dark and bright solitons are distinguished by the density conditions which result from Eq. (29), where we've added the subscripts DS and BS to the total density to indicate dark and bright solitons, respectively. Note the inclusion of the average particle densityn. The densities are defined in terms of the real spinor spatial , with an analogous definition for the bright soliton. Note that the particular cases where ρ =nµ/U or ρ = 2nµ/U correspond to fixed points of the NLDE, where only the trivial constant solutions exists. In the case where ρ > 2nµ/U , no isolated single-soliton solutions exist.
Thus, solitons exist on a 3-dimensional sub-manifold of parameters (defined by the condition µ = |C| 1/2 U ) of the 4-dimensional parameter manifold determined by the low-energy oscillations  chemical potential, interaction, density, and kinetic energy with coordinates denoted as (µ, U, ρ, C). Moreover, the density conditions in Eq. (39) further partition the 3D parameter subspace into the two types of solitons along the boundary ρ =n|C| 1/2 = nµ/U . The concept of a soliton boundary is useful in order to visualize the transition from oscillating solutions at weak nonlinearity into oscillating solutions at strong nonlinearity, with single isolated solitons appearing at the boundary between the two oscillating regimes. In particular, for C = −1 the soliton boundary occurs at (µ/U ) SB = 1. Tuning µ/U → 1 while keeping ρ <n forces the solution to the dark soliton, whereas tuning µ/U → 1 while keepingn < ρ < 2n converges on the bright soliton. Note that the type of soliton obtained is independent of whether µ/U approaches 1 from above or below. Conversely, if we maintain the condition µ/U = 1 while tuning ρ through the critical value ρ c =n we induce a transition between the dark soliton and bright soliton. Thus, tuning µ/U moves the system between oscillating regimes across the soliton boundary, while tuning ρ/n along (µ/U ) SB moves the system between the dark and bright solitons. Quantum phase transitions across the soliton boundary are summarized in Fig. 4. This kind of phase transition in the mean-field theory indicates a possible corresponding quantum phase transition in the underlying microscopic theory [50,49,51,52]. A more exact phase diagram could be calculated from the many body theory via the RLSE, forming a subject for future work. For contemporary works on quantum phase transitions see Ref. [53].
To make our analysis more concrete, we solve the numerical initial value problem defined by Eqs. The values of f A (0) were chosen to pick out branch 2 with f B (0) determined by inverting and solving Eq. 6. Analytical solution methods

Oscillating multi-soliton solutions
Large amplitude oscillating solutions for strong nonlinearity such as those in Figs. 5(g)-(i) may be obtained analytically. Such solutions are bright solitons in the total density over a nonzero background. We start by writing the two-spinor wavefunction in the form of a product of an envelope function η(x) and the internal spinor degrees of freedom parameterized by the function ϕ(x) where we have assumed only that the wavefunction is real, i.e., choosing to work from the zigzag NLDE. Note that there is an arbitrary overall phase constant and translation symmetry x → x 0 , which may included in the final solution. Substituting Eq. (41) into Eqs. (8)- (9), multiplying by cosϕ and sinϕ, respectively, then adding the resulting   equations gives To obtain a second equation we multiply Eqs. (8)-(9) by cosϕ and sinϕ, respectively, then subtract the resulting equations which yields Note that we have divided through by η to arrive at Eqs. (42)- (43). Equations (42)- (43) can be combined by back substitution to get sin(4ϕ) where the prime notation indicates differentiation with respect to x. A formal expression for η is obtained from Eq. (44): To solve Eq. (45), we assume the linear form ϕ(x) = κx, obtain an explicit form for η(x), which we then substitute into Eq. (42) to determine the constant κ and obtain a relation for the chemical potential µ and the interaction U . Equation (45) becomes which, upon integration, yields the result where A is the integration constant. Substituting this result into Eq. (42) and using the linear assumption, gives the expression In order for this expression to hold true it must be that the exponent of the spatial functions is identically zero. Equation (48) then gives the two conditions which may be solved to give The corresponding solution is then where κ = 2U/( c l ). The spinor components in Eq. (53) are plotted in Fig. 7(a) and the corresponding density in Fig. 7(b). Although Eq. (53) was obtained for the zigzag NLDE, a direct transformation to get the associated armchair solution is obtained by  The spatial dependence of the total particle density. The spatial pattern of the density reveals the oscillating solution as a series of equally spaced bright solitons over a nonzero background.
A distinguishing feature of the solitons in Fig. 7 is that the spinor component functions oscillate about zero but the total density does not exhibit nodes. From the density standpoint, the peaks in Fig. 7(b) are bright solitons on a nonzero background.

Dark and bright solitons by parametric transformation
In this section we isolate single dark and bright soliton solutions analytically using a parametric transformation motivated by the form of the invariance relation Eq. (29). A preliminary step requires that we solve Eq. (29) for f B and then back substitute into Eq. (26). This gives us the first order nonlinear equation Equation (55) is separable in the variables f A and x. In the special case where µ = 0, the problem simplifies and yields the integral where D is an integration constant. For |f A | < C 1/4 and C > 0, Eq. (56) shows that f A will remain positive for all x; thus a monotonically increasing and bounded solution exists between −C 1/4 < f A < +C 1/4 , for −∞ < x < +∞. The integral in Eq. (56) is given in terms of the hypergeometric function so that The second spinor component f B is obtained by back substitution into the invariance relation where we find f B = (C − f 4 A ) 1/4 . This type of solution exhibits a form of selfconfinement similar to that studied in [9], where f B remains localized near the region where the slope of f A is steepest and has a shape similar to a tanh function.
. Graphical representation of dark and bright solitons using the parametric transformation. Cast in terms of the angular parameter θ, the dark soliton emerges as an interpolation between θ = π/2 and θ = π, in contrast to the bright soliton which interpolates between θ = π/2 and θ = −π.
For general values of µ, it is helpful to cast Eq. (55) in a more enlightening form by an appropriate parametric transformation. First, we make the substitution Note that in order to keep f A real we must enforce the constraint a > b. Cast in terms of the new variable θ, we have the symmetric result b cosθ = f 2 A − a and b sinθ = a − f 2 B . Equation (58) has an interesting graphical representation which we have depicted in Fig. 8 where we show how the dark and bright solitons emerge when the problem is cast in terms of the angular parameter θ.
The dark and bright solitons may be found in the limit that a → 1, b → 1, and with the choice C = −1 as in Sec. 5. In this case dθ/dx = 0 at θ = π/2 and θ = π for the negative sign under the radical in Eq. (58), and dθ/dx = 0 at θ = −π/2 and θ = π for the positive sign under the radical. In particular, for the negative sign under the radical, the dark soliton and bright soliton interpolate between θ = π/2 and θ = π by following anti-clockwise and clockwise paths, respectively, i.e., for the positive and negative outer signs in Eq. (58).
Equation (58) may be integrated exactly by separation of variables whereby one obtains where D is the integration constant which adds a spatial translation to the soliton core. Equation (59) cannot be directly inverted to find explicit forms for f A (x) and f B (x) (via θ(x)), but we may study it graphically by plotting U x/ c l as a function θ. In Fig. 9 we have plotted Eq. (59) where we have shown singular points occurring at repeating intervals: ..., π/2, π, 5π/2, 3π, ..., alternating between the two types of soliton solutions. Alternatively, Eq. (58) may be solved numerically to obtain an inversion of the solution plotted in Fig. 9, i.e., θ(x). Figure 10 shows this numerical result for the dark soliton in (a) and bright soliton in (b) with the densities for each soliton type shown in the lower panels. Note that the spinor components are nonzero within the soliton cores. Also, there is a clear signature associated with the densities of each soliton type: the dark soliton density dips to a factor of 0.6 of the asymptotic background density whereas the bright soliton peaks at 3.5 times the background.
The invariant form of Eq. (36) says that asymptotically flat soliton solutions must have the same topological invariant quantity Q given by Eq. (37). We verify that this is the case by explicitly computing Q using Eq. (58) and expressing the integral in Eq. (37) in terms of the variable θ. The resulting indefinite integral for Q is When evaluated at the appropriate limits for either of the two types of solitons, the integral in Eq. (60) yields Q = 0. For the dark soliton θ i = π/2 and θ f = π, whereas θ i = π/2 and θ f = −π for the bright soliton. We point out that  This analysis demonstrates that our solitons are indeed of two distinct types, unrelated by any continuous transformation in parameter space.
The character of the parametrization angle θ in Fig. 8 illustrates this point. For example, we are free to parametrize our problem more generally in terms of some variable τ : b → b(τ ) and θ → θ(τ ) in Eq. (58). And, since the singular point (b = 0, θ) is a fixed point of Eqs. (25)- (26), any continuous deformation of the functions b(τ ) and θ(τ ) is allowed provided we avoid the origin where b = 0. Thus, just the fact that the dark and bright soliton paths encircle the origin in opposite directions suggests distinct homotopy classes for the two types of solitons.
It is important to emphasize that we classify our solitons as dark or bright based on the density profiles shown in Fig. (10)(c) and (d). For our single isolated solitons the densities occur as either a suppression (notch) or elevation (peak) with respect to a nonzero condensate background: no nodes are observed in the density profiles. Nevertheless, the underlying spinor components asymptotically approach zero in one direction. This is in marked contrast to the usual case of a single component BEC. For instance, bright solitons occur in presence of an attractive interaction, as in general NLSE based theories with focusing nonlinearity (attractive), while a defocusing nonlinearity (repulsive) leads to a dark soliton profile. It is worth noting that dark and bright solitons may occur in the same system such as in extended Bose-Hubbard models with strong repulsive on-site interactions [54,55,56]. A key feature of this particular model is that the binary nature of the on-site occupation number allows for a mapping to a spin-1/2 system. The particle-hole asymmetry in the spin-1/2 model characterizes the crossover from a nonlinear Schrödinger order parameter to that of a strongly repulsive system and gives rise to dark and bright multi-species setting. In spite of these similarities, our solutions are multi-component bright and dark solitons directly relating to the relativistic context.

Soliton series expansions
A fourth approach to obtain solutions of the NLDE is through a series expansion. This approach allows for general values of the ratio µ/U , interpolating through the soliton boundary and connecting the two regimes shown in Fig. 3. The series expansion uses the same ansatz as that in Sec. 6.1, i.e., Eq. (41), whereby we find that η(x) will take the form of a power series in the quantity cos 4 ϕ + sin 4 ϕ. The function ϕ(x) can then be obtained by integrating η(x) term by term. This offers a convenient method since powers of cos 4 ϕ + sin 4 ϕ are exactly integrable.
We begin by substituting the form Ψ(x) = η(x) (cosϕ(x), sinϕ(x)) into the invariance relation Eq. (29). Upon this substitution, Eq. (29) becomes where the expansion coefficients are the generalized binomial coefficients The terms in the expansion for ϕ(x) can be integrated and expressed in terms of elliptic integrals. The leading order terms in the asymptotic expansion are where the solution for ϕ(x) is expressed in terms of the inverse elliptic integral of the first kind [57]. Conversely, at the other extreme limit |C(U/µ)| 2 1 the following expansion is valid To leading order, only one nonzero solution exists coming from the positive sign inside the square brackets in Eq. (61) η(x) ≈ (U/2µ) 1/2 We point out that Eq. (71) is just the series of bright solitons first obtained in Eq. (53) which we now identify as the large µ/U limit solution in Fig. 5(i). A second solution can be obtained associated with the negative sign inside the brackets in Eq. (61) and coming from the next to leading order term in Eq. (68) which describes plane wave solutions identified with the oscillating solutions in panel A solution for ϕ(x) can be obtained by the method of separation of variables, whereby we find sin(2ϕ) ln [tan(2ϕ)] where D is the integration constant which spatially shifts the solution. One may invert Eq. (76) graphically and substitute this result into Eq. (75) which leads to the dark and bright solitons, respectively, for the positive and negative signs in Eq. (75).

Solitons by the method of numerical shooting
So far we have focused our attention on solving the NLDE using analytical methods. In Sec. 7, we provide a detailed analysis of solutions by the method of numerical shooting. Numerics provide a versatile angle of attack and prove especially convenient for solitons confined by an external trap. Our numerical approach in this section is similar to that used for studying trapped BECs in the absence of a lattice background [58]. The most direct approach is to express Eqs. (8)- (9) in terms of the dimensionless spatial variable χ ≡ (U/ c l )x. The spinor functions f A (χ) and f B (χ) are then expanded in a power series around χ = 0 with a j and b j the expansion coefficients to be determined. Since we are solving two coupled first order equations, we require the initial conditions f A (0) and f B (0). Substituting Eq. (77) into Eqs. (8)- (9) gives us the behavior of the solution at the origin By examination we find that for j ≥ 0 the NLDE implies a recursion relation for the expansion coefficients a j and b j . Recursion relations for the lowest values of the index j are resulting from Eq. (8), and which come from Eq. (9). To obtain soliton solutions using the shooting method, we first fix either a 0 or b 0 , and then vary the other until we obtain convergence to the desired precision. Although it seems that both a 0 and b 0 are free parameters, fixing one to a different value before shooting results in a spatially translated final solution which we verify numerically. This is true as long as the first parameter is fixed to a value between zero and one, for the dark soliton, or between one and the value of the peak for the bright soliton. The second "shooting" parameter can then be tuned to find the stable soliton solution. Taking b 0 = 0 and iterating to obtain the value for a 0 = a soliton 0 to the desired precision gives the soliton configuration. Figure 11 shows the shooting process as we tune a 0 through the dark soliton which is shown in Fig. 11 At larger values of a 0 , the nonlinearity becomes dominant and we start to pick up some of the excited soliton states which can be seen in Figs. 11(a) and (b). For values a 0 < a soliton 0 , the effect of the interaction is reduced, Figs. 11(d) and (e), until finally we see the free particle sine and cosine forms appearing in Fig. 11(f). The particular values of the constant a 0 in Fig. 11 are: (a) a 0 = 1.1, (b) a 0 = 0.9992922, (c) a 0 = 0.999292150145, (d) a 0 = 0.99, (e) a 0 = 0.9, (f) a 0 = 0.5. The bright soliton can be obtained using a similar shooting process. A similar method has been used to study convergent ring dark and bright solitons [59].

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. In this section, we compute the linear spectrum for soliton solutions of the zigzag and armchair NLDE.
To compute soliton lifetimes we must solve the relativistic linear stability equations (RLSE) modified for our quasi-one-dimensional problem [20]. 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 coherence 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 quasi-particle coherence functions. This substitution gives a set of first-order coupled ODEs in one independent variable to be solved consistently for the quasiparticle energies E k and amplitudes u k and v k , where the subscript denotes the mode with momentum p = |k|. Since we are perturbing from a spin-1/2 BEC background, , v kB (x)] T 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 (91)-(94) 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. (91)-(94) are renormalized upon dimensional reduction from 2D to quasi-1D as described in Sec. 3 specializing to the zigzag or armchair geometries. We point out that Eqs. (91)-(94) pertain to the zigzag NLDE, whereas in the armchair version the derivative terms have identical coefficients of −i c l , with i ≡ √ −1. The RLSE respect the same symmetry which converts between the zigzag and armchair NLDE forms via the two-spinor Pauli transformation discussed in Sec. 2: 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 1 = ±0.1862 and E BS 1 = ±0.1902, in units of the quasi-1D renormalized interaction, where the superscripts DS and BS refer to the dark and bright solitons, respectively, for both the zigzag and armchair NLDE. Figure 12 shows the associated quasi-particle functions which are bound states at the defect point of the soliton, i.e., near the region where the density transitions from the A to the B sublattice. The bound states shown in Fig. 12 decay far from the soliton core where the continuum of scattering states is dominant. The negative 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 [60].

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. (91)-(92) and Eqs. (93)-(94). 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 , E = E/U ,x = U x/ c l , we obtain where we have defined new variables Next, we make the substitution Substituting Eqs. (106)-(107) into Eqs. (102)-(103), 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. (102)-(103) are reduced to where Following a similar path yields the reduced forms for Eqs. (104)-(105) where Equations (114)-(115) and (118)-(119) 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. and , (127) 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 (114)-(115), (118)-(119) using approximate forms for the dark soliton spinor components The approximate forms in Eqs.
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 two regimes: (i) all potentials are negative on the interval −∞ < x < +∞ when |Ẽ| < 2 − √ 3, which lead to exponentially growing and decaying solutions, and (ii) two of the potential functions, 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. 10.
We can interpret the regimes (i) and (ii) in light of the density and phase fluctuations, Eqs. (126)-(127). In regime (i), the fluctuations reduce to and A zero-mode solution (E = 0) is supported for which 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, dx [|u A(B) | 2 − |v A(B) | 2 ] = 0, and must instead be normalized according to dx |u A(B) | 2 = dx |v A(B) | 2 . 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 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.  which yields the spectrum where here we work in the renormalized parametersẼ ≡ E/U ,μ ≡ µ/U ,p ≡ p/U . The long wavelength limit of Eq. (137) 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. (137) is plotted in Fig. 13. 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
It is instructive to view the NLDE from a more elegant perspective by recasting the problem 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 [61]. 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π [62]. In spin-1 BECs, domain walls have been studied extensively as boundaries between regions of pseudospin polarization S z = ±1 [63,64], in addition to investigations into the quasi-particle transmission and reflection properties of such boundaries [65].
Domain walls also play an important role in high energy physics, for example as extended supersymmetric objects which isolate different vacua [66,67]. It is thus not surprising that solitonic objects appear at the interface between condensed matter and particle physics. A particular example which highlights this fact is the recent simulation of tachyon condensation using two-component BECs [68]. 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 ground states of the system. The key feature of the deformation is that it is localized; it occurs over a finite region in one of the spatial dimensions.
In the case of NLDE solitons the asymptotic values of the spin components correspond to two different ground states of the NLDE 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 twospinors. 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 and with the additional Pauli matrix σ 1 chosen for the specific case of the zigzag or armchair geometry A thorough description of this formalism may be found in Ref. [20].
To make the connection between superfluidity and relativistic current apparent in the pseudospin formalism, we use the Gordon decomposition approach [69] which lets us separate the phase gradient from the magnitude gradient in Eq. (140). The former describes conventional superfluidity and the latter contains the pseudospin dependent contribution to the overall current. In our previous work [20] we found that the NLDE has the associated current whereψ ≡ ψ † σ 0 . Solving Eq. (140) for ψ and the conjugate of Eq. (140) forψ allows us to express Eq. (144) as where N −1 is the inverse matrix of Eq. (141). 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. (145) can now be expressed as where we use the ∂ ± ν notation to mean the antisymmetric differentiationψ M µν ∂ ± ν ψ ≡ ∂ νψ M µν ψ −ψM µν (∂ µ ψ). The first term in Eq. (147) 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. (147) 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 Ψ(x, t) = e iµt/ [f A (x), f B (x)] T for which we find Thus, for the ordinary dark soliton and bright solitons shown in Fig. 10, 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 the 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 [70] but more closely resembling the boundary between pseudospin domains in two-component BECs. The analogy extends only loosely to spin-1 BECs [65] 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.
The measurable physical quantity is the average pseudospin or spin density defined asF(x, t) =ψ † (x, t)Sψ(x, t). In the BEC limit, a pure state with occupation of the A sublattice only, Ψ(x) = [ψ A (x), 0] T , has whereas occupation of the B sublattice has F(x) = (− /2)|Ψ B (x)| 2 (0, 0, 1). Note that F is defined in terms of single-particles states and so does not have the "hat" operator notation. The total particle density of our solitons varies through the core (the dark soliton dips, the bright soliton peaks), but returns to the constant value lim x/ξ→±∞ ρ = +µ/U with µ/U = 1 for NLDE solitons. In contrast, the spin density has asymptotic values lim x/ξ→−∞ F = ( /2)(0, 0, 1) T , lim x/ξ→∞ F = ( /2)(0, 0, −1) T . This allows for the interpretation of NLDE solitons as boundaries or domain walls separating regions of spin +1/2 and spin -1/2, with the transition between the two asymptotic regions taking place within the soliton core.

Spin waves and Nambu-Goldstone modes
In this section we will extend the discussion from Sec. 11 to the asymptotic linear spectrum far from the soliton core. 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 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. [65]). 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. 9, 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, δF SW = ( /2)(φ B + φ * B ), and δF NG = ( /2)(φ A + φ * A ). In contrast, in the other limit x/ξ → +∞, F z = − /2, δF SW = ( /2)(φ A + φ * A ), and δF NG = ( /2)(φ B + φ * B ). 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. 14.

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 Sec. 3 with the soliton in the xdirection, we encounter only the effects resulting form quantization of spatial modes along the soliton direction in either the zigzag or armchair geometries shown in Fig. 1. We then take the 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. [59]). 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. (160), 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. (163) is thus a direct result of the relativistic linear dispersion of the NLDE. Length scales are defined in Table 2 along with associated momentum and energy scales for quasi-1D NLDE solitons in a harmonic trap. Note that the last three scales in Table 2 are key in our calculations since they contain information about the first two scales.

Physical scale
Length Momentum Energy 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. (164) 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 the parameter values previously stated in Sec. 2 and a trap oscillator frequency ω = 2π × 0.039 Hz, we obtain latt ≈ 2.3 µm, ξ 1D ≈ 10 µm, and trap ≈ 55 µm. § Condition (164) can be overcome by turning to a discrete model but still working in the mean field approximation [20].
To connect with the solutions we have already found, namely the dark and bright solitons in Fig. 10(a)-(b), 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 2. Expressed in terms of the lengths in Table 2 one finds that the derivative, interaction, and chemical potential terms in Eqs. (161)-(162) 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.
Our numerical shooting method then follows the discussion in Sec. 7 and confirms the analytical results displayed in Fig. 10(a)-(b). Here we see the advantage of using numerics over analytical methods in that we are able to isolate any multiple dark or bright soliton solutions over a nonzero background. Panels (a),(c) and (e) in Fig. 15 show spinor components for single, double, and triple dark solitons with corresponding densities in panels (b), (d) and (f). Analogous plots for the bright soliton are shown in Fig. 16. Nodes only appear in the spinor component functions for the case of multiple solitons but not for single ones, but in every case the total density never drops to zero.
The following data was used to obtain the dark solitons in Fig. 15: a 0 = 0.9949684287783,μ = 1, for the single soliton; a 0 = 0.99496892372588591202,μ = 1.00000103, for the double soliton; and, a 0 = 0.993,μ = 1.001, for the triple soliton. The bright solitons in Fig. 16 are associated with the following data: a 0 = 0.010,μ = 1, for the single soliton; a 0 = 0.11,μ = 1.04, for the double soliton; and, a 0 = 0.1, µ = 1.15, for the triple soliton. The solutions are converged to the last digit in the numerical values for a 0 . Greater precision in the value of a 0 is required in the case of the single soliton in order to push excitations out to larger values of x. See also [58] for a study of precision issues in shooting methods related to BEC in a trap.
Next, we solve Eqs. (161)-(162) 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 = 0.001, we find the free parameter a 0 for the single dark soliton at a 0 = 0.94640402384 and for the two and three soliton states a 0 = 0.89882708125 and a 0 = 0.8523151, respectively. The lowest multiple dark solitons in a trap are plotted in Fig. 17 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 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.
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. (161)-(162), 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. 19. The plots in Fig. 19(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. 19(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. 19(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 3. The density peaks and notch depths for both soliton types were computed for the solitons in Figs. 17-18. In Fig. 20 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.

Conclusion
In this article we have presented a variety of methods for solving the zigzag and armchair nonlinear Dirac equation NLDE and mapped out the soliton landscape. The two types of NLDEs are associated with the two orthogonal directions of a quasi-one-dimensional reduction in the plane of a honeycomb optical lattice. The order parameter for bosons propagating along the zigzag pattern in the lattice relates to that along the armchair direction through a transformation expressed in terms of Pauli matrices. In particular, we have found a dark and bright solitons for nonzero chemical potential and a self-confined soliton for the case µ = 0. Both solitons show highcontrast density fringes and are therefore clearly observable in experiments. The spinor component functions for the dark soliton are approximated by the forms ψ A (x) ∼ (1/2)[1 + tanh(x)] and ψ B (x) ∼ (1/2)[1 − tanh(x)], while the bright soliton components resemble the forms ψ A (x) ∼ (1/2)[1 + tanh(x) + sech(x)] and ψ B (x) ∼ (1/2)[1 − tanh(x) + sech(x)], where ψ A (x) and ψ B (x) are the spatial parts of the upper and lower two-spinor components. The crossing point where ψ A = ψ B lies below the fixed point ψ A = ψ B = µ/U , for the dark soliton, and above the fixed point in the case of the bright soliton. We have found that a continuous deformation between the dark soliton and bright soliton forces the solution to flatten out everywhere when pushing through the fixed point which violates conservation of total particle number. Thus, the two solution regimes, large and small local densities, are topologically distinct. In the case where µ = 0, a third soliton solution exists which resembles a flattened sech function in ψ B localized at the center of a tanh form in ψ A .
Solitons in both the zigzag and armchair configurations are stable with positive or negative real eigenvalues, the negative values allowing for dissipation into lower energy Bloch states. We have analyzed the quasi-particle spectrum for modes localized in the soliton core and found an anomalous mode (Nambu-Goldstone mode) associated with phase fluctuations of the core. The quasi-particle spectrum in regions far from the core resemble excitations in other systems where a domain wall is present. 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 [71].
We have computed the discrete spectra for dark and bright solitons bound in a weak harmonic trap which clearly show two 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 . Density contrasts and energies offer vital comparative experimental predictions.
A distinguishing feature of our results is that our dark soliton does not exhibit a node in the total density, even though there are asymptotic zeros in the individual spinor components. An analogous effect occurs in dark solitons in Fermi gases, where a node appears in the gap function ∆(x), but not the total density. In the case of the Fermi gas, a dark soliton may be observed in the gap function ∆(x), the order parameter that encodes pairing of fermion particles and holes in the Bardeen-Cooper-Schrieffer (BCS) to Bose-Einstein condensate crossover. There the notch depth in ∆(x) varies with the translational speed of the soliton. The speed is characterized by an overall phase φ(x) which itself varies across the gap notch. Moreover, the density of Bogoliubov modes inherits a similar gray soliton profile and has been investigated in recent theoretical and experimental work [72,73]. The soliton notch in both the gap and the density becomes deepest and narrowest at unitarity: the regime between BCS and BEC where the pairing length is on the order of the atomic spacing.
Similarities to the Fermi gas brings up several questions worth addressing. First, it is interesting to consider how the combined tuning of a complex gap and atom-atom interactions in the NLDE affect the notch depth and overall phase through the core of our solitons. Second, does a finite velocity boost induce Friedel-like oscillations in NLDE solitons similar to those studied in [72]? Third, what effect does added motion of the soliton and the presence of a gap have on solutions of the RLSE? Moreover, is there a direct relationship between the RLSE and the Andreev equation, a relativistic form of Bogoliubov's equations that emerges in the limit of weak attractive interfermion forces [74]? All such questions relate to the issue of Cooper pairing between NLDE quasiparticles, the relationship between composite particle-hole pairs and the fundamental lattice bosons, and the role of the NLDE in unitary Fermi gas analogs in general.
Alexander von Humboldt foundation and the Heidelberg Center for Quantum Dynamics for additional support. We acknowledge useful discussions with Ken O'Hara at Pennsylvania State University.
Appendix A. Convergence of numerical solutions of the quasi-1D reduction of the NLDE To compute the error of our oscillating numerical solutions, we first calculate the average difference between values of ψ A and ψ B at positions separated by one period. This tells us how the error in the periodicity of our solutions propagates with increasing position. The formula we use for the error ε(x) is where L is the periodicity for the particular solution and the error is computed for both two-spinor component functions ψ A and ψ B . In Fig. A1, we have plotted log 10 |ε(x)|, (a)-(d), and ε(x), (e)-(h), for the solution depicted in Fig. 5(a) using grid sizes N = 10 4 , 10 5 , 10 6 , 10 7 . The solution in Fig. 5(a) is obtained by using a forward stepping finite difference scheme. The initial conditions are obtained by choosing a value for f A which lies on the desired solution branch then substituting into the invariance relation Eq. (40) to obtain f B (0). Next, we check for convergence of the three lowest soliton solutions depicted in Fig. 15. The single soliton was obtained by finite differencing using a shooting method to tune the precision of the initial value of f A near f A ≈ 1 with higher precision forcing oscillations to x grid size. The double and triple soliton solutions are then found by tuning the chemical potential µ. For convergence at a single point we compute the solution at x i = x/ξ Dirac = 15 for several values of the grid size N = 10 1 , 10 2 , 10 3 , 10 4 , 10 5 , 10 6 . We use the formula for the error as a function of the number of grid points . In Fig. A2 we have plotted log 10 |ε(N)| versus log 10 N. The propagation error in the periodicity of the density shown in Fig. 15(f) is computed using Eq. (A.1) with results shown in Fig. A3.