Two-Component Gap Solitons in Self-Defocusing Photonic Crystals Two-Component Gap Solitons in Self-Defocusing Photonic Crystals

Studies of solitons in spatially periodic (lattice) potentials have grown into a vast area of research, with profoundly important applications to nonlinear optics, plasmonics, and matter waves in quantum gases, as outlined in recent reviews [1]-[4]. In ultracold bosonic and fermionic gases, periodic potentials can be created, in the form of optical lattices (OLs), by coherent laser beams illuminating the gas in opposite directions [5]-[7]. Effective lattice potentials for optical waves are induced by photonic crystals (PhCs), which are built as permanent structures by means of various techniques [2, 8, 9], or as laser-induced virtual structures in photorefractive crystals [10]. Reconfigurable PhCs can be also based on liquid crystals [11]. Parallel to the progress in the experiments, the study of the interplay between the nonlinearity and periodic potentials has been an incentive for the rapid developments of theoretical methods [12, 13]. Both the experimental and theoretical results reveal that solitons can be created in lattice potentials, if they do not exist in the uniform space [this is the case of gap solitons (GSs) supported by the self-defocusing nonlinearity, see original works [14]-[17] and reviews [7, 18]], and solitons may be stabilized, if they are unstable without the lattice (multidimensional solitons in the case of self-focusing, as shown in Refs. [17], [19]-[25], see also reviews [1, 3, 4]). The stability of GSs has been studied in detail too—chiefly, close to edges of the corresponding bandgaps—in one [27]-[29] and two [30] dimensions alike.

induced by the PhC makes it possible to create GSs in focusing and defocusing materials alike. Combining assets of PhCs and regular solitons, GSs offer a considerable potential for applications to nonlinear photonics. On the other hand, the use of GSs is limited by the modulational [27] and oscillatory instabilities [35,36]. One of solutions of this problem is the enhancement of the stability (and also mobility of the GSs) in nonlocal nonlinear media [37][38][39][40][41].
In addition to optics, GSs of matter waves have also been theoretically studied [42] and experimentally created [43] in Bose-Einstein condensates formed by atoms with repulsive interactions, loaded into OL potentials. In fact, the OLs controling the dynamics of matter waves may be considered as counterparts of PhCs for coherent atomic wave.
An essential extension of the theme is the study of two-component solitons in lattice potentials. In particular, if both the self-phase-modulation and cross-phase-modulation (SPM and XPM) nonlinearities, i.e., intra-and inter-species interactions, are repulsive, one can construct two-component GSs of intra-gap and inter-gap types, with chemical potentials of the components (or propagation constants, in terms of optical media) falling, respectively, into the same or different bandgaps of the underlying linear spectrum [44,45]. In the case of the attractive SPM, a family of stable semi-gap solitons was found too, with one component residing in the infinite gap, while the other stays in a finite bandgap [45]. The GSs supported by the XPM repulsion dominating over the intrinsic (SPM-mediated) attraction may be regarded as an example of symbiotic solitons. In the free space (without the lattice potential), symbiotic solitons are supported by the XPM attraction between their two components, despite the action of the repulsive SPM in each one [46]- [48]. This mechanism may be additionally enhanced by the linear coupling (interconversion) between the components [49]. Another case of the "symbiosis" was reported in Ref. [50], where the action of the lattice potential on a single component was sufficient for the stabilization of two-dimensional (2D) two-component solitons against the collapse, the stabilizing effect of the lattice on the second component being mediated by the XPM interaction. In addition, the attraction between the components, competing with the intrinsic repulsion, may cause spatial splitting between two components of the GS, as for these components, whose effective masses are negative [16], the attractive interaction potential gives rise to a repulsion force [45,51].
The ultimate form of the model which gives rise to two-component GSs of the symbiotic type is the one with no intra-species nonlinearity, the formation of the GSs being accounted for by the interplay of the repulsion between the components and the lattice potential acting on both of them. In optics, the setting with the XPM-only interactions is known in the form of the "holographic nonlinearity", which can be induced in photorefractive crystals for a pair of coherent beams with a small angle between their wave vectors, giving rise to single- [52,53] and double-peak [54] solitons. Both beams are made by splitting a single laser signal, hence the power ratio between them (which is essential for the analysis reported below) can be varied by changing the splitting conditions. The creation of 2D spatial "holographic solitons" in a photorefractive-photovoltaic crystal with the self-focusing nonlinearity was demonstrated in Ref. [55] (such solitons are stable, as the collapse is arrested by the saturation of the self-focusing). To implement the situation considered here, the sign of the nonlinearity may be switched to self-defocusing by the reversal of the bias voltage, and the effective lattice potential may be induced by implanting appropriate dopants, with the concentration periodically modulated in one direction, which will render the setting quasi-one-dimensional.
In binary bosonic gases, a similar setting may be realized by switching off the SPM nonlinearity with the help of the Feshbach resonance, although one may need to apply two different spatially uniform control fields (one magnetic and one optical) to do it simultaneously in both components. On the other hand, the same setting is natural for a mixture of two fermionic components with the repulsive interaction between them, which may represent two states of the same atomic species, with different values of the total atomic spin (F). If spins of both components are polarized by an external magnetic field, the SPM nonlinearity will be completely suppressed by the Pauli blockade while the inter-component interaction remains active [6], hence the setting may be described by a pair of Schrödinger equations for the two wave functions, coupled by XPM terms.
The objective of this work is to present basic families of one-dimensional symbiotic GSs, supported solely by the repulsive XPM nonlinearity in the combination with the lattice potential, and analyze their stability, via the computation of eigenvalues for small perturbations and direct simulations of the perturbed evolution. The difference from the previously studied models of symbiotic solitons [44,45] is that the solitons where created there by the SPM nonlinearity separately in each component, while the XPM interaction determined the interaction between them and a possibility of creating two-component bound states. Here, the two-component GSs may exist solely due to the repulsive XPM interactions between the components.
We conclude that the symmetric solitons, built of equal components, are destabilized by symmetry-breaking perturbations above a certain critical value of the soliton's power. The analysis is chiefly focused on asymmetric symbiotic GSs, and on breathers into which unstable solitons are transformed. The model is introduced in Section II, which is followed by the analytical approximation presented in Section III. It is an extended version of the Thomas-Fermi approximation, TFA, which may be applied to other models too. In Section IV, we report systematic numerical results obtained for fundamental solitons of both the intra-and inter-gap types, hosted by the first two finite bandgaps of the system's spectrum. The most essential findings are summarized in the form of plots showing the change of the GS stability region with the variation of the degree of asymmetry of the two-component symbiotic solitons, which is a new feature exhibited by the present system. In particular, the stability area of intra-gap solitons shrinks with the increase of the asymmetry, while inter-gap solitons may be stable only if the asymmetry is large enough, in favor of the first-bandgap component, and intra-gap solitons in the second bandgap are completely unstable. The paper is concluded by Section V.

