---
abstract: '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.'
address: |
    $^1$Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Stra[ß]{}e 38, 01187 Dresden, Germany\
    $^2$State Key Laboratory of Transient Optics and Photonics, Xi’an Institute of Optics and Precision Mechanics of Chinese Academy of Sciences, 710119 Xi’an, China\
    $^3$Institute of Condensed Matter Theory and Solid State Optics, Friedrich Schiller University, Max-Wien-Platz 1, 07743 Jena, Germany\
    $^4$Laser Physics Center, Research School of Physics and Engineering, Australian National University, Canberra, ACT 0200, Australia
author:
- 'Yiqi Zhang$^{1,2}$, Stefan Skupin$^{1,3}$, Fabian Maucher$^1$, Arpa Galestian Pour$^3$, Keqing Lu$^{2}$, Wieslaw Królikowski$^4$'
title: Azimuthons in weakly nonlinear waveguides of different symmetries
---

[10]{}

G. I. Stegeman and M. Segev, “Optical spatial solitons and their interactions: Universality and diversity,” Science **286**, 1518–1523 (1999).

A. S. Desyatnikov, A. A. Sukhorukov, and Y. S. Kivshar, “Azimuthons: Spatially modulated vortex solitons,” Phys. Rev. Lett. **95**, 203904 (2005).

P. Coullet, L. Gil, and F. Rocca, “Optical vortices,” Opt. Commun. **73**, 403–408 (1989).

V. M. Lashkin, “Two-dimensional multisolitons and azimuthons in [B]{}ose-[E]{}instein condensates,” Phys. Rev. A **77**, 025602 (2008).

D. Buccoliero, A. S. Desyatnikov, W. Kr[ó]{}likowski, and Y. S. Kivshar, “Spiraling multivortex solitons in nonlocal nonlinear media,” Opt. Lett. **33**, 198–200 (2008).

S. Lopez-Aguayo, A. S. Desyatnikov, and Y. S. Kivshar, “Azimuthons in nonlocal nonlinear media,” Opt. Express **14**, 7903–7908 (2006).

D. Buccoliero, A. S. Desyatnikov, W. Kr[ó]{}likowski, and Y. S. Kivshar, “Laguerre and [H]{}ermite soliton clusters in nonlocal nonlinear media,” Phys. Rev. Lett. **98**, 053901 (2007).

A. Minovich, D. N. Neshev, A. S. Desyatnikov, W. Kr[ó]{}likowski, and Y. S. Kivshar, “Observation of optical azimuthons,” Opt. Express **17**, 23610–23616 (2009).

V. I. Kruglov, Y. A. Logvin, and V. M. Volkov, “The theory of spitral laser beams in nonlinear media,” J. Mod. Opt. **39**, 2277–2291 (1992).

D. V. Skryabin and W. J. Firth, “Dynamics of self-trapped beams with phase dislocation in saturable [K]{}err and quadratic nonlinear media,” Phys. Rev. E **58**, 3916–3930 (1998).

D. Suter and T. Blasberg, “Stabilization of transverse solitary waves by a nonlocal response of the nonlinear medium,” Phys. Rev. A **48**, 4583–4587 (1993).

O. Bang, W. Kr[ó]{}likowski, J. Wyller, and J. J. Rasmussen, “Collapse arrest and soliton stabilization in nonlocal nonlinear media,” Phys. Rev. E **66**, 046619 (2002).

W. Kr[ó]{}likowski, O. Bang, N. I. Nikolov, D. Neshev, J. Wyller, J. J. Rasmussen, and D. Edmundson, “Modulational instability, solitons and beam propagation in spatially nonlocal nonlinear media,” J. Opt. B: Quantum Semiclass. Opt. **6**, S288–S294 (2004).

D. Briedis, D. Petersen, D. Edmundson, W. Kr[ó]{}likowski, and O. Bang, “Ring vortex solitons in nonlocal nonlinear media,” Opt. Express **13**, 435–443 (2005).

S. Lopez-Aguayo, A. S. Desyatnikov, Y. S. Kivshar, S. Skupin, W. Kr[ó]{}likowski, and O. Bang, “Stable rotating dipole solitons in nonlocal optical media,” Opt. Lett. **31**, 1100–1102 (2006).

