A PT-symmetric dual-core system with the sine-Gordon nonlinearity and derivative coupling

We propose a system of sine-Gordon equations, with the $\mathcal{PT}$ symmetry represented by balanced gain and loss in them. The equations are coupled by sine-field terms and first-order derivatives. The sinusoidal coupling stems from local interaction between adjacent particles in the coupled Frenkel-Kontorova (FK) chains and related sine-lattices, while the cross-derivative coupling, which was not considered before, is induced by \emph{three-particle} interactions, provided that the particles in the parallel FK\ chains move in different directions. Nonlinear wave structures are then studied in this model. In particular, kink-kink (KK) and kink-antikink (KA) complexes are explored by means of analytical and numerical methods. It is predicted analytically and confirmed numerically that the complexes are unstable for one sign of the sinusoidal coupling, and stable for another. Stability regions are delineated in the underlying parameter space. Unstable complexes split into free kinks/antikinks that may propagate or become stationary, depending on whether they are subject to gain or loss, respectively.


Introduction
Dual-core waveguides with intrinsic nonlinearity carried by each core offer a convenient setting for the creation of stable dissipative solitons, by application of linear gain to one core, and leaving the parallel-coupled mate one lossy. This possibility was first proposed in the context of nonlinear fiber optics in Ref. [1], see also a review in Ref. [2]. More recently, a similar scheme was elaborated for the application of gain and stabilization of solitons in plasmonics [3], as well as for the creation of stable two-dimensional dissipative solitons and solitary vortices in dual laser cavities [4]. Commonly adopted models of dual-core nonlinear waveguides are based on linearly coupled systems of nonlinear Schrödinger (NLS) equations, which include gain and loss terms [2]. One of the advantages provided by these systems for the theoretical analysis is the availability of exact analytical solutions for stable solitons [5].
A crucial difference between dissipative solitons and their counterparts in conservative media is the fact that the former ones exist as isolated attractors, selected by the balance between gain and loss [6]. On the contrary, nonlinear conservative models, including those originating from optics [7], give rise to continuous families of soliton solutions, rather than isolated ones. and kink-antikink (KA) complexes (in addition to them, nontopological small-amplitude breathers are briefly considered too). Analytical and numerical results for the KK and KA states are reported in Secs. III and IV. The most essential results are existence and stability boundaries for the KK and KA states, delineated in the underlying parameter space by means of analytical and numerical methods. The paper is concluded by Sec. V, suggesting also some potential extensions to future work.

