Vortex Structures in Saturable Media Vortex Structures in Saturable Media

In this chapter, dipole and vortex solitons are computed using spectral renormalization method in the focusing two-dimensional saturable nonlinear Schrödinger (SNLS) equation with periodic and quasicrystal lattice potentials. The nonlinear stability of these multi-humped solitons is investigated using direct simulations of the SNLS equation. It is shown that multiple vortex structures on quasicrystal lattices can be nonlinearly stable as the saturation and the external lattice may prevent the collapse. These results may have applications to investigations of localized structures in nonlinear optics and Bose- Einstein condensates.

Although these structures have not been truly observed in χ ð3Þ cubic media due to collapse, they appear as special numerical solutions of the focusing (2+1)D cubic nonlinear Schrödinger (NLS) equation with an external potential (lattice).
The stability of these solitons is important to applications. Rigorous stability theory of fundamental (positive) solitons is well established [9][10][11][12]. It is known that self-focusing media allow the stable propagation of one-dimensional (1+1)D spatial solitons, due to the balance between the linear diffraction and self focusing. However, two-dimensional (2+1)D optical solitons in Kerr media with the instantaneous cubic nonlinearity are unstable, due to the self-focusing (collapse) at high powers. Nevertheless, saturation of the nonlinearity may prevent the collapse, securing stable soliton propagation.
It has been previously shown that multi-humped soliton solutions of (2+1)D cubic NLS equation with a periodic external lattice can be stable in a certain parameter regime [17][18][19] and also in saturable (e.g. photorefractive) media [20][21][22]. However, the (computational) stability properties of multiple vortex solitons on background quasicrystal lattices have remained relatively unexplored. In Kerr media, in Ref. [23], it is computationally demonstrated that the twodimensional cubic nonlinear Schrödinger (NLS) equation with a quasicrystal lattice potential admits multiple dipole and vortex solitons. In aforementioned work, it is numerically shown that certain multiple vortex structures on quasicrystal lattices tend to be nonlinearly stable if the vortex/dipole humps are located on lattice minima but they are nonlinearly unstable if the humps are located on lattice maxima.
The analytic stability theory of optical modes on lattices that possess a higher topological complexity, for example, lattices containing defects, dislocations and/or quasicrystal lattices, is not sufficiently developed. On the other hand, this issue is studied numerically in recent years [24] and it is reported that dipole and vortex solitons for the two-dimensional nonlinear Schrödinger (NLS) equation with external potentials that possess strong irregularities, that is, edge dislocations and a vacancy defects, exist and it is numerically observed that these multi-humped modes in the defect lattices can be stable or unstable depending on the lattice properties.
In this study, we compute soliton solutions of the focusing (2+1)D saturable nonlinear Schrödinger equation with periodic and quasicrystal (e.g. Penrose) potentials and study their nonlinear stability. The governing equation is where Δu≡u xx þ u yy is the diffraction term, z plays the role of time (or distance). In optics, e 0 corresponds to the applied DC bias field inducing the saturable nonlinearity, uðx, y, zÞ is the complex-valued, slowly varying amplitude of the electric field in the xy-plane propagating in the z direction, Vðx, yÞ is an external optical potential (lattice) that can be formulated as the intensity of a sum of N phase-modulated plane waves [19], that is, where V 0 > 0 is constant and corresponds to the peak depth of the potential, that is, n is a wave vector defined by ðk n x , k n y Þ ¼ ½Kcosð2πn=NÞ, Ksinð2πn=NÞ. The external potential defined in Eq. (2) with N ¼ 2; 3; 4; 6 corresponds to periodic structures (crystals) and with N ¼ 5; 7 correspond to quasicrystals, which have a local symmetry around the origin but, different than that of periodic crystals, these types of structures are not invariant under spatial translations.
In this work, as the external potential, we consider the periodic and quasicrystal lattices corresponding to N ¼ 4 and N ¼ 5 in Eq. (2), respectively. In particular, the quasicrystal with N ¼ 5 is often called the Penrose tiling [25]. Contour images of these potentials are displayed all with V 0 ¼ 1 and K ¼ 1 in Figure 1. Freedman et al. observed solitons in Penrose and other quasicrystal lattices generated by the optical induction method [6].
In order to compute localized solutions (i.e., soliton solutions) to nonlinear evolution equations arising in optics, various techniques have been used, for example, shooting and relaxation techniques and the self-consistency method. A different method is introduced by Petviashvili [26] to construct localized solutions in the two-dimensional Korteweg-de Vries equation (usually referred to as the Kadomtsev-Petviashvili equation).
Petviashvili's method aims to transform the governing equation to Fourier space and determines a convergence factor based upon the homogeneity of a single nonlinear term. This method has been used to find localized solutions in a wide range of nonlinear systems in many works. This method can be successfully applied to nonlinear systems only if the degree of the nonlinearity is fixed in the related evolution equation. It is a well-known fact that in nonlinear optics, many equations involve nonlinearities with different homogeneities, such as cubicquintic, or even lack of homogeneity, as in saturable nonlinearity. A novel numerical scheme is proposed by Ablowitz and Musslimani [27] for computing solitons in nonlinear wave guides. This method is called spectral renormalization method and the idea behind the method is to transform the governing equation into Fourier space and find a nonlinear nonlocal integral equation coupled to an algebraic equation. The coupling prevents the numerical scheme from diverging. By this iteration scheme, the initial condition rapidly converges to an optical mode which is the numerical solution of the governing equation. This method can efficiently be applied to a large class of nonlinear wave problems including higher-order nonlinear terms with different homogeneities.
In this work, we use spectral renormalization method as explained below.

