One- and two-dimensional solitons in PT-symmetric systems emulating the spin-orbit coupling

We introduce a two-dimensional (2D) system, which can be implemented in dual-core planar optical couplers with the Kerr nonlinearity in its cores, making it possible to blend effects of the PT symmetry, represented by the balanced linear gain and loss in the two cores, and spin-orbit coupling (SOC), emulated by a spatially biased coupling between the cores. Families of 1D and 2D solitons and their stability boundaries are identified. In the 1D setting, the addition of the SOC terms leads, at first, to shrinkage of the stability area for PT-symmetric solitons, which is followed by its rapid expansion. 2D solitons have their stability region too, in spite of the simultaneous action of two major destabilizing factors, viz., the collapse driven by the Kerr nonlinearity, and a trend towards spontaneous breakup of the gain-loss balance. In the limit of the SOC terms dominating over the intrinsic diffraction, the 1D system gives rise to a new model for gap solitons, which admits exact analytical solutions.


I. INTRODUCTION
The recent progress in the experimental and theoretical work with engineered optical media has made it possible to emulate, by means of the optical-beam propagation, a wide range of physical effects which were originally predicted or experimentally discovered in other areas of physics. In many cases, the optical emulation offers a possibility to study the effect in question in a pure form, which is often too difficult in the original setting. In other cases, specially designed optical configurations open a way to demonstrate realizations of the effects in forms which are impossible in the source systems. Belonging to this class of optically emulated phenomena are the parity-time (PT ) symmetry [1], based on the paraxial beam propagation in optical media with symmetrically placed gain and loss elements, as predicted theoretically [2] and realized experimentally [3]; Anderson localization in random photonic media [4]; and the photonic emulation of topological insulators [5] and graphene [6]. A recent addition to this topic is a proposal to implement the mechanism of the spin-orbit coupling (SOC), previously elaborated in terms of Bose-Einstein condensates (BECs) [7], in dual-core optical waveguides [8,9]. Basic schemes of "grafting" new physical effects to optics make use of linear beam-propagation regimes. However, the ubiquitous Kerr nonlinearity of optical materials, as well as other types of the optical nonlinear response, suggest one to consider implementation of the new effects in a nonlinear form. A well-known example is provided by PTsymmetric solitons, which have been studied in detail theoretically [10] and created in the experiment [11]. The optical emulation of the SOC was also developed in the context of the nonlinear propagation and formation of solitons [9].
As a further development in these directions, it may be interesting to design optical settings which combine the above-mentioned effects, that would be difficult to achieve in their original realizations. The objective of the present work is to elaborate one-and two-dimensional (1D and 2D) optical systems, based on dual-core waveguides, which make it possible to blend the SOC and PT -symmetry mechanisms. The use of couplers for this purpose is natural, as they provide photonic platforms for the emulation of both the SOC and PT symmetry in 1D [8,12] and 2D [9,13] geometries alike. As an additional result of the analysis, a new 1D conservative model for two-component gap solitons is produced, which admits exact solutions, and demonstrates a nontrivial internal stability boundary in the soliton family. New dynamical features are revealed by the analysis of the combined models, such as a nonmonotonous dependence of the stability area of 1D PT -symmetric solitons on the SOC strength, δ: the area originally shrinks but then strongly expands with the increase of δ. Another noteworthy finding is that the SOC terms stabilize 2D PT -symmetric solitons under the action of the cubic self-attraction, in spite of the simultaneous presence of two mechanisms driving the catastrophic instability in the 2D system: the possibility of the breakup of the PT symmetry, i.e., failure of the gain-loss balance [2,10], and the onset of the critical collapse induced by the Kerr nonlinearity [14] (the latter may also take place in the presence of SOC [15]). The latter result implies extension of the SOC-induced stabilization of 2D solitons (semi-vortices and mixed modes) in the free space with the cubic self-attraction, that was demonstrated recently [16].
The paper is organized as follows. The basic model is introduced in Section II. Results for 1D solitons and their stability are collected in Section III. The new 1D model for gap solitons and its PT -symmetric version are presented in Section IV. Results for 2D solitons are reported in Section V, and the paper is concluded by Section VI.

