Azimuthons in weakly nonlinear waveguides of different symmetries

We show that weakly guiding nonlinear waveguides support stable propagation of rotating spatial solitons (azimuthons). We investigate the role of waveguide symmetry on the soliton rotation. We find that azimuthons in circular waveguides always rotate rigidly during propagation and the analytically predicted rotation frequency is in excellent agreement with numerical simulations. On the other hand, azimuthons in square waveguides may experience spatial deformation during propagation. Moreover, we show that there is a critical value for the modulation depth of azimuthons above which solitons just wobble back and forth, and below which they rotate continuously. We explain these dynamics using the concept of energy difference between different orientations of the azimuthon.


Introduction
Spatial solitons are nonlinear localized states that keep their form during propagation due to the balance between diffraction and self-induced nonlinear potential [1]. Recently, there has been a lot of interest in a generalized type of spatial solitons, the so-called azimuthons. These are azimuthally modulated beams, that exhibit steady angular rotation upon propagation [2]. They can be considered as azimuthally perturbed optical vortices, i.e. beams with singular phase structure [2][3][4]. Theoretical studies demonstrated both, stable and unstable propagation of azimuthons [5][6][7], and the first experimental observation of optical azimuthons was recently achieved in rubidium vapors [8].
It appears that higher order solitonic structures and optical vortices are generally unstable in typical nonlinear media with local (Kerr-like) response [9,10]. On the other hand, it has been shown that a spatially nonlocal nonlinear response provides stabilization of various complex solitonic structures including vortices [11][12][13][14]. Consequently, azimuthons and their dynamics have been studied almost exclusively in the context of spatially nonlocal nonlinear media [7,[15][16][17]. In spite of the fact that there are various physical settings exhibiting nonlocality such as nematic liquid crystals [18][19][20], Bose-Einstein condensates [4,21,22], plasmas [23], thermo-optical materials [24] etc., experimental realization of such systems is always quite involved. Moreover, from the theoretical point of view, nonlocal media are quite challenging for numerical modeling and analytical treatment.
In this paper we propose a much simpler and experimentally accessible optical system to study the propagation of azimuthons: a weakly nonlinear optical multi-mode waveguide. Here, weakly nonlinear means that the nonlinear induced index change, which is proportional to the intensity of the optical beam, is small compared to the index profile (or trapping potential) of the waveguide. Following [25], we can expect that weakly nonlinear azimuthons are stable in multi-mode waveguides.
The paper is organized as follows: In Sec. 2, we briefly introduce the general model equation for beam propagation in weakly nonlinear waveguides, and then we discuss in detail the properties of dipole azimuthons in circular and square waveguides in Sec. 3, respectively. In Sec. 4, higher order azimuthons are investigated, and in Sec. 5 we conclude.

Mathematical modeling
We consider the evolution of a continuous wave (CW) optical beam with amplitude E (ξ , η, ζ ), where (ξ , η) and ζ denote the transverse and longitudinal coordinates, respectively. Then the propagation of this beam in a weakly-guiding waveguide with Kerr nonlinearity in the scalar, slowly varying envelope approximation is described by the following equation [26]: where k 0 = 2πn b /λ 0 refers to the carrier central wave number in the medium; n 2 is the Kerr nonlinear coefficient, n(ξ , η) the linear refractive index distribution, n b the background index, and lim ξ ,η→∞ n(ξ , η) = n b . We consider a weakly-guiding waveguide, so both the linear and nonlinear induced index change |n − n b | and n 2 |E | 2 are small compared to the mean index n b , and at the same time n 2 |E | 2 ≪ |n − n b |, to guarantee a weak nonlinearity. From the mathematical point of view, Eq. (1) is a (2+1)-dimensional nonlinear Schrödinger (NLS) equation with linear potential representing the waveguide profile. For technical convenience, we rescale Eq. (1) to dimensionless quantities with x = ξ /r 0 , y = η/r 0 , z = ζ /(2k 0 r 2 0 ), and σ = sgn(n 2 ), where r 0 represents the transverse spatial extent of the waveguide. Then, the two-dimensional (2D) NLS equation for the scaled wave function ψ = k 0 r 0 2|n 2 |/n b E reads with an "attractive" bounded potential V = 2k 2 0 r 2 0 (n − n b )/n b , given by the spatial profile of the refractive index of the waveguide. Here we will consider propagation in step index waveguides with parameter V (x, y) given as follows V 0 is the height of the potential and determines the number of linear waveguide modes.
Considering the specific case of silica-made waveguides, typical values for the parameters are n b = 1.4, |n − n b | ≤ 9 × 10 −3 , n 2 = 3 × 10 −16 cm 2 /W at a vacuum wavelength of λ 0 = 790 nm, and those values do not vary much when we choose any λ 0 in the range from visible to near-infrared. The laser-induced-damage-threshold (LIDT) for synthetic fused silica is about 10 J/cm 2 for nanosecond pulses (according to Eq. (1) and Table 2 in Ref. [27]) corresponding to ∼ 10 10 W/cm 2 , up to which the model should be valid. For shorter, femtosecond pulses this damage threshold is 1000 times higher, but since we consider "CW" beams the lower value for nanosecond pulses is relevant. It is justified to neglect temporal effects on beam propagation, because the dispersion length is already of the order of kilometers for pulse durations of a few tens of picoseconds, much longer than propagation distances considered in this work. Throughout this paper, we choose V 0 = 1000, and find that r 0 ≈ 25.0 µm , which nicely corresponds to the radius of a standard circular multi-mode fiber [28]. In addition, in order to guarantee intensities smaller than the LIDT the condition |ψ| 2 < 1/3 should be fulfilled.
Note that in Eq. (2) σ is the sign of nonlinearity: σ = 1 (σ = −1) represents focusing (defocusing) nonlinearity. The main difference between weakly nonlinear waveguides with focusing and defocusing nonlinearity is that defocusing nonlinearity supports higher amplitudes of the wave function, and in general leads to more stable and robust configurations. In this paper, however, we consider the experimentally relevant case of a focusing nonlinearity.

Dipole azimuthons in circular waveguides
Azimuthons are a straightforward generalization of the usual ansatz for stationary solutions [2]. They represent spatially rotating structures and hence involve an additional parameter, the rotation frequency ω (see also [29]), so we seek approximate solutions of the form where r = x 2 + y 2 and φ the azimuthal angle in the transverse plane (x, y), U is the stationary profile, ω the rotation frequency, and κ the propagation constant. For ω = 0, azimuthons become ordinary (non-rotating) solitons. The simplest family of azimuthons is the one connecting the dipole soliton with the single charged vortex soliton [15]. A single charged vortex consists of two equal-amplitude dipole-shaped structures representing real and imaginary part of U. If these two components differ in amplitude, the resulting structure forms a "rotating dipole" azimuthon. If one of the components is zero we deal with the dipole soliton, which consists of two out-of-phase humps and does not rotate for symmetry reasons. In a first attempt, let us assume we know the radial shape of the linear vortex mode F(r) which is normalized according to π r|F(r)| 2 dr = 1. Then, using separation of variables, we consider the simplest so-called "rotating dipole" azimuthon with ansatz [16] where A is an amplitude factor, and 1 − B the azimuthal modulation depth of the resulting ringlike structure. Because we are operating in the weakly nonlinear regime, using linear waveguide modes as initial conditions for nonlinear (soliton) solutions is a quite good approximation. After plugging Eq. (3) into Eq. (2), multiplying by U * and ∂U * /∂ φ respectively, and integrating over the transverse coordinates we end up with [16] − κP + ωL z + I + N = 0, (5a) This system relates the propagation constant κ and the rotation frequency ω of the azimuthons to integrals over their stationary amplitude profiles: The first two quantities (P and L z ) have straightforward physical meanings, namely power and angular momentum. The integrals I and I ′ are related to the diffraction mechanism of the system, whereas N and N ′ account for waveguide and nonlinearity. Thus we can formally solve for the rotation frequency with these quantities and obtain After inserting Eq. (4) into Eq. (7), it turns out that only the nonlinear term contributes to the rotation frequency ω (see also Eq. (10) in Ref. [30]). In order to give an estimate for the rotation frequency, we use the linear stationary modal profile F(r) of a circular waveguide expressed in terms of Bessel functions of first kind J 1 and modified Bessel function of second kind K 1 : where C is a normalization factor such that π |F(r)| 2 rdr = 1. Thus, the analytically predicted rotation frequency is From the above equation, one can see that sense of rotation is opposite for focusing and defocusing nonlinearities. Moreover, as ω = 2k 0 r 2 0 · 2π/ℓ, one can find from Eq. (9) an expression for the physical distance ℓ over which the azimuthon in a circular waveguide performs one full rotation: For example, assuming propagation in a silica fiber with the parameters discussed earlier (see Sec. 2) and taking A = 0.4, B = 0.5, we obtain ℓ ≈ 2.4 m. Figure 1 shows the dependence of the azimuthon rotation frequency as a function of its amplitude A (left panel) and the modulation parameter B (right panel). Blue solid lines represent the analytical predictions and, in excellent agreement, red dots are obtained from numerical simulations. An exemplary propagation of the dipole azimuthon in circular waveguide rotating with ω ≈ 0.037 is illustrated in Fig. 2. The top iso-surface plot displays 3D rigid and continuous rotation of the azimuthon during propagation. The row below depicts the transverse intensity (left) and phase (right) distribution of the dipole azimuthon during propagation after rotating by π/2.

Azimuthon-like dipoles in square waveguides
Because of the lack of circular symmetry, azimuthons in the strict sense of Eq. (4) (preservation of shape during rotation with constant frequency) cannot exist in a square waveguide. However, as we will show below, the azimuthon-like behavior is still possible even though the beam propagation is accompanied by variation of the beam transverse intensity distribution. To set the initial condition, we use two linear degenerated orthogonal dipole modes D 1 , D 2 (as in the case of the circular waveguide), which are normalized according to |D 1,2 | 2 dx dy = 1, and superpose them as before to form the azimuthon-like object The field Eq. (11) is then used as an initial condition to the nonlinear Schrödinger equation Eq. (2). However, in contrast to the circular waveguide, the orientation of the dipoles D 1 and D 2 is important in nonlinear regime. If the two orthogonal dipoles D 1 and D 2 are oriented along the diagonals of the waveguide cross-section (see first subplot of Fig. 3), for a given amplitude A rotation occurs only if the modulation parameter B exceeds a critical value B = B cr . Moreover, the rotation is no longer constant as in the case of cylindrical waveguide but fluctuates and hence in what follows we use its average value (termed "average frequency")ω = 1/L L 0 ω(z) dz with L being the propagation distance corresponding to one full 2π rotation. This threshold value B cr decreases if we choose different initial dipole orientations 1 , and appears to be zero (within our numerical accuracy) for parallel orientation with respect to the waveguide boundaries. Let us focus on the case where the initial dipole-like field structure is oriented along the diagonal of the waveguide cross-section, and from now on the notation D 1 and D 2 stands for this orientation. Figure 3 shows an exemplary propagation of the resulting azimuthon-like solution with B > B cr , and an average rotation frequency ofω ≈ 0.0262. The top row depicts intensity and phase distribution at different propagation distance while the bottom plots illustrate the full 3D evolution of the "soliton". The rotation (counter-rotating w.r.t. phase) accompanied by beam deformation is clearly visible. For an amplitude ratio B less than the critical value (B < B cr ) the azimuthon no longer rotates but swings back and forth and hence its average rotation frequency is zero. The propagation of this wobbling dipole with an amplitude ratio B smaller than the critical value is displayed in Fig. 4. The right panels show the maximum angle which the dipole attains during propagation. The 3D surface plot in the bottom row clearly illustrates the swinging movement of the "soliton" upon propagation. Due to the azimuthon profile deformation, it is not possible to find an analytical expression ofω as a function of A and B as in the case of circular waveguide (i.e. Eq. (9)). Therefore, we need to resort to numerical simulations. In Fig. 5 we show the numerically determined relation betweenω and the azimuthon parameters A (left panel) and B (right panel). The threshold-like behavior is evident in the right graph of this figure: the region B < B cr , in which the soliton wobbles, is depicted by a green line.
Let us have a closer look at the boundary between domains of azimuthon rotation (B > B cr ) and wobbling (B < B cr ). We find numerically that the critical value B cr , which separates those two domains, depends very weakly on the azimuthon amplitude A: we are in weakly nonlinear limit, and may use linear waveguide modes to elucidate nonlinear propagation properties of the azimuthons. The appearance of a threshold value B cr , above which the azimuthon rotates, can be explained by considering the Hamiltonian associated with propagation equation Eq. (2), which is a conserved quantity upon propagation. Depending on their orientation, stationary dipoles (B = 0) have a slightly different value of H in nonlinear regime. If we take our diagonal (w.r.t. waveguide cross-section) dipoles D 1 and D 2 from above, we can construct a parallel dipole D p in the following way: The intensity distribution of this dipole D p is very close to that shown in the snapshot of the rotating azimuthon shown in the right panel in Fig. 3. One can show using the definition of the Hamiltonian Eq. (12) that H (D p ) > H (D 1 ) = H (D 2 ). In fact, it turns out that those two dipole orientations (parallel and diagonal) correspond to the extremal values of the Hamiltonian. Now, if the azimuthon rotates it has to pass through all possible orientations, i.e., its value of H has to be larger than H (D p ). We can use this reasoning to estimate the value of B cr : The left panel of Fig. 6 depicts the dependence of the Hamiltonian of superposed diagonal dipoles as a function of B (red line) for constant power P. The blue line represents the value of the Hamiltonian of the parallel dipole ∼ D p at the same power level. As the red curve monotonically increases with modulation parameter B, the intersection of the two curves gives an estimate of the critical value of the azimuthon modulation B cr , which is shown in the right panel in Fig. 5 by a black square. In other words, in order to rotate from the diagonal to vertical position the azimuthon must overcome some kind of energy barrier represented by the difference of the values of the Hamiltonian in those two basic states. As Fig. 6 shows this can be achieved by increasing the value of B in the diagonal state to B 1 = 0.342, which is very close to the value found numerically (see right panel in Fig. 5). Similar phenomena have been identified, for instance, in discrete nonlinear systems where the discreteness introduces the so called Peierls-Nabarro potential which has to be overcome by a soliton to become mobile [31]. It should be emphasized that the crucial difference between azimuthon dynamics in circular and square waveguide originates from the fact that angular momentum of the beam is conserved in the circular waveguide, whereas it changes considerably in the square waveguide. In the right panel of Fig. 6 we display the angular momentum of the beam in the square waveguide for four different values of modulation parameter B. If B is above the critical value B cr , the angular momentum does not change its sign as the blue curve shows. If B < B cr , the angular momentum changes periodically its sign in propagation (red curve), which indicates the swinging motion of the soliton. For values B close to the critical value the variation of the angular momentum  Fig. 3, blue curve), and the twisting ones (corresponding to Fig. 4, red curve). The other curves depict the cases of B slightly above (black) and bellow (green) B cr . All curves are obtained for A = 0.4. are largest as shown by the green and black curves. At the same time, the appearance of broad plateaus indicates slow propagation dynamics.
It is worth mentioning that we also investigated other waveguide symmetries, such as hexagonal waveguides, and found the behavior of the dipole azimuthons to be qualitatively similar to that of the square waveguide, i.e., they either rotate or wobble, depending on modulation depth and initial orientation. However, the difference of the values of the Hamiltonian for the two extremal dipoles is now smaller, resulting in a lower threshold value B cr . We believe that this threshold behavior in the propagation dynamics of azimuthons is generic for any non-circularsymmetric waveguide structure which supports degenerated dipole modes.

