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

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–8]. Recently, analytical solutions of the massive NLDE were obtained for the case of Kerr nonlinearity [9]. Dirac-like spin–orbit couplings for interacting cold atoms have also been investigated, simulating some features of quark confinement [10]. Moreover, solitons appear in systems with Dirac points such as quasi-one-dimensional (quasi-1D) nonlinear optical structures [11–21], acoustic physics [22], and electron propagation in graphene [23–26]. 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, 27]. We note that the (1+1)-dimensional NLDE has also been obtained starting from the nonlinear Schrödinger equation in a periodic potential using an asymptotic multi-scale expansion method [28]. Our own recent work has placed the NLDE in the context of a Bose–Einstein condensate (BEC) [29]. Significantly, our particular form of the NLDE has opened up research in other fields of physics [30–38]. 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.

For graphene nanoribbons the single-particle spectrum associated with the zigzag direction contains edge states as well as states confined to the interior of the ribbon, but none of these states have a linear degenerate band with Dirac-like dispersion. In contrast, for certain nanoribbon widths the spectrum of the armchair ribbon does contain Dirac points [23, 39, 40]. The spectra and single-particle states for armchair and zigzag ribbons have been obtained using both a tight-binding Schrödinger calculation [41, 42], and by starting from the Dirac Hamiltonian for the long-wavelength limit of the 2D lattice [23, 40]. Both methods agree quite well in the intermediate to large ribbon width range, i.e., for $N\gtrsim 10$ where N is the number of zigzag or armchair lines. However, Dirac points only occur in the spectrum for the armchair case and here only for particular values of N. By imposing quantization conditions at the edges of the ribbon one obtains a restriction on N to integer multiples of 3. Thus, for fixed values of the lattice constant one identifies specific ribbon widths that possess the dispersion of interest to us. This is similar to the case of graphene nanotubes where different $(m,n)$ combinations have different properties due to different 1D cuts through the Dirac cone (some are semiconducting, others semimetals). In practice, the quasi-1D NLDE is obtained by isolating the armchair direction in the full two-dimensional (2D) honeycomb lattice theory [29]. 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. The Dirac theory of the full 2D lattice is modified by introducing armchair boundary conditions for the short direction. This leaves a translationally invariant theory in the long direction described by either a massive or massless 1D NLDE (armchair NLDE) determined by the particular subband [40, 41, 43]. In our case, we focus on confinement in the zigzag direction, i.e., transverse to the armchair pattern of the 2D lattice. A schematic of the harmonic magnetic trap, interfering lasers, and BEC required to realize our set up is shown in figure 1. Details of the experimental construction can be found in [44].

Zoom In Zoom Out Reset image size

There are two equivalent forms of the quasi-1D reduction of the NLDE corresponding to a real or complex projection of the Dirac operator in one spatial dimension. Consequently, spinor solutions associated with these two projections are related by a complex Pauli matrix rotation. To obtain soliton solutions we first integrate the armchair NLDE to obtain an invariance relation which describes solutions at fixed values of the diagonal spatial element of the energy–momentum, 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, 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 $\mu /U={(\mu /U)}_{\mathrm{SB}}$, where μ is the chemical potential of the system and U is the quasi-1D renormalized interaction. Tuning $\mu /U$ towards ${(\mu /U)}_{\mathrm{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, as we will explain in section 5.

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 two-spinor. 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 armchair NLDE with real Dirac operator but extension of our results to the complex form via the aforementioned Pauli rotation is straightforward. 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.

The work that we present here is related to a number of parallel studies in condensed matter and cold atomic gases, among other contexts. The dimensional reduction of the quasi-2D honeycomb lattice to a quasi-1D lattice provides a novel way to study BECs. 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 [10]. 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 elativistic models such as the massive Thirring and Gross–Neveu models [10, 45–48]. 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 ten orders of magnitude slower than the speed of light [44].

This article is organized as follows. In section 2, we provide an introduction to the NLDE and discuss the key physical parameters. In section 3, we explain how the NLDE armchair geometry is 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 section 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 is retained as a subject of future work. Focusing on the time-independent armchair 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 section 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 section 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 section 7, we obtain solitons using a numerical shooting method. Finally, in section 8 we conclude.

In this section we introduce the NLDE, a nonlinear extension of the massless Dirac equation, and discuss the experimentally relevant physical parameters. The NLDE for two inequivalent Dirac points describes the dynamics of a Dirac four-spinor of the form $\Psi \equiv {\left({\Psi }_{+},{\Psi }_{-}\right)}^{{\mathsf{T}}}$, with the upper (+) and lower (−) two-spinors relating to opposite ${\bf{K}}$ and ${{\bf{K}}}^{\prime }$ points of the honeycomb lattice (see [25, 26, 29, 49]). 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., ${\rm{E}}({\rm{k}})\approx {\rm{\hslash }}v{\rm{k}}$, a consequence of the underlying symmetry of the honeycomb lattice. For graphene the proportionality constant v is just the Fermi velocity v_F. In BECs v is the quasi-particle group velocity c_l, and required to be less than the speed of sound in order to satisfy the Landau criterion. 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 ${\Psi }_{+}\equiv {\left({\psi }_{A+},\ {\psi }_{B+}\right)}^{{\mathsf{T}}}$ and ${\Psi }_{-}\equiv {\left({\psi }_{B-},\ {\psi }_{A-}\right)}^{{\mathsf{T}}}$. The full NLDE in this case is

$$\begin{eqnarray}&&{\rm{i}}{\rm{\hslash }}{\gamma }^{\mu }{\partial }_{\mu }\Psi +{U}_{2{\rm{D}}}\underset{i=1}{\overset{4}{\displaystyle \sum }}({\Psi }^{\dagger }{{\rm{M}}}_{i})(\Psi {{\rm{M}}}_{i})\;\Psi =0.\end{eqnarray} \tag{ 1 }$$

The matrices ${\gamma }^{\mu }$ are the usual Dirac matrices and the interaction terms are encapsulated in the summation with the matrices ${{\rm{M}}}_{i}$ constructed to give the correct cubic nonlinearities, local to each spinor component [29]. Explicitly, the interaction matrices are

$$\begin{eqnarray}{{\rm{M}}}_{1}=\left(\begin{array}{cccc}1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\end{array}\right),\ \ \ {{\rm{M}}}_{2}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\end{array}\right),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}{{\rm{M}}}_{3}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\end{array}\right),\ \ \ {{\rm{M}}}_{4}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1\end{array}\right).\end{eqnarray} \tag{ 3 }$$

In this simplest version of the NLDE for a BEC in a honeycomb lattice, the interactions do not couple different spinor components, which are only coupled through the kinetic term.³ Thus, in full generality we can focus on the equations for a two-spinor in rectangular coordinates while omitting the Dirac point subscript

$$\begin{eqnarray}&&-{\rm{i}}{\rm{\hslash }}{c}_{l}\left({\partial }_{x}-{\rm{i}}{\partial }_{y}\right){\psi }_{B}+{U}_{2{\rm{D}}}{\left|{\psi }_{A}\right|}^{2}{\psi }_{A}={\rm{i}}{\rm{\hslash }}\;{\partial }_{t}{\psi }_{A},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&-{\rm{i}}{\rm{\hslash }}{c}_{l}\left({\partial }_{x}+{\rm{i}}{\partial }_{y}\right){\psi }_{A}+{U}_{2{\rm{D}}}{\left|{\psi }_{B}\right|}^{2}{\psi }_{B}={\rm{i}}{\rm{\hslash }}\;{\partial }_{t}{\psi }_{B},\end{eqnarray} \tag{ 5 }$$

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}_{2{\rm{D}}}$.

Equations (4) and (5) allow for quasi-1D solutions by confining the BEC in one of the planar directions. From here on we will assume confinement in the y-direction. The experimental construction is shown in figure 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}^{{\rm{o}}}$ in the x y plane. Confinement in the y-direction produces an optical lattice with ribbon geometry wherein ${\psi }_{A}$ and ${\psi }_{B}$ are effectively functions only of x, for energies small compared to the characteristic energy associated with the width. A full derivation of the dimensional reduction from the 2D NLDE to the quasi-1D form is presented in section 3. The stationary states of interest are then obtained by taking the time-dependence to be the usual exponential factor with the chemical potential μ as the frequency: ${\psi }_{A}(x,t)=\;\mathrm{exp}(-{\rm{i}}\mu t/{\rm{\hslash }}){f}_{A}(x)$ and ${\psi }_{B}(x,t)=\;\mathrm{exp}(-{\rm{i}}\mu t/{\rm{\hslash }}){f}_{B}(x)$. Equations (4) and (5) become

$$\begin{eqnarray}&&-{\rm{i}}{\rm{\hslash }}{c}_{l}\;{\partial }_{x}{f}_{B}(x)+{U}_{1{\rm{D}}}\;{| {f}_{A}(x)| }^{2}{f}_{A}(x)=\mu {f}_{A}(x),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&-{\rm{i}}{\rm{\hslash }}{c}_{l}\;{\partial }_{x}{f}_{A}(x)+{U}_{1{\rm{D}}}\;{| {f}_{B}(x)| }^{2}{f}_{B}(x)=\mu {f}_{B}(x).\end{eqnarray} \tag{ 7 }$$

The system, equations (6) and (7), is the time-independent quasi-1D NLDE for the complex form of the Dirac operator. Notice that taking ${f}_{A}\to {\rm{i}}{f}_{B}$ and ${f}_{B}\to {f}_{A}$ converts equations (6) and (7) to the form associated with the real Dirac operator, as can also be obtained by choosing the y-direction in equations (4) and (5). This transformation is equivalent to multiplication by a linear combination of Pauli matrices with an overall complex factor, $[(1+{\rm{i}})/2]({\sigma }_{x}-{\sigma }_{y})$. It is natural to retain the x notation when discussing either form of the NLDE; thus we write the real form as

$$\begin{eqnarray}&&-{\rm{\hslash }}{c}_{l}\;{\partial }_{x}{f}_{B}(x)+{U}_{1{\rm{D}}}\;{| {f}_{A}(x)| }^{2}{f}_{A}(x)=\mu {f}_{A}(x),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\rm{\hslash }}{c}_{l}\;{\partial }_{x}{f}_{A}(x)+{U}_{1{\rm{D}}}\;{| {f}_{B}(x)| }^{2}{f}_{B}(x)=\mu {f}_{B}(x).\end{eqnarray} \tag{ 9 }$$

In section 3, we will see that the interaction ${U}_{1{\rm{D}}}$ in equations (6)–(9) is the renormalized version of the corresponding 2D interaction ${U}_{2{\rm{D}}}$ in equations (4) and (5).

