{'fraction_non_alphanumeric': 0.05356427713255907, 'fraction_numerical': 0.019102923508209853, 'mean_word_length': 3.948082829683939, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 31, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear stationary states in PT-symmetric lattices

In the present work we examine both the linear and nonlinear properties of two related PT-symmetric systems of the discrete nonlinear Schrodinger (dNLS) type. First, we examine the parameter range for which the finite PT-dNLS chains have real eigenvalues and PT-symmetric linear eigenstates. We develop a systematic way of analyzing the nonlinear stationary states with the implicit function theorem at an analogue of the anti-continuum limit for the dNLS equation. Secondly, we consider the case when a finite PT-dNLS chain is embedded as a defect in the infinite dNLS lattice. We show that the stability intervals of the infinite PT-dNLS lattice are wider than in the case of a finite PT-dNLS chain. We also prove existence of localized stationary states (discrete solitons) in the analogue of the anti-continuum limit for the dNLS equation. Numerical computations illustrate the existence of nonlinear stationary states, as well as the stability and saddle-center bifurcations of discrete solitons.


Introduction
The subject of PT-symmetry and its physical implications has gained a tremendous momentum over the past few years. This field was initiated by the original proposal of C. Bender [10] who suggested that the linear Schrödinger operator with a complex-valued potential, which is symmetric with respect to combined parity (P) and time-reversal (T) transformations, is guaranteed to have real spectrum at a certain parametric regime. Thus, this was proposed as a viable alternative for the standard Hermitian quantum mechanics. Yet, it was the pioneering work in the group of D. Christodoulides both at the theoretical [23,28] and experimental [29] levels that showcased nonlinear optics as a fertile ground for the physical implementation of the PT-symmetric potentials. These efforts have motivated a wealth of recent works, especially on the physical side, addressing various aspects of continuous and discrete PT-symmetric systems. These include among others the study of the fragility of PT-symmetry in linear problems [9,26], nonlinear stationary states of few site configurations (also referred to as oligomers, or plaquettes in two-dimensional lattices) [20,21,30,32,34], as well as solitary waves and breathers in infinite systems both continuous [1,2,8,13,14,24] and discrete [12,19,31].
While the number of studies of such PT-symmetric systems both in optics [18] and in atomic physics [15,16] is rapidly growing, the volume of related mathematical works is rather limited and mostly constrained to linear problems [6,7,22,33]. It is the purpose of this paper to provide a number of rigorous results on nonlinear stationary states in PT-symmetric discrete systems. Our emphasis will be two-fold.
First, we will consider finite PT-symmetric chains of the discrete nonlinear Schrödinger (dNLS) type [17]. We examine their phase transitions (from a PT-symmetric oscillatory phase to the exponentially growing phase) when the gain and loss parameter is increased. The nonlinear stationary states bifurcate from the linear PT-symmetric states by means of a standard local bifurcation. On the other hand, we will also consider large-amplitude stationary states in an analogue of the well-known anti-continuum limit for the dNLS equation [25,27], through a suitable rescaling of the PT-dNLS equation. This rescaling enables us to use the implicit function theorem to continue stationary states from the limit, where they are effectively uncoupled and the gain and loss parameter is negligibly small.
Second, we consider the case where the finite PT-symmetric chains are embedded in the infinite nonlinear lattice of the dNLS type. Again, we will examine phase transitions of such systems and will prove that the infinite PT-dNLS lattice has a wider stability interval compared to the isolated PT-dNLS chains. We also develop a proof of the existence of localized stationary states (discrete solitons) in the PT-dNLS equation. Numerical computations illustrate the theoretical results on existence of nonlinear stationary states, as well as the stability and saddle-node bifurcations of discrete solitons.
Note that our technique allows us to prove existence of discrete solitons in the infinite PT-dNLS equation, but such discrete solitons are unstable because the phase transition in this infinite lattice occurs already at the zero value of the gain and loss parameter [26]. Earlier, existence of such discrete solitons was observed in numerical continuations from the diatomic PT-dNLS lattice [19].
The article is structured as follows. Section 2 covers fundamentals of the PT-dNLS equaiton. Section 3 is devoted to finite dNLS chains with four sections on eigenvalues of the linear PT-dNLS equation, local bifurcations of stationary states, bifurcations of large-amplitude stationary states, and numerical results. Section 4 is concerned with the PT-symmetric defects in infinite dNLS lattices and contains three sections on eigenvalues of the linear PT-dNLS equation, bifurcations of discrete solitons from the anti-continuum limit, and numerical results. Section 5 concludes the article with a summary and a discussion of future directions Acknowledgments: The authors thank James Dowdall for help at an early stage of the project during his NSERC USRA work and Dimitri Frantzeskakis for discussions on the subject of PT-symmetry. The work of P.K. is partially supported by the US National Science Foundation under grants NSF-DMS-0806762, and NSF-CMMI-1000337, from the Alexander von Humboldt Foundation and from the US AFOSR under grant FA9550-12-1-0332. The work of D.P. is supported in part by NSERC and by the ministry of education and science of Russian Federation (Project 14.B37.21.0868).

Formalism of the PT-dNLS equation
We consider the discrete nonlinear Schrödinger (dNLS) equation with non-conservative terms that introduce gains and losses of nonlinear oscillators. When gains and losses are combined in a compensated network, the model referred to as the PT-dNLS equation takes the form i du n dt = u n+1 − 2u n + u n−1 + iγ(−1) n u n + |u n | 2 u n , where parameter γ stands for the gain and loss coefficient. The finite PT-dNLS chain is defined for n ∈ {1, 2, ..., 2N } for a positive integer N subject to Dirichlet boundary conditions u 0 = u 2N +1 = 0, whereas the infinite PT-dNLS lattice is defined for all integers n on Z subject to the decay of u n to zero as |n| → ∞. The amplitudes u n for all admissible values of n are complex-valued functions of time t.
For notational consistency, we denote the sequence {u n } n∈Z of complex-valued amplitudes u n by the vector notation u. These vectors are considered in Hilbert space l 2 (Z) equipped with the inner product u, v := n∈Zū n v n and the induced norm u := n∈Z |u n | 2 1/2 . Let us formulate the evolution problem (1) in the complex Hamiltonian form where u is a collection of amplitudes u n for all admissible values of n (denoted by S) in some function space (denoted by X), the bar denotes complex conjugation, and the complex-valued Hamiltonian functional H : X → C takes the form The dynamical system (2) is said to be P T -symmetric if there is a linear real-valued t-independent operator P : X → X such that where N(u,ū) := ∇ūH(u,ū) and I : X → X is an identity operator. If u(t) is a solution of the PT-symmetric dynamical system (2) for t in a symmetric interval J := (−t 0 , t 0 ) ⊂ R with some positive t 0 , then v(t) := Pū(−t) is another solution of the same system for t ∈ J. This statement can be checked by direct substitution. This symmetry suggests the following definition of the operator T : C(J, X) → C(J, X): Note that the operator T is sesquilinear in u and nonlocal in t. The letters P and T stand for parity and time reversal transformations, which correspond to fundamental symmetries in physics. When the vector field N is linear and given by N(u,ū) := Hu associated with a linear complexvalued bounded operator H : X → X, then the P T -symmetry is expressed in the standard form Our first result is to show that the dNLS equation with compensated gain and loss terms (1) is a PT-symmetric dynamical system both for finite and infinite chains. Lemma 1. Define S := {1, 2, ..., 2N } for a positive integer N and X := C 2N . Then, the dynamical system (2) with complex-valued Hamiltonian (3) is PT-symmetric with respect to the operator P : Proof. We verify the statement with the explicit computation. Given the definition of P , we obtain Applying P again, we obtain which recovers the second identity (4).
Remark 1. The symmetry of Lemma 1 can be proven by simple reflection arguments. If the chain of oscillators has the damped site at the left end and the gained site at the right end, then since P reflects all oscillators about the middle point, the reflected chain has now the gained site at the left end and the damped site at the right end, that is, the reflected chain is equivalent to the complex conjugate chain.
Remark 2. Although P in Lemma 1 represents the fundamental physical symmetry, other choices of operator P are possible for the linear terms of the dynamical system (2)- (3). For instance, if N = 2, there exists another operator P a such that P 2 a = id andH = P a HP a , where with either a = 1 √ 5 or a = − 1 √ 5 (this statement can be easily checked by means of symbolic software). Nevertheless, the operator P a does not represent the PT-symmetry of the full nonlinear system (2)-(3) because the nonlinear term N non (u,ū) := N(u,ū) − Hu does not satisfy the second identity (4). For instance, we have hence the second identity (4) is not satisfied. Corollary 1. Let S = Z and X = l 2 (Z, C). For any fixed n 0 ∈ Z, the dynamical system (2) with complex-valued Hamiltonian (3) is PT-symmetric with respect to the operator P : l 2 (Z) → l 2 (Z) given by [P u] n = u n 0 −n , n ∈ Z.
Proof. The proof follows from the proof of Lemma 1 when S = {1, 2, ..., 2N } is replaced by S = Z and the value of n 0 ∈ Z is arbitrary.

Finite PT-dNLS lattices
We shall now consider the PT-dNLS equation (1) for the finite chain S N := {1, 2, ..., 2N }, where N is a positive integer, subject to the Dirichlet boundary conditions u 0 = u 2N +1 = 0. We study 2N eigenvalues of the linear PT-dNLS equation to find the phase transition threshold γ N , which separates the neutral stability of the zero solution for γ ∈ (−γ N , γ N ) and the linear instability of the zero solution for |γ| > γ N . We show that γ N is a monotonically decreasing sequence of N such that γ 1 = 1 and γ N → 0 as N → ∞.
We consider local bifurcations of nonlinear stationary states of the PT-dNLS equation (1) from the linear limit and prove that every simple eigenvalue of the linearized PT-dNLS equation generates a unique (up to a gauge transformation) family of the PT-symmetric stationary states in the parameter space. For γ inside the stability interval (−γ N , γ N ), this yields the existence of 2N branches of stationary states. These 2N branches are extended towards a large-amplitude limit with some intermediate bifurcations.
We characterize the number and properties of the branches of the stationary states in the largeamplitude limit and show that there exist 2 N distinct branches for any γ ∈ (−γ 1 , γ 1 ) = (−1, 1), for which |u n | 2 is large for all n ∈ S N . We also discuss existence of other branches of the stationary states, which are centered at the middle sites of S N and for which |u 1 | 2 is small in the largeamplitude limit.
These analytical results are illustrated with numerical approximations of the nonlinear stationary states of the PT-dNLS equation (1) for the finite chain with N = 1, 2, 3.

Eigenvalues of the linear PT-dNLS equation
We consider the linear stationary PT-dNLS equation on a finite chain S N := {1, 2, ..., 2N }: subject to the Dirichlet boundary conditions w 0 = w 2N +1 = 0. Compared to the PT-dNLS equation (1), the diagonal term of the discrete Laplacian operator has been included in the definition of the parameter E (see Remark 3 below). We shall find all 2N eigenvalues of the linear stationary dNLS equation (7) in explicit form, a result from which the phase transition threshold γ N is computed also explicitly.
Theorem 1. Eigenvalues of the linear eigenvalue problem (7) are found explicitly from the set of quadratic equations: In particular, all eigenvalues are simple and real for γ ∈ (−γ N , γ N ), where Proof. By writing we can rewrite the linear eigenvalue problem (7) in the equivalent form: where the boundary conditions are now y 0 = 0 and x N +1 = 0. Expressing y k from the second equation of the system (10) and substituting it to the first equation of the system, we obtain a second-order difference equation where the boundary conditions are now x 0 = −x 1 and x N +1 = 0. Using the discrete Fourier transform, we represent the eigenvector satisfying the boundary condition x N +1 = 0 in the form Parameter θ in the fundamental interval [0, π] defines uniquely the spectral parameter z := γ 2 + E 2 from the dispersion relation z := γ 2 + E 2 = 2 + 2 cos θ = 4 cos 2 θ 2 .
From the remaining boundary condition x 0 + x 1 = 0, we obtain where cos θ 2 = 0 (since {x k } N k=1 must not be identically zero). From the roots of sin θ(1+2N ) 2 , we obtain the admissible values of θ as follows: which yields the result by (11).

Remark 3.
For each eigenvalue E of the linear stationary dNLS equation (7) with the eigenvector w, there exists another eigenvalueĒ with the eigenvector Pw. This is an elementary consequence of the PT-symmetry, which produces a new solution v(t) = Pū(−t) = (Pw)e −i(E−2)t of the time-dependent dNLS equation (1) from the solution u(t) = we −i(E−2)t of the same equation. In particular, if E is a simple real eigenvalue (as in Theorem 1), then the eigenvector w can be chosen to satisfy the PT-symmetry We list some numerical values of the phase transition thresholds: Note that lim N →∞ γ N = 0.

Stationary states: local bifurcations
We shall now consider nonlinear stationary states on a finite chain S N , which satisfy the nonlinear stationary PT-dNLS equation: Ew n = w n+1 + w n−1 + iγ(−1) n w n + |w n | 2 w n , n ∈ S N , subject to the Dirichlet boundary conditions w 0 = w 2N +1 = 0. We shall work in the space X = C 2N . Assuming that the linear stationary PT-dNLS equation (7) admits a simple real eigenvalue E 0 with the eigenvector w 0 ∈ X, we shall prove the existence of a branch of the PT-symmetric stationary states w ∈ X satisfying the nonlinear stationary PT-dNLS equation (13) for E in a onesided neighborhood of E 0 . The solution branch is unique up to a gauge transformation: w → e iα w, where α ∈ R. This result corresponds to the standard local bifurcation of the nonlinear state w from the linear eigenstate w 0 , which is complicated here due to the presence of the PT-symmetry.
The local bifurcation results were considered with formal perturbation expansions by Zezyulin & Konotop [34]. Here we give a rigorous version of the same result.
Theorem 2. Assume that E 0 is a simple real eigenvalue of the linear stationary PT-dNLS equation (7) with the PT-symmetric eigenvector w 0 = Pw 0 in X = C 2N . Then, there exists a unique (up to a gauge transformation) PT-symmetric solution w = Pw of the nonlinear stationary PT-dNLS equation (13) for real E > E 0 . Moreover, the solution branch is parametrized by a small parameter a such that the map R ∋ a → (E, w) ∈ R × X is C ∞ and for sufficiently small a, there is a positive constant C such that Proof. We write the nonlinear stationary PT-dNLS equation (13) in the abstract form where H : X → X is the linear (matrix) operator associated with the right-hand side of the linearized stationary PT-dNLS equation (7) and N non (w) : X → X is the cubic nonlinear part. We note that according to our assumptions, we have Using the standard Lyapunov-Schmidt method, we write where (∆, a, u) ∈ C×C×X are determined from the nonlinear equations (15) The projection to Ker(E 0 − H) + is written in the scalar form: By the implicit function theorem, the projection to Ran(E 0 − H) + (not written here) guarantees the existence and uniqueness of a smooth (C ∞ ) map from (∆, a) ∈ C 2 to u ∈ Ran(E 0 − H) ⊂ X.
Moreover, for small values of ∆ and a, there is a positive constant C such that For a = 0, we have a unique zero solution u = 0 and the equation (17) is satisfied identically. In what follows, we assume a = 0.
We claim that P w 0 , w 0 = 0 under the assumption that E 0 is a simple eigenvalue of H. Indeed, if P w 0 , w 0 = 0, there exists a generalized eigenvector w 1 ∈ X for the same eigenvalue E 0 from a solution of the inhomogeneous equation which is a contradiction to the assumption that E 0 is a simple eigenvalue of H.
Therefore, P w 0 , w 0 = 0. Then, there exists a unique smooth map from a ∈ C to ∆ ∈ C solving the bifurcation equation (17). Moreover, for small values of a, there is a positive constant C such that Note that both P w 0 , w 0 and P w 0 , N non (w 0 ) are real because of the PT symmetry of the eigenvector w 0 = Pw 0 and the nonlinear field N non satisfies the second identity (4). Indeed, we have Therefore, ∆ is real at the leading order O(|a| 2 ). To exclude the gauge transformation, let us consider the real values of a. Because the nonlinear vector field N non preserves the PT-symmetry, the unique solution for u and ∆ is P T -symmetric, so that u = Pū and ∆ is real. The bound (14) follows from (16), (18), and (19). To be precise, we obtain It remains to prove that ∆ 2 > 0. However, using the explicit representation from Theorem 1, for the eigenvalue with θ = 2πj 1+2N , 1 ≤ j ≤ N , we obtain the eigenvector w 0 with components Therefore, > 0, and the proof of the theorem is complete.
Remark 4. The local bifurcation results do not apply in the limit N → ∞ because of two reasons. First, the spectrum of the linear stationary dNLS equation (7) becomes continuous as N → ∞. Second, for any γ = 0, the spectrum includes complex (purely imaginary) points of E because γ N → 0 as N → ∞.
Let us consider the simplest example N = 1 when Theorem 2 works. The two simple eigenvalues are E 0 = ± 1 − γ 2 and the eigenvectors w 0 = Pw 0 are given by the same expression In this case, ∆ 2 = 3 4 , so that we have the expansion In fact, it follows from the exact solution (30) below that the error term O(A 4 ) is identically zero.

Stationary states: bifurcation from infinity
We shall now consider the stationary states of the nonlinear stationary PT-dNLS equation (13) in the limit of large values of E. This corresponds to the anti-continuum limit of weak couplings in the PT-dNLS lattice after a suitable scaling transformation (which is also discussed in [19]). Note that the standard anti-continuum limit arising when the coupling parameter in front of the discrete Laplacian operator vanishes fails to generate any solutions of the stationary dNLS equation (13) for real values of E and γ = 0.
We shall develop methods to analyze a bifurcation from infinity for solution branches. In particular, we shall prove the existence of 2 N branches of the PT-symmetric stationary states w of the nonlinear stationary PT-dNLS equation (13) for γ ∈ (−1, 1) and for large values of E, for which |w n | 2 is large for all n ∈ S N . The solution branches are unique up to the gauge transformation w → e iα w with α ∈ R. The complication of proving this result is caused by the degeneracy of asymptotic solutions of the nonlinear algebraic system (13) in the limit E → ∞. Indeed, setting w = √ EW and taking the limit E → ∞, we obtain an uncoupled set of algebraic equations with N PT-symmetric solutions where ϕ k ∈ R is an arbitrary parameter and e k is a unit vector on the finite chain S N . However, the space of solutions of the nonlinear algebraic system (13) in the limit E → ∞ does not enjoy the linear superposition principle and parameters {ϕ k } N k=1 must be fixed from O(1) conditions as E → ∞. To prove persistence of continuations of the limiting roots for large but finite values of E, we have to unfold the degeneracy of the nonlinear system by a special transformation, after which the result is guaranteed by the implicit function theorem. Along these lines, we prove the following main result.
Theorem 3. For any γ ∈ (−1, 1), the nonlinear stationary PT-dNLS equation (13) in the limit of large real E admits 2 N PT-symmetric solutions w = Pw (unique up to a gauge transformation) such that, for sufficiently large E, the map E → w is C ∞ at each solution and there is a positive Proof. We set E = 1 δ and w = W √ δ for small positive δ and write the stationary dNLS equation (13) in the equivalent form: subject to the Dirichlet boundary conditions W 0 = W 2N +1 = 0. We consider a PT-symmetric solution with W = PW such that the system can be closed at N algebraic equations for 1 ≤ n ≤ N subject to the reflection boundary condition W N +1 =W N . Note that parameter γ ∈ R is fixed.
Case N = 1: In this case, we only have one nonlinear algebraic equation to solve: Setting 1 e iϕ 1 , we separate the real and imaginary parts of equation (22) as follows: For any γ ∈ (−1, 1), there exist two solutions for ϕ in [0, π] from the second equation written as sin(2ϕ 1 ) = −γ. For each ϕ, we have a unique solution of the first equation written as Case N ≥ 2: Let us now unfold the degeneracy of the nonlinear algebraic system (21) in the limit δ → 0 by using the transformation: where amplitudes A 1 ,A 2 ,...,A N and phases ϕ 1 , ϕ 2 , ..., ϕ N are all real. After substitution and separation of real and imaginary parts, we obtain N equations for phases and N equations for amplitudes For δ = 0, the system of amplitude equations (25) has a unique solution at the point A 1 = A 2 = · · · = A N = 1. The vector field of the nonlinear system is smooth with respect to (A 1 , A 2 , . . . , A N ) and δ near this point for all (ϕ 1 , ϕ 2 , . . . , ϕ N ) ∈ T N , where T denotes the fundamental interval [0, 2π] subject to the periodic boundary conditions. The Jacobian matrix with respect to (A 1 , A 2 , . . . , A N ) at this point has eigenvalue 1 of geometric multiplicity one and algebraic multiplicity N . By the Implicit Function Theorem, for all (ϕ 1 , ϕ 2 , . . . , ϕ N ) ∈ T N and small δ ∈ R, there is a unique solution of the nonlinear system (25) such that the map (ϕ 1 , ϕ 2 , . . . , ϕ N , δ) → (A 1 , A 2 , . . . , A N ) is C ∞ and there is a positive δ-independent constant C such that Bound (20) follows from this bound and the scaling transformation. Now we consider the system of phase equations (24), which is δ independent. Nevertheless, it depends on δ via amplitudes (A 1 , A 2 , . . . , A N ). For δ = 0, the nonlinear system (24) can be written in the explicit form: Denote ψ := 2(ϕ 1 + ϕ 2 + · · · + ϕ N ). For any γ ∈ (−1, 1), there are 2 N possible solutions of (27) for (ψ, ϕ 2 , . . . , ϕ N ) ∈ T N , depending on the binary choice of the roots of the sinusoidal functions on the fundamental period. Because ϕ 1 = ψ 2 − ϕ 2 − · · · − ϕ N , there are actually four solutions for ϕ 1 in T, however, the solutions with ϕ 1 ∈ (π, 2π] are reducible to the solutions with ϕ 1 ∈ (0, π] by the transformation W → −W, which is a particular case of the gauge transformation. In what follows, we only consider the two possible solutions for ϕ 1 in [0, π]. The vector field of the nonlinear system (24) with (A 1 , A 2 , . . . , A N ) obtained from the nonlinear system (25) is smooth with respect to (ϕ 1 , ϕ 2 , . . . , ϕ N ) and δ. The Jacobian matrix with respect to (ϕ 1 , ϕ 2 , . . . , ϕ N ) for δ = 0 is given by the matrix Now it is clear that cos(ϕ n ) = 0 for all 2 ≤ n ≤ N if γ ∈ (−1, 1). In addition, the last equation in the system (27) is given by either sin(ψ) = 0 if N is even or sin(ψ) = γ if N is odd. In either case, cos(ψ) = 0 if γ ∈ (−1, 1). Hence, the Jacobian matrix is invertible if γ ∈ (−1, 1). By the Implicit Function Theorem, for all small δ ∈ R, there is a unique continuation of any of the 2 N possible solutions (ϕ * 1 , ϕ * 2 , . . . , ϕ * N ) of the nonlinear system (27) as a solution of the nonlinear system (24) such that the map δ → (ϕ 1 , ϕ 2 , . . . , ϕ N ) is C ∞ and there is a positive δ-independent constant C such that This completes the proof of the theorem.
Remark 5. The number of solution branches grows as N → ∞ for any fixed value of γ in the interval (−1, 1). However, all these solution branches are delocalized in the sense that |w n | 2 ≈ E as E → ∞ for all n in S N . Therefore, none of the solution branches of Theorem 3 approach to a localized state (discrete soliton) as N → ∞. . If N → ∞, such stationary states approach to a localized state (discrete soliton). Note that the discrete solitons are unstable on the unbounded lattice because the continuous spectrum of the linearized dNLS equation (7) is complex for any γ = 0, recall that γ N → 0 as N → ∞.
Remark 7. The arguments of the implicit function theorem can not be applied to construct solution branches which are centered anywhere but at the central sites in the set S N . Indeed, the numerical results below show that no such solution branches exist for large values of E.

Numerical results
We shall construct here the simplest nonlinear stationary states for N = 1, 2, 3. For N = 1, this corresponds to the nonlinear dimer, where the solution branches can be obtained analytically, as in [20,28,32]. For N = 2, this corresponds to the nonlinear quadrimer and the solution branches can be at best approximated numerically [20,34]. For N = 3, the numerical approximations of the nonlinear stationary states are added here for the first time.
For N = 1, we use the reduction w 2 =w 1 and write w 1 = Ae −iϕ with real A and ϕ. Then, the stationary PT-dNLS equation (13) yields two equations With two solutions of the first equation for ϕ ∈ [0, π], we obtain two solution branches where γ 1 = 1. The two solution branches coalesce into one branch for γ = γ 1 and disappear via a saddle-center bifurcation for γ > γ 1 .
Positivity of A 2 ± shows that E > E ± = ± 1 − γ 2 , where E ± are the simple eigenvalues of the linear stationary PT-dNLS equation for γ ∈ (−γ 1 , γ 1 ). We note that A 2 ± → 0 as E → E ± and that A 2 ± ∼ E as E → ∞. These analytical results clearly illustrate the bifurcation results in Theorems 2 and 3.
For N = 2, we use the reduction w 4 =w 1 and w 3 =w 2 . Writing w 1 = Ae −iϕ−iψ and w 2 = Be −iψ with real A, B, ϕ, and ψ in the nonlinear stationary dNLS equation (13), we obtain the system of nonlinear equations: and The roots of the algebraic system (31) and (32) can be investigated numerically and the results depend on the value of γ. Figure 1 shows the solution branches on the (E, A 2 )-plane (top) and the (E, B 2 )-plane (bottom) for γ = 0.5 < γ 2 ≈ 0.618 (left), γ = 0.75 (middle), and γ = 1.1. Note that no solution branches exist for γ > γ * 2 := 2 cos π 5 ≈ 1.618, because no simple real eigenvalues occur in the linearized dNLS equation (7) for these values of γ.
According to Theorem 2, we count exactly four (2N = 4) solution branches for small amplitudes A and B for γ < γ 2 and exactly two small solution branches for γ 2 < γ < γ * 2 . According to  For N = 3, we write w 1 = Ae −i(ϕ+ψ+θ) =w 6 , w 2 = Be −i(ψ+θ) =w 5 , and w 3 = Ce −iθ =w 4 . The roots of the resulting algebraic system are investigated numerically by a homotopy method and the results are shown on Figure 2. We count six (2N = 6) branches in the small-amplitude limit if γ < γ 3 ≈ 0.445 and four branches if γ 3 < γ < γ * 3 ; = 2 cos 2π 7 ≈ 1.247. We also count eight (2 N =3 = 8) branches for large amplitudes A 2 if γ < 1 and no branches for large amplitudes A 2 if γ > 1. More branches are counted for large amplitudes B 2 and even more branches for large amplitudes C 2 . Overall, the results for N = 3 are similar to the results for N = 2.

PT-symmetric defects embedded in infinite PT-dNLS lattices
We shall now consider an infinite PT-dNLS lattice, where the particular emphasis is on the existence and stability of localized stationary states (discrete solitons). Because the phase transition in the PT-dNLS equation (1) on the infinite lattice occurs already at γ N →∞ = 0, there is no way to obtain stable discrete solitons in such systems with extended gain and loss [26]. Therefore, we modify the PT-dNLS lattice by considering the PT-symmetric potential as a finite-size defect. Such defects were considered recently in the physical literature [31] and [4,5].
Let N be a positive integer and S N := {1, 2, ..., 2N } be the sites of the lattice, where the PT-symmetric defects are placed. The model takes the form i du n dt = u n+1 − 2u n + u n−1 + iγ(−1) n χ n∈S N u n + |u n | 2 u n , n ∈ Z, where χ n∈S N is a characteristic function for the set S N . When N = 1, the PT-dNLS equation (33) corresponds to the embedded dimer in the infinite PT-dNLS lattice. When N = 2, it corresponds to the embedded quadrimer, and so on. We study the linearized PT-dNLS equation and find the phase transition thresholdγ N . It is quite remarkable thatγ N > γ N for any N ∈ N, in particular,γ 1 = √ 2. Nevertheless,γ N is still a monotonically decreasing sequence of N such thatγ N → 0 as N → ∞.
Then, we employ the large-amplitude (anti-continuum) limit of the PT-dNLS equation (33) to study the existence of discrete solitons supported at the PT-symmetric defect S N . For recent results on existence of discrete solitons in the anti-continuum limit for the regular dNLS equation (in the absence of PT-symmetry), see e.g. [3,11]. We find that for all γ ∈ (−1, 1), 2 N branches of the discrete solitons exist in this limit, for which |u n | 2 is large for all n ∈ S N .
The existence and stability of discrete solitons is illustrated numerically and we show that the stable branches of the discrete solitons for γ = 0 originate from the stable branches in the Hamiltonian version (γ = 0) of the dNLS equation [25,27].

Eigenvalues of the linear PT-dNLS equation
We consider the linear stationary PT-dNLS equation: Because the PT-symmetric potential is compact, the continuous spectrum of the linear PT-dNLS equation (34) is located for E ∈ [−2, 2]. Besides the continuous spectrum, isolated eigenvalues may exist outside the continuous spectrum. To characterize isolated eigenvalues, we introduce a parametrization E := 2 cos θ, Re(θ) ∈ [−π, π], Im(θ) > 0, and look for exponentially decaying solutions of the linear PT-dNLS equation (34) in the form: which still leaves a set of 2N unknown variables {w n } n∈S N . To find {w n } n∈S N , we close the linear eigenvalue problem at the algebraic system subject to the boundary conditions Note that each eigenvalue E is complex if Im(θ) > 0 and there exists a complex conjugate eigenvalueĒ by the PT-symmetry (see also Remark 3). The following result is similar to the result of Theorem 1.
and persists for |γ| > γ N,k except possibly finitely many points on any compact interval of γ, where the pair coalesces into a double (semi-simple) pair of real eigenvalues. In particular, no complex eigenvalues exist for γ ∈ (−γ N ,γ N ), wherẽ Proof. We set and rewrite the linear eigenvalue problem (37) in the equivalent form: subject to the modified boundary conditions y 0 = x 1 e iθ and x N +1 = y N e iθ . Expressing y k from the second equation of the system (40) by and substituting these expressions to the first equation of the system, we obtain a second-order difference equation where the boundary conditions are now Note that equations (41) are obtained by multiplying every term by 2 cos θ − iγ, which is hence supposed to be non-zero. The second-order difference equation (41) admits an exact solution where α is defined from the transcendental equation and (C + , C − ) are non-zero solutions of the linear system following from the boundary conditions (42): Note that for fixed γ, the values of θ are obtained from the characteristic equation for this linear homogeneous system, after the values of α are excluded from the transcendental equation (44). Also note that only the values of θ with Im(θ) > 0 determine isolated eigenvalues of the linear stationary PT-dNLS equation (34) by means of the representation (35). After some algebraic manipulations, the characteristic equation for the linear system takes the form (cos(αN ) sin θ sin α − i sin(αN ) cos θ(1 − cos α)) (γ + 2i cos θ) = 0.
The equation γ + 2i cos θ = 0 gives an artificial root corresponding to the value α = π, because the second-order difference equation (41) is obtained by multiplying every term by 2 cos θ − iγ. Therefore, we drop this nonzero factor from the characteristic equation and reduce it to the transcendental equation Case N = 1: Equation (45) yields e 2iθ = 2 cos α − 1. When this constraint is used in equation (44), we obtain e −2iθ = 1 − γ 2 . Since γ ∈ R, we obtain the existence of a simple eigenvalue with Re(θ) > 0 for γ 2 > 2 and the bifurcation occurs for γ = √ 2 and corresponds to θ = ± π 2 , when cos θ = 0.
Case N ≥ 2: In a general case, we first consider bifurcations of complex values of θ from real values of θ. Hence, we set θ ∈ R and realize from (44) that Im(cos α) = 0, which implies either α ∈ R or α ∈ πk + iR for any k ∈ Z. In both cases, the characteristic equation (45) with real θ implies that either θ = 0 or θ = ± π 2 . If θ = 0, then the characteristic equation (45) reduces to the equation sin(αN ) = 0, whereas equation (44) implies that γ 2 = −2(1 − cos α) ≤ 0, which is outside of the parameter range we are interested in. (Recall here that no eigenvalues with Imθ > 0 exist in the self-adjoint case with γ = 0.) On the other hand, if θ = ± π 2 , then the characteristic equation (45) reduces to the equation cos(αN ) = 0 (recall here that the artificial root α = π is neglected) or equivalently, From equation (44), we obtain that the bifurcation occurs at which corresponds to the values (38). Next, we show that the values with Im(θ) > 0 correspond to the range γ 2 > γ 2 k . This implies that a new isolated eigenvalue E of the linear stationary PT-dNLS equation (34) with Im(E) = 0 bifurcates from the value E = 0 that corresponds to θ = ± π 2 at γ = ±γ k and persists for all |γ| > γ k . A symmetric complex-conjugate eigenvalueĒ exists by the PT symmetry of the stationary PT-dNLS equation (34).
To show the above claim, we use the parametrization z := e iα , ζ := e 2iθ , and s := γ 2 . The system of transcendental equations (44) and (45) becomes the system of algebraic equations: We know that for s = s 0 := γ 2 k , there exists a solution of this algebraic system for ζ = −1 and z 2N = −1. Therefore, we consider the continuation of this complex-valued root in real values of s. By computing the derivative in s, we obtain For s = s 0 , we obtain from this linear system that and since ζ = e 2iθ , this proves that d ds Im(θ)| θ=± π 2 > 0. By continuity of roots of the algebraic system above, Im(θ) remains positive for s > s 0 , that is, for γ 2 > γ 2 k , near the bifurcation point. We have shown that roots with Im(θ) > 0 bifurcate from the points θ = ± π 2 and remain in the upper half-plane for all γ 2 > γ 2 k . These roots correspond to complex eigenvalues E if Re(θ) = 0 (in which case E > 2) or if Re(θ) = ±π (in which case E < 2). Both situations can not be a priori excluded, however, we have here two facts: • When Re(θ) = 0 or Re(θ) = ±π, a pair of complex-conjugate eigenvalues E coalesce at the real line into a double (semi-simple) eigenvalue, because of the PT-symmetry with Im(θ) = 0 generates two linearly independent eigenvectors for the same real eigenvalue.
• Roots θ are analytic with respect to parameter γ by analytic dependence of roots of the algebraic system.
Combining these two facts together, we realize that that the double (semi-simple) eigenvalues E can not split along the real axis (as each real eigenvalue after splitting would then become a double eigenvalue by the PT-symmetry). And they can not persist on the real axis as a double eigenvalues because of analyticity of the parameter continuation of roots θ with respect to γ. Therefore, these double real roots split back to the complex domain. In addition, analyticity of the parameter continuations guarantees that there are finitely many points on any compact interval in γ, where the pairs of complex-conjugate eigenvalues can coalesce at the real axis. The proof of the theorem is complete.
We list some numerical values of the phase transition thresholds:

Stationary states: bifurcations from the anti-continuum limit
We now consider the stationary states, which satisfy the nonlinear stationary PT-dNLS equation We explore the large-amplitude limit similarly to Section 3.3. Hence, we set E = 1 δ and w = W √ δ , where δ is a small positive number. The stationary PT-dNLS equation (46) is rewritten in the equivalent form: We consider the PT-symmetric solutions with W = PW, where P is given by [P W] n = W 2N +1−n . Note that the choice of n 0 in P given by Corollary 1 is adjusted to the center of the PT-symmetric defect. Therefore, the existence of the PT-symmetric solutions can be considered in the framework of the following two subsystems: and subject to the boundary condition W N +1 =W N . The following result is similar to the result of Theorem 3.
Theorem 5. For any γ ∈ (−1, 1), the nonlinear stationary PT-dNLS equation (47) in the limit of small positive δ admits 2 N PT-symmetric solutions W = PW ∈ l 2 (Z) (unique up to a gauge transformation) such that, for sufficiently small δ, the map δ → W is C ∞ at each solution and there is a positive δ-independent constant C such that where W 0 is a solution of Theorem 3 (after rescaling).
Proof. We consider small solutions of the subsystem (49) for a given W 1 ∈ C and small δ ∈ R.
Parameter γ ∈ R is fixed. The nonlinear system represents a bounded C ∞ map from ( where Z − is the set of negative integers including zero. For δ = 0 and arbitrary W 1 ∈ C, W − = 0 is a root of the nonlinear map and the Jacobian with respect to W − at W − = 0 is invertible. By the Implicit Function Theorem, for all W 1 ∈ C and small δ ∈ R, there is a unique solution of the nonlinear system (49) such that the map (W 1 , δ) → W − is C ∞ and there is a positive δ-independent constant C such that W − l 2 (Z − ) ≤ C|δ|. Substituting W 0 from the map constructed above to the first equation of the subsystem (48), we close the system at N nonlinear equations for {W n } 1≤n≤N . The only difference from the N nonlinear equations considered in the proof of Theorem 3 is the boundary condition for given W 0 , however, W 0 is small as W 0 = O(δ). The two applications of the Implicit Function Theorem developed in the proof of Theorem 3 apply directly to our case and yield the assertion of this theorem. The bound (50) follows from bounds (26) and (28). We give details of the perturbative expansions for the two (most fundamental) discrete solitons supported by the dimer defect for N = 1. The subsystems (48) and (49) are rewritten explicitly as follows: The perturbation expansion allows us to compute at the leading order W (0) n = e −iθ δ n,1 , where θ ∈ [0, π] is arbitrary at this point. At the O(δ) order, we obtain the equations: 1 := Ze −iθ and rewrite (55) as −(Z +Z) = e 2iθ −iγ. The solvability condition is sin(2θ) = γ and it gives exactly two values for θ ∈ [0, π] for any γ ∈ (−1, 1). Then, the first-order correction term is found explicitly as follows: The perturbation expansion (53) can be continued to higher orders of δ thanks to C ∞ smoothness in Theorem 5. Hence, we obtain two branches of soliton states supported by the dimer defect.
Remark 9. The phase transition threshold for S 1 isγ 1 = √ 2, whereas the soliton states of Theorem 5 are only constructed for γ ∈ (−1, 1) in the limit E → ∞. Numerical studies (see Figure 5) show that the solution states exist for γ ∈ (−1, 1) for any fixed value of E.

Numerical results
We now test these analytical results for the PT-symmetric chain with embedded defects.
First, we consider the linear limit of such chains and examine the corresponding PT-symmetry phase transitions expected to occur atγ N . In particular, in Fig. 3, we consider the case of N = 1, i.e., a single embedded dimer which has been predicted to have a PT-phase transition atγ 1 = √ 2, for the infinite lattice. In the figure, we can clearly discern the relevant transition (the corresponding vertical dashed line shows the theoretical prediction ofγ 1 ). Nevertheless, for the finite lattice considered (here 800 sites are used), an additional bifurcation is observed at γ = 1 (see also the works of [5,26]). This bifurcation is suppressed at the infinite lattice limit, but for a finite chain a "bubble" of complex eigenvalues arises around E = 0 (which is the middle point of the spectral band). This bubble keeps expanding and slowly increasing in imaginary part between γ = 1 and γ = √ 2 (notice that this growth is barely visible in the linear scale of the figure but it is noticeable in the logarithmic scale of the inset). At the latter critical point, the rapid growth of the isolated unstable eigenvalue pair becomes dominant for the instability of the lattice with an embedded dimer.
Similar conclusions can be drawn for the case of embedded quadrimer (N = 2) and hexamer (N = 3) from Fig. 4. The figures reveal, however, that in this case in addition to the actual (infinite chain) critical points ofγ 2 ≈ 0.765 andγ 3 ≈ 0.518, respectively, there are multiple additional points of weak instability emergence due to finite size effects. Such features are noticeable due to bifurcations of instability bubbles at the edges of the spectral band (at E = 2 and E = −2), at  γ ≈ 0.46 for N = 2 and at γ ≈ 0.30 for N = 3. Additional bubbles emerge near the middle point of the spectral band (at E = 0), at γ ≈ 0.61 for N = 2 and at γ ≈ 0.45 for N = 3.
Next, we turn to the existence of nonlinear stationary solutions in the PT-symmetric dimer (N = 1) embedded in the nonlinear chain, as shown in Fig. 5. We have identified two branches of localized states starting from the Hamiltonian limit where 2θ = 0 and 2θ = π (i.e., in-phase and anti-phase discrete solitons, respectively). Indeed, the left panel shows the phase difference between w 1 and w 2 for the PT-symmetric embedded dimer. These results clearly illustrate that at the Hamiltonian limit of γ = 0, the system starts from the two well-known in phase and outof-phase solutions [25]. The former (lower norm one) is unstable for our focusing nonlinearity, while the latter one is spectrally stable. Interestingly, in accordance with what is known for an isolated dimer [20,28,31], these two solutions disappear in a saddle-center bifurcation at γ = 1. In fact, both this increased proximity and the eventual collision and disappearance are captured very accurately by the solvability condition sin(2θ) = γ. The resulting angle from the numerical computation and from the essentially coincident analytical prediction are shown in the left panel of Fig. 5. Note that while the results are shown in Fig. 5 for δ = 0.01, they remain similar not only for smaller values of δ (such as e.g. 0.001), but even for larger values up to δ = 1 that were considered. In particular, the point of the saddle-center bifurcation γ = 1 has been found to be identical for other values of δ. This confirms the observation made in Remark 9.
In the Hamiltonian case of γ = 0, our linear stability results for these branches fall back on the analysis of [25]. In fact, we retrieve the exact same condition, for the leading order eigenvalue correction, namely λ 2 /δ = 4 cos(2θ). This expression reveals the instability of the in-phase mode and the stability (for our focusing nonlinearity) of the out-of-phase one. In the presence of gain/loss (i.e., for finite γ), the fundamental difference lies in the existence condition which mandates that sin(2θ) = γ and therefore the eigenvalues of the in-phase unstable state approach the origin (from the real axis), as do the ones of the out-of-phase marginally stable state (from the side of the imaginary axis) as γ is increased towards 1. These eigenvalue pairs end up colliding at the origin at the point of the PT-phase transition at γ = 1. The trajectories of the two sets of eigenvalues are shown in the right panel of Fig. 5. It can be seen that similarly to the Hamiltonian case of [25], the agreement is better for the unstable branch of real eigenvalues. Nevertheless, for both branches the comparison of computation and analysis is highly favorable.
For the case of N = 1, we now turn to numerical simulations, in order to briefly discuss the difference between the manifestation of the instability of the in-phase localized states for different values of γ. Two prototypical examples are shown in Fig. 6. Both panels illustrate the evolution of the two most central, maximum amplitude sites of the solution. In the Hamiltonian case of γ = 0 shown in the left panel, it can be seen that small perturbations give rise to an amplified symmetry breaking in the dynamics. Nevertheless, the conservative nature of the ensuing dynamics ascertains the rapid saturation of this symmetry breaking and the eventual oscillations that arise lead to a (nearly) periodic alternation between a symmetric and a symmetry broken state. On the other hand, for γ = 0.5 we observe a drastically different manifestation of this instability shown in the right panel. More specifically, the site associated with gain grows indefinitely in a nearly exponential form, as illustrated in the inset. At the same time, the site associated with damping decreases in amplitude in a similar fashion.
Lastly, we address the bifurcation and stability results for the case of N = 2 i.e., for an embedded quadrimer. In this case, the eigenvalues of the linearization problem are shown in Fig. 7. We have  examined four principal configurations with all four sites excited (there exist also configurations with two sites excited as for N = 1, in accordance with Remark 8). In agreement with Theorem 3, these four configurations are coded as + + ++ (when all sites excited in phase), + − +− (when all sites are out of phase with their immediate neighbors), + + −−, and + − −+. As can be seen in the figure and as is known from the Hamiltonian limit of γ = 0, the only spectrally stable among these configurations is the out of phase configuration + − +− with three imaginary eigenvalue pairs, while the in-phase configuration + + ++ is the most unstable with three real pairs. Upon continuation over γ, for both of these configurations, two of the eigenvalue pairs move towards zero (which they reach as γ → 1), while one remains real for the in-phase, and imaginary for the out of phase. Remarkably, pairwise these configurations collide and disappear in the limit of γ → 1. More specifically, the in-phase configuration + + ++ collides with the configuration + − −+. Similarly the out-of-phase configuration + − +− collides with the configuration + + −−. It should also be noted that we weren't able to continue any asymmetric mixed phase configurations (such as e.g. + + +−, + − −−, − + −−, − − +− or their opposite parity variants) past the Hamiltonian limit of γ = 0.

Conclusion
In this paper, we examined two distinct scenaria for PT-symmetric dynamical lattices. In the first one, N -site PT-symmetric chains were considered as a finite-dimensional dynamical system. In the second one, we have considered the embedding of the finite PT-symmetric system as a defect in an infinite dNLS lattice. In both cases, we have examined the linear problem, explicitly computing the corresponding eigenvalues and identifying the strengths of γ beyond which instabilities (and the phase transition breaking the PT-symmetry) arise due to real eigenvalues. We have also considered the nonlinear states when they emerge from the linear limit, as well as when they arise from a highly nonlinear limit under suitable rescaling (analogous to the anti-continuum limit of the standard dNLS lattice). In that case, we argued about the disparity of the branch counts in these two limits (for general N ), which suggests the existence of a number of bifurcations, such as saddle-center ones, at intermediate values of the corresponding parameter.
In the case of the infinite PT-symmetric PT-dNLS equation, i du n dt = u n+1 − 2u n + u n−1 + iγ(−1) n u n + |u n | 2 u n , n ∈ Z, we note that the phase transition threshold is now set at γ N →∞ = 0. Therefore, for any γ = 0, the linear PT-dNLS equation is unstable with a complex-valued continuous spectrum. Nevertheless, we can still obtain existence of stationary localized states (discrete solitons) in the large-amplitude limit for γ ∈ (−1, 1) for any of the configurations described in Theorems 3 and 5. Moreover, the discrete soliton can be centered at any site n 0 ∈ Z because of the shift invariance of the PT-dNLS equation (56).
There are other numerous directions in which one can envision generalizations of the present study. In the present work, we considered the case where there is a single parameter γ, for each of the N pairs of sites with gain and loss. However, it is also relevant to generalize such considerations to the case of many independent parameters for such sites [34]. On the other hand, one can consider generalizations of the present setting that aim towards the case of higher dimensionality. Arguably, the simplest such generalization concerns the setting of two one-dimensional coupled (across each of their sites) chains in the form of a railroad track as in [30] and the consideration of multi-site excitations therein. However, the genuinely higher dimensional problem and the examination of generalizations of plaquette configurations [21], whereby the potential of vortices exists in the Hamiltonian limit is of particular interest in its own right. These themes will await further consideration.