Semidiscrete vortex solitons

We demonstrate a possibility of creation of stable optical solitons combining one continuous and one discrete coordinates, with embedded vorticity, in an array of planar waveguides with intrinsic cubic-quintic nonlinearity. The same system may be realized in terms of the spatiotemporal light propagation in an array of tunnel-coupled optical fibers with the cubic-quintic nonlinearity. In contrast with zero-vorticity states, semidiscrete vortex solitons cannot exist without the quintic term in the nonlinearity. Two types of the solitons, \emph{viz.}, intersite- and onsite-centered ones (IC and OC, respectively), with even and odd numbers $N$ of actually excited sites in the discrete direction, are identified. We consider the modes carrying the embedded vorticity $S=1$ and $2$. In accordance with their symmetry, the vortex solitons of the OC type exhibit an intrinsic core, while the IC solitons with small $N$ may have a coreless structure. Facilitating their creation in the experiment, the modes reported in the present work may be much more compact states than their counterparts considered in other systems, and they feature strong anisotropy. They can be set in motion in the discrete direction, provided that the coupling constant exceeds a certain minimum value. Collisions between moving vortex solitons are considered too.


I. INTRODUCTION
The creation of stable solitons in multidimensional geometry is a subject of intensive ongoing research in nonlinear optics, Bose-Einstein condensates (BECs), and other fields [1,2]. A fundamental problem which impedes straightforward making of such solitons is that the ubiquitous cubic self-focusing nonlinearity gives rise to the critical and supercritical collapse in the two-and three-dimensional (2D and 3D) space [3,4], which destabilizes the solitons, on the contrary to stable ones supported by the cubic nonlinearity in 1D.
Thus, stabilization of multidimensional solitons is the central problem in this field. One possibility is the use of spatially periodic potentials, induced by photonic crystals in optics [5][6][7] and by optical lattices in BEC [8,9]. In the limit of very deep periodic potentials, the nonlinear medium creates effectively discrete solitons [10][11][12][13][14][15]. Recently, another stabilization method was elaborated for two-component solitons, based on the use of the linear spin-orbit coupling in binary BEC [16][17][18], or its counterpart in the bimodal light propagation in planar waveguides [19].
The above-mentioned means for the stabilization of 2D and 3D solitons rely on the use of linear effects. The stability may also be provided by replacement of cubic nonlinearity by other terms, which do not lead to the collapse. In optics, stable 2D spatial solitons and spatiotemporal light bullets were created in quadratic nonlinear crystals [20,21], saturable [22] and nonlocal [23] nonlinear media, and, finally, in a bulk material which features competing cubic and quintic (CQ) focusing and defocusing interactions [24]. It is relevant to mention that the CQ nonlinearity, with negligible corrections from higher-order terms, was experimentally identified in optical materials such as CS 2 (for power densities up to hundreds of GW/cm 2 [25]) and colloidal suspensions of metallic nanoparticles [26,27]. In BEC, stable 2D matter-wave solitons were predicted in dipolar BECs [28,29] and microwave-coupled binary condensates [30], where the nonlinearity is cubic but effectively nonlocal. Recently, experiments have revealed soliton-like multidimensional matter-wave states, in the form of quantum droplets (QDs), filled by an incompressible ultradilute quantum fluid [31], in dipolar BECs [32][33][34], and in binary condensates with contact interactions [35][36][37][38][39]. The latter experiments followed the prediction of the stabilization of the QDs by the Lee-Huang-Yang (LHY) correction to the mean-field interactions [40], reported in Refs. [41,42]. The correction originates from quantum fluctuations around the mean-field states.
A still more challenging issue is the creation of stable 2D and 3D bright solitons with embedded vorticity, because such a state is subject to the instability against splitting into fragments by azimuthal perturbations, which is a more destructive (faster) factor than the collapse [2]. Similar to their fundamental (zero-vorticity) counterparts, bright vortex solitons can be stabilized by means of linear effects, such as lattice potentials [8,9,43], and with the help of modified nonlinearity. In particular, effectively two-dimensional stable discrete and lattice optical vortex solitons were predicted [10] and created in the experiment [11][12][13]. Stable 2D solitons with intrinsic vorticity (topological charge) S = 1 in the uniform CQ medium were first predicted in Ref. [44], which was then extended to S ≥ 2 [45,46]. The existence of stable 3D solitons with S = 1 in the same model was reported in Ref. [47]. In binary BECs, stable 2D vortex solitons with topological charges up to S = 5 have been predicted in the above-mentioned microwave-coupled condensates, as well as in QDs with embedded vorticity [48,49]. Stable 3D QDs filled by the "swirling" condensate with S = 1 and 2 were also revealed by the analysis [50].
A relevant extension of the analysis of topological solitons aims to predict stable vortex modes with an anisotropic shape. In optics, stable 2D discrete vortex states were predicted in 2D anisotropic discrete lattices [51,53]. In BEC, vortex solitons were predicted in anisotropic lattice potentials [52], as well as in one component of the spin-orbitcoupled binary BEC with the anisotropic dipole-dipole interaction between atoms [54][55][56]. Very recently, creation of a novel form of 2D anisotropic vortices, viz., semidiscrete QDs with imprinted vorticity, was elaborated in a system with one continuous and one discrete coordinates, realized as a binary condensate loaded in an array of parallel tunnelcoupled quasi-1D traps [57]. In that model, the intrinsic nonlinearity of the traps is represented by the combination of self-attractive quadratic LHY and repulsive cubic mean-field terms. A nontrivial peculiarity of the setting is identification of the vorticity, which is defined in terms of the global phase pattern carried by the semidiscrete state.
FIG. 1: The setting based on the array of planar optical waveguides (blue slabs), separated by gray isolating layers, with the continuous transverse coordinate, x, and the discrete one, n. Light is coupled into the array along the z direction.
Semidiscrete vortices were not yet considered in nonlinear optics, which is the subject of the present work. Because arrays of coupled nonlinear waveguides are ubiquitous systems in optics [14,15], the creation of nontrivial self-trapped states in the arrays is a very relevant objective. The present work aims to predict stable semidiscrete optical vortex solitons in a effectively 2D setting built as an array of parallel tunnel-coupled planar waveguides with the intrinsic CQ nonlinearity, as sketched in Fig. 1. This setting may be naturally considered as a semidiscrete one. The respective model is somewhat similar to the above-mentioned one developed for the array of BEC traps with the intrinsic quadratic-cubic nonlinearity, which also supports stable vortex solitons [58,59]. However, results concerning the shape and stability of semidiscrete vortex solitons in the present system are essentially different. In particular, they may feature strong anisotropy, unlike quasi-isotropic modes reported in Ref. [57], and they may be built with much tighter shapes.
The rest of the paper is structured as follows. The model is introduced in Sec. II, including two different realizations in optics, viz., in terms of the light propagation in the spatial domain, as outlined in Fig. 1, and in a spatiotemporal one, realized as an array of tunnel-coupled fiber-like waveguides. Systematic results for the existence and stability of vortex solitons with S = 1 and S = 2 are summarized in Sec. III. The paper is concluded by Sec. IV.

