<figure><img src="image_1.png" /></figure>

HAL
open science

# Efficient multimode vectorial nonlinear propagation solver beyond the weak guidance approximation

Pierre Béjot

► To cite this version:

Pierre Béjot. Efficient multimode vectorial nonlinear propagation solver beyond the weak guidance approximation. Journal of the Optical Society of America B, 2024, 41 (5), pp.1160-1169. 10.48550/arXiv.2402.16492 . hal-04498495

HAL Id: hal-04498495

https://hal.science/hal-04498495v1

Submitted on 11 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An efficient multimode vectorial nonlinear propagation solver beyond the weak guidance approximation

PIERRE BÉJOT¹, *

¹Laboratoire Interdisciplinaire Carnot de Bourgogne (ICB), UMR6303 CNRS–Université de Bourgogne, 21078 Dijon, France

*Corresponding author: pierre.bejot@u-bourgogne.fr

Compiled February 8, 2024

In this article, we present an efficient numerical model able to solve the vectorial nonlinear pulse propagation equation in circularly symmetric multimode waveguides. The algorithm takes advantage of the conservation of total angular momentum of light upon propagation and takes into account the vectorial nature of the propagating modes, making it particularly relevant for studies in ring-core fibers. While conventional propagation solvers exhibit a computational complexity scaling as $N_{\text{mode}}^4$ where $N_{\text{mode}}$ is the number of considered modes, the present solver scales as $N_{\text{mode}}^{3/2}$. As a first example, it is shown that orbital angular momentum modulation instability processes take place in ring-core fibers in realistic conditions. Finally, it is predicted that the modulation instability process is followed by the appearance of breather-like angular structures. © 2024 Optical Society of America

http://dx.doi.org/10.1364/ao.XX.XXXXXX

## 1. INTRODUCTION

Advancements in optical fiber technology have sparked renewed interest in the potential use of multimode waveguides to address limitations inherently associated with monomode propagation. The abundance of potential waveguide architectures now available opens doors to exploring an extremely rich variety of propagation dynamics [1, 2]. Solitons [3–8] and associated resonant radiations through geometric parametric instabilities [9], invariant space-time wavepackets [10–14], self-cleaning [15, 16] and related condensation/thermalization [17–22] phenomena are some nonlinear fascinating effects, among others, taking place in multimode fibers that have been intensively studied during the last few years. Beside these nonlinear phenomena, the use of multimode fibers for increasing the number of data-transmitting channels lets envision the possibility to overcome the expected data-carrying capacity crunch of single-mode fibers. In this regard, specially designed multimode optical fibers, called ring-core fibers, have been developed so as to minimize spatial-mode couplings [23–25]. The originality of these fibers lies in the modal content they support. First, contrary to conventional multimode fibers, these fibers only support one single kind of ring-shaped radial modes (p-modes) with well separated propagation constants, which has the advantage to strongly limit mode coupling. The other strong advantage of the latter property is that this reduces the initially full 3D+1 problem ([r, θ, t, z], where (r, θ, z) are the cylindrical coordinates, z being the direction along which the light propagates, and t stands for time) into a far simpler 2D+1 (θ, t, z) problem. In a theoretical point of view, this dimensionality reduction lets envision the possibility to study the existence

of new kind of solutions, otherwise impossible to find in more complex waveguides. Another important property of the modes supported by these fibers is that they are also eigenvectors of the total angular momentum (TAM) operator. This is in fact true not only for ring-core fibers but also for any circularly symmetric waveguides. However, the peculiarity of the ring-core fiber modes lies in their polarization content. While modes of conventional fibers are in very good approximation circularly polarized (excepted for the special case of transverse-electric and transverse-magnetic modes), the ones of ring-core fibers are extremely vectorial [26], i.e., their polarization contents strongly vary within their transverse sections. Theoretically speaking, the immediate consequence of this point is that scalar propagation equation does not hold anymore in such a system. Accordingly, if one wants to simulate the (nonlinear) propagation within these fibers, one necessarily has to use vectorial propagation solvers, which obviously increases the complexity of the simulations. The latter can become strongly computationally demanding, in particular when the waveguide is massively multimode. In this context, it is then of paramount importance to develop fast and accurate full vectorial models able to capture as close as possible the nonlinear propagation dynamics in waveguiding structures potentially supporting hundreds of modes.

Recently, we developed a numerical algorithm able to efficiently solve the scalar version of the nonlinear unidirectional pulse propagation equation (UPPE) [10]. This algorithm, valid in the weak guidance approximation, was shown to be not only more accurate since based on UPPE, but also far faster [27] than the classical algorithm solving the multimode generalized nonlin-
ear Schrödinger equation (MM-GNLSE), even if the latter was massively parallelized [28]. This is because the complexity of MM-GNLSE scales approximately as $N^4$ [29], where N is the number of modes to handle while, as shown below, the developed algorithm approximately scales as $N^{3/2}$. Note however, that the superiority of our algorithm on MM-GNLSE ones in terms of calculation speed only holds when the number of modes becomes higher and higher. In this article, we present a vectorial version of our algorithm that remains valid beyond the weak guidance approximation, i.e., that takes into account the vectorial nature of the modes. We restrict this study to waveguides presenting a cylindrical symmetry around the propagation axis. In the first section, after recalling the propagation equation, we discuss the general properties of vectorial modes and their intimate link with the TAM of light. In particular, we will show how the conservation of TAM in circularly symmetric waveguide can be advantageously used for constructing a fast-modal transform, i.e., a transformation allowing to go back and forth between a modal representation of the field and its representation in the direct space. In the second section, we will present the developed algorithm for solving in an efficient way the vectorial nonlinear UPPE equation. Finally, in the last section we will present a numerical example. In particular, it is predicted, for the first time, that a moderately intense continuous field propagating in a ring-core fiber can evolve towards an angular breather-like structure.

## 2. LINEAR PROPAGATION EQUATION AND MODE STRUCTURE OF CIRCULARLY SYMMETRIC WAVEGUIDES

### A. Propagation equation

We consider here an isotropic non-magnetic medium without current and charges, whose linear susceptibility $\chi^{(1)}$ is inhomogeneous in space. At this point, since this section is devoted to study the modal structure of waveguides, we consider a linear propagation regime. The Maxwell's equations read:

$$ 
\begin{aligned}
\vec{\nabla} \times \vec{E} &= -\partial_t \vec{B}, \quad \vec{\nabla} \cdot \vec{D} = 0 \\
\vec{\nabla} \times \vec{B} &= \mu_0 \partial_t \vec{D}, \quad \vec{\nabla} \cdot \vec{B} = 0,
\end{aligned}
\quad (1) $$

where $\mu_0$ is the vacuum permeability, $\vec{E}$ and $\vec{B}$ are the electric and magnetic fields respectively. In the frequency space, the electric displacement field $\vec{D}$ respects

$$ \vec{D} = \epsilon_0 n^2 \vec{E}, \quad (2) $$

where $n = \sqrt{1 + \chi^{(1)}}$ is the frequency- and spatially-dependent refractive index and $\epsilon_0$ is the vacuum permittivity. Since $\vec{D}$ has a divergence equal to zero, it implies:

$$ \vec{\nabla} \cdot \vec{E} = -\vec{\nabla} [\log(n^2)] \cdot \vec{E}. \quad (3) $$

The electric field propagation equation then reads in the frequency space

$$ \vec{\Box} \vec{E} = \vec{0}, \quad (4) $$

with

$$ \vec{\Box} \vec{E} = \nabla^2 \vec{E} + \frac{n^2 \omega^2}{c^2} \vec{E} + \vec{S}[\vec{E}], \quad (5) $$