The model
The model outlined above is represented by the system of XPM-coupled Schrödinger equations for local amplitudes of co-propagating electromagnetic waves in the planar optical waveguide, u(x, z) and v(x, z), where x and z are the transverse coordinate and propagation distance, without the SPM terms, and with the lattice potential of depth 2ε > 0 acting on both components: The variables are scaled so as to make the lattice period equal to π, and the coefficients in front of the diffraction and XPM terms equal to 1. In the case of matter waves, u and v are wave functions of the two components, and z is replaced by time t. Direct simulations of Eqs.
(1) and (2) were performed with the help of the split-step Fourier-transform technique.
Stationary solutions to Eqs. (1), (2) are looked as where the real propagation constants, k and q, are different, in the general case, and real functions U(x) and V(x) obey equations with the prime standing for d/dx. Numerical solutions to Eqs. (4) and (5) were obtained by means of the Newton's method.
Solitons are characterized by the total power, with both P u and P v being dynamical invariants of Eqs. (1), (2), and by the asymmetry ratio, The total power and asymmetry may be naturally considered as functions of the propagation constants, k and q.
The well-known bandgap spectrum of the linearized version of Eqs. (4), (5) (see, e.g., book [12]) is displayed in Fig. 1, the right edge of the first finite bandgap being The location of GSs is identified with respect to bandgaps of the spectrum. In this work, results are reported for composite GSs whose two components belong to the first and second finite bandgaps.
Stability of the stationary solutions can be investigated by means of the linearization against small perturbations [26]- [28]. To this end, perturbed solutions of Eqs. (1) and (2) are looked for as where u 1,2 and v 1,2 are wave functions of infinitesimal perturbations, and λ is the respective instability growth rate, which may be complex (the asterisk stands for the complex conjugate). The instability takes place if there is at least one eigenvalue with Im(λ) > 0. The substitution of ansatz (9) into Eqs. (1), (2) and the linearization with respect to the small perturbations leads to the eigenvalue problem based on the following equations: These equations were solved by means of the fourth-order center-difference numerical scheme.
Results for the shape and stability of GSs of different types are presented below for lattice strength ε = 6, which adequately represents the generic case. Note, in particular, that Fig. 1 was plotted for this value of the lattice-potential's strength.