Finally, we provide a brief consideration of physical implementation of the NLDE in BECs, with more detailed treatment given in [44]. 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\sqrt{3}/2{\rm{\hslash }}$; and the atom–atom binary interaction strength, ${U}_{2{\rm{D}}}={L}_{z}g{\bar{n}}^{2}3\sqrt{3}{a}^{2}/8$, where ${U}_{2{\rm{D}}}$ is the 2D optical lattice renormalized version of the usual interaction $g=4\pi {{\rm{\hslash }}}^{2}{a}_{s}/M$. Appearing in these definitions are the average particle density $\bar{n}=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 ${t}_{h}\equiv 1.861{\left({V}_{0}/{E}_{R}\right)}^{1/4}{E}_{R}\mathrm{exp}\left(-1.582\sqrt{{V}_{0}/{E}_{R}}\right)$ [50], where V₀ 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-1D regime (which we will discuss in detail in section 3) are $T=8.0\;\mathrm{nK}$, $\bar{n}=1.5\times {10}^{9}\;{\mathrm{cm}}^{-3}$, ${L}_{x}\gg {L}_{y}\gg {L}_{z}=3.0\;\mu {\rm{m}}$, ${t}_{h}=4.31\;\mathrm{nK}=1.04\;\mathrm{kHz}$, ${V}_{0}/{E}_{R}=16$, ${U}_{2{\rm{D}}}=0.391\;\mathrm{pK}=0.051\;\mathrm{Hz}$, $a=0.55\;\mu {\rm{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}\mathrm{Rb}$, and the associated scattering length ${a}_{s}=5.77\;\mathrm{nm}$. A complete discussion of NLDE parameters and constraints can be found in [44].

The quasi-1D NLDE is physically realized by starting with the 2D honeycomb lattice, then adding harmonic confinement in the y-direction. For generality we will explain the construction of both the armchair and zigzag geometries. We then focus exclusively on the armchair case since here the spectrum contains the linear degenerate points crucial to our results. The potentials for the honeycomb lattice with harmonic confinement for the armchair and zigzag geometries are given by

$$\begin{eqnarray}{V}_{\mathrm{armchair}}({\bf{r}}) & = & -\displaystyle \frac{\alpha {E}_{0}^{2}}{4}\;\left(3+2\;\mathrm{cos}\left[({{\bf{k}}}_{1}-{{\bf{k}}}_{2})\cdot {\bf{r}}\right]+2\;\mathrm{cos}\left[({{\bf{k}}}_{2}-{{\bf{k}}}_{3})\cdot {\bf{r}}\right]\right.\\ & & +\left.2\;\mathrm{cos}\left[({{\bf{k}}}_{1}-{{\bf{k}}}_{3})\cdot {\bf{r}}\right]\right)+\displaystyle \frac{1}{2}M{\omega }_{y}^{2}{y}^{2},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}{V}_{\mathrm{zigzag}}({\bf{r}}) & = & -\displaystyle \frac{\alpha {E}_{0}^{2}}{4}\;\left(3+2\;\mathrm{cos}\left[({{\bf{k}}}_{2}-{{\bf{k}}}_{1})\cdot {\bf{r}}\right]+2\;\mathrm{cos}\left[({{\bf{k}}}_{1}-{{\bf{k}}}_{3})\cdot {\bf{r}}\right]\right.\\ & & +\left.2\;\mathrm{cos}\left[({{\bf{k}}}_{2}-{{\bf{k}}}_{3})\cdot {\bf{r}}\right]\right)+\displaystyle \frac{1}{2}M{\omega }_{y}^{2}{y}^{2},\end{eqnarray} \tag{ 11 }$$

where α is the polarizability of the atoms, E₀ is the electric field strength, ${\bf{r}}=(x,y)$ is the planar coordinate vector, M is the atomic mass of the constituent bosons, and ${\omega }_{y}$ is the frequency of the harmonic potential which adds the additional confinement. The wave vectors ${{\bf{k}}}_{1}$, ${{\bf{k}}}_{2}$, and ${{\bf{k}}}_{3}$ in equations (10) and (11) are defined as

$$\begin{eqnarray}&&{{\bf{k}}}_{1}=k\;\hat{{\bf{y}}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{\bf{k}}}_{2}=\left(\displaystyle \frac{\sqrt{3}}{2}\hat{{\bf{x}}}-\displaystyle \frac{1}{2}\hat{{\bf{y}}}\right)k,\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\bf{k}}}_{3}=\left(-\displaystyle \frac{\sqrt{3}}{2}\hat{{\bf{x}}}-\displaystyle \frac{1}{2}\hat{{\bf{y}}}\right)k.\end{eqnarray} \tag{ 14 }$$

The different forms of V_armchair and V_zigzag in equations (10) and (11) reflect a uniform rotation of the beams by ${90}^{{\rm{o}}}$. The wavelength of the laser light that forms the lattice is $\lambda =2\pi /k$ which defines the distance between sites on a hexagonal sublattice as $a=2\lambda /3$. The underlying electric field which produces the potentials in equations (10) and (11) is the superposition of the fields from each beam and is given explicitly by

$$\begin{eqnarray}&&{\bf{E}}=\displaystyle \frac{{E}_{0}}{2}\hat{{\bf{z}}}\underset{i=1}{\overset{3}{\displaystyle \sum }}\mathrm{exp}\left[{\rm{i}}\left({{\bf{k}}}_{i}\cdot {\bf{r}}+{\omega }_{i}t\right)\right],\end{eqnarray} \tag{ 15 }$$

where ${\omega }_{i}=c| {{\bf{k}}}_{i}|$ ($c\;=$ speed of light $\approx 2.99\times {10}^{8}\;{\rm{m}}\;{{\rm{s}}}^{-1}$) here must be distinguished from the trapping frequency, and $\hat{{\bf{z}}}$ indicates the polarization of the beams perpendicular to the plane of the lattice. In the AC stark effect relevant to ultracold atoms trapped in optical lattices, the lattice potential is then obtained by taking the average of the square of the electric field

$$\begin{eqnarray}&&V=-\displaystyle \frac{1}{2}\alpha \bar{{| {\bf{E}}| }^{2}},\end{eqnarray} \tag{ 16 }$$

where V here can stand for either potential in equations (10) and (11). Plots of the armchair and zigzag potentials are shown in figures 1(c) and (d). For the remainder of our work we focus exclusively on the armchair configuration.

We now solve equations (4) and (5) in the presence of the armchair potential in equation (10). Stationary states are obtained as in the unconfined case but now with additional armchair boundary conditions in the y-direction for ${\psi }_{A}$ and ${\psi }_{B}$. Respecting translational invariance along the x-direction, separation of variables for the spinor functions gives ${\psi }_{A}(x,y,t)=\;\mathrm{exp}[{\rm{i}}({k}_{x}x-\mu t/{\rm{\hslash }})]\;\eta (y){g}_{A}(y)$ and ${\psi }_{B}(x,y,t)=\;\mathrm{exp}[{\rm{i}}({k}_{x}x-\mu t/{\rm{\hslash }})]$ $\;\eta (y){g}_{B}(y)$, where η is a real envelope function and ${g}_{A(B)}$ are superpositions of left and right plane wave excitations along the width of the ribbon. Equations (4) and (5) become

$$\begin{eqnarray}&&{\rm{\hslash }}{c}_{l}\left({k}_{x}-{\partial }_{y}\right){g}_{B}-{\rm{\hslash }}{c}_{l}{g}_{B}\displaystyle \frac{{\partial }_{y}\eta }{\eta }+{U}_{2{\rm{D}}}\;{\eta }^{2}{| {g}_{A}| }^{2}{g}_{A}+\displaystyle \frac{1}{2}M{\omega }_{y}^{2}{y}^{2}\;{g}_{A}=\mu \;{g}_{A},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\rm{\hslash }}{c}_{l}\left({k}_{x}+{\partial }_{y}\right){g}_{A}+{\rm{\hslash }}{c}_{l}{g}_{A}\displaystyle \frac{{\partial }_{y}\eta }{\eta }+{U}_{2{\rm{D}}}\;{\eta }^{2}{| {g}_{B}| }^{2}{g}_{B}+\displaystyle \frac{1}{2}M{\omega }_{y}^{2}{y}^{2}\;{g}_{B}=\mu \;{g}_{B}.\end{eqnarray} \tag{ 18 }$$

Our approach is to treat the slowly varying envelope function η using the Thomas–Fermi approximation and impose armchair boundary conditions at the ribbon edges for the plane wave functions. To implement this method we separate the chemical potential in terms of the envelope contribution ${f}_{e}\mu$ and the contribution from the plane-wave excitations ${f}_{p}\mu$, where ${f}_{e}+{f}_{p}=1$. The system equations (17) and (18) splits into two pairs of equations

$$\begin{eqnarray}&&-{\rm{\hslash }}{c}_{l}{g}_{B}\displaystyle \frac{{\partial }_{y}\eta }{\eta }+{U}_{2{\rm{D}}}\;{\eta }^{2}{| {g}_{A}| }^{2}{g}_{A}+\displaystyle \frac{1}{2}M{\omega }_{y}^{2}{y}^{2}\;{g}_{A}={f}_{e}\mu \;{g}_{A},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\rm{\hslash }}{c}_{l}{g}_{A}\displaystyle \frac{{\partial }_{y}\eta }{\eta }+{U}_{2{\rm{D}}}\;{\eta }^{2}{| {g}_{B}| }^{2}{g}_{B}+\displaystyle \frac{1}{2}M{\omega }_{y}^{2}{y}^{2}\;{g}_{B}={f}_{e}\mu \;{g}_{B},\end{eqnarray} \tag{ 20 }$$

and

$$\begin{eqnarray}&&{\rm{\hslash }}{c}_{l}\left({k}_{x}-{\partial }_{y}\right){g}_{B}={f}_{p}\mu \;{g}_{A},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{\rm{\hslash }}{c}_{l}\left({k}_{x}+{\partial }_{y}\right){g}_{A}={f}_{p}\mu \;{g}_{B}.\end{eqnarray} \tag{ 22 }$$

Applying the Thomas–Fermi approximation in equations (19) and (20) by neglecting derivatives of the envelope function in favor of the interaction and harmonic terms, we obtain

$$\begin{eqnarray}&&{\eta }^{2}(y)=\displaystyle \frac{{f}_{e}\mu -(1/2)M{\omega }_{y}^{2}{y}^{2}}{{U}_{2{\rm{D}}}{| {g}_{A}| }^{2}},\ \ \mathrm{and}\ \ \ {\eta }^{2}(y)=\displaystyle \frac{{f}_{e}\mu -(1/2)M{\omega }_{y}^{2}{y}^{2}}{{U}_{2{\rm{D}}}{| {g}_{B}| }^{2}},\end{eqnarray} \tag{ 23 }$$

which requires the condition $| {g}_{A}| =\pm | {g}_{B}|$, for consistency. We combine equations (21) and (22) to arrive at

$$\begin{eqnarray}&&{{\rm{\hslash }}}^{2}{c}_{l}^{2}\left({k}_{x}^{2}-{\partial }_{y}^{2}\right){g}_{B}={f}_{p}^{2}{\mu }^{2}\;{g}_{B},\ \ {g}_{A}=\displaystyle \frac{1}{{f}_{p}\mu }{\rm{\hslash }}{c}_{l}\left({k}_{x}-{\partial }_{y}\right){g}_{B},\end{eqnarray} \tag{ 24 }$$