where $\omega$ is the angular frequency, $c$ is the speed of light in vacuum, and

$$ \vec{S}[\vec{E}] = \vec{\nabla} [\vec{\nabla} (\log n^2) \cdot \vec{E}]. \quad (6) $$

For a circularly symmetric waveguide (supposed homogeneous all along z), the refractive index only depends on r in the cylindrical coordinates system $(r, \theta, z)$ and, for simplicity, we define $f(r) = \log n^2$. In this case, one has

$$ \vec{\nabla} f(r) = \partial_r f(r) \vec{e}_r. \quad (7) $$

so that

$$ \vec{S}[\vec{E}] = \vec{\nabla} (\vec{\nabla} f \cdot \vec{E}) = \begin{pmatrix} \partial_r^2 f \tilde{E}_r + \partial_r f \partial_r \tilde{E}_r \\ \frac{1}{r} \partial_r f \partial_\theta \tilde{E}_r \\ \partial_r f \partial_z \tilde{E}_r \end{pmatrix}, \quad (8) $$

where the electric field $\vec{E}$ is decomposed in cylindrical coordinates as

$$ \vec{E} = \tilde{E}_r \vec{e}_r + \tilde{E}_\theta \vec{e}_\theta + \tilde{E}_z \vec{e}_z. \quad (9) $$

Physically speaking, the operator $\vec{S}[\vec{E}]$ is directly related to the spin-orbit coupling phenomenon, as we will discuss a bit latter.

### B. Total angular momentum of light: eigenvectors and eigenvalues

Let us consider the operators $L_z$ and $S_z$ acting on the electric field as:

$$ L_z = -i \frac{\partial}{\partial \theta} I_d, \quad S_z = i \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}, \quad (10) $$

where $I_d$ is the two-dimensional identity matrix, and $\theta$ the polar angle. The operator $L_z$ (resp. $S_z$) corresponds to the projection along z of the orbital (resp. spin) angular momentum of light. They can be seen as infinitesimal generators of rotations around the z axis of the field amplitude and polarization, respectively. While the eigenvectors of $S_z$ correspond to the left and right circularly polarized fields ($|\sigma_+\rangle, |\sigma_-\rangle$) with respective eigenvalues $(1, -1)$, the eigenvectors of $L_z$ are fields owing an helical wavefront $e^{i\ell\theta}$ ($\ell \in \mathbb{Z}$) with eigenvalues $\ell$. If one now considers the TAM of light, any eigenvectors of $J_z = L_z + S_z$ write

$$ |\Phi_{j_\ell}\rangle = e^{i\ell\theta} [a_+ |\sigma_+\rangle + a_- e^{2i\theta} |\sigma_-\rangle], \quad (11) $$

with corresponding eigenvalues $j_\ell = \ell + 1$ ($\ell \in \mathbb{Z}$) and $(a_+, a_-)$ the respective complex amplitudes of the two circularly polarized components. In particular, for a given $j_\ell$, one can define the two orthonormal eigenvectors $|\Psi_{\pm j_\ell}\rangle$

$$ |\Psi_{\pm j_\ell}\rangle = \frac{e^{i\ell\theta}}{\sqrt{2}} [ |\sigma_+\rangle \pm e^{2i\theta} |\sigma_-\rangle ], \quad (12) $$

that form a bi-dimensional basis of representation of any eigenvectors of $J_z$ with corresponding eigenvalue $j_\ell = \ell + 1$. In particular, one has:

$$ |\Phi_{j_\ell}\rangle = \left( \frac{a_+ + a_-}{\sqrt{2}} \right) |\Psi_{+,j_\ell}\rangle + \left( \frac{a_+ - a_-}{\sqrt{2}} \right) |\Psi_{-,j_\ell}\rangle. \quad (13) $$

It is also interesting for the following to express $|\Psi_{\pm j_\ell}\rangle$ in cylindrical coordinates ($\vec{e}_r, \vec{e}_\theta$):

$$ |\Psi_{+,j_\ell}\rangle = e^{ij_\ell\theta} \vec{e}_r, \quad |\Psi_{-,j_\ell}\rangle = ie^{ij_\ell\theta} \vec{e}_\theta. \quad (14) $$
### C. Total angular momentum eigenvectors: an orthonormal basis for representing the electric field

First, one has to emphasize that analysing only the transverse part of the field $\vec{E} \perp (r, \theta)$ is sufficient at this stage. This is because, first, the longitudinal part of the field $\vec{E}_z$ can be retrieved using the divergence constraint (Eq. 3) and secondly, because this component is generally speaking extremely weak in fibers. For simplicity, we will then use the shorthand notation $\vec{E}$ for the transverse part. Moreover, note that all derivations made below apply in the frequency domain. However, we will omit when not necessary the notation specifying that the field is expressed in the frequency domain.

Since the electric field is intrinsically angularly $2\pi$ periodic, both circularly polarized field components can be decomposed as a Fourier series in polar angle:

$$ \vec{E}(r, \theta) = \sum_{\ell} a_{+,\ell}(r) e^{i\ell\theta} |\sigma_+\rangle + a_{-,\ell}(r) e^{i\ell\theta} |\sigma_-\rangle. \quad (15) $$

In cylindrical coordinates, it writes

$$ \vec{E}(r, \theta) = \sum_{\ell} \frac{a_{+,\ell}(r)}{\sqrt{2}} e^{i(\ell+1)\theta} (\vec{e}_r + i \vec{e}_\theta) + \sum_{\ell} \frac{a_{-,\ell}(r)}{\sqrt{2}} e^{i(\ell-1)\theta} (\vec{e}_r - i \vec{e}_\theta), \quad (16) $$

which can be recast as

$$ \vec{E}(r, \theta) = \sum_{j\ell} A_{+,j\ell}(r) |\Psi_{+,j\ell}\rangle + A_{-,j\ell}(r) |\Psi_{-,j\ell}\rangle, \quad (17) $$

with $A_{+,j\ell}(r) = \frac{a_{+,\ell}(r) \pm a_{-,\ell+2}(r)}{\sqrt{2}}$. The eigenvectors $|\Psi_{+,j\ell}\rangle$ hence form an orthonormal basis, any electric field being represented as a unique weighted sum of TAM eigenvectors.

### D. Spin-orbit coupling

Spin-orbit coupling [30] can be highlighted by looking at the action of the operator $\vec{S}$ on a circularly polarized field embedding a well defined orbital angular momentum. Let us consider a circularly polarized electric field $\vec{E}_{(|\sigma_{\pm}\rangle, \ell)}$ embedding an orbital angular momentum $\ell$.

$$ \vec{E}_{(|\sigma_{\pm}\rangle, \ell)} = A(r) e^{i\ell\theta} |\sigma_{\pm}\rangle \quad (18) $$

In cylindrical coordinates, the electric field components ($E_r, E_{\theta}$) write

$$ \begin{aligned} E_r (|\sigma_{\pm}\rangle, \ell) &= A(r) e^{i(\ell \pm 1)\theta}, \\ E_{\theta} (|\sigma_{\pm}\rangle, \ell) &= \pm i A(r) e^{i(\ell \pm 1)\theta}. \end{aligned} \quad (19) $$

Injecting Eq. 19 in Eq. 8, one has

$$ \vec{S} [\vec{E}_{(|\sigma_{\pm}\rangle, \ell)}] = \begin{pmatrix} \partial_r^2 f A(r) + \partial_r f \partial_r A(r) \\ \frac{i(\ell \pm 1)}{r} \partial_r f A(r) \end{pmatrix} e^{i(\ell \pm 1)\theta}. \quad (20) $$

The vector $\vec{S}$ can be finally projected back in the circular polarization basis:

