Discrete solitons and vortices on two-dimensional lattices of $\mathcal{PT}$-symmetric couplers

We introduce a 2D network built of $\mathcal{PT}$-symmetric dimers with on-site cubic nonlinearity, the gain and loss elements of the dimers being linked by parallel square-shaped lattices. The system may be realized as a set of $\mathcal{PT}$-symmetric dual-core waveguides embedded into a photonic crystal. The system supports $\mathcal{PT}$-symmetric and antisymmetric fundamental solitons (FSs) and on-site-centered solitary vortices (OnVs). Stability of these discrete solitons is the central topic of the consideration. Their stability regions in the underlying parameter space are identified through the computation of stability eigenvalues, and verified by direct simulations. Symmetric FSs represent the system's ground state, being stable at lowest values of the power, while anti-symmetric FSs and OnVs are stable at higher powers. Symmetric OnVs, which are also stable at lower powers, are remarkably robust modes: on the contrary to other PT-symmetric states, unstable OnVs do not blow up, but spontaneously rebuild themselves into stable FSs.


I. INTRODUCTION
Dynamics of light fields in dissipative media keeps drawing much interest in nonlinear optics [1,2]. A special kind of such systems with exactly balanced spatially separated gain and loss realize the general concept of the parity-time (PT )-symmetry [3,4]. PT -symmetric systems are represented by Hamiltonians containing a complex potential, V (r), which is subject to the symmetry constraint, V (r) = V * (−r). In optics this condition can be realized by guiding light through a medium with a complex refractive index satisfying condition n(r) = n * (−r) [5][6][7][8][9][10][11][12], hence its real and imaginary part, the latter one representing the combination of spatially distributed gain and loss, must be, respectively, even and odd functions of spatial coordinates.
Unlike generic dissipative systems, in which the balance between gain and loss selects parameters of isolated stable modes (attractors), PT -symmetric settings support continuous families of modes, sharing this property with conservative media, which was a motivation for many experimental [13,14] and theoretical studies, including the consideration of a great diversity of nonlinear PT -symmetric systems . In the latter case, solitons are a major subject of the studies.
An important variety of PT -symmetric systems in optics are discrete ones. The simplest among them, which was actually the first PT -symmetric system created experimentally [37], is realized as a pair of linearly coupled waveguides carrying gain and loss with equal strengths (a dimer [20,21,26,27,29,31]). An extension to various more complex discrete systems is provided by arrays of linearly coupled dimers [38][39][40][41][42][43][44][45][46][47][48]. In particular, robust one-dimensional (1D) discrete solitons [39] and scattering states [41] are supported by a chain of dimers with the cubic nonlinearity, in which active (gain-carrying) and passive (lossy) elements are coupled horizontally, while at each site the active and passive poles are coupled vertically. Accordingly, the chain is described by a system of linearly coupled discrete nonlinear Schrödinger (DNLS) equations with equal strengths of linear gain and loss in the two components. Stability, mobility and interactions between fundamental discrete solitons in this system were studied by means of numerical methods and analytical approximations.
There are fewer publications which reported PT -symmetric solitons in 2D settings [23,30,35]. The objective of the present work is to introduce a 2D network of horizontally coupled nonlinear dimers, which is a generalization of the 1D system from [39], and can also be realized in optics, using a photonic-crystal matrix into which nonlinear dual-core waveguides with the balanced gain and loss are embedded (see Fig. 1 below). We construct families of fundamental and vortical two-component discrete solitons in this 2D system and analyze their stability.
FIG. 1: The transverse cross-section of the waveguiding system built as a juxtaposition of networks of active (red) and passive (blue) elements, which are coupled into PT -symmetric dimers (the couplings are designated by short diagonal segments). The elements are connected into the networks by thin plates forming the double photonic crystal. The propagation axis (z) is directed perpendicular to the plane of the drawing.
The paper is structured as follows. The model is introduced in Section 2. Fundamental and vortex solitons are produced in Sections 3 and 4, respectively. The paper is concluded by Section 5.