which we then solve by the plane wave decomposition

$$\begin{eqnarray}&&{g}_{B}={b}_{1}{{\rm{e}}}^{{\rm{i}}{k}_{n}y}+{b}_{2}{{\rm{e}}}^{-{\rm{i}}{k}_{n}y},\end{eqnarray} \tag{ 25 }$$

with g_A determined by

$$\begin{eqnarray}&&{g}_{A}=\displaystyle \frac{{\rm{\hslash }}{c}_{l}}{{f}_{p}\mu }\sqrt{{k}_{x}^{2}+{k}_{n}^{2}}\left[{b}_{1}\;{{\rm{e}}}^{-{\rm{i}}\;\alpha }\;{{\rm{e}}}^{{\rm{i}}{k}_{n}y}+{b}_{2}\;{{\rm{e}}}^{{\rm{i}}\;\alpha }\;{{\rm{e}}}^{-{\rm{i}}{k}_{n}y}\right],\end{eqnarray} \tag{ 26 }$$

where the phase in equation (26) is defined as $\alpha \equiv {\mathrm{tan}}^{-1}({k}_{n}/{k}_{x})$. We have anticipated spatially quantized states along the width of the ribbon by including the integer subscript k_n, $n\in {\mathbb{Z}}$ in equation (26). Imposing armchair boundary conditions on ${g}_{A(B)}$ and their counterparts at the ${\bf{K}}\prime$ point leads to the conditions

$$\begin{eqnarray}&&{b}_{1}={b}_{2}^{\prime },\ \ \ \ {b}_{2}={b}_{1}^{\prime }=0,\end{eqnarray} \tag{ 27 }$$

and for the wavenumber along the width of the ribbon

$$\begin{eqnarray}&&\mathrm{sin}\left[\left({k}_{n}+K\right){L}_{y}\right]=0,\end{eqnarray} \tag{ 28 }$$

which gives the spectral condition

$$\begin{eqnarray}&&{k}_{n}=\displaystyle \frac{n\pi }{{L}_{y}}-\displaystyle \frac{4\pi }{3{a}_{0}}.\end{eqnarray} \tag{ 29 }$$

In terms of the number of dimers N (adjacent pairs between the ribbon edges) and the lattice spacing a₀, the width L_y of the armchair nanoribbon is ${L}_{y}={{Na}}_{0}/4$. Focusing attention on the Dirac points in the spectrum, which occur for ${k}_{n}=(4\pi /{a}_{0})[(n/N)-1/3]=0$, we find the constraint $N=3n$ relating the width integer N and the subband quantization number n [40]. Thus, Dirac points occur when N is an integer multiple of 3. Combining the constraint derived from equation (23), namely $| {g}_{A}| =\pm | {g}_{B}|$, with equations (25) and (26) and equation (27) gives us the dispersion relation

$$\begin{eqnarray}&&{f}_{p}\mu ={\rm{\hslash }}{c}_{l}\sqrt{{k}_{x}^{2}+{k}_{n}^{2}}.\end{eqnarray} \tag{ 30 }$$

To organize our results so far we first note that the construction of the optical lattice nanoribbon requires specifying the atomic polarizability of the atoms α, the electric field strength of the optical lattice E₀, and the harmonic trap frequency ${\omega }_{y}$. For a particular choice of ribbon width L_y (equivalently N) the envelope function equation (23) must vanish at the edges $\eta (\pm {L}_{y}/2)=0$, with the edges defined along the line where the lattice potential well depth equals the magnitude of the harmonic potential: ${V}_{\mathrm{lattice}}(\pm {L}_{y}/2)={V}_{\mathrm{trap}}(\pm {L}_{y}/2)$. The condition for vanishing envelope implies that

$$\begin{eqnarray}&&{\eta }^{2}(\pm {L}_{y}/2)=\displaystyle \frac{[(1-{f}_{p})/{f}_{p}]\;{\rm{\hslash }}{c}_{l}\;\sqrt{{k}_{x}^{2}+{k}_{n}^{2}}-(1/8)M{\omega }_{y}^{2}{L}_{y}^{2}}{{U}_{2{\rm{D}}}{| {b}_{1}| }^{2}}=0,\end{eqnarray} \tag{ 31 }$$

and the condition for the optical and harmonic potentials is

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}M{\omega }_{y}^{2}{L}_{y}^{2}=\alpha {E}_{0}^{2}.\end{eqnarray} \tag{ 32 }$$

Combining equations (31) and (32) determines the envelope and plane wave fractions f_e and f_p

$$\begin{eqnarray}&&{f}_{p}=1-{f}_{e}={\left(\displaystyle \frac{\alpha {E}_{0}^{2}}{4{\rm{\hslash }}{c}_{l}\sqrt{{k}_{x}^{2}+{k}_{n}^{2}}}+1\right)}^{-1},\end{eqnarray} \tag{ 33 }$$

which combines with equation (30) to give the full expression for the dispersion

$$\begin{eqnarray}&&\mu ={\rm{\hslash }}{c}_{l}\sqrt{{k}_{x}^{2}+{k}_{n}^{2}}+\displaystyle \frac{M{\omega }_{y}^{2}{L}_{y}^{2}}{8}.\end{eqnarray} \tag{ 34 }$$

The first term in equation (34) arises from excitations in the long and short directions of the ribbon and the second term accounts for the finite width of the ribbon, which in practicality may be subtracted off as an overall constant energy by defining the shifted chemical potential $\mu \to \mu -M{\omega }_{y}^{2}{L}_{y}^{2}/8$. Upon inclusion of the ${\bf{K}}\prime$-point contribution, the 2-component spinor wavefunction is given by

$$\begin{eqnarray}&&\Psi (x,y,t)\;=\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{{\rm{e}}}^{{\rm{i}}\left({k}_{x}x-\mu \;t/{\rm{\hslash }}\right)}\;\sqrt{\displaystyle \frac{M{\omega }_{y}^{2}{L}_{y}^{2}}{2{U}_{2{\rm{D}}}}\left(1-4\displaystyle \frac{{y}^{2}}{{L}_{y}^{2}}\right)}\left(\begin{array}{c}\mathrm{cos}\left({k}_{n}\;y\right)\\ \pm \displaystyle \frac{{\rm{i}}}{\sqrt{{k}_{x}^{2}+{k}_{n}^{2}}}\left[{k}_{x}\mathrm{cos}\left({k}_{n}\;y\right)+{k}_{n}\mathrm{sin}\left({k}_{n}\;y\right)\right]\end{array}\right),\end{eqnarray} \tag{ 36 }$$

defined over the width $-{L}_{y}/2\lt y\lt +{L}_{y}/2$. In particular, for the linear degenerate subband where ${k}_{n}=0$ the wavefunction reduces to

$$\begin{eqnarray}&&\Psi (x,y,t)={{\rm{e}}}^{{\rm{i}}\left({k}_{x}x-\mu \;t/{\rm{\hslash }}\right)}\;\sqrt{\displaystyle \frac{M{\omega }_{y}^{2}{L}_{y}^{2}}{2{U}_{2{\rm{D}}}}\left(1-4\displaystyle \frac{{y}^{2}}{{L}_{y}^{2}}\right)}\left(\begin{array}{c}1\\ \pm {\rm{i}}\end{array}\right).\end{eqnarray} \tag{ 37 }$$

In the remainder of this section we use our results up to this point to derive the quasi-1D NLDE, equations (6) and (7) and equations (8) and (9), starting from the 2D NLDE in equations (4) and (5). Transforming to the quasi-1D regime requires that ${L}_{y}\ll {L}_{x}$, which ensures that excitations along the x-direction have much lower energy than those in the y-direction. In particular, this condition must be satisfied for the armchair nanoribbon geometry discussed thus far. Equivalently, these constraints may be expressed in terms of the trap frequencies and the atom-atom interaction: ${\omega }_{x}\ll {\omega }_{y}$, ${\rm{\hslash }}{\omega }_{x}\ll U$. For the moment we neglect the harmonic trap in x. A modified renormalized interaction U and careful consideration of phase coherence in quasi-1D BECs should be sufficient to account for tight confinement, but are not necessary here [51–54].

The quasi-1D NLDE is obtained by separating ${\psi }_{A}({\bf{r}},t)$ and ${\psi }_{B}({\bf{r}},t)$ into longitudinal and transverse modes following similar arguments as in [55]:

$$\begin{eqnarray}&&{\psi }_{A}({\bf{r}},t)={\left(\displaystyle \frac{M{\omega }_{y}^{2}{L}_{y}^{3}}{3{U}_{2D}}\right)}^{-1/2}h(y){f}_{A}(x){{\rm{e}}}^{-{\rm{i}}\mu t/{\rm{\hslash }}},\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{\psi }_{B}({\bf{r}},t)={\left(\displaystyle \frac{M{\omega }_{y}^{2}{L}_{y}^{3}}{3{U}_{2D}}\right)}^{-1/2}h(y){f}_{B}(x){{\rm{e}}}^{-{\rm{i}}\mu t/{\rm{\hslash }}},\end{eqnarray} \tag{ 39 }$$

where ${f}_{A(B)}(x)$ contains the longitudinal x-dependence and $h(y)$ is the dimensionless function that describes the transverse part of the wavefunction from equation (37), i.e.,

$$\begin{eqnarray}&&h(y)=\sqrt{\displaystyle \frac{M{\omega }_{y}^{2}{L}_{y}^{2}}{2{U}_{2{\rm{D}}}}\left(1-4\displaystyle \frac{{y}^{2}}{{L}_{y}^{2}}\right)}.\end{eqnarray} \tag{ 40 }$$

Note that we have included a normalization prefactor in equations (38) and (39). Substituting equations (38) and (39) into equations (4) and (5) using the expression for $h(y)$ in equation (40) and integrating over the y-direction leads directly to the quasi-1D NLDE with complex coefficients in equations (6) and (7). Consequently, we find the quasi-1D renormalized interaction to be

$$\begin{eqnarray}&&{U}_{1{\rm{D}}}\equiv {U}_{2{\rm{D}}}{\left(\displaystyle \frac{3}{2{L}_{y}}\right)}^{2}{\displaystyle \int }_{-{L}_{y}/2}^{+{L}_{y}/2}{\rm{d}}y\;\left(1-4\displaystyle \frac{{y}^{2}}{{L}_{y}^{2}}\right)\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&=\;\left(\displaystyle \frac{6}{5{L}_{y}}\right){U}_{2{\rm{D}}}.\end{eqnarray} \tag{ 42 }$$

A key result here is that the chemical potential is not modified by dimensional reduction since the Dirac equation is first order in the spatial derivatives and the extra term proportional to dh/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-1D renormalized interaction.

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 equations (8) and (9) we look for real solutions and write the NLDE as a derivative field (or direction field)

$$\begin{eqnarray}&&\quad {f}_{B}^{\prime }=-\displaystyle \frac{U}{{\rm{\hslash }}{c}_{l}}{f}_{A}\left(\displaystyle \frac{\mu }{U}-{f}_{A}^{2}\right),\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{f}_{A}^{\prime }=\displaystyle \frac{U}{{\rm{\hslash }}{c}_{l}}{f}_{B}\left(\displaystyle \frac{\mu }{U}-{f}_{B}^{2}\right),\end{eqnarray} \tag{ 44 }$$