$$ \begin{aligned} S_{|\sigma_+\rangle} &= \frac{S_r - iS_\theta}{\sqrt{2}} e^{-i\theta} \\ S_{|\sigma_-\rangle} &= \frac{S_r + iS_\theta}{\sqrt{2}} e^{i\theta}, \end{aligned} \quad (21) $$

<figure><img src="image_3.png" /><figcaption>Fig. 1. Refractive index profile of the air-core fiber used in the present example (a) and propagation constants as a function of the TAM $j\ell$ of the vectorial modes (b). Only one group of radial modes is guided within the ring (blue squares and red circles) while the second (and higher) order radial modes are not (orange diamonds and green stars). The modes depicted as blue squares and red circles differ by their distinct polarization patterns. Panel (c) displays the second-order derivative of the guided modes propagation constants with respect to $j\ell$.</figcaption></figure>

where $S_r$ (resp. $S_\theta$) is the component along $\vec{e}_r$ (resp. $\vec{e}_\theta$) of $\vec{S}$ given in Eq. 20. One then notices that the term $\vec{S}$ induces a spin-orbit coupling through the refractive index gradient. More particularly, for a circularly polarized electric field ($|\sigma_{\pm}\rangle$), the component of opposite helicity (i.e., along $|\sigma_{\mp}\rangle$) of $\vec{S}$ is

$$ S_{|\sigma_{\mp}\rangle} = \frac{\partial_r^2 f A(r) + \partial_r f \partial_r A(r) \mp \frac{\ell \pm 1}{r} \partial_r f A(r)}{\sqrt{2}} e^{i(\ell \pm 2)\theta}. \quad (22) $$

One can readily deduce the following selection rules that apply in circularly symmetric media:

$$ \begin{aligned} (|\sigma_+\rangle, \ell) &\xrightarrow{\vec{S}} (|\sigma_-\rangle, \ell + 2), \\ (|\sigma_-\rangle, \ell) &\xrightarrow{\vec{S}} (|\sigma_+\rangle, \ell - 2), \end{aligned} \quad (23) $$

The operator $\vec{S}$ then does not conserve neither the spin nor the orbital angular momentum on their own but nevertheless conserves the projection along the z axis of the TAM operator $\mathcal{J} = \mathcal{L} + \mathcal{S}$ as we discuss now.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Electric field transverse profiles of the 7 first modes of the first group. Each column corresponds to a mode of well-defined TAM. The first row depicts the radial intensity profile of the circular left (solid blue) and right (dashed red) polarization component of the mode. The middle (resp. bottom) row depicts the spatial phase profile of the circular left (resp. right) component of the mode.</figcaption></figure>

### E. Total angular momentum conservation in circularly symmetric media

The conservation of TAM can be highlighted by looking at the action of the propagation operator $\vec{\vec{E}}$ (Eq. 5) on an electromagnetic field embedding a well defined TAM $j_\ell$. In this regard, it is sufficient to express the propagation equation for the two orthonormal eigenvectors of the TAM (Eq. 12). As described previously (Eq. 4), the propagation equation of the electric field spectrum $\vec{\vec{E}}$ in a non-magnetic, inhomogeneous, isotropic and linear medium reads

$$ \vec{\vec{E}} \vec{\vec{E}} = \vec{0}. \qquad (24) $$

Let us now consider the following two fields $\vec{\vec{E}}_{\pm, j_\ell}$

$$ \vec{\vec{E}}_{\pm, j_\ell} = A_{\pm, j_\ell}(r) |\Psi_{\pm, j_\ell}\rangle e^{iK_z z}, \qquad (25) $$

that own the same TAM $j_\ell$. Injecting them into Eq. 4, one obtains

$$ \nabla^2 \vec{\vec{E}}_{\pm, j_\ell} + \frac{n(r)^2 \omega^2}{c^2} \vec{\vec{E}}_{\pm, j_\ell} = \mathcal{D}_{+, j_\ell} \left[ A_{\pm, j_\ell}(r) \right] e^{iK_z z} |\Psi_{\pm, j_\ell}\rangle + \mathcal{D}_{-, j_\ell} \left[ A_{\pm, j_\ell}(r) \right] e^{iK_z z} |\Psi_{\mp, j_\ell}\rangle, \qquad (26) $$

with

$$ \begin{aligned} \mathcal{D}_{+, j_\ell} &= \left( \partial_r^2 + \frac{1}{r} \partial_r - \frac{j_\ell^2 + 1}{r^2} + \frac{n(r)^2 \omega^2}{c^2} - K_z^2 \right), \\ \mathcal{D}_{-, j_\ell} &= \frac{2j_\ell}{r^2}. \end{aligned} \qquad (27) $$

and, in the case of circularly symmetric inhomogeneous media,

$$ \begin{aligned} \vec{\vec{S}}[\vec{\vec{E}}_{+, j_\ell}] &= S_{+, j_\ell}[A_{\pm, j_\ell}(r)]e^{iK_z z}|\Psi_{\pm, j_\ell}\rangle \\ &\quad + S_{-, j_\ell}[A_{\pm, j_\ell}(r)]e^{iK_z z}|\Psi_{\mp, j_\ell}\rangle, \qquad (28) \\ \vec{\vec{S}}[\vec{\vec{E}}_{-, j_\ell}] &= \vec{0}. \end{aligned} $$

with

$$ \begin{aligned} S_{+, j_\ell} &= \partial_r^2 f + \partial_r f \partial_r, \\ S_{-, j_\ell} &= \frac{\ell+1}{r} \partial_r f. \end{aligned} \qquad (29) $$

The action of the propagation operator hence couples only states owing the same TAM, confirming thereby that the propagation operator $\vec{\vec{E}}$ preserves the TAM of light. Other said, the electromagnetic fields that are invariant at a given absolute phase over the propagation (i.e., optical modes) in circularly symmetric inhomogeneous media are also necessarily eigenvectors of TAM. This remark then implies that optical modes necessarily take the following functional form

$$ \vec{\vec{E}}_{\text{mode}}(r, \theta, z) = \left( A_{+,j_\ell}(r) |\Psi_{+,j_\ell}\rangle + A_{-,j_\ell}(r) |\Psi_{-,j_\ell}\rangle \right) e^{iK_z z}, \qquad (30) $$

which can be recast in cylindrical coordinates as

$$ \vec{\vec{E}}_{\text{mode}}(r, \theta, z) = e^{ij_\ell \theta} [A_{r,\ell}(r) \vec{\vec{e}}_r + iA_{\theta,\ell}(r) \vec{\vec{e}}_\theta] e^{iK_z z}. \qquad (31) $$

Before presenting how the radial part of the modes can be evaluated, it is interesting to have a close look on the structure of the matrix coupling $|\Psi_{+,j_\ell}\rangle$ and $|\Psi_{-,j_\ell}\rangle$. Schematically, according to Eqs. (27-28), the coupling matrix $C_{j_\ell}$ has the following form:

$$ C_{j_\ell} = \begin{pmatrix} d_+ + S_+ - \frac{j_\ell^2}{r^2} & j_\ell d_- \\ j_\ell (d_- + s_-) & d_+ - \frac{j_\ell^2}{r^2} \end{pmatrix}, \qquad (32) $$

