Solitary waves in a two-dimensional nonlinear Dirac equation: from discrete to continuum

TL;DR

This paper investigates solitary waves in a two-dimensional nonlinear Dirac equation, analyzing their existence, stability, and evolution from discrete to continuum models in waveguide arrays.

Contribution

It demonstrates the continuation of continuum Dirac solitons across discretization levels and explores their stability and dynamics in discrete and continuum regimes.

Findings

Continuum Dirac solitons persist for all discretization parameters.

Solutions with 1- or 2-sites at the anti-continuum limit can be continued to large couplings.

Some solutions are found to be spectrally unstable and exhibit specific dynamical behaviors.

Abstract

In the present work, we explore a nonlinear Dirac equation motivated as the continuum limit of a binary waveguide array model. We approach the problem both from a near-continuum perspective as well as from a highly discrete one. Starting from the former, we see that the continuum Dirac solitons can be continued for all values of the discretization (coupling) parameter, down to the uncoupled (so-called anti-continuum) limit where they result in a 9-site configuration. We also consider configurations with 1- or 2-sites at the anti-continuum limit and continue them to large couplings, finding that they also persist. For all the obtained solutions, we examine not only the existence, but also the spectral stability through a linearization analysis and finally consider prototypical examples of the dynamics for a selected number of cases for which the solutions are found to be unstable.

Equations

$$i{\partial }_t{\psi }_1$$

$$i{\partial }_t{\psi }_2$$

$$i{\partial }_t{\psi }_1$$

$$i{\partial }_t{\psi }_2$$

$$\phi(\vec{r})=\left[\begin{array}[]{c}u(r)\mathrm{e}^{iS\theta}\\ i\,v(r)\mathrm{e}^{i(S+1)\theta}\end{array}\right],$$

$$\left({\partial }_r+\frac{S+1}{r}\right)v+(m-\omega -g{u}^2)u$$

$$-\left({\partial }_r-\frac{S}{r}\right)u-(m+\omega +g{v}^2)v$$

$$\psi(\vec{r},t)=\left[\begin{array}[]{c}\left\{u(r)+\delta\left[a_{1}(r)\mathrm{e}^{iq\theta}\mathrm{e}^{\lambda t}+b^{*}_{1}(r)\mathrm{e}^{-iq\theta}\mathrm{e}^{\lambda^{*}t}\right]\right\}\mathrm{e}^{iS\theta}\\[8.61108pt] \!i\left\{v(r)+\delta\left[a_{2}(r)\mathrm{e}^{iq\theta}\mathrm{e}^{\lambda t}+b_{2}^{*}(r)\mathrm{e}^{-iq\theta}\mathrm{e}^{\lambda^{*}t}\right]\right\}\mathrm{e}^{i(S+1)\theta}\end{array}\!\right]\mathrm{e}^{-i\omega t},$$

$$\mathcal{M}_{q}=i\left(\begin{array}[]{cc}L_{1}&L_{2}\\ \\ -L_{2}^{*}&-L_{1}^{*}\end{array}\right)-\frac{iq}{r}\left(\begin{array}[]{cc}\sigma_{1}&0\\ \\ 0&\sigma_{1}\end{array}\right),$$

$$L_{1}=\left(\begin{array}[]{cc}\omega-m+2gu^{2}&-\left(\partial_{r}+\frac{S+1}{r}\right)\\ \\ \left(\partial_{r}-\frac{S}{r}\right)&\omega+m+2gv^{2}\end{array}\right)\ ,\quad L_{2}=\left(\begin{array}[]{cc}gu^{2}&0\\ \\ 0&gv^{2}\end{array}\right)\ ,\quad\sigma_{1}=\left(\begin{array}[]{cc}0&I\\ \\ I&0\end{array}\right)\ .$$

$$\sigma \left(\mathcal{M}\right)=\bigcup _{q\in \mathbb{Z}}\sigma \left({\mathcal{M}}_q\right).$$

$$i{\partial }_t{V}_{n,m}$$

$$i{\partial }_t{U}_{n,m}$$

$$Q=\sum _n\sum _m{\rho }_{n,m},{\rho }_{n,m}=|U_{n,m}{|}^2+|V_{n,m}{|}^2,$$

