Discrete solitons and scattering of lattice waves in guiding arrays with a nonlinear $\mathcal{PT}$-symmetric defect

Discrete fundamental and dipole solitons are constructed, in an exact analytical form, in an array of linear waveguides with an embedded $\mathcal{PT}$-symmetric dimer, which is composed of two nonlinear waveguides carrying equal gain and loss. Fundamental solitons in tightly knit lattices, as well as all dipole modes, exist above a finite threshold value of the total power. However, the threshold vanishes for fundamental solitons in loosely knit lattices. The stability of the discrete solitons is investigated analytically by means of the Vakhitov-Kolokolov (VK) criterion, and, in the full form, via the computation of eigenvalues for perturbation modes. Fundamental and dipole solitons tend to be stable at smaller and larger values of the total power (norm), respectively. The increase of the strength of the coupling between the two defect-forming sites leads to strong expansion of the stability areas. The scattering problem for linear lattice waves impinging upon the defect is considered too.


I. INTRODUCTION
Arrays of evanescently coupled waveguides made of nonlinear materials is the fundamental model of discrete nonlinear optics [1][2][3][4]. Guided propagation of light in such arrays emulates electronic wave functions in fundamental periodic and disordered potentials of solid state physics. Therefore, discrete arrays of optical waveguides may be used to implement photonic counterparts of semiconductor devices used in electronic circuits [2,5]. On the other hand, the flexibility in the creation of virtual (photoinduced) [6,7] and permanent [8][9][10] arrayed waveguides opens the way to exploration of phenomena which are difficult to directly observe or control in solid-state settings, such as the Anderson localization [11,12] and eigenmodes of quasicrystal potentials [13]. Furthermore, the use of the material Kerr nonlinearity makes it possible to predict [14][15][16][17][18] and experimentally create [19][20][21] nonlinear modes in the waveguiding lattices -notably, various species of discrete vortex solitons.
The light propagation in arrays with embedded defects has also drawn a great deal of interest, as the use of defects can enhance the functionally of such transmission schemes [22][23][24][25][26][27]. In particular, defects carrying the balanced gain and loss, which give rise to the optical realization of the PT symmetry [28][29][30][31][32], were recently introduced too [33][34][35]. Such systems, although governed by discrete nonlinear Schrödinger equations corresponding to non-Hermitian Hamiltonians [36], give rise to entirely real propagation spectra, provided that the common strength of the gain and loss terms does not exceed a critical level, past which the PT -symmetry of the propagating eigenmodes is destroyed. The simplest version of discrete PT -symmetric nonlinear systems amounts to dimers, i.e., systems of two linearlycoupled sites, which carry the balanced gain and loss, combined with the onsite nonlinearity [37][38][39][40].
In this work, we consider one-dimensional arrays of linear waveguides with an embedded defect, formed by a pair of Kerr-nonlinear guiding cores with a tunable strength of the coupling between them. The two defect-forming waveguides carry mutually balanced linear gain and loss, thus realizing the PT symmetry in this context (a similar configuration built in continuous medium is called a PT -symmetric dipole [41]). This setting, as well as a similar one, with a pair of nonlinear waveguides side-coupled to the linear array [42][43][44], which can be readily implemented in the experiment, offer an advantage of finding symmetric [45,46] and asymmetric [47] solutions for discrete solitons in an exact analytical form. However, the PT -symmetric extension of the dual-core nonlinear defect embedded into a linear array was not addressed in previous works, although it can also be realized experimentally, by means of the methods recently developed in optics [28][29][30][31][32]. We produce both exact analytical and numerical solutions for stable fundamental and dipole-mode discrete solitons in this system (asymmetric solitons cannot be supported by the PTsymmetric defect, as they cannot maintain the equilibrium between the gain and loss). In particular, we conclude that the stability and instability regions for the discrete solitons is accurately identified by the Vakhitov-Kolokolov (VK) criterion. The rest of the paper is structured as follows. The system is introduced in Sec. II. Analytical and numerical results are presented in Sec. III. The scattering of lattice waves on the PT -symmetric defect is studied by means of direct simulations in Sec. IV The paper is concluded by Sec. V.

