Rabi oscillation, an interband oscillation, describes periodic motion between two states that belong to different energy levels, in the presence of an oscillatory driving field. In photonics, Rabi oscillations can be mimicked by applying a weak longitudinal periodic modulation to the refractive index. However, the Rabi oscillations of nonlinear states have yet to be introduced. We report the Rabi oscillations of azimuthons—spatially modulated vortex solitons—in weakly nonlinear waveguides with different symmetries. The period of the Rabi oscillations can be determined by applying the coupled mode theory, which largely depends on the modulation strength. Whether the Rabi oscillations between two states can be obtained or not is determined by the spatial symmetry of the azimuthons and the modulating potential. Our results not only deepen the understanding of the Rabi oscillation phenomena, but also provide a new avenue in the study of pattern formation and spatial field manipulation in nonlinear optical systems.

## 1.

## Introduction

Rabi oscillations were introduced in quantum mechanics,^{1} but by now are widely investigated in a variety of optical and photonic systems that include fibers,^{2}^{,}^{3} multimode waveguides,^{4}^{–}^{6} coupled waveguides,^{7} waveguide arrays,^{8}^{–}^{10} and two-dimensional modal structures.^{11}^{,}^{12} Recently, Rabi oscillations of topological edge states^{13}^{,}^{14} and modes in the fractional Schrödinger equation^{15} were also reported. Rabi oscillations are interband oscillations that require an AC field to be applied as an external periodic potential. In optics, the longitudinal periodic modulation of the refractive index plays the role of an AC field in temporal quantum systems, and Rabi oscillations are indicated by the resonant mode conversion. As far as we know, the investigation of optical Rabi oscillations thus far has been limited to the linear regime only, and the Rabi oscillation in nonlinear systems is still an open problem that needs to be explored.

The aim of this work is to investigate Rabi oscillations of azimuthons in weakly nonlinear waveguides that is accomplished by applying a weak longitudinally modulated periodic potential. Azimuthons are a special type of spatial soliton; they are azimuthally modulated vortex beams that exhibit steady rotation upon propagation.^{16} Generally, azimuthons, especially the ones with higher-order angular momentum structures, are unstable in media with local Kerr or saturable nonlinearities. To overcome the instability drawback, a nonlocal nonlinearity is introduced, and recently published reports demonstrate that the stable propagation of azimuthons can indeed be obtained.^{17}^{–}^{19} In addition, it was also reported that the spin-orbit-coupled Bose–Einstein condensates can support stable azimuthons as well.^{20} However, the treatment of nonlocal nonlinearity and spin-orbit-coupled Bose–Einstein condensates is challenging in both theoretical modeling and experimental demonstration. Nevertheless, it has been confirmed that weakly nonlinear waveguides^{21} represent an ideal platform for the investigation of stable azimuthons,^{22}^{,}^{23} even with higher-order modal structures.

Following this path of inquiry, we first investigate Rabi oscillations of azimuthons in a circular waveguide and then in a square waveguide. Since in this nonlinear three-dimensional wave propagation problem no analytical solutions are known, the mode of inquiry will necessarily be predominantly numerical with some theoretical background. In the circular waveguide, the azimuthons will exhibit Rabi oscillation while rotating during propagation. In the square waveguide, the behavior of the azimuthons is different in two aspects:^{15} (i) azimuthons will rotate only if the corresponding Hamiltonian (energy) is bigger than a certain threshold value and (ii) azimuthons will be deformed during propagation. Hence, in this work, we choose azimuthons with large enough energies to avoid wobbling motions in the square waveguide during propagation.

## 2.

## Results

## 2.1.

### Theoretical Analysis

The propagation of a light beam in a photonic waveguide can be described by the following Schrödinger-like paraxial wave equation:

## Eq. (1)

$$i\frac{\partial}{\partial Z}\mathrm{\Psi}+\frac{1}{2{k}_{0}}(\frac{{\partial}^{2}}{\partial {X}^{2}}+\frac{{\partial}^{2}}{\partial {Y}^{2}})\mathrm{\Psi}+{k}_{0}\frac{{n}_{2}}{{n}_{b}}{|\mathrm{\Psi}|}^{2}\mathrm{\Psi}+{k}_{0}\frac{n(X,Y)-{n}_{b}}{{n}_{b}}[1+\mu \text{\hspace{0.17em}}\mathrm{cos}(\delta Z)]\mathrm{\Psi}=0,$$## Eq. (2)

