**4. Higher-charged vortices in mixed linear-nonlinear circular arrays**

### **4.1. 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–64] or in purely nonlinear lattices [19, 65–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.

**Figure 9.** Stable and unstable propagation of azimuthons. (a, b) *n* = 5, *m* = 2, *p* = 30, *b* = 0.4. (c, d) *n* = 6, *m* = 4, *p* = 35, *b* = 0.4. White noise was added into (a).

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.

#### **4.2. Theoretical model**

14 Will-be-set-by-IN-TECH

The azimuthons aforementioned are restricted to the particular cases of *n* = *nJ*. In fact, azimuthon solutions can also be found when *n* � *nJ*. Contrary to the intuition and the cases in nonlocal media [53], the charge *m* of available azimuthon solutions is independent of the lattice order *nJ* 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

**4. Higher-charged vortices in mixed linear-nonlinear circular arrays**

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

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–64] or in purely nonlinear lattices [19, 65–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

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.

**Figure 9.** Stable and unstable propagation of azimuthons. (a, b) *n* = 5, *m* = 2, *p* = 30, *b* = 0.4. (c, d)

*n* = 6, *m* = 4, *p* = 35, *b* = 0.4. White noise was added into (a).

response in our model [31].

femtosecond laser pulses [61].

**4.1. Introduction**

of the power.

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*:

$$i\frac{\partial A}{\partial z} = -\frac{1}{2} \left( \frac{\partial^2 A}{\partial \mathbf{x}^2} + \frac{\partial^2 A}{\partial y^2} \right) - [1 - \delta \mathcal{R}(\mathbf{x}, y)] |A|^2 A - p\mathcal{R}(\mathbf{x}, y)A = 0 \tag{7}$$

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*) = Σ*n <sup>k</sup>*=<sup>1</sup> exp[−(*<sup>x</sup>* <sup>−</sup> *xk* )2/*d*<sup>2</sup> <sup>−</sup> (*<sup>y</sup>* <sup>−</sup> *yk*)2/*d*2], which means a circular array composed by *<sup>n</sup>* Gaussian waveguides with widths *d* = 1/2 and centers (*xk*, *yk*) 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*<sup>0</sup> = 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) admits several conserved quantities, including the power or energy flow *U* = <sup>∞</sup> <sup>−</sup><sup>∞</sup> <sup>|</sup>*A*<sup>|</sup> <sup>2</sup>*dxdy* and the Hamiltonian *H* = <sup>∞</sup> −∞[ <sup>1</sup> <sup>2</sup> <sup>|</sup> *<sup>∂</sup><sup>A</sup> ∂x* | <sup>2</sup> + <sup>1</sup> <sup>2</sup> <sup>|</sup> *<sup>∂</sup><sup>A</sup> ∂y* | <sup>2</sup> <sup>−</sup> *pR*|*A*<sup>|</sup> <sup>2</sup> <sup>−</sup> <sup>1</sup> <sup>2</sup>*γ*|*A*| <sup>4</sup>]*dxdy*.

We search for vortex solutions of Eq. (7) in the form of *A*(*x*, *y*, *z*)=[*qr*(*x*, *y*) + *iqi*(*x*, *y*)] exp(*ibz*), where *qr*(*x*, *y*) and *qi*(*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 *qr* and *qi*, 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).

**Figure 10.** (a) Power of vortex solitons with unit charge versus *b* for different *δ* values. (b) Hamiltonian versus power corresponding to (a). (c-f) Soliton profiles marked in (a) at *b* = 2.1 (c) and 3.5, 4.12, 4.60 in (d-f). In all panels *n* = 4, *p* = 6.

#### **4.3. 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 2*mπ* 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 2*mπ*/*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* → *blow* 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, due to the concentration of the light energy in regions where the nonlinearity is stronger. Interestingly, there exists a transition point of propagation constant *btr* at which the vortex soliton features a ring-shaped profile [Fig. 10(e)], which is analogous to the cases in lattice-free media or in media with imprinted Bessel lattices [17]. The nonlinearity modulation gives rise to unusual power-controlled soliton shape transformations [Figs. 10(d)-10(f)]. Note that the four lobes of vortex solitons at *b* > *btr* for *δ* = 1 reside in the regions between neighboring waveguide sites of the circular array [Fig. 10(f)], while the lobes of vortex solitons at *b* < *btr* reside in the waveguide sites [Fig. 10(d)]. Similar to the cases in harmonic lattices [56], vortex solitons may be termed as "on-site" vortices at *b* < *btr* and "off-site" vortices at *b* > *btr*.

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)].