II. THE MODEL
The system considered in this work is displayed in Fig. 1, that features the one-dimensional linear waveguide array with the embedded pair of cores, which, in addition to the Kerr nonlinearity, supply equal amounts of local gain and loss. The nonlinearity can be induced in the inserted cores by means of appropriate dopants [48], or by making the light confinement in them much tighter than in the (quasi-)linear waveguides. The gain and loss may be imposed by means of doping too, with the external pump focused solely on the core which is selected to carry the gain.
The separation between the embedded cores is considered as a control parameter of the system. In other words, outside of the defect, the inter-core coupling coefficient is constant C 0 , while for the two cores forming the defect it is C d . It is implied that C d /C 0 > 1 and C d /C 0 < 1 correspond, respectively, to the distance between the two defect-building waveguides which is smaller or larger than the separation between the cores in the linear array.
The propagation of light in the present system is governed by the discrete nonlinear Schrödinger equation, which can be derived by means of well-known methods (the tight-binding approximation) [1]: = − (C n,n−1 u n−1 + C n+1,n u n+1 ) − γ n |u n | 2 u n + iκ n u n .
Here z is the propagation distance, u n (z) is the amplitude of the electromagnetic field in the n-th core, real C n,n±1 are the above-mentioned coupling constants, while real coefficients γ n and κ n > 0/κ n < 0 account for the nonlinearity, and for the gain/loss constants in the defect-forming cores: C n,n−1 = C d , n = N/2 C 0 , n = N/2 , γ n = γ, n = N/2, N/2 + 1 0, n = N/2, N/2 + 1 , where even N is the total number of the waveguides in the system. The total field power is Stationary solutions to Eq. (1) are looked for as where U n is the distribution of the local amplitudes, and K the propagation constant. Stability of the localized stationary modes was investigated numerically by means of computing eigenvalues for small perturbations, and the results were verified by means of direct simulations of the perturbed evolution. The perturbed solution was taken as where the asterisk stands for the complex conjugate. The substitution of this expression into Eq. (1) and linearization leads to the eigenvalue problem for the perturbation frequency, G ≡ G r + iG i , and the eigenmodes, {w n , v n }: Solution U n is stable if all the eigenvalues G are real. In the next section, we carry out the analytical and numerical study of the stationary discrete solitons in the present system. By means of rescaling, we set C 0 = 0.5, γ = 1, and produce numerical results for the system of size N = 64, while P [the total power (3)], C d and κ are left as control parameters.