Spectral renormalization method
Spectral renormalization method is essentially a Fourier iteration method. The idea of this method was proposed by Petviashvili. Later, this method is improved by Ablowitz et al. and applied to (2+1)D NLS equation. In this subsection, the method is configured so that it can be applied to the saturable NLS equation.
Consider the saturable (2+1)D nonlinear Schrödinger equation with a potential in three-dimensional space: iu z ðx, y, zÞ þ u xx ðx, y, zÞ þ u yy ðx, y, zÞ − e 0 uðx, y, zÞ Using the ansatz uðx, y, zÞ Multiplying both sides of this equation by e iμz results in By applying Fourier transformation, one obtains To prevent singularities in the denominator in future calculations, the term rf is added to and subtracted from Eq. (6): Solving for thef in the second term above yieldŝ Making the substitution f ðx, yÞ ¼ λwðx, yÞ where λ is a non-zero constant gives When indexed, Eq. (10) can be utilized in an iterative method in order to find w. For this purpose,ŵ can be calculated using the following iteration scheme: with the initial condition taken as a Gaussian-type function (for the fundamental soliton) and the stopping criteria jw n −w n−1 j < 10 −8 . Here, the values of x 0 and y 0 define the location of the initial condition. In order to center the initial condition on the lattice maximum (the one that appears at the center of the lattice), one should, for example, take x 0 ¼ y 0 ¼ 0 and to center the initial condition on a lattice minimum (usually taken as one of the closest minima to the central maximum), one should take x 0 ¼ π, y 0 ¼ 0, for the periodic lattice N ¼ 4.
However, λ is unknown and hence must be calculated for each iteration Multiplying Eq. (10) by k 2 x þ k 2 y þ r leads to After moving all terms to the left side, one has Multiplying by the conjugate ofŵ, that is, byŵ Ã results in Finally, by integrating this equation, one gets This is nothing but an equation of the form FðλÞ ¼ 0 which can be solved for λ by employing the Newton-Raphson method, for instance. The Newton-Raphson method is a numerical method for finding roots of equations and is given by the following iteration scheme: where F ′ ¼ dF dλ and λ 0 are the initial guess for the root. Here λ 0 ¼ 1 and the stopping criterion is jλ n −λ n−1 j < 10 −8 .
In this work, we numerically find both multiple dipole and vortex solitons on periodic (N ¼ 4) and quasicrystal "Penrose" (N ¼ 5), background lattices. The nonlinear (in)stabilities are also examined for these localized structures by direct computations of Eq. (1). The initial conditions are taken to be a dipole or a vortex.

Numerical investigation of dipole and vortex solitons
In this section, we show the existence of both dipole and vortex solitons centered at the lattice maxima for both periodic and the Penrose potentials. Hereafter, the potential depth is set to V 0 ¼ 1, the saturation parameter is set to e 0 ¼ 8 and the propagation constant for all dipole and vortex structures is taken to be as μ ¼ 4.
For the spectral renormalization, we used the following initial condition, centered at the lattice maxima for both periodic and the Penrose potentials where x n , y n represent the location of vortex solitons, θ n is the phase difference, M corresponds to the number of humps and A is a positive integer.
Since the lattice solitons located on lattice maxima are found to be nonlinearly unstable for cubic nonlinearity case in Ref. [23], in this work, we deal with the dipoles and vortex solitons located on lattice maxima to demonstrate the stabilization effect of the saturation.

Dipole solitons
A dipole or a two-phase localized vortex is found for the periodic and Penrose lattices.
The dipole profile, its phase structure and the contour plot of the dipole humps superimposed on the underlying periodic potential are shown in Figure 2.
Dipoles at the lattice maxima for quasicrystal (Penrose) potentials are also found numerically; here A ¼ 3 and r ¼ 6:2504. The dipole profile, its phase structure and the contour plot of the dipole humps superimposed on the underlying Penrose lattice are shown in Figure 3.