$$i\frac{\partial}{\partial z}\psi +\frac{1}{2}(\frac{{\partial}^{2}}{\partial {x}^{2}}+\frac{{\partial}^{2}}{\partial {y}^{2}})\psi +\sigma |\psi {|}^{2}\psi +V[1+\mu \text{\hspace{0.17em}}\mathrm{cos}(dz)]\psi =0,$$Different materials are used to produce waveguides, and silica is one of the most popular, with the typical parameters ${n}_{b}=1.4$, $|n-{n}_{b}|\le 9\times {10}^{-3}$, and ${n}_{2}=3\times {10}^{-16}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{2}/\mathrm{W}$ for light beams with the wavelengths ranging from the visible to the near-infrared. Without loss of generality, we choose ${\lambda}_{0}=800\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$ in this work. Therefore, if we choose ${V}_{0}=500$, the value of ${r}_{0}\approx 25.0\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ is obtained, according to the relation adopted in Eq. (2). Indeed, such a value is reasonable for a multimode fiber.^{24} According to the wavelength and ${r}_{0}$, one learns that the diffraction length is $\sim 7\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$. Considering that the group velocity dispersion coefficient is $\sim 35\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{fs}}^{2}/\mathrm{mm}$ at the wavelength of ${\lambda}_{0}=800\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$, the dispersion length is on the order of kilometers for a picosecond light beam, which is much longer than the propagation distance taken in this work. As a result, it is safe to neglect temporal effects.

To start with, we consider the modes supported by the deep potential alone; therefore, the nonlinear term and the longitudinal modulation in Eq. (2) are initially neglected. The corresponding solution of the reduced linear Eq. (2) can be written as $\psi (x,y,z)=u(x,y)\mathrm{exp}(i\beta z)$, where $u(x,y)$ is the stationary profile of the mode and $\beta $ is the propagation constant. Plugging this solution into the reduced Eq. (2), one obtains

## Eq. (3)

$$\beta u=\frac{1}{2}(\frac{{\partial}^{2}}{\partial {x}^{2}}+\frac{{\partial}^{2}}{\partial {y}^{2}})u+Vu,$$Therefore, to seek approximate azimuthons in weakly nonlinear waveguides, one takes the degenerate modes and makes a superposition of them as the initial wave

where $A$ is an amplitude factor, $1-B$ is the azimuthal modulation depth, and ${u}_{1,2}(x,y)$ are the degenerate linear modes (see examples in Fig. 1). Thus we set the transverse profile of the azimuthon at the initial place as and numerically propagate it to obtain an output mode at an arbitrary $z$. We should note that the inputs depicted by Eq. (5) do not rotate in a linear medium, since the modes are degenerate. Rotation appears only when the nonlinearity is added to the model.In Fig. 2, we display such approximate azimuthons with $A=0.4$ and $B=0.5$. One finds that the phase of azimuthons is nontrivial, displaying angular momentum and topological charge. For dipole azimuthons the topological charge is $\pm 1$, whereas for quadrupole, hexapole, and octopole azimuthons the values are $\pm 2$, $\pm 3$, and $\pm 4$, respectively. As expected, these azimuthons will rotate with a constant angular frequency $\omega $ during propagation when the nonlinear term in Eq. (2) is included. Therefore, the wave $U(x,y,z)$ can be rewritten as $U(r,\theta -\omega z)$ in polar coordinates, where $r=\sqrt{{x}^{2}+{y}^{2}}$ and $\theta $ is the azimuthal angle in the transverse plane $(x,y)$. This fact allows for a bit of theoretical analysis.

After plugging Eq. (4) into Eq. (2) with $\mu =0$, multiplying by ${U}^{*}$ and ${\partial}_{\theta}{U}^{*}$, respectively, and integrating over the transverse coordinates, one ends up with a linear system of equations:

## Eq. (6)

$$-\beta P+\omega {L}_{z}+I+N=0,\phantom{\rule{0ex}{0ex}}-\beta {L}_{z}+\omega {P}^{\prime}+{I}^{\prime}+{N}^{\prime}=0,$$After these preliminaries, we are ready to address the Rabi oscillation of azimuthons. To this end, we adopt the superposition of two azimuthons ${U}_{m,n}(x,y)\mathrm{exp}(i{\beta}_{m,n}z)$ as an input