II. MODEL
The network model, outlined above as the 2D extension of the chain of dimers introduced in [39], is based on a system of two linearly coupled normalized DNLS equations for amplitudes ψ m,n (z) and ϕ m,n (z) of the electromagnetic field propagating in the active and passive cores of dual-core waveguides (that represent the PT -symmetric dimers), which, as schematically shown in Fig. 1, are embedded into a quasi-2D photonic crystal, built of thin plates connecting the active and passive cores into respective sub-networks: Here (m, n) are discrete coordinates in the transverse plane, z is the propagation distance, the coefficient of the coupling between the active and passive cores is scaled to be 1, γ > 0 is the gain-loss coefficient, and C is coupling constant in the networks. Localized states produced by Eq. (1) are characterized by the total power, To estimate physical parameters of the setting, we assume that it is built of AlGaAs (the nonlinear coefficient of this material is n 2 ≃ 1.5 × 10 −13 cm 2 /W [49]), with core area 4 µm 2 of the lattice site, and refractive indices of the core and thin plates 3.263 and 3.256, respectively (the corresponding refractive-index difference can be produced by varying the share of aluminium in the composition of those parts). The standard coupled-mode theory [50,51] shows that, if we fix the distance between the adjacent active and passive lattice sites as 4 µm, the coupling coefficient between them is 3.6 mm −1 . Because we scale this coefficient to be 1 in Eq. (1), we can estimate the accordingly normalized values of the coupling and gain-loss coefficients, C and γ, and the unit of the total power, P 0 = 0.55 kW. For example, the lattice spacing 8 µm and the physical value of the gain-loss coefficient, 1.8 mm −1 , imply C = 0.3 and γ = 0.5. These values are realistic for the experimental realization, and, on the other hand, their normalized values give rise to nontrivial results, as shown below.
Stationary modes with real propagation constant k are looked for as solutions to Eq.
(1) as usual, {ψ m,n (z), ϕ m,n (z)} = e ikz {u m,n , v m,n }, where stationary amplitudes obey equations in the form of (k + iγ) u m,n = C 2 (u m,n+1 + u m,n−1 + u m−1,n + u m+1,n − 4u m,n ) + |u m,n | 2 u m,n + v m,n , Because our study is focused on localized modes, we apply zero boundary conditions at infinity (in fact, at edges of the numerical domain).
In the numerical simulations, it is convenient to combine u m,n and v m,n into a single "long" vector, U , of length 2N 2 , where N × N is the size of the numerical domain. Thus, Eq. (3) is rewritten in the vectorial form: where diag(|U | 2 ) means a diagonal matrix with the respective elements, the other matrices beingL = D −I −I D andΓ = γ I 0 0 −I , where I and D are, respectively, the N 2 × N 2 unit and linear-coupling matrices, the latter one The stability of the stationary modes is investigated in a numerical form, computing eigenvalues for small perturbations and verifying the results by direct simulations. Linearized equations for perturbation eigenmodes are produced by writing the perturbed solution, in terms of the "long vector", U , as {ψ m,n (z), ϕ m,n (z)} = e ikz (U + αe −iλz + β * e iλ * z ), where α and β are the respective vectors combining amplitudes of the eigenmodes, and * stands for the complex conjugation: The underlying solution, U , is stable if all eigenvalues λ are real.