The extended Thomas-Fermi approximation
It is well known that, close to edges of the bandgap, GSs feature an undulating shape, which may be approximated by a Bloch wave function modulated by a slowly varying envelope [14,16]. On the other hand, deeper inside the bandgap, the GSs are strongly localized (see, e.g., Figs. 4 and 7 below), which suggests to approximate them by means of the variational method based on the Gaussian ansatz [45,56]. This approximation was quite efficient for the description of GSs in single-component models [56], while for two-component systems it becomes cumbersome [44,45]. Explicit analytical results for well-localized patterns can be obtained by means of the TFA [5], which, in the simplest case, neglects the kinetic-energy terms, U ′′ and V ′′ , in Eqs. (4), (5). Assuming, for the sake of the definiteness, q < k and also |k| < ε (the TFA is irrelevant for |k| > ε), the approximation yields the fields inside the inner layer of the solution: Thus, the TFA predicts the core part of the solution in the form of peaks in the two components with the same width, 2x 0 , but different heights, Further, the expansion of expressions (14) around the soliton's center (x = 0), yields The substitution of the second derivatives of the fields at x = 0, calculated as per Eq. (15), into Eqs. (4), (5) yields a corrected expression for the soliton's amplitudes: along with the condition for the applicability of the TFA: For the symmetric-GS families in the two first finite bandgaps, the amplitude predicted by the improved TFA in the form of Eq. (16) is displayed, as a function of k = q, in Fig. 2(a) and compared to its numerically found counterpart. It is worthy to note that the correction terms in Eq. (16) essentially improve the agreement of the TFA prediction with the numerical findings.
At |x| > x 0 , the TFA gives V(x) = 0, which is a continuous extension of the respective expression (14) in the V component, while the continuity of fields U(x) and U ′ (x) makes it necessary to match the respective expression (14) to "tails", which, in the lowest approximation, satisfy equation U ′′ = 0. The continuity is provided by the following tail solution:  Total power P for the same soliton families, whose stable and unstable portions are designated by the bold dotted and dashed lines, respectively. The latter one is destabilized by symmetry-breaking perturbations, while the entire family is stable in the framework of the single-component model. Coordinates of the stability/instability border are given by Eq. (21). The chain of squares shows the analytical dependence produced by the TFA, see Eq. (20).
The integration of expressions (14) and (18) yields the following approximation for the powers of the two components: The substitution of approximation (19) into definitions (6) and (7) of the total power and asymmetry demonstrates an agreement with numerical results. For instance, the slope of the curve R(q) for the intra-gap GSs at fixed k (see Fig. 8 below) at the symmetry point (k = q), as predicted by Eq. (19) for ε = 6 and k = 1, is (∂R/∂q) | q=k ≈ −0.155, while its numerically found counterpart is ≈ −0.160. Further, the analysis of Eq. (19) readily demonstrates that the strongly asymmetric solitons may exist up to the limit of R → 1, which is corroborated by the existence area for the intra-gap solitons shown below in Fig. 11(b).
Another corollary of Eqs. (19) is the prediction for the total power for the symmetric solitons, which is plotted in Fig. 2(b), along with its numerically found counterpart. Although the TFA does not predict edges of the bandgaps, the overall analytical prediction for P(k), as well as the soliton's amplitude shown in Fig. 2(a), run quite close to the numerical curves, except for near the right edge, where, indeed, condition (17) does not hold for ε = 6 and k = 3.75, see Eq. (8). Note that these simple analytical approximations were not derived before in numerous works dealing with single-component GSs.