II. THE MODEL
The starting point is a model for the dual-core planar optical waveguide which governs the spatiotemporal evolution of local amplitudes of the electromagnetic waves in two cores, U 1 (x, t, z) and U 2 (x, t, z), under the action of the anomalous group-velocity dispersion (GVD) and Kerr nonlinearity: In these coupled 2D nonlinear Schrödinger equations (NLSEs) z is the propagation distance, x the transverse coordinate, t the reduced temporal variable, the second x derivatives account for the paraxial diffraction, while the coupling coefficient, together with the GVD and Kerr coefficients, are scaled to be 1. The 1D reduction of Eqs. (1) and (2) amounts to the commonly known 1D model for dual-core optical fibers or the spatial propagation in the twin-core planar waveguides [17]. By itself, the system of Eqs. (1) and (2) gives rise solely to unstable 2D solitons, due to the occurrence of the collapse in the same setting. The stabilization of the spatiotemporal solitons against the collapse may be provided by the optical counterpart of the SOC, which amounts to taking into account the temporal dispersion of the coupling coefficient [9], represented by terms δ (U 2 ) t and δ (U 1 ) t to be added to Eqs. (1) and (2), respectively, with a real dispersion coefficient δ [18]. The accordingly modified coupled equations, which also include the balanced gain and loss terms, with real coefficient Γ > 0 [12], that introduce the PT symmetry, take the form of All the ingredients of the present system, including the separated gain and loss [3], can be realized experimentally in optics. Nevertheless, the model based on Eqs. (3) and (4) turns out to be irrelevant, because the linearized version of the equations for excitations, in the usual form of U 1,2 ∼ exp (ikz − iωt + iqx), gives rise to an unstable dispersion relation between the propagation constant, k, real temporal frequency, ω, and real transverse wavenumber, k: Indeed, it follows from Eq. (5) that, in the presence of the gain and loss, the zero solution is unstable, which is accounted for by an imaginary part of k, against perturbations with (1 + δ · ω) 2 < Γ 2 . In other words, the PT symmetry is always broken by the SOC terms in this system, while it gives rise to stable 2D solutions at Γ = 0 [9].
A system which maintains the PT symmetry in the presence of the SOC is based on the following coupled NLSEs: This is a model of a planar coupler, in which the temporal dispersion of the coupling coefficient is disregarded, while terms ∓δ (U 2,1 ) x account for "skewness" of the coupling in the transverse direction, assuming that the layer between the guiding cores has an oblique structure. Roughly speaking, the latter means that light tunnels between point x in the first (second) core and point x + δ (respectively, x − δ) in the mate core. The dispersion relation for the linearized version of Eqs. (6) and (7) is Unlike Eq. (5), it demonstrates that the PT symmetry holds at Γ ≤ 1 (in fact, under the condition that the gain-loss coefficient is smaller than the inter-core coupling constant, which is scaled to be 1).
The system may produce solitons in its bandgap, i.e., at values of the propagation constant, k, which cannot be obtained from Eq. (8). Simple analysis demonstrates that the dispersion relation (8) gives rise to a semi-infinite bandgap, in which solitons are expected to exist:

