Stable Higher-Charged Vortex Solitons in Optically-Induced Lattices

way for


Introduction
Vortices are fundamental objects which appear in many branches of physics [1] such as optics [2,3] and Bose-Einstein condensates [4]. In nonlinear optics, vortex solitons are associated with the phase dislocations (or phase singularities) carried by the nondiffracting optical beams [5], and share many common properties with the vortices observed in other systems, e.g., superfluids and Bose-Einstein condensates [6,7]. In a homogeneous medium, stable vortex solitons were proposed to exist in the so-called cubic-quintic or other similar nonlinear media, for example, combination of χ (2) and χ (3) nonlinear media, based on competing self-focusing and self-defocusing nonlinearities [8][9][10][11]. However, the experimental realization of vortex solitons in such media is hard, as the requirement of very high energy flow of light usually excites other higher-order nonlinearities, which may be dominant and suppress the occurrence of competing nonlinearities.
Thus far, dynamics of higher-charged vortex solitons are still poorly understood. The optical settings allowing stable higher-charged vortex solitons are rare. Main efforts in bulk or lattice-modulated nonlinear media were devoted to the analysis of vortex solitons with charges less than or equal to two. In the following sections, three different schemes for the realization of vortex solitons with higher charges will be addressed.
In section 2, the dynamics of vortex solitons in a radial lattice with a lower-index defect covering several lattice rings is revealed. The defect scale can be utilized to control the energy flow of vortices. Vortex solitons with various charges are stable in a region near the upper cutoffs of propagation constant. Although higher-charged vortices at higher energy flow suffer oscillatory instability, they can survive very long distances without visible distortions. Vortex solitons at lower or moderate energy flow are completely stable under appropriate conditions. The variation of topological charges slightly influences the existence and stability domains of vortex solitons. This property provides an effective way for the experimental realization of vortex solitons with higher charges in an optical setting with fixed parameters.
In section 3, the existence, stability and propagation dynamics of vortex solitons in a defocusing Kerr medium with an imprinted azimuthally modulated Bessel lattice are investigated. Since the special amplitude distribution of vortex soliton resembles to the azimuthons stated in [31] and the phase distribution is also a staircase function of the polar angle, such vortex solitons can also be termed as "azimuthons". The azimuthal refractive index modulation admits stable vortex solitons with lower or higher topological charges. The "stability rule" of azimuthons in defocusing cubic media is exactly opposite to that of vortex solitons in focusing media with the same transverse refractive index modulation. It is the first example of stable azimuthons in local nonlinear media.
In section 4, the stability of vortex solitons supported by a circular waveguide array with out-of-phase modulation of linear and nonlinear refractive indices is studied. The out-of-phase competition between two effects substantially modifies the stability properties of vortex solitons. Vortex solitons undergo remarkable power-dependent shape transformations. They expand or shrink radially with the propagation constant, depending on the phase difference between the neighboring lobes. In particular, we revealed that increasing waveguide number of circular array can stabilize vortex solitons with higher topological charges.

Introduction
Defects and defect states exist in a variety of linear and nonlinear systems, including solid state physics, photonic crystals, and Bose-Einstein condensates. When lights propagate in an optical lattice with a local defect, the band-gap guidance results in the formation of linear or nonlinear defect modes [32,33]. Recently, defect guiding phenomena of light in diverse settings, such as photonic crystals [13], fabricated waveguide arrays [34,35], and optically induced photonic lattices [36][37][38][39][40][41][42], have been predicted theoretically and observed experimentally. Ye and his coworkers proposed that stable nonlinear modes can be trapped in a lower-index defect sandwiched between two optical lattices, or in the cylindrical core of a radial lattice [43]. The variation of defect scales, depths and shapes can be used to stabilize and reshape the fundamental, dipole and vortex solitons [20].
In this section, we reveal that the defocusing media with an imprinted radially symmetric lattice with a lower-index defect covering several lattice rings can support stable vortex solitons with higher charges under appropriate conditions. In contrast to the cases in competing media, vortex solitons can propagate stably at lower or moderate energy flow. In lattices with fixed depth and defect scale, vortex solitons are completely stable provided that the propagation constant exceeds a critical value. In particular, we illustrate that the variation of topological charges slightly influences the existence and stability domains of vortex solitons. This is in contrast to all previous studies and allows the experimental realization of vortex solitons with higher charges without changing the parameters of an optical setting.