$$H=-\frac{1}{2}\sum _n\left[C{V}_{n,m}^{\ast }(i{\nabla }_x{U}_{n,m}+{\nabla }_y{U}_{n,m})+C{U}_{n,m}^{\ast }(i{\nabla }_x{V}_{n,m}-{\nabla }_y{V}_{n,m})-\frac{g}{k+1}(|U_{n,m}{|}^4+|V_{n,m}{|}^4)+m(|U_{n,m}{|}^2-|V_{n,m}{|}^2)\right].$$

$$\mathrm{i}{\dot{U}}_{n,m}=\frac{\delta H}{\delta U_{n,m}^{\ast }},\mathrm{i}{\dot{V}}_{n,m}=\frac{\delta H}{\delta V_{n,m}^{\ast }}.$$

$$(\omega -m+g|v_{n,m}{|}^2)v_{n,m}+C[i{\nabla }_x{u}_{n,m}+{\nabla }_y{u}_{n,m}]$$

$$(\omega +m+g|u_{n,m}{|}^2)u_{n,m}+C[i{\nabla }_x{v}_{n,m}-{\nabla }_y{v}_{n,m}]$$

$$U_{n,m}(t)=e^{-i\omega t}\left[u_{n,m}+\delta (a_{n,m}{e}^{\lambda t}+c_{n,m}^{\ast }{e}^{{\lambda }^{\ast }t})\right],V_{n,m}(t)=e^{-i\omega t}\left[v_{n,m}+\delta (b_{n,m}{e}^{\lambda t}+d_{n,m}^{\ast }{e}^{{\lambda }^{\ast }t})\right]$$

$$\mathcal{M}=i\left(\begin{array}[]{cc}L_{1}&L_{2}\\ \\ -L_{2}^{*}&-L_{1}^{*}\end{array}\right)$$

$$L_{1}=\left(\begin{array}[]{cc}\omega-m+2g|u|^{2}&C(i\nabla_{x}+\nabla_{y})\\ \\ C(i\nabla_{x}-\nabla_{y})&\omega+m+2g|v|^{2}\end{array}\right)\ ,\quad L_{2}=\left(\begin{array}[]{cc}u^{2}&0\\ \\ 0&v^{2}\end{array}\right)\ .$$

$$\epsilon (t)=\frac{|Q(t)-Q(0)|}{Q(0)},$$

Full text

I Introduction

Optical waveguide arrays dnc ; moti constitute one of the settings that have led to numerous recent experimental and theoretical developments as regards the analysis of wave phenomena in Hamiltonian lattices. Both in this and in the related context of photorefractive crystals, features such as discrete diffraction yaron and its management yaron1 , Talbot revivals christo2 , $\mathcal{PT}$-symmetry and its breaking kip , as well as discrete solitons yaron ; yaron2 and vortices neshev ; fleischer were not only theoretically predicted but also experimentally observed. Variants of the theme of optical waveguide arrays have involved multi-component models bearing multiple polarizations dncm ; rudy , waveguides featuring quadratic (so-called $\chi^{2}$) nonlinearities stege ; rudy2 , and the examination of dark-solitonic states shandarov ; hadi . Another theme where extensive related studies have been conducted is the atomic physics realm of Bose-Einstein condensates (BECs) in optical lattices ober ; konotop1 .

Recently, research on binary waveguide arrays has been gaining momentum Tran ; akyl ; jesuskip ; Aceves ; yannan . Part of the reason for this is that under suitable limiting conditions, this system can lead to Dirac-like nonlinear equations that are of increasing prominence and wide relevance in diverse physical contexts. These include, among others, bosonic evolution in honeycomb lattices Carr1 ; Carr2 and a growing class of atomically thin two-dimensional (2D) Dirac materials diracmat such as graphene, silicene, germanene and transition metal dichalcogenides tmdc . They also arise when studying light propagation in honeycomb photorefractive lattices (the so-called photonic graphene) peleg ; ablowitz3 ; ablowitz4 . These Dirac settings have been argued to present fundamental differences, e.g., with respect to their more well known, non-relativistic limits of the nonlinear Schrödinger equation, such as the absence of collapse for an extended interval of frequencies in two spatial dimensions; see, e.g., the recent work of PRL and the discussion therein.