S. Skupin, M. Grech, and W. Kr[ó]{}likowski, “Rotating soliton solutions in nonlocal nonlinear media,” Opt. Express **16**, 9118–9131 (2008).

F. Maucher, D. Buccoliero, S. Skupin, M. Grech, A. S. Desyatnikov, and W. Krolikowski, “[Tracking azimuthons in nonlocal nonlinear media]{},” [Opt. Quant. Electron.]{} **[41]{}**, [337–348]{} ([2009]{}).

M. Peccianti, K. A. Brzdkiewicz, and G. Assanto, “Nonlocal spatial soliton interactions in nematic liquid crystals,” Opt. Lett. **27**, 1460–1462 (2002).

C. Conti, M. Peccianti, and G. Assanto, “Route to nonlocality and observation of accessible solitons,” Phys. Rev. Lett. **91**, 073901 (2003).

C. Conti, M. Peccianti, and G. Assanto, “Observation of optical spatial solitons in a highly nonlocal medium,” Phys. Rev. Lett. **92**, 113902 (2004).

R. Nath, P. Pedri, and L. Santos, “Soliton-soliton scattering in dipolar [B]{}ose-[E]{}instein condensates,” Phys. Rev. A **76**, 013606 (2007).

V. M. Lashkin, E. A. Ostrovskaya, A. S. Desyatnikov, and Y. S. Kivshar, “Vector azimuthons in two-component [B]{}ose-[E]{}instein condensates,” Phys. Rev. A **80**, 013615 (2009).

A. G. Litvak, V. Mironov, G. Fraiman, and A. Yunakovskii, “Thermal self-effect of wave beams in plasma with a nonlocal nonlinearity,” Sov. J. Plasma Phys. **1**, 31–37 (1975).

C. Rotschild, O. Cohen, O. Manela, M. Segev, and T. Carmon, “Solitons in nonlinear media with an infinite range of nonlocality: First observation of coherent elliptic solitons and of vortex-ring solitons,” Phys. Rev. Lett. **95**, 213904 (2005).

S. Skupin, U. Peschel, L. Berg[é]{}, and F. Lederer, “Stability of weakly nonlinear localized states in attractive potentials,” Phys. Rev. E **70**, 016614 (2004).

R. W. Boyd, *Nonlinear Optics* (Academic Press, Amsterdam, 2008), 3rd ed.

N. Kuzuu, K. Yoshida, H. Yoshida, T. Kamimura, and N. Kamisugi, “Laser-induced bulk damage in various types of vitreous silica at 1064, 532, 355, and 266 nm: Evidence of different damage mechanisms between 266-nm and longer wavelengths,” Appl. Opt. **38**, 2510–2515 (1999).

G. P. Agrawal, *Nonlinear Fiber Optics* (Academic Press, Singapore, 2009), 4th ed.

D. V. Skryabin, J. M. McSloy, and W. J. Firth, “Stability of spiralling solitary waves in [H]{}amiltonian systems,” Phys. Rev. E **66**, 055602 (2002).

F. Maucher, S. Skupin, M. Shen, and W. Kr[ó]{}likowski, “Rotating three-dimensional solitons in [B]{}ose-[E]{}instein condensates with gravitylike attractive nonlocal interaction,” Phys. Rev. A **81**, 063617 (2010).

Y. S. Kivshar and D. K. Campbell, “[P]{}eierls-[N]{}abarro potential barrier for highly localized nonlinear modes,” Phys. Rev. E **48**, 3077–3081 (1993).

Introduction
============

Spatial solitons are nonlinear localized states that keep their form during propagation due to the balance between diffraction and self-induced nonlinear potential [@Stegeman_science_1999]. 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 [@desyatnikov_prl_2005]. They can be considered as azimuthally perturbed optical vortices, i.e. beams with singular phase structure [@coullet_oc_1989; @desyatnikov_prl_2005; @lashkin_pra_2008]. Theoretical studies demonstrated both, stable and unstable propagation of azimuthons [@buccoliero_ol_2008; @lopez-aguayo_oe_2006; @buccoliero_prl_2007], and the first experimental observation of optical azimuthons was recently achieved in rubidium vapors [@minovich_oe_2009].