Theoretical model
We consider light propagation along the z axis of a defocusing Kerr medium with an imprinted transverse modulation of the refractive index. Dynamics of the beam can be described by the nonlinear Schrödinger equation for the dimensionless complex field amplitude A: Here, the longitudinal z and transverse x, y coordinates are scaled in the terms of diffraction length and beam width, respectively; p denotes the lattice depth; the refractive-index profile is given by R(x, y) = cos 2 (Ωr) for r ≥ (2N − 1)π/(2Ω) and R(x, y) = 0 otherwise, where r = (x 2 + y 2 ) 1/2 is the radial distance, Ω is the frequency, and N = 1, 2... is the number of rings removed from the lattice and characterizes the defect scale. Thus, the transverse modulation of refractive index features a lower-index guiding core. By comparing the defocusing bulk media without external potentials, the radial lattices with defects can confine the beams in a local region. An example of such refractive-index landscapes is shown in Fig. 1. Although there are defects in radial lattices, the wings of nonlinear modes still penetrate into the bulk of lattices. Thus, the existence of nonlinear modes strongly depends on the transverse lattices. Since the term 1/rd/dr in Laplacian can be neglected at r → ∞, the band-gap structure of a radially symmetric lattice is slightly different from that of 1D periodic lattice [43]. Thus, it is convenient to use the band-gap structure of 1D periodic lattice to approximately analyze the existence of solitons. Due to the fact that nonlinear modes in defocusing Kerr media can only be found in the finite gaps, we are interested in the solitons residing in the first finite gap. Equation (1) conserves several quantities, including the energy We search for stationary solutions of Eq. (1) by assuming A(x, y, z) = q(r) exp(ibz + imφ), where q is a r-dependent real function depicting the profile of stationary solution, b is a propagation constant associating with the energy flow, and m is an integer known as the topological charge of vortex soliton. The nonlinear mode degenerates to a fundamental radially symmetric mode when m = 0. The substitution of the light field into Eq. (1) yields: which can be solved numerically by means of a Newton iterative method. Mathematically, various families of stationary solutions are determined by the propagation constant b, lattice depth p, modulation frequency Ω and defect scale N. We vary b, p, N and fix Ω ≡ 2 in this section.
The stability of solitons can be analyzed by considering the perturbed solution in the form of: , here the perturbation components u, v could grow with a complex rate λ during propagation, and n is an integer representing the angle dependence of the perturbation and is termed as an azimuthal index. The substitution of the perturbed solution into Eq. (1) results in a system of eigenvalue equations: The coupled equations can be solved by a finite-difference method. In Cartesian coordinates, the square of the above linearization operator is self-adjoint if the stationary solutions are angle independent (fundamental solitons). Thus, the discrete eigenvalue is either purely real or purely imaginary. The instability growth rates with purely real parts correspond to the Vakhitov-Kolokolov (V-K) instability [44]. When the stationary solutions are angle dependent (vortex solitons), the eigenvalues may have both real and imaginary parts associating with an oscillatory instability. Stationary solutions are completely stable provided that all real parts of eigenvalues equal zero.