Vortex solitons
Vortex solitons on lattice maxima including four-hump vortex solitons on a periodic lattice and both five-and 10-hump vortex solitons on a Penrose lattice were also investigated.   First, we consider the lattice which is periodic along the x and y directions which corresponds to N ¼ 4 in Eq. (2). Four-humped vortex solitons were obtained by taking the initial maxima of the periodic lattice, x n , y n , θ n to be Using the above initial condition, we numerically found a four-hump vortex centered at the periodic lattice maxima. The final vortex profile, the phase structure and the contour plot of the vortex solitons superimposed on the lattice maxima are shown in Figure 4.
Similarly, a five-hump vortex centered at the maxima of the Penrose potential shown in Figure 5 is obtained when θ n ¼ 2πn=5. Both the five-hump and the 10-hump vortices are found.
The 10-hump vortex located at the Penrose lattice maxima is shown in Figure 6; in this case θ n ¼ πn=5.

Nonlinear stability of vortex and dipole solitons
An important issue is the nonlinear stability of these vortex and dipole solitons. In order to examine the nonlinear stability of the vortex and dipole solitons found above, we directly compute Eq. (1), over a long distance (z ¼ 30 is typically found to be sufficient) for both types of potentials. The initial conditions were taken to be a multi-phased optical mode, namely a dipole or a vortex.
The nonlinear stability of the dipoles and the vortex structures are investigated by monitoring the maximum amplitude versus the propagation distance, the change in the location of centers of mass and the phase of the dipole/vortex solitons.
A dipole/vortex is assumed to be stable if it preserves 1. Its peak amplitude, as opposed undergoing self-focusing and/or finite-distance collapse.
2. Its position on the lattice, that is, be drift-stable (drift-unstable solitons are typically characterized by "humps" that drift from lattice maxima toward nearby minima).
The center of mass is calculated as juj 2 dxdy is the soliton power.
One should note that, for each case, the average displacement of the center of mass in the xdirection may change from negative to positive values depending on the initial locations of the dipole/vortex humps.  The nonlinear stability of vortex and dipole solitons is investigated in the following subsections separately.

Nonlinear stability of dipole solitons
In this subsection, the nonlinear stability properties of dipole solitons on periodic and Penrose lattice maxima are investigated.
The nonlinear stability properties of the dipole solitons on periodic lattice maxima are demonstrated in Figure 7. As the initial condition, we took the dipole solitons at periodic lattice maxima, as shown in Figure 2. The maximum amplitude, the location of the centers of mass of dipole solitons on periodic lattice maxima versus the propagation distance z ¼ 30 are plotted in Figure 7. The dipole profile at the final point z ¼ 30 is also depicted in the same figure. It is clearly seen that the dipole solitons on periodic lattice maxima are nonlinearly stable since the maximum amplitude oscillates with relatively small amplitude and the dipole solitons preserve their shape and location during the evolution.
The contour plots and the complex phase structures of this dipole are also depicted along the propagation distance z in Figure 8. This figure reveals that the dipole located on periodic lattice maxima preserves its shape and complex phase along the evolution distance.
Next, we investigate the nonlinear stability of dipole solitons on the Penrose lattice maxima that is previously shown in Figure 3. Similar to the dipole on the periodic lattice maxima, these structures are also found to be nonlinearly stable due to the conservation of the location and the phase structure and at the same time, only showing small oscillations in the maximum amplitude during the evolution (see Figures 9 and 10).
We should note that, for cubic nonlinearity case, in Ref. [23], the nonlinear stability properties of the dipole solitons on both periodic and the Penrose lattice maxima are demonstrated and both dipole solitons are found to be nonlinearly unstable since they exhibit strong localization after a few diffraction lengths and breakup in their phase structures. The Penrose lattice dipole solitons also suffer from the drift instability since the dipole humps both move from the lattice maxima toward nearby lattice minima immediately.