A. The analytical consideration
Following [47,49] (and a similar approach for the continuous system developed in [41]), an exact solution for stationary amplitudes U n with propagation constant K in Eq. (4) can be sought for as where K and λ are linked by the dispersion relation (hence, for given C 0 , the propagation constant may only take values K > 2C 0 ), amplitude A may be assumed real, and φ is, in this form, an arbitrary phase shift. The substitution of Eqs. (4) and (6) into Eq. (1) at n = N/2 and n = N/2 + 1 leads to the final system of nonlinear equations: where, as said above, γ = 1 is fixed. The balance of imaginary terms in Eq. (8) yields two solutions for phase shift φ, hence the solution exists under the constraint of κ < C d . As usual, this means that PT -symmetric solutions exist provided that the gain/loss coefficient does not exceed a certain critical value. Then, the balance of real terms in Eq. (8) yields an expression for the squared amplitude: where the upper and lower signs correspond, respectively, to φ + and φ − in Eq. (9), and dispersion relation (7) was used to simplify this expression. Further, solution (10) with the lower sign exists at all values of λ ≥ 0, the limit case of λ = 0 corresponding not to a localized mode, but to a flat one with a constant amplitude, As concerns solution (10) with the lower sign, they exist up to λ = 0, also going over into a flat mode with a constant amplitude if the linear lattice is a tightly knit one, with the coupling constant which large enough, (obviously, the lattice with C 0 ≥ C d is always a tightly knit one). In fact, both flat states, given by Eqs. (11) and (12), are unstable because, as shown below, they are limit forms of broad discrete solitons which are unstable according to the VK criterion.
On the other hand, for a loosely knit linear lattice, with C 0 < C (0) 0 , solutions (10) with the upper sign exist above a finite value of λ, In other words, width W of the respective discrete soliton pinned to the PT -symmetric nonlinear defect takes values smaller than the respective maximum value: W ≡ 1/λ ≤ W max ≡ 1/λ min . Note that, unlike amplitudes (11) and (12) of the flat states, which remain finite at amplitude at λ = 0, A 2 + (λ) vanishes at λ = λ min . Figures 2 and 3 show typical examples of stable and unstable solutions produced, severally, by Eq. (6) with φ + and φ − (blue and red lines display Re {U n } and Im {U n }, respectively). The profiles of the solutions in Fig. 2 show that the soliton corresponding to φ + and A 2 + in Eqs. (9) and (10) may be identified as a fundamental one, with an even shape of the taller (real) component, while the profiles in Fig. 3 demonstrate that the soliton corresponding to φ − and A 2 − represents an excited state in the form of a dipole, whose taller (imaginary) component is odd, with respect to the lattice coordinate. Unstable solitons undergo a blowup in the real-time propagation, see Figs. 2(f) and 3(f).
Two solutions (10) give rise to the following expression for the total power of the soliton, calculated as per Eq. (3): In particular, the total power diverges for the dipole modes in the limit of λ → 0, when the soliton degenerates into the flat state with constant amplitude (11), and similarly for the fundamental mode in the tightly knit lattice, with C 0 > C 0 , see Eq. (13). In these cases, the power attains a finite minimum (threshold), given by Eq. (15) at values λ = λ 0 determined by a cubic equation, On the other hand, for the fundamental mode in the loosely knit lattice, with C 0 < C 0 , where the solution family starts from point (14), the P (λ) dependence commences from P (λ = λ min ) = 0.
It is reasonable to expect that a necessary stability condition for the fundamental and dipole-mode solution families may be given by the Vakhitov-Kolokolov (VK) criterion [50][51][52], dP/dK > 0, or, which is more technically convenient, Recall that P+(λ) and P−(λ) pertain, respectively, to the fundamental and dipole-mode discrete solitons. dP/d(e λ ) > 0, in terms of Eqs. (6) and (7). The application of the VK criterion in this form to expression (15) leads to the following stability condition: cf. Eq. (16). Straightforward analysis demonstrates that the family of the fundamental solitons in the loosely knit lattice, with C 0 < C (0) 0 , satisfies the VK criterion at all values of λ where such solitons exist, i.e., in the region defined by Eq. (14). On the other hand, for the family of dipole solitons, as well as for the fundamental modes in the tightly knit lattice [see Eq. (13)], the P (λ) dependence always has a minimum at λ = λ 0 , corresponding to ∆ VK (λ 0 ) = 0, see Eqs. (17) and (16). The slope of the dependence is positive, suggesting the VK stability, at λ > λ 0 , and negative, i.e.,VK-unstable, at λ < λ 0 .
The above conclusions are illustrated by Fig. 4, which displays P ± (λ) for the fundamental solitons and dipole modes in the tightly and loosely knit lattice, respectively. Below, we verify the predictions of the VK criterion for the fundamental (φ + ) and dipole (φ − ) discrete PT solitons by means of computing stability eigenvalues, and using direct simulations of the perturbed evolution.