Discussions
First, we consider vortex solitons with unit charge supported by the defocusing Kerr media with an imprinted radial lattice with a defect. Without loss of generality, we set the defect scale N = 10 in the following discussions. In contrast to the fundamental solitons, the energy flow U of vortices is always a monotonically decreasing function of propagation constant b [ Fig. 2(a)]. We stress that the vortex solitons here are bright or ring-profile ones whose amplitudes decay to zero at infinity. Such vortices cannot be bifurcated from the dark vortices supported by the bulk defocusing media in the vanishing lattice. In other words, vortex soliton only exists when the lattice depth exceeds a critical value. For example, as shown in Fig. 2(b), the threshold value of lattice depth for the appearances of vortices with unit charge is p th ≈ 1.18, below which no vortex solutions can be found. Therefore, the vortex in the present lattice system is not a continuum of dark vortex in the vanishing lattice case, and it belongs to a different soliton family.  We summarize the linear-stability analysis results in Fig. 2(b). We show the critical value of propagation constant b n=1 cr above which no perturbations with the azimuthal index n and nonzero real part of growth rate were found. Vortex solitons are dynamically stable in a broad region near the upper cutoffs of propagation constant. It is the combination of defocusing nonlinearity and confining potential who affords the stability of vortex solitons. The precise structure of instability regions (patched) is rather complicated. There may exist multiple narrow stability windows. Now, we focus on the vortex solitons with higher topological charges in radially lattices with defects imprinted in a defocusing Kerr medium. Figure 4  The instability of vortex solitons with higher charges usually depends on the azimuthal index n [17,26,45]. Linear-stability analysis results reveal that for vortices with m = 3, the instability area associating with n = 2 is always dominant. For vortex solitons with m = 3 in a lattice with p = 5, the widths of instability windows associating with azimuthal indices n = 1, 2 and 3 occupy ≈ 26.77%, ≈ 38.69% and ≈ 16.93% of the width of the whole existence domain, respectively. An example of instability growth rate corresponding to azimuthal index n = 2 versus propagation constant is illustrated in Fig. 4(d). It indicates that vortex solitons suffer weak azimuthal instability, which allows them to propagate without obvious shape distortion over large propagation distances. Vortex solitons will be completely stable provided that the propagation constant exceeds a critical value.  vortices with m = 1, which constitutes one of our central results. That is to say, the stability area is slightly affected by the growth of topological charge, which allows one to realize stable vortex solitons with even higher charges. Since vortices with different charges share a collective stability area, one can input beams with different charges to excite vortex solitons with corresponding charges in certain parameter windows without changing the lattice depth, defect scale, modulation frequency, etc. It should be noted again that vortex solitons with different charges can propagate stably at lower or moderate energy flow, which is in sharp contrast to the cases in competing media, where very high energy flow is needed to stabilize the vortices [8,10]. Thus, in addition to the Bessel lattice [17], the radial lattice with a defect is another effective alternative for the realization of stable vortex solitons at lower or moderate energy flow, especially for vortices with higher charges.  photons decreases with the increase of rotational energy. The conclusion may be generalized to vortex solitons with continuous intensity distributions in other models.
The above discussions can also explain the decrease of the thickness of vortex solitons shown in Fig. 5(b). With the growth of topological charge, the decrease of effective mass of the beam is in companion with the increase of effectively rotational radius and angular velocity. For fixed propagation constant, the delocalization of vortex soliton weakens with the growth of topological charge. A representative propagation example of unstable vortex solitons is illustrated in Fig. 5(c). The vortex can propagate without visible shape distortion over hundreds of diffraction lengths.