16 Will-be-set-by-IN-TECH

δ=0.5

1

(b)

<sup>0</sup> <sup>13</sup> <sup>26</sup> <sup>39</sup> −420

**Figure 10.** (a) Power of vortex solitons with unit charge versus *b* for different *δ* values. (b) Hamiltonian versus power corresponding to (a). (c-f) Soliton profiles marked in (a) at *b* = 2.1 (c) and 3.5, 4.12, 4.60 in

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 2*mπ* 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 2*mπ*/*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* → *blow* and

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, due to the concentration of the light energy in regions where the nonlinearity is stronger. Interestingly, there exists a transition point of propagation constant *btr* at which the vortex soliton features a ring-shaped profile [Fig. 10(e)], which is analogous to the cases in lattice-free media or in media with imprinted Bessel lattices [17]. The nonlinearity modulation gives rise to unusual power-controlled soliton shape transformations [Figs. 10(d)-10(f)]. Note that the four lobes of vortex solitons at *b* > *btr* for *δ* = 1 reside in the regions between neighboring waveguide sites of the circular array [Fig. 10(f)], while the lobes of vortex solitons at *b* < *btr* reside in the waveguide sites [Fig. 10(d)]. Similar to the cases in harmonic lattices [56], vortex solitons may

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)].

be termed as "on-site" vortices at *b* < *btr* and "off-site" vortices at *b* > *btr*.

U

−280

−140

H

0

<sup>1</sup> 6.8 12.6 18.4 <sup>0</sup>

0.5

(a)

δ=1

b

(d-f). In all panels *n* = 4, *p* = 6.

*b* → ∞.

**4.3. Single-charged vortices**

13

26

U

39

**Figure 11.** (a) Power of vortex solitons with unit charge versus *b*. Soliton profiles at *b* = 2.50 (b), 5.70 (c) and 7.57 (d). Rings in (b-d) denote circular array of waveguides. In all panels *n* = 8, *δ* = 0.7, *p* = 6.

It increases firstly for small and moderate *b* and reaches a maximum, afterwards, the power decreases gradually and approaches to a constant ultimately. The initial increase of power indicates that the linear refractive index modulation plays a dominant role. However, the nonlinear contribution to the refractive index dominates when the power decreases with *b*. Thus, nonlinearity modulation imposes a restriction on the maximal power of vortex solitons. Hamiltonian of vortex solitons has inflexion points corresponding to the propagation constant value *b* = *bin* where *dU*/*db* = 0 [Fig. 10(b)]. Despite the fact that power is a decreasing function when *b* > *bin*, the peaks of vortices still increase with the propagation constant, which implies that the width of vortex lobes decreases rapidly for solitons at *b* > *bin*.