The coupled sine-Gordon system
The P T -symmetric system of coupled SG equations for a real amplified field, φ (x, t), and an attenuated one, ψ (x, t), is adopted in the following form: where coefficient α represents the balanced gain and loss, and is the coefficient of the inter-chain sinusoidal coupling in the double FK chain [23] (a similar coupling appears in a triangular system of three long JJs with a trapped magnetic flux [34]). Such sinusoidal coupling has also been considered extensively in terms of the so-called sine-lattices [36], which represent, e.g., base-rotator models of the DNA double helix [37]. Coefficient β in Eqs. (1) and (2) represents a new type of the anti-symmetric cross-derivative coupling between the two SG equations (note that reflection x → −x makes it possible to fix β > 0). It is different from the usual magnetic coupling between stacked JJs, which would be represented by symmetric second-order cross-derivatives [32]. The derivation of this coupling in terms of coupled FK chains is outlined below.
The Hamiltonian corresponding to the conservative version of Eqs. (1) and (2), with α = 0, is (the last term in the integrand may be replaced by its symmetrized form, − (β/2) (φψ x − ψφ x )). The P T transformation for the system of Eqs. (1) and (2) is defined as follows: It includes the swap of ψ and φ as the P transformation in the direction transverse to x, as in usual P T -symmetric couplers [15,16]. Obviously, the system is invariant with respect to transformation (4). The anti-symmetric cross-derivative coupling emerges in a "triangular" dual FK system schematically displayed in Fig. 1, where a and h are the spacing in each FK chain, and the separation between the parallel chains, respectively. It is assumed that particles with coordinates v n (t) belonging to the bottom chain may move in the horizontal direction, while particles with coordinates u n (t), which belong to the top chain, move along a different direction, under fixed angle θ with respect to the horizontal axis. The inner energies of the two chains (the interaction between them is considered below) can be written as where m is the mass of each particle, κ is the strength of the elastic coupling along each chain (equal coefficients in front of (u n − u n−1 ) 2 and (v n − v n−1 ) 2 is a consequence of identity cos 2 θ + sin 2 θ ≡ 1), W 0 is the depth on the onsite potential, and b is its period. The model may have b a, if the FK chains are built as superlattices on top of an underlying lattice potential; then, kinks, which are considered below, represent a relatively weak deformation of the chains. The equations of motion are generated by the energy as usual, In the continuum approximation, which corresponds to inner energy (5) generates the corresponding terms in Eqs.
Proceeding to the consideration of the coupling between the top and bottom chains, we note that the usual local coupling may be interpreted as produced by energies of diagonal springs linking adjacent particles. These energies are, in turn, proportional to squared lengths of these links. In particular, the sum of the squared lengths for the pair of links connecting the n-th particle in the top chain to its neighbors in the bottom chain, with numbers n and n + 1, is A straightforward consideration of the continuum limit (7) for the corresponding FK Hamiltonian, demonstrates that the last term in expression (8) represents the local-coupling energy, which indeed gives rise to the linearized form [sin (φ − ψ) ≈ φ − ψ] of terms ∼ in Eqs. (1) and (2), so that cos θ ∼ − (the most essential results are obtained below for < 0, which thus corresponds to cos θ > 0). Other terms in Eq. (8) do not represent the inter-chain coupling, and may be absorbed into an appropriate definition of the Hamiltonian of each chain. In particular, term a (v n+1 − v n ) becomes a derivative ∼ ψ x in the continuum limit, hence this term drops out from the continuum Hamiltonian, as mentioned above.
On the other hand, the cross-derivative couplings ∼ β in Eqs. (1) and (2) may be induced by three-particle interactions (TPIs) in the underlying FK system, while usual binary interactions cannot give rise to this coupling (FK models with TPIs were considered in few previous works [38]). In the simplest case, the energy of the relevant TPI can be defined to be proportional to the sum of areas of the respective three-body triangles, as shown in Fig. 1. (i.e., the TPI is induced by the "surface tension" of the triangles). Interactions of this type may be realized in heterogeneous structures with FK chains attached to nanolayers which provide the surface tension, such as graphene (by dint of a The energy of the local two-particle coupling between the chains, which gives rise to terms ∼ in Eqs. (1) and (2), is determined by squared lengths of the red links, see Eq. (8). The total energy of the three-particle coupling, which generates terms ∼ β in Eqs.
(1) and (2), is proportional to the combined area of the shaded triangles, see Eq. (9). The polarization angle θ of the motion in the top chain determines the relative strength of the two couplings in Eqs. (1) and (2) as per Eq. (10). Other details of the setting are explained in the text.
technique similar to that reported in Ref. [39]) or other materials (see, e.g., Ref. [40]). This type of the TPI may also be adopted as a relatively simple phenomenological model. To derive the cross-derivative-coupling term in the Hamiltonian induced by the TPI, we note that the area of the triangle, confined by the links whose lengths are given by Eq. (8), is The transition to the continuum limit as per Eq. (7) transforms the last term in Eq. (9) into the last term in the Hamiltonian density corresponding to Eq. (9), similar to its above-mentioned counterpart a (v n+1 − v n ) in Eq. (8), becomes a full derivative, ψ x , in the continuum limit, hence it may be dropped from the Hamiltonian. Also, term (sin θ) au n may be absorbed into the definition of the inner-chain FK Hamiltonian. Lastly, the above considerations demonstrate that the "polarization angle", θ, in the dual-FK chain (see Fig. 1) determines the relative strength of the two couplings: This relation shows that the cross-derivative coupling, β = 0, emerges when the motion directions in the coupled FK chains are not parallel (θ = 0), while the usual coupling, with = 0, acts unless the two directions are mutually perpendicular, θ = π/2.