Rotating higher order localized modes
In analogy to dipole-azimuthons discussed so far, it is possible to consider dynamics of higher order azimuthons constructed by using pairs of degenerated higher order waveguide modes. In particular, in a circular waveguide degenerate higher order modes have the form of quadrupoles, hexapoles, octapoles, decapoles, etc. and other modes with even number of lobes (optical necklaces). In a square waveguide, one can identify pairs of degenerate modes in a form of hexapoles, octapoles, decapoles, dodecapoles, etc. (optical matrices). Interestingly, the quadrupoles found in square waveguides are not degenerate and hence the propagation dynamics of their corresponding superposed states in weakly nonlinear regime is dominated by linear mode beating. Fig. 7 displays an example of hexapole azimuthon in a circular waveguide. The left panel represents the initial intensity and phase distribution of the azimuthon. The 3D surface plot in the right panel depicts stable and smooth rotation of the hexapole azimuthon. As far as the higher order azimuthons of the square waveguide are concerned, they share the dynamical properties of dipole azimuthons discussed before. As an example, Fig. 8 shows the rotating (top row) and swinging (middle row) hexapoles. Due to the beam deformation, the hexapole with B = 0.5 transforms into a (2 × 3) matrix while rotating. The left two panels in the third row of Fig. 8 display the rotating dodecapole azimuthon and its deformation into a (3 × 4) solitonic matrix; the right panels show rotating icosapole azimuthon and its deformation into a (4 × 5) solitonic matrix.

Conclusion
In conclusion, we demonstrated numerically stable propagation of azimuthons, i.e. localized rotating nonlinear waves in weakly nonlinear waveguides. Depending on the waveguide profile, different nonlinearity induced propagation dynamics can be observed. We showed that azimuthons in circular waveguides rotate continuously. The analytically predicted dependence of rotation frequency ω as a function of soliton parameters was found to be in excellent agreement with numerical simulations. Further on, we discussed propagation of azimuthon-like structures in square waveguides. We showed that their dynamics critically depends on the initial conditions. In particular, we found a threshold-like behavior in the propagation dynamics, separating rotating and wobbling solutions. We showed that this effect is associated with different values of the Hamiltonian for different azimuthon orientations. As our analysis relies on physical parameters of actual multi-mode waveguides, our findings may open a relatively easy route to experimental observations of stable azimuthons.

Acknowledgement
Numerical simulations were partly performed on the SGI XE Cluster and the Sun Constellation VAYU Cluster of the Australian Partnership for Advanced Computing (APAC). This research was supported by the Australian Research Council.