Quiescent Gap Solitons in Coupled Nonuniform Bragg Gratings with Cubic-Quintic Nonlinearity

: We study the stability characteristics of zero-velocity gap solitons in dual-core Bragg gratings with cubic-quintic nonlinearity and dispersive reﬂectivity. The model supports two disjointed families of gap solitons (Type 1 and Type 2). Additionally, asymmetric and symmetric solitons exist in both Type 1 and Type 2 families. A comprehensive numerical stability analysis is performed to analyze the stability of solitons. It is found that dispersive reﬂectivity improves the stability of both types of solitons. Nontrivial stability boundaries have been identiﬁed within the bandgap for each family of solitons. The effects and interplay of dispersive reﬂectivity and the coupling coefﬁcient on the stability regions are also analyzed.


Introduction
In the past few decades, significant attention has been given to fiber Bragg gratings (FBGs) owing to their applications in filtering, data format conversion, and dispersion compensation [1][2][3][4][5][6][7]. These applications have been realized, in part, due to a distinctive property of FBGs, namely the existence of a stopband (or photonic bandgap) in their linear spectrum [3,8]. Additionally, FBGs possess a strong dispersion that can be a million times larger than the dispersion of silica [8,9].
In the case of FBGs written in nonlinear media (e.g., Kerr nonlinearity), the nonlinearity can be counterbalanced by the FBG's effective dispersion giving rise to the generation of gap solitons (GSs) [9][10][11]. Gap solitons have received much attention both theoretically and experimentally due to their potential applications in slow light applications, optical memory elements [12,13], and optical logic operations [14,15].
Bragg grating solitons have extensively been studied in the cubic nonlinearity both theoretically [16][17][18][19] and experimentally [20][21][22][23]. It has been established that gap solitons in uniform Bragg gratings and Kerr nonlinearity exist throughout the bandgap and can be characterized by two parameters. Moreover, these gap solitons can travel at any velocity between zero (quiescent) and c/n where c is the speed of light in vacuum and n is the refractive index of the medium. As for stability of the gap solitons, it has been shown analytically and numerically that about 50% of the soliton family is stable [24,25]. Experimentally, gap solitons with a speed of 23% of the speed of light have been observed [23].
The standard model of gap solitons assumes that the FBG is uniform and weak. However, for specialized gratings such as deep gratings, gratings based on photonic nanowires [37] and Bragg superstructures [38,39], due to the inhomogeneity of the bandgap spectra, the effect of dispersion of reflectivity needs to be incorporated into the model [40]. Numerical simulations have demonstrated that dispersive reflectivity improves the stability of the moving and quiescent gap solitons and results in the enlargement of the stable region [40,41].
Soliton switching and multicomponent solitons in uniform optical media and structures such as nonlinear couplers are potential building blocks for optical signal processing devices [42][43][44][45]. Indeed, the ability to perform optical switching and/or signal processing with gap solitons is an intriguing one. Such optical devices can be constructed by writing Bragg gratings in one or both cores of a fiber coupler [46]. As a result, gap solitons in grating-assisted couplers with identical and dissimilar Bragg gratings are of great interest [30,31,47]. It has been shown that generalization of these systems to include dispersive reflectivity results in a much richer dynamics which can be instrumental in designing novel optical devices [48].
In Reference [47], the stability of GSs in coupled identical Bragg gratings with cubicquintic nonlinearity was analyzed. In the present work, a generalization of the model of Reference [47] is considered which includes the effect of dispersive reflectivity in both gratings. We analyze the characteristics and dynamical stability of zero-velocity GSs in the generalized model. In Section 2, the generalized model and the features of the bandgap are discussed. In Section 3, the stationary quiescent soliton solutions and their characteristics are analyzed. The analysis of the stability of the quiescent GSs is presented in Section 4. A summary of the main results is provided in Section 5.

The Model and Its Linear Spectrum
To derive the mathematical model for the system, we start from Maxwell's equations and following the steps for the derivation of GMTM [9], the coupled mode equations [49], and the approaches for the inclusion of dispersive reflectivity and cubic-quintic nonlinearity described, respectively, in References [36,40] the following system of normalized partial differential equations can be derived: In Equation (1), forward-propagating waves in core 1 and 2 are represented by u 1,2 (x, t) and the backward propagating ones are denoted by v 1,2 (x, t). m > 0 represents the strength of the dispersive reflectivity, η > 0 is the normalized quintic nonlinearity coefficient, and λ > 0 is the coefficient of the linear coupling between the cores. Based on the experimental values measured for the cubic-quintic nonlinearity for chalcogenide glasses [50,51] and some organic materials [52,53], it is found that η can range between 0.05 and 0.62. As a result, our analysis is confined to the range 0 ≤ η ≤ 1. In addition, as was noted in Reference [40], m > 0.5 may not be physically significant. Therefore, the analysis will be limited to 0 ≤ m ≤ 0.5.
To determine the linear spectrum of the system, u 1,2 , v 1,2 ∼exp(i(kx − ωt)) are inserted into Equation (1). Linearizing the resulting equations and some straightforward algebraic manipulations lead to the following dispersion relation: The negative(positive) sign in (2) refers to the lower(upper) branch of the dispersion relation. When m = 0, the dispersion diagram relation reduces to the one in Reference [47]. The size of the bandgap can be determined from the roots of dω 2 dk 2 = 0. This leads to the following bandgap: where ω 2 0 is the minimum value of the ω 2 (k 2 ). For m ≤ 0.5, ω 0 is attained at k 2 0 = 0. On the other hand, for m > 0.5, ω 0 occurs at It is worth noting that for 0 ≤ m ≤ 0.5, the bandgap does not depend on m and it closes at λ = 1. Figure 1 shows the dispersion curves for different values of λ and m.