Conditions for stability of the flat states
As the main objective of the work is to produce KK and KA complexes, which interpolate between flat states with φ, ψ = 0 (mod 2π), and explore stability of the complexes, a necessary preliminary condition is the stability of the flat states against small perturbations. To address this issue, we use the linearized version of Eqs. (1) and (2), which governs the evolution of small perturbations added to the flat states: The substitution of the usual ansatz for eigenmodes of small perturbations, {φ, ψ} = {φ 0 , ψ 0 } exp (ikx − iωt) , in Eqs. (11) and (12) yields the biquadratic dispersion equation for ω: Straightforward algebraic manipulations demonstrate that Eq. (13) gives rise to a purely real, i.e., instability-free spectrum, with ω 2 (k 2 ) ≥ 0 for all k 2 ≥ 0, under the following conditions: Note that, if satisfies constraint (14), then the expression on the right-hand side of Eq. (17) is always non-negative, i.e., the corresponding stability interval for α 2 exists.

The small-amplitude limit: coupled NLS equations
Small-amplitude solutions to Eqs. (1) and (2) for oscillatory nontopological solitons (breathers) may be looked for as where U and V are small-amplitude slowly varying complex functions, and c.c. stands for the complex conjugate. In the lowest nontrivial approximation, the complex amplitudes obey a system of coupled NLS equations, which are derived from Eqs. (1) and (2) by the substitution of ansatz (18): The dispersion relation for the linearized version of Eqs. (19) and (20) yields cf. Eq. (13), hence the zero-background solution is stable under the condition of which is not affected by coefficient β, cf. Eqs. (14)- (17). Note that the expansion of Eq. (17) for small leads to the same condition, illustrating the consistency of the analysis.
The system of coupled NLS equations (19) and (20) with β = 0 is identical to the above-mentioned model of the P T -symmetric coupler, which may be realized in terms of nonlinear fiber optics, admitting an exact analytical solution for solitons and their stability [15]. The additional terms ∼ β in the coupler model may represent the temporal dispersion of the coupling strength in fiber optics [41] (in that case, t and x are replaced, respectively, by the propagation distance and reduced time [7]). Solitons and their stability in the framework of Eqs. (19) and (20) with β = 0 constitute a separate problem which will be considered elsewhere. Here, we mention that a broad class of exact solutions to Eqs. (19), (20) can be found in the case of supersymmetry [42], α = (note that it is precisely the edge of the stability region (22)), when both the gain and loss coefficients in the two cores are exactly equal to the coefficient of the coupling between them. Indeed, in this case, the substitution of reduces Eqs. (19), (20) to a single integrable NLS equation, Equation (24) gives rise to the commonly known vast set of single-and multisoliton solutions [7], which generate the respective two-component solutions via Eq. (23). The analysis of the stability of this solutions is a subject for a separate work.

Stationary equations
Quiescent solutions to Eqs. (1) and (2), φ(x) and ψ(x), satisfy the stationary equations, which may be considered (with x formally replaced by time) as equations of two-dimensional motion of a mechanical particle of unit mass in the plane with coordinates (φ, ψ), under the action of the Lorentz force ∼ β and a force produced by an effective potential, (3) for the Hamiltonian of the SG system. As indicated above, we are interested in solutions for KK and KA complexes interpolating between different flat states, i.e., fixed points (FPs) of Eqs. (25) and (26), φ 0 and ψ 0 . The FPs are determined by equations Equation (27) gives rise to three sets of the FPs, with arbitrary integer n. In addition to that, at | | > 1/2 there also exist FPs with It is easy to see that the Hamiltonian density defined by Eq. (3) has a minimum only at FP (28), while FPs (29) and (30), (31) correspond to a maximum or saddle points, respectively, hence stable FPs may be produced solely by Eq. (28) [for this reason, the detailed stability analysis, which produces conditions (14)- (17), was presented above only for this type of the FP]. The KK and KA complexes should connect the FPs with different values of n, hence these complexes represent heteroclinic trajectories of the dynamical system based on Eqs. (25) and (26).
As concerns FP (32), it is straightforward to check that, at < −1/2, they correspond to an absolute maximum of the Hamiltonian density, therefore they are unstable. On the other hand, at > 1/2, when FP (28) is unstable, according to Eq. (14), FP (32) realizes an absolute minimum of the Hamiltonian density, hence it may represent a stable flat solution. Heteroclinic solutions connecting such FPs can be constructed, but they are different from 2π kinks. In particular, for 0 < 2 − 1 1 and β = 0, an approximate solutions connecting two FPs (32) with n = 0 and opposite signs chosen for ± is Detailed consideration of such heteroclinic solutions at > 1/2 is beyond the scope of the present work.
The gain-loss coefficient, α, does not appear in Eqs. (25) and (26), therefore it has no bearing on the shape of stationary states. Nevertheless, α does affect stability of KK and KA complexes, as shown below. This is similar to what has been shown earlier in P T -symmetric SG models in Ref. [20].