II. THE MODEL
The propagation of light in the waveguiding array displayed in Fig. 1 is modeled, in the paraxial approximation, by coupled nonlinear Schrödinger equations (NLSEs) [60,61] with the CQ terms [62]: where A n is the envelope of the electromagnetic field in the n-th guiding core, with carrier wavenumber k 0 = 2n 0 π/λ, κ is the strength of the tunnel coupling between adjacent cores, n 0 , n 2 and n 4 represent, respectively, the linear, thirdand fifth-order refractive indices of the medium, and ζ is loss coefficient. Equation (1) does not include multiphoton absorption, which would be represented by nonlinear losses, as we do not consider media with resonant interactions or ionization of the optical material, which would induce such losses [63]. The coupling constant is determined by the difference of the refractive index between the guiding cores and dielectric material which separates them, the carrier wavelength, and the thickness of the separating layer, exponentially decaying with the increase of the latter. In practical terms, experimentally relevant values of the thickness are tantamount to a few wavelengths (which corresponds to several microns). This choice leads to the propagation length, for which the coupling effects become essential, measured in several millimeters [14]. As concerns the losses, in available optical materials they take values dB/m or still smaller, which makes the losses negligible for experimentally available propagation lengths 10 cm. For this reason, the dissipative term is neglected in the following consideration.
The semidiscrete systems adequately modeled by Eq. (1) with the cubic-only nonlinearity (n 4 = 0), were first introduced, and soliton-like states in them were investigated, in Refs. [60,[64][65][66][67][68], in terms of the spatiotemporal propagation of light in arrays of optical fibers. On the other hand, the 2D discrete model with the CQ onsite nonlinearity, and discrete soliton modes in them, were considered in Ref. [69]. A 2D quasi-discrete model, combining a deep checkerboard potential and the CQ nonlinearity, was also addressed, along with its vortex-soliton modes, in Ref. [70].
By applying rescaling, Eq. (1) is cast in the normalized form, with the effective intersite coupling constant, C. The total power of the field in the scaled form, which is a dynamical invariant of Eq. (3), is defined as The model also conserves its Hamiltonian, To estimate physical parameters of the model, we assume that, for instance, the cubic-quintic material is CS 2 , which was used for the creation of stable fundamental 2D solitons in Ref. [24]. At the carrier wavelength λ = 800 nm, its parameters are n 0 = 1.61, n 2 = 3.1 × 10 −19 m 2 /W, and n 4 = 5.2 × 10 −35 m 4 /W 2 [71,72]. Then, relations between scaled units in Eq. (3) and physical ones are estimated as follows: x = 1 ⇐⇒ 0.235 µm, z = 1 ⇐⇒ 0.07 mm, and P = 1 ⇐⇒ 70 kW, if the thickness of the single planar waveguide is 0.5 µm. As mentioned above, the losses in the medium may be neglected for experimentally relevant propagation distances [24]. With respect to these estimates, the characteristic transverse size of self-trapped modes considered below, ∆x ∼ 20 − 40, corresponds to the physical width ∼ 5 − 10 µm, and the characteristic propagation distance, which is z ∼ 500 − 3000, amounts to 3.5 ∼ 20 cm. These values are realistic for feasible experiments, see Refs. [11][12][13][14]. Note also that ∆x is sufficiently large in comparison to λ, which justifies the use of the paraxial propagation equation (3).
The same model applies to the temporal-domain light propagation, with carrier group velocity V gr , in an array of tunnel-coupled optical fibers with the intrinsic CQ nonlinearity, if coordinate x in Eq. (3) is replaced by the temporal coordinate, τ = t − z/V gr . In that case, the solitons represent semidiscrete "light bullets", and Eq. (4) defines the total energy of the spatiotemporal optical signal.
Stationary semidiscrete vortex solitons with topological charge S, of the onsite-centered (OC) and intersite-centered (IC) types, with the pivot located, respectively, at a lattice site or between two sites, were produced by means of the power-conserved-square numerical method [73], initiated by inputs where A and α are positive real constants. For the OC and IC solitons, we set in Eq. (7) respectively. Examples of the inputs used for the generation of the OC and IC solitons are displayed in Fig. 2. To make the pictures clearer, in this and other figures values of local intensities, |U n (x)| 2 , are shown, at discrete values of n, by means of finite-width stripes separated by empty ones ["data" and "empty" stripes in Fig. 2(a)].
Then, stability of the stationary solutions was verified by direct simulations of Eq. (3) for propagation distance z = 1000. The soliton is stable if its intensity profile stays unchanged throughout the simulation. According to the above-mentioned estimate of parameters for CS 2 , z = 1000 corresponds to Z = 7 cm, which is sufficiently long for observing a stable spatial soliton in the experiment. In the numerical computations, control parameters are total power P (4) and coupling constant C in Eq. (3). In terms of vortex modes displayed below in Fig. 3, z = 1000 corresponds to 10 characteristic Rayleigh (diffraction) lengths with respect to the continuous and discrete coordinates, hence this distance is sufficient to make conclusions concerning the stability of the modes.
Typical examples of the vortex solitons of the OC and IC types are displayed in Fig. 3 [the sign of the vorticity, S = 1, is determined by the comparison of phase patterns in the figure with the standard expression, u ∼ exp (iSΘ), in the continuum coordinate plane with angular coordinate Θ]. The solitons are characterized by the number of effectively excited sites, N (individual waveguides in which the field takes non-negligible values). Another important characteristic is number N core of sites in the core of the semidiscrete vortex, i.e., waveguides in which stationary fields u n (x) [see Eq. (6)] cross zero (vanish) at x = 0, thus having opposite signs at x > 0 and x < 0, while zero crossing is absent in waveguides which do not belong to the core. Although fundamental semidiscrete solitons, with zero vorticity, are also characterized by finite N [60, [64][65][66][67][68], only vortices may feature the core. Examples displayed in Fig. 3 exhibit values of the excited and in-core sites 2 ≤ N ≤ 9 and 0 ≤ N core ≤ 3, with odd and even numbers N and N core pertaining, severally, to the vortex modes of the OC and IC types. Note that the states of the former type always have N core ≥ 1, and, due to their intrinsic symmetry, any vortex soliton of the OC type with odd S has a real odd modal function u 0 (x) in the central waveguide (n = 0), see Fig. 5(a) below. On the other hand, the IC vortex states may exist both with N core = 0 and N core ≥ 2, see panels (a,b) and (c,d) in Fig. 3. In the case of N core = 0, the vorticity is accounted for by opposite signs of fields u 0 (x) and u 1 (x) in Eq. (6), with the "virtual pivot" of the vortex set between n = 0 and n = 1. It is also worthy to note that the semidiscrete vortex solitons may be strongly anisotropic, elongated in the continuous direction, in the case of small C, see, e.g., Figs. 3(a,e).
Stability areas for the OC and IC species of semidiscrete vortex solitons, which are characterized by numbers N , are displayed in the (P, C) plane in Figs. 4(a,b). A remarkable feature of the stability chart is bistability, allowing coexistence of the vortices with equal powers P and different numbers of sites, and even tristability -in particular, overlap of the stability regions for the vortices of the IC type with N = 4, 6, 8, at P > 83 [see point X in Fig. 4(a)]. In fact, the multistability area is even broader, as Fig. 4(b) demonstrates that, in addition to these three IC modes, two stable vortices of the OC type exist too, with N = 5 and 7, at the same values of C and P . In the limit of P → ∞, which corresponds to the 2D continuum space, the multistability agrees with the known fact that 2D vortex solitons of the nonlinear Schrödinger equation with CQ nonlinearities become stable as the solitons expand to accommodate indefinitely growing values of the norm [45]. The calculation of values of the Hamiltonian, according to Eq. (5), demonstrates that the minimum of the Hamiltonian, i.e., the ground state of the system, is realized, in the case of bi-or tristability, by the mode with the smallest number of excited discrete sites. An example of this calculation is shown in Fig. 4(c,d). Typically, the tristability area for the IC solitons with N = 4, 6 and 8 is shaded in Fig. 4(c). Outside of the stability area but close to its border, unstable vortex solitons can also be found; in direct simulations, they split in fragments. Typical examples of the evolution of stable and unstable vortex solitons are displayed in Fig.  5. Far from the stability boundary of the vortices, only fundamental (zero-vorticity) solitons are produced by the numerical calculations. Note that the vortex solitons of the IC and OC types with only two and three excited sites (and, respectively, N core = 0 and N core = 1) realize two types of the solitons with the minimum size in the discrete direction. Such tightly localized semidiscrete modes were not reported in the previously studied semidiscrete system [57]. In fully discrete 2D models, the possibility that the smallest vortex soliton includes four sites was theoretically predicted [10] and experimentally demonstrated [12,13]. Figure 6(a) displays the propagation constant, β, for two types of compact solitons [IC and OC with (N, N core ) = (4, 0) and (N, N core ) = (5, 1), respectively] as a function of the total power, P , for fixed C. These results indicate that the β(P ) curves satisfy the Vakhitov-Kolokolov criterion, dβ/dP > 0, which is a well-known necessary stability condition [3,74]. Figure 6(b) shows the β(C) dependence for stable vortex solitons with different numbers of excited sites N and fixed P , showing that β gradually decreases with C.
To characterize anisotropy of the vortex soliton, we define parameter This definition is adopted from Ref. [57], where ε = 1 implied that solitons were effectively isotropic modes, while ε > 1 and ε < 1 indicated that their shape was anisotropic, namely, elongated in the continuous or discrete direction, respectively. getting less compact. Moreover, passage of ε(P ) and ε(C) through the level of ε = 1 implies that the vortex solitons may be tuned to the isotropic shape by selecting specific values of P or C. A natural question is how the competing self-focusing cubic and defocusing quintic onsite terms in Eq. (3) affect the existence and stability of semidiscrete vortex solitons. To address this issue, we here rescale the equation, to make the intersite coupling constant equal to 1 and admit the presence of a free coefficient in front of the quintic term as a free parameter:Ũ hence the rescaled total norm isP = √ CP . The substitution of rescaling (12) in Eq. (3) casts it in the following form, where we omit the tilde, and use symbol σ ≡ C, to stress that C appears here as a newly defined parameter controlling the relative strength of the quintic term: To characterize effects of the competition between the cubic and quintic nonlinearities, Fig. 7 displays anisotropy parameter ε [see Eq. (10)] as a function of σ for the IC and OC vortex solitons, as obtained from the numerical solution of Eq. (13) with a fixed total power, P = 250. The figure shows that, at σ > 0.2, ε stays very close to ε = 1, which indicates that the semidiscrete vortex solitons keep an effectively isotropic profile (hence, they have a broad quasi-continuum shape) in this case. At σ < 0.2, the vortex solitons feature a compact shape, adequately characterized by the number of excited sites, N . The segments, which are characterized by the number of excited sites 5 ≤ N ≤ 9, terminate in Figs. 7(a,b) at their top and bottom points because solutions could not be found above or below them. At still smaller values of σ the vortex solitons are unstable, and no vortex modes were found at σ → 0, when the CQ nonlinearity turn into the cubic form. Thus, the inclusion of the self-defocusing quintic term is necessary for the existence and stability of the semidiscrete vortex solitons. A challenging issue is whether stable vortex solitons can be found for the double topological charge, S = 2. Numerical results demonstrate that such solitons indeed exist -naturally, with larger powers than their counterparts with S = 1. They can be found, at least, at P > 200. Here, we fix P = 250 and consider characteristics of the vortex solitons with S = 2, varying intersite coupling C. The minimum numbers of sites necessary for constructing the double-vortex modes of the IC and OC types are N = 6 (with N core = 2) and N = 7 (with N core = 3), respectively. Typical examples for 6 ≤ N ≤ 13 and 2 ≤ N core ≤ 5 (even and odd N and N core for the solitons of the IC and OC types, respectively) are displayed in Fig. 8, and dependences β(C) and ε(C) for stable vortex solitons with S = 2 and different values of N are displayed in Fig. 9.
The intrinsic symmetry of the vortex solitons of the OC type, with even vorticity S ≥ 2, predicts that its modal wave function in the central waveguide (n = 0), u 0 (x), is a real even function vanishing at x = 0. Indeed, an example displayed in Fig. 10(a) clearly corroborates the prediction. This feature may be compared to its counterpart in the case of S = 1, i.e., the odd real wave function u 0 (x), see Fig. 5(a) above. It is also worthy to stress that, while u 0 (x) virtually vanishes in a broad core area in Fig. 10(a), u 0 (x) keeps small positive values at all x = 0, i.e., it does not cross zero.
Similar to what is demonstrated above for S = 1, Fig. 9(a) shows that β gradually decreases with the increase of C, while ε(C) may strongly deviate from ε = 1, indicating strong anisotropy of vortex solitons with S = 2 [see in Fig. 9(b)]. The branches shown in Fig. 8 are completely stable, whereas their extensions, not shown in the figure, carry solitons which are unstable against splitting in the discrete direction. Typical examples of the evolution of the stable and unstable vortex solitons with S = 2 are displayed in Figs. 10(b,c). A comparison of values of the Hamiltonian between the vortex solitons with S = 1 and 2, which have equal powers, P , are shown in Figs. 9(c,d). The comparison indicates that, for the same N , the Hamiltonian is smaller for S = 1. Furthermore, the existence region for the vortices with S = 1 is larger (in some cases, much larger) than that for their counterparts with S = 2. Interestingly, two different stable branches are found for N = 8 with smaller and larger values of C, which are shown by the two red-color segments in Fig. 9(c). This result implies that there are two types of vortex soliton with S = 2 for N = 8. To illustrate this conclusion, Fig. 8(b) shows an example of the vortex state belonging to the branch with larger C and N core = 2. An example of the stable vortex soliton from the branch with smaller C is displayed in Fig. 10(d), with N core = 4.

D. Mobility of the semidiscrete vortex solitons
A nontrivial issue is mobility of the vortex solitons in the discrete direction, initialized by the application of a kick, with strength η, to them: where U (0) n (x) represents a quiescent soliton. Here, we address this issue for U n (x) taken as the stable vortex of the OC type with S = 1 and P = 100. At C < 0.2, the kick cannot set the soliton in progressive motion, just destroying it if η is too large, as shown in Fig. 11(a), or causing a finite leap, as shown in Fig. 11(b). The kicked vortex soliton demonstrates mobility at moderate discreteness, with C > 0.2, see an example in Fig. 11(c).
Collisions between two vortex solitons moving in the opposite discrete directions can be initialized by taking where 2n 0 is the initial separation between the solitons. Elastic collisions occur if η is small. In this case, the two colliding vortex solitons retain their vorticities after the collision. Inelastic collisions occur between the vortex solitons kicked with larger values of η. In this case, the collision splits each soliton in fragments. Figure 12 shows typical examples of elastic and inelastic collisions for the vortex solitons with (P, C, ±n 0 ) = (100, 0.3, ±16). In this case, the collision remains elastic at η < 0.016π.

IV. CONCLUSION
We have introduced the spatial-domain model for stacked set of tunnel-coupled planar waveguides with the combination of intrinsic self-focusing and defocusing cubic and quintic nonlinearities. The model applies as well to arrays of tunnel-coupled fiber waveguides with the same nonlinearity and anomalous group-velocity dispersion, which is a temporal-domain counterpart of the paraxial diffraction in planar guiding cores. Unlike fundamental semidiscrete solitons, previously studied in a similar model with the cubic-only nonlinearity, we here aimed to construct solitons with embedded vorticities, S = 1 and 2. It is found that such vortex solitons of the IC and OC (intersite-and onsite-centered) types, composed of N excited sites carrying non-negligible amplitudes in the transverse direction, form stable families, starting from minimum values (N IC ) min = 2 and (N OC ) min = 3, respectively. The semidiscrete vortex solitons of the OC type always feature an internal core, while the IC solitons with small N may have a coreless structure. The system admits multistability, i.e., coexistence of stable vortex solitons with equal norms and different values of N . The existence and stability of such modes are only possible when the quintic defocusing term is not too small. The vortex solitons feature mobility in the discrete direction, provided that the coefficient of the intersite coupling exceeds a certain minimum value. Semidiscrete vortex solitons moving in opposite directions collide elastically or inelastically if they are set in motion by relatively weak or strong kicks, respectively.