where $d_+ = \partial_r^2 + \frac{1}{r}\partial_r - \frac{1}{r^2} + \frac{n(r)^2\omega^2}{c^2}$, $d_- = \frac{2}{r^2}$, and $s_- = \frac{1}{r}\partial_r f$. Calculating the optical modes amounts to solve the equation system $C_{j_\ell}[\phi_{j_\ell}] = K_z^2[\phi_{j_\ell}]$, where $\phi_{j_\ell}$ is an eigenvector of $C_{j_\ell}$ (and consequently a mode of the waveguide) and $K_z^2$ is its associated eigenvalue that corresponds to the square of the modal propagation constant. First, it is clear that the determinant of $C_{j_\ell}$ only depends on $j_\ell^2$. As a consequence, the modal propagation constants will not depend on the sign of $j_\ell$. This essential degeneracy, linked to the system symmetry, reflects that there is no preferential rotation direction for the field. Moreover, for not too strong refractive index gradient, one has $S \ll D$. Looking at Eq. 32 for $j_\ell \neq 0$, in the limiting case $S \to 0$, optical modes will tend to the following functional forms:

$$ \begin{aligned} \vec{\vec{E}}_{\text{mode}}(r, \theta, z) &\approx A(r)e^{ij_\ell \theta} [\vec{\vec{e}}_r \pm i \vec{\vec{e}}_\theta] e^{iK_z z}, \\ &\approx A(r)e^{i(j_\ell \mp 1)\theta} e^{iK_z z} |\sigma_\pm\rangle. \end{aligned} \qquad (33) $$
The modes $j_\ell \neq 0$ are then in good approximation circularly polarized owing a well defined orbital angular momentum. Note that this last remark only holds for fiber exhibiting a not too strong refractive index gradient, such as conventional fibers. This is not the case, for instance, for ring-core fibers, as we will discuss a bit later. On the contrary, when $j_\ell = 0$, $|\Psi_{+,0}\rangle$ and $|\Psi_{-,0}\rangle$ are not coupled anymore. They are thus eigenvectors of the propagation equation but do not share the same eigenvalue, excepted in the bulk case. The modes for $j_\ell = 0$ will then be either radially (i.e., along $\vec{e}_r$) or azimuthally (i.e., along $\vec{e}_\theta$) polarized:

$$ \begin{aligned} \vec{E}_{+,j_\ell=0}(r, \theta, z) &= A_r(r)e^{iK_r z} \vec{e}_r, \\ \vec{E}_{-,j_\ell=0}(r, \theta, z) &= A_\theta(r)e^{iK_\theta z} \vec{e}_\theta, \end{aligned} \qquad (34) $$

with $K_r \neq K_\theta$. The next section is now devoted to describe how the modal transverse profile can be calculated efficiently.

## F. Determination of the radial distribution of the modes

As discussed before, optical modes of circularly symmetric waveguides are necessarily eigenvectors of the TAM operator so that their general functional forms follow Eq. 31. One has to emphasize here that the refractive index derivatives appearing in $C_{j\ell}$ has to be understood in the sense of distributions that make possible to differentiate functions whose derivatives do not exist in the classical sense. This point is of particular importance when one deals with step index fibers whose refractive index profiles exhibit discontinuities. As far as the radial distribution of the modes is concerned, it can be numerically evaluated for any refractive index distribution by using a convenient representation basis for the radial dimension. Note again that this procedure has to be repeated, if needed, for every frequencies independently. As it was done in our previous scalar algorithm [10], we decided here to use a radial basis set composed of zero order Bessel functions $J_0(k_\perp r)$. The main advantage of such a choice is that it is extremely easy to go back and forth between a direct space representation of the field and its representation in the Bessel basis (through a fast Hankel transform with its related unitary transformation matrix $H$ [31]). In the present case, the basis set used for the calculation is then composed of vector field functions $B$ that follow the general functional form given in Eq. 31