Next, we investigate the properties of vortex solitons with unit charge in circular waveguide array with *n* = 8. The vortex solitons bifurcate from the linear modes at *b* = *blow* (*U* = 0) and stop to exist at *b* = *bupp* [Fig. 11(a)]. The reason is that, for *b* > *bupp*, where the nonlinearity modulation absolutely dominates, the strong focusing nonlinearity results in vortex solitons with high peaks and narrow widths, which in return destroy the global staircase-like phase connection between the neighboring lobes. In other words, the vortices become eight independent components for *b* > *bupp*. Vortex solitons tend to shrink with the propagation constant [Figs. 11(b)-11(d)], which is quite different from the necklace solitons in the same regime, where necklace beams expand with the propagation constant [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.

**Figure 12.** (a) Power of vortex solitons with *m* = 1, *p* = 6, *n* = 8 versus *b* at different *δ* values. (b) Power of vortex solitons with different charges. (c) Instability growth rate of vortices with *m* = 1, *n* = 4 versus *b*. (d) Instability growth rate of vortices with different *m* versus *b*. Field modulus distributions of vortex solitons with *m* = 3 at *b* = 2.5 (e) and 7.17 (f). *n* = 8, *p* = 7 in (b-f) and *δ* = 0.7 except for (a, c).

#### **4.4. 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(a) displays the power curves of vortex solitons with unit charge in the waveguide array with *n* = 8 for varying nonlinearity modulation depth *δ*. Vortex solitons cease to exist at *b* = *bupp* if the nonlinearity modulation depth exceeds a certain value, which is in contrast to the vortices in linear lattices. The power curves at different *δ* share a common lower propagation constant cutoff. It means that all vortices bifurcate from a common linear mode of linearized version of Eq. (7). The existence domains of vortex solitons become narrower with the increase of *δ*.

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

18 Will-be-set-by-IN-TECH

**Figure 12.** (a) Power of vortex solitons with *m* = 1, *p* = 6, *n* = 8 versus *b* at different *δ* values. (b) Power of vortex solitons with different charges. (c) Instability growth rate of vortices with *m* = 1, *n* = 4 versus *b*. (d) Instability growth rate of vortices with different *m* versus *b*. Field modulus distributions of vortex

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(a) displays the power curves of vortex solitons with unit charge in the waveguide array with *n* = 8 for varying nonlinearity modulation depth *δ*. Vortex solitons cease to exist at *b* = *bupp* if the nonlinearity modulation depth exceeds a certain value, which is in contrast to the vortices in linear lattices. The power curves at different *δ* share a common lower propagation constant cutoff. It means that all vortices bifurcate from a common linear mode of linearized version of Eq. (7). The existence domains of vortex solitons become narrower with the increase of *δ*.

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

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

than those of *m* = 1 and 3 which may also be attributed to the group symmetry.

solitons with *m* = 3 at *b* = 2.5 (e) and 7.17 (f). *n* = 8, *p* = 7 in (b-f) and *δ* = 0.7 except for (a, c).

**4.4. Higher-charged vortices and their stability**

**Figure 13.** Stability (white) and instability (shaded) areas on the (*δ*, *b*) plane for vortex solitons at *m* = 1, *n* = 4 (a) and *m* = 3, *n* = 8 (b). Vortices in regions above *bin* suffer V-K instability. In all panels *p* = 7.

conclude that waveguide array with larger *n* tends to stabilize the vortices with charges close to *n*/2.

Typical examples of vortices with *m* = 3 at different *b* are displayed in Figs. 12(e) and 12(f). Contrary to the vortices with *m* = 1 shown in Figs. 11(b)-11(d), vortices with *m* = 3 expand with the increase of *b*. 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

**Figure 14.** Field modulus distributions of vortex solitons at *m* = 3, *b* = 2.1, *δ* = 0.5 (a) and *m* = 5, *b* = 5, *δ* = 0.7 (b). (c) Phase structure of (b). (d) Instability growth rate of solitons with different *m*. In all panels *p* = 7, *n* = 12.

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* > *bin*), 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* < *bin*.

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 *blow* 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].

**Figure 15.** Propagation simulations of stable (b-d, f) and unstable (a, e) vortex solitons. (a) Difference between a vortex soliton at *z* = 120 and *z* = 0. *m* = 1, *b* = 4.5, *n* = 4. Propagation results of vortex solitons at *m* = 2, *b* = 3.6, *n* = 8 (b), *m* = 3, *b* = 5.16, *n* = 8 (c), *m* = 3, *b* = 2.1, *n* = 12, *δ* = 0.5, *z* = 69 (e), and *m* = 5, *b* = 5, *n* = 12 (f). (d) Phase structure of (c). *δ* = 0.7 except for (e) and *z* = 512 except for (a, e). In all cases *p* = 7.

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 *blow* 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.