Quiescent Gap Soliton Solutions
For m = 0, Equation (1) simplify to the model of Reference [47] which admits exact analytical soliton solutions. However, for m = 0, there is no analytical solution to Equation (1) and the soliton solutions have to be obtained via numerical methods. To this end, we assume that the quiescent soliton solutions can be written as Substituting these expressions into Equation (1) and simplifying leads to the following system of differential equations: The relaxation algorithm is used to solve Equation (3). The relaxation algorithm is a numerical technique for solving boundary value problems. The algorithm employed in our analysis utilizes a deferred correction technique together with Newton iteration (for further details on the relaxation algorithm, see References [54,55]). Similar to the model of Reference [47], two disjointed families of zero-velocity solitons are found in the model of Equation (1), namely, Type 1 and Type 2 solitons. Furthermore, both asymmetric (U 1 = U 2 , V 1 = V 2 ) and symmetric (U 1 = U 2 , V 1 = V 2 ) solitons exist in these soliton families. Type 1 and Type 2 solitons differ in their phase structure and the amplitudes of the solitons. More specifically, for both symmetric and asymmetric Type 1 solitons, the real(imaginary) parts of the forward and backward waves are even(odd) functions of x. However, for symmetric and asymmetric Type 2 solitons, the real and imaginary parts of the the forward and backward waves have opposite parities compared with those of the Type 1 ones. The boundary that separates the soliton families has to be found numerically. It is found that the boundary is approximated very well by the expression . It should be noted that no soliton solutions exist on the boundary.
Typical profiles for Type 1 and Type 2 symmetric solitons are displayed in Figure 2. In the single core and dual-core gap soliton models without dispersive reflectivity (see References [36,47]), the Type 2 solitons have a sharp (but nonsingular) peak and close to the boundary between the soliton families, they possess a two-tier structure. However, our analysis shows that the profiles of Type 2 asymmetric and symmetric solitons in the present model are drastically different in that their peaks are not sharp and the two-tier structure is absent. Another observation is that, similar to the other models with dispersive reflectivity, there exists a critical value of m beyond which the solitons develop sidebands. More specifically, the sidebands appear in the profiles of the solitons if the following inequality is satisfied [56] m > 1 Examples of soliton profiles with sidebands are shown in Figure 3. Similar to the models without dispersive reflectivity [30,47], the analysis of Type 1 symmetric and asymmetric soliton solutions shows that at any ω, there exists a certain coupling coefficient, i.e., ±λ c , at which the symmetric branch bifurcates into three branches. Two of these branches that are mirror images of one another represent the stable Type 1 asymmetric solitons and the remaining branch which is the extension of the original symmetric branch corresponds to the unstable symmetric Type 1 solitons.  Bifurcation diagrams can be constructed by plotting λ versus the effective asymmetry parameter θ which is given by where U 1max and U 2max represent the peak amplitudes of the u-components in core 1 and 2, respectively. Bifurcation plots for various values of m, ω, and η are shown in Figure 4. It is worth noting that bifurcation does not occur for Type 2 solitons. As a result, for all values of ω within the bandgap, Type 2 asymmetric and symmetric solitons always exist together.