III. SOLITONS IN THE 1D SYSTEM
A family of stationary soliton solutions of the 1D version of Eqs. (6) and (7), in which the temporal dependence is dropped, are looked for as with complex functions u 1,2 (x) satisfying equations with the prime standing for d/dx. The PT symmetry amounts to the cross-symmetry constraint for the two components, which all the soliton solutions obey: (the asterisk stands for the complex conjugation). The solutions are characterized by their norm, i.e., the total power, in terms of the underlying optics model: In the absence of the SOC terms (δ = 0), exact PT -symmetric soliton solutions to Eqs. (11) and (12) are well known [12], where φ is a constant phase. The exact stability condition for these solutions is known too [12]. Combined with the obvious existence condition, k > √ 1 − Γ 2 , it produces the following interval filled by the stable PT -symmetric solitons at δ = 0 (in the absence of the SOC terms): It is more relevant to write this in terms of the norm of soliton (15), In the presence of the SOC (δ > 0), we have built stationary localized solutions of the system by means of the well-known imaginary-time-integration method [19]. A notable SOC effect is the splitting between peaks of the two components. A simple perturbative analysis of Eqs. (11) and (12) demonstrates that, for small δ, the distance between the peaks is    (17) is virtually exact at δ ≤ 0.4. It is worthy to note that the profiles of |u 1 (x)| and |u 2 (x)| are nearly identical for different values of the gain-loss coefficient, Γ = 0.3 and 0.8, with the same δ = 2, in Figs. 1(b,c). This observation is in qualitative agreement with the fact that, for given N , the profiles of exact solitons (15) do not depend on Γ either.
An obvious effect of the increase of δ is a transition from the smooth soliton profile to a multi-lobe structure, as seen from the comparison of panels (a) and (b,c) in Fig. 1. This fact can be explained by the inversion of the dispersion relation (8), to express q in terms of k, at ω = 0. The result of a simple algebra is that, precisely at k > √ 1 − Γ 2 , i.e., when the exact soliton (15) exists for δ = 0, wavenumber q becomes complex, which implies the wavy profile of the soliton's tails, for The numerical results demonstrate that the transition occurs at δ ≈ 0.6 for Γ = 0.3, and at δ ≈ 0.7 for Γ = 0.8, which is consistent with Eq. (18) in which numerically obtained values of k are substituted. The stable soliton family is further characterized by Fig. 2(a), which shows the largest value of the field |u 1 (x)| as a function of N , for δ = 2 and Γ = 3. Naturally, |u 1 (x = 0)| increases monotonously with N , up to a critical value, N c , at which the destabilization of the solitons happens via the breakup of the PT symmetry (the solitons exist at N > N c as unstable solutions, which cannot be found by means of the imaginary-integration method). In the conservative system with Γ = 0, stable asymmetric solitons, with u 1 (x) = u 2 (−x), may exist at N > N c , but in the presence of Γ > 0 this is impossible, as asymmetric solitons would not maintain the balance between the gain and loss. Figure 2(d) shows an example of a stable asymmetric solution found at Γ = 0 for N = 5 and δ = 2 In fact, N c is the most important characteristic of the soliton families. It is shown, as a function of the SOC strength, δ, for fixed Γ = 0, 0.4, and 0.8, in Fig. 2(b). In particular, N c is given by the exact expression (16) for δ = 0. The initial decrease of N c , which is observed in Fig. 2(b) at relatively small δ, can be easily understood: as shown by Eq. (17), the two components get separated by distance ∆x ≈ δ, hence the effective linear coupling between the components weakens with the increase of δ, making a smaller strength of the Kerr nonlinearity (measured by the norm) sufficient to initiate the symmetry breaking between the two components. Obviously, this effect should scale as δ 2 , which is consistent with Fig. 2(b) at δ small enough. However, a nontrivial finding is that, at δ exceeding values corresponding to the minimum of N c in Fig. 2(b), N c features rapid increase with δ, which becomes asymptotically linear at large values of δ. This property is explained in the next section which addresses the limit form of system (6), (7) with the SOC terms dominating over the paraxial diffraction.
Further, Fig. 2(c) shows N c as a function of Γ for fixed SOC strengths, δ = 0, 1, and 2, the dashed curve representing the exact result given by Eq. (15) for δ = 0. This figure shows too that the SOC terms tend to essentially expand the stability area for the PT -symmetric solitons, by increasing N c . However, even if N c keeps a nonzero value up to Γ = 1 at δ > 0, there is no stability region at Γ > 1, as the zero solution (the background of the solitons) is unstable in the latter case, as follows from Eq. (8).