Symmetric solitons
Obviously, the shape of symmetric solitons, built of two equal components [with k = q and U(x) = V(x), see Eq. (3)], is identical to that of GSs in the single-component model. However, there is a drastic difference in the stability of the symmetric GSs between the single-and two-component systems. Almost the entire symmetric family is stable against symmetric perturbations, i.e., it is stable in the framework of the single-component equation (in agreement with previously known results [12]), except for a weak oscillatory instability, accounted for by quartets of complex-conjugate eigenvalues, in the form of λ = ±iIm (λ) ± Re (λ) (with two mutually independent signs ±), which appears near the left edge of the second bandgap-namely, at k < k min ≈ −3.45. An example of the development of the latter instability is displayed in Fig. 3. On the other hand, Fig. 2 demonstrates that a considerable part of the family in the first finite bandgap, and the entire family in the second bandgap are unstable against symmetry-breaking perturbations in the two-component system. The boundary separating the stable and unstable subfamilies of the fundamental symmetric GSs in the first finite bandgap corresponds to the power and propagation constants is found at the symmetric solitons being stable in the intervals of 0 < P < 7.19, 1.05 < k < k max  Typical examples of stable and unstable fundamental symmetric GSs, found in the first and second bandgaps (not too close to their edges), are displayed in Fig. 4. Further, direct simulations demonstrate that the evolution transforms the unstable symmetric solitons into persistent localized breathers, as shown in Figs. 5 and 6, in accordance with the fact that the corresponding instability eigenvalues are complex. Although the emerging breather keeps the value of R = 0, see Eq. (7), the u-and v-components of the breather generated by the symmetry-breaking instability are no longer mutually identical. This manifestation of the symmetry-breaking instability is illustrated by Fig. 5(c), which displays the evolution of the difference between the peak powers of the two components, and the separation between their centers. The latter is defined as It is relevant to mention that the second finite bandgap also contains a branch of the so-called subfundamental solitons, whose power is smaller than that of the fundamental GSs [29], [57], [58]. These are odd modes, squeezed, essentially, into a single cell of the underlying lattice potential. The subfundamental solitons are unstable, tending to rearrange themselves into fundamental ones belonging to the first finite bandgap, therefore they are not considered below.

Asymmetric solitons of the intra-gap type
As said above, two-component asymmetric fundamental GSs, with different propagation constants, k = q, may be naturally classified as solitons of the intra-and inter-gap types if k and q belong to the same or different finite bandgaps [44]. In this subsection, we report results for asymmetric intra-gap solitons with both k and q falling into the first finite bandgap, as well as for asymmetric breathers developing from such solitons when they are unstable.
Examples of stable asymmetric GSs of the intra-gap type are displayed in Fig. 7, for a fixed asymmetry ratio, R = −0.5, defined as per Eq. (7). The GS family, along with the family of persistent breathers into which unstable solitons are spontaneously transformed, is The top and bottom plots display, respectively, the evolution of the peak-power difference, max(|u (x, z)| 2 ) − max(|u (x, z)| 2 ), and the separation between centers of the two components, X u and X v , which is defined as per Eq. (22).  Fig. 5(a,b), but for an unstable soliton in the second bandgap, with k = q = −3.5. It is possible to explain the fact that all the R(q) curves converge to R = −1, as q approaches the right edge of the bandgap in Fig. 8. In this case, the V component turns into the delocalized Bloch wave function with a diverging power, P v , that corresponds to P u /P v → 0 [it is tantamount to R → −1, as per Eq. (7)]. An example of a stable GS, close to this limit, with k = −0.5, q = 3.65 and R = −0.9811, is shown in Fig. 9. The central core of the V-component is described by the TFA, based on Eq. (14), as the corresponding necessary condition (17) holds in this case, while the TFA does not apply to the U-component. The presence of undulating tails, which are close to the Bloch functions, rather than the simple approximation (18), which is valid far from the edge of the bandgap, is also visible in Fig. 8.
Those asymmetric intra-gap GSs which form unstable subfamilies in Fig. 8 are destabilized by oscillatory perturbations. The instability transform the solitons into breathers, see a typical V(x), which pertain to propagation constants k and q, are shown, respectively, by the magenta (lower) and blue (taller) profiles.
example in Fig. 10 [cf. the examples of the destabilization of the symmetric GSs shown in Fig. 5(b,c)]. The emerging breathers keep values of the asymmetry ratio (7) almost identical to those of their parent GSs; for instance, in the case displayed in Fig. 10, the unstable soliton with R initial = −0.3617 evolves into the breather with R final = −0.3623.
It is relevant to stress that the transformation of unstable stationary GSs into the breathers gives rise to little radiation loss of the total power, P. On the other hand, in the general case a given unstable gap soliton does not have a stable counterpart with a close value of P, hence this unstable soliton cannot transform itself into a slightly excited state of another stable GS. Thus, the breathers represent a distinct species of localized modes.
The most essential results of the stability analysis for the asymmetric solitons of the intra-gap type, and for breathers replacing unstable solitons, are summarized by diagrams in the planes of (k, q) and (P, R), which are displayed in Fig. 11. The predictions of the analysis based on the computation of the stability eigenvalues for the stationary solitons, as per Eqs. (10)-(13), always comply with stability tests provided by direct simulations of Eqs. (1) and (2).
As mentioned above, the instability of a part of the branch of the symmetric solitons along the line of R = 0 in Fig. 11(b) implies that the symmetry-breaking perturbations destabilize the symmetric solitons in the first finite bandgap at P > P cr , see Eq. (21), while their counterparts are stable in the single-component system. Another clear conclusion is that the stability region gradually shrinks with the increase of the asymmetry.

