{'fraction_non_alphanumeric': 0.04282688838340262, 'fraction_numerical': 0.012419274714356682, 'mean_word_length': 4.419098894871068, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlocal gap solitons in PT-symmetric periodic potentials with defocusing nonlinearity

Existence and stability of PT-symmetric gap solitons in a periodic structure with defocusing nonlocal nonlinearity are studied both theoretically and numerically. We find that, for any degree of nonlocality, gap solitons are always unstable in the presence of an imaginary potential. The instability manifests itself as a lateral drift of solitons due to an unbalanced particle flux. We also demonstrate that the perturbation growth rate is proportional to the amount of gain (loss), thus predicting the observability of stable gap solitons for small imaginary potentials.


I. INTRODUCTION
Systems with non-Hermitian Hamiltonian possessing a real eigenvalue spectrum gained interest after the early work by Bender et al., [1], who showed that the condition of Hermiticity to obtain real eigenvalues can be replaced if the system satisfies an additional condition of parity-time invariance in the form PTH = HPT, with the additional constraint that H and PT share a common set of eigenvectors, with P and T parity and time reversal operators, respectively. Moreover, it was demonstrated that quantum systems encompassing a potential such that V (x) = V * (−x) satisfy the PT symmetry. In fact, the eigenvalues remain purely real in a given subset of system parameters: a phase transition towards complex eigenvalues connected with a spontaneous symmetry breaking can be observed [2,3]. Such systems are appealing from both fundamental and practical perspectives, as evident from the vast literature in diverse areas, including quantum mechanics, integrated optics, plasmonics and metamaterials among others [4][5][6][7][8][9][10].
In optics, where the attention paid to PT theory has been higher, the interest relies mainly on the fact that the Maxwell equations can be often recast as a Schrödinger equation; the biggest advantage is the ease in realizing complex Hamiltonians both in the spatial [11] and in the temporal domains [12,13] and, consequently, the accessible experimental demonstration of the theoretical findings. Singular light properties such as nonreciprocal photon propagation have been demonstrated, leading to new technological achievements like the realization of an all-optical diode in Silicon [14]. The role of nonlinearity has been addressed as well, with the prediction of PT -symmetric solitons [15]. In optical systems, PTsymmetric potentials can be obtained by suitable tailoring of the real and imaginary parts of the refractive indices, such that the real part is a symmetric function with position and the imaginary part encompassing gain and loss terms is anti-symmetric in space (including metal insertions [3,14] and pumped system [11,16]); in matter wave systems PT -symmetric potentials can be achieved by suitable gain and loss mechanisms introduced through coupling with external reservoirs.
On the other hand, in the last decade attention has been paid to waveguide arrays and the formation of discrete solitons [17,18], including spatial gap solitons [19,20], that is, self-localized nonlinear waves with a propagation constant within the linear bandgap of the periodic structure. Light propagation through PTsymmetric periodic potentials, including exotic phenomena such as PT -symmetry breaking, nonreciprocal behavior and double refraction were investigated from the beginning, both in the linear and nonlinear regimes [15,21,22]. The existence of PT -symmetric solitons and breathers was also discussed in the context of couplers and array geometries [23]. The existence and stability of two-dimensional gap solitons in PT linear lattices with a local Kerr response was discussed in [24].
While stable gap solitons can exist in PT -symmetric linear periodic potentials in nonlinear defocusing local media, the role of nonlocality on the stability and mobility of localized solutions is still an open issue. Recently, studies were reported on solitons in PT -symmetric periodic potentials in nonlocal self-focusing media [25] and in the presence of linear defects [26], as well as on nonlocal bright solitons in defocusing media [27], including geometries with localized PT -symmetric potentials [28]. Nevertheless, the effect of nonlocality on gap solitons in PTsymmetric periodic potentials with self-defocusing was never addressed to date.
In the present work, we study the existence and stability of PT -symmetric gap solitons in defocusing nonlocal media in the presence of a periodic potential and with a diffusive-like nonlinear response, such that exhibited by, e.g., nematic liquid crystals [18,29], thermo-optic materials [30] and atomic vapors [31]. While stationary modes are obtained only on-site (i.e., corresponding to potential minima) due to lack of parity symmetry for the overall system, we demonstrate that gap solitons in PTpotentials are always unstable for any degree of nonlocality, generally undergoing oscillatory instabilities (OI). Multi-component and multi-parameter systems as well as dissipative systems are commonly associated with OI, the latter originally discussed in the context of parametrically driven damped Kerr media [32] and later extended to the generalized Thirring model [33,34] and other systems [35,36]. It is to be noted that the existence of OI in Hermitian systems results in the eventual decay or blow-up of the stationary solution [33,36]. Here we demonstrate that the existence of OI in PT -symmetric systems leads to unidirectional energy transfer from one lattice site to the other, specifically from the gain region towards the loss region.

II. DEFOCUSING NONLOCAL GAP SOLITONS IN THE BANDGAP REGION
We consider a generic field Ψ describing either the particle distribution for matter waves or the electric field for light (electromagnetic) waves. We further assume that the material is linearly inhomogeneous and can be modeled by a linear periodic potential V (x) with a nonvanishing imaginary part [21]. Finally, we study a nonlocal nonlinear medium with a diffusive response. The evolution of the field Ψ versus propagation z (time for matter waves, propagation distance for light) obeys the following system where x is the transverse coordinate and V NL is the nonlinear portion of the complete potential V + V NL . Equation (1) is the well-known Gross-Pitaevskii equation (i.e. nonlinear Schrödinger equation in optics), whereas Eq. (2) is in the form of a screened Poisson or Yukawa equation. The parameter σ (defined as a strictly positive quantity) determines the range of nonlocality: a large σ corresponds to a highly nonlocal response [37]; in particular, the width d of the Green In writing Eq. (2) we assumed a self-defocusing Kerr medium as V NL is proportional to the intensity |Ψ| 2 : a bell-shaped Ψ gives rise to a repulsive potential. For the potential to be PTsymmetric, a necessary (but not sufficient) condition is that V (x) = V * (−x) [11]: hereby we will make the ansatz V (x) = V r (x) + iV i (x) = V r sin 2 (x) − iV i sin(2x), setting to π the period of the linear potential without any loss of generality.
We look for propagation invariant solutions of Eqs. the propagation constant for light waves or µ the chemical potential for matter waves. The substitution into (1)-(2) provides the nonlinear eigenvalue problem A. Linear eigenmodes The linearized version of Eq.
(3) (i.e., with V NL = 0) supports Bloch waves of the form φ K (x) = u K (x) exp(i K x) [38] where u K (x + π) = u K (x), i.e. u is a π-periodic function of x. Figure 1 shows the behavior of the linear eigenvalues µ versus the real potential V r for purely real (V i = 0) and for PT -symmetric (V i = 0) potentials [15]. For V i = 0 the number and size of the bandgaps increase with V r due to the stronger confinement in each lattice site. For V i = 0 the two lowest bands merge for small V r (inside the shaded region), corresponding to the break-up of the PT -symmetry: the eigenvalues µ belonging to this band are complex. When V r overcomes a threshold, dependent on V i and explicitly given by V r = 2V i , the two lowest bands split, the PT -symmetry is fulfilled and the eigenvalues become purely real [15,21], corresponding to the unshaded region in Fig. 1.
The eigenvalue spectrum for the free particle case, i.e. in the absence of any potential, follows µ = K 2 /2. As well known, for V i = 0 the system spectrum exhibits gaps; in particular, for V r = 4, there are two finite bandgaps for µ < 8 (thickest line in Fig. 2(a)). Figure 2 graphs the dependence of the dispersion relation on the imaginary part V i of the potential when V r is arbitrarily fixed to 4, together with the free particle case (lowest curve). Starting from V i = 0 and then increasing V i (see Fig. 2 (c-d)), the eigenvalue spectrum remains real (PTsymmetry is conserved) until potential reduces to V (x) = V r [1 − exp(2ix)]/2: thus, for V i = V c the energy bandgaps disappear, resulting in a spectrum equivalent to that of a free particle, but shifted above by V c , i.e. µ = K 2 /2 + V c [21,27]. Conversely, when V i > V c (broken PT -symmetry) the spectrum becomes complex; furthermore, there is no energy bandgap in the first Brillouin zone, as the first bandgap shifts to the second Brillouin zone [see Fig. 2 (b)]. Such symmetry breaking in the spectrum was discussed in detail in [15,21,39].
For V i = V c and neglecting the nonlinearity, the Floquet mode u K (x) obeys (up to the end of the section we will consider a potential of period T in order to generalize our result) with T = π in our case. Equation (5) can be solved by expressing u K (x) in its Fourier series as u K (x) = ∞ n=−∞ a n (K)e 2πinx T . After defining G = K + 2πn/T , a direct substitution into Eq. (5) provides the recursive relation − µ a n = V r 2 a n−1 .
FIG. 3. Gap soliton profiles for µ = 3 versus x. The symmetric and anti-symmetric curves correspond to the real and imaginary parts, respectively. The dotted lines show the corresponding nonlinear potential VNL for two values of the nonlocality σ = 1 (left) and σ = 10 (right). Each row corresponds to Vi = 0.1 and 1.5, respectively. The yellow/white boxes sketche the qualitative trend of the linear real potential Vr(x) along x, with mid-points in the yellow and white regions corresponding to minima and maxima of the real potential, respectively; otherwise stated, the yellow regions correspond to the guiding channels.

By looking at Eq. (6) it is apparent that, if
Vr +G 2 2 − µ = 0, when one of the coefficients a n vanishes then all a n go to zero, leading to the trivial solution u K = 0. Hence, non-trivial Bloch modes exist only if G 2 = 2µ − V r , with G 2 spanning all the positive half of the real axis; thus, the eigenvalue µ has to satisfy the necessary condition For a fixed pair (µ, K) we have a non-zero u K (x) if there is a given n(µ, K) such that 0.5 (K + 2πn/T ) 2 = µ − V r /2; then the expansion coefficients a n (µ, K) are zero for n < n − 1, whereas for n > n − 1 they can be computed via Eq. (6). Summarizing, the system shows a continuous spectrum with a cutoff at µ = V c , in agreement with our numerics [27]. For the sake of simplicity and without loss of generality, hereafter we take V r = 4: Fig. 1 (b) shows that the transition from a real to a complex spectrum occurs for V i = 2.

B. Nonlinear case
We now concentrate on calculating the form assumed by shape-preserving wavepackets in the presence of nonlinearity. As is well known, gap solitons exist only in the linear bandgaps (hence the name): we will look for fundamental PT -symmetric gap solitons embedded in the first linear bandgap by choosing V i < V c to avoid symmetry breaking. We also note that, due to the presence of the odd imaginary potential V i , off-site gap solitons do not exist in our case [40]. Typical profiles of on-site gap solitons, obtained using a pseudo-spectral technique based on Chebyshev polynomials, are shown in Fig. 3 for various degrees of nonlocality and for µ close to the upper band-edge. We set a numerical window much wider than the nonlocal length d to avoid spurious effects from the boundaries [37,41]. Analogously to the case of purely real lattices [40,42], the real component of the fundamental gap soliton is mainly localized in a given guide, with tails extending towards adjacent guides and shaped so that the linear modes in these guides are excited in phase opposition with respect to the main lobe [17]. Moreover, regardless of V i , the larger the nonlocality is the lower the tails are (compare left and right columns in Fig. 3), similar to the case V i = 0 [40]. A physical explanation of this behavior relies on coupled mode theory: to exist, gap solitons require out-of-phase excitation of neighboring guides, with the nonlinear response providing the needed difference in propagation constants (energies for matter waves) between core and side channels. In the local Kerr case, a strong excitation is needed in the adjacent sites to reach the necessary nonlinear phase modulation, whereas in the nonlocal case the nonlinear perturbation induced by the mode in the core guide, spreading outwards owing to the finite size of the Green function G(x), provides the required modulation.
At variance with the case V i = 0, gap solitons have a symmetric real part and an anti-symmetric imaginary part in order to fulfill power (number of particles) conservation. In fact, in the presence of a complex potential the particle conservation for the imaginary Schrödinger equation (that is, containing a complex potential) reads where ρ = |φ| 2 is the wave intensity and the particle flux j is given by According to Eq. (9), an even ρ corresponds to an odd χ, that is, an even real part and an odd imaginary part, respectively. Physically, particles are created within the gain regions and then they diffuse (transversely) towards the loss regions in order to keep the overall (i.e., integrated along x) particle number constant. Noteworthy, the flux j is an even function, that is, a unidirectional flow takes place [15]. Equation (9) also states that the larger V i is the larger the anti-symmetric part is (see the rightmost panel in Fig. 4), due to the increase in the   Fig. 4), with the tails diminishing as the eigenvalue µ approaches the upper band-edge (dotted lines without symbols in Fig. 4), such dynamics being fully analogous to what happens in a purely real potential [40]. The phases associated with these solutions are illustrated in the last panel of Fig. 4: the phase has hyperbolic tangent profile in the central region and deviates from it when overlapping with the adjacent guides. Figure 5 plots the soliton power (number of particles for matter waves), defined as P = |φ| 2 dx. First, the soliton power increases with the nonlinear eigenvalue µ, consistently with the self-defocusing nature of the nonlinear response. The power P r,i = |φ r,i | 2 dx carried by the real and the imaginary parts of the solution is also graphed. With increasing magnitude of the complex potential, the power carried by the imaginary part also increases, in turn affecting the stability of the solutions as we will demonstrate in the following. We also note that the power carried by the soliton increases with the nonlocality σ due to the lower nonlinear effect (see the form of the Green function G and Ref. [40]), with a ratio P r /P i roughly proportional to V i but slightly dependent on σ.