IV. THE LIMIT CASE OF NEGLIGIBLE INTRINSIC DIFFRACTION
The results displayed above in Fig. 2(b) for large δ suggest to consider the limit case of the system in which the SOC terms dominate over the paraxial diffraction. Dropping the second derivatives in the 1D version of Eqs. (3) and (4), one thus arrives at the simplified system, whose dispersion relation is cf. Eq. (15). It gives rise to a finite bandgap, k 2 < 1 − Γ 2 [unlike the semi-infinite bandgap given by Eq. (9)], hence the corresponding localized states may be considered as gap solitons. Actually, the system of Eqs. (19) and (20) with Γ = 0 is a new conservative model generating gap solitons, therefore it makes sense to consider its solutions, along with its PT -symmetric version, corresponding to 0 < Γ < 1. An essential difference from the standard gap solitons generated by the Bragg-grating model [20] is the separation between peaks of the two components, and, on the other hand, gap-soliton solutions are real in the present model. Looking for stationary solutions to Eqs. (19) and (20) as per Eq. (10), in the absence of the gain and loss, Γ = 0, functions u 1 (x) and u 2 (x) satisfy a system of real equations whose solutions obey the cross-symmetry constraint, cf. Eq. (13): u 1 (−x) = u 2 (x). In fact, δ may be easily absorbed into rescaled coordinate x, therefore numerical results are presented below for δ = 1. The total norm (14), defined according to the original coordinate, scales as δ, which explains the asymptotically linear N c (δ) dependence in Fig.  2

(b).
It is possible to find exact soliton solutions to Eqs. (22) and (23), noting that the evolution of u 1 and u 2 along x conserves the corresponding formal Hamiltonian: Next, it is convenient to represent solutions in the "polar" form, For soliton solutions which vanish at x → ±∞, one should set h = 0, hence Eq. (24) makes it possible to eliminate u 2 in favor of θ, after substituting expressions (25): Finally, combining two equations (22) and (23) and Eq. (26), it is easy to derive a single equation for θ(x): A solution to Eq. (27), which generates a soliton after the substitution in Eq. (26), exists for k 2 < 1: The total norm of the solitons calculated as per Eqs. (14) and (25)- (28), can be written as Thus, with the increase of k from −1 to +1, N (k) monotonously grows from N (k = −1) = 0 to N (k = +1) = 2 √ 2πδ. The analytical solution makes it possible to find the distance between peaks of the two components, cf. Eq. (17). In the general case, the analytical expression for ∆x is cumbersome. It takes a relatively simple form for k = 0, which corresponds to a stable gap soliton (see below): Figure 3(a) displays an example of the exact solution for k = −0.5 and δ = 1, whose total norm is N = 3.35, in agreement with Eq. (29). Further, Fig. 3(b) shows a result of the test of the stability of this gap soliton, produced by direct simulations of Eqs. (19) and (20) with Γ = 0 and small random perturbations added to the initial conditions. Figure 3(b) clearly shows that the gap soliton is stable. Figure 3(c) displays another example of the exact soliton,for k = 0.6 and δ = 1, with the total norm N = 5.82, which also agrees with Eq. (29). In this case, the simulations, the result of which is shown in the top plot of Fig.  3(d) at z = 800, demonstrate weak instability of the soliton, which leads to generation of an undulating tail. The simulations for other values of k, that are displayed too in Fig. 3(d), suggest that the intrinsic boundary between stable and unstable gap solitons is located at k ≈ 0.5, cf. a qualitatively similar results for the gap solitons in the Bragg-grating model [21].
The PT -symmetric version of Eqs. (19) and (20) with Γ > 0 was solved numerically, which also produced solitons. Two examples, for δ = 1 and Γ = 0.5, but different propagation constants, k = −0.5 and k = 0.6, with the norms, respectively, N = 3. 20 and 6.31,are displayed in Figs. 4(a) and (c). Simulations of the perturbed evolution of the former soliton demonstrate its stability in Fig. 4(b). On the other hand, a set of results of the simulations of the solitons with several values of k, presented in Fig. 4(d) at z = 900, make it evident that the latter soliton is unstable, the stability boundary being located at k ≈ 0.2.