Our aim in the present work is to revisit the context of binary waveguide arrays, motivated by the realistic models studied in Tran ; akyl for two-dimensional geometries. In an earlier work, we considered the phenomenology of a 1D discrete nonlinear Dirac equation JPA . However, this was done for a discretization of the less straightforwardly applicable (from a physical perspective) Soler model. Here, we turn our attention to the more realistic setting of the waveguide arrays with onsite nonlinearity, aiming to explore existence, spectral stability and dynamics of nonlinear modes. We do this in a complementary way between the (highly) discrete and the continuum limits. On the one hand, we explore the configurations at strong coupling (the continuum soliton) and subsequently reduce the coupling going towards the discrete limit. Here, we eventually find that a configuration bearing nine sites turns out to be the limiting highly discrete analogue of the continuum solitary wave. On the other hand, we also start our search for model solutions from the highly discrete limit of vanishing coupling (the so called anti-continuum limit of macaub ) with one- or two-site configurations and continue them into the strong coupling regime. Utilizing a spectral stability analysis, we identify the regimes of lattice coupling as well as of solution frequency for which the relevant waveforms are dynamically stable. When instability is identified, some prototypical examples of the configuration’s unstable evolution are given.

Our exposition will be structured as follows: in Sec. II we provide an overview of the theoretical properties of both the discrete and the continuous model. In Sec. III, we numerically explore the existence, stability and dynamical features of the models, while in Sec. IV we summarize our findings and present our conclusions as well as some interesting directions for future work.

II Theoretical Setup

Following the setup of akyl ; Tran , the two-dimensional (continuum) Dirac model of relevance to the binary waveguide problem is of the form:

[TABLE]

where $(\psi_{1},\psi_{2})$ denotes the spinor field (mode amplitude), while $m$ represents the mass associated with the propagation mismatch between the two different types of waveguides. The cubic nonlinearity stems from the Kerr effect and breaks the Lorentz symmetry (contrary, e.g., to what is the case with the Soler model; cf. JPA ; PRL ). We consider this model both at the discrete and at the continuum level.

II.1 Continuous model

The analysis of the continuum model of Eq. (II) can be performed in radial coordinates, following a procedure similar to the one proposed in PRL . This results in the form:

[TABLE]

The form of this equation suggests that we look for stationary solutions as $\psi(\vec{r},t)=\exp(-i\omega t)\phi(\vec{r})$ with

[TABLE]

where $u(r)$ and $v(r)$ are real-valued. The value $S\in\mathbb{Z}$ can be cast as the vorticity of the first spinor component (the second component in this formulation has vorticity $S+1$). The equation fulfilled by stationary profiles then only depends on $r$, casting the problem into a 1D one:

[TABLE]

with $r>0$. As we discuss in the next section, stationary solutions are sought by numerical means.

In order to capture the linear stability of the stationary solutions, we introduce the following ansatz into (II.1):

[TABLE]

and subsequently solve the ensuing [to O$(\delta)$] eigenvalue problem: $\lambda(a_{1},a_{2},b_{1},b_{2})^{T}=\mathcal{M}_{q}(a_{1},a_{2},b_{1},b_{2})^{T}$ with $\mathcal{M}_{q}$ being

[TABLE]

and

[TABLE]

The key observation which facilitates a computation of the spectrum is that the explicit form of Eq. (6) for $\mathcal{M}_{q}$ contains $r$ and $\partial_{r}$, but not $\theta$. This allows us to compute the full 2D stability spectrum $\mathcal{M}$ as the union of spectra of the one-dimensional spectral problems:

[TABLE]

In what follows, for concreteness we will set $m=g=1$ (as this choice can be made by renormalizing the time and the wave function). Now, the only free parameter that will be considered in the continuum limit is the frequency $\omega$.

II.2 Discrete model

The discrete version of Eq. (II) will be based on a discretization similar to that considered, e.g., in akyl (cf. the discussion around Eq. (2) therein). From a numerical approximation perspective, this is tantamount to a centered difference discretization of the first derivative in Eq. (II) and leads to:

[TABLE]