## Eq. (8)

$$\psi ={c}_{m}(z){U}_{m}(x,y)\mathrm{exp}(i{\beta}_{m}z)+{c}_{n}(z){U}_{n}(x,y)\mathrm{exp}(i{\beta}_{n}z),$$## Eq. (9)

$$i\frac{\partial {c}_{m}}{\partial z}{U}_{m}\text{\hspace{0.17em}}\mathrm{exp}(idz)+\frac{1}{2}\mu {c}_{m}V{U}_{m}[1+\mathrm{exp}(2idz)]\phantom{\rule{0ex}{0ex}}+i\frac{\partial {c}_{n}}{\partial z}{U}_{n}+\frac{1}{2}\mu {c}_{n}V{U}_{n}[\mathrm{exp}(idz)+\mathrm{exp}(-idz)]=0.$$Since the azimuthons are constructed based on Eq. (5), they satisfy the relation $\u27e8{U}_{m},{U}_{n}\u27e9\ne 0$ if $m=n$ and $\u27e8{U}_{m},{U}_{n}\u27e9=0$ if $m\ne n$, thus forming a complete set of eigenstates. Here we borrowed the bra-ket notation from quantum mechanics. Note that the orthogonality of azimuthon shapes is only valid in the weakly nonlinear regime. As a result, one obtains two coupled equations based on Eq. (9):

## Eq. (10)

$$i\frac{\partial {c}_{m}}{\partial z}+\frac{1}{2}\mu \frac{\u27e8{U}_{m}V{U}_{n}\u27e9}{\u27e8{U}_{m}{U}_{m}\u27e9}{c}_{n}=0,\phantom{\rule{0ex}{0ex}}i\frac{\partial {c}_{n}}{\partial z}+\frac{1}{2}\mu \frac{\u27e8{U}_{n}V{U}_{m}\u27e9}{\u27e8{U}_{n}{U}_{n}\u27e9}{c}_{m}=0,$$## Eq. (12)

$${\mathrm{\Omega}}_{R}=\frac{\mu}{2}\frac{\u27e8{U}_{m}V{U}_{n}\u27e9}{\sqrt{\u27e8{U}_{m}{U}_{m}\u27e9\u27e8{U}_{n}{U}_{n}\u27e9}}.$$We note that the azimuthon conversion happens at half of the period, i.e., at ${z}_{R}/2$. Note also that the Rabi spatial frequency directly depends on the modulation strength $\mu $.

## 2.2.

### Circular Waveguide

We investigate the propagation of azimuthons in the circular weakly nonlinear waveguide by including the longitudinal modulation, and the results are displayed in Fig. 3. Without loss of generality, we choose the dipole and hexapole azimuthons, which are shown in Figs. 3(a) and 3(b), respectively. By taking the dipole azimuthon as an example [Fig. 3(a)], we want to see whether the Rabi oscillation between the dipole azimuthon [Fig. 2(a)] and its corresponding second-order dipole azimuthon [Fig. 2(e)] can be established. So, to induce resonance, we set the modulation frequency to be the difference between the eigenvalues of the two modes, which is $d\approx 25.2$. As shown by Eq. (12), the period of the Rabi oscillation is expected to depend on the modulation strength $\mu $, and here we set it to be $\mu \approx 0.031$ to make the period ${z}_{R}\sim 60$. As a consequence, one expects to see the second-order dipole azimuthon at a distance $\sim 30$ after turning on the longitudinal modulation.

In Fig. 3(a), the propagation of the dipole azimuthon is exhibited as a three-dimensional (3-D) isosurface plot, in which the longitudinal modulation exists only in the interval $30\le z\le 90$. When the propagation distance is smaller than $z\le 30$, one in fact observes the stable rotational propagation of the dipole azimuthon. The selected amplitude distributions at $z=0$ and $z=30$ are shown above the 3-D isosurface plots. In the interval $30\le z\le 90$, which is about one period of the Rabi oscillation, the oscillation between the dipole azimuthon and the second-order azimuthon is displayed, in which the dipole azimuthon completely switches to the second-order azimuthon at $z=60$. Indeed, the corresponding amplitude distribution is the same as that in Fig. 2(e) except for rotation, and the reason is quite obvious—azimuthons rotate steadily during propagation. When the propagation distance reaches $z=90$, the dipole azimuthon is recovered and the longitudinal modulation is also lifted at the same time. Therefore, one observes a stable rotating dipole azimuthon in the interval $90\le z\le 120$ and the amplitude distributions at $z=90$ and $z=120$, which are dipole azimuthons explicitly, and are shown above the isosurface plot. The analogous propagation dynamics of the hexapole azimuthon is shown in Fig. 3(b), the setup of which is the same as that of Fig. 3(a); it also clearly displays the Rabi oscillation of a higher-order azimuthon.