A. Symmetric modes
As suggested by the analysis of the PT -symmetric coupler [19,22], the parity-time symmetry of Eqs. (1) remains unbroken at γ < 1, in which case the system of two coupled equations (3) can be reduced to a single one, by substitution v m,n = e iχ± u m,n , χ + = arcsin γ, χ − = π − arcsin γ, with ± in Eq. (6) corresponding to χ ± in Eq. (7). Because in the limit of γ = 0 relation (7) with χ + and χ − represents, respectively, symmetric and antisymmetric states in the two-component conservative system, these solutions may be naturally named PT -symmetric and antisymmetric ones [19]. Discrete 2D fundamental solitons (FSs) in the conservative version of Eqs. (1) and (3) with γ = 0 were previously studied in the context of the spontaneous formation of asymmetric solitons, with u m,n = v m,n , from the symmetric ones [52]. The transition to asymmetric solitons is not relevant for PT -symmetric systems, where the spontaneous symmetry breaking causes blow-up, rather than emergence of asymmetric modes [19,22]) (recently, it was demonstrated that asymmetric solitons may exist in a specially designed 1D PT -symmetric model [34]). However, the instability may lead, instead of the blow-up, to transformation into a stable soliton of another type, if different solitons species coexist in the system. In particular, it is demonstrated below that unstable PT -symmetric solitary vortices do not blow up, but rather rebuild themselves into stable symmetric FSs.
Equation (6) is the usual stationary form of the 2D DNLS equation, which gives rise to discrete FSs [53] and solitary vortices [54] whose total power (2) grows with the increase of the effective propagation constant, so that these modes satisfy the Vakhitov-Kolokolov criterion [55,56], dP/dk > 0, which is a necessary, but not sufficient, condition for their stability. Therefore, for given k, term ∓ 1 − γ 2 in Eq. (8) makes the total power of the PT -symmetric FS smaller than that of its antisymmetric counterpart. The latter fact, in turn, implies that the symmetric FS may represent the system's ground state (it is demonstrated below that this conjecture is true), and it should be more relevant to the experiment, as the creation of the mode with a smaller power is, obviously, easier.
Stationary solutions for symmetric FSs were produced, first, by means of the imaginary-time-propagation method (ITM) [57], applied to Eq. (1). For these computations, the set of control parameters includes the gain-loss and coupling strengths, (γ, C), and the total power of the soliton, P . It was checked that precisely the same family of symmetric solutions is generated by the Newton's method applied to stationary equations (4). The fact that these solutions are produced by the ITM method strongly suggests that they indeed represent the system's ground state [58]. It is also relevant to stress that Eq. (4) produces exactly the same solutions as the reduced equation (6), although reduction (7) was not applied to Eq. (4). A typical example of a stable PT -symmetric FS is displayed in Fig. 2.
For γ = 0, the results corroborate well-known properties of the 2D FSs in the DNLS equation [53]. In particular, symmetric on-site-centered FSs (their off-site-centered counterparts are completely unstable, as usual [53], and are not considered here) exist above a power threshold, Because stationary equations (3) for PT -symmetric modes reduce to the single equation (6), it is obvious that P (1) th does not depend on γ. Further, if the P (k) dependence for the FS in the single-component DNLS equation with γ = 0 is P 1 (k), which is well known in the numerical form [53], then an exact result for the two-component PT -symmetric and antisymmetric FSs, following from Eq. (6), is Numerical solutions for both types of the solitons completely agree with Eq. (10) (not shown here in detail). In fact, for both the PT -symmetric and antisymmetric FSs the P (k) dependence can be rather accurately fitted to an empiric relation An essential finding is that the symmetric FS becomes unstable above a higher threshold value of the total power, P = P (2) th . This stability boundary may be understood as the one at which the symmetric soliton is destabilized by the spontaneous symmetry breaking, modified by the presence of γ > 0 in comparison with the similar effect in the conservative system with γ = 0 [52]. Because asymmetric solitons cannot exist in the system with the balanced gain and loss, this symmetry breaking always leads to blow-up of the FS in the active component [ψ m,n in Eq. (1)]. The blow-up is quite similar to that shown below in Fig. 7(b) for an antisymmetric vortex.  The most essential results of the numerical analysis of the PT -symmetric FSs are summarized in Fig. 3, in the form of stability diagrams in parameter planes (γ, C), (γ, P ), and (C, P ). The horizontal and diagonal cutoffs in the (γ, C), (γ, P ) and (C, P ) planes exactly represent the threshold power, as given by Eq. (9). A natural trend is shrinkage of the stability area with the increase of γ, and its full disappearance close to γ = 1. It is natural as well that the stability area expands with the increase of C at fixed P , because the large coupling constant makes the modes, effectively, "less nonlinear", for the given total power, and this, in turn, suppresses the trend to the nonlinearity-driven spontaneous symmetry breaking which destabilizes the soliton.