with $U_{n,m}$ and $V_{n,m}$ ($-N/2+1\leq(n,m)\leq N/2$) being the components of the spinor (amplitude modes) $\Psi_{n,m}\equiv(U_{n,m},V_{n,m})$, and $\nabla_{x}\Psi_{n,m}\equiv(\Psi_{n+1,m}-\Psi_{n-1,m})$, $\nabla_{y}\Psi_{n,m}\equiv(\Psi_{n,m+1}-\Psi_{n,m-1})$ being the $x$ and $y$ components of the discrete gradient. The connection to the corresponding continuum limit can be assigned by selecting $C=1/(2h)$ with $h$ being the lattice spacing (discretization parameter).

The dynamical system of Eq. (II.2) presents a number of conserved quantities, such as the charge (squared $\ell^{2}$ norm):

[TABLE]

with $\rho_{n,m}$ being the charge density, and the Hamiltonian:

[TABLE]

Equations (II.2) can be derived from the Hamiltonian (11) by means of Hamilton’s equations:

[TABLE]

Our main focus hereafter will be on stationary solutions and their stability as well as their dynamics. Such solutions can be found by using $U_{n,m}(t)=\exp(-i\omega t)u_{n,m}$, $V_{n,m}(t)=\exp(-i\omega t)v_{n,m}$, when they possess frequency $\omega$ and satisfy the coupled algebraic equations:

[TABLE]

Once stationary solutions of the algebraic system of Eqs. (II.2) are calculated (by, e.g., fixed point methods as discussed below), their linear stability is considered by means of a linearized stability analysis. More specifically, considering small perturbations [of order ${\rm O}(\delta)$, with $0<\delta\ll 1$] of the stationary solutions, we substitute the ansatz

[TABLE]

into Eqs. (II.2), and then solve the ensuing [to O$(\delta)$] eigenvalue problem: $\lambda(a_{n,m},b_{n,m},c_{n,m},d_{n,m})^{T}=\mathcal{M}(a_{n,m},b_{n,m},c_{n,m},d_{n,m})^{T}$ with $\mathcal{M}$ being

[TABLE]

and

[TABLE]

The potential existence of an eigenvalue with non-vanishing real part suggests the existence of a dynamical instability. If all the eigenvalues are imaginary, then the solution is spectrally (neutrally) stable.

As in the continuum limit, we will set $m=g=1$ and vary $\omega$ as well as $C$ as our relevant parameters in order to characterize the behavior of the solution and the variation of its stability properties.

III Numerical results

III.1 Continuous model

We start our exposition by showing the numerical results regarding fundamental solutions ($S=0$ solitary waves, for which we nevertheless note that their second component $\psi_{2}$ bears a vortex of charge $1$) and $S=1$ vortices in the continuous setting. Numerical analysis has been performed in a similar fashion as in Ref. PRL , using spectral methods for dealing with spatial derivatives. Figure 1 shows several examples of the profiles for $S=0$ and $S=1$ stationary solutions. To assess stability, Fig. 2 shows the dependence on the frequency of the real and imaginary parts of the eigenvalues for $S=0$ and $S=1$ solitary waves. We observe that, similar to the Soler model PRL , the continuum $S=0$ solitons are unstable below a critical frequency ($\omega=0.388$) for this model. However, it is interesting to observe that contrary to the Soler model, only $q=0$ instabilities are present for the $S=0$ case; these instabilities are of exponential nature and, consequently, can be predicted by the Vakhitov–Kolokolov criterion, as the curve representing charge versus frequency presents a maximum at the bifurcation point (see bottom panel of Figs. 2). In addition, solutions get more localized with decreasing frequency and tend to be localized at $r=0$ as $\omega\rightarrow 0$.

We can see from the stability analysis of Fig. 2 that for the fundamental solution, there are wide intervals of stability/instability for $S=0$, while the $S=1$ waveform is unstable for all values of $\omega$ (and coupling $C$); cf. with the similar results of PRL . This prompts us to turn to the dynamical evolution of the instability. The study of the dynamics of the unstable solutions for $S=0$ shows that they feature collapse, as is depicted in Fig. 3. As a numerical diagnostic for the accuracy of the simulation, we have monitored the relative norm error $\varepsilon$ defined as

[TABLE]

with $Q(t)$ being the soliton’s charge. We have observed that the norm is preserved within a factor $\sim 10^{-7}$, despite an obvious departure of the norm from its initial value as collapse occurs.

In movie1 one can find a movie with the soliton evolution. In every movie, the top panels correspond to the density of each spinor component and bottom ones to their phase.