Introduction
Besides the harmonic lattice, there is another important optical lattice with unique symmetry, the Bessel lattice, which can be created by nondiffracting Bessel beams with cylindrical symmetry. Kartashov et al. systematically investigated the dynamics of various types of solitons supported by the Bessel lattices, including multipole-mode solitons [46], ring-profile vortex solitons [17], spatiotemporal solitons [27] etc. Necklace [47], broken ring solitons [48] can also be trapped stably in different order Bessel lattices. For a review of the early works, see [23].
Interestingly, Bessel lattices with azimuthal modulation are also possible [25,49]. Such lattices resemble highly nonlinear micro-structured fibres [50] and may be realized in experiment by several incoherent Bessel beams with different intensities and orders [49]. The complex lattices can also be created in photorefractive crystals by the phase-imprinting technique [50,51]. The azimuthally modulated lattices exhibit several discrete guiding channels of linear refractive index. Stable soliton complexes and azimuthal switching in focusing cubic media with modulated Bessel lattices were reported in [49]. Neighboring components in soliton complex are out-of-phase. Ring-shaped and single-site solitons were observed in azimuthally modulated lattices [50,51]. Especially, by using group-theory techniques, Kartashov et al. derived a general "charge/stability rule" for vortex solitons supported by the azimuthal Bessel lattice [25].
In Ref. [31], Desyathikov and his coworkers introduced a novel class of spatially localized self-trapped ringlike singular optical beams in focusing cubic and saturable media, the so-called "azimuthons". The amplitude of such states is a spatially localized ring modulated azimuthally, and the phase of the azimuthon is a staircase function of the polar angle. This concept provided an important missing link between the radially symmetric vortices and rotating soliton clusters [52]. Following this work, stable azimuthons in nonlocal nonlinear media were found when the nonlocality parameter exceeds a certain threshold value [53,54].
However, stable azimuthons are only found in media with nonlocal responses. Azimuthons in local nonlinear media unavoidably experience azimuthal modulation instability upon propagation. In this section, we elucidate the existence and stability properties of azimuthons (vortex solitons with special amplitude distribution) supported by the azimuthally modulated Bessel lattices. It is the combination of nontrivial phase and lattice confinement who affords the existence of azimuthons. Thus, the azimuthons we obtained provide a missing link between the radially symmetric vortices and nonrotating soliton clusters although they break the radial symmetry due to the potential we used. Similar to the "azimuthons" stated in [31], the nonlinear localized modes we discussed can also be attributed to the two contributions induced by the internal energy flow and the modulated beam respectively. In sharp contrast to the cases in focusing cubic media [25], we reveal that the "stability rule" in defocusing cubic media is quite the reverse. The result is in good agreement with the conclusion given by [55] where the stability of discrete vortex solitons supported by hexagonal photonic lattices in focusing media is opposite to the stability in the defocusing one, though the discussions were limited to the single-charged and double-charged vortex solitons.

Theoretical model
We consider beam propagation along the z axis in defocusing cubic media with an imprinted transverse refractive index modulation. The dynamics of the nonlinear modes supported by such a scheme can be described by Eq. (1) Here, the parameter p describes the lattice depth. The profile of the modulated lattice is given by R(x, y) = J 2 n J [(2b lin ) 1/2 r] cos 2 (nφ), where n J denotes the order of the Bessel function, φ is the azimuthal angle, n stands for the azimuthal index and b lin defines the transverse lattice scale. Typical transverse linear refractive index modulation induced by the first order Bessel lattice with azimuthal index n = 2 is shown in Fig. 6(a). The local lattice maxima situated closer to the lattice center are more pronounced than others. The number of guiding channels in the main ring is given by 2n. Experimentally, Eq. (1) can be realized by launching a modulated Bessel beam into a photorefractive crystal in the ordinary polarization direction and a soliton beam in the extraordinary polarization direction [46]. In the particular case of optical lattice induction in SBN crystal biased with dc electric field ∼ 10 5 V/m, for laser beams with 10μm the propagation distance z ∼ 1 corresponds to 1mm of actual crystal length, while amplitude q ∼ 1 corresponds to peak intensity about 50 mW/cm 2 [25]. Note that Eq. (1) can also be treated as Gross-Pitaevskii equation for a 2D Bose-Einstein condensate with repulsive interatomic interactions trapped in an optical lattice created by an azimuthally modulated Bessel beam.
We search for stationary solutions of Eq. (1) in the form of A(x, y, z) = [q r (x, y) + iq i (x, y)] exp(ibz), where q r and q i are real and imaginary parts of the solution profiles and b is a nonlinear propagation constant. The twisted phase structure of the stationary solutions can be defined by m = arctan(q i /q r )/2π, where m is the so-called "topological charge" of vortex solitons. Substituting the expression into Eq. (1), we obtain: We fix b lin ≡ 2 and vary b, p and n without loss of generality.
To elucidate the stability properties of solitons, we search for perturbed solutions of Eq. (1) in the form A(x, y, z) = [q r + iq i + (u r + iu i ) exp(λz)] exp(ibz), where u r and u i are the real and imaginary parts of the perturbations, respectively. Substituting the perturbed solution into Eq. (1) and linearizing u r,i around q r,i yield a system of coupled Schrödinger-type equations for perturbation components u r,i : where u r,i may grow with a complex rate λ during the propagation of solitons. The eigenfunctions u r,i and eigenvalues λ can be solved numerically. The solitons are stable only when all real parts of λ equal zero.