B. Numerical results for fundamental solitons
The numerical analysis was started with search for stationary discrete solitons by means of the imaginary-time propagation method [53,54], with the purpose to check that there are no other solutions but those found above in  6) and (7), κ is the gain/loss coefficient of the PT -symmetric defect, and C d is the coefficient of the coupling between the two cores forming the defect]. The solitons do not exist in the white regions, viz., exactly as follows from the analytical solution, at κ > C d (at the top of the panels) and at λ < λmin(κ) [see Eq. (14)], i.e., in the left bottom corners of panels (a3) and (a4). In the gray area, the solitons exist but do not satisfy the VK criterion [see Eq. (17)], i.e., they are definitely unstable. The computation of the stability eigenvalues, based on Eq. (5), demonstrates that the solitons are completely stable in the red area, but unstable against perturbation modes not covered by the VK criterion in the yellow area. The white, gray and red/yellow colors have the same meaning in existence and stability diagrams displayed in other figures below. The horizontal dotted lines in panel (a3) and (a4) designate the border between regions of the tightly and loosely knit linear lattices, defined as per Eq. (13), the lattices being "loose" [i.e., with λmin > 0, see Eq. (14)] beneath these lines. In panels (a1) and (a2), which correspond to C0 ≥ C d , the linear lattices may only be tightly knit ones. (b1)-(b4) The stability and instability areas from (a1)-(a4) for the fundamental solitons, replotted in the plane of (P, κ) at same values of C d .(c1)-(c4) The same stability diagrams as in (b1)-(b4), but produced by direct simulations. In the white area, the imaginary-time simulations do not converge to stationary solitons.
the analytical form, see Eqs. (6), (7), (9), and (10). In this way, no additional solutions have been found, while the actually produced ones are fully identical to their analytical counterparts.
For the fundamental solitons, which correspond to φ + and A 2 + in Eqs. (9) and (10), the VK-stable region, as predicted by Eq. (17), is displayed in the plane of (κ, λ), for different values of C d , in Fig. 5(a1)-(a4). In these figures, the VK criterion does not hold [i.e., ∆ VK is negative, see Eq. (17)] in gray areas. As this criterion is only necessary, but not sufficient, for the stability, we have used numerically generated solutions of the eigenvalue problem, based on Eq. (5), to identified fully stable areas [red in Fig. 5(a1)-(a4)], and those (yellow ones) where the fundamental solitons are unstable, although they obey the VK criterion.
The stability regions from Fig. 5(a1)-(a4) are replotted in the (P, κ) parameter plane in Fig. 5(b1)-(b4), as total power P is a physically relevant characteristic of soliton modes. Further, Fig. 5(c1)-(c4) displays a counterpart of  FIG. 6: (a) Areas of the stability (red) and instability (yellow) for the fundamental discrete solitons in the (P, C d /C0) plane at a fixed value of the gain-loss coefficient, κ = 0.01. (b) Regions of the existence of stable symmetric (red) and asymmetric (yellow) fundamental discrete solitons in the same system with κ = 0 (no gain and loss).
the same stability diagram, produced by summarizing results of direct simulations of the perturbed evolution. The comparison clearly shows that the computation of the eigenvalues, in the combination with the VK criterion, and, on the other hand, direct simulations lead to almost identical stability areas.
An evident trend featured by these results is destabilization of the fundamental solitons with the increase of the total power. This can be explained by the fact that, at large values of P , the nonlinearity, which acts only at the two defect sites, becomes a dominant term, which makes the system close to an isolated PT -symmetric dimer [37][38][39][40], and suppresses the stabilizing effect of the linear lattice attached to the dimer.
In Fig. 5, the increase of the strength of the coupling between the two defect-forming sites, C d , makes the stability area much larger. This feature can be readily understood, as, at small C d , the action of the gain at the amplified site is weakly checked by the power flow to the adjacent dissipative site, hence perturbations can easily initiate blowup in the gain-carrying core.
The enhancement of the stability with the increase of C d is additionally illustrated by Fig. 8(a), which shows the stability areas of the fundamental solitons in the plane of (P, C d /C 0 ) for a small fixed value of the gain/loss coefficient of the PT -symmetric defect, κ = 0.01. From this figure, we conclude that the stable fundamental solitons exist above the corresponding minimum value of the coupling constant, C d /C 0 > 0.74. At P > 1.5, the stability boundary may be fitted by a simple linear expression, P = C d /C 0 . Further, in Fig. 8(b), we depict similar results for κ = 0, when the system reduces to the usual nondissipative one, with a double nonlinear defect embedded into the linear array. Recall that both symmetric [45,46] and asymmetric [47] discrete solitons can be supported in the linear lattice by the double nonlinear defect. It is seen in Fig. 6 that the stability and instability areas for the fundamental solitons in the PTsymmetric system smoothly carry over, respectively, into stability regions for symmetric and asymmetric solitons in the nondissipative system. Indeed, because stationary asymmetric solitons cannot exist in the PT -symmetric system (as the balance between the gain and loss cannot be maintained by them), the stability region for asymmetric solitons in the conservative system, where symmetric solitons exist too but are unstable [47], corresponds to the instability region in the PT -symmetric counterpart of the conservative lattice.