III.2 Discrete model

We now turn to the existence, stability and dynamics of solitons in the discrete 2D nonlinear Dirac equation (NLDE) in the form of Eq. (II.2). The solutions are obtained by making use of fixed point methods based on the anti-continuous (AC) limit macaub , that allows to solve Eq. (II.2). The principal challenge in this case is to identify a suitable solution in the AC limit given that there are many solutions at $C=0$ that can be extended to the continuous $C\rightarrow\infty$ limit. In fact, all the solutions we have analyzed can be traced upon increasing the coupling strength $C$ towards the continuum limit.

Among all the solutions in the AC limit, remarkably we find that the one that leads to the $S=0$ solitary waves of the continuum limit – discussed in the previous section – is the 9-site soliton, i.e. $u_{0,0}=u_{0,\pm 1}=u_{\pm 1,0}=u_{\pm 1,\pm 1}=\sqrt{1-\omega}$. The $v$ field must be vanishing at $C=0$, as can be seen from Eq. (II.2). The left panel of Fig. 4 shows the profile of a typical such solution at finite coupling. Notice that in the left panels of Fig. 4, one can observe that, contrary to the soliton in the continuum, the imaginary part of $u_{n,m}$ is not null; however, $||\mathrm{Im}(u_{n,m})||_{\infty}$ tends asymptotically to 0 when $C\rightarrow\infty$. An additional advantage of the AC limit is that the decoupled nature of the lattice enables us to analytically calculate the spectrum at $C=0$. This consists of 9 pairs of eigenvalues at 0, $N^{2}-9$ pairs at $\lambda=\pm i(1+\omega)$ and $N^{2}$ pairs at $\lambda=\pm i(1-\omega)$. As in the 1D case, some of the eigenvalues become real when the coupling is switched on and the soliton is unstable for finite coupling. In particular, 7 eigenvalue pairs detach from zero whereas 2 additional pairs remain at zero for every coupling. Of the seven remaining pairs, two become imaginary and five pairs acquire a real part; for very low coupling $C\lesssim 0.01$, four among these real pairs experience Hamiltonian Hopf bifurcations yielding complex quartets. This scenario persists for $C<0.80$. Beyond this value, a rather complicated scenario arises, as can be seen in Fig. 5, where the dependence of the stability eigenvalues as a function of the coupling constant $C$ for $\omega=0.7$ is shown. One can observe that for large $C$ there are only two sources of oscillatory instabilities: (1) for $C\gtrsim 2.7$, one of the quartets that exists for small $C$ persists even for high $C$; (2) at $C\approx 8.8$, an additional Hamiltonian Hopf bifurcation takes place. Contrary to the Soler model JPA , the $2i\omega$ eigenvalue is not present either in the discrete or in the continuum limit.

It is intriguing to note that for the discretization considered at least one of the Hamiltonian Hopf bifurcations persists for large values of $C$. Such a feature has been previously encountered in the finite difference analysis of the Soler model in JPA . Another interesting feature that occurs below a certain critical frequency is the existence of two branches of solitons, one starting at $C=0$ and another one finishing at $C\rightarrow\infty$. These branches are connected by an intermediate one through two saddle-center bifurcations i.e., the dependence of $Q$ on $C$ has an S-shaped form. This type of bifurcation point arises for higher values of $C$ when $\omega$ decreases. Figure 6 shows the relevant dependence of the charge with respect to $C$ for $\omega=0.7$ and $\omega=0.6$. We can observe in the latter case the existence of this intermediate branch, as well as the associated fold points.

In a similar fashion to the 1D case JPA , 1-site solitons also exist and can be continued to the continuum limit. In the AC limit, they are given by $u_{0,0}=\sqrt{1-\omega}$, while $v_{0,0}=0$. The spectrum at $C=0$ comprises $N^{2}-1$ pairs at $\lambda=\pm i(1+\omega)$, $N^{2}$ pairs at $\lambda=\pm i(1-\omega)$ and a sole pair at $\lambda=0$. As the single pair at 0 must remain there because of the U(1) symmetry, no eigenvalue can depart from 0 and the soliton is stable, at least for small values of the coupling parameter $C$.

