Robust PT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{P}}{\bf{T}}$$\end{document} symmetry of two-dimensional fundamental and vortex solitons supported by spatially modulated nonlinearity

The real spectrum of bound states produced by PT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{P}}{\bf{T}}$$\end{document}-symmetric Hamiltonians usually suffers breakup at a critical value of the strength of gain-loss terms, i.e., imaginary part of the complex potential. The breakup essentially impedes the use of PT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{P}}{\bf{T}}$$\end{document}-symmetric systems for various applications. On the other hand, it is known that the PT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{P}}{\bf{T}}$$\end{document} symmetry can be made unbreakable in a one-dimensional (1D) model with self-defocusing nonlinearity whose strength grows fast enough from the center to periphery. The model is nonlinearizable, i.e., it does not have a linear spectrum, while the (unbreakable) PT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{P}}{\bf{T}}$$\end{document} symmetry in it is defined by spectra of continuous families of nonlinear self-trapped states (solitons). Here we report results for a 2D nonlinearizable model whose PT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{P}}{\bf{T}}$$\end{document} symmetry remains unbroken for arbitrarily large values of the gain-loss coefficient. Further, we introduce an extended 2D model with the imaginary part of potential ~xy in the Cartesian coordinates. The latter model is not a PT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{P}}{\bf{T}}$$\end{document}-symmetric one, but it also supports continuous families of self-trapped states, thus suggesting an extension of the concept of the PT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{P}}{\bf{T}}$$\end{document} symmetry. For both models, universal analytical forms are found for nonlinearizable tails of the 2D modes, and full exact solutions are produced for particular solitons, including ones with the unbreakable PT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{P}}{\bf{T}}$$\end{document} symmetry, while generic soliton families are found in a numerical form. The PT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{P}}{\bf{T}}$$\end{document}-symmetric system gives rise to generic families of stable single- and double-peak 2D solitons (including higher-order radial states of the single-peak solitons), as well as families of stable vortex solitons with m = 1, 2, and 3. In the model with imaginary potential ~xy, families of single- and multi-peak solitons and vortices are stable if the imaginary potential is subject to spatial confinement. In an elliptically deformed version of the latter model, an exact solution is found for vortex solitons with m = 1.

While wave functions of quantum systems may be complex, spectra of their energy eigenvalues must be real, which is usually secured by restricting the underlying Hamiltonian to be Hermitian 1 . However, the condition of the reality of the energy spectrum does not necessarily imply that it is generated by an Hermitian Hamiltonian. Indeed, it is well known that non-Hermitian Hamiltonians obeying the parity-time (PT ) symmetry may also produce entirely real spectra [2][3][4][5][6][7] . In terms of the single-particle complex potential, the PT symmetry requires its real and imaginary parts to be even and odd functions of coordinates 2 : where γ > 0 0 and β ≥ 0 are constants, features the unbreakable or nearly unbreakable 2D PT symmetry, represented by several species of families of stable solitons: single-and double-peak ones, as well as 2D solitons with embedded integer vorticity (topological charge), = m 1, 2, 3. The second model, with is not, strictly speaking, a PT -symmetric one, but it is equally relevant for the realization in optics, and it shares basic manifestations of the PT symmetry, maintaining families of single-and multi-peak solitons [featuring up to five peaks, in accordance with the structure of W x y ( , )] and solitary vortices, also with = m 1, 2, 3. The latter result is a contribution to the general topic of constructing models more general than the PT -symmetric ones with similar properties(including the case of the partial PT symmetry 51 ), which has been addressed in various settings [52][53][54]56,[70][71][72][73][74][75] , see also review 47 .
In both models, universal analytical forms are obtained for tails of solitons, and full exact solutions are produced for particular species of single-peak solitons, with β = 0 in Eqs (3) and (4). In the former case, the existence of the exact solitons at arbitrarily large values of γ 0 in Eq. (3) explicitly demonstrates the unbreakability of the PT symmetry. In the latter case two different families of exact solutions are found, which, however, exist only for γ ≤ 2 0 in Eq. (4) with β = 0. In addition, an anisotropic version of the latter model gives rise to particular exact solutions for vortex solitons with topological charge = m 1. Generic soliton families with = m 0, 1, 2, 3, which include the exact single-peak solutions as particular ones, are constructed in a numerical form in both models, and their stability is investigated numerically-both through computation of eigenvalues for small perturbations and by means of direct simulations.