Here we would like to note that the Rabi oscillation is not feasible between two arbitrary azimuthons. Only azimuthons with similar structures (e.g., the dipole and the higher-order dipole azimuthons) can switch between each other, and azimuthons with different symmetries (e.g., the dipole and the quadrupole azimuthons) will not, because the overlap integrals in general are zero, $\u27e8{U}_{m}V{U}_{n}\u27e9=0$. Note also that the Rabi oscillation between two modes with opposite symmetry is also possible if the potential is antisymmetrically modulated in the transverse plane.^{25}

Generally, there must exist some frequency detuning $\ell =d-{d}^{\prime}$ between the real modulation frequency ${d}^{\prime}$ and the resonant frequency $d$. Therefore, it is reasonable to have a look at the efficiency of the azimuthon conversion versus the detuning $\ell $. Unfortunately, one cannot obtain a direct measure of the efficiency via the projection of the field amplitude $\psi $ on the targeting azimuthons due to the rotation of azimuthons during propagation. However, the efficiency is reflected on the Rabi oscillation period ${z}_{R}$—the bigger the value of ${z}_{R}$, the bigger the efficiency.^{13}^{–}^{15} The dependence of ${z}_{R}$ on the frequency detuning $\ell $ is shown in Fig. 3(c). As expected, one finds that the efficiency of azimuthon conversion is the biggest at the resonant frequency, and it reduces with the growth of the frequency detuning $\ell $.

## 2.3.

### Square Waveguide

Now we investigate the azimuthon transition in a square waveguide, which is generated by the potential in Eq. (2) in the form $V(x,y)={V}_{0}\text{\hspace{0.17em}}\mathrm{exp}[-({x}^{10}+{y}^{10})/{w}^{10}]$. Again, we solve for the linear eigenmodes supported by the deep square waveguide using the plane-wave expansion method. Connected with the geometry of the potential, the amplitude distributions of the linear modes are more complex than those in the regular circular waveguide; therefore, we denote them as the deformed modes. In Fig. 4, we display three kinds of deformed modes and the corresponding azimuthons, which will transform mutually because of the relation $\u27e8{U}_{m}V{U}_{n}\u27e9\ne 0$.

Different from the azimuthons in circular waveguides, the azimuthons in square waveguides rotate conditionally. Due to the symmetry of the square waveguide, a rotating azimuthon will deform, i.e., its profile will change. Without considering nonlinearity, the linear superposition of two degenerate modes [e.g., dipoles in Fig. 4(a), labeled as ${u}_{1,2}(x,y)$] is also a dipole solution in the square potential, as ${u}_{1,2}^{\prime}(x,y)=[{u}_{1}(x,y)\pm {u}_{2}(x,y)]/\sqrt{\iint {|{u}_{1}(x,y)\pm {u}_{2}(x,y)|}^{2}\mathrm{d}x\text{\hspace{0.17em}}\mathrm{d}y}$. We found that an azimuthon in a square waveguide will rotate when the condition $\mathcal{H}(\psi )>\mathcal{H}({u}_{1,2}^{\prime})$ is met. Here $\mathcal{H}$ represents the Hamiltonian of the model.^{22} Even though the wobbling azimuthons can be established in our numerical simulations, we are more interested in the rotating azimuthons, so we set again $A=0.4$ and $B=0.5$ in this section to guarantee the rotation of azimuthons during propagation.