where the dependence on x is implied. Together, equations (43) and (44) comprise a vector field $({f}_{A}^{\prime },{f}_{B}^{\prime })$ which describes the flow of two-spinor solutions. The various combinations of conditions on the signs of ${f}_{A}^{\prime }$ and ${f}_{B}^{\prime }$ partition the $({f}_{A},{f}_{B})$ solution space into 16 regions. However, analysis of symmetries of equations (43) and (44) shows that only eight combinations lead to distinct solution types. In particular, the transformation ${f}_{A}\to -{f}_{B},\ {f}_{B}\to {f}_{A}$ leaves our equations invariant. We have listed these regions in table 1 along with the corresponding signs for the derivatives ${f}_{A}^{\prime }$ and ${f}_{B}^{\prime }$. Nine fixed points exist in the solution space $({f}_{A},{f}_{B})$: $(0,0),(\pm \sqrt{\mu /U},0),(\pm \sqrt{\mu /U},\pm \sqrt{\mu /U}),(0,\pm \sqrt{\mu /U})$.

Table 1.  Distinct solution regions of the NLDE. The sign of each derivative is determined by the values of both spinor functions. Note that in each region the signs of ${f}_{A}^{\prime }$ and ${f}_{B}^{\prime }$ remain fixed, where a sign change occurs across a boundary. The boundary between two regions is defined by ${f}_{A}^{\prime }=0$ or ${f}_{B}^{\prime }=0$, which corresponds to replacing an inequality by equality in any of the conditions listed above. The far right column lists the three types of soliton solutions and their associated regions.

Region	Condition for fA	Condition for fB	${f}_{A}^{\prime }$	${f}_{B}^{\prime }$	Solution type
I	$\sqrt{\mu /U}\lt {f}_{A}$	$\sqrt{\mu /U}\lt {f}_{B}$	−	+	Bright, multi-soliton
II	$\sqrt{\mu /U}\lt {f}_{A}$	$0\lt {f}_{B}\lt \sqrt{\mu /U}$	+	+	Bright, multi-soliton
III	$\sqrt{\mu /U}\lt {f}_{A}$	$-\sqrt{\mu /U}\lt {f}_{B}\lt 0$	−	+	Multi-soliton
IV	$\sqrt{\mu /U}\lt {f}_{A}$	${f}_{B}\lt -\sqrt{\mu /U}$	+	+	Multi-soliton
V	$0\lt {f}_{A}\lt \sqrt{\mu /U}$	$\sqrt{\mu /U}\lt {f}_{B}$	−	−	Bright, multi-soliton
VI	$0\lt {f}_{A}\lt \sqrt{\mu /U}$	$0\lt {f}_{B}\lt \sqrt{\mu /U}$	+	−	Dark soliton
VII	$0\lt {f}_{A}\lt \sqrt{\mu /U}$	$-\sqrt{\mu /U}\lt {f}_{B}\lt 0$	−	−	—
VIII	$0\lt {f}_{A}\lt \sqrt{\mu /U}$	${f}_{B}\lt -\sqrt{\mu /U}$	+	−	Multi-soliton

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 (periodic) solutions exist for which f_A and f_B both have amplitudes greater than $\sqrt{\mu /U}$, in which case each component has three critical points (${f}_{A(B)}^{\prime }=0$) during a half period. Such solutions cover all regions in table 1 except regions VI and VII. We give an example of this type of solution in this section, which we obtain analytically. Other oscillating solutions occur around each of the fixed points $(\pm \sqrt{\mu /U},0)$, which may be discerned from equations (43) and (44) using linear perturbation theory by substituting ${f}_{A}(x)=\epsilon \;\mathrm{cos}{kx}\pm \sqrt{\mu /U}$ or ${f}_{A}(x)=\epsilon \;\mathrm{cos}{kx}$, with the small amplitude such that $\epsilon /\sqrt{\mu /U}\ll 1$. These solutions cross four region boundaries. For example, regions II, III, VI, and VII in the case where ${f}_{A}(x)=\epsilon \;\mathrm{cos}{kx}\pm \sqrt{\mu /U}$, with each derivative ${f}_{A}^{\prime }$ and ${f}_{B}^{\prime }$ changing sign once per half period. We also find non-oscillating asymptotically flat solutions having a similar form to the function $(1/2)[1\pm \mathrm{tanh}(x)]$, which we study in section 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 which crosses regions I, II, and V. 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.

Zoom In Zoom Out Reset image size

To deepen our analysis, a detailed map of the solution space of equations (43) and (44) can be arrived at by uncovering spatially invariant quantities. To obtain the first invariant quantity we multiply equation (43) by the right hand side of equation (44) and vice versa, then add the resulting equations to obtain

$$\begin{eqnarray}&&\left(\mu {f}_{A}-U\;{f}_{A}^{3}\right){f}_{A}^{\prime }+\left(\mu {f}_{B}-{{Uf}}_{B}^{3}\right){f}_{B}^{\prime }=0,\end{eqnarray} \tag{ 45 }$$

which then gives

$$\begin{eqnarray}&&{\left[\left(2\mu {f}_{A}^{2}-{{Uf}}_{A}^{4}\right)+\left(2\mu {f}_{B}^{2}-{{Uf}}_{B}^{4}\right)\right]}^{\prime }=0,\end{eqnarray} \tag{ 46 }$$

where the prime notation indicates differentiation with respect to x. Integrating equation (46) gives a relation between the functions f_A and f_B

$$\begin{eqnarray}&&{\left({f}_{A}^{2}-\displaystyle \frac{\mu }{U}\right)}^{2}+{\left({f}_{B}^{2}-\displaystyle \frac{\mu }{U}\right)}^{2}-2{\left(\displaystyle \frac{\mu }{U}\right)}^{2}=C,\end{eqnarray} \tag{ 47 }$$

where we have simplified the expression by completing the squares in f_A and f_B, and C is the integration constant. The meaning of equation (46) is seen by multiplying equation (43) by f_B and equation (44) by f_A and then adding the resulting equations, which gives

$$\begin{eqnarray}&&-\displaystyle \frac{U}{2}\left[\left({f}_{A}^{4}+{f}_{B}^{4}\right)-2\displaystyle \frac{\mu }{U}\left({f}_{A}^{2}+{f}_{B}^{2}\right)\right]={\rm{\hslash }}{c}_{l}\left({f}_{A}^{\prime }{f}_{B}-{f}_{A}{f}_{B}^{\prime }\right)+\displaystyle \frac{U}{2}\left({f}_{A}^{4}+{f}_{B}^{4}\right).\end{eqnarray} \tag{ 48 }$$

The expression in brackets on the left of equation (48) is the same as the left hand side of equation (47), while the right side of equation (48) is the ${T}^{11}$ element of the energy–momentum tensor ${T}^{\mu \nu }$ [56]. Thus, we find

$$\begin{eqnarray}&&{T}^{11}=-\displaystyle \frac{{UC}}{2},\end{eqnarray} \tag{ 49 }$$

so that equation (47) is a statement of uniformity of pressure ${T}^{11}/{c}_{l}$ in the x-direction. The other elements of the energy–momentum tensor can also be computed, whereby we obtain the total energy density ${T}^{00}=-({\rm{\hslash }}{c}_{l}/2)({f}_{A}^{\prime }{f}_{B}-{f}_{A}{f}_{B}^{\prime })+(U/2)({f}_{A}^{4}+{f}_{B}^{4})$ and current ${T}^{01}={T}^{10}=({\rm{i}}/2)({\Psi }^{\dagger }{\Psi }^{\prime }-{{\Psi }^{\dagger }}^{\prime }\Psi )=0$. Note that the current vanishes but the energy density is spatially dependent. Furthermore, the energy density can be expressed in terms of the pressure ${T}^{11}$ and the interaction energy as ${T}^{00}=-(1/2){T}^{11}+(3/4)U({f}_{A}^{4}+{f}_{B}^{4})$.

A second quantity which characterizes solutions of the NLDE is arrived at by expressing equations (43) and (44) in the form

$$\begin{eqnarray}&&{f}_{B}{f}_{B}^{\prime }=-\displaystyle \frac{U}{{\rm{\hslash }}{c}_{l}}{f}_{B}{f}_{A}\left(\displaystyle \frac{\mu }{U}-{f}_{A}^{2}\right),\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{f}_{A}{f}_{A}^{\prime }=\displaystyle \frac{U}{{\rm{\hslash }}{c}_{l}}{f}_{A}{f}_{B}\left(\displaystyle \frac{\mu }{U}-{f}_{B}^{2}\right),\end{eqnarray} \tag{ 51 }$$

then adding, integrating, and combining terms to get a total derivative on the left hand side:

$$\begin{eqnarray}&&{\left({f}_{A}^{2}+{f}_{B}^{2}\right)}^{\prime }=\displaystyle \frac{2U}{{\rm{\hslash }}{c}_{l}}{f}_{B}{f}_{A}\left({f}_{A}^{2}-{f}_{B}^{2}\right).\end{eqnarray} \tag{ 52 }$$

Equation (52) can be simplified by introducing the average spin components ${S}_{x}\equiv \bar{\Psi }{\sigma }_{x}\Psi =2{f}_{A}{f}_{B}$, ${S}_{z}\equiv \bar{\Psi }{\sigma }_{z}\Psi =2{f}_{A}^{2}-{f}_{B}^{2}$, and the total density $\rho \equiv {f}_{A}^{2}+{f}_{B}^{2}$, so that equation (52) becomes ${\rho }^{\prime }=\left(U/{\rm{\hslash }}{c}_{l}\right){S}_{x}{S}_{z}$. Equation (52) states that the total density varies most where the wavefunction lies between a pure chiral state and a highly mixed state, as S_z is a measure of chirality and S_x a measure of the degree of chiral mixing within a particular state. Integrating equation (52) over an interval $a\lt x\lt b$ gives

$$\begin{eqnarray}&&{\displaystyle \int }_{a}^{b}{\rm{d}}x\;{S}_{x}{S}_{z}=\displaystyle \frac{{\rm{\hslash }}{c}_{l}}{U}[\rho (b)-\rho (a)].\end{eqnarray} \tag{ 53 }$$

In particular, equation (53) allows for solutions which asymptotically approach a constant value, i.e.,

$$\begin{eqnarray}&&{\displaystyle \int }_{-\infty }^{+\infty }{\rm{d}}x\;{S}_{x}{S}_{z}={\delta }_{\pm \infty },\end{eqnarray} \tag{ 54 }$$

where the constant ${\delta }_{\pm \infty }\equiv ({\rm{\hslash }}{c}_{l}/U)[\rho (+\infty )-\rho (-\infty )]$ is a global parameter that depends only on the difference between the asymptotic values of the total density. Hence, stationary solitons are topologically stable in that local fluctuations in f_A and f_B will cancel out, with ${\delta }_{\pm \infty }$ remaining fixed.