Discussions
Before we discuss the dynamics of vortex solitons (azimuthons), it is important to understand the origin of such nonlinear modes. After removing the nonlinear term in Eq. (4), the linear equation has infinite eigenvalues and the corresponding linear eigen-modes. Nonlinear modes bifurcate from these linear modes while the nonlinearity cannot be ignored. Mathematically, the refractive index modulation contributed by the modulated Bessel lattice increases linearly with the growth of lattice depth. However, this relationship cannot hold for practical crystal when the lattice is modulated very deep. Thus, the practical realization of stable vortex solitons with higher topological charges becomes infeasible by solely increasing the lattice depth of the first-order lattice to a very large value. Fortunately, the higher-order Bessel lattice can suppress the azimuthal instability of vortex solitons effectively [26]. To study the properties of vortex solitons with higher charges, one must consider the higher order modulated Bessel lattices.
The following discussion will focus on azimuthons (vortex solitons) carrying different topological charges supported by azimuthally modulated Bessel lattices imprinted in defocusing cubic media. For the convenience of comparing with the results of Ref. [25], we assume R(x, y) = J 2 n [(2b lin ) 1/2 r] cos 2 (nφ), where the order of lattice equals to the azimuthal index. We also search for stationary solutions of azimuthons by a relaxation method. A Gauss beam multiplying a phase dislocation with charge m was selected as an initial iterative guess solution. Figure 7 displays some instances of azimuthons supported by the azimuthally modulated lattices with n = 4 and 6. The azimuthons exhibit spatially modulated patterns which are in contrast to the vortices in unmodulated Bessel lattices [17], where the vortices are ring-shaped. The number of amplitude peaks is determined by the azimuthal index n. Like the vortices in focusing media [25], azimuthons with similar amplitude distributions allow different topological charges. In the fourth-order Bessel lattices with azimuthal index n = 4, azimuthons can be found only for m = 1, 2 and 3. In sufficiently deep lattices with fixed b and p, the discreteness of azimuthons increases with the growth of the topological charge m, while the "radii" of the azimuthons are almost the same. Such properties are similar to the vortex solitons in focusing media [25]. For fixed p and n, the azimuthons will expand to the outer lattice rings at small b and shrink to the main guiding lattice ring at larger b. The local minima of azimuthons around the lattice ring approaches to zero when b → b co .
We also find azimuthons supported by lattices with different azimuthal indices. Numerical study reveals that azimuthon solutions can be found only when the relation 0 < m < n is satisfied. The relation also holds for the vortex solitons in focusing cubic media [25]. The phase difference between the neighboring components is mπ/n, which differs from the vortex solitons in harmonic lattices [56] or necklace solitons in Bessel lattices [47].
The properties of azimuthons supported by azimuthally modulated Bessel lattice are summarized in Fig. 8. The power of azimuthons is a descending curve due to the defocusing nonlinearity [ Fig. 8(a)]. Azimuthon solutions cannot be found when the propagation constant exceeds a certain value which corresponds to an eigen-value of the linearized equation of Eq. (4). The upper propagation constant cutoffs of vortex solutions with m = 1 and 3 are displayed in Fig. 8(b) and 8(c). The existence areas expand with the growth of lattice depth for fixed topological charge and shrink with the growth of topological charges for fixed lattice depth. There is a lower threshold lattice depth for supporting azimuthons. Comparing the points of b → 0 in Fig. 8(b) and 8(c), we find that the threshold lattice depth grows with the increase of topological charge m.
To comprehensively understand the stability properties of azimuthons supported by lattices with different depths and azimuthal indices, we performed the linear stability analysis on azimuthons in lattices with order (azimuthal index) n up to 10 and lattice depths p ≤ 80. We numerically derived an important "stability rule" for azimuthons in the azimuthally modulated Bessel lattice imprinted in defocusing media. That is, azimuthons might be stable only when the topological charge satisfies the condition:  Available charges and stability status  n=2  m=1 unstable  n=3  m=1 stable  m=2 unstable  n=4  m=1 stable  m=2 stable  m=3 unstable  n=5  m=1 stable  m=2 stable  m=3 unstable m=4 unstable  n=6 m=1 stable m=2 stable m=3 stable m=4 unstable m=5 unstable where n > 2. Azimuthons with m = n/2 for even n may be stable or unstable depending on the lattice parameters. There exists a narrow instability area near b → 0 when the lattice is modulated shallow (near its lower threshold value). For deeper lattices, completely stable azimuthons are possible. A summary of "stability rule" is presented in Table 1. The table shows the stability status of azimuthons for different lattice orders. It is exactly opposite to the Table I in [25], which was derived by the group-theory and is valid in the focusing cubic media. This finding also verifies the very recent reports [55,57] in which the stability of discrete vortex solitons supported by hexagonal photonic lattices in focusing media is proved to be opposite to the stability in the defocusing ones. We note that our conclusion is more general since the above two studies are restricted to the single-and double-charge discrete vortex solitons.