B. Antisymmetric solitons
It has been found that the system supports stable and unstable antisymmetric FSs too [see a typical examples of a stable one in Fig. 2(b)]. The antisymmetric FSs cannot be produced by the ITM, which suggests that, as argued above, they do not represent the system's ground state. Nevertheless, they have been generated by the Newton's method applied to Eq. (4). For this reason, the propagation constant, k, rather than the total power, P , plays the role of the control parameter for the antisymmetric solitons. Their stability was identified through a numerical solution of Eq. (5), and verified by means of direct simulations of Eq. (1). Results of the stability analysis are collected in Fig. 4, in which panels (a-c) display the stability area in the (k, C) and (γ, C) panels. In addition, the stability area is mapped into the (C, P ) plane in plot 4(d), using the P (k) dependence [see Eq. (11)]. In the latter panel, the boundary of the white non-existence area is given by Eq. (9), i.e., it is exactly the same as in Fig. 3(c1) and 3(c2), once the existence of both symmetric and antisymmetric FSs is determined by the same equation, Eq. (6).
We have found that, at γ < 0.95, the stability area shown in panel 4(d) tends to extend to P → ∞ [or, in other words, to k → ∞ in panel 4(a)], up to the accuracy limits of the numerical calculations. This conclusion is in drastic contrast with the situation shown for the symmetric FSs in Fig. 3, where the stability is always bounded from above, in terms of P . The explanation for the stability of the antisymmetric solitons at large P is straightforward: unlike their symmetric counterparts, they do not undergo the destabilization via the symmetry breaking at large P . However, at γ > 0.95, the stability area becomes finite in terms of k, see Fig. 4(c) for γ = 0.999. Although the area is not vanishingly small at γ so close to 1, it exactly collapses to nil at γ = 1, as it should, as the PT symmetry gets completely broken at this point.
A general conclusion about the antisymmetric FSs (which pertains equally well to antisymmetric vortex solitons, see below) is that, although they tend to be stable at large values of the power, where their symmetric counterparts are unstable, they are less suitable for the experimental realization just because they demand large powers for the stability, while the symmetric FS, which represents the ground state of the system, may be stable at lower levels of the power.