We now combine the qualitative information and the invariance relation (equation (47)) from section 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, equations (43) and (44) 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 equations (43) and (44) by plotting equation (47) for a particular value of ${T}^{11}$. This allows us to see how solutions evolve as we vary the ratio $\mu /U$. For the choice ${T}^{11}={\mu }^{2}/2U$, which corresponds to $C=-{(\mu /U)}^{2}$ in equation (47), we have plotted the solution space for $\mu /U=0.8$, 0.9, 0.98, 1.0, 1.02, 1.25, 1.5, 1.75, 2.0, in figures 3(a)–(i). In figures 3(a)–(c), solutions oscillate about the fixed points $({f}_{A},{f}_{B})=(\pm \sqrt{\mu /U},\pm \sqrt{\mu /U})$ with the orbits beginning to coalesce in figure 3(c). In figure 3(d), $\mu /U=1$ and solution paths begin and end at the saddle fixed points $(0,\pm \sqrt{\mu /U}),(\pm \sqrt{\mu /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 $\mu /U$ is increased from unity, orbits bifurcate into solutions which oscillate about $(0,0)$, with small and large amplitudes indicated by the paths 1 and 2 in figure 3(e). These orbits continue to smooth out through figure 3(i), at which point $\mu /U=2$. Soliton solutions such as those depicted in figure 3(d) correspond to the case ${T}^{11}={\mu }^{2}/2U$. In general though, isolated dark and bright single or multi-solitons are distinguished by the density conditions

$$\begin{eqnarray}&&0\lt {\rho }_{\mathrm{DS}}\lt \displaystyle \frac{\mu }{U}\;\bar{n},\displaystyle \frac{\mu }{U}\;\bar{n}\lt {\rho }_{\mathrm{BS}}\lt \displaystyle \frac{\mu }{U}\left[1+\sqrt{2\left(1-\displaystyle \frac{{T}^{11}U}{{\mu }^{2}}\right)}\right]\bar{n}.\end{eqnarray} \tag{ 55 }$$

The upper and lower bounds $\mu \bar{n}/U$ result from equations (43) and (44), and the upper bound on the second inequality comes from equations (47) and (49). We have added the subscripts DS and BS to the total density in equation (55) to indicate dark and bright solitons, respectively. Note the inclusion of the average particle density $\bar{n}$. The densities are defined in terms of the real spinor spatial functions ${\rho }_{\mathrm{DS}}(x)=\bar{n}\left[{f}_{A,\mathrm{DS}}^{2}(x)+{f}_{B,\mathrm{DS}}^{2}(x)\right]$, with an analogous definition for the bright soliton.

Zoom In Zoom Out Reset image size

The bounds for the inequalities in equation (55) can be proved as follows. For the proof of the upper bound of the first inequality, we start by assuming that ${\rho }_{\mathrm{DS}}({x}_{1})\geqslant \bar{n}\mu /U$ for at least one element in the domain ${x}_{1}\in {\mathbb{R}}$, and where $\bar{n}\mu /U\gt 0$. Since we are considering isolated dark solitons, there exist an infinite number of points ${x}_{2}\in {\mathbb{R}}$ such that ${\rho }_{\mathrm{DS}}({x}_{1})\lt {\rho }_{\mathrm{DS}}({x}_{2})$. Moreover, we have ${\mathrm{lim}}_{x\to \pm \infty }\;{\rho }_{\mathrm{DS}}(x)={\rho }_{\mathrm{DS}}^{\mathrm{sup}}\equiv \mathrm{sup}\left\{{\rho }_{\mathrm{DS}}(x)\in {\mathbb{R}}:x\in {\mathbb{R}})\right\}$, where ${\rho }_{\mathrm{DS}}^{\mathrm{sup}}$ denotes the supremum of the dark soliton density. Thus, $\bar{n}\mu /U\lt {\rho }_{\mathrm{DS}}^{\mathrm{sup}}$ and ${\mathrm{lim}}_{\;x\to \pm \infty }\;{\rho }_{\mathrm{DS}}^{\prime }(x)=0$. It follows that ${\mathrm{lim}}_{\;x\to \pm \infty }\left({f}_{A,\mathrm{DS}}^{2}+{f}_{B,\mathrm{DS}}^{2}\right)={\rho }_{\mathrm{DS}}^{\mathrm{sup}}/\bar{n}$ $\Rightarrow$ $-\sqrt{{\rho }_{\mathrm{DS}}^{\mathrm{sup}}/\bar{n}}\lt {f}_{A(B),\mathrm{DS}}\lt +\sqrt{{\rho }_{\mathrm{DS}}^{\mathrm{sup}}/\bar{n}}$, and that ${\mathrm{lim}}_{\;x\to \pm \infty }{f}_{A(B),\mathrm{DS}}^{\prime }(x)=0$. Here the subscript indicates that the condition applies to both f_A and f_B. By equations (43) and (44), asymptotically vanishing derivatives imply four possible combinations: ${\mathrm{lim}}_{\;x\to \pm \infty }\left\{{f}_{A,\mathrm{DS}}^{2}(x),{f}_{B,\mathrm{DS}}^{2}(x)\right\}=\{0,0\},\{\mu /U,0\},\{0,\mu /U\},\{\mu /U,\mu /U\}$. The first three cases lead to ${\rho }_{\mathrm{DS}}^{\mathrm{sup}}=0,\bar{n}\mu /U,\bar{n}\mu /U$, respectively, which contradict our earlier result that $0\lt \bar{n}\mu /U\lt {\rho }_{\mathrm{DS}}^{\mathrm{sup}}$. The fourth combination cannot occur as one may deduce from equations (43) and (44) and the analysis in table 1, proving by contradiction the upper bound of the first inequality in equation (55). Proving the lower bound in the second inequality proceeds by the reverse argument, i.e., using the infimum ${\rho }_{\mathrm{BS}}^{\mathrm{inf}}$ and the initial assumption that ${\rho }_{\mathrm{DS}}({x}_{1})\leqslant \bar{n}\mu /U$, but otherwise the steps are similar to those in the first proof. Finally we address the upper bound in the second inequality in equation (55). This bound comes from the invariance relation equation (47) where one sees that the sum of squared terms places an upper bound on ${f}_{A,\mathrm{BS}}(x)$ and ${f}_{B,\mathrm{BS}}(x)$, and thus on the total density ${\rho }_{\mathrm{BS}}(x)$. We obtain this result by setting ${f}_{A,\mathrm{BS}}(x)$ or ${f}_{B,\mathrm{BS}}(x)$ equal to zero in equation (47) and using equation (49) to replace C by $-{T}^{11}U/2$. Note that the radical in equation (55) implies the upper bound ${T}_{\mathrm{max}}^{11}={\mu }^{2}/U$.

Thus, solitons exist on a three-dimensional sub-manifold of parameters (defined by the condition $\mu ={| 2{T}^{11}U| }^{1/2}$) of the four-dimensional parameter manifold determined by the chemical potential, interaction, density, and energy–momentum tensor with coordinates denoted as $(\mu ,U,\rho ,{T}^{11})$. Moreover, the density conditions in equation (55) further partition the 3D parameter subspace into the two types of solitons along the boundary $\rho =\bar{n}{| 2{T}^{11}/U| }^{1/2}=\bar{n}\mu /U$, where the last equality pertains to the spacial case ${T}^{11}={\mu }^{2}/2U$.

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 ${T}^{11}=\mu /2$ the soliton boundary occurs at ${(\mu /U)}_{\mathrm{SB}}=1$. Tuning $\mu /U\to 1$ while keeping $\rho \lt \bar{n}$ forces the solution to the dark soliton, whereas tuning $\mu /U\to 1$ while keeping $\bar{n}\lt \rho \lt 2\bar{n}$ converges on the bright soliton. It is important to keep in mind that the upper bound here, $\rho \lt 2\bar{n}$, comes from choosing a particular value for ${T}^{11}$ and that bright solitons exist at higher densities but are associated with a different choice of ${T}^{11}$. Note that the type of soliton obtained is independent of whether $\mu /U$ approaches 1 from above or below. Conversely, if we maintain the condition $\mu /U=1$ while tuning ρ through the critical value ${\rho }_{c}=\bar{n}$ we induce a transition between the dark soliton and bright soliton. Thus, tuning $\mu /U$ moves the system between oscillating regimes across the soliton boundary, while tuning $\rho /\bar{n}$ along ${(\mu /U)}_{\mathrm{SB}}$ moves the system between the dark and bright solitons. Quantum phase transitions across the soliton boundary are summarized in figure 4. This kind of phase transition in the mean-field theory indicates a possible corresponding quantum phase transition in the underlying microscopic theory [57–60]. 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 [61].

Zoom In Zoom Out Reset image size

To make our analysis more concrete, we solve the numerical initial value problem defined by equations (43) and (44) with the initial conditions taken from equation (47), which relates f_A and f_B at x = 0. The value of ${f}_{A}(0)$ is chosen to coincide with a particular branch and a three-point balanced finite difference scheme for the first derivative is implemented. In figures 5(a)–(i), we have plotted the evolution of solutions corresponding to the upper branch (labeled as branch 2) of the orbits in figure 3 using the initial values $({f}_{A}(0),{f}_{B}(0))=(0.8000,1.1421)$, $(0.8000,1.2819)$, $(0.8000,1.3702)$, $(1.0010,0.0014)$, $(1.0100,0.0000)$, $(1.1000,1.6454)$, $(1.1000,1.8298)$, $(1.1000,1.9871)$, $(1.1000,2.1272)$. The values of ${f}_{A}(0)$ were chosen to pick out branch 2 with ${f}_{B}(0)$ determined by inverting and solving equation (47) for f_B

$$\begin{eqnarray}&&{f}_{B}(0)=\pm {\left[\mu /U\pm \sqrt{2{T}^{11}/U+2{(\mu /U)}^{2}-{\left({f}_{A}{(0)}^{2}-\mu /U\right)}^{2}}\right]}^{1/2}.\end{eqnarray} \tag{ 56 }$$

Similar plots focusing on branch 1 are shown in figure 6 using the initial values: $({f}_{A}(0),{f}_{B}(0))=(0.8000,1.1421)$, $(0.8000,1.2819)$, $(0.8000,1.3702)$, $(0.001,0.9993)$, $(0.001,0.9050)$, $(0.001,0.7071)$, $(0.001,0.6180)$, $(0.001,0.5176)$. The two types of soliton solutions can be seen in figures 5(d) and 6(d), respectively, corresponding to branch 1 and 2 in figure 3(d).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1. Oscillating multi-soliton solutions

Large amplitude oscillating solutions for strong nonlinearity such as those in figures 5(g)–(i) may be obtained analytically. Such solutions are bright solitons in the total density over a nonzero background. We note that analytical solutions in the massive case have recently been studied in detail by Khawaja [9]. We start by writing the two-spinor wavefunction in the form of a product of an envelope function $\eta (x)$ and the internal spinor degrees of freedom parameterized by the function $\varphi (x)$

$$\begin{eqnarray}&&{\Psi }_{\mathrm{zigzag}}(x)=\eta (x)\left(\begin{array}{c}\mathrm{cos}\varphi (x)\\ \mathrm{sin}\varphi (x)\end{array}\right),\end{eqnarray} \tag{ 57 }$$

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\to {x}_{0}$, which may included in the final solution. Substituting equation (57) into equations (8) and (9), multiplying by $\mathrm{cos}\varphi$ and $\mathrm{sin}\varphi$, respectively, then adding the resulting equations gives

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\varphi }{{\rm{d}}x}=-\displaystyle \frac{\mu }{{\rm{\hslash }}{c}_{l}}\left[1-{\eta }^{2}\;(U/\mu )({\mathrm{cos}}^{4}\varphi +{\mathrm{sin}}^{4}\varphi )\right].\end{eqnarray} \tag{ 58 }$$