Lattice order
Linear instability analysis results of some unstable azimuthons supported by the fourth-, fifth-, and sixth-order lattices are shown in Fig. 8  The azimuthons aforementioned are restricted to the particular cases of n = n J . In fact, azimuthon solutions can also be found when n = n J . Contrary to the intuition and the cases in nonlocal media [53], the charge m of available azimuthon solutions is independent of the lattice order n J but less than the azimuthal index n. The initial input guess solutions with m ≥ n may converge to the nonlinear modes of the following three different categories: 1. an azimuthon with charge m < n; 2. a multipole mode or necklace soliton with neighboring components out-of phase; 3. a multipole mode or necklace soliton embedded into a global skew phase whose charge m = m − n. Thus, we conclude that azimuthon solutions can only be found for m < n. The reason may be attributed to the Kerr media with a local nonlinear response in our model [31].

Introduction
Recent studies demonstrated that the transverse nonlinearity modulation of optical materials can substantially affect the existence conditions and stability properties of spatial solitons [58]. Current technologies allow one to realize both the linear refractive index and the nonlinearity modulation of materials, for example, in photonic crystals with holes infiltrated with highly nonlinear materials [59]. Nonlinearity can also be modulated by changing the concentration of dopants upon fabrication [60] or be achieved in arrays written in glass by high-intensity femtosecond laser pulses [61].
Very recently, various types of one-dimensional solitons in the form of odd, even, dipole, triple, vector solitons, and two-dimensional solitons in the form of fundamental, multi-pole, vortex, surface solitons, and three-dimensional solitons in the form of light-bullets and optical tandems were predicted in competing linear and nonlinear lattices [18,[62][63][64] or in purely nonlinear lattices [19,[65][66][67]. On the other hand, spatially modulated nonlinearity may result in controllable soliton shape transformations [18,62]. Bound states with an arbitrary number of solitons can be found in media with spatially inhomogeneous nonlinearities [68]. The power-dependent location of stationary solitons and their stability in linear-nonlinear lattices were analyzed in Refs. [69], where the position and stability of the solitons become functions of the power.
Besides models with periodically modulated nonlinearity, settings with local nonlinearity modulation were also suggested [70]. In particular, Kartashov et al. demonstrated that necklace solitons can be supported by circular waveguide arrays with out-of-phase modulation of nonlinearity and linear refractive index [63]. They revealed that the stability domain of necklace solitons expands with the increase of the number of necklace spots. Although vortex solitons in linear-nonlinear or in purely nonlinear lattices were studied in [18,19], the discussions were limited to the vortices with unit charge. In fact, stable localized vortex solitons with charges higher than two were only predicted in media modulated by Bessel-like lattices [25,26]. In addition, the stability domain of vortex solitons in all previous studies shrinks rapidly with the growth of topological charge.
In this section, following the theoretical model proposed by Kartashov and his coworkers [63], we address the existence and stability properties of vortex solitons with higher charges supported by a circular array with mixed out-of-phase linear-nonlinear refractive index modulation. We find that vortex solution can be found only when its topological charge is less than half of the number of waveguides. Vortex solitons expand radially with the propagation constant when the phase difference between neighboring lobes is greater than π/2 and shrink radially when the phase difference between neighboring lobes is smaller than π/2. For vortex solitons with fixed charges, the stability weakens with the growth of waveguide sites. In particular, we demonstrate that vortex solitons with higher charges tend to be more stable than those with lower charges, which in sharp contrast to the cases in uniform or periodical lattice modulated media.

Theoretical model
We consider a beam propagation along the z axis in a circular waveguide array with out-of-phase modulation of linear refractive index and nonlinearity coefficient. Evolution of the nonlinear waves is governed by the nonlinear Schröinger equation for the dimensionless field amplitude A: where p and δ are the depths of modulation of linear refractive index and nonlinearity. We adopt modulation function of the linear refractive index and nonlinearity as R(x, y) , which means a circular array composed by n Gaussian waveguides with widths d = 1/2 and centers (x k , y k ) located on a ring of radius nr 0 /2 [63]. To guarantee that waveguides do not overlap and the soliton components residing in neighboring waveguides can interact with each other, we select r 0 = 0.6. The nonlinear coefficient γ = 1 − δR attains its minima where the linear refractive index attains maxima. Thus, the nonlinearity modulation and linear refractive index modulation are out-of-phase. The length of the arc between adjacent waveguides remains the same for any n since the radius of the array increases linearly with n. Equation (1) We search for vortex solutions of Eq. (7) in the form of A(x, y, z) = [q r (x, y) + iq i (x, y)] exp(ibz), where q r (x, y) and q i (x, y) are real and imaginary parts of the solution profiles and b is a nonlinear propagation constant. Substituting the expression into Eq. (7) leads to two coupled ordinary differential equations in terms of q r and q i , which can be solved numerically by a two-dimensional relaxation algorithm. The stability characteristics of vortex solitons can be understood by numerically solving the eigenvalue of Eqs. (5).

Single-charged vortices
First, we address the properties of vortex solitons supported by a circular waveguide array with n = 4. Typical vortex solitons with unit charge are illustrated in Fig. 10. Vortex solitons reflect their distinctive modulated profiles. The field modulus distribution of such state is a spatially localized ring modulated azimuthally, and the phase changes 2mπ along a closed trajectory including the beam centre [Figs. 10(c)-10(f)]. Thus, they look like the so-called "azimuthons" in local and nonlocal nonlinear media [31,53]. Vortex solitons are composed of n lobes located in different waveguides of the array with phase difference 2mπ/n between the neighboring lobes. This phase difference determines the net force contributed by the interactions of vortex lobes. There are n local intensity minima connecting the n lobes for small and moderate propagation constant. The discreteness of solitons strengthens for b → b low and b → ∞.
If the modulation depth of the local nonlinearity is small (δ ≤ 0.5), the lobes of vortices stay around sites of linear waveguide array. However, the competition between the linear and nonlinear refraction may result in a remarkable shape transformation at δ > 0. 5 Fig. 10(d)]. Similar to the cases in harmonic lattices [56], vortex solitons may be termed as "on-site" vortices at b < b tr and "off-site" vortices at b > b tr .
The power of vortex solutions is a nonmonotonic function of propagation constant, due to the competition between the linear and nonlinear refractive index contributions [ Fig. 10(a)].  [63]. It can be interpreted in physics that the phase difference between the neighboring lobes is π/4 which results in an attractive force between the neighboring components of vortices. The net attractive force strengthens with the propagation constant which leads to the radial shrinkage of vortices. In fact, the positions of four lobes of vortex solitons in waveguide array with n = 4 do not vary with the propagation constant because the phase difference between the neighboring lobes is π/2.

Higher-charged vortices and their stability
The stability of vortex solitons, especially for higher-charged vortices, is always an important problem. The following discussions will focus on the stability of vortex solitons with different charges in circular waveguide arrays with different n and δ values. Figure 12 The power curves of vortex solitons with different charges at fixed δ are shown in Fig. 12(b). Due to the discrete symmetry group of the circular waveguide array, vortex solutions can only be found when the condition 0 < m < n/2 is satisfied, just similar to the cases in azimuthally modulated Bessel lattices [25]. The existence domain of vortex solitons with m = 2 is narrower than those of m = 1 and 3 which may also be attributed to the group symmetry.
For vortex solitons with unit charge in a waveguide array with n = 4, the stability area shrinks quickly with the growth of nonlinearity modulation depth δ [ Fig. 12(c)]. However, vortex solitons with unit charge in waveguide array with n = 8 are unstable in their entire existence domain. Interestingly enough, vortex solitons with higher topological charges can propagate stably in the same configuration [ Fig. 12(d)]. The stability analysis results shown in Fig. 12(d) suggest that vortex solitons in a circular array can be stable only when the condition n/4 ≤ m < n/2 is satisfied. Furthermore, the stability area of vortex solitons with m = 3 is broader than that of vortices with m = 2. That is to say, vortex solitons with higher charges are more stable than those with lower charges which constitutes one of our central findings. We As can also be explained by the phase difference between the neighboring lobes. The phase difference between the neighboring lobes of vortices with m = 1, 2 and 3 is π/4, π/2 and 3π/4 respectively. The net force contributed by the eight lobes is attractive for vortices with m = 1, zero for vortices with m = 2 and repulsive for  vortices with m = 3. Since the peaks of lobes grow with the propagation constant, the net force becomes stronger and the lobes cannot be trapped on the waveguide sites, which leads to the shrinkage of vortices with m = 1 and expansion of vortices with m = 3. The lobes of vortices with m = 2 do not change their positions with the variation of b, just similar to the vortices with m = 1 in waveguide array with n = 4 [ Fig. 10]. Thus, we draw a conclusion that the expansion or shrinkage of vortices is solely determined by the radio between the topological charge and the number of waveguides.
To comprehensively understand the stability characteristics of vortex solitons, we conduct linear stability analysis on the stationary solutions with different m in circular waveguide array with different n for varying nonlinearity modulation depth δ. Representative stability and instability domains on the (δ, b) plane are illustrated in Fig. 13. At larger b where dU/db < 0 (b > b in ), vortices are expected to be linearly unstable according to the Vakhitov-Kolokolov (V-K) criterion [44]. This is in good agreement to the numerical analysis results. However, other types of instability (e.g. oscillatory instability with complex growth rates) may arise when b < b in .
For vortex solitons with m = 1, n = 4 and m = 3, n = 8, the stability domains shrink with the growth of nonlinearity modulation depth. Note again vortex solitons with m = 1, n = 8 are completely unstable. Comparing Fig. 13(a) with Fig. 13(b), one can immediately find that the stability area of vortices with m = 3, n = 8 is broader than that of m = 1, n = 4. Thus, increasing waveguide number can significantly suppress the instability of vortex solitons with higher charges. There is an instability band near b low for vortex solitons with m = 2, n = 8 [ Fig. 12(d)]. We stress that vortex solitons with higher topological charges are stable in a wide parameter window which is hardly realized in other settings except in Bessel-like lattice modulated media [25,26]. To shed more light on the vortex solitons in circular waveguide array, we explore the stationary solutions, stability and propagation dynamics of vortices in a circular array with larger n (e. g. n = 12, 18 etc.). Figures 14(a) and 14(b) show the field modulus distributions of vortex solitons at m = 3, δ = 0.5, b = 2.1 and m = 5, δ = 0.7, b = 5 in waveguide array with n = 12 respectively. As mentioned above, the stability of vortex solitons weakens with the increase of nonlinearity modulation depth. Yet, vortex solitons with m = 3, δ = 0.5 are stable in a very narrow area near b low while vortex solitons with m = 5, δ = 0.7 are stable in a wide parameter window [ Fig. 14(d)]. It reveals again that vortex solitons tend to be more stable when m approaches to n/2.
Finally, we perform extensive propagation simulations of vortex solitons by the split-step Fourier method to verify the linear stability analysis results. Typical stable and unstable propagation examples are shown in Fig. 15. Figure 15(a) plots the difference of vortex profiles at z = 120 and 0. One can infer from the simultaneous increase of four lobes that the soliton experiences a V-K instability, that is, the unstable eigenvalues of λ are purely real, which do not destroy the phase distribution. On the other hand, vortex soliton at m = 3, b = 2.1, n = 12, δ = 0.5 suffers oscillatory instability (growth rate with both real and imaginary parts) [ Fig. 15(e)]. In the numerical simulations, we add white noises into the initial inputs for stable vortices and add no noises for unstable vortices.