{'fraction_non_alphanumeric': 0.046564164812192325, 'fraction_numerical': 0.015307842740575898, 'mean_word_length': 4.280220758411963, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solitary Waves in Discrete Media with Four Wave Mixing

In this paper, we examine in detail the principal branches of solutions that arise in vector discrete models with nonlinear inter-component coupling and four wave mixing. The relevant four branches of solutions consist of two single mode branches (transverse electric and transverse magnetic) and two mixed mode branches, involving both components (linearly polarized and elliptically polarized). These solutions are obtained explicitly and their stability is analyzed completely in the anti-continuum limit (where the nodes of the lattice are uncoupled), illustrating the supercritical pitchfork nature of the bifurcations that give rise to the latter two, respectively, from the former two. Then the branches are continued for finite coupling constructing a full two-parameter numerical bifurcation diagram of their existence. Relevant stability ranges and instability regimes are highlighted and, whenever unstable, the solutions are dynamically evolved through direct computations to monitor the development of the corresponding instabilities. Direct connections to the earlier experimental work of Meier et al. [Phys. Rev. Lett. {\bf 91}, 143907 (2003)] that motivated the present work are given.


I. INTRODUCTION
Recently, nonlinear Hamiltonian lattice dynamical systems with a large number of degrees of freedom have become a focal point for a variety of application areas [1]. Several diverse physical contexts in which such models (typically discrete in space and continuous in the evolution variable) arise are (i) spatial dynamics of optical beams in coupled waveguide arrays arising in nonlinear optics [2], (ii) temporal evolution of Bose-Einstein condensates (BECs) in optical lattices in soft-condensed matter physics [3] and (iii) DNA double strand in biophysics [4].
A prototypical model, applicable to different degrees of approximation, in all of the above contexts is the so-called discrete nonlinear Schrödinger (DNLS) equation. This model was first proposed in [5] and implemented for the first time experimentally in Refs. [6,7]. A systematic presentation of the experimental results in optical systems is given in a recent work [8]; for BEC-related experiments, see [3,9]. In this model, self-localized excitations (discrete solitons) are possible as a result of the interplay between the Kerr nonlinearity and discrete linear coupling. Many properties of optical discrete spatial solitons have been systematically explored in theory and experiment, including generalizations to diffraction management [10,11], diffraction-managed solitons [12], and soliton transport and gating [13,14,15,16].
On the other hand, a topic that has received considerably less attention has been the study of vector analogs of the DNLS equation. A number of studies have addressed the existence and stability of diverse families of solitary waves/localized states in the DNLS equation (see [17,18,19]). These issues have also been studied for both cubic and quadratic nonlinearities [18,20] in one and two dimensions [17,21,22]. However, the first experimental realization for such a system occurred rather recently [23] (see also the very recent paper [24]). This work sheds new light on the theoretical investigations of this topic, in that the relevant model for the corresponding AlGaAs waveguide array experiment included four wave mixing (FWM) terms that were rarely included in previous theoretical studies. Subsequent studies (see [25,26]) addressed a number of specific properties of the model system put forward in [23]. Such properties related to switching, instability-induced amplification, modulational instability as well as an analysis of the energy barrier between stable discrete solitary waves centered on a lattice site versus two lattice sites.
The focus of the present study is to extend the analysis of the vector NLS lattice system incorporating FWM which is relevant to the experiments of [23]. More specifically, we provide a systematic analysis of the existence and stability requirements of the principal solitary wave modes of the system. In so doing, we have found that the concept of the anti-continuum limit (in which the coupling between waveguides is absent) provides a powerful tool for the study of the existence and stability concerning these solutions. Subsequently, continuations in the coupling parameter are used to construct a full two-parameter bifurcation diagram that gives the complete picture of stability for both one-component as well as "mixed" mode solutions (involving both components). We note that all of the results obtained via our dimensionless model can be directly placed in the context of the experimental results obtained in Refs. [23,24].
The remainder of the paper is organized in the following way. In section II, we present the theoretical model of interest as well as proposing a setup that allows for analysis of the existence and stability of the model's solutions. In section III, we fully analyze the relevant solutions and their stability in the anti-continuum limit (i.e., zero coupling between adjacent waveguides). In section IV, we proceed numerically to determine solutions to the model for nonzero coupling and observe how the bifurcation points and stability features are modified for this case. The instability regions obtained in section IV are then monitored through direct numerical experiments in section V. Finally, in section VI, we summarize our findings and present our conclusions.

II. MODEL AND SETUP
Following [23], we will use the dimensionless model iȧ n = −a n − ǫ (a n+1 + a n−1 ) − |a n | 2 + A|b n | 2 a n − Bb 2 n a ⋆ n (1) In these equations, a n and b n are the appropriately normalized, slowly-varying, complex field envelopes for the transverse electric (TE) and transverse magnetic (TM) polarized waves respectively. The constants A and B are respectively associated with the cross-phase modulation (XPM) and four-wave mixing (FWM) and were evaluated in [23] to be approximately equal to A ≃ 1 and B ≃ 1/2. Wherever possible, the results obtained in this analysis will be given for general A and B, even though the above values have been used in the specifics of our numerical computations. The overdot in Eqs. (1)-(2) denotes differentiation with respect to the evolution variable (which is z in this case) and the ⋆ denotes the complex conjugate. For the properties of the waveguide array and incident light used in the experiment (see [23] for details), the dimensionless power is connected with its dimensional analog P d (measured in Watts) according to P d ≃ 56.4P . Furthermore, the dimensionless coupling ǫ used in [23] was ǫ ≃ 0.921. This is crucial since it serves to highlight differences between the results obtained in the computational analysis of [23] and the generalization of the corresponding results presented herein. More generally, we take ǫ to be a free parameter in our analysis. We note that ǫ is related to the dimensional coupling constant, κ, given in [23], via ǫ = 2.741κ. Here, κ is measured in mm −1 .
We look for stationary solutions in the form: a n =ã n e iqz and b n =b n e iqz and subsequently dropping the tildes, we obtain the corresponding stationary equations: (q − 1)a n − ǫ (a n+1 + a n−1 ) − |a n | 2 + A|b n | 2 a n − Bb 2 n a ⋆ n = 0 (4) In this context, q is an additional free parameter, namely the propagation constant (q is also used as a free parameter for the computations presented in [23]). Our parameter q is related to the corresponding dimensional propagation constant, δ, given in [23], by δ = 0.3648q. Here, δ is measured in mm −1 .
It is worth pointing out that for Fig. 1 in Ref. [23], the range of parameters used in their computations appear to correspond to the case of q > 3. This is verified by our analysis below. Again, this is highlighted to account for the disparities of some of the conclusions obtained in [23] with the ones their generalizations presented in the following sections.
To perform linear stability computations, provided that an exact solution (a 0 n , b 0 n ) of the stationary Eqs. (4)-(5) has been obtained (either analytically in the anti-continuum limit, or numerically for ǫ = 0), we use the perturbation ansatz: In this ansatz, δ is a (small) formal expansion parameter and ω is the eigenfrequency of the perturbation eigenvector (c n , d ⋆ n , f n , g ⋆ n ) T (where the superscript denotes the transpose). Then, upon lengthy but straightforward calculation, one obtains the linearized matrix eigenvalue problem to O(δ). These dynamical equations have the form: The N × N (where N is the size of the lattice) blocks of the linearization matrix are given by: In the above, the short-hand notation (∆ 2 + 2)z n = z n+1 + z n−1 , invoking the concept of the discrete Laplacian ∆ 2 . With this linearized matrix eigenvalue problem, we can establish a criterion for the existence and stability of solutions to this system by examination of our system at the anti-continuum limit (i.e., when ǫ = 0).

A. Existence
We begin by examining eqs. (4) − (5) for the case ǫ = 0. Setting a n = r n e iθn and b n = s n e iφn in this case, one finds: From the above, we infer that for solutions to exist in this limit it is necessary for 2(θ n − φ n ) to be an integer multiple of π, hence, we obtain that with k ∈ Z.
The simplest possible solutions are the ones that involve only one of the two branches and were hence termed TE and TM modes respectively in [23]. The TE solution of Eqs. (24)-(25) has the form (in the present limit) and exists only for q > 1. On the other hand, the TM mode features and is only present for q > −1. Also, there are two possible mixed mode solutions permitted by the quantization condition of eq. (26), allowing the exponential of Eqs. (24)-(25) to be +1 or −1 respectively. The first one (e 2i(θn−φn) = 1) was characterized as a linearly polarized (LP) branch in [23], involving in-phase contributions from both the TE and TM components. In this case, the linear system of Eqs. (24)- (25) has the general solution: If (A + B) 2 > 1 (which is the case in the experimentally relevant setting), this branch only exists for (A + B)(q + 1) > (q − 1) and (A + B)(q − 1) > (q + 1) [the sign of the two above inequalities needs to be reversed for existence conditions in the case of (A + B) 2 < 1]. Among the two conditions, in the present setting, the second one is the most "stringent" for the case A = 2B = 1, which yields the constraint q ≥ 5 (while the first condition requires for the same parameters q ≥ −5).
The second possibility for a mixed mode (stemming from e 2i(θn−φn) = −1 in Eqs. (24)-(25)) yields the, so-called, elliptically polarized mode (EP) [23], involving a π/2 phase shift between the TE and TM components. The relevant amplitudes in this case are: The condition for this branch to exist if . In this case, the first condition is more constraining than the second, imposing for our special case of interest herein (A = 2B = 1) q ≥ 3 (while the second condition only requires q ≥ −3).

B. Stability
We can now turn to the study of the stability of the corresponding branches via the anti-continuum limit (ǫ = 0). The major advantage of this limit is that all the blocks of the linearization matrix L now become diagonal. Furthermore, the structures that we consider, per the branches presented above, involve the excitation of a single site (all other sites are inert) with amplitudes corresponding to the appropriate TE, TM, LP or EP modes. The interesting byproduct of this formulation is that all inert sites end up yielding purely diagonal entries in the stability matrix equal to (q − 1) in the first (2N-1)×(2N-1) blocks and (q + 1) in the second set of (2N-1)×(2N-1) blocks. Hence the matrix will have 2N − 1 eigenfrequencies equal to (q − 1) and equally many eigenfrequencies of (q + 1). The remaining 4 × 4 matrix will correspond to the single excited site of the lattice and its eigenvalues will yield the eigenfrequencies pertaining to the corresponding branch of solutions. Hence, by solving this considerably simpler eigenvalue problem, we can infer the full stability picture for ǫ = 0 from the anti-continuum limit. This, we believe, is a major advantage in using this approach.
We note that due to the fact that the power is conserved for our system (see previous section), two of the relevant eigenvalues of this 4 × 4 problem should be identically 0. This is confirmed in the computations presented here. Now, taking the case for the "TE" mode in eqs.(8) − (23), we determine the remaining pair of eigenvalues in this case: Examining our model for the experimental case, A = 1, B = 1/2 and q > 0, we find that this eigenfrequency is real for q < 5, while it is imaginary for q > 5 (hence, implying the presence of an instability, since the corresponding eigenvalue λ = iω becomes real and exponential growth occurs ∝ e λz ). Indeed, we can show that for the TE mode, stability is ensured for 1 < q < 5. For the TM mode, the corresponding expression can be obtained as: and for the case of interest, stability is ensured for −1 < q < 3, while instability ensues (due to a real eigenvalue) for q > 3. For the LP mode, the expression for the corresponding eigenfrequency is too cumbersome to provide in analytic form (even though such a form has been obtained) for general A and B. Hence, we will restrict ourselves to the case of A = 2B = 1, where the relevant eigenfrequency reads: implying stability for q ≥ 5 (which coincides with the branch's existence region). Finally, for the EP mode, the expression similarly reads: implying, in principle, stability for q ≥ 3 (again, coinciding with the solutions domain of existence); see, however, the discussion below for this branch. From the above observations, one can put together the complete bifurcation picture of the four branches and their corresponding stability in the anti-continuum limit (again, we present this here for concreteness for A = 2B = 1). More specifically, the TE branch exists for q ≥ 1 and is stable for 1 ≤ q ≤ 5. For q > 5, the branch is destabilized as a new branch emerges, namely the LP branch, through a pitchfork bifurcation; notice that the TM component of this branch, per Eq. (30) is exactly zero at q = 5, hence it directly bifurcates from the TE branch. The two branches of this super-critical pitchfork correspond to the two signs of s n in Eq. (30). The bifurcating branch "inherits" the stability of the TE branch for all larger values of q, while the latter branch remains unstable thereafter.
In a similar development, the TM branch exists for q ≥ −1 and is stable in the interval −1 ≤ q ≤ 3. However, at q = 3, a new branch (in fact, a pair of branches) bifurcates acquiring non-zero r n , beyond the bifurcation point, as per Eq. (31). This is accompanied by the destabilization of the TM branch (due to a real eigenvalue) and the apparent stability of the ensuing EP branch for all values of q > 3. Fig. 1 shows how the two supercritical pitchfork bifurcations are generated as a function of the coordinates (r n , s n , q). We can now turn on the coupling and see how the existence and stability of each branch is affected by varying q and ǫ.

IV. FINITE COUPLING
One of the key features that will help us in our analysis of the finite coupling case consists of determining where the bands of the continuous spectrum of the linearized problem occurs. This is obtained analytically using the plane wave ansatz a n ∼ e i(kn−ωz) and b n ∼ e i(kn−ωz) . Plugging our ansatz into the coupled linearized problem (essentially, using eqs. (1) − (2) and ignoring the nonlinear terms), the resulting dispersion relations are Hence the continuous spectrum (in the infinite lattice case) will consist of the frequency intervals [q − 1 − 2ǫ, q − 1 + 2ǫ] and [q + 1 − 2ǫ, q + 1 + 2ǫ]. As may be expected, for ǫ = 0, these intervals degenerate to two single points (ω = q − 1 and ω = q + 1), in consonance with the discussion of the previous section. But perhaps more importantly, for q − 1 < 2ǫ (equivalently for ǫ > (q − 1)/2), the continuous spectrum branch will be crossing the origin, leading to the collision of the eigenvalues with their mirror symmetric opposites, that will in turn lead to instabilities. Hence, we need only consider couplings in the interval ǫ ∈ [0, (q − 1)/2]. We now numerically examine each of the main four branches, using one-parameter continuation in the abovementioned interval of ǫ's for different values of q and thus constructing a two-parameter bifurcation diagram for each branch. Notice that the solutions are numerically obtained for each pair of (q, ǫ) by means of a fixed point iteration. Once convergent to within a pre-set tolerance (10 −8 typically for the computations presented herein), the procedure is followed by numerical linear stability analysis, obtaining the eigenfrequencies for the linearization matrix L.

A. TE branch
The continuation of the TE branch is detailed in Fig. 2. It turns out that especially in the case of this branch, solutions cannot be obtained for ǫ > (q−1)/2, i.e., the branch terminates at that point with its amplitude going to zero at this critical point. Within its region of existence, the branch has a domain of stability and one of instability. The point of separation between the two in the anti-continuum limit, studied previously, was the critical point of q = 5. For ǫ = 0, the separatrix curve is shown in Fig. 2 and can be well approximated by the curve ǫ c T E ≈ (4 √ 2/5) √ q − 5. Hence for q ≤ 5, the solution is stable for all values of ǫ in its range of existence (0 < ǫ < (q − 1)/2), while for q ≥ 5, the solution is only stable for ǫ c T E < ǫ < (q − 1)/2 and unstable (due to a real eigenvalue pair) for 0 < ǫ < ǫ c T E .

B. TM branch
The TM branch is considerably more complicated, structurally, than its TE counterpart. Firstly, it does not disappear beyond the critical ǫ = (q − 1)/2, however, it does become unstable as predicted previously, hence we will, once again, restrict ourselves to this parameter range. Furthermore, similarly to the TE branch case, there is an ǫ c T M below which the branch is always unstable, whereas for ǫ > ǫ c T M , the branch may be stable. At the anti-continuum limit, the critical point for the instability is q = 3, as discussed previously; for q > 3, the critical point is obtained numerically in Fig. 3. It can be well approximated (close to q = 3) by ǫ c T M ≈ (9/10) √ q − 3. However, within the range of potential stability (0 ≤ ǫ ≤ (q − 1)/2 for q ≤ 3, and ǫ c T M ≤ ǫ ≤ (q − 1)/2 for q ≥ 3), we observe an additional large region of instability in the two parameter bifurcation diagram of Fig. 3, due to a complex quartet of eigenvalues. This instability appears to "emanate" from the point with (q, ǫ) = (2.2, 0) in the anti-continuum limit, and to linearly expand its range as ǫ increases. Hence, we will try to understand it at the level of ǫ = 0. What happens at the value of q = 2.2 in the ǫ = 0 limit is that the point spectrum eigenfrequency of Eq. (34) "collides" with the continuous spectrum point of concentration, corresponding to ω = q − 1. However, the eigenvector of this eigenvalue [34] has a Krein signature (see e.g. [27,28,29] for relevant definitions of this signature) opposite to that of the continuous band at ω = q − 1. Hence, according to Krein theory [30], the resulting collision leads to the formation of a quartet of eigenvalues emerging in the complex plane and, in turn, implying the instability of the TM configuration. As ǫ is varied from 0, the continuous spectrum band grows linearly in ǫ, hence the corresponding interval of q's, where this instability is present (due to the collision with the opposite Krein sign eigenvalue) also grows at the same rate. Along the same vein, it is worth pointing out that the line of this instability threshold and that of ǫ = (q − 1)/2 are parallel.

C. LP branch
The mixed LP branch, involving in-phase contributions of the TE and TM has been continued for various values of q > 5, as a function of ǫ. Recall that this branch emerges through a supercritical pitchfork as q is varied for for fixed ǫ (see Fig. 2). In an exactly analogous way, if ǫ is varied for a fixed q > 5, the branch appears to terminate at ǫ = ǫ c T E through a subcritical pitchfork. Hence, the LP branch only arises for 0 < ǫ < ǫ c T E and q > 5 and it is stable throughout its interval of existence. For ǫ > ǫ c T E , this branch "degenerates" back into the stable (for such ǫ's) TE mode. This is clearly illustrated in Fig. 4 that shows the continuation of the LP branch as a function of ǫ for a fixed q = 5.5. The top panels show the norms of the two components illustrating that the second component is absent for ǫ > ǫ c T E = 0.81 in this case. The middle and bottom left panels show a case for ǫ = 0.5 < ǫ c T E where a genuine (and stable) TE solution exists, while the corresponding right panels are for ǫ = 1 > ǫ c T E , past the critical point where the LP solution degenerates into a TE mode.

D. EP branch
The case of the EP branch, analyzed in Fig. 5, is somewhat analogous to that of the LP branch. For fixed q close to (and larger than) 3 and varying ǫ, the branch exists and is stable for 0 < ǫ < ǫ c T M . For larger values of ǫ, the EP mode degenerates into a TM mode, as is shown in the top left panel of Fig. 5 for q = 3.5 (cf. with the analogous for the LP case of Fig. 4).  Fig. 3, but now for the TM branch. The top left set of panels shows the continuation in ǫ for q = 2.5, along with the special cases of ǫ = 0.1 (stable) and ǫ = 0.5 (unstable due to quartet of eigenvalues) in the left and right middle and bottom panels. The top right set of panels shows the q = 3.5 (unstable) case with the middle/bottom panels corresponding to ǫ = 0.5 (unstable due to real eigenvalue pair)/ǫ = 0.65 (stable) respectively.
However, for q > 3.62, this phenomenology appears to change and an expanding (for increasing q) interval of oscillatory instability within the range of existence of the EP branch (0 < ǫ < ǫ c T M ) seems to arise. To rationalize this feature, we again turn to the anti-continuum limit. Just as in the case of the (TM) branch from which it originates, the EP branch has a point spectrum eigenfrequency given by Eq. (36) which has a negative Krein signature and upon collision with the continuous spectrum band of eigenfrequencies will give rise to an instability. Setting the frequency of Eq. (36) equal to q − 1, we obtain that this collision occurs at q = 9. For lower values of q, this "collision" will occur for a finite (nonzero) interval of values of ǫ -see e.g. the top right panels of Fig. 5 for q = 4.5-, which, in turn, yields the region of oscillatory instability of the EP mode, clearly separated from its region of stability by a dashed line in the bottom panel of Fig. 5.
In conclusion, the EP solutions exist in the regime where the TM mode is unstable (cf. bottom panels of Figs. 3 and 5), illustrated by the curve originating from the point (q, ǫ) = (3, 0) in the anti-continuum limit. However, the EP solution, as is shown in the top panels of the figure, may (top left) or may not (top right) necessarily be stable in its region of existence due to the appearance of an oscillatory instability for sufficiently large q.

V. DYNAMICAL EVOLUTION OF UNSTABLE BRANCHES
We now examine a series of different instability scenarios pertaining to the branches of solutions discussed above. In so doing, we solve Eqs. (1)-(2) using a fourth-order Runge-Kutta algorithm. To ensure the validity of our numerical results, we have monitored the conservation of the power of Eq. (3) and of the Hamiltonian (energy) of the system: The power is conserved to numerical precision, while the Hamiltonian is conserved to 1 part in 10 8 . In all cases, 150 spatial sites are used in our computations.
Our first example, shown in Fig. 6, illustrates the case of q = 5.5 and ǫ = 0.5 for the TE mode. This is a case representative of the region of parameter space (cf. Fig. 2), where the TE mode is unstable due to a real eigenvalue pair (while the LP mode is stable). The solution actually coincides with the steady TE branch until approximately time 10. The configuration then deviates from the TE branch and subsequently develops a breathing oscillation around a stable LP state.
The next three cases are concerned with the different scenarios for the instability of the TM mode. Fig. 7, with q = 3.5 and ǫ = 0.5, is representative of the region of parameter space where the TM mode is unstable due to a real eigenvalue pair, while the EP mode is stable. In this case, As in the previous case, the solution coincides with the steady TM mode up until approximately time 10. It is then observed that for short times the configuration deviates from the TM branch and oscillates around the (stable) EP configuration. However, for longer times, the system appears to wander away from the latter configuration and into a breathing state, where the two components have roughly equal power and exchange a small fraction of it (almost) periodically.
The second instability scenario of the TM branch, shown in Fig. 8, involves the case of an eigenvalue quartet as e.g. for q = 2.5 and ǫ = 0.5 (cf. Fig. 3). In this case, we observe the configuration departing from the unstable TM steady state and resulting into a breathing solution involving both components. In this and the remaining cases, the departure from the steady state occurs later than in the first two cases.
Next, we examine the evolution of the instability for the TM mode in the presence, also, of continuous spectrum instabilities for q = 2.5 and ǫ = 1 (cf. Fig. 3). In this case a different phenomenology occurs. The mode is completely destroyed in favor of small amplitude excitations.
Finally, we examine the evolution of the instability for an elliptically polarized mixed mode case within its regime of oscillatory instability (cf. Fig. 5), for ǫ = 0.5 and q = 5. We observe the configuration departing the steady state and giving rise to a persistent breathing state involving both components.

VI. CONCLUSIONS
In this paper, we have examined in detail the principal modes arising in the discrete vector lattices of AlGaAs waveguides, motivated by the recent experimental investigations of [23,24]. We have constructed the full twoparameter bifurcation diagrams of all four principal branches of solutions. These include the transverse electric and transverse magnetic modes (occupying respectively each of the two components of the vector lattice), but also the linearly polarized (two components in or out of phase) and elliptically polarized modes (two components out of phase by odd multiples of π/2). The latter two modes, emerge from the former two in pairs through a pitchfork bifurcation that has been completely elucidated through the study of the system in the anti-continuum limit. In fact, the anticontinuum limit has provided with a very effective tool to analytically highlight existence and stability properties of the different branches as well as potential instabilities. We have continued this monoparametric picture (as a In the latter case, the solution has already degenerated into a TM mode. Top right panels: q=4.5. The middle and bottom panels correspond to the stable case of ǫ = 0.1 (left panels) and the unstable case of ǫ = 1 (right panels). The bottom panel shows the two-parameter diagram for this mode. The region for q > 3 and 0 < ǫ < ǫ c T M , where this modes exists is separated by the dashed line into a stable and an unstable regime. For comparison the extension of the regions of stability/instability of the TM mode from Fig. 3 are also included. We note that the choice of star symbols (instead of circles as e.g. in Figs. 2, 4) in the spatial profile of an is used to indicate that the imaginary part of an is shown. It is noteworthy that the process of obtaining such a multiparametric bifurcation picture has allowed us to reveal a considerably deeper understanding of the system's parameter space, in comparison to earlier publications (such as e.g. the theoretical conclusions of the works of [23,24]). In particular, for instance, the limited parameter space examined in the latter setting (where a monoparametric continuation was used for ǫ ≃ 0.921 and apparently varying q only for q > 3) led the authors to infer that the TM mode is always unstable and similarly that the bifurcating EP mode is also unstable. Here, we have illustrated explicitly regions of stability of each of these modes in the parameter space of the system (which definitely exist for lower values of the coupling).
Finally, while the recent investigations of the various principal branches have offered a rather substantial understanding of their properties, there are multiple directions of future interest in this topic that also appear within experimental reach. Such examples consist e.g. of the study of two-dimensional waveguide arrays and of discrete solitons [31] and vortices [32] in them. Another such example is the study of more elaborate, multi-site constructions, such as twisted modes [33] even in one spatial dimension. While the above directions of examination of stationary states are of interest in their own right, an additional dimension of relevance revealed by our direct instability simulations of section V pertains to the systematic study of the existence and stability of breathing excitations due to their principal role in the system's dynamics. Such investigations are currently in progress and will be reported in future publications.