III. VARIATIONAL APPROACH
Before studying the stability of the solutions, we adopt a variational analysis to look for analytical solutions and verify how well this semi-analytical approach predicts the solutions of a PT -symmetric system. We believe this is the first time that variational analysis is applied to optical PT -symmetric systems. The semi-analytical results using a variational approximation [40,43,44] are obtained from the Lagrangian of the system Eq. (1) for solutions of the form Ψ(x, z) = φ(x, z)exp(−iµz); for the conservative part the Lagrangian reads As Fig. 4 shows, the width of the intensity profile does not appreciably vary for increasing µ, the main portion of the wavepacket being confined in a single channel. At the same time, the exact soliton profiles have tails which, for µ close to the upper band, become negligibly small. This allows us to choose very general trial functions of the form where A corresponds to the field amplitude, θ is the amplitude of the phase profile and f (x) is its spatial distribution along x; A NL is the amplitude of the nonlinear potential and, finally, w b and w NL are the widths of the soliton and of the nonlinear perturbation V NL , respectively. Given that the real part of the solution is symmetric and that the imaginary part is anti-symmetric, we assume f (x) to be an odd function of x. The nonlocal nonlinear potential well is parabolic-like in the proximity to the soliton and proportional to the Green function of Eq.
(2) away from it (i.e., exp(−|x|/d)); for simplicity we take a Gaussian profile. The phase factor is accounted for in the θ term and, as shown below, can predict the observed behavior of the solutions to a reasonable extent. The standard variational approach for systems with dissipative terms can be modified as [43,45] d dz where β stands for all the parameters free to vary in our variational computation (A, A NL , P , w b , w NL , θ); we also defined Q = iV i sin(2x)φ and L C = ∞ −∞ L C dx. Using Eqs. (11)(12), the reduced Lagrangian is The first term on the RHS of Eq.
Equation (15) clearly shows that the phase associated with the solution depends on the imaginary part of the refractive index V i ; it vanishes when the potential is purely real, eventually providing a flat-phase soliton. Similarly, the variation of Eq. (13) with respect to A NL and w NL yields in full analogy with a conservative system [40]. Figure 6 (a) compares Eq. (16) with numerical solutions: the two approaches are in good agreement (flat behavior versus µ), except for a constant factor due to the different shapes of the nonlinear perturbation. Clearly, the variational approach is not able to model the broadening of V NL near the band-edge due to the presence of non-negligible tails in the actual soliton. The effects of the complex potential are apparent in the expression for w b and µ Figure 4 suggests a hyperbolic tangent profile for the transverse phase distribution, i.e. f (x). Therefore we take f (x) = tanh(x) and plot the beam width using Eq. (17) for various P and V i in Fig. 6(c). The propagation constant µ is evaluated using Eq. (18), as well, and is compared with numerical results in Fig. 6(b). The dependence on soliton power P and imaginary potential V i match in both cases: in fact, w b increases with power owing to self-defocusing, whereas the soliton broadens with V i because of the larger flux of particle to accommodate. Fig. 6 (d) plots the numerical and variational phase profiles for various V i , with differences ascribable to the erroneous evaluation of w b in the variational method: the substitution of w b from numerics into Eq. (15) provides a nearly perfect agreement.
In the highly nonlocal limit we can get another set of simplified analytical formulae for w b and µ. Following Ref. [46] we can set V NL ≈ V and the corresponding results are plotted in Fig. 6(e)-(f). The observed discrepancy in phase profile can be attributed to the fact that w b calculated using this approach is shifted by a small constant value (Fig. 6 (f)).
Using w b from numerics in Eq. (21) yields a perfect agreement with the numerically evaluated phase profile.

