Coupled backward- and forward-propagating solitons in a composite right- and left-handed transmission line

PHYSICAL REVIEW E 88, 013203 (2013) Coupled backward- and forward-propagating solitons in a composite right- and left-handed transmission line G. P. Veldes,1,2J. Cuevas,3P. G. Kevrekidis,4and D. J. Frantzeskakis1 1Department of Physics, University of Athens, Panepistimiopolis, Zografos, Athens 15784, Greece 2Department of Electronics, Technological Educational Institute of Lamia, Lamia 35100,Greece 3Grupo de F´ ısica No Lineal, Universidad de Sevilla. Departamento de F´ ısica Aplicada I, Escuela Polit´ ecnica Superior, C/ Virgen de ´ Africa, 7, 41011 Sevilla, Spain 4Department of Mathematics and Statistics, University of Massachusetts, Amherst Massachusetts 01003-4515, USA (Received 13 February 2013; revised manuscript received 28 May 2013; published 12 July 2013) We study the coupling between backward- and forward-propagating wave modes, with the same group velocity, in a composite right- and left-handed nonlinear transmission line. Using an asymptotic multiscale expansion technique, we derive a system of two coupled nonlinear Schr¨ odinger equations governing the evolution of the envelopes of these modes. We show that this system supports a variety of backward- and forward-propagating vector solitons of the bright-bright, bright-dark, and dark-bright type. Performing systematic numerical simulations in the framework of the original lattice that models the transmission line, we study the propagation properties of the derived vector soliton solutions. We show that all types of the predicted solitons exist, but differ on their robustness: Only bright-bright solitons propagate undistorted for long times, while the other types are less robust, featuring shorter lifetimes. In all cases, our analytical predictions are in very good agreement with the results of the simulations, at least up to times of the order of the solitons’ lifetimes. DOI: 10.1103/PhysRevE.88.013203 PACS number(s): 41.20.Jb, 42.65.Tg, 78.20.Ci I. INTRODUCTION Left-handed (LH) metamaterials are artificial, effectively homogeneous structures, featuring negative refractive index at specific frequency bands where the effective permittivity  and permeability μare simultaneously negative [1–3]. In fact, all known realizations of LH metamaterials rely on the use of common right-handed (RH) elements and thus in a realistic situation such a composite material features both LH and RH behavior in certain frequency bands. Physically speaking, the difference between the two is that in the LH (RH) regime, the energy and the wave fronts of the electromagnetic (em) waves propagate in the opposite (same) directions, giving rise to backward- (forward-) propagating waves. Transmission line (TL) theory constitutes a convenient framework for the analysis of LH metamaterials. Such an analysis relies on the connection of the em properties of the medium (and μ) with the electric elements of the TL’s unit cell, namely, the serial and shunt impedance. As mentioned above, in practice, composite right- and left-handed (CRLH) structures are quite relevant, giving rise to pertinent CRLH-TL models. These models are in fact dynamical lattices that can be used for the description of a variety of metamaterial-based devices and systems, such as resonators, directional couplers, and antennas [1–4]. Nonlinear CRLH-TL models, with a serial and/or shunt impedance depending on voltages or currents, have also attracted attention. Such structures may be realized by inserting diodes, which mimic voltage-controlled nonlinear capacitors, into resonant conductive elements (such as split-ring resonators) [5–7]. Such nonlinear CRLH-TL models have been used in various works dealing, e.g., with the parametric shielding of em fields [8], the long-wave–short-wave interaction [9], or soliton formation [10–12]. Experiments in nonlinear CRLH-TL models have also been performed (see [13]) and the formation of bright [14,15] and dark [15,16] envelope solitons, described by an effective nonlinear Schr¨ odinger (NLS) equation, was reported. Notice that in earlier studies on RH-TL models it was shown that two (or more) solitons propagating with the same group velocity can be described by a system of two (or more) NLS equations [17](seealso[18]for theoretical as well as experimental results). Such coupled NLS equations have been studied extensively in nonlinear optics and mathematical physics; see, e.g., Refs. [19–21] and references therein. They are well known to give rise to a variety of vector solitons, including bright-bright (BB), bright-dark (BD), and dark-dark (DD) ones. In this work we study analytically and numerically the interaction between backward- and forward-propagating solitons in a nonlinear CRLH-TL model. Our model is a nonlinear version of a generic CRLH-TL model, particularly relevant to the context of LH metamaterials (see, e.g., Refs. [2,4]). The considered nonlinear element in the unit cell of the TL is the shunt capacitor, which simulates the presence of a heterostructure barrier varactor (HBV) diode [6](the capacitance of the HBV diode depends on the applied voltage). Starting from the discrete lump element model of the CRLH TL, we derive a nonlinear lattice equation. First, we study the linear regime and show that for certain frequency bands, RH and LH modes can propagate with the same group velocity. Then we treat the nonlinear lattice equation in the framework of the quasidiscrete (or quasicontinuum) approximation (see, e.g., [15,19,22,23] for a review): we thus seek envelope soliton solutions of the nonlinear lattice model, characterized by a discrete carrier and a continuum envelope and employ an asymptotic multiscale expansion method to derive a system of two coupled NLS equations. Each of these equations describes the evolution of the envelope of a backward- (LH) and a forward-propagating (RH) mode. 013203-1 1539-3755/2013/88(1)/013203(14) ©2013 American Physical Society VELDES, CUEVAS, KEVREKIDIS, AND FRANTZESKAKIS PHYSICAL REVIEW E 88, 013203 (2013) A systematic analysis of the system of the NLS equations reveals the existence, in certain frequency bands, of three different types of vector solitons: (a) a backward-propagating bright soliton coupled with a forward-propagating bright soliton, (b) a backward-propagating bright soliton coupled with a forward-propagating dark soliton, and (c) a backwardpropagating dark soliton coupled with a forward-propagating bright soliton. This way we propose and study in detail all possible vector solitons in this nonlinear CRLH-TL setting, highlighting how individual soliton components interact via the nonlinearity induced by the insertion of the HBV diodes. Importantly, our analysis suggests how the characteristics (amplitude, width, etc.) of each solitonic mode can control the ones of the other mode; this possibility arises from the fact that the relevant parameters of each soliton are connected with the ones of the other soliton, thus providing a means to adjust, e.g., the amplitude or width of the LH mode by changing the parameters of the RH mode. The above analytical predictions are then tested against direct numerical simulations, which are performed in the framework of the original nonlinear lattice model. The results of the simulations verify the existence of the aforementioned types of vector solitons in the full TL model, but also offer important information regarding their robustness. In particular, the results of direct simulations performed for long times indicate that BB solitons are the most robust among the members of the vector soliton family. Indeed, the mixed [dark-bright (DB) or BD] types are found to be less robust; however, the DB solitons in a specific frequency band, although deformed during their evolution, are found to be more robust than those in other bands, as well as the bright-dark solitons, which are destroyed for the same propagation time. The observed long-time behavior of the above solitons of the mixed type may be qualitatively explained, in the framework of the effective coupled NLS description, by the following fact: in the less robust cases, the continuous-wave background carrying the dark soliton is prone to modulational instability [18,24], while in the more robust case it is not. In any case, our results indicate the existence of all three types: robustness of BB solitons and partial or substantial deformation of the other types. We can thus conclude that bright-bright (LH-RH) solitons, as well as dark-bright (LH-RH) solitons in certain frequency bands, have a better chance of being observed in experiments. The paper is organized as follows. In Sec. II we introduce the nonlinear CRLH-TL model and the pertinent lattice equation and derive the system of the two coupled NLS equations (relevant details are also appended in the Appendices). In Sec. III we present analytical and numerical results for each type of vector soliton. Finally, in Sec. IV we summarize and give our conclusions. II. THE MODEL AND ITS ANALYTICAL CONSIDERATIONS A. The nonlinear CRLH-TL model We consider a generic CRLH-TL model, composed by both right- and left-handed elements, as shown in its unit-cell circuit in Fig. 1[2,4]. The (RH) elements of this TL are the inductance FIG. 1. Unit-cell circuit of the nonlinear CRLH model. LRand capacitance CR, while the LH ones are the inductance LLand capacitance CL. We assume that the TL is loaded with a nonlinear capacitance (CR, while the capacitance CL will be assumed to be fixed and voltage independent). This can be implemented by proper insertion of diodes in the TL (see, e.g., pertinent experiments as well as theoretical work in Refs. [10–16]); in other words, we assume that the shunt capacitor CRis nonlinear (see details below). Let us now consider Kirchhoff’s voltage and current laws for the unit-cell circuit of Fig. 1, which respectively read Vn−1=Vn+LR dIn dt +Un,(1) In=In+1+IL+d dt (CRVn),(2) where Unis the voltage across the capacitance CLand IL is the current across the inductor LL. The above equations, together with the auxiliary equations Vn=LLdIL/dt and In=CLdUn/dt, lead to the following system: LRLLCL d4 dt4(CRVn)+LL d2 dt2(CRVn)+LRCL d2Vn dt2 −LLCL d2 dt2(Vn+1+Vn−1−2Vn)+Vn=0.(3) To proceed further, we now consider a specific voltage dependence for the nonlinear capacitance CR.Herewewill assume that for sufficiently small values of the voltage Vn,the function CR(Vn) can be approximated as follows, via a Taylor expansion: CR(Vn)≈CR0+C R0(Vn−V0)+1 2C R0(Vn−V0)2,(4) where CR0≡CR(V0) is a constant capacitance corresponding to the bias voltage V0, while C R0and C R0also assume constant values, depending on the particular form of CR(V). Below we will further discuss this approximation, in connection with the HBV diode, used in the experiments described in Ref. [13] (similar varactor-type diodes were also used in the experiments of Ref. [7]). Next, substituting Eq. (4) into Eq. (3) and using the scale transformations t→ωsht[where ω2 sh =(LLCR0)−1] and Vn→[C R0(2CR0)−1]Vn, we obtain d4Vn dt4−β2d2 dt2(Vn+1+Vn−1−2Vn)+(1 +δ2)d2Vn dt2 +δ2Vn+δ2d2V2 n dt2+δ2μd2V3 n dt2+d4V2 n dt4+μd4V3 n dt4=0, (5) where the constant parameters δ,β, and μare given by δ=fse fsh ,β=fRH fsh ,μ=2C R0 3C2 R0 CR0.(6) 013203-2 COUPLED BACKWARD- AND FORWARD-PROPAGATING ... PHYSICAL REVIEW E 88, 013203 (2013) In the above expressions, fse and fsh denote series and shunt frequencies, respectively, and fRH denotes the characteristic frequency related to the RH part of the unit-cell circuit; the above frequencies are defined as fse =1 2π√LRCL ,f sh =1 2π√LLCR0 , fRH =1 2π√LRCR0 .(7) Note that if fse/fsh =1, i.e., δ=1, then the CRLH-TL model is usually referred to as balanced, in the sense that the characteristic impedances of the purely LH and RH TLs, defined as ZL=√LL/CLand ZR=√LR/CR0, are equal, i.e., ZL=ZR[2]. In contrast, if fse/fsh >1, i.e., δ>1, the LH part of the TL dominates, in the sense that the TL has a more pronounced LH behavior (the serial branch features a capacitive character, while the shunt branch an inductive one). In the opposite case, fse/fsh <1, i.e., δ<1, the RH part of the TL dominates and the TL has a more pronounced RH behavior (the serial branch features an inductive character, while the shunt branch shows a capacitive one). It is now relevant to adopt physically relevant parameter values for Eq. (5). For applications in the microwave frequency range (e.g., for microstrip lines [2] or coplanar waveguide structures loaded with split-ring resonators [3]; see also Ref. [23] for recent work), typical values of the capacitances and inductances involved in the CRLH structure are of the order of pF and nH, respectively. Here we will use the values LR=1nH,CL=0.1 pF, and LL=0.12 nH; thus the frequencies in Eqs. (7) take the values fse =15.92 GHz, fsh = 14.53 GHz, and fRH =5.03 GHz. In contrast, as concerns the parameters involved with the nonlinear capacitor CR,we assume that the pertinent capacitance corresponds to a HBV diode, which is characterized by the following equation [6] [see also [7], where the same form of C(V) is used, but different parameter values]: C(V)=Cj0Ada 1+|V| Vbr −m ,(8) where Cj0=1.53 fF/μm2is the capacitance corresponding to bias voltage V0=0.2V,Ada =650 μm2is the device area, Vbr =12 V is the breakdown potential, and the exponent m= 2.7 results from fitting experimental data. It is clear that, for sufficiently small V, by Taylor expanding Eq. (8) one obtains Eq. (4), where the constant parameter values involved are CR0=1pF,C R0=−0.24 pF/V, and C R0=−0.08 pF/V2. To this end, the values of the normalized parameters δ,β, and μappearing in Eq. (5) take the following values: δ≈1.1,β≈0.35,μ≈−0.9.(9) Below we will use these values for the purposes of our analytical and numerical considerations (we have checked that other values lead to qualitatively similar results). Notice that our choice leads to δ>1, i.e., we consider the case where the TL has a more pronounced LH character; however, when considering the linear setting (see next section), this parameter will also assume other values, corresponding to the balanced and RH-dominated behavior as well. B. Linear analysis We now assume plane wave solutions of Eq. (5) of the form Vn=V0exp[i(kn −ωt)], where kand ωdenote the wave number and angular frequency, respectively, and the amplitude of the wave is V01. Substituting the above ansatz into Eq. (5) and keeping only the linear terms in V0, we obtain the following linear dispersion relation: ω4−1+δ2+4β2sin2k 2ω2+δ2=0.(10) The above result is illustrated in Fig. 2, where we plot the frequency f/fsh as a function of the wave number k(in rad/cell) for three different values of δ(note that here we consider one period of k, i.e., −π⩽kj⩽π). It is clear that for δ= 1.0954 [Fig. 2(a)] there exist two frequency bands where em wave propagation is possible: the RH band [high-frequency band depicted by the dashed (blue) line] for 1.0954 <f < 1.4535 and the LH band [low-frequency band depicted by the solid (red) line] for 0.7538 <f <1. In the same case Wave number k (rad/cell) f/fsh (a) −3 −2 −1 0 1 2 3 0.8 1 1.2 1.4 RH−RH+ LH−LH+ GAP δ=1.0954 −3 −2 −1 0 1 2 3 0.8 1 1.2 1.4 Wave number k (rad/cell) f/fsh (b) RH+ RH− LH−LH+ δ=1 −3 −2 −1 0 1 2 3 0.6 0.8 1 1.2 Wave number k (rad/cell) f/fsh (c) RH−RH+ LH−LH+ GAP δ=0.7746 FIG. 2. (Color online) Dispersion relation showing the normalized frequency f/fsh as a function of the wave number k(in rad/cell) for different values of δ:(a)δ=1.0954, (b) δ=1, and (c) δ=0.7746. The solid (red) and dashed (blue) lines show the dispersion relation in the LH low- and RH high-frequency regions, respectively; RH±and LH±denote branches with k>0ork<0. If δ= 1 a gap is formed; the width of the gap is |δ−1|for (a) δ>1or (c) δ<1. 013203-3 VELDES, CUEVAS, KEVREKIDIS, AND FRANTZESKAKIS PHYSICAL REVIEW E 88, 013203 (2013) (δ=1.0954), there exists a gap for 1 <f/f sh <δwhere em wave propagation is not possible. In the case where δ=1 (corresponding, e.g., to the value CL=0.12 pF) the gap vanishes [see Fig. 2(b)] and the TL is balanced. In the balanced case, em wave propagation is possible in two frequency bands as well: the RH band [high-frequency band represented by the dashed (blue) line] with 1 <f/f sh <1.405 and the LH band [low-frequency band represented by the solid (red) line] with 0.7117 <f/f sh <1. Finally, for δ=0.7746 (corresponding, e.g., to CL= 0.2 pF), a gap appears again for δ<f/f sh <1 [Fig. 2(c)]. In this case too there exist a RH high-frequency band and a LH low-frequency band for 0.588 <f/f sh <0.7746 and 1<f/f sh <1.317, respectively. Note that in all cases, the RH±and LH±branches correspond to positive or negative k, respectively. Thus, generally, in the linear setting and for a given frequency the em waves may propagate in either the RH region (forward wave propagation) or the LH region (backward wave propagation). However, in the nonlinear setting, coupling between modes propagating in the LH and RH regimes is possible (see, e.g., relevant earlier work in Refs. [17,18]). Below we will demonstrate that this is the case indeed and study the coupling (interaction) between LH and RH modes with equal group velocities. Since the latter are tangents in the dispersion curves, inspection of Fig. 2shows that it is possible to identify domains, belonging to the RH±and LH∓branches, exhibiting parallel tangents, i.e., equal group velocities. To further elaborate on this, we may use Eq. (10) to obtain the group velocity vg≡∂ω/∂k: vg=ω3β2sin k ω4−δ2.(11) In Fig. 3we show the dependence of the group velocity vg on the normalized frequency f/fsh for the values of δused in Fig. 2. Notice that the figure depicts only the group-velocity branches with vg>0 [see solid (red) and dashed (blue) lines] corresponding, respectively, to the LH−and RH+branches of the dispersion curves; the branches with vg<0 (pertinent to the LH+and RH−branches of the dispersion curve) are mirror symmetric with respect to the ones shown in the figure, due to the parity of the dispersion relation. Considering a horizontal cut of the group-velocity curves, say, at vg=0.1or0.075 [see horizontal lines in Figs. 3(a) and 3(c)], it is readily observed that, indeed, a LH−and a RH+ mode can share a common group velocity (and interact in the nonlinear regime, as mentioned above). In fact, inspection of the group-velocity curves, say, in Fig. 3(a), shows that the maximum possible common vgis given by vgmax =0.1339, the local maximum of vg, occurring at f=0.9391, in the (shorter in height) LH low-frequency band. Then one can divide each of the LH and RH group-velocity curves into two subregions, depending on the sign of the group-velocity dispersion (GVD) ∂vg/∂ω, where such coupling with equal group velocities may occur. These subregions are the subbands I (0.7538 < f/fsh <0.9391) and II (0.9391 <f/f sh <1) for the LH lowfrequency band, characterized by positive and negative GVD, respectively, and the subbands III (1.0954 <f/f sh <1.1195) and IV (1.356 <f/f sh <1.4535) for the RH high-frequency band, again characterized by positive and negative GVD, f/fsh Group Velocity υg (a) 0.8 1 1.2 1.4 0 0.05 0.1 0.15 0.2 III III IV GAP δ=1.0954 0.8 1 1.2 1.4 0 0.05 0.1 0.15 0.2 f/fsh Group Velocity υg (b) LH RH δ=1 (c) f/fsh Group Velocity υg 0.6 0.8 1 1.2 0 0.05 0.1 0.15 0.2 III II GAP IIV δ=0.7746 FIG. 3. (Color online) Group velocity vgas a function of the normalized frequency f/fsh (for the values of δused in Fig. 2). The solid (red) and dashed (blue) lines indicate branches corresponding to the LH−and RH+regimes, respectively. The intersection of the group velocity curves with the horizontal solid (black) line depicts frequencies of modes with the same group velocity vg=0.1. Regions I–IV indicate possible interactions between LH−and RH+modes with the same vgbut different signs of GVD. respectively. Thus nonlinear LH and RH modes of equal vg can feature the following four different possible interactions: (i) the LH mode in band II and RH mode in band IV, both featuring negative GVD; (ii) the LH mode in band I and RH mode in band IV, where the LH (RH) mode features positive (negative) GVD; (iii) the LH mode in band I and RH mode in band III, both featuring positive GVD; and (iv) the LH mode in band II and RH mode in band III, where the LH (RH) mode features negative (positive) GVD. It is clear that the above set of possibilities arises from the existence of the gap in the considered case with δ=1.0954. A similar situation also occurs for δ<1, e.g., for δ=0.7746, as in Figs. 2(c) and 3(c). In contrast, for δ=1 the gap no longer exists and thus the only possible interaction is between a LH mode with positive GVD and a RH mode with negative GVD; this interaction can occur for group velocities vg⩽0.175, i.e., beneath the dashed horizontal line in Fig. 3(b). This possibility, however, is already taken into consideration [see case (ii) above]; furthermore, soliton formation in the balanced CRLHTL model (δ=1) has already been studied in the literature 013203-4 COUPLED BACKWARD- AND FORWARD-PROPAGATING ... PHYSICAL REVIEW E 88, 013203 (2013) [11]. For these reasons, below we will proceed by analyzing the case corresponding to δ=1.0954, which offers all possible scenarios. It is clear that the case of δ=0.7746 shares similar qualitative features; this similarity extends beyond the linear wave case and into the nonlinear solitonic one. Although, as explained above, we are not going to analyze soliton formation and solitons in the special case of the balanced CRLH-TL model with δ=1, it is worth mentioning that in the case of δ=1, the dispersion relation exhibits a Dirac point, namely, it is approximately linear in the vicinity of k=0, i.e., ω≈±[1 +(β/2)k] [see Fig. 2(b)]. The emergence of Dirac points is particularly interesting in the two-dimensional (2D) setting of triangular and hexagonal lattices arising in different contexts, such as optics [25], atomic Bose-Einstein condensates [26], and the so-called photonic graphene [27]. This has also led to an interest in this subject from a rigorous mathematical perspective [28]. It is thus quite intriguing that, in principle, 2D balanced CRLH-TL models may host a variety of fundamental effects, such as conical diffraction, formation of topological defects, and even phase transitions, as in Refs. [25–27]. C. Nonlinear analysis: Coupled NLS equations and solitons To describe the coupling between a RH and a LH nonlinear mode with equal group velocities, we will use the quasidiscrete approximation, which takes into account the inherent discreteness of the system (see, e.g., Ref. [19]fora review and Refs. [15,23] for relevant recent work). Generally, this approach allows for the description of quasidiscrete envelope solitons (usually satisfying an effective NLS model), characterized by a discrete carrier and a slowly varying continuum envelope. In our case, since we are interested in the description of two different modes, we seek a solution of Eq. (5) in the form Vn= 2  j=1 Vj(X,T )exp(iθj)+c.c., (12) where c.c. denotes complex conjugate. In Eq. (12), subscripts j=1,2 correspond to the LH and RH modes, Vj(X,T ) are unknown (continuous) slowly varying envelope functions depending on the slow scales X=(n−vgt) (where vgis the common group velocity) and T=2t, and exp(iθj), with θj=kjn−ωjt, are the (discrete) carriers of frequencies ωj and wave numbers kj. Finally, is a formal small parameter setting the field amplitude and the slow scales of the envelope functions. At this point, we should note that the field Vnas expressed in Eq. (12) is in fact the leading-order form of a more general ansatz employing multiple time and space scales. In this context, use of a formal multiscale expansion method leads to a hierarchy of equations at various powers of , which are solved up to third order. We will present the main results here and provide further details in Appendix A. In particular, from the first- and second-order problems [i.e., at O() (linear limit) and O(2), respectively] we derive the dispersion relation (10) and the group velocity (11). Finally, at the next order, O(3), we obtain the following coupled NLS equations: i∂TV1+1 2D1∂2 XV1+(g11|V1|2+g12|V2|2)V1=0,(13) f/fsh g11, g22 (a) 0.8 1 1.2 1.4 −2.5 −2 −1.5 −1 −0.5 0 II GAP IV g11 g22 ** f/fsh g11, g22 (b) 0.8 1 1.2 1.4 −2.5 −2 −1.5 −1 −0.5 0 IGAP IV g11 g22 ** f/fsh g11, g22 (c) 0.8 1 1.2 1.4 −2.5 −2 −1.5 −1 −0.5 0 IGAP III g11 g22 * * f/fsh g11, g22 (d) 0.8 1 1.2 1.4 −2.5 −2 −1.5 −1 −0.5 0 II GAP III g11 g22 * * FIG. 4. (Color online) Nonlinearity coefficients g11 [solid (red) lines] and g22 [dashed (blue) lines] as functions of f/fsh for interactions in bands (a) II and IV, (b) I and IV, (c) I and III, and (d) II and III. The parameter δtakes the value δ=1.0954. Stars (in black) in each panel depict parameter values used in relevant simulations; see Figs. 6and 7[corresponding to (a)], Figs. 9and 10 [corresponding to (b)], Figs. 13 and 14 [corresponding to (c)], and Figs. 15 and 16 [corresponding to (d)] below. i∂TV2+1 2D2∂2 XV2+(g21|V1|2+g22|V2|2)V2=0,(14) where the normalized GVD coefficients Dj, the self-phase modulation (SPM) coefficients gjj , and the cross-phase modulation (CPM) coefficients gj,3−j(with j=1,2) are respectively given by Dj≡∂2ωj ∂k2 j=vgcot kj−ω4 j+3δ2 ωjω4 j−δ2vg,(15) gjj =ω3 jω2 j−δ2 2ω4 j−δ2(3μ−Aj),(16) gj,3−j=ω3 jω2 j−δ2 2ω4 j−δ2(6μ−B3−j),(17) where the coefficients Ajand B3−jare defined in Appendix A. It is now useful to remark that for all four possible wave interactions mentioned in the preceding section, the SPM coefficients gjj are negative (see Fig. 4). Then, measuring normalized time Tand densities |Vj|2in units of |D1|−1and |D1/gjj |, respectively, we express Eqs. (13) and (14) in the following form: i∂TV1+s 2∂2 XV1+(λ1|V2|2−|V1|2)V1=0,(18) i∂TV2+d 2∂2 XV2+(λ2|V1|2−|V2|2)V2=0,(19) where s=sgn(D1),d=D2 |D1|,λ 1=g12 |g22|,λ 2=g21 |g11|.(20) As seen from Eqs. (18) and (19), in the absence of CPM coupling (λj=0) the evolution of either the LH mode V1 or the RH mode V2is described by a single NLS equation. The latter supports soliton solutions of the bright or the dark 013203-5 VELDES, CUEVAS, KEVREKIDIS, AND FRANTZESKAKIS PHYSICAL REVIEW E 88, 013203 (2013) type, depending on the sign of dispersion coefficient (see, e.g., Ref. [20]): In particular, the mode V1(V2) supports bright solitons for s<0(d<0) or dark solitons for s>0(d> 0). These conditions, however, are modified for λj= 0 and various types of coupled (vector) solitons can be found in the full version of Eqs. (18) and (19). Below we will present these types of coupled backward- and forward-propagating solitons belonging, respectively, to the LH low- and RH high-frequency bands. In principle, four types of vector solitons are possible: bright-bright solitons in the form V1(X,T )=a1sech(bX)exp(−iν1T),(21) V2(X,T )=a2sech(bX)exp(−iν2T),(22) BD solitons in the form V1(X,T )=a1sech(bX)exp(−iν1T),(23) V2(X,T )=a2tanh(bX)exp(−iν2T),(24) DB solitons in the form V1(X,T )=a1tanh(bX)exp(−iν1T),(25) V2(X,T )=a2sech(bX)exp(−iν2T),(26) and DD solitons in the form V1(X,T )=a1tanh(bX)exp(−iν1T),(27) V2(X,T )=a2tanh(bX)exp(−iν2T).(28) In the above equations, a1,2and ν1,2denote the amplitudes and frequencies of each soliton and bis the (common) inverse width of the solitons. Now, each of the above ans¨ atze is substituted into Eqs. (18) and (19), leading to a set of equations connecting the soliton parameters. In particular, the equations connecting the amplitudes a1,2with the inverse width bare of the form (aj/b)2=−αj,(29) (aj/b)2=(−1)jαj,(30) (aj/b)2=(−1)3−jαj,(31) (aj/b)2=αj(32) for the BB, BD, DB, and DD solitons, respectively, where parameters αj(j=1,2) are given by α1=dλ1+s 1−λ1λ2 ,α 2=sλ2+d 1−λ1λ2 .(33) In contrast, the frequencies of the BB, BD, DB, and DD solitons are respectively given by ν1=−s 2b2,ν 2=−d 2b2,(34) ν1=−s 2b2−λ1a2 1,ν 2=a2 2,(35) ν1=a2 1,ν 2=−d 2b2−λ2a2 2,(36) ν1=a2 1−λ1a2 2,ν 2=a2 2−λ2a2 1.(37) 0.94 0.96 0.98 1 −4.0 −2.0 0 2.0 4 .0 f1/fsh α1,α2 (a) α1α2 BD BB DB * II and IV 0.94 0.96 0.98 1 −2 −1.5 −1 −0.5 0 f1/fsh λ1,λ2,d (b) λ1λ2d * s=−1 FIG. 5. (Color online) Parameters for soliton interactions in bands II and IV. (a) Dependence of the coefficients α1[solid (blue) line] and α2[dashed (red) line] on the normalized frequency f1/fsh. (b) Dependence of the parameters λ1[thin solid (blue) line], λ2 [dashed (blue) line], and d[bold solid (red) line] on the normalized frequency f1/fsh. It is clear that each type of vector soliton is characterized by five independent parameters, connected by a set of four equations; thus each of the above vector solitons is characterized by one free parameter. Furthermore, the fact that the parameters of each soliton component depends on the ones of the other component clearly highlights an important possibility arising from the nonlinear coupling of the solitonic modes: One can control, e.g., the characteristics (amplitude, width, etc.) of a LH soliton by means of the parameters of its RH pair and vice versa. Requiring that the right-hand sides of Eqs. (29)–(32) are positive, it can be concluded that the existence of each type of vector soliton is determined by the signs of parameters αj: in particular, BB solitons exist for α1<0 and α2<0, BD solitons exist for α1<0 and α2>0, DB solitons exist for α1>0 and α2<0, and DD solitons exist for α1>0 and α2>0. In contrast, we should mention that the signs of αj depend on frequency, through the frequency dependence of the parameters involved in Eqs. (33). Thus, for each different type of mode interaction [see cases (i)–(iv) in Sec. II B], below we will present results for the sign (and magnitude) of αj in bands II and IV, I and IV, I and III, and II and III, and accordingly determine which type of soliton is possible. This way we will show that BB, BD, and DB solitons can exist in certain frequency bands, while DD solitons do not exist: this is due to the fact that α1and α2are either of opposite signs or both negative, as attested to by Figs. 5(a),8,11, and 12. Before proceeding with the presentation of the coupled soliton solutions we make the following comments. First, vector solitons will be found in a stationary form; however, using these stationary solutions, one can also find traveling 013203-6 COUPLED BACKWARD- AND FORWARD-PROPAGATING ... PHYSICAL REVIEW E 88, 013203 (2013) soliton solutions, with an additional free parameter, i.e., the velocity C, by means of the following Galilean boost: V1(X,T )→V1(X−CT,T) ×exp i sCX +C2 2T,(38) V2(X,T )→V2(X−CT,T) ×exp i dCX +C2 2T.(39) Second, it is interesting to note that, contrary to what is often the case in the mathematically studied multicomponent variants of the NLS equation [21], the model of Eqs. (18) and (19) does not necessarily respect the condition λ1=λ2. The latter condition ensures the existence of an underlying Hamiltonian structure and is customary in other physical applications (such as atomic physics [29]). Nevertheless, as we will see below, this is not a necessary condition for the existence of the exact soliton solutions considered below. III. SOLITON INTERACTIONS IN DIFFERENT FREQUENCY BANDS: NUMERICAL RESULTS A. Numerical procedure Let us now proceed to study numerically the evolution of the coupled solitons presented in the previous section in the framework of the fully discrete model of Eq. (5). In order to compare the analytical approximations with the results of numerical simulations, we will make use of two diagnostic quantities: The first one is the evolution of the center of mass defined as X(t)=n=N n=−NnV 2 n n=N n=−NV2 n (40) and the second one is a powerlike quantity defined as P(t)= n=N  n=−N V2 n,(41) with 2N+1 being the lattice size. The above quantities can readily be determined for each type of vector soliton that is predicted analytically in the framework of the coupled NLS equations. In all simulations, which have been performed by means of a fixed-step fourth-order Runge-Kutta scheme with a time step equal to 0.01, we have fixed the value of the small parameter as =0.02 and we have used periodic boundary conditions. Use of the latter leads to the requirement that the wave number kof a dark soliton component must be equal to 2πq/p, with q,p ∈Zand qalso being odd. In all figures below (Figs. 6–16), unless stated otherwise, we show the density plots of Vn, the spatial profile of Vnat t=2000, the time evolution of the center of mass X(t) and the quantity P(t). Regarding the evolution time of the simulations, we should note the following. Most of our simulations are performed for relatively large normalized times, typically up to t∼107 in some cases. However, given our time normalization, the physical unit time (set by the frequency fsh =14.529 GHz) n t (a) −400 −200 0200 400 0 500 1000 1500 2000 −0.01 −0.005 0 0.005 0.01 200 300 400 −0.01 0 0.01 0.02 n Vn (b) Numerical Analytical 0 500 1000 1500 −100 0 100 200 300 t X(t) (c) Numerical Analytical 0 500 1000 1500 2.4 2.42 2.44 2.46 2.48x 10−3 t P(t) (d) Numerical Analytical FIG. 6. (Color online) Bright-bright solitons in regions II–IV. (a) Density plot of the space-time evolution of Vnobtained numerically. (b) Comparison of the analytical and numerical profiles of Vnat t=2000. Also shown is the time evolution of (c) the center of mass and (d) the power diagnostic. The parameters used are f1/fsh =0.965 45 and f2/fsh =1.365 35, which give k1=−0.4061 and k2=1.8576, i.e., a bright-bright soliton in bands II and IV (this particular choice corresponds to the points depicted by stars in Fig. 5). The difference in the powers can be attributed to the approximate nature of our solution. is very small, namely, t0=(2πfsh)−1≈11 ps (see Sec. II A). Actually, since all characteristic frequencies of the system [see Eq. (7)] are in the microwave regime, all characteristic times are less than a nanosecond and thus, obviously, simulations for time teven of the order of 109are extremely time consuming. Nevertheless, our results for normalized times up to t=107 n t (a) −500 0 500 0 5 10x 104 −0.01 −0.005 0 0.005 0.01 −500 0 500 −0.01 0 0.01 t=0 (b) −500 0 500 −0.01 0 0.01 t=5×106 Vn −500 0 500 −0.01 0 0.01 t=107 n 0 5 10 15 x 104 −500 0 500 t X(t) (c) 0 5 10 15 x 104 2.4 2.42 2.44 2.46 2.48 x 10−3 t P(t) (d) Numerical Analytical FIG. 7. (Color online) The bright-bright soliton of Fig. 6evolves until t=105. All the panels are similar to those of Fig. 6except for (b), in which snapshots of the soliton at t=5×106and t=107are compared to the initial condition of the simulation in order to examine its robustness under a very long evolution time. 013203-7 VELDES, CUEVAS, KEVREKIDIS, AND FRANTZESKAKIS PHYSICAL REVIEW E 88, 013203 (2013) 0.8 0.85 0.9 −10 −5 0 5 10 f1/fsh α1,α2 α1α2 BD * * I and IV FIG. 8. (Color online) Parameters for solitons in bands I and IV. Shown is the dependence of the coefficients α1[solid (blue) line] and α2[dashed (red) line] on the normalized frequency f1/fsh.Stars depict parameter values used for the simulations shown in Figs. 9and 10 below. (corresponding to a physical time of the order of one-tenth of a millisecond) demonstrate good agreement with our analytical predictions in suitable cases (see below). Furthermore, the results of such long simulations can also be used as a reliable indication of the solitons’ robustness. Hence, in the case where the solitary waves are found to be very robust, we expect that they would survive for the longer time scales that would render them experimentally observable. B. BD, BB, and DB solitons in bands II and IV First, we consider the interaction between a backwardpropagating soliton, with a frequency lying in band II, and a forward-propagating soliton, with a frequency lying in band IV (for δ=1.0954). In this case, s=−1 (see Fig. 3) and, following the analysis of the previous section, we find three different types of vector solitons in certain subbands. n t (a) −500 0 500 0 500 1000 1500 2000 −5 0 5 x 10−3 0 200 400 600 −0.01 −0.005 0 0.005 0.01 n Vn (b) Numerical Analytical 0 500 1000 1500 −200 0 200 400 t Xd(t) (c) Numerical Analytical 0 500 1000 1500 0.048 0.0485 0.049 0.0495 0.05 t P(t) (d) Numerical Analytical FIG. 9. (Color online) Bright-dark solitons in regions I and IV. (a) Density plot of the time evolution of Vnobtained numerically. (b) Comparison of the analytical and numerical profiles of Vnat t= 2000. Also shown is the time evolution of (c) the center of mass and (d) the power diagnostics. The parameters used are f1/fsh =0.8831 and k2=5π/8≈1.9625, which give k1=−1.0404 and f2/fsh =1.3748 (see the corresponding points depicted by stars in Fig. 8), i.e., a bright-dark soliton in bands I and IV. n t (a) −1000 0 1000 0 5 10 15 x 104 −5 0 5 x 10−3 −1000 −500 0 500 1000 −0.01 0 0.01 t=0 (b) −1000 −500 0 500 1000 −0.01 0 0.01 t=2×105 Vn −1000 −500 0 500 1000 −0.01 0 0.01 t=5×105 n 0 5 10 15 x 104 −4000 −2000 0 2000 4000 t Xd(t) (c) 0 5 10 15 x 104 0.048 0.0485 0.049 0.0495 0.05 t P(t) (d) Numerical Analytical FIG. 10. (Color online) The bright-dark soliton of Fig. 9evolves until t=2×105. All the panels are similar to those of Fig. 9except for (b), in which snapshots of the soliton at t=2×105and t= 5×105are compared to the initial condition of the simulation. Notice that the center of mass is not bounded in [−N,N]; this can be caused by the soliton splitting. In this setting, the modulational instability of the background (see the discussion in Appendix B) appears to be responsible for the breakup of the wave packet. (i) If 0.94 <f 1/fsh <0.96 then α1<0 and α2>0; thus BD solitons exist, with the frequency of the bright (dark) soliton component being in the LH low- (RH high-) frequency band. (ii) If 0.96 <f 1/fsh <0.98 then α1<0 and α2<0; thus BB solitons exist, with the soliton frequencies being one in the LH low- and one in the RH high-frequency band. (iii) If 0.98 <f 1/fsh <1 then α1>0 and α2<0; thus DB solitons exist, with the frequency of the dark (bright) soliton being in the LH low- (RH high-) frequency band. Here we will focus on the BB soliton case (BD and DB solitons will be studied below). In particular, we will investigate the dynamics of a BB soliton corresponding to the frequency value f1/fsh =0.965: For this value, the dispersion coefficient dbecomes d≃−1=s, the CPM coefficients λ1,2 take the values λ1=λ2=λ=−1.7, and α1=α2=0.487 (see the intersection point of the relevant curves depicted by a star in Fig. 5). In this case, the BB soliton has a (common 0.8 0.85 0.9 −100 −50 0 50 100 f1/fsh α1,α2 α1α2 DB * * I and III FIG. 11. (Color online) Same as Fig. 8, but for soliton interactions in bands I and III. Stars depict parameter values used for the simulations shown in Figs. 13 and 14 below. 013203-8 COUPLED BACKWARD- AND FORWARD-PROPAGATING ... PHYSICAL REVIEW E 88, 013203 (2013) 0.94 0.96 0.98 1 −40 −20 0 20 40 f1/fsh α1,α2 α1α2 DB * * II and III FIG. 12. (Color online) Same as Fig. 8, but for soliton interactions in bands II and III. Stars depict parameter values used for the simulations shown in Figs. 15 and 16 below. for both components) group velocity vg=0.1288, which occurs when the (normalized) carrier frequencies for the modes V1and V2take, respectively, the values f1/fsh =0.965 (as mentioned above) and f2/fsh =1.365. Notice that other BB soliton solutions not obeying the above particular conditions (λ1=λ2and d=−1) exist as well, as per our analysis of the previous section. At this point, we should mention that this choice of the parameter values leads to symmetric coupled NLS equations (18) and (19) (and also to symmetric BB solitons of equal amplitudes; see, e.g., Ref. [30]). A very interesting subcase, corresponding to λ=−1 (i.e., equal SPM and CPM coefficients), would render Eqs. (18) and (19) the completely integrable Manakov system [31]; in such a case, in Fig. 5(b), the curve for parameter dwould intersect the curves for λ1,2, with all parameters taking the value −1. However, we have checked that other, physically relevant, choices of the values of the parameter δdo not lead to the Manakov case [the chosen n t (a) −500 0 500 0 500 1000 1500 2000 −0.02 −0.01 0 0.01 0.02 −500 0500 1000 −0.02 0 0.02 n Vn (b) Numerical Analytical 0500 1000 1500 −100 0 100 200 300 t Xd(t) (c) Numerical Analytical 0 500 1000 1500 0.25 0.26 0.27 0.28 0.29 t P(t) (d) Numerical Analytical FIG. 13. (Color online) Dark-bright solitons in regions I and III. (a) Density plot of the time evolution of Vnobtained numerically. (b) Comparison of the analytical and numerical profiles of Vnat t=2000. Also shown is the time evolution of (c) the center of mass and (d) the width diagnostic. The parameters used are k1=−6π/5≈ −1.884 and f2/fsh =1.1002, which give f1/fsh =0.8003 and k2= 0.1232 (cf. points depicted by stars in Fig. 11), i.e., a dark-bright soliton in bands I and III. n t (a) −2000 0 2000 0 5 10 15 x 104 −0.02 −0.01 0 0.01 0.02 −3000 −2000 −1000 01000 2000 3000 −0.02 0 0.02 t=0 (b) −3000 −2000 −1000 01000 2000 3000 −0.02 0 0.02 t=2×105 Vn −3000 −2000 −1000 01000 2000 3000 −0.02 0 0.02 t=5×105 n 0 5 10 15 x 104 −4000 −2000 0 2000 4000 t Xd(t) (c) 0 5 10 15 x 104 0.71 0.72 0.73 0.74 0.75 t P(t) (d) Numerical Analytical FIG. 14. (Color online) The dark-bright soliton of Fig. 13 evolves until t=2×105. All the panels are similar to those of Fig. 13 except for (b), in which snapshots of the soliton at t=2×105and t= 5×105are compared to the initial condition of the simulation in order to examine its robustness under a very lengthy time evolution. value δ=1.0954 only leads to symmetric NLS equations (18) and (19)]. Utilizing the above-mentioned BB soliton solution, we can express the voltage Vn(t)inEq.(5) in terms of the original coordinates nand tas follows: Vn(t)≈V0[R1(n,t) cos(k1n−1t) +R2(n,t) cos(k2n−2t)],(42) n t (a) −500 0 500 0 500 1000 1500 2000 −0.02 −0.01 0 0.01 0.02 −500 0 500 1000 −0.02 0 0.02 n Vn (b) Numerical Analytical 0500 1000 1500 −100 0 100 200 t Xd(t) (c) Numerical Analytical 0500 1000 1500 0.71 0.72 0.73 0.74 t P(t) (d) Numerical Analytical FIG. 15. (Color online) Dark-bright solitons in regions II and III. (a) Density plot of the time evolution of Vnobtained numerically. (b) Comparison of the analytical and numerical profiles of Vnat t= 2000. Also shown is the time evolution of (c) the center of mass and (d) the width diagnostic. The parameters used are k1=−3π/23 ≈ −0.4095 and f2/fsh =1.1162, which give f1/fsh =0.965 and k2= 0.2758 (cf. stars in Fig. 12), i.e., a dark-bright soliton in bands II and III. 013203-9