Backward- and forward-wave soliton coexistence due to second-neighbor coupling in a left-handed transmission line

We study analytically and numerically the impact of the second-neighbor interactions on the propagation of an initial wave-packet through a coupled nonlinear left-handed transmission line. We start with a detailed analysis of the dispersion curve and group velocity. We then use the quasi-discrete approximation with respect to long wavelength transverse perturbations to reduce the nonlinear discrete model of the lattice to a nonlinear Schrödinger equation. A detailed analysis of the product of the group velocity dispersion (P) and the nonlinear term (Q) is presented with respect to the relevant parameters of the system. Here we see the important role of the second-neighbor coupling, which makes the dispersion curve non-monotonic. Thus, for a given frequency, the group velocity can be both positive and negative, while for each choice the product, PQ, can be either positive or negative. Our numerical results confirm that for sufficiently strong second-neighbor coupling, bright backward-wave pulses and dark solitons can coexist at a common frequency, and these travel in opposite directions through the lattice.


Introduction
Metamaterials (MMs) represent an active area of research that continues to attract substantial attention in the scientific community. This is because of their novel properties that could bring technological and scientific advancements [1]. Often called left-handed materials (LHMs) whose effective permittivity, , and magnetic permeability, μ, are simultaneously negative, MMs possess many extraordinary qualities, such as negative refraction in the microwave spectrum, backwardwave propagation and backward Cherenkov radiation [2,3]. In these materials with negative index, energy and wavefronts would travel in opposite directions [4][5][6]. The applications of such materials include superlenses and cloaking [7]. One requirement is that electromagnetic waves propagate with negative group velocity, and this leads us to left-handed transmission lines (LH TL).
Today the LH TL is well known to provide an alternative representation of LHMs. Since its first realization by [8,9], a diverse set of works has explored various aspects from a physics as well as a mathematical outlook. The advantage of using LH TL is that it provides a useful way to study the wave properties of the LHMs, to check the role of nonlinearity, and even to model the exotic properties of new systems. Nonlinearity may arise from embedding an array of wires and split ring resonators into a nonlinear dielectric or by inserting diodes into resonant conductive elements in LH TL [1,[10][11][12], and this allows for soliton formation and pulse propagation in such media. Several works have focused on soliton formation in LH TL and showed that the dynamics was governed by the wellknown nonlinear Schrödinger equation (NLSE) [13][14][15][16]; these studies also looked at modulational instability (MI) and found exact analytical solutions. Other studies elucidated interesting consequences of secondneighbor and non-local interactions on nonlinear excitations in discrete systems more broadly [17][18][19][20][21][22][23]. All of these studies highlight that second-neighbor interactions can substantially change the dynamics of the system.
We devote the present study to soliton formation in a LH TL with second-neighbor interactions. In this physical context, we briefly highlight two recent references. The first [24] examined the highly nonlinear regime, where the discreteness of the lattice is essential; the second [25] investigated the more weakly nonlinear limit, where a semi-discrete approximation can be used. The former found resonant ILMs (in addition to the more standard ILMs). The latter study explored solitons (which are less sharply localized than ILMs) and found the simultaneous presence of two types of solitonic excitations, namely either two bright solitons or one bright and one dark soliton, coexisting at the same frequency.
The system investigated in this paper centrally differs from those two studies in that the lattice is intrinsically left-handed before adding second-neighbor interactions. To our knowledge, this combination of nonlocal interactions and left-handedness is novel. Furthermore, in addition to an analytical treatment, we explore the governing circuit equations numerically-a nontrivial endeavor, given the combination of non-locality, left-handedness and nonlinearity.
Our detailed analysis shows that the nature of the lattice excitations can be either backward-or forward-propagating, depending on the strength of the second-neighbor interaction. This indicates that second-neighbor coupling does more than just act as a perturbation for solutions that would exist in its absence; it allows for qualitatively new wave phenomena to emerge. The framework for analysis is the quasi-discrete approximation with respect to long wave-length perturbations, which yields a nonlinear Schrödinger (NLS) equation. A close examination of the P and Q coefficients, as well as their product, allows us to classify the nature of the analytical solutions; such an approach proves fruitful in many physical context (see, for instance, [26]). Intriguingly, for sufficiently strong second-neighbor coupling, the lattice exhibits a non-monotonicity in its dispersion curve, and this opens up regimes where both backward-wave and forward-wave solitonic excitations can propagate at the same frequency. We confirm this prediction numerically. The fact that these two excitations can exist at a common frequency raises the prospect of them being supported by a single driver in real lattices that feature dissipation.
The presentation of our results and the structure of the paper are as follows. In Sect. 2, we describe our model and establish the nonlinear network equations including the second-neighbor coupling. From there, we derive the dispersion relation of plane waves as well as the group velocity. Section 3 outlines the semidiscrete approximation and the perturbation analysis (in the limit of small amplitude) to obtain the NLS equation. We focus on evaluating the product P Q to see how the second-neighbor interactions affect the solutions and their propagation behavior. In Sect. 4, we outline the nontrivial numerical procedure used to integrate the governing equations. We show the numerical results and compare to the analytical predictions.
Finally, in Sect. 5, we summarize our conclusion and present some possibilities for future studies.

Model and equations of motion
2.1 Presentation of the model Figure 1 shows the coupled LH TL under consideration in the present study. We consider a coupled model of the unit cell studied by [13] with second-neighbor interactions. This model consists of N identical cells, and each elementary cell is comprised of a reverse polarized diode, where the capacitance C L (V ) varies with voltage, and an inductor L L = 470 µH in parallel with a capacitor C R = 0.05C L . In order to take into account second-neighbor interactions, we introduce a second linear inductor with inductance L. By applying Kirchhoff's loop rule for two loops involving I 1 (the right one additionally involving I 3 and traversing capacitors C n+1 and C n+2 , and the left involving I 2 and traversing capacitors C n−1 and C n ), we obtain the following equations fairly quickly: where U n = V n−1 −V n is the voltage across the nonlinear capacitance C L (V ), and I L L is the current through inductor L L . Differentiating Eq. (1) with respect to time, while substituting Eq. (2) in the latter, we obtain: where the nonlinear capacitance C L (V ) is chosen as in Ref. [13], C L (U n ) = C 0 − a 1 U n + a 2 U 2 n , and the constant coefficients take the values C 0 = 800 pF, a 1 = 490 pF/V and a 2 = 156 pF/V 2 . This equation can be reduced to the following dimensionless form: where time τ is in units of √ (L L C 0 ) and voltage in units of C 0 a 1 ; δ = L L L , g = C R C 0 and γ = a 2 C 0 (a 1 ) 2 .

Dispersion relation and the effect of the second-neighbor interactions
By considering the solutions of Eq. (4) in the form of plane waves of low amplitudes defined by, where V 0 1, ω and k denote the angular frequency and wave number, respectively, we obtain the following dispersion relation: The dispersion curve as a function of the wavenumber for different values of the coupling inductance L is shown in Fig. 2. Each color of the curve corresponds to a specific value of L. Looking at the plots, we observe that the frequency increases with decreasing values of the coupling inductance. For small values of L, the upper cutoff frequency moves higher, and the dispersion curve loses its monotonicity, where it first increases and then decreases as k is incremented. In the figure, this happens for the coupling inductances of 0.066L L , 0.1L L , and 0.2L L (which translates to δ = 15, 10, 5, respectively). In contrast, for the large values of the coupling inductance the cutoff frequency does not change, and only the intermediary frequencies rise-see the green and red curves around k ∈ [0.5, 2.5] when L = L L or L = 2L L . Note, for instance, that for k = 1 rad/cell the frequency increased significantly from 261 kHz (for the ideal line, δ = 0) to 513 kHz when δ = 1.
From these curves, we see that below a certain wavenumber, and for sufficiently low L (or high δ), the corresponding lattice would be expected to behave like a right-handed transmission line for long wavelengths. We can use Eq. (6) to derive the group velocity, v g , as a function of wavenumber: Equation (7) is plotted in Fig. 3 to visualize the dependance of the group velocity on wavenumber. We see that v g is always negative for sufficiently large coupling inductance (L = +∞, L L , 2L L , which correspond to the solid blue, green and red curves, respectively); for these traces, v g < 0 for all values of wavenumbers and frequencies. For smaller coupling inductances, v g can be both positive and negative (i.e., for 0.066L L , 0.1L L , 0.2L L ). This confirms the loss of monotonicity of the dispersion curve.
Another relevant feature worth mentioning is that we get two values of k corresponding to a unique value of frequency in the region above the upper cutoff frequency for the ideal lattice (without second-neighbor coupling). This means that two different waves can simultaneously propagate at the same frequency. The influence of the second-neighbor interactions on the dispersion curve and group velocity illustrates that L substantially affects the propagation of plane-waves, and it hints at its role in soliton formation as well. Let us now investigate this relationship analytically to explore the dynamics of the system in more detail.

The nonlinear Schrödinger equation
In order to theoretically understand the dynamics of the system and to showcase the effects of the secondneighbor interactions, we employ the semi-discrete approximation and the reductive perturbation method in the quasi-discrete limit [27]. Here we consider waves with a slowly varying envelope in time and space with respect to a given carrier with angular frequency ω and wavenumber k. Thus, let us look for solutions to Eq. (4) of the form: where V l (X, T ) is the envelope of the wave; exp(il(ωτ − kn)) represents its carrier; l is an integer that designates the order of magnitude of the perturbation; ω and k denote the angular frequency and the wavenumber, respectively; and c.c. represents the conjugate complex. The dominant motion is in the n direction, and we define two stretched variables in the wave frame: Here V g is the the group velocity, and ε is a small parameter related to the wave amplitude. By substituting Eq. (8) into Eq. (4), we get the resulting equations from the multi-scale expansion, and in this way we arrive at the dependent equations of order ε and ε exp iθ .
In the linear limit (order O(ε)), we obtain the dispersion relation for plane-waves, Eq. (6). The next order ( ε 2 , exp iθ ) returns the group velocity, given by Eq. (7). Finally, to order ( ε 3 , exp iθ ), the well-known NLSE for the unknown voltage V l = V l (X, T ) is obtained: with and Q = −8ω 5 (2 sin(k) − sin(2k)) (cos(k) − 1) sin(k) Here P = 1 2 ∂ V g ∂k and Q are the coefficients of the dispersion and nonlinearity, respectively. Physically, P measures the linear dispersion of the wave and Q determines how the nonlinear effects modulate the frequency. Equation (10) is in the form of the NLSE equation, and this equation is known to admit bright and dark soliton solutions. It has been studied extensively in various fields of physics and has proved useful in characterizing bright and dark solitons observed in many experiments.

Bright and dark soliton-like solutions to the NLSE
In mathematical physics, it has been demonstrated that the NLSE equation has bright and dark soliton-like ana-lytical solutions. These solutions depend in general on the sign of product P Q. Bright soliton solutions are, of course, the well-known pulse excitations whose shape does not deform during the propagation and arises from a stable competition between the nonlinear and dispersive effects in the system [28].
Let us now examine the behavior of the product P Q as a function of frequency. The sign of the product P Q for Eq. (10) depends directly on the wavenumber, and indirectly on the frequency, and it determines the nature of the solutions. The panels of Fig. 4 show the development of the sign of P Q for various couplinginductance values, L. In each panel, we compare the product P Q to what was found in Ref. [13] in the case of L = ∞ (blue curve) where the second neighbor is neglected. This information is also summarized in Table 1.
We now highlight the most salient features contained in Table 1. At larger values of δ, given in the last two rows of the table and shown in panels Fig. 4c and d, there are special zones where we obtain two values of the product P Q at a given frequency. These two values occur at different k-values, but they map to the same frequency. For instance, when L = 0.1L L and δ = 10, the special zone extends from 1100 kHz to 1512 kHz, creating the distinctive loops we see in the figure panels. In fact, we can further sub-divide this special zone into (i) 1100 kHz < f < 1226 kHz, where a unique frequency yields two values of the product P Q: one negative and the other positive; (ii) 1226 kHz < f < 1352 kHz, with two negative values of P Q; and (iii) 1352 kHz < f < 1512 kHz, with two positive values of P Q. In the case of opposite signs for the product (positive and negative values), we expect the degeneration of the initial wave into both bright soliton (for the positive value) and dark soliton (for the negative value). Conversely, in the case where both values are positive, two bright solitons should be able to propagate through the network. Finally in the case of two negative values of P Q, we expect the existence of two different dark solitons.
We can conclude by noting that the emergence of the simultaneous positive and negative value of the product P Q, or two negative, or two positive values for a unique frequency could result in the emergence of a pair of two bright solitons, a bright and a dark soliton, or two dark solitons. This prediction will be interesting to confirm by either numerical simulations or experimental observations.

Numerical analysis approach
In this section, we present the numerical simulations performed on the discrete Eq. (4), governing the propagation of the waves along the discrete LH TL comprised of 500 cells as illustrated in Fig. 1.
The structure of Eq. (4) precludes a direct numerical integration using the more standard techniques, such as Euler or Runge-Kutta; its second time-derivatives are entangled with the discrete spatial Laplacians, while nonlinear terms also feature prominently. This combination prevents even the use of a straightforward implicit algorithm. The numerical scheme we employ on this equation needed to combine implicit and explicit numerical elements, as we now describe in more detail.
The first step involves casting the governing equations, Eq. (4), in the form of a system of first-order differential equations. Here it proves slightly advantageous to move the the variables U n = V n−1 − V n that measure the voltages across the diodes, as introduced earlier, yielding: Now let us introduce W n =U n . While it is not necessary to do so, for simplicity of presentation we now assume that γ is small and drop the term proportional to γ . Using implicit differentiation for the time derivatives and some additional steps, Eq. (13) can be written in the following form: Equation 14 can then be cast in matrix form as: whereẆ is the N-dimensional column-vector comprised of the entries in the set {Ẇ n }, similarly for U, and W 2 represents the column vectors with entries from {W 2 n }. Furthermore, A, B, and C are N × N matrices. Assuming U n 1, the three matrices simplify to: Notice also how periodic boundary conditions were incorporated in these matrices. We can now multiply both sides of Eq. (15) on the left by A −1 , thereby isolating the derivative column-vector on the left side of Eq. (15). Remembering thatU = W, we now have a system of first-order equations that can be numerically integrated using the standard algorithms. For simplicity, here we choose the Euler-Cromer method with a small time-step. Performance could be improved by incorporating more sophisticated numerical recipes with adaptive step-sizes [29,30]. The final considerations pertain to the initial conditions. Since this is a second-order system, we need to specify both U n (0) and W n (0). When we test for bright soliton solutions, we start with initial conditions given by: where P and Q are given by Eqs. (11) and (12), and k represents the wavenumber. If W n (0) is choosen to be zero, we will get two equal-amplitude pulses traveling in opposite directions in the ring. We can also steer the soliton to the right or left by choosing a W n (0) consistent with U n (0) and the dispersion curve, Eq. (6), i.e., W n (0) = A sech A Q 2P n ω k sin(kn). For dark solitons, we chose the following form: For more information, "Appendix" lists the pseudo code used.

Numerical results
We used the software package Igor Pro [31] for performing the simulations. The pseudo-code used is given in "Appendix." Typical run-times for N = 500 are several minutes on a standard PC. When δ is small, the second-neighbor coupling is suppressed and we expect to see backward-wave soliton propagation. This prediction is confirmed numerically in Fig. 5. The soliton is seen to travel to the left, whereas the phase-velocity is to the right, as indicated by the orientation of any particular trough or crest.
We now move to test one of the central findings of this paper, namely the possibility of both forwardand backward-wave solitons existing at the same frequency in this system when the second-neighbor effect is increased. Thus, we increase the δ to 10-the resulting linear dispersion curve is shown in Fig. 6. We will test two wavenumbers that share the same frequency: k 1 = 0.12 and k 2 = 1.4 for which ω = 4.8. Notice that the slopes and thus group-velocities are of opposite signs for these two k-values. Moreover, for this frequency of 4.8, we are in the special zone-see Table 1; one of the wavenumbers should yield a bright soliton and the other a dark soliton. The simulations for k = 1.4 are shown in the density graph of Fig. 7a. We see that this wavenumber produces a backward-wave bright soliton that travels to the left. Figure 7b depicts the soliton profile at two times. The lower (black) trace shows the profile close to the beginning of the simulation at t = 25, and the upper (red) trace gives a snapshot at t = 450. Due to initial conditions, there is a small piece of the launched wavepacket that travels to the right, but even this piece represents a backward wave-packet. Most of the initial wavepacket, however, forms the left-moving soliton, which maintains its shape over time.
Next, let us explore the other wavenumber consistent with the frequency ω = 4.8, namely k = 0.12. Here we would expect a stable dark soliton. The results are shown in Fig. 8. Upon close inspection, it is also evident that the dark soliton travels in the same direction as the phase-velocity.
We should note that launching a bright wave-packet at this wavenumber (k = 0.12) does not seem to result in a stable propagating soliton but instead the wavepacket significantly broadens over time.
Finally, we explore the question raised before about whether these two solitons-one bright and backwardwave, and the other dark and forward-wave-can propagate in the same network. For this purpose, we launch the dark soliton (with k = 0.12) centered at N = 200, while simultaneously launching a bright soliton (with k = 1.4) centered at N = 450. The results are shown in Fig. 9. We see that the two solitons travel through each other in opposite directions. The bright discrete soliton does appear to broaden somewhat after the collision, but this could be due to the numerical approximation where we dropped some nonlinear terms in the matrix A.

Conclusion and outlook
We investigated analytically and numerically the formation of solitons in a nonlinear transmission line with left-handed nearest-neighbor coupling but righthanded second-neighbor coupling. It was found that such a lattice has the possibility of supporting four types of soliton-like solution-bright and dark, left-handed and right-handed-depending on the strength of the second-neighbor coupling. After deriving the nonlinear model of the line, we employed the quasi-discrete approximation with respect to long wave-length transverse perturbations, and we reduced the dynamics of the lattice to a NLS equation which supports bright and dark soliton-like solutions. Then the detailed effect of the second neighbor interactions (represented by a linear inductor) on the product of the NLS coefficients was examined. Our results indicate that the novel secondneighbor coupling inductance can generate more solutions than what had been seen by [13], where the line without this coupling was considered. More interestingly, there is now a frequency band where the product P Q is positive and negative for a unique value of frequency. This raises the possibility that, in this zone, both bright and dark solitons could simultaneously propagate through the network.
To verify these predictions, we devised a numerical scheme to integrate the governing equations-a task complicated by their intricate form that entangle spatial and temporal second derivatives and nonlinearities. Our numerical results proved broadly consistent with the analytical findings. In particular, we were able to verify the coexistence of bright and dark solitons propagating on opposite directions.