V. TWO-DIMENSIONAL SOLITONS
Numerical solution of the full 2D system of Eqs. (3) and (4) produces stable solitons too, which are characterized by the respective norm, Examples of the 2D solitons are displayed in Figs. 5 and 6, which, in particular, demonstrate the mirror symmetry of |U 1 (x, t)| and |U 2 (x, t)|. Note that, similar to what was observed above for 1D solitons, the increase of the SOC strength, δ, generates a complex multi-lobe shape of the 2D solitons. As shown in Fig. 7, two different critical values of the norm can be identified for the 2D solitons. The upper one, N (upp) c , is the value at which, similar to what is found above for the 1D system, the PT -symmetric solitons are destabilized by the spontaneous symmetry breaking. Specific to the 2D setting is a lower limit, N (low) c , below which the solitons cannot self-trap. At N < N (low) c simulations demonstrate spreading of the wave fields. The existence of stable 2D solitons in the free space, under the action of the cubic-only attractive nonlinearity, is a nontrivial finding, as it was commonly believed that all solitons in such setting are destabilized by the (critical) collapse [14,15]. Recently [16], it was reported that the linear SOC terms may stabilize 2D solitons, as these terms break the specific scaling invariance of the 2D NLSEs, which accounts for the critical collapse in this case. The stability regions of the PT -symmetric 2D solitons, demonstrated in Fig. 7, are a still stronger result, as the solitons are able to stay stable despite the simultaneous presence of two major destabilization factors: the possibility of the collapse, and the trend to spontaneous breakup of the balance between the gain and loss in the two cores of the coupler. It is pertinent to mention that other types of the SOC terms, which are relevant to two-component BEC, do not give rise to N (low) c , the respective stability region being 0 < N < N (upp) c [16]. However, the PT -symmetry cannot be introduced in a physically relevant form in that setting. Naturally, Fig. 7(a) shows that the stability region found in the present system, N (low) c < N < N (upp) c , shrinks to nothing in the limit of δ → 0, when the 2D modes degenerate into the commonly known unstable Townes solitons [14]. On the other hand, a nontrivial finding is that the stability region attains its largest size at a finite value of the SOC strength, δ ∼ 1, as seen in Fig. 7. Similar to the 1D system [cf. Fig. 2(c)], the stability region keeps a (small) finite width at Γ = 1, completely disappearing at Γ > 1.

VI. CONCLUSION
The objective of this work is to introduce 1D and 2D optical systems which makes it possible to emulate the PT symmetry in combination with effects of the SOC (spin-orbit coupling), thus creating a novel physical setting. The systems are based on dual-core planar optical couplers with the Kerr nonlinearity in their cores. The PT symmetry is represented by equal amounts of the linear gain and loss in the two cores, while the SOC is induced by the skew form of the coupling between the cores. Families of stable 1D and 2D PT -symmetric solitons have been identified by means of numerical and analytical methods. The size of the stability region of the 1D solitons nonmonotonously depends of the SOC strength, δ, originally shrinking and then rapidly expanding with the increase of δ. In the limit of the SOC terms dominating over the paraxial diffraction, the 1D system produces a new model for gap solitons, for which exact analytical solutions have been found. 2D solitons may be stable against the combined action of two major destabilization factors, viz., the critical collapse driven by the Kerr self-focusing, and the trend to spontaneous breakup of the balance between the gain and loss in the coupled cores.
The present analysis may be extended in other directions. In particular, it is interesting to consider the mobility of 1D and 2D solitons and collisions between them, which is a nontrivial problem in the presence of the SOC [16]. Dynamical regimes, such as Josephson-like oscillations of localized modes between the coupled cores, may be interesting