Exact KK and KA solutions for β = 0
In the case of β = 0, stationary equations (25) and (26) admit two obvious types of solutions. One of them is symmetric, with φ 0 (x) being any stationary solution of the usual sine-Gordon equation, such as the 2π kink, antikink, or periodic kink chains.
In the absence of the gain and loss terms, α = 0, the stability of the symmetric solutions with β = 0 can be investigated in the general form. Indeed, eigenmodes of small perturbations added to solution (34) can be looked for as symmetric or antisymmetric ones: where ζ is an infinitesimal amplitude of the perturbation, ω is an eigenfrequency of the perturbation mode, and φ (±) 1 (x) are the eigenmodes themselves. The stability condition is that all eigenfrequencies must be real, ω 2 ≥ 0.
The linearization of the nonstationary coupled SG equations (1) and (2) (with β = α = 0) leads to the following equations for the modal functions, φ Obviously, solutions of Eq. (37) have the same eigenvalues, ω 2 ≡ ω 2 0 , as in the case of the usual single SG equation. In particular, the usual SG 2π kink (or antikink), see Eq. (49) below, is commonly known to be stable, hence it gives rise to ω 2 0 ≥ 0. Thus, the symmetric perturbations cannot destabilize the KK complex in the coupled system.
On the other hand, Eq. (38) for the antisymmetric perturbations gives rise to eigenvalues Note that Eq. (37) has a zero-mode solution, with ω 0 = 0: (because dφ 0 /dx for the 2π kink has no zeros at finite x, the fact that eigenmode (40) corresponds to ω 2 0 = 0 confirms, via the Sturm theorem [43], that all higher-order eigenmodes generated by Eq. (37), that have zero crossings, correspond to ω 2 > 0, which corroborates the stability of the KK complex against the symmetric perturbations). Then, the substitution of ω 2 0 = 0 in Eq. (39) gives rise to ω 2 = −2 . Thus, the symmetric KK solution of the coupled SG system is unstable for > 0, and stable for < 0.
The other type of solutions to stationary equations (25) and (26) with β = 0 is antisymmetric, with φ 0 (x) being any solution of the stationary double SG equation, An exact 2π-kink solution to Eq. (42), which corresponds to the antisymmetric KA complex produced by the coupled SG system, has the known form [44]: In the limit of = 0, this solution is tantamount to the usual 2π kink. Obviously, solution (43) exists at < 1/2, which agrees with condition (14). In the limit of = 1/2, Eq. (43) takes a limit form, which is a relevant solution too, in this special case: At > 1/2, a valid solution can be obtained from Eq. (43) by analytic continuation, in the form of a spatially periodic state (essentially, a KA chain): All KA chains are unstable [22], hence solution (45) is unstable too.
The analysis of the stability of the KA complex, represented by solution (43) in the framework of coupled SG equations (1) and (2), is not analytically tractable, even in the case of β = α = 0. The respective numerical results are presented below -see, in particular, Figs. 3, 7 and 8.