It appears that higher order solitonic structures and optical vortices are generally unstable in typical nonlinear media with local (Kerr-like) response [@kruglov_jmo_1992; @skryabin_pre_1998]. On the other hand, it has been shown that a spatially nonlocal nonlinear response provides stabilization of various complex solitonic structures including vortices [@suter_pra_1993; @bang_pre_2002; @krolikowski_job_2004; @briedis_oe_2005]. Consequently, azimuthons and their dynamics have been studied almost exclusively in the context of spatially nonlocal nonlinear media [@lopez-aguayo_ol_2006; @buccoliero_prl_2007; @stefan_oe_2008; @Fabian:oqe:09]. In spite of the fact that there are various physical settings exhibiting nonlocality such as nematic liquid crystals [@peccianti_ol_2002; @conti_prl_2003; @conti_prl_2004], Bose-Einstein condensates [@lashkin_pra_2008; @nath_pra_2007; @lashkin_pra_2009], plasmas [@litvak_plasmas_1975], thermo-optical materials [@rotschild_prl_2005] 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 [@stefan_pre_2004], we can expect that weakly nonlinear azimuthons are stable in multi-mode waveguides.

The paper is organized as follows: In Sec. \[modelling\], 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. \[dipoles\], respectively. In Sec. \[higherorder\], higher order azimuthons are investigated, and in Sec. \[conclusion\] we conclude.