Thus the azimuthons will be deformed during propagation because of the symmetry of the potential, hence to exhibit the azimuthon conversion in a more clear way one has to properly choose the value of $\mu $ to make the Rabi oscillation period almost equal to one half of the target azimuthon rotation period, since the modulation frequency $d$ is determined beforehand. Numerical simulations reveal that the rotation periods of the hexapole azimuthon and the higher-order dipole azimuthon are $\sim 148.4$ and $\sim 146.4$, respectively. Therefore, according to Eq. (12), the modulation strength $\mu $ for the two cases should be $\sim 0.085$ and $\sim 0.034$, respectively. The numerical demonstration of the Rabi oscillation is shown in Fig. 5. Different from the setting for the case of a circular waveguide, the longitudinal modulation always accompanies the square waveguide.

As shown in Fig. 5(a), the hexapole azimuthon is obtained at $z\sim 37.1$ (one half of the Rabi oscillation period and also a quarter of the rotation period of the hexapole azimuthon), whereas in Fig. 5(b) the higher-order dipole azimuthon is obtained at $z\sim 36.6$. When the propagation distance reaches one Rabi oscillation period, the dipole azimuthon is recovered with a small deformation. To show the azimuthon conversion more transparently, we also display the corresponding phase distributions. Evidently, there is only one phase singularity in the phase of the dipole azimuthon (the topological charge is 1), five singularities for the hexapole azimuthon (the topological charge is 3), and nine for the higher-order dipole azimuthon (the topological charge is again 1). As seen, the phase distributions are in accordance with the expectations and with those displayed in Fig. 4.

## 3.

## Conclusion

We investigated and demonstrated Rabi oscillations of azimuthons in weakly nonlinear waveguides with weak longitudinally periodic modulations. Based on the coupled mode theory, we find the period of Rabi oscillation, which is affected by the modulation strength and also by the spatial symmetry of azimuthons. The analysis is feasible for both circular and square waveguides and can be extended to the waveguides of other symmetries.

Based on the model considered in this work, the switching between a vortex-carrying azimuthon and a multipole that is free-of-vortex will not happen. The reason is that the initial azimuthon is composed of two degenerate modes ${u}_{1,2}$ that will switch into another two degenerate modes ${u}_{3,4}$ during propagation. So the output is a composition of ${u}_{3,4}$, which also carries a vortex. However, if the potential is modulated transversely in a proper manner, such a switch becomes possible once both the longitudinal and transverse phase matching is satisfied.

## Acknowledgments

This work was supported by Guangdong Basic and Applied Basic Research Foundation (No. 2018A0303130057), the National Natural Science Foundation of China (Nos. U1537210, 11534008, and 11804267), and the Fundamental Research Funds for the Central Universities (Nos. xzy012019038 and xzy022019076). M. R. B. acknowledges support from the NPRP 11S-1126-170033 project from the Qatar National Research Fund, whereas K. C. J. and Y. Q. Z. acknowledge the computational resources provided by the HPC platform of Xi’an Jiaotong University. All authors acknowledge anonymous reviewers for their valuable comments. The authors declare no conflicts of interest.

## References

## Biography

**Kaichao Jin** received his BS degree from Wuhan University of Technology in 2018. He is a master student working under the supervision of Yiqi Zhang at the School of Electronic Science and Engineering of Xi’an Jiaotong University. Currently, he is working on topological and waveguide array physics.

**Yongdong Li** received his PhD from Xi’an Jiaotong University, Xi’an, China, in 2005. He is currently a professor at the School of Electronic Science and Engineering of Xi’an Jiaotong University. His research interests include modeling and applications of plasmas.

**Feng Li** received his PhD from the University Nice Sophie Antipolis, France, in 2013. He is currently a professor at the School of Electronic Science and Engineering of Xi’an Jiaotong University. His research interests include quantum and topological effects in photonic micro/nanostructures.

**Milivoj R. Belić** received his PhD from The City College of New York in 1980. He is a professor in physics at Texas A&M University at Qatar. His research interests include topological photonics, nonlinear optics, and nonlinear dynamics.

**Yanpeng Zhang** received his PhD from Xi’an Jiaotong University, Xi’an, China, in 2000. He is currently a professor at the School of Electronic Science and Engineering of Xi’an Jiaotong University. His research interests include quantum optics, topological photonics, and nonlinear optics.

**Yiqi Zhang** received his PhD from Xi’an Institute of Optics and Precision Mechanics, Chinese Academy of Sciences, Xi’an, China, in 2011. He is currently an associate professor at the School of Electronic Science and Engineering of Xi’an Jiaotong University. His research interests include topological photonics and nonlinear optics.