Results
the models and analytical solutions for solitons. The underlying equations. The 1D NLSE for the amplitude of the electromagnetic field, u(x, z), with the local strength of the self-defocusing nonlinearity, Σ x ( ), growing from = x 0 towards = ±∞ x faster than |x| (this condition is necessary for self-trapping imposed by the self-repulsion 64 ), which is capable to maintain bright solitons with unbreakable PT symmetry, is 60 Here z and x are scaled propagation coordinate and transverse coordinate, in terms of the planar optical waveguide. In work 60 , the analysis was presented for a steep 1D modulation profile, with σ ≥ 0, where coefficients equal to 1 may be fixed to these values by means of rescaling. The choice of this profile allows one to obtain a particular exact solution for solitons 64 . Of course, in a real physical medium the local strength of the nonlinearity, defined as per Eq. (6), cannot grow to infinitely large values at | | → ∞ x . However, in reality it is sufficient that it grows according to Eq. (6) to finite values, that correspond to |x| which is essentially larger than the width of the soliton created by this profile. The growth of Σ x ( ) may be safely aborted at still larger |x| 64 .
Further, the spatially-odd imaginary potential, which accounts for the PT -symmetric gain-loss profile (cf. Eq. (1)), was introduced in ref. 60 as with γ > 0 0 and β ≥ 0. In the case of the spatially uniform self-focusing cubic nonlinearity, the 1D imaginary potential in the form given by Eq. (7) was introduced in ref. 76 .
Here, we aim to define a 2D extension of the model, as the NLSE for the propagation of the electromagnetic field with amplitude u(x, y, z) in the bulk waveguide with transverse coordinates (x, y): where ≡ + r x y 2 2 is the radial coordinate, and the nonlinearity-modulation profile is chosen similar to its 1D counterpart (6): with σ ≥ 0. Further, we consider two different versions of the 2D imaginary potential. First, it is a PT -symmetric one given by Eq. (3). The other imaginary potential, defined as per Eq. (4), is not PT -symmetric, because the  transformation, → − − x y x y ( , ) ( , ), does not reverse the sign of W(x, y), in this case. However, in terms of the implementation in optics the gain-loss distribution corresponding to Eq. (4) is as relevant as that defined by Eq. (7), and, as mentioned above, properties of solitons in models which are akin to PT -symmetric ones is a subject of considerable interest.
Stationary states with a real propagation constant, k, are looked for as solutions to Eq. (8) in the form of with complex function U(x, y) satisfying the following equation: Asymptotic solutions. As mentioned above, Eqs (8) and (11) are nonlinearizable, i.e., they cannot be characterized by a linear spectrum. Indeed, straightforward analysis of Eq. (11) demonstrates that it may produce localized solutions (solitons), with tails decaying at → ∞ r according to an asymptotic expression which is determined by the full nonlinear equation, rather than by its linearization. For the PT -symmetric imaginary potential (3) with β = 0, it is provided that σ ≠ 0. In the case case of σ = 0, this asymptotic solution is replaced by www.nature.com/scientificreports www.nature.com/scientificreports/ Note that asymptotic solutions given by Eqs (12) and (13) exist at arbitrarily large γ 0 , suggesting the unbreakability of the PT symmetry in this case, as corroborated by exact solution (19) produced below.
The imaginary potential defined by Eq. (4) with β = 0 produces the following result: asympt 0 2 2 0 for σ ≠ 0, and if σ = 0, the result is On the contrary to the the above asymptotic solutions, given by Eqs (12) and (13), which are available for arbitrarily large γ 0 , their counterparts produced by Eqs (14) and (15) exist only at γ < 2 0 , i.e., if the gain-loss coefficient is not too large.
It is relevant to stress the universal character of all asymptotic approximations given by Eqs (12)(13)(14)(15): they depend solely on coefficients σ and γ 0 of the underlying model, and, unlike the commonly known asymptotic forms of solitons in usual systems, do not depend on the propagation constant, k. The single exception is presented by exact solution Eq. (18) given below, whose asymptotic form (actually coinciding with the exact soliton solution, in that case) explicitly depends on k, but this happens solely for specially chosen parameters given by Eq. (17). In the generic case, a dependence on k appears in the next-order correction to the shape of the asymptotic tail. In particular, the correction to the tails given by Eqs (12) and (13) are Furthermore, for more complex solutions, such as multi-peak solitons and solitary vortices, as well as for higher-order radial states of the single-peak solitons, which are produced below in the numerical form, the asymptotic expression at large r is exactly the same as given by Eqs (12)(13)(14)(15).
Exact solutions for single-peak solitons. Precisely at the above-mentioned critical value γ = 2 0 , the asymptotic solutions (14) and (15) vanish. However, in the special case, 0 the vanishing asymptotic solution given by Eq. (15) is replaced by a different one, which, as can be easily checked, is an exact solution to Eq. (11) (not just an asymptotic approximation valid at large r), It exists, as the continuous family, at all values of < − k 1. Further, Eq. (11) which includes the PT -symmetric imaginary potential Eq. (3), with β = 0, gives rise to an exact solution at a special value k x 0 ( ) of the propagation constant: which exists at all values of coefficients γ 0 and σ, except for σ = 0. In other words, at = k k x 0 ( ) the asymptotic approximation Eq. (12) is tantamount to the exact solution. This solution features the unbreakable PT symmetry, as it persists at arbitrarily large values of the gain-loss coefficient, γ 0 . Moreover, although Eq. (19) yields the exact solution at the single value of the propagation constant, given by Eq. (20), which is embedded in a generic family of numerically found fundamental solitons, as demonstrated below in Figs 1, 2 and 3, the entire family asymptotically shrinks to the exact solution in the limit of large γ 0 . Indeed, it is easy to find that, for γ  1 0 2 and a relatively small deviation of the propagation constant from the special value (20) www.nature.com/scientificreports www.nature.com/scientificreports/ Next, Eq. (11) with the imaginary potential taken as per Eq. (4) with β = 0, and with σ ≠ 0 in the nonlinearity-modulation profile (9), gives rise to the following exact solution, at the respective single value of k:   (3) and (9). Stable fundamental single-peak solitons are marked by green dots. All unstable solitons are marked by red crosses, irrespective of their structure. Exact soliton solutions, given by Eqs (19) and (20), are indicated by green stars (except for one at γ = 2 0 , which is designated by the red cross, as the exact solutions are unstable at γ ≥ 2 0 ). Green numbers ≥2 in this figure and below denote stable solitons with the same number of peaks. Further, green numbers 1 label stable single-peak solitons with the higher-order radial structure, as in Fig. 1(b). Green numbers 1 or 2, placed close to green dots, imply bistability, i.e., coexistence of stable fundamental single-peak solitons and stable higher-order or double-peak ones. Red crosses placed on top of green dots imply coexistence of fundamental single-peak solitons with some unstable mode. Soliton solutions were not found in white areas.
www.nature.com/scientificreports www.nature.com/scientificreports/ In this case too, the asymptotic approximation Eq. (14) becomes identical to the exact solution at = k k xy 0 ( ) , both existing at γ < 2 0 , on the contrary to exact solution (19), which exists at all values of γ 0 . Thus, the models considered here do not have the linear spectrum. Instead of it, they are characterized by spectra (families) of self-trapped nonlinear solutions (solitons). The radical change of the concept of the system's spectrum implies a respective change in the concept of the PT symmetry, which now applies not to the set of eigenvalues of the linearized system, but directly to families of nonlinear states. Lastly, it is worthy to note that all the asymptotic and exact solutions produced above, including the first correction (16) to the asymptotic tails, feature isotropic shapes of |U(x, y)|, although the imaginary potentials Eqs (3) and (4) are obviously anisotropic.
Exact solutions for elliptic vortices in an anisotropic model. In addition to 2D fundamental solitons, similar to the exact ones presented here, we also address below, by means of numerical methods, solitons with embedded vorticities, = … m 1, 2, 3 . A challenging issue is to seek for exact solutions for vortex solitons. Such solutions can be found in the case of imaginary potential Eq. (4) with β = 0, in a more general anisotropic version of the nonlinearity-modulation profile in Eq. (8) with σ = 0, namely, where positive ≠ g 1 accounts for the ellipticity of the modulation profile. Then, an exact solution for elliptically deformed vortex solitons with = m 1 is given by the following ansatz [cf. Eq. (22)]: where real ≠ b 1 accounts for the ellipticity of the soliton's phase field, and a is another real constant. The substitution of this ansatz and expressions Eqs (4) and (24) (with β = 0) in the accordingly modified Eq. (11) leads to the following relations between parameters of the ansatz: ,  www.nature.com/scientificreports www.nature.com/scientificreports/ The system of three equations (26) for two free parameters a and b demonstrates that the exact vortex solution is a nongeneric one, as it may exist only if an additional constraint, which can be derived by eliminating a and b in Eq. (27), is imposed on parameters g and γ 0 : In the isotropic model, with = g 1, Eq. (26) has no nontrivial solutions. However, they can be found for ≠ g 1. A particular example is Numerical results for zero-vorticity solitons. The PT -symmetric imaginary potential (3): single-and double-peak solitons. The isolated exact solution of the model with the PT -symmetric gain-loss distribution, given by Eqs (19) and (20), can be embedded in a continuous family of solitons, produced by a numerical solution of Eq. (11), with Σ r ( ) and γ(x) taken as per Eqs (3) and (9). The appropriate numerical algorithm is the Newton conjugate gradient method 77 , which is briefly outlined in section Method below. The stability of the stationary states was identified by numerical computation of eigenvalues of small perturbations, using linearized Eq. (32) for perturbations around the stationary solitons. Finally, the stability predictions, produced by the eigenvalues, were verified by simulations of the perturbed evolution of the solitons (some technical details are reported elsewhere 63 ).
It is relevant to stress that the convergence of the algorithm which produces stationary states depends on appropriate choice of the initial guess. While stationary modes were not found in " holes" appearing in stability charts which are displayed below in Figs 2, 3, 7, 10, 11, 12, 16 and 17, it is plausible that stationary solutions exist in the holes too, being, however, especially sensitive to the choice of the input. On the other hand, the intricate alternation of stability and instability spots, which is also observed in the charts, is a true peculiarity of the present www.nature.com/scientificreports www.nature.com/scientificreports/ model. Moreover, genuine structure of the stability charts may be fractal, but analysis of this possibility is beyond the scope of the present work.
Generic examples of numerically found stable solitons with single-and double-peak shapes are displayed in Fig. 1. Note that the double-peak modes have their two maxima separated in the direction of x, in accordance with the anisotropic shape of the imaginary potential in Eq. (3). As concerns single-peak modes, two different varieties of stable ones were found: fundamental solitons, with the shape similar to that of the exact solution given by Eqs (19) and (20) [see Fig. 1(a)], and higher-order states with a radial ring surrounding the central peak, see  www.nature.com/scientificreports www.nature.com/scientificreports/ Fig. 1(b). It is worthy to note that, unlike many other models, where higher-order radial states are unstable [78][79][80][81][82][83] , they are stable in the present case. Note also that shapes of both species of the single-peak solitons, fundamental and higher-order ones, seem isotropic in terms of |U(x, y)|, similar to exact solution (19). The isotropy is obviously broken by double-peak modes, see Fig. 1(c).
Results of the stability analysis, based on the computation of perturbation eigenvalues, are summarized in the stability map in the plane of (k, γ 0 ) [the soliton's propagation constant and strength of the gain-loss term in Eq.
0 2 in Figs 2 and 3, respectively. Several noteworthy features are revealed by these plots. First, it is worthy to note significant stability areas for both the double-peak and higher-order single-peak PT  9), in panels (a and b) panels, respectively. Green circles and red crosses denote stable and unstable vortex solitons, respectively. The same notation is used below in other stability charts for vortex solitons. www.nature.com/scientificreports www.nature.com/scientificreports/ -symmetric solitons in Figs 2 and 3. Further, bistability is observed at many points, in the form of coexisting stable fundamental and double-peak solitons, or fundamental and higher-order radial states of single-peak ones. As concerns the possibility of maintaining the unbreakable PT symmetry, Fig. 2 demonstrates shrinkage of the existence and stability regions of the modes with the increase of γ 0 at β = 0 to the exact soliton solution given by Eqs (19) and (20), in agreement with the trend represented by approximate solution (21). Eventually, the exact solution loses its stability at γ ≥ 2 0 . On the other hand, the introduction of a relatively weak confinement of the gain-loss term, with β = .
0 2 in Eq. (3), demonstrates that the PT symmetry remains unbreakable in Fig. 3, where  . The same as in Fig. 6, but for a case when the stable vortex soliton with = m 1, featuring a complex shape, is created, the parameters in Eqs (3) and (9)  www.nature.com/scientificreports www.nature.com/scientificreports/ both the existence and stability regions extend in the direction of large values of −k and γ 0 , without featuring any boundary.
As concerns unstable solitons, they typically blow up in the course of the evolution, see an example below in Fig. 18. Although it shows the blowup of a vortex soliton, the instability development of zero-vorticity ones is quite similar.
The stability charts, drawn in Figs 2 and 3 for σ = 1 in Eq. (9), are similar to their counterparts produced at other values of σ, including σ = 0, when the exact solution given by Eqs (19) and (20) does not exist, while the asymptotic form of the solitons' tails is given by Eq. (13).
The imaginary potential (4): single-and multi-peak solitons. A drastic difference revealed by the stability analysis of the model based on Eqs (4), (8) and (9) is that the respective exact solutions, given by Eq. (18) for the special case (17), and by Eqs (22) and (23)     www.nature.com/scientificreports www.nature.com/scientificreports/ As mentioned above, the steep growth of Σ r ( ) in Eq. (9) cannot extend to infinity, it being sufficient to maintain the adopted profile of Σ r ( ) on a scale which is essentially larger than a characteristic size of solitons supported by this profile. The same pertains to the linear growth of the imaginary potential at large |x| in Eq. (3): in reality, it should not continue at distances much larger than the size of the stable solitons considered in the previous section. However, the presence of β min implies that the corresponding "tacit" confinement of γ(x, y) in Eq. (4) is not sufficient to produce stable 2D solitons. At β β > min , the numerical solution generates stable fundamental single-peak solitons and their higher-order radial counterparts with isotropic shapes of |U(x, y)|, as shown in Fig. 4(a,b). Further, stable multi-peak solitons are found too. Due to the 2D structure of the imaginary potential (4), they feature a four-or five-peak structure, built along both the x and y axes, as shown in Fig. 4(c,d), instead of the uniaxial double-peak modes supported by the quasi-1D imaginary potential (3), cf. Fig. 1(c).
A typical stability chart for the 2D solitons generated by the model with β β > min is displayed in Fig. 5. It features bistability between the fundamental single-peak solitons and the higher-order ones, or four-and five-peak complexes, in a relatively small region of the (k, γ 0 ) plane, at sufficiently small values of γ 0 . Figure 5 clearly shows that no solitons were found at γ ≥ 2 0 , this restriction coinciding with that for the exact solution given by Eqs (22) and (23). Thus, unlike the PT -symmetric imaginary potential (3), the model based on potential (4) does not produce unbreakable soliton families.
Vortex solitons. Soliton solutions of Eq. (11) with embedded vorticity were found numerically by means of the above-mentioned Newton conjugate gradient method, initialized by the ansatz with integer vorticity ≥ m 1 added to the previously found 2D stationary solutions of Eq. (11): Vortex solitons in the case of the PT -symmetric imaginary potential. In the framework of the model with imaginary potential (3), stable vortex solitons were found in the case of β = 0 (no gain-loss confinement) with = m 1, while vortices with ≥ m 2 do not exist or are unstable. An example of stable vortices is shown in Fig. 6, and the respective stability charts for different values of σ in Eq. (9) are presented in Fig. 7. The strongly anisotropic shape of the vortex is a consequence of the anisotropy of the underlying imaginary potential (3). www.nature.com/scientificreports www.nature.com/scientificreports/ The introduction of the confinement of the gain and loss in Eq. (3) (in particular, setting β = . 0 5) makes it possible to construct stable vortex solitons with higher vorticities, corresponding to > m 1 in Eq. (28). An example of a stable vortex with = m 3 is shown in Fig. 8. In most cases, stable vortices generated by input (28) from double-peak stationary solutions have the same shape as those originating from their single-peak counterparts. However, in few cases the application of the lowest vorticity, with = m 1 in Eq. (28), to the double-peak input leads to the creation of stable vortex solitons with a complex shape, see an example in Fig. 9.
Stability charts for the vortex solitons with m = 1, 2, and 3, supported by the PT -symmetric imaginary potential which is subject to the spatial confinement, with β = .
0 5 in Eq. (3), are shown in Figs 10, 11 and 12. While the stability area shrinks with the increase of m, a few stable isolated modes were found even for = m 4 (not shown here). The comparison of Figs 7 and 10 shows that the introduction of the spatial confinement of the gain-loss profile helps to expand the stability area for = m 1 towards larger values of γ 0 , thus upholding the trend to observe the unbreakable PT symmetry in this 2D model. In direct simulations, the evolution of unstable vortex modes leads towards the blowup, via their fusion into a single peak, similar to what is displayed below in Fig. 18.
Vortex solitons in the model with imaginary potential (4). Starting from input Eq. (28), stable vortices can be constructed in the model with the gain-loss profile given by Eq. (4) only if it is subject to the spatial confinement (recall the same is reported above for zero-vorticity solitons). Examples of stable solitons with vorticities m = 1, 2 and 3 found in this model are shown in Figs 13,14 and 15. Note that higher-order states with ≥ m 2 are actually compound modes built of m unitary vortices, whose pivots do not merge into a single one, remaining separated, although with a small distance between them, as can be seen for = m 2 in Fig. 14 (cf. a similar effect recently reported in ref. 84 ). The separated pivots form arrays along axes x or y, the particular direction being randomly chosen by the initial conditions, as is clearly seen in Fig. 15. Nevertheless, the overall shapes of the unitary and higher-order vortices are nearly isotropic, due to the structure of the gain-loss term in Eq. (4) (cf. strongly anisotropic shapes of vortices in Figs 6, 8 and 9, supported by the imaginary potential (3)).
Stability charts obtained in this model for the solitons with embedded vorticities m = 1 and 2 are shown in Figs 16 and 17. Only few examples of stable vortices with m = 3 have been found in this case (for instance, the one shown in Fig. 15, as well as at σ = 0, γ = . 0 4 0 , = − . k 1 2). Finally, a generic example of the evolution of an unstable vortex soliton is shown in Fig. 18. The strong difference between vertical scales in different panels of the figure clearly suggests that the instability leads to the blowup of the unstable mode, in the course of which the original vortex tends to fuse into a single peak. In fact, all unstable solitons considered in this work tend to develop the blowup in direct simulations.