Nonlinear stability of vortex solitons
The nonlinear stability properties of the vortex solitons on periodic and Penrose lattice maxima are also investigated in this study.
Following the same order with that of dipoles on lattice maxima case, we start with the fourhump vortex solitons on periodic lattice maxima. Taking the four-hump vortex on periodic lattice maxima (see Figure 4) as the initial condition, the nonlinear stability properties are examined. In Figure 11, we plot the maximum amplitude and the location of the centers of   mass of the vortex humps versus the propagation distance z ¼ 30 and the vortex profile at the final point z ¼ 30. It is seen that the maximum amplitude of the vortex solitons oscillates with small amplitude. Furthermore, during the evolution, the vortex solitons do not move from their locations, meaning that there is no drift instability either.
The contour plots and the complex phase structures of this four-hump vortex are also depicted along the propagation distance z in Figure 12. This figure reveals that the four-hump vortex located on periodic lattice maxima perfectly preserves its shape and complex phase along the evolution distance. As a result of these facts, four-hump vortex solitons at periodic lattice maxima are found to be nonlinearly stable.  Previously in Ref. [23], it is numerically demonstrated that, for the cubic nonlinearity case, the maximum amplitude of the four-hump vortex solitons centered at the periodic lattice maxima increases sharply, after z ¼ 0:3 and due to this instant blow up in the maximum amplitude these vortex solitons are found to be nonlinearly unstable. It is shown that the phase structure also breaks up for the same vortex structure.
We investigated the nonlinear stability properties of the five-and 10-hump vortex solitons at the lattice maxima of Penrose potential by taking the five-( Figure 5) and 10-hump ( Figure 6) vortices on Penrose lattice maxima as initial conditions, respectively.
As can be seen from Figure 13, the maximum amplitude of these five-hump vortex solitons oscillate with relatively small amplitude and vortex humps stay at the same place during the direct simulation (no drift instability). The contour plots and the complex phase structures of this five-hump vortex are also depicted along the propagation distance z in Figure 14. This figure reveals that the five-hump vortex located on Penrose lattice maxima preserves its shape and complex phase along the evolution distance. As a result of these facts, vortex solitons on the Penrose lattice maxima are found to be nonlinearly stable.  Next, the nonlinear stability properties of vortex solitons on a Penrose lattice with 10 humps at lattice maxima are examined.
In Figure 15, we plot the maximum amplitude and the location of the centers of mass of a 10hump vortex on Penrose lattice maxima versus the propagation distance z ¼ 15. The contour plots and the complex phase structures of this 10-hump vortex are also depicted along the propagation distance z in Figure 16.  It is seen that the maximum amplitude of these vortex solitons increases somewhat until z ¼ 10 and after this value, the maximum amplitude increases sharply. On the other hand, although the vortex humps preserve their locations, the phase structure starts to break up and amalgamation of four-vortex humps (two by two) is observed at around z ¼ 10. At z ¼ 13, the phase structure is entangled and amalgamation of vortex humps is very clear resulting in the disappearance of four-vortex humps. At z ¼ 15, the phase structure is totally entangled and only four-vortex humps keep fairly acceptable amplitude and the others seem to disappear (see Figure 16). Therefore, as a result, we found the 10-hump vortex on Penrose lattice maxima as unstable.  In this case, it is clearly observed that the instability occurs due to amalgamation of vortex humps. In the Penrose lattice, as the number of the vortex humps increases, the humps get closer to each other which increases the possibility of amalgamation and this results in the instability. But still, the effect of saturation is very clear that the instability occurs later than that of cubic nonlinearity case. Again, one should note that, in previous cited work [23], for both cases, it is shown that the maximum amplitudes increase quite sharply after a short propagation distance (around z ¼ 0:8 and z ¼ 0:75 for five-and 10-hump cases, respectively) indicating nonlinear instability. The phase also breaks up these cases and moreover, during the evolution, both five-hump and 10-hump vortex solitons move from their initial locations (lattice maxima) to the lattice minima. The drift instability of lattice solitons, which is characterized by humps moving from the lattice maxima to the nearby lattice minima, was previously observed numerically for fundamental lattice solitons by Ablowitz et al. [28] (see also [14][15][16]).

Conclusion
We have numerically investigated dipole and multi-vortex structures associated with both periodic (N ¼ 4) and Penrose (N ¼ 5) lattices. First, we showed the numerical existence of vortex and dipole solitons on periodic and Penrose lattices. We investigated the nonlinear stability properties by simulations of saturable NLS equation. The simulations of the NLS equation showed that i. both the dipole and the vortex solitons centered at the periodic lattice maxima are found to be nonlinearly stable; ii. both the dipole and the five-hump vortex solitons centered at the Penrose lattice maxima are found to be nonlinearly stable; iii. 10-hump vortex solitons situated at Penrose lattice maxima are nonlinearly unstable. In this case, the nonlinear instability occurs as a result of the amalgamation of two or more vortex humps.
In conclusion, if one compares the cubic nonlinearity case to saturable nonlinearity case, saturation helps to suppress collapse for both dipole and vortex solitons located on lattice maxima. Even for a complex vortex structure such as 10-hump vortex on Penrose lattice maxima, collapse is suppressed and amalgamation of vortex humps and phase entanglement is delayed.

Author details
İlkay Bakırtaş Address all correspondence to: ilkayb@itu.edu.tr İstanbul Technical University, Faculty of Science and Letters, Department of Mathematics, Maslak, Istanbul, Turkey