Three-dimensional topological solitons in PT-symmetric optical lattices

We address the properties of fully three-dimensional solitons in complex parity-time (PT)-symmetric periodic lattices with focusing Kerr nonlinearity, and uncover that such lattices can stabilize both, fundamental and vortex-carrying soliton states. The imaginary part of the lattice induces internal currents in the solitons that strongly affect their domains of existence and stability. The domain of stability for fundamental solitons can extend nearly up to the PT-symmetry breaking point, where the linear lattice spectrum becomes complex. Vortex solitons feature spatially asymmetric profiles in the PT-symmetric lattices, but they are found to still exist as stable states within narrow regions. Our results provide the first example of continuous families of stable three-dimensional propagating solitons supported by complex potentials.

Generation of three-dimensional (3D) solitons has been a salient problem of fundamental importance since the birth of nonlinear physics. The problem critically hinges in the elucidation of physical settings that allow existence of three-dimensional self-sustained excitations that, in addition, are stable upon propagation. The latter requirement is particularly challenging because the common cubic (Kerr) nonlinearity present in most potentially suitable materials leads to supercritical collapse and thus cannot support stable higher-dimensional solitons in uniform media [1,2]. A number of approaches to stabilize 3D solitons have been suggested over the years [3,4]. Most of them suggest using nonlinearities that are different from pure cubic nonlinearity, introduce various higher-order effects, or rely on spatial modulations of system parameters. Thus, it has been predicted that stable 3D solitons (termed light bullets in optics) can be supported by media with saturable [5] and quadratic [6][7][8] nonlinearities, materials with competing [9][10][11] or nonlocal [12][13][14][15] nonlinearities, offresonant two-level systems [16], etc. Stable evolution of 3D solitons is shown to be possible under the action of higher-order effects arising upon filamentation [17] or in the presence of higher-order dispersion or nonparaxial corrections [18]. 3D solitons were studied in dissipative settings too, where in addition to competing conservative nonlinearities, higherorder absorption is usually present [19][20][21][22][23][24][25].
A powerful strategy for realization of stable 3D bullets (as well as of topological 2D states [26,27]) relies on the spatial modulations of material parameters. Thus, longitudinal tandems [28,29], graded-index fibers [30], and the structures with combined transverse and longitudinal modulations [31], such as shaken optical lattices [32], can be used to suppress collapse. Stabilization of 3D solitons in a transversally periodic medium was predicted in discrete waveguide arrays [33][34][35]. This idea was trans-ferred to continuous optical systems [36][37][38][39] and also to Bose-Einstein condensates (BECs) [40][41][42]. Importantly, this strategy afforded the first experimental observation of fundamental light bullets in fiber arrays [43] that was later extended to vortex solitons [44]. In all previous works, the transverse linear potential stabilizing the 3D solitons was real, i.e. it was created only by the refractive index modulation (in optics) or by using standing optical lattices (in BECs). Nevertheless, the recent discovery [45] of a new class of so-called  -symmetric complex potentials, obeying the condition ( ) ( ) * = r r   and possessing purely real spectra when Im ( ) r  is below certain critical level, has opened important new horizons for the control of linear and nonlinear light propagation.
The concept of  -symmetry introduced initially in quantum mechanics upon construction of non-Hermitian operators with real spectra has been extended to several areas of science. Among them are lasers with single-mode operation and enhanced tunability [46,47], resonant systems such as atomic gases [48], design of metamaterials and metasurfaces where symmetry breaking may manifest itself in changes of the state of polarization of reflected light [49], formation of light-matter excitation in dissipative exciton-polariton condensates [50], micro-cavities [51], or acoustic devices [52], to name a few. Importantly, many of these systems may operate in the nonlinear regime where under appropriate conditions self-sustained nonlinear states may form. Elucidation of suitable conditions for the existence of such stable excitations, especially multidimensional ones, is a topic of continuously renewed interest. Optics provides the ideal setting to study the concept.
It should be mentioned that optical waveguiding systems offer particularly convenient platforms for the realization of  -potentials [53,54], because they allow the simultaneous modulation of the refractive index and the gain and losses.  -symmetry has been experimentally observed in several settings [55,56]. There is considerable interest in the investigation of solitons in  -symmetric potentials. The properties of 1D and 2D solitons in periodic  -symmetric lattices are known [57][58][59][60][61][62][63][64][65][66], but 3D solitons have not been achieved. Note that previous works [67,68] about 3D solitons used complex potentials localized in three dimensions and employed symmetries absent in a periodic medium. Thus, the question of whether stable 3D solitons are possible in complex potentials remains open. In this paper we show for the first time that stable fundamental and vortex 3D solitons can form in 2D  -symmetric periodic optical lattices imprinted in focusing cubic media. To the best of our knowledge, this is the first known example of stable 3D solitons supported by complex periodic lattices that do not require higher-order dissipation for stable propagation. We address the propagation of a spatiotemporal wavepacket along the x -axis of a cubic nonlinear optical medium with anomalous dispersion and imprinted transverse modulation of the refractive index and gain/losses. The evolution of the wavepacket is governed by the dimensionless nonlinear Schrödinger equation: is the dimensionless field amplitude; 2 n is the nonlinear coefficient; is the diffraction length; 2 / k n p l = is the wavenumber; n is the refractive index;   Here we assume that the lattice profile does not change in time. Note that temporal lattices have been created experimentally [69], but the extension of the technique to realize the model studied here is not straightforward. In contrast, since doping of fibers with active or absorbing centers is a well-established technology, properly designed arrays of dissipative fibers, similar to conservative arrays from [43,44], may be potentially readily used for the experimental observation of the phenomena described here. Indeed, spatial modulation of gain and losses may be realized by packing absorbing and amplifying fibers into one complex periodic array and using properly shaped pump beams. By assuming similar nonlinear and dispersive properties of the material as in the experiments that addressed the formation of conservative light bullets that were performed in waveguide arrays with silica cores [43,44], for a characteristic transverse scale of Experiments with pulsed excitation in  -symmetric systems may also be performed with other materials featuring much stronger nonlinearities, see e.g. [55,56].
On rigorous grounds, once dispersion of the dielectric permittivity of the medium is taken into account, the exact  -symmetry may be realized only for a discrete set of optical frequencies, as follows from the Kramers-Kronig relations [70]. Nevertheless, as long as gain and losses are provided by dopants in such low concentrations that do not notably affect the dispersion of the real part of the permittivity of the host material, as it is assumed here, and, in turn, as long as the imaginary part of the permittivity is determined exclusively by the dopants, the deviations from the exact  -symmetric conditions may not have a significant effect within the frequency range occupied by the propagating soliton pulses. Moreover, gain and losses may be provided by different mechanisms, including nonlinear ones [56], which may further reduce or even vanish the impact of the issue. As a specific setting for the experimental implementation of the model (1) we also suggest  -symmetric refractive index landscapes imprinted in a cold gas of two atomic isotopes in a Λ-type configuration (e.g., 87 Rb and 85 Rb isotopes) loaded in an atomic cell, as described in [48]. Due to the interference of two Raman resonances, the required spatial distribution of the refractive index may potentially be realized by using a proper combination of a control laser field and a far-off-resonance Stark laser field. Since the refractive index of the atomic vapor is determined by the two external laser fields whose intensities could be at microwatt levels, one may use them for fine-tuning the  -symmetric potential. In addition, such a system is characterized by a large Kerr nonlinearity (sometimes exceeding usual non-resonant nonlinearities by many orders of magnitude [71]) that may favor the formation of solitons at remarkably low power levels. At the same time, inducing the required dispersion in such atomic system may be a challenge.
The existence domains of 3D solitons are closely connected to the existence domains of Substitution of the field in such a form into the linear version of Eq. (1) without the temporal derivative yields the linear eigenvalue problem:  Second, in a nonlinear dispersive medium, exponential localization of the tails of solitons is guaranteed only when b remains above the edge of the continuous spectrum both in the temporal and spatial domains. Thus, the lower cutoff for the existence of the 3D solitons is given by . Figure 2 summarizes the main properties of the stationary 3D solitons supported by  -symmetric lattices, and Fig. 3 illustrates representa-tive profiles of such states. Stationary soliton solutions were obtained from Eq. (1) in the form , where the function a describes the soliton shape that remains unchanged upon propagation in the absence of perturbations. Here a is complex, since even for a fundamental solution the internal energy currents from amplifying to absorbing domains appear due to nonzero Im ( , ) h z  , but the propagation constant b is always real, since we are searching for nonlinear states for which gain and losses are integrally balanced. Solutions were obtained with the method of squared operator described in Ref. [72] using fast Fourier transforms similarly to spectral renormalization methods [73]. The method consists in the iterative calculation of the soliton shape using the expression †  Eq. (1) is intrinsically dissipative, therefore the total energy is not a rigorous conserved quantity anymore (there exist inputs for which U grows or decays upon propagation).
Instead, conservation of the quantity ( , , ) ( , , ) is guaranteed by the symmetry of ( , ) h z  . Nevertheless, it is convenient to use the ( ) U b dependencies to characterize the soliton families, since the energy U is a directly measurable experimental parameter and it does not change upon evolution when a stationary unperturbed state is considered. Note that the existence of entire soliton families parameterized by the propagation constant b despite the presence of gain and losses is a unique property of  -symmetric potentials.  p > is required for appearance of nonmonotonic ( ) U b curves for i 2 p = ]. As shown in the plot, we found that there is no energy threshold for the existence of fundamental 3D solitons in the  -symmetric lattice (the energy tends to zero when b ¥ ). Indeed, when solitons become strongly localized at large b values, the lattice affects their shapes only weakly and the ( ) U b curves for different r,i p values approach the monotonically decreasing dependence as in uniform medium [10,11]. Scaling arguments that can be used to obtain the latter dependence nevertheless become inapplicable when the lattice strongly affects soliton shapes, hence at intermediate b values one observes considerable deviations from the Importantly, note that Fig. 2(a) shows that stability can potentially be achieved only within a finite interval of energies and propagation constants.
When propagation constant approaches the cutoff value co b , the soliton energy increases rapidly. In the regime i r p p < with unbroken symmetry, this cutoff coincides with the upper edge of the first band max re b in the spectrum from Fig. 1. Hence, close-to-cutoff solitons strongly expand both in space and time (see representative shapes in Fig. 3, left column), featuring multiple spatial oscillations. Increasing b results in a higher localization of the solitons and in the gradual equilibration of its spatial and temporal widths (Fig. 3, central column). Figure 3 also   where r Re =   and i Im =  . The magnitude of the current is directly determined by the imaginary part of the potential. Moreover, the current j strongly impacts the soliton shape. It should be noted that the soliton phase varies not only in space, but also in time t . This is illustrated in the last row of Fig. 3. However, the temporal variation of the phase is rather slow in comparison with its spatial variation, so that the corresponding temporal current component ¶ ¶ is small and therefore affects only weakly the soliton shape in comparison with , j h z . Figure 2(b) illustrates the impact of the imaginary part of the lattice on the properties of the 3D fundamental soliton families. Increasing i p results in a decrease of the cutoff for soliton existence. The cutoff vanishes for i r p p ³ . When the  -symmetry is unbroken, both spatial w h and temporal w t integral soliton widths, defined according to Fig. 2(c)]. The ratio of spatial and temporal widths A qualitatively similar picture was obtained for 3D solitons carrying a nested off-site topological dislocation, or vortex light bullets. Such solitons are composed of four optical peaks located in the vicinity of four neighboring lattice maxima. Traces of a canonical, staircase vortex phase distribution and phase singularity are visible in Fig. 3 (right), but the entire phase distribution y is rather complex due to the presence of internal currents from amplifying to absorbing domains. Snapshots of the intensity distributions carried by the spatio-temporal vortex solitons at different t values reveal their spatial asymmetry. On physical grounds, such asymmetry is a consequence of the interference of the global current associated with vorticity and the local currents arising inside waveguides due to inhomogeneous gain and losses. Namely, in some waveguides the directions of global and local currents coincide, leading to an increase of total current, but in other waveguides the currents take opposite directions, which results in current reduction. As a consequence, the ( , ) h z j distributions at any t do not feature a four-fold rotation symmetry as in conservative lattices, and such spatial asymmetry directly translates into an asymmetry of the field distribution q as it follows from Eq. (3). The effect is readily visible in Fig. 3, where we compare the field distributions (modulus) of the 3D spatio-temporal vortex solitons with positive and negative topological charges l . The asymmetry in the shape of the vortex solitons increases when the imaginary part of the lattice grows (compare 0 t = cross-sections for vortex solutions shown in right column of Fig. 3, obtained for i 6 p = , and the solution shown in Fig. 4, obtained for i 0.55 p = ). Typical ( ) U b curves for vortex solitons are presented in Fig.  2(d). Curves are qualitatively similar to those featured by fundamental solitons, but the energy carried by a vortex solitons is approximately four times higher than that of a fundamental soliton. It should also be noted that the 3D vortex solitons do not bifurcate from the upper edge of the band, i.e. their cutoff is slightly higher than that of fundamental solitons.  Fig. 2(b). In all cases r 6 p = . Distributions in the secondfifth rows are shown in the , , [ 4,4] h z t Î -+ windows.
The central result of this paper is summarized in Fig. 5, which depicts the stability properties of both, fundamental and vortex 3D solitons in  -symmetric lattices. In order to check stability we first resort to direct integration of Eq. (1) with inputs in the form of solitons solutions perturbed with weak small-scale noise. Such perturbed states were allowed to propagate for large distances, up to 3 10 x = , which allow to identify the parameter domains where 3D solitons are dynamically stable. This method of stability analysis tests the presence of a wide range of potential instabilities, because the broadband input noise excites a correspondingly broad spectrum of perturbations. Due to the large propagation distances considered, exponentially growing modes are readily detectable. At the same time, because noise is introduced with small amplitude, the perturbation introduces only weak deformations into the input soliton that thus do not lead to pronounced longitudinal oscillations of its amplitude/width. We found that fundamental solitons are stable in a considerable part of their existence domain on the i ( , ) p b plane, nearly up to the  -symmetry breaking point [in the region between red dots in Fig. 5(a)]. The Vakhitov-Kolokolov criterion correctly predicts the lower border low cr b of the stability domain for the fundamental solitons that is very close to the cutoff for soliton existence, but it does not predict its upper border upp cr b for i 0 p > . This is because at upp cr b b = oscillatory rather than exponential instabilities come into play. These are accompanied by transverse shift/oscillations of the soliton center. The stability domain shrinks at first with the increase of i p , but then it expands again. It completely disappears slightly below the symmetry breaking point i r p p = , where the spatial background becomes unstable. We also found that the  -symmetric lattices can stabilize 3D vortex solitons, even though they are spatially asymmetric. The domain of stability for the vortex solitons [ Fig. 5(b)] is much narrower than that of fundamental solitons in terms of i p , but it is more extended in terms of b . The 3D vortex solitons are also be stable within the interval low of propagation constants. The domain of stability shrinks completely at i 0.6 p » for r 6 p = , but its width increases with growing r p . The development of instability for the 3D vortex solitons at is initiated by the concentration of energy in one of the constituent sub-peaks and its subsequent collapse. In contrast, fundamental 3D solitons with upp cr b b > oscillate and radiate energy away until their parameters reach a stability interval and a new stable state is formed.
where 2 2 a a * -= - , 2 2 a a * + = +  . Standard eigenvalue solvers can give solution of this 3D problem only for small number of transverse points (even though the corresponding matrix is sparse), a procedure that however does not guarantee sufficient accuracy. Instead, here we used the method put forward in [74] and obtained the eigenvalues and eigenfunctions of  by solving the equation accelerating convergence that is analogous to the operator used to find soliton solutions above. This equation is solved using a fourth-order Runge-Kutta method. At each step in x the eigenvalue d is computed from Y as 1 1 , is a standard inner product. It was shown in Ref. [74] that for sufficiently broad class of input conditions, at large x this method converges to the perturbation mode Y whose ei-genvalue d has the largest real part max r d , provided that the corresponding solution is dynamically unstable. A typical dependence max r ( ) b d calculated with the described approach is shown in Fig. 5(c) for the case of 3D vortex solutions. The plot confirms that linear stability analysis predicts stable vortex soliton solutions (i.e., with max r 0 d = ) in the same interval of propagation constants that was obtained by direct propagation of the perturbed solutions presented above. The typical normalized decay distance for unstable soliton solutions can be estimated as max r 1/ x d  , which corresponds to an actual length of max dif r / L d . Summarizing, it is readily apparent that inclusion of the dissipative part of the potential may drastically modify the properties of the system. As it becomes non-Hamiltonian, strong internal currents emerge that substantially modify soliton shapes that are especially pronounced in the case of vortex solitons, and instabilities come into play that may destroy even fundamental solitons in the regions where they are otherwise stable in conventional conservative systems.
Nevertheless, we have revealed that complex (i.e., with real and imaginary parts arising from gain and losses)  -symmetric twodimensional lattices are capable of supporting robust, stable fully threedimensional fundamental and vortex solitons even in Kerr nonlinear media. We showed that the strength of gain/losses acting in the system affects strongly the soliton parameters and the corresponding stability domains. The states discovered here in the  -symmetric systems are the first known examples of stable, three-dimensional, propagating soliton families supported by complex potentials. It is readily apparent that dissipative  -symmetric lattices offer additional tools for controlling the shapes, domains of existence, and stability domains of threedimensional excitations as compared to conservative lattices. In particular, the dissipative effects introduce features in the soliton shapes and stability properties that are not present in comparable conservative settings.