To obtain a second equation we multiply equations (8) and (9) by $\mathrm{cos}\varphi$ and $\mathrm{sin}\varphi$, respectively, then subtract the resulting equations which yields

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}(\mathrm{ln}\eta )}{{\rm{d}}x}=\displaystyle \frac{U}{4{\rm{\hslash }}{c}_{l}}\;{\eta }^{2}\;\mathrm{sin}(4\varphi ).\end{eqnarray} \tag{ 59 }$$

Note that we have divided through by η to arrive at equations (58) and (59). Equations (58) and (59) can be combined by back substitution to get

$$\begin{eqnarray}&&\left[\displaystyle \frac{\mathrm{sin}(4\varphi )}{{\mathrm{cos}}^{4}\varphi +{\mathrm{sin}}^{4}\varphi }\right]{\left(\varphi +\displaystyle \frac{\mu x}{{\rm{\hslash }}{c}_{l}}\right)}^{\prime }=4\;(\mathrm{ln}\eta )\prime \;\end{eqnarray} \tag{ 60 }$$

where the prime notation indicates differentiation with respect to x. A formal expression for η is obtained from equation (60)

$$\begin{eqnarray}&&{\eta }^{4}=\mathrm{exp}\left[\int {\rm{d}}x\;\displaystyle \frac{\mathrm{sin}(4\varphi )}{{\mathrm{cos}}^{4}\varphi +{\mathrm{sin}}^{4}\varphi }{\left(\varphi +\displaystyle \frac{\mu x}{{\rm{\hslash }}{c}_{l}}\right)}^{\prime }\right].\end{eqnarray} \tag{ 61 }$$

To solve equation (61), we assume the linear form $\varphi (x)=\kappa x$, obtain an explicit form for $\eta (x)$, which we then substitute into equation (58) to determine the constant κ and obtain a relation for the chemical potential μ and the interaction U. Equation (61) becomes

$$\begin{eqnarray}&&{\eta }^{4}=\mathrm{exp}\left[(\kappa +\mu /{\rm{\hslash }}{c}_{l})\int {\rm{d}}x\;\displaystyle \frac{\mathrm{sin}(4\kappa x)}{{\mathrm{cos}}^{4}\kappa x+{\mathrm{sin}}^{4}\kappa x}\right],\end{eqnarray} \tag{ 62 }$$

which, upon integration, yields the result

$$\begin{eqnarray}&&\eta (x)=A{\left({\mathrm{cos}}^{4}\kappa x+{\mathrm{sin}}^{4}\kappa x\right)}^{-(1+\mu /\kappa {\rm{\hslash }}{c}_{l})/4},\end{eqnarray} \tag{ 63 }$$

where A is the integration constant. Substituting this result into equation (58) and using the linear assumption, gives the expression

$$\begin{eqnarray}&&\kappa =-\displaystyle \frac{\mu }{{\rm{\hslash }}{c}_{l}}\left[1-\displaystyle \frac{U}{\mu }{\left({\mathrm{cos}}^{4}\kappa x+{\mathrm{sin}}^{4}\kappa x\right)}^{1-(1+\mu /\kappa {\rm{\hslash }}{c}_{l})/2}\right].\end{eqnarray} \tag{ 64 }$$

Since κ is constant, it must be that the exponent of the spatial functions is identically zero. Equation (64) then gives the two conditions

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}-\displaystyle \frac{\mu }{2\kappa {\rm{\hslash }}{c}_{l}}=0,\end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray}&&-\displaystyle \frac{\mu }{{\rm{\hslash }}{c}_{l}}\left(1-\displaystyle \frac{U}{\mu }\right)=\kappa ,\end{eqnarray} \tag{ 66 }$$

which may be solved to give

$$\begin{eqnarray}&&\kappa =\displaystyle \frac{\mu }{{\rm{\hslash }}{c}_{l}},\end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray}&&\mu =2\;U.\end{eqnarray} \tag{ 68 }$$

The corresponding solution is then

$$\begin{eqnarray}&&\Psi (x)=A{\left({\mathrm{cos}}^{4}\kappa x+{\mathrm{sin}}^{4}\kappa x\right)}^{-1/2}\left(\begin{array}{c}\mathrm{cos}\;\kappa x\\ \mathrm{sin}\;\kappa x\end{array}\right),\end{eqnarray} \tag{ 69 }$$

where $\kappa =2U/({\rm{\hslash }}{c}_{l})$. The spinor components in equation (69) are plotted in figure 7(a) and the corresponding density in figure 7(b). Although equation (69) was obtained for the NLDE with real coefficients, a direct transformation to get the associated solution for the case of complex coefficients is obtained by taking ${\psi }_{A}\to {\rm{i}}{\psi }_{B}$, ${\psi }_{B}\to {\psi }_{A}$ to get

$$\begin{eqnarray}&&\Psi (x)=A{\left({\mathrm{cos}}^{4}\kappa x+{\mathrm{sin}}^{4}\kappa x\right)}^{-1/2}\left(\begin{array}{c}{\rm{i}}\;\mathrm{sin}\;\kappa x\\ \mathrm{cos}\;\kappa x\end{array}\right).\end{eqnarray} \tag{ 70 }$$

Zoom In Zoom Out Reset image size

A distinguishing feature of the solitons in figure 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 figure 7(b) are bright solitons on a nonzero background.

6.2. 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 equation (47). A preliminary step requires that we solve equation (47) for f_B and then back substitute into equation (44). This gives us the first order nonlinear equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{f}_{A}}{{\rm{d}}x}=\mp \displaystyle \frac{U}{{\rm{\hslash }}{c}_{l}}\left\{\left[\displaystyle \frac{\mu }{U}\pm \sqrt{2\displaystyle \frac{{T}^{11}}{U}+2{\left(\displaystyle \frac{\mu }{U}\right)}^{2}-{\left({f}_{A}^{2}-\displaystyle \frac{\mu }{U}\right)}^{2}}\right]\right.\end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray}&&\left.\times \left[2\displaystyle \frac{{T}^{11}}{U}+2{\left(\displaystyle \frac{\mu }{U}\right)}^{2}-{\left({f}_{A}^{2}-\displaystyle \frac{\mu }{U}\right)}^{2}\right]\right\}.\end{eqnarray} \tag{ 72 }$$

Equation (72) is separable in the variables f_A and x. In the special case where $\mu =0$, the problem simplifies and yields the integral

$$\begin{eqnarray}&&\int {\rm{d}}{f}_{A}{\left(2\displaystyle \frac{{T}^{11}}{U}-{f}_{A}^{4}\right)}^{-3/4}=\pm \displaystyle \frac{U}{{\rm{\hslash }}{c}_{l}}x+D,\end{eqnarray} \tag{ 73 }$$

where D is an integration constant. For $| {f}_{A}| \lt {\left(2{T}^{11}/U\right)}^{1/4}$ and ${T}^{11}\gt 0$, equation (73) shows that ${f}_{A}^{\prime }$ will remain positive for all x; thus a monotonically increasing and bounded solution exists between $-{\left(2{T}^{11}/U\right)}^{1/4}\lt {f}_{A}\lt +{\left(2{T}^{11}/U\right)}^{1/4}$, for $-\infty \lt x\lt +\infty$. The integral in equation (73) is given in terms of the hypergeometric function so that

$$\begin{eqnarray}&&\displaystyle \frac{U}{{\rm{\hslash }}{c}_{l}}x=\displaystyle \frac{{f}_{A}}{{\left({T}^{11}/U\right)}^{3/4}}{\ }_{2}{F}_{1}\left(\displaystyle \frac{1}{4},\displaystyle \frac{3}{4};\displaystyle \frac{5}{4};\displaystyle \frac{{{Uf}}_{A}^{4}}{{T}^{11}}\right)+D.\end{eqnarray} \tag{ 74 }$$

The second spinor component f_B is obtained by back substitution into the invariance relation where we find ${f}_{B}={\left(2{T}^{11}/U-{f}_{A}^{4}\right)}^{1/4}$. This type of solution exhibits a form of self-confinement similar to that studied in [10], where f_B remains localized near the region where the slope of f_A is steepest and has a shape similar to a $\mathrm{tanh}$ function.

For general values of μ, it is helpful to cast equation (72) in a more enlightening form by an appropriate parametric transformation. First, we make the substitution ${f}_{A}^{2}=a+b\;\mathrm{cos}\theta$, with $a\equiv \mu /U$ and ${b}^{2}\equiv 2{T}^{11}/U+2{(\mu /U)}^{2}$, which transforms equation (72) to

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\theta }{{\rm{d}}x}=\pm \;\displaystyle \frac{U}{{\rm{\hslash }}{c}_{l}}\sqrt{a+b\;\mathrm{cos}\theta }\ \sqrt{a\pm b\;\mathrm{sin}\theta }.\end{eqnarray} \tag{ 75 }$$

Note that in order to keep f_A real we must enforce the constraint $a\gt b$. Cast in terms of the new variable θ, we have the symmetric result $b\;\mathrm{cos}\theta ={f}_{A}^{2}-a$ and $b\;\mathrm{sin}\theta =a-{f}_{B}^{2}$. Equation (75) has an interesting graphical representation which we have depicted in figure 8 where we show how the dark and bright solitons emerge when the problem is cast in terms of the angular parameter θ.

Zoom In Zoom Out Reset image size

The dark and bright solitons may be found in the limit that $a\to 1$, $b\to 1$, and with the choice ${T}^{11}={\mu }^{2}/2U$ as in section 5. In this case ${\rm{d}}\theta /{\rm{d}}x=0$ at $\theta =\pi /2$ and $\theta =\pi$ for the negative sign under the radical in equation (75), and ${\rm{d}}\theta /{\rm{d}}x=0$ at $\theta =-\pi /2$ and $\theta =\pi$ for the positive sign under the radical. In particular, for the negative sign under the radical, the dark soliton and bright soliton interpolate between $\theta =\pi /2$ and $\theta =\pi$ by following anti-clockwise and clockwise paths, respectively, i.e., for the positive and negative outer signs in equation (75).

Equation (75) may be integrated exactly by separation of variables whereby one obtains