Discussion
The objective of this work is to elaborate 2D models with the spatially modulated self-defocusing nonlinearity and gain-loss distributions [imaginary potentials, iW(x, y)] which give rise to families of stable single-peak, multi-peak, and vortical solitons, including ones which may persist and remain stable ("unbreakable") at arbitrarily large values of strengths γ 0 of the imaginary potential. The unbreakability is possible in the case of the PT -symmetric imaginary potential, which is given by Eq. (3). An asset of the models, which can be implemented in bulk nonlinear optical waveguides with embedded gain and loss elements, is that they produce universal asymptotic solutions for solitons' tails, along with full exact solutions for selected species of 2D fundamental and vortex solitons (the latter one is available in the elliptically deformed version of the model). In particular, in the limit of large γ 0 , the unbreakable family of fundamental solitons tends to shrink towards the exact solution. Generic www.nature.com/scientificreports www.nature.com/scientificreports/ families of zero-vorticity solitons, including single-and multi-peak ones and higher-order radial states of single-peak solitons, as well as families of self-trapped modes with embedded vorticity m = 1, 2, and 3, are constructed in the numerical form, and their stability is identified by means of the numerical computation of eigenvalues for small perturbations, and verified by direct simulations. In the case of the PT -symmetric imaginary potential (3) the solitons are stable in vast parameter regions, and feature a trend towards maintaining the unbreakable PT symmetry. Under the action of the imaginary potential (4), families of stable fundamental and vortex solitons exist too, provided that the imaginary potential is subject to spatial confinement.
A relevant extension of the analysis may be to address the elliptically deformed model, which is considered in the present work in a brief form. A challenging problem is the possibility of the fractal structure of the stability patterns in the models' parameter planes.

Methods
The Newton conjugate gradient method for the 2D robust PT-symmetry model. Solutions of the stationary Eq. (11) were constructed by means of the Newton conjugate gradient method, which is presented in detail in book 77 . In terms of this method, the stationary-solution operator L 0 is defined by Eq. (11), while the respective linearization operator L 1 is defined as where the nonlinearity coefficient, Σ r ( ), and imaginary potential, W(x, y) are defined, respectively, by Eqs (3) or (4) and (9).
Simulations of the evolution of the wave fields. Direct    The solutions were numerically constructed in the 2D spatial domain, |x, y| ≤ 9, which was covered by a discrete grid of size × = × N N 512 512 x y . The direct simulations were carried out with step ∆ = − z 10 5 . This small step was selected to provide sufficient accuracy of the numerical solutions obtained in the presence of the " exotic" nonlinearity-modulation and gain-loss profiles (9)