Perturbative solutions for small β
If the cross-derivative coupling constant β in Eqs. (25) and (26) is treated as a small perturbation, it is easy to see that the first correction to the symmetric KK solution (34), which was obtained above for β = 0, is antisymmetric. Thus, the full (approximate) solution becomes asymmetric at β = 0: with perturbation φ 1 (x) determined by the linearized equation An exact solution to Eq. (47) can be found, making use of the zero mode (40): In particular, the unperturbed 2π kink/antikink is with polarity σ = +1/ − 1. The respective perturbed KK solution (46) is Below, these approximate analytical results are compared to their numerically found counterparts in Fig. 2.
Similarly, the first correction to the antisymmetric KA solution (41), induced by small β, must be symmetric, cf. Eq. (46), hence the corresponding approximate solution is again asymmetric: with perturbation φ 1 determined by the linearized equation (cf. Eq. (47))

Stationary KK and KA solutions and stability equations
In this section, we report results for KK and KA solutions of the full coupled system of SG equations (1) and (2) and their stability, obtained by means of numerical methods and complementing the analytical results of the previous section. Computations were performed with the help of a finite-difference scheme for the spatial derivatives, using a central-difference scheme for the first-order ones, and free-end (Neumann) boundary conditions. Figure 2 shows profiles of solutions of both the KK and KA types at > 0 and < 0. For the KK complexes, we display the central part of the φ and ψ components, and compare them to the perturbative prediction given by Eq. (50); for the KA modes, we display both components only in the numerical form, as an analytical solution of perturbative equation (51) is not available.
To study the spectral stability of the solutions, we proceed by adding small perturbations to the stationary solutions as follows: where λ is a (generally, complex) (in)stability eigenvalue, u 1,2 (x) and v 1,2 (x) being the corresponding eigenmodes of the small perturbations. The spectral stability condition is that there should not exist eigenvalues with Re(λ) > 0. The substitution of expression (52) into Eqs. (1)-(2) and the subsequent linearization to O(δ) leads to the following eigenvalue problem: where I is the identity operator. As said above, the KK and KA profiles do not depend on the gain-loss coefficient α, although their stability does depend on α, through its explicit inclusion inside the matrix of Eq. (53).
Instabilities are not only caused by the existence of localized KK or KA complexes, but can also emerge from the background if conditions (14)- (17) are not fulfilled. In particular, if only Eqs. (14) and/or (15) are violated, background eigenvalues with nonzero real parts possess a zero imaginary part. On the other hand, if solely Eqs. (16) and/or (17) are violated, the respective unstable background eigenvalues have nonzero imaginary parts, with a slight difference between the two cases: if condition (16) is violated, the eigenvalues with nonzero real parts are those which possess the largest imaginary part, whereas, if condition (17) fails, unstable eigenvalues arise with the smallest imaginary part (i.e., close to wavenumbers k = 0).

Instability of the KK and KA complexes at > 0
Both KK and KA complexes are found to be exponentially unstable at > 0, in agreement with the exact analytical result obtained above for the KK complexes in the case of β = α = 0 [see Eq. (39)]. Further, the complexes of both types exist, at > 0, below a critical point, β < β c , which depends on . For given , the critical values β c are different for the KK and KA solutions. At β = β c , the real eigenvalue responsible for the kink's instability vanishes. Actually, β c satisfies inequality (15), i.e., the existence range is smaller than the range of the background stability whenever Eqs. (16) and (17) hold. Figure 3 shows the dependence of the real part of the instability eigenvalues on β for fixed = 0.25. Figure 4 shows the dependence of β c on for both KK and KA complexes.
This instability leads to motion of the kinks, with the different components moving in opposite directions, i.e., the instability splits the KK and KA complexes, as shown in Figs. 5 and 6. Note that, at α = 0, the kinks move at constant velocities, with equal absolute values of the velocities in the two components, φ and ψ. However, at α = 0, the kinks move with monotonously varying velocities, whose absolute values in the two components are different. In the latter case, the kink in the The critical value of the cross-derivative coupling, β c , above which the KK and KA complexes do not exist, versus . Note that both complexes exist beyond the right edge (at ≥ 0.5); however, as condition (14) does not hold in that area, the flat states are unstable in it. The picture is independent of α, as only stability and dynamical properties depend of this parameter. gain component accelerates, while that of the lossy component -in a way reminiscent to observations in [20]-drops its speed and finally come to a halt.

