## 1 Introduction  
The study of \(\mathcal{PT}\) (parity-time) symmetric systems was initiated through the works of [PERSON] and collaborators [1, 2]. Originally, it was proposed as an alternative to the standard quantum theory, where the Hamiltonian is postulated to be Hermitian. In these works, it was instead found that Hamiltonians invariant under \(\mathcal{PT}\)-symmetry, which are not necessarily Hermitian, may still give rise to completely real spectra. Thus, the proposal of [PERSON] and co-authors was that these Hamiltonians are appropriate for the description of physical settings. In the important case of Schrodinger-type Hamiltonians, which include the usual kinetic-energy operator and the potential term, \(V(x)\), the \(\mathcal{PT}\)-invariance is consonant with complex potentials, subject to the constraint that \(V^{*}(x)=V(-x)\). A decade later, it was realized (and since then it has led to a decade of particularly fruitful research efforts) that this idea can find fertile ground for its experimental realization although not in quantum mechanics where it was originally conceived. In this vein, numerous experimental realizations sprang up in the areas of linear and nonlinear optics [3, 4, 5, 6, 7, 8], electronic circuits [9, 10, 11], and mechanical systems [12], among others. Very recently, this now mature field of research has been summarized in two comprehensive reviews [7, 8]. One of the particularly relevant playgrounds for the exploration of the implications of \(\mathcal{PT}\)-symmetry is that of nonlinear optics, especially because it can controllably involve the interplay of \(\mathcal{PT}\)-symmetry and nonlinearity. The relevant efforts have gradually progressed from the simpler setting of coupled waveguides (bearing gain and loss) to entire \(\mathcal{PT}\)-symmetric lattices and the identification of optical solitons in them. These developments hold considerable promise for the potential realization of more complex settings, such as the one proposed herein. In this context,the propagation of light (in systems such as optical fibers or waveguides [7; 8]) is modeled by the nonlinear Schrodinger equation of the form:  
\[i\Psi_{z}+\Psi_{xx}+\Psi_{yy}+U(x,y)\Psi+\sigma|\Psi|^{2}\Psi=0. \tag{1}\]  
In the optics notation that we use here, the evolution direction is denoted by \(z\), the propagation distance. Here, we restrict our considerations to two spatial dimensions and assume that the potential \(U(x,y)\) is complex valued, representing gain and loss in the optical medium, depending on the sign of the imaginary part (negative for gain, positive for loss) of the potential. In this two-dimensional setting, the condition of full \(\mathcal{PT}\)-symmetry in two dimensions is that \(U^{*}(x,y)=U(-x,-y)\). Potentials with full \(\mathcal{PT}\) symmetry have been shown to support continuous families of soliton solutions [13; 14; 15; 16; 17; 18]. However, an important recent development was the fact that the condition of (full) \(\mathcal{PT}\) symmetry can be relaxed. That is, either the condition \(U^{*}(x,y)=U(-x,y)\) or \(U^{*}(x,y)=U(x,-y)\) of, so-called, partial \(\mathcal{PT}\) symmetry can be imposed, yet the system will still maintain all real spectra and continuous families of soliton solutions [19]. Other subsequent results include stable vortex structures in cubic nonlinear media under the influence of partially \(\mathcal{PT}\)-symmetric azimuthal potentials [20]. In the original contribution of [19], only the focusing nonlinearity case was considered for two select branches of solutions and the stability of these branches was presented for isolated parametric cases (of the frequency parameter of the solution). Our aim in the present work is to provide a considerably more \"spherical\" perspective of the problem. The corresponding physical scenario that we have in mind involves a medium featuring a cubic nonlinearity. Given that the potential of [19] was principally featuring four \"nodes\", we are envisioning an implementation involving a four waveguide system.