$$ B = e^{ij_\ell \theta} [J_0(k_\perp r) \vec{e}_r + iJ_0(k'_\perp r) \vec{e}_\theta] . \qquad (35) $$

Numerically speaking, the theoretically infinite and dense Bessel basis becomes finite and discrete. This is because a finite and discrete number of radial points is used. As a result, sampling the radial part of each field polarization component with $N_r$ spatial points implies the use of a Bessel radial basis containing $N_r$ vectors. The numerical procedure for calculating the optical modes consequently amounts, for any given TAM $j_\ell$, to diagonalize $C_{j\ell}$ expressed in the basis $B$ and whose dimension is $4N_r^2$ (the factor 4 coming from the fact that the subspace associated to a given $j_\ell$ is bi-dimensional). The diagonalization returns $2N_r$ orthonormal eigenvectors $v_{j\ell}$ that form an orthonormal basis $B_{eig}(j_\ell)$. The matrix $V_{j\ell}$ filled with the coordinates of the eigenvectors $v_{j\ell}$ in the basis $B$, is then the transformation matrix from $B$ to $B_{eig}(j_\ell)$. As a result, once an electric field is decomposed in the basis $B$ (which corresponds to the Hankel transform of the radial profiles of both polarization components), it is extremely easy to describe it in the modal basis. In the following section, we will see in more details how this can be done numerically in an efficient way.

## G. Fast modal transform

Let us assume that the electric field $\vec{E}$ to be decomposed in the eigenbasis is initially written in the direct space (in polar coordinates) and in the circular polarization basis

$$ \vec{E}(r, \theta) = a_+(r, \theta)|\sigma_+\rangle + a_-(r, \theta)|\sigma_-\rangle. \qquad (36) $$

Note that, numerically speaking, the coordinates $(r$ and $\theta)$ are sampled with a number of points $N_r$ and $N_\theta$ respectively. For simplicity, the field is written as a matrix of dimension $2N_r \times N_\theta$, where the first (resp. last) $N_r$ rows correspond to the left (resp. right) circular polarization component. Note that the angular sampling implies that only a finite number of TAM values is considered ($j_\ell \in [-\frac{N_\theta}{2}, \frac{N_\theta}{2} - 1]$). The field can be easily recast so as to be written along the vectors $\vec{e}_r$ and $\vec{e}_\theta$ with the following transformation $T$

$$ \vec{E}(r, \theta) = a_r(r, \theta) \vec{e}_r + ia_\theta(r, \theta) \vec{e}_\theta, \qquad (37) $$

with

$$ a_r(r, \theta) = \frac{e^{i\theta} a_+ + e^{-i\theta} a_-}{\sqrt{2}}, \quad a_\theta(r, \theta) = \frac{e^{i\theta} a_+ - e^{-i\theta} a_-}{\sqrt{2}}. \qquad (38) $$

One can now decompose the two field components as an angular Fourier series

$$ \begin{aligned} \vec{E} &= \sum_{j\ell} e^{j\ell\theta} [a_{r,j\ell}(r) \vec{e}_r + ia_{\theta,j\ell}(r) \vec{e}_\theta] \\ \vec{E} &= \sum_{j\ell} a_{r,j\ell}(r) |\Psi_{+,j\ell}\rangle + a_{\theta,j\ell}(r) |\Psi_{-,j\ell}\rangle, \end{aligned} \qquad (39) $$

where the coefficients $a_{r,j\ell}$ and $a_{\theta,j\ell}$ are easily obtained by making a fast Fourier transform $FFT_\theta$ of the field matrix along the second dimension (i.e., along $\theta$). After this operation, the field is represented, at each $r$, in the basis of the TAM operator (Eq. 14). Numerically speaking, the $J^{th}$ column (of length $2N_r$) of the matrix representing the field is then now $[a_{r,j\ell}(J), a_{\theta,j\ell}(J)]$. If one now wants to express the field in the modal basis, the last remaining operation is to multiply each column by the matrix $M_{j\ell} = V_{j\ell}H$, where $H$ is the Hankel transform unitary matrix [31]. Note that the matrix $M_{j\ell}$ to use is not the same for each column (i.e., for each $j_\ell$). Numerically speaking, it is then advantageous to reshape at this point the field matrix as a column vector (of length $2N_r N_\theta$) and define a sparse block diagonal matrix $M$ (of dimension $4N_r^2 N_\theta^2$) that is defined as

$$ M = \begin{pmatrix} M_{j_1} & (0) & (0) \\ (0) & \ddots & (0) \\ (0) & (0) & M_{j_{N_\theta}} \end{pmatrix} \qquad (40) $$

Doing so avoids a numerically heavy loop over $j_\ell$ and allows to express the field in the modal basis for every $j_\ell$ at once. After the operation, the field is expressed in the modal basis and can be reshaped back into a $2N_r \times N_\theta$ matrix. The global operation for going from the real space to the modal basis finally reads

$$ \bar{E}(p, j_\ell) = \mathcal{R}_2^{-1} M \mathcal{R}_2 (FFT_\theta [T \vec{E}(r, \theta)]), \qquad (41) $$

where $\mathcal{R}_2$ is the numerical operation for transforming the 2D matrix into a column vector and $\bar{E}(p, j_\ell)$ is the electric field expressed in the modal basis. Equation 41 constitutes the fast
<figure><img src="image_5.png" /><figcaption>Fig. 3. Electric field transverse profiles of the 7 first modes of the second group. Each column corresponds to a mode of well-defined TAM. The first row depicts the radial intensity profile of the circular left (solid blue) and right (dashed red) polarization component of the mode. The middle (resp. bottom) row depicts the spatial phase profile of the circular left (resp. right) component of the mode.</figcaption></figure>

modal transform that allows to go back and forth from the real space to the modal space. Note that, we considered here that the initial field was monochromatic. When dealing with pulsed fields, the transformation described above as to be performed in the frequency domain since the transformation matrix (Eq. 40) obtained during the diagonalization is now frequency dependent and denoted hereafter $M_\omega$. In this case, the modal transform of a field $\vec{E}(r, \theta, t)$ reads

$$ \vec{E}(p, j_l, \omega) = \mathcal{R}_3^{-1} M_\omega \mathcal{R}_3 \left( \text{FFT}_\theta \left[ \text{FFT}_t \left( T \vec{E}(r, \theta, t) \right) \right] \right), \quad (42) $$

where FFT$_t$ is the one-dimensional fast Fourier transform operated along the time dimension, $\mathcal{R}_3$ stands for the numerical operation transforming the 3D matrix representing the field into a column vector, and $M_\omega$ is the following sparse block diagonal transformation matrix (of dimension $4N_r^2 N_\theta^2 N_l^2$, where $N_l$ is the number of sampling points used in the time domain):

$$ M_{\omega} = \begin{pmatrix} M_{\omega_1} & (0) & (0) \\ (0) & \cdots & (0) \\ (0) & (0) & M_{\omega_{N_l}} \end{pmatrix} \quad (43) $$

The modal transform described above can be used for projecting any vector field distributions in the modal basis. As we will see later, it will be particularly useful for solving the nonlinear unidirectional pulse propagation equation.

## H. Numerical example

For the present example, we consider an air-core fiber whose (circularly symmetric) refractive index profile is depicted Fig. 1(a). Such kind of fiber is now well-known for allowing the propagation of modes embedding orbital angular momenta with extremely low cross-talks over long distances. We consider for this example a wavelength of 1.035µm. We used 250 (resp. 64) points for sampling the radial (resp. angular) coordinate, meaning that a total of 32000 modes (guided in the ring-core or not) have been calculated. As shown in Fig. 1(b), among all modes calculated, the considered fiber supports a total of 62 modes guided within the ring-core, all belonging to the same group of

radial modes. The guided modes can be divided into two distinct groups of modes, whose (spatially-dependent) polarization patterns are orthogonal. Each mode belonging to a particular group of modes has a distinct TAM $j_l$. Note that modes of the same group but with opposite TAM own exactly the same effective refractive index. Figure 2 (resp. Fig. 3) depicts the seven first guided modes of the first (resp. second) polarization group. As shown, all modes share almost the same radial distribution. Moreover, they are characterized by their strong vectorial nature, i.e., their spatially-dependent polarization pattern. However the polarization of the modes gradually approaches to a pure circular polarization as the TAM value increases. This can be seen by noticing that the amplitude of one polarization component, relatively to the other, becomes weaker and weaker as the TAM $|j_l|$ increases. As far as the phase spatial distribution is concerned (middle and bottom rows), one retrieves the expected signature of fields owing a well-defined TAM (Eq. 11). After having presented the algorithm for calculating the modes, their associated propagation constants and the modal transform, one can now describe how to solve the unidirectional pulse propagation equation.

## 3. NUMERICAL ALGORITHM FOR SOLVING THE NON-LINEAR UNIDIRECTIONAL PULSE PROPAGATION EQUATION

### A. Nonlinear unidirectional pulse propagation equation

When nonlinear effects are taken into account, the electric displacement field $\vec{D}$ reads in the frequency domain

$$ \vec{D} = \epsilon_0 (1 + \chi^{(1)}) \vec{E} + \vec{P}_{\text{NL}}, \quad (44) $$

where $\vec{P}_{\text{NL}}$ is the nonlinear polarization. The divergence of the electric field then now verifies

$$ \vec{\nabla} \cdot \vec{E} = - \left[ \vec{\nabla} (\log n^2) \cdot \vec{E} + \frac{\vec{\nabla} \cdot \vec{P}_{\text{NL}}}{\epsilon_0 n^2} \right], \quad (45) $$

where $n = \sqrt{1 + \chi^{(1)}}$ is the linear refractive index. As a result, the nonlinear propagation equation now reads in the frequency
<figure><img src="image_6.png" /><figcaption>Fig. 4. Computation time per propagation loop as a function of the number of radial (solid blue) and angular (dashed red) points. The solid blue (resp. dashed red) curve has been calculated with $N_θ=16$ (resp. $N_r=250$).</figcaption></figure>

domain

$$ \vec{\overleftrightarrow{E}} = -\frac{\omega^2 \vec{P}_{\text{NL}}}{\epsilon_0 c^2} - \frac{\vec{\nabla} (\vec{\nabla} \cdot \vec{P}_{\text{NL}})}{\epsilon_0 n^2}, \quad (46) $$

with $\vec{\overleftrightarrow{E}}$ the operator defined in Eq. 5. If one projects this equation on the modal space (i.e., on the basis composed of the vectorial modes of the linear propagation equation) and neglecting the backward field, one obtains the nonlinear UPPE equation

$$ \partial_z \vec{E} = i K_z \vec{E} + \frac{i}{2\epsilon_0 K_z} \left( \frac{\omega^2 \vec{P}_{\text{NL}}}{c^2} + \left[ \frac{\vec{\nabla} (\vec{\nabla} \cdot \vec{P}_{\text{NL}})}{n^2} \right] \right), \quad (47) $$

where $\vec{A}$ stands for the decomposition of a three-dimensional vector field $\vec{A}$ on the modal basis and $K_z$ is the associated propagation constant. Note here that Eq. 47 takes into account both the transverse field and its associated longitudinal component (i.e., the field component polarized along $z$). For a more detailed derivation of the nonlinear unidirectional pulse propagation equation, we invite the reader to refer to [32]. In the context of propagation within an optical fiber, a fair approximation is to neglect the last term in Eq. 47, and also the longitudinal component of the field since this part remains generally small as compared to its transverse part. As a result, its contribution to the nonlinear polarization will be neglected. As a result, we will now consider that the field is in fair approximation purely transverse (i.e., contained in the $[r, \theta]$ plane). However, one has to keep in mind that these approximations break when the field embeds modes that have propagation constants far smaller than the propagation constant in bulk ($K_z \ll \frac{n\omega}{c}$).

## B. Complex envelope equation

Most of the time, it is more convenient to deal with a propagation equation acting on the complex envelope of the field. This is done by defining the complex field $\vec{\epsilon}$ as

$$ \vec{E} = \frac{\vec{\epsilon} + \vec{\epsilon}^*}{2}. \quad (48) $$

For convenience, we also introduce the vector field $\vec{\xi}$ as

$$ \vec{\xi} = \sqrt{\frac{2}{\epsilon_0 c n}} \vec{\epsilon}, \quad (49) $$

so that the total intensity $I$ reads

$$ I = \vec{\xi} \cdot \vec{\xi}^*. \quad (50) $$

In the context of optical fiber propagation, we consider that the medium is in good approximation isotropic, at least, in a microscopic point of view. In this regard, the number of independent elements of the third-order nonlinear susceptibility tensor $\chi^{(3)}$ reduces to one. We can consequently define an unique nonlinear refractive index $n_2$ as

$$ n_2 = \frac{3\chi_{xxxx}^{(3)}}{4\epsilon_0 c n_0^2}, \quad (51) $$

where $n_0$ is the bulk refractive index taken at the central frequency of the considered field. Neglecting the term responsible for third harmonic generation and considering at this point a pure electronic instantaneous Kerr effect, the nonlinear polarization reads in the circular polarization basis

$$ \vec{P}_{\text{NL}} = \frac{2n_0 n_2}{3} \left( \frac{[|\vec{\xi}_+|^2 + 2|\vec{\xi}_-|^2] \vec{\xi}_+}{[2|\vec{\xi}_+|^2 + |\vec{\xi}_-|^2] \vec{\xi}_-} \right), \quad (52) $$

where $\vec{\xi}_+$ (resp. $\vec{\xi}_-$) is the left (resp. right) circular component of $\vec{\xi}$. Note that, if needed, the Raman-induced contribution to the nonlinear polarization can also be added. As one can remark, the nonlinear polarization term is easily calculated in the direct space (and in the circular polarization basis). Finally, it is of common use to consider a sliding time origin (corresponding, at any propagation distance $z$, to the arrival time of a pulse propagating at a group velocity $v_{g0}$). The complex unidirectional pulse propagation finally reads in the modal basis

$$ \partial_z \vec{\xi} = i \left( K_z - \frac{\omega}{v_{g0}} \right) \vec{\xi} + i \frac{n_0 \omega^2}{K_z c^2} \vec{P}_{\text{NL}}. \quad (53) $$

In the next section, we will see how Eq. 53 can be numerically resolved in an efficient way using the modal transform presented in section G.

## C. Numerical algorithm

The starting point for the numerical algorithm presented here is to note that the linear term can be easily evaluated in the modal space while it is far more convenient to calculate the nonlinear one in the direct space. As a result, the strategy for solving Eq. 53 is to use a split-step algorithm during which both terms (linear and nonlinear) are evaluated in a different representation space. In this regard, one can see the present method as a split-step modal algorithm, in analogy with the well-known split-step Fourier and split-step Hankel-Fourier algorithms. The proposed algorithm is then extremely close to the one we presented for solving the scalar version of UPPE [10]. The key point of the algorithm lies in the construction of the modal transform, described in section G (Eq. 42), for being able to efficiently go from the direct space to the modal space (and vice versa). Knowing the expression of a field $\vec{\xi}$ ($r, \theta, z, t$) at a given $z$, one can calculate the field at a close distance $z + dz$ as follow. First, one evaluates its expression $\vec{\xi}$ in the modal basis using Eq. 42. Then, a linear propagation step over $dz/2$ is performed so that

$$ \vec{\xi} \rightarrow e^{i \left( K_z - \frac{\omega}{v_{g0}} \right) \frac{dz}{2}} \vec{\xi}. \quad (54) $$

After this, the nonlinear polarization is evaluated in the direct space. This is done by evaluating the new field $\vec{\xi}$ in the direct space. Once calculated in the direct space, the nonlinear polarization is expressed in the modal basis with the modal
transform so as to solve the second term in Eq. 53 by a fourth-order Runge-Kutta algorithm. Finally, a second linear step over $dz/2$ is performed in the modal basis. Doing so avoids to deal with multiple summations of overlap integrals as it is usually done in simulations solving the multimode-generalized nonlinear Schrodinger equation. In order to assess the efficiency of the algorithm, we evaluated how the computation time increases with the number of modes involved in the code [which is equal to $2N_r N_\theta$, where $N_r$ (resp. $N_\theta$) is the number of points in the radial (resp. angular) dimension]. For this, we timed one hundred nonlinear propagation loops as a function of $N_r$ and $N_\theta$ in the monochromatic case ($N_t=1$). The calculations have been performed in a Dell Optiplex 5260 computer (32 Go RAM, processor Intel I7 8th generation). The results are summarized in Fig. 4. As shown, the propagation algorithm complexity is approximately $O(Nr^2N_\theta)$. The fact that the algorithm complexity is higher in $N_r$ than in $N_\theta$ comes from the fact that only fast Fourier transforms are performed along $\theta$ while the modal transform needs to deal with matrix vector multiplications in the radial domain. If one wants to evaluate how the code behaves with the number of modes involved $N_{\text{mode}}$, this can be done by setting $N_\theta = N_r$, which leads to a complexity $O(N_{\text{mode}}^{3/2})$, which is far more efficient than algorithms solving the MM-GNLSE, known to behave as $O(N_{\text{mode}}^4)$ [28, 29].

## 4. NUMERICAL EXAMPLE: ANGULAR MODULATION INSTABILITY AND ANGULAR BREATHER-LIKE GENERATION IN RING-CORE FIBER

Nonlinear effects in ring-core fibers have become a recent subject of active study, leading to effects never observed in conventional fibers [33, 34]. Here, using the code presented in this paper, we predict that, first, vector modulation instability takes place in ring-core-fibers, and second, that angular breather-like solutions can emerge in appropriate conditions. For this purpose, we consider the propagation of a continuous ($\lambda_0 = 1.035\mu$m) laser field propagating in the ring-core fiber considered in section H. The effective mode area for the transverse electric mode is approximately $150 \mu m^2$ and the associated effective nonlinear coefficient is $\gamma \simeq 1.3 W^{-1}km^{-1}$. The initial field $\vec{\zeta}(r, \theta, z = 0)$ is coupled in the transverse electric (fundamental) mode ($J = 0$) with a power $P_0=200$ kW.

$$ \vec{\zeta}(r, \theta, z = 0) = \sqrt{I_0 F(r)} \vec{e}_{\theta}, \qquad (55) $$

where $F(r)$ is the normalized radial profile of the transverse electric mode of the ring-core fiber. The peak intensity $I_0$ is defined so that

$$ \iint r \vec{\zeta} \cdot \vec{\zeta}^* drd\theta = P_0. \qquad (56) $$

and approximately equals to $185$ GW/cm$^2$. A weak random noise (both in phase and amplitude) is superimposed in the modal space to the initial field. Note that such an initial condition is experimentally realistic with femtosecond, picosecond and even nanosecond Ytterbium doped fiber amplifiers [35] now commercially available. In particular, the use of long pulses will limit effects, such as self-phase modulation and modal group velocity dispersion, that could be potentially detrimental for observing the effect presented below. In order to limit the computational time, we discarded the time dependence of the field. At the early stage of propagation, only two particular modes start to exponentially grow up. These modes belong to the pump beam vectorial group and have a TAM $j_\ell = \pm 1$

<figure><img src="image_7.png" /><figcaption>Fig. 5. Energy gain of the modes $j_\ell = \pm 1$ as a function of the propagation distance. The blue (resp. red) solid line corresponds to the mode $j_\ell = 1$ (resp. $j_\ell = -1$) of the first group, i.e., belonging to the vectorial group of the pump. The green (resp. yellow) dashed line depicts the mode $j_\ell = 1$ (resp. $j_\ell = -1$) belonging to the second vectorial group of modes.</figcaption></figure>

(see Fig. 5). This indicates that only TAM-conserved processes of kind $2J_{\text{pump}} \rightarrow J_1 + J_2$ take place at the very first stage of the propagation. Quickly after ($\approx 1$ cm), the two modes with same total angular momenta ($j_\ell = \pm 1$) but belonging to the second vectorial group starts to grow with more or less the same gain ($\approx 22$ dB/m) than the one experienced by the two first modes. During the 7 first centimeters of propagation, only these four modes are populated. If one now decomposes the field in left and right circularly polarized components and looking at their orbital angular momentum spectra [see Figs. 6(a,b)], only four contributions experiences gain ($|\ell = 0, \sigma_-\rangle, |\ell = 2, \sigma_-\rangle, |\ell = 0, \sigma_+\rangle$, and $|\ell = -2, \sigma_+\rangle$) during the first 7 cm. During this first step, the angularly-resolved power distribution remains flat [see Figs 6(c,d)]. Even if a rigorous theoretical study of this process is beyond the scope of the present paper, all these observations nevertheless well correspond to an angular analogue of the (time) modulation instability process taking place in single mode fibers in the anomalous dispersion regime [36]. The analogy is corroborated by noticing that the second (discrete) derivative of the propagation constant with respect to the TAM of the modes [Fig. 1(c)] is negative. In this regard, it is particularly interesting to note that the numerical gain experienced by the two first modes very well matches the analytical formula

$$ G(j_\ell) = \frac{20}{\log 10} \sqrt{-\frac{k_{j\ell}^{(2)}}{2} j_\ell^2 \left( \frac{k_{j\ell}^{(2)}}{2} j_\ell^2 + 2\gamma P_0 \right)}, \qquad (57) $$

where $k_{j\ell}^{(2)}$ is the second-order (discrete) derivative of the propagation constant with respect to $j_\ell$ [see Fig. 1(c)] evaluated for the transverse electric fundamental mode and where G is the energy gain evaluated in dB/m. Equation 57 is the angular analogue of the gain formula for modulation instability in the frequency domain. However, note that this formula is established only empirically. Further studies will be obviously needed so as to confirm or infirm its validity. The fact that angular modulation instability seems to take place suggests that angular analogue of solitons and breathers [37, 38] could be also observed. Note that such a possibility is made possible because there is only one single radial guided mode in such kind of fibers, which extremely simplifies the physical process. In this regard, the propagation of the initial field has been continued up to 20 cm. From approximately 10 cm, the orbital angular momentum spectrum of the field starts to strongly broaden for becoming maximal at a
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Orbital angular momentum spectrum (a) [resp. (b)] of the left (resp. right) circularly polarized component of the field as a function of the propagation distance. Angular distribution of the intensity profile (c) [resp. (d)] of the left (resp. right) circularly polarized component of the field as a function of the propagation distance. Energy gain of the fourth components growing during the first stage of propagation (e). Angular spectrum (f) of the left (blue squares) and right (red circles) circular components at the distance where the angular compression is maximal (z ≃ 13 cm) and associated angular distribution (g).</figcaption></figure>

distance ≈ 13 cm. At this particular distance, the gain saturates and the OAM spectrum of both circularly polarized components exhibits a triangular shape in logarithmic scale [Fig. 6(f)]. In the direct space, this *spectral* broadening is associated to an extreme angular localization of the power [Fig. 6(g)], whose shape turns out to be very similar to breather structures already observed in fibers in the time domain. Afterwards, the angular dynamics of both polarization components exhibits multiple periodic breathing. Here again, this latter dynamics is extremely similar to those already studied in time [39]. However, it is interesting to note that the periodic and finite nature of the angular dimension over which the dynamics takes place nevertheless implies some differences in behavior with their time-domain counterparts. For instance, the fact that the orbital angular momentum spectrum is discrete by nature implies that the modulation instability process has a power threshold below which the effect vanishes. This discrete nature also explained why angular breathers appear in the above example without the use of a deterministic and coherent seeding as needed for their time-domain analogues.

**5. CONCLUSION**

To conclude, an efficient numerical algorithm for solving the vectorial nonlinear unidirectional pulse propagation equation has been presented. The algorithm, valid beyond the weak guidance approximation for circularly symmetric waveguides, takes advantage of the TAM conservation for being able to calculate the vectorial optical modes. The modal calculation is used in turn so as to build a fast modal transform, which allows to easily switch from a representation of the field in the direct space to its decomposition in the modal space (and *vice versa*). The presented algorithm has a complexity that roughly behaves as

$N_{\text{mode}}^{3/2}$, which is far more efficient than conventional algorithms that scale as $N_{\text{mode}}^4$. As a first numerical example, we considered the nonlinear propagation of a moderately intense transverse electric field within a ring-core fiber. It was first shown that angular modulation instability is expected to take place, suggesting in turn the possibility to observe angular analogue of temporal solitons and breathers. This expectation has been then corroborated by exhibiting a realistic case where vectorial breather-like structures seem to emerge. While a rigorous study of such a phenomenon is out of the scope of the present paper, the perspective to observe angular breathers in ring-core fibers could nevertheless open the door to the study of new kind of nonlinear solutions. Beyond the presented monochromatic case, an extremely rich family of spatio-temporal (time+angle) dynamics could be also studied in such kind of systems. In particular, it could be interesting to see if bi-dimensional breather solutions can be found, either in the normal or anomalous dispersion (frequency) regimes. Note that since the number of angular modes is finite, we anticipate that the quest to find analytical bi-dimensional breather solutions in such a case will be far more tractable than in the general 3D+1 situation where both radial and angular dimensions have to be taken into account. Note, however, that the fact that optical modes are strongly vectorial could nevertheless complicate a bit the theoretical study.

**ACKNOWLEDGMENTS**

The author acknowledges financial support of French programs "Investments for the Future" operated by the National Research Agency (ISITE-BFC, contract ANR-15-IDEX-03; EIPHI Graduate School, contract ANR-17-EURE-0002; EQUIPEX+ Smartlight, contract ANR-21-ESRE-0040), from Bourgogne Franche-Comté
region and European Regional Development Fund. The author thanks the CRI-CCUB for CPU loan on the multiprocessor server.

## DISCLOSURES

The author declares no conflicts of interest.

## DATA AVAILABILITY

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## REFERENCES

1. K. Krupa, A. Tonello, A. Barthélémy, T. Mansuryan, V. Couderc, G. Millot, P. Grelu, D. Modotto, S. A. Babin, and S. Wabnitz, Multimode nonlinear fiber optics, a spatiotemporal avenue, APL Photonics 4, 110901 (2019).

2. L. G. Wright, F. O. Wu, D. N. Christodoulides, and F. W. Wise, Physics of highly multimode nonlinear optical systems, Nat. Physics 18, 1018-1030 (2022).

3. B. Crosignani and P. D. Porto, Soliton propagation in multimode optical fibers, Opt. Lett. 6, 329-330 (1981).

4. A. Hasegawa, Self-confinement of multimode optical pulse in a glass fiber, Opt. Lett. 5, 416-417 (1980).

5. B. Crosignani, A. Cutolo, P. D. Porto, Coupled-mode theory of nonlinear propagation in multimode and single-mode fibers: envelope solitons and self-confinement, J. Opt. Soc. Am. 72, 1136-1141 (1982).

6. S. S. Yu, C. H. Chien, Y. Lai, and J. Wang, Spatio-temporal solitary pulses in graded-index materials with Kerr nonlinearity, Opt. Commun. 119, 167-170 (1995).

7. S. Raghavan, and G. P. Agrawal, Spatiotemporal solitons in inhomogeneous nonlinear media, Opt. Commun. 180, 377-382 (2000).

8. W. H. Renninger, and F. W. Wise, Optical solitons in graded-index multimode fibres, Nat. Commun., 4, 1719 (2013).

9. K. Krupa, A. Tonello, A. Barthélémy, V. Couderc, B. M. Shalaby, A. Bedahmane, G. Millot, and S. Wabnitz, Observation of Geometric Parabolic Instability Induced by the Periodic Spatial Self-Imaging of Multimode Waves, Phys. Rev. Lett. 116, 183901 (2016).

10. P. Béjot, Multimodal unidirectional pulse propagation equation, Phys. Rev. E 99, 032217 (2019).

11. B. Kibler and P. Béjot, Discretized conical waves in multimode optical fibers, Phys. Rev. Lett. 126, 023902 (2021).

12. P. Béjot and B. Kibler, Spatiotemporal Helicon Wavepackets, ACS Photonics, 8(8), 2345-2354 (2021).

13. P. Béjot and B. Kibler, Quadrics for Structuring Invariant Space-Time Wavepackets, ACS Photonics 9(6), 2066-2072 (2022).

14. K. Stefanska, P. Béjot, K. Tarnowski, and B. Kibler, Experimental Observation of the Spontaneous Emission of a Space-Time Wavepacket in a Multimode Optical Fiber, ACS Phot. 10(3), 727-732 (2023).

15. Z. Liu, L. G. Wright, D. N. Christodoulides, and F. W. Wise, Kerr self-cleaning of femtosecond-pulsed beams in graded-index multimode fiber, Opt. Lett. Vol. 41(16), 3675-3678 (2016).

16. K. Krupa, A. Tonello, B. M. Shalaby, M. Fabert, A. Barthélémy, G. Millot, S. Wabnitz, and V. Couderc, Spatial beam self-cleaning in multimode fibres, Nat. Photonics 11, 237-241 (2017).

17. P. Aschieri, J. Garnier, C. Michel, V. Doya, and A. Picozzi, Condensation and thermalization of classsical optical waves in a waveguide, Phys. Rev. A 83, 033838 (2011).

18. A. Picozzi, J. Garnier, T. Hansson, P. Suret, S. Randoux, G. Millot, and D. N. Christodoulides, Optical wave turbulence: Towards a unified nonequilibrium thermodynamic formulation of statistical nonlinear optics, Phys. Reports 542, 1-132 (2014).

19. F. O. Wu, A. U. Hassan, and D. N. Christodoulides, Thermodynamic theory of highly multimoded nonlinear optical systems, Nat. Photonics 13, 776-782 (2019).

20. H. Pourbeyram, P. Sidorenko, F. O. Wu, N. Bender, L. Wright, D. N. Christodoulides, and F. Wise, Direct observations of thermalization to a Rayleigh-Jeans distribution in multimode optical fibres, Nat. Physics 18, 685-690 (2022).

21. F. O. Wu, Q. Zhong, H. Ren, P. S. Jung, K. G. Makris, and D. N. Christodoulides, Thermalization of Light's Orbital Angular Momentum in Nonlinear Multimode Waveguide Systems, Phys. Rev. Lett. 128, 123901 (2022).

22. E. V. Podivilov, F. Mangini, O. S. Sidelnikov, M. Ferraro, M. Gervaziev, D. S. Kharenko, M. Zitelli, M. P. Fedoruk, S. A. Babin, and S. Wabnitz, Thermalization of Orbital Angular Momentum Beams in Multimode Optical Fibers, Phys. Rev. Lett. 128, 243901 (2022).

23. N. Bozinovic, Y. Yue, Y. Ren, M. Tur, P. Kristensen, H. Huang, A. E. Willner, and S. Ramachandran, Terabit-Scale Orbital Angular Momentum Mode DivisionMultiplexing in Fibers, Science 340(6140), 1545-1548 (2013).

24. C. Brunet, P. Vaity, Y. Messaddeq, S. LaRochelle, and L. A. Rusch, Design, fabrication and validation of an OAM fiber supporting 36 states, Opt. Express 22(21), 26117-26127 (2014).

25. S. Ramachandran, P. Gregg, P. Kristensen and S. E. Golowich, On the scalability of ring fiber designs for OAM multiplexing, Opt. Express 3 (2015).

26. P. Gregg, P. Kristensen, and S. Ramachandran, Conservation of orbital angular momentum in air-core optical fibers, Optica 2(3), 267-270 (2015).

27. K. Tarnowski, S. Majchrowska, P. Béjot, and B. Kibler, Numerical modelings of ultrashort pulse propagation and conical emission in multimode optical fibers, J. Opt. Soc. Am. B 38(3), 732-742 (2021).

28. L. G. Wright, Z. M. Ziegler, P. M. Lushnikov, Z. Zhu, M. A. Eftekhar, D. N. Christodoulides, F. W. Wise, Multimode nonlinear fiber optics: massively parallel numerical solver, tutorial and outlook, IEEE J. Sel. Top. Quantum Electron. 24, 1-16 (2017).

29. F. Poletti and P. Horak, Description of ultrashort pulse propagation in multimode optical fibers, J. Opt. Soc. Am. B 25, 1645-1653 (2008).

30. K. Y. Bliokh, F. J. Rodriguez-Fortuno, F. Nori, and A. V. Zayats, Spin-orbit interactions of light, Nature Phot. 9, 796-808 (2015).

31. M. Guizar-Sicairos and J. C. Gutierrez-Vega, Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields, J. Opt. Soc. Am. A 21, 53 (2004).

32. J. Andreasen and M. Kolesik, Nonlinear propagation of light in structured media: Generalized unidirectional pulse propagation equations, Phys. Rev. E 86, 036706 (2012).

33. X. Liu, E. N. Christensen, K. Rottwitt, and S. Ramachandran, Nonlinear four-wave mixing with enhanced diversity and selectivity via spin and orbital angular momentum conservation, APL Photonics 5, 010802 (2020).

34. S. K. Dacha, W. Zhu, A. Agrawal, K. J. Ritter, and T. E. Murphy, Nonlinear rotation of spin-orbit coupled states in hollow ring-core fibers, Opt. Express 30(11), 18481-18495 (2022).

35. W. Torruellas, Y. Chen, B. McIntosh, J. Farroni, K. Tankala, S. Webster, D. Hagan, M. J. Soileau, M. Messerly, and J. Dawson, High peak power ytterbium doped fiber amplifiers, Proceedings Fiber Lasers III: Technology, Systems, and Applications 6102 61020N (2006).

36. P. Béjot, B. Kibler, E. Hertz, B. Lavorel, and O. Faucher, General approach to spatiotemporal modulational instability processes, Phys. Rev. A 83, 013830 (2011).

37. B. Kibler, J. Fatome, C. Finot, G. Millot, F. Dias, G. Genty, N. Akhmediev, and J. M. Dudley, The Peregrine soliton in nonlinear fibre optics, Nature Phys. 6,790-795 (2010).

38. J. M. Dudley, F. Dias, M. Erkintalo, and G. Genty, Instabilities, breathers and rogue waves in optics, Nat. Photonics. 8(1), 755-764 (2014).

39. M. Erkintalo, K. Hammani, B. Kibler, C. Finot, N. Akhmediev, J. M. Dudley, and G. Genty, Higher-Order Modulation Instability in Nonlinear Fiber Optics, Phys. Rev. Lett. 107, 253901 (2011).