Stable KK and KA complexes at < 0
We have also considered the complexes in the case of ≤ 0. Similar to the situation at > 0, bound states of both KK and KA types exist at β < β c . Their existence limits are also shown in Fig. 4, whereas Fig. 7 showcases the dependence of the stability eigenvalues on β for given < 0 and α = 0. It can be seen here that the imaginary eigenvalues approach 0 as β approaches β c . A drastic difference from the case of > 0 is that the KK complexes are spectrally stable whenever they exist at ≤ 0, once again in the full agreement with the analytical result for α = 0, given by Eq. (39). The KA modes are stable too if the condition ≤ 0 is supplemented by α ≤ α c , for some appropriate finite critical value of the gain-loss coefficient, α c ≤ β [i.e., α c satisfies condition (16)]. Exactly at = 0, the KK and KA only exist for β = 0 (notice that this is the uncoupled limit where φ and ψ are independent) and they are stable only for α = 0 as, if α = 0 for β = 0 the complexes cannot be stable because conditions (16)-(17) are violated. In this context, it is perhaps more precise to say that one component is subject to damping, while the other subject to pumping, with the latter featuring unstable dynamics.
The KA solutions become exponentially unstable at α > α c , as shown in Figs. 8 and 9. The latter figure shows that the instability splits the KA complex into two components. Naturally, the kink in the component (ψ) which is subject to the action of the dissipation comes to a halt, while its counterpart in the amplified component (φ) becomes a traveling one, accelerating over time. Each of the separated kinks creates its "shadow" in the mate component, in the form of a dip.
The dependence of α c on β is presented in Fig. 10. A noteworthy feature of this dependence is its non-monotonous form, with a maximum of α c attained at a particular value of β, which depends on . This maximum is caused by the fact that, on the one hand, if β falls below a critical value, then α c > β and condition (16) does not hold, hence the background (flat-state's) instability masks the exponential instability; on the other hand, if β is above that critical value, then β > α c and the latter instability manifests itself.
From a technical standpoint, it is relevant to note in passing that the finite-difference discretization of the first-order spatial and temporal derivatives introduces a number of additional, yet spurious numerical instabilities (associated with complex eigenvalues stemming from the continuous spectrum), disappearing as one approaches the continuum, infinite-domain limit; see Ref. [20] for similar examples with temporal first derivatives in P T -symmetric sine-Gordon and related systems, and Ref. [45] for such examples with spatial first derivatives in problems featuring Dirac   operators. In what we have discussed above, we have not considered these instabilities, focusing on the true dynamical features of the continuum problem.

Conclusions
We have introduced a system of coupled sine-Gordon equations, with mutually balanced gain and loss in them, which represents the P T symmetry in the system. The consideration of this system helps to understand possibilities for the implementation of the P T symmetry, which was elaborated in detail in optics, in other physical settings. Two types of coupling were included: sinusoidal terms, and the cross-derivative coupling. The former coupling corresponds to a commonly utilized interaction term between the corresponding FK (Frenkel-Kontorova) chains. The cross-derivative coupling was not considered in previously studied models. It may be generated by three-body interactions, assuming that the particles belonging to parallel chains move along different directions. The 2π KK (kink-antikink) and KA (kink-antikink) complexes were constructed in the system, and their stability was studied, by means of analytical and numerical methods. It has been found that the complexes are stable or unstable, depending on the sign of the sinusoidal coupling term. Stability regions for the complexes of both types were identified in the parameter space of the temporal gain/loss and the cross-derivative coupling strength. Simulations reveal splitting of unstable complexes into separating kinks and antikinks. The latter move at constant speed in the absence of gain and loss, and follow the dynamics imposed by the gain (acceleration) or loss (deceleration) in their presence.
There remain numerous interesting issues to address as a continuation of the present work. In particular, the possibility of the existence of stable traveling KK and KA complexes is a challenging problem. Moreover, a potentially analytical or semi-analytical study in the spirit of Ref. [46] could provide a set of guidelines on the expected motion (and stability) of the kinks in the presence of the additional terms considered herein. It is interesting too to identify stability boundaries for P T -symmetric solitons of the NLS type in the framework of Eqs. (19) and (20) with β = 0.