The middle panel of Fig. 4 shows the profile of a 1-site soliton for $\omega=0.7$ and $C=1$. Notice that such solitons, for every coupling, possess the following properties:

• •

$\mathop{\mathrm{Re}}(u_{n,m})=0$ if $n$ and $m$ are odd,

• •

$\mathop{\mathrm{Im}}(u_{n,m})=0$ if $n$ and $m$ are even,

• •

$\mathop{\mathrm{Re}}(v_{n,m})=0$ if $n$ is odd and $m$ is even,

• •

$\mathop{\mathrm{Im}}(v_{n,m})=0$ if $n$ is even and $m$ is odd,

resembling the properties of the soliton of Fig. 2 of Tran in the large $C$ limit. This also endows the solitons’ real and imaginary parts with a “staggered” structure with alternating rows missing, as seen in the middle panels of Fig. 4. Figure 7 shows the dependence of the complex eigenvalues on $C$. One can see that the soliton is stable (see the real part in the bottom left panel of the figure) for $C$ below $1.04$. At this critical point there is a bifurcation caused by a mode that destabilizes, after bifurcating from the linear modes band. Above this point, the soliton is exponentially unstable, becoming stable again at $C=1.631$. However, at $C=1.73$ it experiences a similar instability anew, while the structure finally stabilizes $\forall C\geq 3.0$ in this case of $\omega=0.7$.

Another interesting structure is the 2-site soliton, which, in the continuum limit is reminiscent of the soliton of Fig. 5 in Tran . In the AC limit, such a wave structure is given by $u_{0,0}=u_{0,1}=\sqrt{1-\omega}$, while once again the $v$ field is vanishing. The spectrum at $C=0$ is composed by $N^{2}-2$ pairs at $\lambda=\pm i(1+\omega)$, $N^{2}$ pairs at $\lambda=\pm i(1-\omega)$ and two pairs at $\lambda=0$. As in the 9-site case, the two pairs remain invariant at [math], and thus no eigenvalue departs from 0, making the structure spectrally stable for small $C$.

The right panel of Fig. 4 shows the profile of a 2-site soliton for $\omega=0.7$ and $C=1$. Here too, similarly to the 1-site soliton, we observe a staggered pattern with the following features:

• •

$\mathop{\mathrm{Re}}(u_{n,m})=0$ if $n$ is odd,

• •

$\mathop{\mathrm{Im}}(u_{n,m})=0$ if $n$ is even,

• •

$\mathop{\mathrm{Re}}(v_{n,m})=0$ if $n$ is even,

• •

$\mathop{\mathrm{Im}}(v_{n,m})=0$ if $n$ is odd.

Figure 7 shows the dependence of the linearization eigenvalues with $C$. One can see that the behavior is essentially the same as in the case of the 1-site soliton, with the bifurcations taking place at the same points as in the 1-site case. In fact, the spectrum is almost identical in both 1-site and 2-site case, as the panels of Fig. 7 show.

We have also considered the stability of the above mentioned configurations for fixed coupling near the continuum limit and variable frequency. Figure 8 shows the profile of such solutions for $\omega=0.5$ and $C=5$.

Figure 9 shows, for $C=2.5$ and $C=5$, the dependence of the complex eigenvalues with $\omega$ for the 9-site solitons. Notice that solitons only exist above a critical value of $\omega(C)$; this critical value decreases with $C$. For instance, solitons only exist for $\omega\geq 0.47$ ($\omega\geq 0.22$) if $C=2.5$ ($C=5$). This is caused by a bifurcation similar to that shown in Fig. 6. This bifurcation is not found (at least for the considered coupling $C=5$) in the 1-site and 2-site solitons (see Fig. 10). In the latter case, we observe that while generally dynamical instabilities may arise for small values of $\omega$ as well as in a narrow interval in the vicinity of $\omega=1$, a wide parametric interval of frequencies exists where the solitary waves are dynamically stable.

Finally, we have considered the dynamics of prototypical unstable solutions. To this aim, we have included the spectral plane of the corresponding solution and a link to the movie with the evolution.