Solitons of the inter-gap type
All The asymmetry measure, R(q), for families of the inter-gap solitons is plotted in Fig. 14 The (in)stability of the inter-gap solitons is summarized by the diagrams in the planes of (k, q) and (P, R) presented in Fig. 15, cf. similar diagrams for intra-gap solitons shown above in Fig. 11. As well as in that case, unstable inter-gap solitons are spontaneously replaced by robust localized breathers. The spontaneous transformation increases the initial degree of the asymmetry: For instance, an unstable inter-gap soliton with R = 0.019 is converted into a breather with R = 0.028.
In the present case too, the existence region of stable modes shrinks with the increase of the asymmetry; note also that the stationary inter-gap solitons may be stable solely at sufficiently large values of the asymmetry, R ≥ R min ≈ 0.5. The asymmetric shape of the stability diagram in Fig. 15(b) with respect to R > 0 and R < 0 [unlike the symmetry of the diagram for the intra-gap solitons implied in Fig. 11(b)] is explained by the fact that, in definition (7), R > 0 implies that the dominant component resides in the first finite bandgap, where it is more robust than in the second bandgap.
Finally, out additional analysis has demonstrated that all the stationary GSs-not only the symmetric ones [see Fig. 2], but also all the asymmetric solitons of the intra-gap type-are completely unstable in the second finite bandgap. They too tend to spontaneously rearrange themselves into breathers, which is not shown here in detail.    Fig. 11, but for inter-gap solitons. In (a), the dotted rectangle delineates the region occupied by wavenumbers k and q belonging to the first and second finite bandgaps, respectively.

Conclusion
We have introduced the model of symbiotic two-component GSs (gap solitons), based on two nonlinear Schrödinger equations coupled by the repulsive XPM terms and including the lattice potential acting on both components, in the absence of the SPM nonlinearity. The model has a realization in optics, in terms of "holographic solitons" in photonic crystals, and as a model of binary quantum gases (in particular, a fully polarized fermionic one) loaded into the optical-lattice potential. Families of fundamental asymmetric GSs have been constructed in the two lowest finite bandgaps, including the modes of both the intra-gap and inter-gap types, i.e., those with the propagation constants of the two components belonging to the same or different bandgaps, respectively. The existence and stability regions of the symbiotic GSs and breathers, into which unstable solitons are transformed, have been identified. A noteworthy finding is that symmetry-breaking perturbations destabilize the symmetric GSs in the first finite bandgap, if their total power exceeds the critical value given by Eq. (21), along with all the symmetric solitons in the second bandgaps. It was demonstrated too that the stability area for the intra-gap GSs shrinks with the increase of the asymmetry ratio, R. On the other hand, inter-gap GSs may be stable only for sufficiently large ratio, R > 0.5. The intra-gap solitons are completely unstable in the second bandgap. Some features of the GS families were explained by means of the extended TFA (Thomas-Fermi approximation), augmented by the tails attached to the taller component, in the case of asymmetric solitons.
A natural extension of the analysis may deal with 2D symbiotic gap solitons, supported by the square-or radial-lattice potentials. In that case, it may be interesting to consider two-component solitary vortices too.