C. Numerical results for dipole solitons
Recall that dipole solitons correspond to φ − and A 2 − in Eqs. (9) and (10), and they always have a nonzero existence threshold in terms of the total power (unlike the fundamental solitons, whose existence threshold vanishes in the loosely-knit lattices, as shown above). Stability regions for the dipole modes in the (P, κ) plane are displayed in Fig.  7. The stability and instability in this figure are identified through the computation of eigenvalues in the framework of Eq. (5). Comparing this to similar results for the fundamental solitons, displayed in Fig. 5(b1)-(b4), we conclude that, as well as in the that case, the stability areas expands with the increase of the coupling constant C d which links the two defect-forming cores. On the other hand, the minimum (threshold) value of the total power, which is necessary for the existence of the dipole solitons, is higher than it was found above for their fundamental counterparts (in the tightly-knit lattice). This difference is quite natural, as fundamental solitons, being more compact than dipoles, need less power to build themselves. Further, on the contrary to the fundamental modes, the dipoles tend to be more stable at large values of the total power, P (although Fig. 7 exhibits alternation of stability and instability areas). This trend may be understood via the comparison of the present system with an isolated PT -symmetric dimer [37][38][39][40], to which it becomes closer with the enhancement of the nonlinear terms acting solely at the defect sites.

IV. SCATTERING OF INCIDENT WAVES ON THE DEFECT
In the framework of Eq. (1), the scattering problem for linear lattice waves with wavenumber q hitting the embedded defect corresponds to stationary solutions (4) with field U n looked for as Here, K and q are related by the dispersion equation for propagating waves, cf. Eq. (7), I, R and T being amplitudes of the incident, transmitted, and reflected waves, respectively. It is possible to fix I as a real positive constant, while T = T r + iT i and R = R r + iR i are complex. Then, substituting the ansatz based on Eqs. (18) and (19) into Eqs. (1) at n = N/2 and n = N/2 + 1, the following two complex equations for amplitudes T and R are obtained: A numerical solution of these equations produces the reflection and transmission coefficients, i.e., |R| 2 /I 2 and |T | 2 /I 2 , as well as their sum, (|R| 2 ) + |T | 2 )/I 2 , as a function of wavenumber q > 0 and PT coefficient κ, are displayed by the (q, κ) plane in Fig. 8 for different values of I. Note that the plane comprises both κ > 0 and κ < 0, which correspond, respectively, to the incident wave originally hitting the amplifying or dissipative element, according to Eq. (2). Naturally, there is no symmetry between outcomes of the scattering for κ > 0 and κ < 0 (a similar asymmetry explains unidirectional transmission in PT -symmetric systems [55][56][57]).
From the figures, we conclude that, for I = 0.1, which correspond to an almost linear system, the solutions features a highly nonconservative behavior, with the fields being strongly attenuated or, sometimes, amplified, in different parts of the (q, κ) plane. With the increase of I, the nonlinear term becomes a dominant one, effectively suppressing the nonconservative effects, which are generated by the linear gain and loss terms. In particular, at I = 10, which correspond to a strong nonlinearity, the incident wave is almost entirely transmit across the defect.

V. CONCLUSION
The objective of this work is to propose a simple experimentally relevant system, based on the linear waveguide array, with an embedded nonlinear defect in the form of the PT dipole. The system makes it possible to obtain the full family of exact solutions for fundamental and dipole solitons pinned to the defect. The fundamental solitons in tightly knit lattices, as well as the dipole modes in all cases, exist above a finite threshold value of the total power, while for the fundamental solitons in loosely knit lattices the threshold is absent. The stability of the solitons is partly predicted by the VK criterion, and is obtained in the full form via the numerical computation of eigenvalues for perturbation modes. The fundamental and dipole solitons tend to be stable at lower and higher values of the total power (norm), respectively. The stability areas strongly expand with the increase of the strength of the coupling between the two waveguides forming the defect, which can be readily explained. The scattering problem for lattice waves impinging upon the defect was solved numerically. It shows a highly non-conservative behavior in the quasi-linear limit, while the strong nonlinearity suppresses the nonconservative features in the scattering states.
It may be interesting to extend the analysis for two-dimensional arrays with an inserted defect in the form of a PT -symmetric nonlinear quadrimer, as the corresponding generalization of the one-dimensional dimer. In this case, however, analytical solutions for modes pinned to the defect are not available, hence the entire analysis should be performed in a numerical form.