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 