IV. STABILITY AND DYNAMICAL EVOLUTION
A. Linear stability analysis (LSA) The stability of the calculated gap soliton can be addressed by considering the effect of small perturbations in the form where V µ NL is the nonlinear potential computed via Eq. (4), corresponding to the soliton φ(x) once the propagation constant µ, the degree of nonlocality σ and the linear properties (i.e., V r and V i ) of the structure are selected. In the perturbative limit, i.e., neglecting nonlinear terms stemming from the added perturbation, we obtain from Eqs. (3)-(4) the following linear eigenvalue problem where we defined the operator L = 0.5∂ 2 x − V µ NL (x) − V r (x) + µ. Moreover, we introduced D = (1 − σ∂ 2 x ) −1 , capable to compute the field intensity from the nonlinear potential V NL ; in other words, D is the convolution FIG. 7. Eigenvalues λg for various µ (i.e., soliton power) and imaginary potential Vi in the complex plane (the horizontal axis corresponds to the real part, the vertical axis to the imaginary one). Blue + and red △ correspond to σ = 1 and σ = 10, respectively.
between the Green function of Eq. (2) and the intensity profile.
The solution is stable if Im(λ g ) = 0 holds for all the eigenvalues, i.e., if the system has only real eigenvalues. The existence of complex eigenvalues corresponds to an oscillatory instability with Im(λ g ) > 0 (< 0), the lat-ter implying exponentially decaying (growing) modes together with intensity oscillations while the wave evolves along z. We solved the system of Eqs. (24-25) using pseudo-spectral techniques based on Chebyshev polynomials to compute both the diffraction operator and the operator D. We chose a grid extending to 40 along x (much larger than the maximum degree of nonlocality used, i.e. σ = 10) to avoid artifacts. The grid consisted of 1001 points, the latter ensuring independence from the numerical resolution for both the eigenvalues and the eigenfunctions in the range of interest (see below).
As a first attempt, we solved system (24)-(25) for V i = 0 and any degree of nonlocality: in agreement with previous works [40], gap solitons are stable (that is, |Im(λ g )| is less than 10 −13 , well below our numerical accuracy) for every σ and µ. Next, we considered V i = 0. In the Kerr case (σ=0) we find solitons are unstable, as briefly described in Ref. [27]; here we want to address the role played by nonlocality in soliton stability. Figure 7 illustrates the behavior of λ g versus the gain/loss coefficient V i and the soliton propagation constant µ. The values of Re(λ g ) are limited to the interval [− 6 6], as values outside it are associated with high frequency noise. First, regardless of the soliton and structure parameters, there is an oscillatory instability due to the ubiquitous presence of eigenvalues with non-vanishing real as well as imaginary parts. Moreover, all the eigenvalues responsible for  Fig. 8. The first two rows correspond to λ1, the other two rows to λ2. In the first two rows, red solid and black dashed lines correspond to |p(x)| (q(x) is negligibly small compared with p(x)) for µ = 1.9 and µ = 2.5, respectively; in the plots we took only the modes with Re(λ1) > 0 and Im(λ1) < 0. In the third and fourth rows, for Vi = 0.1 red solid and black dashed lines correspond to |p(x)| and |q(x)|, respectively; for larger Vi the red solid lines with and without symbols correspond to |q(x)| for stable perturbations (|p(x)| is not plotted as it has the same behavior), whereas the black dashed lines with and without symbols correspond to growing modes. Lines without and with symbols correspond to µ = 1.9 and µ = 2.5, respectively. the OI appears in quartets, featuring the same |Im(λ g )| and |Re(λ g )| [47]. The eigenvalue distribution is quite complex, nevertheless some general trends can be observed. For V i = 0.1 the OI eigenvalues are located close to the origin of the complex plane, whereas for higher V i new branches of OI eigenvalues appear, stemming from the broken degeneracy of the purely real eigenvalues. We also note that a new quartet of OI eigenvalues appears for V i ≥ 0.5, with |Re(λ g )| close to 4.
The trend of the instability versus V i and µ can be assessed by looking at Fig. 8. To completely study the instability we computed the eigenvalues having the largest imaginary part (i.e., the maximum growth rate) in the range |Re(λ g )| ∈ [0 3] (we name it λ 1 ) and in the interval |Re(λ g )| ∈ [3 8] (we name it λ 2 ). Let us start from λ 1 : in general, the growth rate increases with µ, and at the same time |Re(λ 1 )| gets larger, that is, the leading unstable mode moves away from the origin. The former statement always holds true except for σ = 1 and V i = 0.1 or V i = 1.5: in both these cases the real part of λ 1 vanishes, i.e., there is no OI. By looking at the growth rate, for V i = 1.5 we see a rapid noise amplification, whereas for V i = 0.1 |Im(λ 1 )| is negligible, that is, the instability should be appreciable only over lengths much longer than the beam Rayleigh distance. Finally, we also note that the growth rate tends to increases with V i . The behavior of λ 2 is different. For low V i the growth rate is zero for every nonlocality σ. The real part of λ 2 does not depend on V i , whereas it drops off as µ increases. Interestingly, the growth rate associated with λ 2 is dominant with respect to λ 1 for large V i (in Fig. 8 this occurs for both V i = 1 and V i = 1.5).
We can summarize our findings as follows: except for the cases discussed above, the instability rate grows with both V i and µ. This relates to a break up of the PTsymmetry when a perturbation is added to the soliton, analogously to Ref. [21]. In fact, large V i induce the appearance of a complex spectrum, see Fig. 1 for example; in a similar way, large µ correspond to lower real potentials trapping the wave, thus leading to a reduced effective PT -breaking threshold. To validate our interpretation of soliton stability, Fig.  9 graphs the eigenfunctions corresponding to the eigenvalues λ 1 and λ 2 . For λ 1 the eigenfunctions p(x) and q(x) are strongly asymmetric with respect to x = 0: the eigenfunctions featuring a positive growth rate (i.e., Im(λ 1 ) < 0) are centered around x = −π, i.e., on the adjacent channel with respect to the gap soliton; eigenfunctions featuring a negative growth rate corresponding to an exponential decay on propagation (i.e., for Im(λ 1 ) > 0, not shown in Fig. 9) can be found by a mirror reflection with respect to x = 0. Moreover, the transverse phase profiles of p(x) and q(x) are non uniform. All these properties confirm that the instability is due to the presence of perturbation modes breaking the PT symmetry: otherwise stated, the gap soliton propagates in a sea of unstable modes, excited by differences between the actual field profile and the exact soliton shape. The LSA allows us to predict that the soliton instability takes place as an asymmetric transfer of power towards the gain regions (in our geometry x < 0), manifesting itself together with longitudinal oscillations in the field intensity due to the non-vanishing real part of λ 1 . For λ 2 the eigenfunctions are more complicated: for V i > 0.1 they appear as delocalized Bloch waves (in Fig.  9 a zoom around x = 0 is plotted) with spatial frequency dictated by Re(λ 2 ) and conserving the main properties found for λ 1 : asymmetry around x = 0, positive and negative growth rates, a non-flat phase profile. A difference exists due to their periodic profile along x: their shape is such that it is impossible to easily determine the drift direction, or its occurrence at all. Before proceeding with the dynamical stability analysis of the stationary solution, we investigate the tradeoffs between diffraction, linear lattice and self-defocusing nonlinearity on wave propagation by numerically integrating Eq. (1) with a standard Beam Propagation Method (BPM); we use splitting of the propagation operator and a Crank-Nicolson scheme for the diffraction term. We launch in the nonlinear lattice a wavepacket with spatial distribution corresponding to a gap soliton, but varying its input amplitude. This approach corresponds to adding low frequency noise to the soliton. We start with an input excitation which is one quarter of the soliton power; Fig. 10 (top row) plots the correspond- ing numerical results. Discrete diffraction is observed, with the linear PT -symmetric potential introducing a left/right asymmetry in the intensity distribution owing to the non-reciprocity of the Bloch-Floquet modes [21].
Accordingly, the asymmetry becomes prominent as the imaginary part of the refractive index V i gets larger (compare different columns in Fig. 10). The complex potential also results in longitudinal oscillations of the beam intensity versus propagation, in agreement with the presence of OI. A further increase in power to half that of the soliton reduces diffraction, but for increasing V i the symmetry in propagation is broken and the beam diffracts only in one direction (second row of Fig. 10). Increments in excitation reduce discrete diffraction until, when the power corresponds to the exact soliton, the wavepacket profile along z becomes invariant; noteworthy, in agreement with Fig. 4, the larger V i is the larger is the amount of power coupled to adjacent guides owing to a lower linear confinement.
For powers above soliton generation (Fig. 10, third row) self-defocusing becomes comparable with the linear trapping potential V r , until eventually a powerdependent breaking of PT -symmetry (see Fig. 1) occurs and, consequently, a strong left/right asymmetry [21]. The net effect is a transverse motion of the nonlinear wavepacket from one channel to the other, the direction dictated by V i : the wave is attracted towards the gain region (in our case (see section II) negative x): this is in perfect agreement with the LSA carried out in Sec. IV A.
Even though the LSA is formally valid only for small perturbations, its results can describe other minor features of the wave evolution: the oscillation period increases with V i (LSA predicts a period of about 60 for V i = 1.5, in agreement with BPM simulations) and the exponential growth gets larger with V i as well (for V i = 1.5 LSA predicts an increase equal to e over a length of 50, in agreement with BPM); the LSA can also account for the small excitation of the adjacent guide on the right side (i.e., towards the loss region) for large V i (see the third panel in the second row of Fig. 9). Finally, Fig. 10 graphs the quantity Q, defined as the quasi-power Q = Ψ(x, z)Ψ * (−x, z)dx which is conserved in the linear regime [21] together with the real power P : the power is almost conserved everywhere whereas Q changes strongly when the PT -symmetry is broken due to the nonlinear response. It is also important to understand the role played by nonlocality. To this extent, we repeated the simulations in Fig. 10 for a lower nonlocality, in particular for σ = 1, as displayed in Fig. 11: in agreement with the LSA, in this case the gap soliton is slightly more stable, with a trend opposite to that of solitons for V i = 0, both in the discrete [17] and in the continuous cases [46].
The wavepacket behavior versus input power can be explained by resorting to particle conservation, as expressed by Eq. (8). After recasting the divergence of the flux j as (∂ρ/∂x)∂χ/∂x + ρ∂ 2 χ/∂x 2 , the particle conservation expressed by Eq. (8) provides ∂ρ/∂z > 0 when According to Fig. 4, corresponding to the soliton we can set χ(x) ≈ c 0 x, i.e., the phase follows a linear trend across x, with c 0 a constant proportional to V i ; from the figure we find c 0 > 0. Equation (26) turns into (c 0 /ρ)∂ρ/∂x < 2V i (x). First, we note that the amplitude change cancels out due to the linearity of the conservation equation in ρ. Second, due to self-defocusing nonlinearity, powers exceeding the soliton case yield a broadening of the soliton intensity profile ρ (see Fig. 6), with a diminished absolute value of the derivative of ρ. This means that condition (26) is satisfied for x < 0 ( gain region), whereas it is broken for x > 0 (loss region). In other words, particles undergo accumulation in the gain region and depletion in the loss region due to the imbalance of the flux j.
Let us now consider how this affects power coupling between adjacent guides. The accumulation of particles on the left of the core guide increases the net number of particles tunneling to the next channel, this being enhanced by the defocusing character of the nonlinearity, i.e., by lowering of the Peierls-Nabarro barrier; conversely, in self-focusing media the particle flux outwards is reduced by the nonlinearity; hence, gap solitons are stable [15]. The opposite phenomenon takes place on the right side, with particles moving to the core guide from the lateral. Summarizing, the net effect is a particle motion towards negative x. Noteworthy, the field increase towards the lateral guide progressively reduces the nonlinear dephasing between them, allowing a partial back-coupling of power towards the input guide and thus inducing the oscillatory instability. The simulations demonstrate that, over several cycles, the particles tend to acquire a net motion towards negative x. Analogous phenomena occur for powers below soliton formation, but with beam broadening due to diffractive spreading.
The explanation just provided agrees with both LSA and BPM results: the instability of gap solitons is related with the transverse flux of particle, thus it is very small for low V i (gap solitons for V i = 0 are stable) whereas is enhanced as the size of the gain/loss ratio is bigger. Additionally, flux considerations allow us to understand the dependence of soliton stability on nonlocality: a higher nonlocality yields narrower solitons, that is, larger ∂ρ/∂x, in turn increasing the flux j.
C. BPM study of noise effects on soliton propagation Next we analyze the dynamical stability of propagating gap solitons by adding Gaussian noise at the input section z = 0 with a magnitude corresponding to the 0.01% of the soliton amplitude. For low V i LSA predicts |Im(λ g )| on the order of 10 −3 , that is, the instability should appear for propagation lengths larger than 10 3 (blue solid lines in Fig. 8) (for the sake of comparison with the linear regime, discrete diffraction induces an appreciable spreading at z = 20). The numerical simulations do not show appreciable changes in soliton profile up to z = 200, thus confirming the LSA results, regardless of the value of µ (first two columns in Figs. 12 -13). For larger V i LSA predicts a much higher growth rate, comprising oscillations in the intensity evolution versus z and growing modes for all σ: the numerical results in the last two columns of Figs. 12-13 confirm the predictions. Consistently with the LSA, the simulations demonstrate that, for µ closer to the edge of the bottom band, the soliton stability improves more than for µ closer to the edge of the top band, i.e., the growth rate increases with µ. The LSA predicts similar instability lengths for solutions with different σ and large V i , with instability in general increasing as nonlocality becomes larger, the latter finding being confirmed in presence of low-frequency noise, see Sec. IV B. Conversely, BPM simulations shown in Figs. 12-13 indicate that the instability drops off when nonlocality is increased, but confirming the small dependence of the growth rate from nonlocality. The observed behavior can be explained in the context of LSA: when we add high-frequency noise, we are exciting unstable modes with large Re(λ g ) (larger than 8, thus out of the range plotted in Fig. 7), encompassing a complicated distribution of the eigenvalues versus nonlocality. The last statement is confirmed by the high-frequency variations in the intensity distribution in Figs. 12-13 in comparison with the smooth behavior followed in Figs. 10-11. Finally, Figure 14 provides an estimation of the instability behavior plotting the value of V i at which the instability manifests on a distance lower than 200, when a soliton perturbed with a noise of constant amplitude (0.01), regardless of the nonlocality parameter σ and of the propagation constant µ. As discussed previously,  instability increases both with V i and µ; at the same time, broader nonlinear response helps in inhibiting soliton blow-up.