Mathematical modeling {#modelling}
=====================

We consider the evolution of a continuous wave (CW) optical beam with amplitude $\mathcal{E} (\xi, \eta, \zeta)$, where $(\xi,\eta)$ and $\zeta$ 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 [@boyd_book_2008]: $$i \frac{\partial}{\partial \zeta} \mathcal{E} +
\frac{1}{2k_0} \left( \frac{\partial ^2}{\partial \xi^2} + \frac{\partial ^2}{\partial \eta^2} \right) \mathcal{E} +
k_0 \frac{n_2}{n_b} |\mathcal{E}|^2 \mathcal{E} + k_0 \frac{n(\xi,\eta)-n_b}{n_b} \mathcal{E}=0,
\label{physicaleq}$$ where $k_0=2 \pi n_b / \lambda_0$ refers to the carrier central wave number in the medium; $n_2$ is the Kerr nonlinear coefficient, $n(\xi,\eta)$ the linear refractive index distribution, $n_b$ the background index, and $\lim_{\xi,\eta \rightarrow \infty}n(\xi,\eta)=n_b$. We consider a weakly-guiding waveguide, so both the linear and nonlinear induced index change $|n-n_b|$ and $n_2|\mathcal{E}|^2$ are small compared to the mean index $n_b$, and at the same time $n_2|\mathcal{E}|^2 \ll |n-n_b|$, to guarantee a weak nonlinearity.

From the mathematical point of view, Eq. (\[physicaleq\]) is a (2+1)-dimensional nonlinear Schr[ö]{}dinger (NLS) equation with linear potential representing the waveguide profile. For technical convenience, we rescale Eq. (\[physicaleq\]) to dimensionless quantities with $x=\xi/r_0$, $y=\eta/r_0$, $z=\zeta/(2 k_0 r_0^2)$, and $\sigma=\mathrm{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 $\psi = k_0 r_0 \sqrt{2|n_2|/n_b} \mathcal{E}$ reads $$i\frac{\partial}{\partial z} \psi + \left( \frac{\partial ^2}{\partial x^2}
 + \frac{\partial ^2}{\partial y^2} \right) \psi +
\sigma |\psi|^2 \psi + V \psi=0,
\label{partial}$$ with an “attractive” bounded potential $V=2k_0^2 r_0^2 (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(x,y)=
\begin{cases}
V_0 & \mathrm{~where~}
\begin{cases}
\sqrt{x^2+y^2} \leq 1 & \mathrm{~for~circular~waveguide}
\\
|x| \leq 1 ~\&~ |y| \leq 1 & \mathrm{~for~square~waveguide}
\end{cases}
\\
0 & \mathrm{elsewhere}
\end{cases}.$$ $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| \leq 9 \times 10^{-3}$, $n_2 = 3\times 10^{-16} ~\mathrm{cm}^2/\mathrm{W}$ at a vacuum wavelength of $\lambda_0=790$ nm, and those values do not vary much when we choose any $\lambda_0$ in the range from visible to near-infrared. The laser-induced-damage-threshold (LIDT) for synthetic fused silica is about $10 ~\mathrm{J}/\mathrm{cm}^2$ for nanosecond pulses (according to Eq. (1) and Table 2 in Ref. [@kuzuu_ao_1999]) corresponding to $\sim 10^{10} ~\mathrm{W}/\mathrm{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 \approx 25.0 ~\mu \mathrm{ m}$ through $r_0=\sqrt{\left|n_bV/\left[2k_0^2(n-n_b)\right]\right|}$, which nicely corresponds to the radius of a standard circular multi-mode fiber [@agrawal_book_2009]. In addition, in order to guarantee intensities smaller than the LIDT the condition $|\psi|^2<1/3$ should be fulfilled.

Note that in Eq. (\[partial\]) $\sigma$ is the sign of nonlinearity: $\sigma =1$ ($\sigma=-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.

Rotating localized dipoles {#dipoles}
==========================

Dipole azimuthons in circular waveguides
----------------------------------------

Azimuthons are a straightforward generalization of the usual ansatz for stationary solutions [@desyatnikov_prl_2005]. They represent spatially rotating structures and hence involve an additional parameter, the rotation frequency $\omega$ (see also [@skryabin_pre_2002]), so we seek approximate solutions of the form $$\psi (r,\phi,z) = U (r, \phi - \omega z) \exp (i \kappa z),
\label{ansatz}$$ where $r=\sqrt{x^2+y^2}$ and $\phi$ the azimuthal angle in the transverse plane $(x,y)$, $U$ is the stationary profile, $\omega$ the rotation frequency, and $\kappa$ the propagation constant. For $\omega=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 [@lopez-aguayo_ol_2006]. 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 $\pi \int r|F(r)|^2 \mathrm{d}r =1$. Then, using separation of variables, we consider the simplest so-called “rotating dipole” azimuthon with ansatz [@stefan_oe_2008] $$U(r,\phi-\omega z)=AF(r)\left[\cos(\phi-\omega z)+i B \sin(\phi-\omega z) \right],
\label{complexamplitude}$$ where $A$ is an amplitude factor, and $1-B$ the azimuthal modulation depth of the resulting ring-like 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. (\[ansatz\]) into Eq. (\[partial\]), multiplying by $U^*$ and $\partial U^* / \partial \phi$ respectively, and integrating over the transverse coordinates we end up with [@stefan_oe_2008]

$$-\kappa P + \omega L_z + I + N = 0,$$

$$-\kappa L_z + \omega P' + I' + N' = 0.$$

This system relates the propagation constant $\kappa$ and the rotation frequency $\omega$ of the azimuthons to integrals over their stationary amplitude profiles:

$$P = \iint r|U(r)|^2 ~\mathrm{d}r~\mathrm{d}\phi,$$

$$L_z = -i\iint r\frac{\partial U(r)}{\partial \phi} U^*(r) ~\mathrm{d}r~\mathrm{d}\phi,$$

$$I = \iint rU^*(r) \Delta_{\perp} U(r) ~\mathrm{d}r~\mathrm{d}\phi,$$

$$N = \iint r\left[\sigma|U(r)|^2+V \right] |U(r)|^2 ~\mathrm{d}r~\mathrm{d}\phi,$$

$$P' = \iint r\left| \frac{\partial U(r)}{\partial \phi} \right|^2 ~\mathrm{d}r~\mathrm{d}\phi,$$

$$I' = i\iint r\frac{\partial U^*(r)}{\partial \phi} \Delta_{\perp} U(r) ~\mathrm{d}r~\mathrm{d}\phi,$$

$$N' = i\iint r\left[\sigma|U(r)|^2+V \right] \frac{\partial U^*(r)}{\partial \phi} U(r) ~\mathrm{d}r~\mathrm{d}\phi.$$

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 $$\omega=\frac{P(I'+N')-L_z(I+N)}{L_z^2-PP'}.
\label{angularfrequency}$$

After inserting Eq. (\[complexamplitude\]) into Eq. (\[angularfrequency\]), it turns out that only the nonlinear term contributes to the rotation frequency $\omega$ (see also Eq. (10) in Ref. [@fabian_pra_2010]). 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$: $$F(r) = C \times
\begin{cases}
J_1 ( \sqrt{V_0-\kappa} r) & \mathrm{for}~ 0\leq r\leq1
 \\
\dfrac{J_1(\sqrt{V_0-\kappa})}{K_1(\sqrt{\kappa})} \cdot K_1 (\sqrt{\kappa}r) & \mathrm{for}~ r>1
\end{cases},$$ where $C$ is a normalization factor such that $\pi\int|F(r)|^2rdr=1$. Thus, the analytically predicted rotation frequency is $$\omega=\frac{\pi \sigma}{2} \int r |F(r)|^4 \mathrm{d}r \cdot A^2 B.
\label{analyticalangular}$$ From the above equation, one can see that sense of rotation is opposite for focusing and defocusing nonlinearities. Moreover, as $\omega = 2 k_0 r_0^2 \cdot 2 \pi/ \ell$, one can find from Eq. (\[analyticalangular\]) an expression for the physical distance $\ell$ over which the azimuthon in a circular waveguide performs one full rotation: $$\ell= \frac{\displaystyle 8 k_0 r_0^2}{\displaystyle A^2 B \int r |F(r)|^4 \mathrm{d}r}.
\label{physicaldistance}$$ For example, assuming propagation in a silica fiber with the parameters discussed earlier (see Sec. \[modelling\]) and taking $A=0.4$, $B=0.5$, we obtain $\ell \approx 2.4 ~\mathrm{m}$.

Figure \[circularrelations\] 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 $\omega \approx 0.037$ is illustrated in Fig. \[circularevolution\]. 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 $\pi/2$.

![Azimuthon rotation frequency $\omega$ versus amplitude factor $A$ (left panel) and amplitude ratio $B$ (right panel). Blue solid lines show analytical predictions from Eq. (\[analyticalangular\]). Red dots denote results obtained from numerical simulations of Eq. (\[partial\]). The values of $A$ and $B$ are shown next to the lines.[]{data-label="circularrelations"}](relation_ca "fig:"){width="40.00000%"} ![Azimuthon rotation frequency $\omega$ versus amplitude factor $A$ (left panel) and amplitude ratio $B$ (right panel). Blue solid lines show analytical predictions from Eq. (\[analyticalangular\]). Red dots denote results obtained from numerical simulations of Eq. (\[partial\]). The values of $A$ and $B$ are shown next to the lines.[]{data-label="circularrelations"}](relation_cb "fig:"){width="40.00000%"}

![The propagation of a dipole azimuthon with $A=0.4,~B=0.5$ in the circular waveguide. The iso-surface plot on the top clearly displays the spiraling of the azimuthon during propagation. Figures in the row below depict the azimuthon’s transverse intensity distribution at input, and after rotating $\pi/2$, respectively. The white line indicates the waveguide boundaries.[]{data-label="circularevolution"}](cd_prop "fig:"){width="55.00000%"}\
![The propagation of a dipole azimuthon with $A=0.4,~B=0.5$ in the circular waveguide. The iso-surface plot on the top clearly displays the spiraling of the azimuthon during propagation. Figures in the row below depict the azimuthon’s transverse intensity distribution at input, and after rotating $\pi/2$, respectively. The white line indicates the waveguide boundaries.[]{data-label="circularevolution"}](cdipole0 "fig:"){width="35.00000%"} ![The propagation of a dipole azimuthon with $A=0.4,~B=0.5$ in the circular waveguide. The iso-surface plot on the top clearly displays the spiraling of the azimuthon during propagation. Figures in the row below depict the azimuthon’s transverse intensity distribution at input, and after rotating $\pi/2$, respectively. The white line indicates the waveguide boundaries.[]{data-label="circularevolution"}](cdipole85 "fig:"){width="35.00000%"}

Azimuthon-like dipoles in square waveguides
-------------------------------------------

Because of the lack of circular symmetry, azimuthons in the strict sense of Eq. (\[complexamplitude\]) (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 $\iint |D_{1,2}|^2~\mathrm{d}x~\mathrm{d}y=1$, and superpose them as before to form the azimuthon-like object $$U(x,y,z=0)=A\left(D_1+iBD_2\right).
\label{squareansatz}$$ The field Eq. (\[squareansatz\]) is then used as an initial condition to the nonlinear Schr[ö]{}dinger equation Eq. (\[partial\]). 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. \[square\_r\_evolution\]), for a given amplitude $A$ rotation occurs only if the modulation parameter $B$ exceeds a critical value $B=B_{\rm 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”) $\bar \omega=1/L\int_0^L\omega(z)~\mathrm{d}z$ with $L$ being the propagation distance corresponding to one full $2\pi$ rotation. This threshold value $B_{\rm 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.

![The propagation of the azimuthon-like dipole, that rotates in the square waveguide. Top row: intensity (left) and phase (right) distribution corresponding to the input and soliton rotation by $\pi/4$, respectively. The white line indicates the waveguide boundaries. The iso-surface plot at the bottom displays the rotation and deformation during propagation. The initial amplitude and modulation parameters are $A=0.4,~B=0.5$.[]{data-label="square_r_evolution"}](sdr0 "fig:"){width="35.00000%"} ![The propagation of the azimuthon-like dipole, that rotates in the square waveguide. Top row: intensity (left) and phase (right) distribution corresponding to the input and soliton rotation by $\pi/4$, respectively. The white line indicates the waveguide boundaries. The iso-surface plot at the bottom displays the rotation and deformation during propagation. The initial amplitude and modulation parameters are $A=0.4,~B=0.5$.[]{data-label="square_r_evolution"}](sdr30 "fig:"){width="35.00000%"}\
![The propagation of the azimuthon-like dipole, that rotates in the square waveguide. Top row: intensity (left) and phase (right) distribution corresponding to the input and soliton rotation by $\pi/4$, respectively. The white line indicates the waveguide boundaries. The iso-surface plot at the bottom displays the rotation and deformation during propagation. The initial amplitude and modulation parameters are $A=0.4,~B=0.5$.[]{data-label="square_r_evolution"}](sdr_prop "fig:"){width="55.00000%"}

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 \[square\_r\_evolution\] shows an exemplary propagation of the resulting azimuthon-like solution with $B>B_{\rm cr}$, and an average rotation frequency of $\bar \omega \approx 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_{\rm 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. \[square\_t\_evolution\]. 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.

![The propagation of a dipole azimuthon that wobbles in the square waveguide. The first row depicts the dipole at input and maximum displacement, respectively. The white line indicates the waveguide boundaries. The iso-surface plot below displays the twist and deformation during propagation. The chosen amplitude factor and ratio are $A=0.4,~B=0.2$.[]{data-label="square_t_evolution"}](sdt0 "fig:"){width="35.00000%"} ![The propagation of a dipole azimuthon that wobbles in the square waveguide. The first row depicts the dipole at input and maximum displacement, respectively. The white line indicates the waveguide boundaries. The iso-surface plot below displays the twist and deformation during propagation. The chosen amplitude factor and ratio are $A=0.4,~B=0.2$.[]{data-label="square_t_evolution"}](sdt40 "fig:"){width="35.00000%"}\
![The propagation of a dipole azimuthon that wobbles in the square waveguide. The first row depicts the dipole at input and maximum displacement, respectively. The white line indicates the waveguide boundaries. The iso-surface plot below displays the twist and deformation during propagation. The chosen amplitude factor and ratio are $A=0.4,~B=0.2$.[]{data-label="square_t_evolution"}](sdt_prop "fig:"){width="55.00000%"}

Due to the azimuthon profile deformation, it is not possible to find an analytical expression of $\bar \omega$ as a function of $A$ and $B$ as in the case of circular waveguide (i.e. Eq. (\[analyticalangular\])). Therefore, we need to resort to numerical simulations. In Fig. \[squarerelations\] we show the numerically determined relation between $\bar \omega$ 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_{\rm cr}$, in which the soliton wobbles, is depicted by a green line.

![Averaged rotation frequency $\bar \omega$ versus amplitude factor $A$ (left panel) and amplitude ratio $B$ (right panel) in a square waveguide. Red dots are numerical results obtained from Eq. (\[partial\]), blue curves ($B<0.4$) and blue lines ($B>0.4$) are the fitting results to the numerical simulations. The green line represents the superposed dipoles twist during propagation. The values of $A$ and $B$ are shown next to the curves. The black square represents the analytical estimate for $B_{\rm cr}$.[]{data-label="squarerelations"}](relation_sa "fig:"){width="40.00000%"} ![Averaged rotation frequency $\bar \omega$ versus amplitude factor $A$ (left panel) and amplitude ratio $B$ (right panel) in a square waveguide. Red dots are numerical results obtained from Eq. (\[partial\]), blue curves ($B<0.4$) and blue lines ($B>0.4$) are the fitting results to the numerical simulations. The green line represents the superposed dipoles twist during propagation. The values of $A$ and $B$ are shown next to the curves. The black square represents the analytical estimate for $B_{\rm cr}$.[]{data-label="squarerelations"}](relation_sb "fig:"){width="40.00000%"}

Let us have a closer look at the boundary between domains of azimuthon rotation ($B>B_{\rm cr}$) and wobbling ($B<B_{\rm cr}$). We find numerically that the critical value $B_{\rm 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_{\rm cr}$, above which the azimuthon rotates, can be explained by considering the Hamiltonian associated with propagation equation Eq. (\[partial\]), $$\mathcal{H}=\iint \left(|\nabla \psi|^2-\frac{1}{2}|\psi|^4-V|\psi|^2 \right)~\mathrm{d}x ~\mathrm{d}y,
\label{hamiltonian}$$ which is a conserved quantity upon propagation. Depending on their orientation, stationary dipoles ($B=0$) have a slightly different value of $\mathcal{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: $$D_p=\frac{D_1+D_2}{\displaystyle \sqrt{\iint |D_1+D_2|^2 ~\mathrm{d}x ~\mathrm{d}y}}.$$ 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. \[square\_r\_evolution\]. One can show using the definition of the Hamiltonian Eq. (\[hamiltonian\]) that $\mathcal{H}(D_p)>\mathcal{H}(D_1)=\mathcal{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 $\mathcal{H}$ has to be larger than $\mathcal{H}(D_p)$. We can use this reasoning to estimate the value of $B_{\rm cr}$: The left panel of Fig. \[angularmomentumcurve\] 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 $\sim 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_{\rm cr}$, which is shown in the right panel in Fig. \[squarerelations\] 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. \[angularmomentumcurve\] 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. \[squarerelations\]). 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 [@Kivshar_Campbell].

![ Left panel: Hamiltonian of $A(D_1+iBD_2)$ (red curve) as a function of $B$ and Hamiltonian of $D_3$ (blue line) (for the same power $P$). The intersection of the two curves defines $B_{\rm cr}$ (depicted by the black square). Right panel: The angular momenta of the rotating superposed dipoles (corresponding to Fig. \[square\_r\_evolution\], blue curve), and the twisting ones (corresponding to Fig. \[square\_t\_evolution\], red curve). The other curves depict the cases of $B$ slightly above (black) and bellow (green) $B_{\rm cr}$. All curves are obtained for $A=0.4$. []{data-label="angularmomentumcurve"}](hamil "fig:"){width="45.00000%"} ![ Left panel: Hamiltonian of $A(D_1+iBD_2)$ (red curve) as a function of $B$ and Hamiltonian of $D_3$ (blue line) (for the same power $P$). The intersection of the two curves defines $B_{\rm cr}$ (depicted by the black square). Right panel: The angular momenta of the rotating superposed dipoles (corresponding to Fig. \[square\_r\_evolution\], blue curve), and the twisting ones (corresponding to Fig. \[square\_t\_evolution\], red curve). The other curves depict the cases of $B$ slightly above (black) and bellow (green) $B_{\rm cr}$. All curves are obtained for $A=0.4$. []{data-label="angularmomentumcurve"}](sangular "fig:"){width="45.00000%"}

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. \[angularmomentumcurve\] 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_{\rm cr}$, the angular momentum does not change its sign as the blue curve shows. If $B<B_{\rm 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 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_{\rm cr}$. We believe that this threshold behavior in the propagation dynamics of azimuthons is generic for any non-circular-symmetric waveguide structure which supports degenerated dipole modes.

Rotating higher order localized modes {#higherorder}
=====================================

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.

![The propagation of a hexapole azimuthon with $A=0.4,~B=0.5$ in the circular waveguide. Left panel shows the original hexapole azimuthon and the corresponding phase; right panel shows the iso-surface plot of the propagation.[]{data-label="cquadr_evolution"}](chexa0 "fig:"){width="45.00000%"} ![The propagation of a hexapole azimuthon with $A=0.4,~B=0.5$ in the circular waveguide. Left panel shows the original hexapole azimuthon and the corresponding phase; right panel shows the iso-surface plot of the propagation.[]{data-label="cquadr_evolution"}](chexa_prop "fig:"){width="50.00000%"}

Fig. \[cquadr\_evolution\] 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. \[hexa\_evolution\] 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 \times 3$) matrix while rotating. The left two panels in the third row of Fig. \[hexa\_evolution\] display the rotating dodecapole azimuthon and its deformation into a ($3 \times 4$) solitonic matrix; the right panels show rotating icosapole azimuthon and its deformation into a ($4 \times 5$) solitonic matrix.

![The propagation of superposed hexapoles in the square waveguide. The first row shows a rotating azimuthon with $A=0.4,~B=0.5>B_{\rm cr}$ (Media 1), and the second row shows a twisting azimuthon with $A=0.4,~B=0.2<B_{\rm cr}$ (Media 2). The superposed dodecapoles (left two panels), and superposed icosapoles (right two panels) as well as their deformations are shown in the third row; in both cases we choose $A=0.4$ and $B=0.5$.[]{data-label="hexa_evolution"}](shexa_r0 "fig:"){width="35.00000%"} ![The propagation of superposed hexapoles in the square waveguide. The first row shows a rotating azimuthon with $A=0.4,~B=0.5>B_{\rm cr}$ (Media 1), and the second row shows a twisting azimuthon with $A=0.4,~B=0.2<B_{\rm cr}$ (Media 2). The superposed dodecapoles (left two panels), and superposed icosapoles (right two panels) as well as their deformations are shown in the third row; in both cases we choose $A=0.4$ and $B=0.5$.[]{data-label="hexa_evolution"}](shexa_r30 "fig:"){width="35.00000%"}\
![The propagation of superposed hexapoles in the square waveguide. The first row shows a rotating azimuthon with $A=0.4,~B=0.5>B_{\rm cr}$ (Media 1), and the second row shows a twisting azimuthon with $A=0.4,~B=0.2<B_{\rm cr}$ (Media 2). The superposed dodecapoles (left two panels), and superposed icosapoles (right two panels) as well as their deformations are shown in the third row; in both cases we choose $A=0.4$ and $B=0.5$.[]{data-label="hexa_evolution"}](shexa_t0 "fig:"){width="35.00000%"} ![The propagation of superposed hexapoles in the square waveguide. The first row shows a rotating azimuthon with $A=0.4,~B=0.5>B_{\rm cr}$ (Media 1), and the second row shows a twisting azimuthon with $A=0.4,~B=0.2<B_{\rm cr}$ (Media 2). The superposed dodecapoles (left two panels), and superposed icosapoles (right two panels) as well as their deformations are shown in the third row; in both cases we choose $A=0.4$ and $B=0.5$.[]{data-label="hexa_evolution"}](shexa_t30 "fig:"){width="35.00000%"}\
![The propagation of superposed hexapoles in the square waveguide. The first row shows a rotating azimuthon with $A=0.4,~B=0.5>B_{\rm cr}$ (Media 1), and the second row shows a twisting azimuthon with $A=0.4,~B=0.2<B_{\rm cr}$ (Media 2). The superposed dodecapoles (left two panels), and superposed icosapoles (right two panels) as well as their deformations are shown in the third row; in both cases we choose $A=0.4$ and $B=0.5$.[]{data-label="hexa_evolution"}](sdodeca "fig:"){width="35.00000%"} ![The propagation of superposed hexapoles in the square waveguide. The first row shows a rotating azimuthon with $A=0.4,~B=0.5>B_{\rm cr}$ (Media 1), and the second row shows a twisting azimuthon with $A=0.4,~B=0.2<B_{\rm cr}$ (Media 2). The superposed dodecapoles (left two panels), and superposed icosapoles (right two panels) as well as their deformations are shown in the third row; in both cases we choose $A=0.4$ and $B=0.5$.[]{data-label="hexa_evolution"}](sicosa "fig:"){width="35.00000%"}

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 $\omega$ 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.

[^1]: In the linear system, any dipole orientation is possible, and can be constructed from orthogonal basis ($D_1,D_2$) by superposition.