$$\begin{eqnarray}&&\displaystyle \frac{U}{{\rm{\hslash }}{c}_{l}}x+D=\pm \displaystyle \frac{\left(\mathrm{cos}\theta -\mathrm{sin}\theta +1\right)\mathrm{ln}\left[\left(\mathrm{cos}\theta /2\right){\left(\mathrm{cos}\theta /2-\mathrm{sin}\theta /2\right)}^{-1}\right]}{\sqrt{1+\mathrm{cos}\theta }\;\sqrt{1-\mathrm{sin}\theta }},\end{eqnarray} \tag{ 76 }$$

where D is the integration constant which adds a spatial translation to the soliton core. Equation (76) cannot be directly inverted to find explicit forms for ${f}_{A}(x)$ and ${f}_{B}(x)$ (via $\theta (x)$), but we may study it graphically by plotting ${Ux}/{\rm{\hslash }}{c}_{l}$ as a function of θ. In figure 9 we have plotted equation (76) where we have shown singular points occurring at repeating intervals: $...,\pi /2,\pi ,5\pi /2,3\pi ,...$, alternating between the two types of soliton solutions.

Zoom In Zoom Out Reset image size

Alternatively, equation (75) may be solved numerically to obtain an inversion of the solution plotted in figure 9, i.e., $\theta (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.

Zoom In Zoom Out Reset image size

The invariant form of equation (53) says that asymptotically flat soliton solutions must have the same topological invariant quantity ${\delta }_{\pm \infty }$ given by equation (54). We verify that this is the case by explicitly computing ${\delta }_{\pm \infty }$ using equation (75) and expressing the integral in equation (54) in terms of the variable θ. The resulting indefinite integral for ${\delta }_{\pm \infty }$ is

$$\begin{eqnarray}{\delta }_{\pm \infty } & = & {\displaystyle \int }_{{\theta }_{i}}^{{\theta }_{f}}{\rm{d}}\theta \sqrt{1+\mathrm{cos}\theta }\;\sqrt{1-\mathrm{sin}\theta }(\mathrm{cos}\theta +\mathrm{sin}\theta )\\ & = & \displaystyle \frac{\sqrt{1+\mathrm{cos}\theta }\;\sqrt{1-\mathrm{sin}\theta }[\mathrm{sin}\theta +(\mathrm{sin}\theta -1)\mathrm{cos}\theta ]}{1+\mathrm{cos}\theta -\mathrm{sin}\theta }{| }_{{\theta }_{i}}^{{\theta }_{f}}.\end{eqnarray} \tag{ 77 }$$

When evaluated at the appropriate limits for either of the two types of solitons, the integral in equation (77) yields ${\delta }_{\pm \infty }=0$. For the dark soliton ${\theta }_{i}=\pi /2$ and ${\theta }_{f}=\pi$, whereas ${\theta }_{i}=\pi /2$ and ${\theta }_{f}=-\pi$ for the bright soliton. We point out that the dark and bright solitons are distinct and cannot be deformed into each other. In the case of the former, the crossing point where ${f}_{A}={f}_{B}$ lies below the fixed point $(\sqrt{\mu /U},\sqrt{\mu /U})$, whereas for the latter the crossing point lies above the fixed point. Any continuous deformation in the NLDE solution space which moves the crossing point towards the fixed point must force f_A and f_B to flatten out everywhere pushing the crossing point out to $x\to \infty$. 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 figure 8 illustrates this point. For example, we are free to parametrize our problem more generally in terms of some variable τ: $b\to b(\tau )$ and $\theta \to \theta (\tau )$ in equation (75). And, since the singular point $(b=0,\theta )$ is a fixed point of equations (43) and (44), any continuous deformation of the functions $b(\tau )$ and $\theta (\tau )$ 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 figures 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 [62–64]. 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.

6.3. 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 $\mu /U$, interpolating through the soliton boundary and connecting the two regimes shown in figure 3. The series expansion uses the same ansatz as that in section 6.1, i.e., equation (57), whereby we find that $\eta (x)$ will take the form of a power series in the quantity ${\mathrm{cos}}^{4}\varphi +{\mathrm{sin}}^{4}\varphi$. The function $\varphi (x)$ can then be obtained by integrating $\eta (x)$ term by term. This offers a convenient method since powers of ${\mathrm{cos}}^{4}\varphi +{\mathrm{sin}}^{4}\varphi$ are exactly integrable.

We begin by substituting the form $\Psi (x)=\eta (x)(\mathrm{cos}\varphi (x),\mathrm{sin}\varphi (x))$ into the invariance relation equation (47). Upon this substitution, equation (47) becomes

$$\begin{eqnarray}&&{\eta }^{2}=\displaystyle \frac{1}{{\mathrm{cos}}^{4}\varphi +{\mathrm{sin}}^{4}\varphi }\left(\displaystyle \frac{\mu }{U}\right)\left[1\pm \sqrt{1+2\left(\displaystyle \frac{{T}^{11}}{U}\right){\left(\displaystyle \frac{U}{\mu }\right)}^{2}({\mathrm{cos}}^{4}\varphi +{\mathrm{sin}}^{4}\varphi )}\right],\end{eqnarray} \tag{ 78 }$$

and the associated equation for $\varphi (x)$ is equation (58) which takes the form

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\varphi }{{\rm{d}}x}=\pm \displaystyle \frac{\mu }{{\rm{\hslash }}{c}_{l}}\sqrt{1+2\left(\displaystyle \frac{{T}^{11}}{U}\right){\left(\displaystyle \frac{U}{\mu }\right)}^{2}({\mathrm{cos}}^{4}\varphi +{\mathrm{sin}}^{4}\varphi )}.\end{eqnarray} \tag{ 79 }$$

Because of the bounds $1/2\leqslant {\mathrm{cos}}^{4}\varphi +{\mathrm{sin}}^{4}\varphi \leqslant 1$, there are two limits for which equations (78) and (79) simplify, defined by ${| (2{T}^{11}/U)(U/\sqrt{2}\mu )| }^{2}\gg 1$ and ${| (2{T}^{11}/U)(U/\mu )| }^{2}\ll 1$. For ${| (2{T}^{11}/U)(U/\sqrt{2}\mu )| }^{2}\gg 1$, equations (78) and (79) may be expanded in an asymptotic series, whereby one finds

$$\begin{eqnarray}\eta {(x)}^{2} & = & {| 2{T}^{11}/U| }^{1/2}{\left({\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x)\right)}^{-1/2}\\ & & \times \underset{n=0}{\overset{\infty }{\displaystyle \sum }}\left(\begin{array}{c}1/2\\ n\end{array}\right){\left[(2{T}^{11}/U){(U/\mu )}^{2}\left({\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x)\right)\right]}^{-n},\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}x}\varphi (x) & = & \pm \displaystyle \frac{\mu }{{\rm{\hslash }}{c}_{l}}{\left({\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x)\right)}^{1/2}\\ & & \times \underset{n=0}{\overset{\infty }{\displaystyle \sum }}\left(\begin{array}{c}1/2\\ n\end{array}\right){\left[(2{T}^{11}/U){(U/\mu )}^{2}\left({\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x)\right)\right]}^{-n},\end{eqnarray} \tag{ 81 }$$

where the expansion coefficients are the generalized binomial coefficients

$$\begin{eqnarray}&&\left(\begin{array}{c}1/2\\ n\end{array}\right) & \equiv & \displaystyle \frac{1/2(1/2-1)(1/2-2)\cdots }{n!\;(1/2-n)(1/2\;n-1)(1/2\;n-2)\cdots }\\ & = & \displaystyle \frac{1/2(1/2-1)\cdots (1/2\;n+1)}{n!}.\end{eqnarray} \tag{ 82 }$$

The terms in the expansion for $\varphi (x)$ can be integrated and expressed in terms of elliptic integrals. The leading order terms in the asymptotic expansion are

$$\begin{eqnarray}&&\eta (x)\approx {| 2{T}^{11}/U| }^{1/4}{\left({\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x)\right)}^{-1/4}+{\mathcal{O}}{\left[(2{T}^{11}/U){(U/\sqrt{2}\mu )}^{2}\right]}^{-1},\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&\varphi (x)\approx \pm {{\rm{E}}}^{-1}\left(\displaystyle \frac{2{T}^{11}{Ux}}{{\rm{\hslash }}{c}_{l}}| \displaystyle \frac{1}{2}\right)+{\mathcal{O}}{\left[(2{T}^{11}/U){(U/\sqrt{2}\mu )}^{2}\right]}^{-1},\phantom{\rule{9em}{0ex}}\end{eqnarray} \tag{ 84 }$$

where the solution for $\varphi (x)$ is expressed in terms of the inverse elliptic integral of the first kind [65]. Conversely, at the other extreme limit ${| (2{T}^{11}/U)(U/\mu )| }^{2}\ll 1$ the following expansion is valid

$$\begin{eqnarray}\eta {(x)}^{2} & = & {\left({\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x)\right)}^{-1}\left(\displaystyle \frac{\mu }{U}\right)\\ & & \times \left[1\pm \underset{n=0}{\overset{\infty }{\displaystyle \sum }}\left(\begin{array}{c}1/2\\ n\end{array}\right){\left[(2{T}^{11}/U){(U/\mu )}^{2}\left({\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x)\right)\right]}^{n}\right],\end{eqnarray} \tag{ 85 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}x}\varphi (x)=\pm \displaystyle \frac{\mu }{{\rm{\hslash }}{c}_{l}}\underset{n=0}{\overset{\infty }{\displaystyle \sum }}\left(\begin{array}{c}1/2\\ n\end{array}\right){\left[(2{T}^{11}/U){(U/\mu )}^{2}\left({\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x)\right)\right]}^{n}.\qquad \quad \end{eqnarray} \tag{ 86 }$$

To leading order, only one nonzero solution exists, coming from the positive sign inside the square brackets in equation (78):

$$\begin{eqnarray}&&\eta (x)\approx \displaystyle \frac{{(U/2\mu )}^{1/2}}{{\left[{\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x)\right]}^{1/2}}+{\mathcal{O}}\left[(2{T}^{11}/U){(U/\mu )}^{2}\right],\end{eqnarray} \tag{ 87 }$$

$$\begin{eqnarray}&&\varphi (x)\approx \pm \displaystyle \frac{\mu x}{{\rm{\hslash }}{c}_{l}}+{\mathcal{O}}\left[(2{T}^{11}/U){(U/\mu )}^{2}\right].\end{eqnarray} \tag{ 88 }$$

We point out that equation (88) is just the series of bright solitons first obtained in equation (69) which we now identify as the large $\mu /U$ limit solution in figure 5(i). A second solution can be obtained associated with the negative sign inside the brackets in equation (78), coming from the next to leading order term in equation (85):

$$\begin{eqnarray}&&\eta (x)\approx \sqrt{\displaystyle \frac{| {T}^{11}| }{2\mu }}+{\mathcal{O}}{\left[(2{T}^{11}/U){(U/\mu )}^{2}\right]}^{2},\end{eqnarray} \tag{ 89 }$$