V. CONCLUSIONS
In conclusion, we investigated shape, existence curve and stability of PT -symmetric gap solitons in a selfdefocusing medium featuring a nonlocal nonlinearity and a linear periodic potential. We showed that fundamental gap solitons exist inside all the first bandgap, with appreciable tails in proximity of the top of the lowest band edge. As in the absence of gain/loss, the magnitude of the tails decreases as the nonlocal range broadens; moreover, for a given propagation constant, a higher nonlocality requires higher powers to form a soliton. We also demonstrated that solitons possess a real and imaginary parts in order to conserve particle number, with the imaginary part proportional to gain/loss terms. A variational approach suitable for the study of PT -symmetric solutions was developed, as well as an analytical method in the highly nonlocal limit, confirming the same dependence of soliton features on system parameters as showed by exact numerical solutions. We found that, in the presence of an imaginary potential, gap solitons become oscillatory unstable. Moreover, using both a linear stability analysis and BPM simulations, we showed that the perturbation growth rate changes dramatically with the imaginary potential: for large gain/loss terms the soliton shape is conserved over much shorter distances than in the case of small imaginary potentials (much shorter than the characteristic discrete diffraction length). The instability manifests mainly as a transverse particle flux (photons in the electromagnetic case) across the periodic lattice, with soliton motion towards the gain region. We also demonstrated that instability slightly changes with the response width of the nonlinearity, with the behavior being strongly dependent on the spectral contents of the applied perturbation. Finally, we demonstrated that solitons closer to the lower band-edge, i.e., with lower µ, are more stable than those with high µ, the former solitons being well within the region where PT -symmetry is fulfilled. Last, our findings can find application, for example, in the power-driven control of optical signals into a waveguide array.