IV. VORTEX SOLITON
To the best of our knowledge, the symmetry breaking of discrete vortex solitons in 2D two-component linearlycoupled lattices was not investigated even in the conservative model, with γ = 0. We have found symmetric, antisymmetric, and asymmetric solitary vortices, and the symmetry-breaking transition, in the model based on Eqs. (1) and (3) with γ = 0. These findings will be reported elsewhere, while here we focus on PT -symmetric and antisymmetric vortices in the system with γ > 0 (in fact, it is shown below that, for symmetric vortices, the stability area only weakly changes with the variation of γ, hence the present results give an idea of the stability of symmetric vortices in the conservative system too). Because the ITM cannot converge to the definitely non-ground-state vortical modes, the solutions were found by means of the Newtons's method applied to Eq. (4), hence the control parameter of the solution families is the propagation constant, k, rather than total power P .
It is well known that two types of vortices are possible in discrete systems, viz., on-and off-site-centered ones [53] (alias "rhombuses" and "squares"). Both can be produced by the Newton's method applied to Eq. (3), starting from the following inputs: where U and V are real constants.
In agreement with the known result for the usual 2D DNLS equation with γ = 0 [54], we have found that the off-site vortices are unstable at practically all values C > 0 (they may be stable at extremely small C), therefore in what follows below we consider only the modes of the on-site types.
As well as in the case of FSs, stationary solutions for on-site-centered vortices (OnVs) are subject to reduction (7), which leads to the single equation (6), hence the OnV may also be classified into PT symmetric or antisymmetric varieties. Examples of stable symmetric and antisymmetric ones are displayed in Figs. 5 and 6.
A unique dynamical property of the symmetric OnVs is that, when they are unstable, they do not blow up, unlike the FSs, and unlike antisymmetric OnVs, which do blow up if being unstable, as shown in Fig. 7(b). Instead, as seen in Fig. 7(a), the instability transforms symmetric vortices into stable symmetric FSs. Thus, the symmetric vortices demonstrate remarkable robustness as self-trapped modes of the PT -symmetric 2D nonlinear system. The most important results for the symmetric and antisymmetric vortices, namely, their existence and stability areas in planes of control parameters (k, C) and (C, P ), are summarized in Fig. 8. For the symmetric OnVs, it is not excluded that, above the triangular existence regions in Fig. 8(a1), and above the existence region in Fig. 8(a2), there exist very strongly unstable vortices which could not be produced by the numerical method.
Numerical data demonstrate that dependence P (k) for the symmetric and antisymmetric vortices alike is almost exactly given by where P FS (k) is the respective dependence for the fundamental solitons, see Eq. (11). This relation is naturally explained by the fact that the on-site-centered vortex may be constructed, essentially, as a set of four weakly interacting FSs, with phase shifts π/2 between adjacent ones [54]. Relation (14) also helps to map the existence and stability regions for the OnVs in the (C, P ) parameter plane in Figs. 8(a2) and (c2). Actually, the fact that the symmetric OnVs (at least those which could be produced by the numerical method) exist only at relatively small values of P [at P < 4 in Fig. 8(a2)] explains that, as stressed above, their instability leads not to the blow-up, but rather to the spontaneous transformation into stable symmetric FSs. It is worthy to note too that the stability region for the symmetric vortices in Fig. 8(a2) is located above the instability area (at larger P ), while for the symmetric FSs the situation was opposite, cf. Figs. 3(c1) and 3(c2). In this respect, the symmetric vortices resemble antisymmetric FSs, as suggested by the comparison with Fig. 4(d). However, the symmetric OnVs share with symmetric FSs the trend to exist and be stable at values of the total power essentially lower that those at which antisymmetric counterparts are stable, hence the symmetric vortices are more appropriate for the experimental creation.
The stability areas for the antisymmetric OnVs, displayed in Figs. 8(b1)-8(c2), resemble their counterparts shown in Fig. 4 for antisymmetric FSs. In particular, the stability zone tends to extend to k → ∞, i.e., P → ∞, except for the case when the gain-loss coefficient, γ, is very close to its critical value, 1, see Fig. 8(c1) for γ = 0.999.
The overall dependence of the stability of the vortices on the gain-loss parameter γ is characterized by total areas of the stability regions in the (k, C) plane [see Figs. 8(a1) and 8(b1)-(c1)], which are shown as functions of γ in Fig.  9. The drastic difference between these dependences for the symmetric and antisymmetric vortices is obvious.
It was argued above that the fundamental PT -symmetric FS are likely candidates to the role of the system's ground state. To further verify this conjecture, in Table 1, we compare values of the chemical potential, µ, for a typical set of the control parameters, (C, γ) = (0.05, 0.1) and P = 1.9. All the species of discrete solitons considered above exist and are stable at this point (symmetric and antisymmetric FSs and OnVs). The table demonstrates that the symmetric FS has a deep minimum of µ, thus definitely representing the ground state, the three other species being excited states. It is worthy to note too that the lowest excited state is the PT -symmetric vortex, but not the antisymmetric FS. Similar relations between chemical potentials of the four soliton species have been found at other values of the parameters. The objective of this work is to introduce the 2D network of PT -symmetric dimers with the cubic nonlinearity, which support stable PT -symmetric and antisymmetric fundamental and vortical discrete solitons. While the shape of these solitons can be reduced to that of their counterparts in the conservative two-component system (although discrete vortices were not systematically studied even in the conservative system with two linearly coupled components), the study of the stability is crucially important for the consideration of the PT -symmetric system. It has been found that the symmetric FS (fundamental soliton) represents the ground state of the system, and may be stable at lowest values of the total power, thus being the most promising target for the experimental realization. On the contrary, anti-symmetric FS and OnV (on-site-centered vortex) tend to be stable at higher powers. The symmetric OnV represents the first metastable mode above the FS ground state. It also tends to be stable at relatively low powers, and demonstrates remarkable robustness: unlike all the other discrete-soliton species found in this system, unstable symmetric OnVs do not suffer the blow-up, but rather spontaneously transform into stable FSs.
The above analysis did not address mobility of the solitons, but, unlike the 1D version of the system, mobility of 2D discrete solitons is not plausible in the case of the cubic nonlinearity [53,59,60]. In this connection, it may be interesting to study 2D networks of PT -symmetric dimers with the quadratic nonlinearity, cf. Refs. [62][63][64], and also Ref. [61], where the mobility was demonstrated for 2D discrete solitons in the conservative systems with the quadratic nonlinearity. Another promising extension may deal with 2D systems built of elements which include nonlinear PT -symmetric gain and loss terms [17,21].