$$\begin{eqnarray}&&\varphi (x)\approx \pm \displaystyle \frac{\mu x}{{\rm{\hslash }}{c}_{l}}+{\mathcal{O}}\left[(2{T}^{11}/U){(U/\mu )}^{2}\right].\end{eqnarray} \tag{ 90 }$$

Equations (89) and (90) describe plane wave solutions identified with the oscillating solutions in panel figure 6(i). To make the connection to the single dark and bright solitons in figures 5(d) and 6(d), we set ${T}^{11}={\mu }^{2}/2U$ then Taylor expand equations (78) and (79) in powers of $| 1-(U/\mu )|$ for $\mu /U\approx 1$. Such an expansion is valid as long as $| 1-(U/\mu )| \lt 1$. Solving equations (78) and (79) to leading order in $| 1-(U/\mu )|$ gives

$$\begin{eqnarray}{\eta }^{2}(x) & \approx & \displaystyle \frac{1}{{\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x)}\\ & & \times \left[1\pm \sqrt{1-({\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x))}\right]+{\mathcal{O}}[| 1-(U/\mu )| ],\end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}x}\varphi (x)\approx \pm \displaystyle \frac{\mu }{{\rm{\hslash }}{c}_{l}}\sqrt{1-({\mathrm{cos}}^{4}\varphi (x)+{\mathrm{sin}}^{4}\varphi (x))}+{\mathcal{O}}[| 1-(U/\mu )| ].\qquad \quad \end{eqnarray} \tag{ 92 }$$

A solution for $\varphi (x)$ can be obtained by the method of separation of variables, whereby we find

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{sin}(2\varphi )\;\mathrm{ln}\left[\mathrm{tan}(2\varphi )\right]}{\sqrt{1-\mathrm{cos}(4\varphi )}}=\pm \displaystyle \frac{\mu }{{\rm{\hslash }}{c}_{l}}x+D,\end{eqnarray} \tag{ 93 }$$

where D is the integration constant which spatially shifts the solution, as found also in equations (73)–(76). One may invert equation (93) graphically and substitute this result into equation (92), leading to the dark and bright solitons, respectively, for the positive and negative signs in equation (92).

So far we have focused our attention on solving the NLDE using analytical methods. In section 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 [66].

The most direct approach is to express equations (8) and (9) in terms of the dimensionless spatial variable $\chi \equiv (U/{\rm{\hslash }}{c}_{l})x$. The spinor functions ${f}_{A}(\chi )$ and ${f}_{B}(\chi )$ are then expanded in a power series around $\chi =0$

$$\begin{eqnarray}&&{f}_{A}(\chi )=\underset{j=0}{\overset{\infty }{\displaystyle \sum }}{a}_{j}{\chi }^{j},\ \ \ \ {f}_{B}(\chi )=\underset{j=0}{\overset{\infty }{\displaystyle \sum }}{b}_{j}{\chi }^{j},\end{eqnarray} \tag{ 94 }$$

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 equation (94) into equations (8) and (9) gives us the behavior of the solution at the origin

$$\begin{eqnarray}&&{f}_{A}^{\prime }(0)\sim \displaystyle \frac{1}{{\rm{\hslash }}{c}_{l}}\left[\mu -{{Uf}}_{B}{(0)}^{2}\right]{f}_{B}(0),\quad \end{eqnarray} \tag{ 95 }$$

$$\begin{eqnarray}&&{f}_{B}^{\prime }(0)\sim -\displaystyle \frac{1}{{\rm{\hslash }}{c}_{l}}\left[\mu -{{Uf}}_{A}{(0)}^{2}\right]{f}_{A}(0).\end{eqnarray} \tag{ 96 }$$

By examination we find that for $j\geqslant 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

$$\begin{eqnarray}&&{b}_{1}-{a}_{0}^{3}=-{a}_{0},\end{eqnarray} \tag{ 97 }$$

$$\begin{eqnarray}&&2{b}_{2}-3{a}_{1}{a}_{0}^{2}=-{a}_{1},\end{eqnarray} \tag{ 98 }$$

$$\begin{eqnarray}&&3{b}_{3}-3{a}_{2}{a}_{0}^{2}-3{a}_{0}{a}_{1}^{2}=-{a}_{2},\end{eqnarray} \tag{ 99 }$$

$$\begin{eqnarray}&&4{b}_{4}-3{a}_{3}{a}_{0}^{2}-3{a}_{2}{a}_{1}^{2}-{a}_{1}^{3}=-{a}_{3},\end{eqnarray} \tag{ 100 }$$

resulting from equation (8), and

$$\begin{eqnarray}&&-{a}_{1}-{b}_{0}^{3}=-{b}_{0},\end{eqnarray} \tag{ 101 }$$

$$\begin{eqnarray}&&-2{a}_{2}-3{b}_{1}{b}_{0}^{2}=-{b}_{1},\end{eqnarray} \tag{ 102 }$$

$$\begin{eqnarray}&&-3{a}_{3}-3{b}_{2}{b}_{0}^{2}-3{b}_{0}{b}_{1}^{2}=-{b}_{2},\end{eqnarray} \tag{ 103 }$$

$$\begin{eqnarray}&&-4{a}_{4}-3{b}_{3}{b}_{0}^{2}-3{b}_{2}{b}_{1}^{2}-{b}_{1}^{3}=-{b}_{3},\phantom{\rule{8em}{0ex}}\end{eqnarray} \tag{ 104 }$$

which come from equation (9). To obtain soliton solutions using the shooting method, we first fix either a₀ or b₀, and then vary the other until we obtain convergence to the desired precision. Although it seems that both a₀ and b₀ are free parameters, fixing one to a different value before shooting results in a spatially translated final solution, as we verified 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}_{0}^{\mathrm{soliton}}$ to the desired precision gives the soliton configuration. Figure 11 shows the shooting process as we tune a₀ through the dark soliton which is shown in figure 11(c) for which ${a}_{0}^{\mathrm{soliton}}=0.999\;292\;150\;145\pm {10}^{-8}$. Note that converging to a value ${a}_{0}^{\mathrm{soliton}}\lt 1$ ensures that the total density remains in the dark soliton regime. A similar process converges on a second value ${a}_{0}^{\mathrm{soliton}}\gt 1$ associated with the bright soliton.

Zoom In Zoom Out Reset image size

As the precision in the value of ${a}_{0}^{\mathrm{soliton}}$ is increased, oscillations are pushed out to larger values of x. In figures 11(a) and (b) we show the solution for values of ${a}_{0}\gt {a}_{0}^{\mathrm{soliton}}$. At larger values of a₀, the nonlinearity becomes dominant and we start to pick up some of the excited soliton states which can be seen in figures 11(a) and (b). For values ${a}_{0}\lt {a}_{0}^{\mathrm{soliton}}$, the effect of the interaction is reduced, figures 11(d) and (e), until finally we see the free particle sine and cosine forms appearing in figure 11(f). The particular values of the constant a₀ in figure 11 are: (a) ${a}_{0}=1.1\pm {10}^{-19}$, (b) ${a}_{0}=0.999\;2922\pm {10}^{-13}$, (c) ${a}_{0}=0.999\;292\;150\;145\pm {10}^{-8}$, (d) ${a}_{0}=0.99\pm {10}^{-18}$, (e) ${a}_{0}=0.9\pm {10}^{-19}$, (f) ${a}_{0}=0.5\pm {10}^{-19}$. 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 [67].

In this article we have presented a variety of methods for solving the armchair NLDE and mapped out the soliton landscape. A discrete symmetry allows for two types of NLDEs associated with the quasi-1D reduction to the armchair geometry in the plane of a honeycomb optical lattice. This discrete symmetry is expressed in terms of Pauli matrices acting on the order parameter for bosons propagating along the length of the armchair pattern in the lattice.

In particular, we have found dark and bright solitons for nonzero chemical potential and a self-confined soliton for the case $\mu =0$. Both dark and bright solitons show high-contrast density fringes and are therefore clearly observable in experiments. The spinor component functions for the dark soliton are approximated by the forms ${\psi }_{A}(x)\sim (1/2)[1+\mathrm{tanh}(x)]$ and ${\psi }_{B}(x)\sim (1/2)[1-\mathrm{tanh}(x)]$, while the bright soliton components resemble the forms ${\psi }_{A}(x)\sim (1/2)[1+\mathrm{tanh}(x)+\mathrm{sech}(x)]$ and ${\psi }_{B}(x)\sim (1/2)[1-\mathrm{tanh}(x)+\mathrm{sech}(x)]$, where ${\psi }_{A}(x)$ and ${\psi }_{B}(x)$ are the spatial parts of the upper and lower two-spinor components. The crossing point where ${\psi }_{A}={\psi }_{B}$ lies below the fixed point ${\psi }_{A}={\psi }_{B}=\sqrt{\mu /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. Thus, the two solution regimes, large and small local densities, are topologically distinct. In the case where $\mu =0$, a third soliton solution exists which resembles a flattened $\mathrm{sech}$ function in ${\psi }_{B}$ localized at the center of a $\mathrm{tanh}$ form in ${\psi }_{A}$.

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 $\Delta (x)$, but not the total density. In the case of the Fermi gas, a dark soliton may be observed in the gap function $\Delta (x)$, the order parameter that encodes pairing of fermion particles and holes in the Bardeen–Cooper–Schrieffer (BCS) to BEC crossover. There the notch depth in $\Delta (x)$ varies with the translational speed of the soliton. The speed is characterized by an overall phase $\phi (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 [68, 69]. 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, suggesting possible future research topics. 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 [68]? All such questions relate to the issue of Cooper pairing between NLDE quasi-particles, 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.

This material is based in part upon work supported by the National Science Foundation under grant number PHY-1306638, AFOSR grant FA9550-14-1-0287, the Alexander von Humboldt foundation, and the Heidelberg Center for Quantum Dynamics. We acknowledge useful discussions with Ken O'Hara at Pennsylvania State University.

To compute the error of our oscillating numerical solutions, we first calculate the average difference between values of ${\psi }_{A}$ and ${\psi }_{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 $\varepsilon (x)$ is

$$\begin{eqnarray}&&{\varepsilon }_{A(B)}(x)\equiv \displaystyle \frac{\mathrm{diff}(x)}{\mathrm{avg}(x)}=2\left[\displaystyle \frac{{\psi }_{A(B)}(x+L)-{\psi }_{A(B)}(x)}{{\psi }_{A(B)}(x+L)+{\psi }_{A(B)}(x)}\right],\end{eqnarray} \tag{ A.1 }$$

where L is the periodicity for the particular solution and the error is computed for both two-spinor component functions ${\psi }_{A}$ and ${\psi }_{B}$. In figure A1 , we have plotted ${{log}}_{10}| \varepsilon (x)|$, (a)–(d), and $\varepsilon (x)$, (e)–(h), for the solution depicted in figure 5(a) using grid sizes ${\rm{N}}={10}^{4},\;{10}^{5},\;{10}^{6},\;{10}^{7}$. The solution in figure 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 equation (56) to obtain ${f}_{B}(0)$.

Zoom In Zoom Out Reset image size