We start with the typical evolution for the 9-site soliton at small coupling, where it possesses several coexisting instabilities (see Fig. 11a). The evolution for $C=0.5$ and $\omega=0.7$, shown in movie2 , leads to the destruction of the structure. In the case of 1-site solitons, we consider two cases corresponding to each exponential instability (Figs. 11b-c); the first one leads to the expansion i.e. dispersion of the solitary wave (movie3 , for $\omega=0.7$ and $C=1.45$), and the second one to slow soliton motion (movie4 , for $\omega=0.7$ and $C=1.9$). For the 2-site soliton, where the spectral dependence on $C$ is qualitatively similar to the 1-site soliton, the dynamics for the first instability is similar to the 1-site case (i.e. the soliton expands with time); however, for the second instability (see spectrum in Fig. 11d) we observe that the 2-site soliton breaks up (movie5 , for $\omega=0.7$ and $C=1.9$), repartitioning its mass into predominantly single site structures.

Another interesting regime for the dynamical observation of the solutions’ instability regards the setting close to the continuum limit. To this end, we have fixed $C=5$ and observed the 9-site soliton dynamics for three cases (see Figs. 11e-g), which are traced in movie6 ($\omega=0.7$), movie7 ($\omega=0.6$) and movie8 ($\omega=0.3$). In the first case, there is only an oscillatory instability that leads to the transformation of the soliton in a pair of precessing (and periodically recombining) solitary structures; in the second case, there is an additional oscillatory instability whose consequence is the eventual destruction of the solitary wave after it is initially converted into a 1-site soliton pair. Finally, in the third case, apart from the two oscillatory instabilities identified in this case, there is a dominant exponential instability which leads to the soliton’s expansion and subsequent pulsation. An example of the instability of the 1-site soliton as we approach the continuum limit (for large $C=5$ and $\omega=0.2$) is shown in movie9 . Here the soliton ends up breathing as a result of its exponential instability.

IV Conclusions & Future Challenges

In the present work, we have given a systematic account of some of the prototypical solitary states that a two-dimensional model of a Dirac type can bear as stationary solutions. Our work was motivated by models of binary waveguide lattices that have recently appeared in the literature Tran ; akyl ; thus, we have considered nonlinearities that are onsite in each component. We have explored the model from two complementary perspectives. We have identified the continuum limit solution and extended it all the way to the anti-continuum limit where somewhat surprisingly we have found it to correspond to a 9-site solution. On the other hand, we have constructed some of the simplest solutions of the anti-continuum limit, such as the 1- and 2-site ones and extended them over all couplings towards the continuum limit. In addition to the existence problem, we have provided a road map towards the corresponding stability properties. The 1- and 2-site solutions with their staggered structure appear to be rather robust and, in the exception of some finite intervals of instability, appear to feature stable dynamics. On the other hand, the 9-site solution contains considerably more directions of potential instability, yet most of these disappear in the large coupling regime. The unstable dynamics of the different waveforms were also considered showing examples of breathing, mobility, and fission as potential manifestations of the instability, depending on the solution of interest and its specific (frequency and coupling) parameters.

Finally, the current study suggests a number of future directions of interest. In the context of the discrete version of the nonlinear Schrödinger equation, a systematic perturbative analysis was developed from the anti-continuum limit that enabled a characterization of the stability features in the vicinity of this limit and the development of an understanding of the conditions under which structures near this limit might be stable dnlsbook . A similar theory seems to be within reach in the case of the Dirac model (see also yannan for a recent analysis in 1D binary waveguide arrays), but has not been developed as of yet. On the other hand, and although it is less relevant to the optical problem per se, extending Dirac-like models and associated consideration to three-dimensional settings would be a particularly challenging theme of work. Here, once again the continuum limit preliminary conclusions of PRL suggest possible existence of stable solutions, which may in principle be possible to continue between the continuum and anti-continuum limit and be spectrally stable in wide parametric intervals between these two limits. Studies along these directions are currently in progress and will be reported in future publications.

Acknowledgements. P.G.K. gratefully acknowledges the support of NSF-PHY-1602994, the Alexander von Humboldt Foundation, the Stavros Niarchos Foundation via the Greek Diaspora Fellowship Program, and the ERC under FP7, Marie Curie Actions, People, International Research Staff Exchange Scheme (IRSES-605096). J.C.-M. thanks the European Regional Development Funds program (EU-FEDER) and the MEIC (project MAT2016-79866-R) for financial support. This work was supported in part by the U.S. Department of Energy.