Results of the Stability Analysis and Discussion
The aforementioned bifurcation diagrams can be utilized to make conjectures regarding the stability of Type 1 symmetric and asymmetric solitons [32,57]. According to the elementary bifurcation theory, for those values of λ for which the asymmetric and symmetric solitons exist together (i.e., when λ is in the range |λ| < λ c ), the symmetric solitons (i.e., solitons corresponding to the dotted lines in Figure 4) should be unstable and the asymmetric ones should be stable. When |λ| > λ c , there are only symmetric solitons in which case the bifurcation theory predicts that the symmetric solitons would be stable. Since bifurcation diagrams cannot be constructed for Type 2 solitons (see Section 3), their stability must be determined numerically.
To determine the stable regions within the bandgap and to check the validity of the conclusions of the bifurcation theory, we have performed a systematic and extensive numerical stability analysis of the symmetric and asymmetric Type 1 and Type 2 solitons. To this end, a symmetrized split-step Fourier scheme has been employed (see Reference [58]) to solve Equation (1). In addition, to prevent reflection of radiation, absorbing boundary conditions have been incorporated into our simulations. The numerically-obtained solutions were used as the initial configuration for the numerical simulations. The evolution of solitons was simulated up to t = 2000. The noise inherent in the numerical scheme was found to be adequate for the onset of the instability development. The results of our numerical simulations fully corroborate the predictions of the bifurcation theory. As is shown in Figure 5, for the parameters where the Type 1 symmetric and asymmetric GSs exist together, the asymmetric soliton is stable (Figure 5a) and the symmetric one is unstable (Figure 5b,c). Additionally, Figure 5b,c demonstrate that the evolution of the unstable Type 1 symmetric solitons is accompanied by energy loss as radiation and eventually they evolve to a Type 1 asymmetric soliton. As regards the stability of Type 2 gap solitons, the simulations show that there is a region within the bandgap where Type 2 asymmetric solitons are stable. An example of the stable Type 2 asymmetric soliton is shown in Figure 6a. Type 2 symmetric solitons, on the other hand, are found to be always unstable. An example of the evolution of a Type 2 symmetric soliton is displayed in Figure 6b,c. It is noteworthy that the symmetric Type 2 soliton evolves to an asymmetric moving Type 2 soliton. This confirms that stable moving Type 2 asymmetric solitons are supported by Equation (1). The propagation of unstable Type 1 and Type 2 asymmetric solitons might result in different outcomes. For example, as is shown in Figure 7a, the asymmetric unstable Type 1 soliton sheds some of its energy during the propagation and then rearranges itself into a Type 1 asymmetric moving stable soliton. On the other hand, the evolution of Type 2 asymmetric solitons might result in either small amount of radiation followed by the formation of another quiescent stable Type 2 asymmetric soliton (Figure 7b) or complete destruction of the soliton (Figure 7c). One of the goals of the stability analysis is to identify the stable regions inside the bandgap for both types of solitons. In addition, it is necessary to analyze how the system parameters (i.e., λ and m) affect the stability of solitons. To this end, we have carried out the stability analysis for various values of m and λ. The results of these stability analyses are displayed in the stability diagrams of Figure 8. Compared with the model of Reference [47] (i.e., m = 0), the stability diagrams in Figure 8 show that the stable regions are significantly altered due to the dispersive reflectivity. In particular, m = 0 gives rise to the enlargement of the stable regions for asymmetric Type 2 solitons and symmetric Type 1 solitons. Another notable feature of Figure 8 is that for a given λ, if m is varied from 0.1 to 0.5, the sizes of the stable regions for asymmetric Type 2 and symmetric Type 1 solitons progressively increase. In the case of λ = 0.1, as m increases in the range 0.1 ≤ m ≤ 0.3, the stable region for Type 1 asymmetric solitons becomes larger. However, for values of m > 0.3 this trend is reversed (see the third diagram in the first row of Figure 8). Figure 8 shows that, for a fixed m, increasing the coupling coefficient leads to significant enlargement of the stable region for Type 1 symmetric solitons. This, however, results in the shrinking of the stable region for Type 1 asymmetric solitons. In addition, as λ increases, the stable region for asymmetric Type 2 solitons is shrunk. In other words, for the asymmetric solitons of both types, stronger coupling coefficient counteracts the stabilizing effect of the dispersive reflectivity.

Conclusions
We have analyzed the stability dynamics of zero-velocity gap solitons in coupled FBGs with dispersion of reflectivity and cubic-quintic nonlinearity. The model supports two disjointed families of GSs which are categorized as Type 1 and Type 2 solitons. Type 1 solitons are the generalizations of the gap solitons in coupled Bragg gratings in the third order (cubic) nonlinearity with dispersive reflectivity. For Type 2 solitons, the effect of quintic nonlinearity is more pronounced. A key finding is that, unlike the cubic-quintic models without dispersive reflectivity, the Type 2 solitons' profiles are smooth and they do not have a sharp peak.
The stability of Type 1 and Type 2 asymmetric and symmetric solitons has been investigated systematically via numerical stability analysis. The stability analysis corroborates the predictions of the bifurcation theory. More specifically, it is found that the symmetric solitons are stable only in the parametric region where they exist by themselves. Through the stability analysis, the stable and unstable regions for Type 1 and Type 2 asymmetric solitons have been identified. A notable finding is that, for a fixed coupling coefficient, increasing the dispersive reflectivity enlarges the stable region for the asymmetric Type 2 solitons. On the other hand, for a given value of dispersive reflectivity, increasing the coupling coefficient makes the stable regions for Type 1 and Type 2 asymmetric solitons smaller.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.