Fast numerical estimations of axisymmetric multimode fibers modes
Estimating the propagation constants and the transverse mode profiles of multimode fibers is not as easy as it sounds. In our recent work we highlighted here, we needed to estimate the mode profiles for a standard graded-index fiber. It turned out that many standard approximations done in the literature to estimate the propagation constants do not give results accurate enough for the mode profile. The general approach we introduced in a previous tutorial to numerically find the fiber modes for any index profile using a 2D scalar finite differences approach is still valid. However, to provide accurate results, it needs a fine discretization of the space that leads to important memory and computational time requirements when the fiber core increases. If we consider an axisymmetric fiber, we can obtain a 1D formulation of the problem, that is unfortunately unstable under naive finite differences approaches. We detail here a stable formulation that leads to accurate and fast estimations of the mode profiles.
This tutorial is mostly based on the Supplementary Information of our recent work [M.W. Matthès et al., arxiv, 2010.14813 (2020)] that is largely inspired by the work from [L.S. Tamil, Appl. Opt., 9 (1991)].
We want to estimate the modes profiles of a perfect straight graded-index fiber under the scalar approximation. Graded-index fiber mode profiles and dispersion relation do not have a closed-form analytical expression. However, approximate analytical expressions can be found, for instance, using perturbation theory or a variational approach (see for instance [A. Sharma and P. Bindal, LAMP Series Report (1992)] for a nice introduction). Arguably the most widely used approximation is the Wentzel-Kramers-Brillouin (WKB) approximation. It leads to an analytical dispersion relation when assuming an infinite quadratic spatial profile of the refractive index. While leading to accurate estimations of the propagation constants, it has limited accuracy for the expression of the spatial mode profiles ( see [A. Gedeon, Opt. Commun., 3, (1974)] and [L. Maksymiuk and G. Stepniak, Opt. Quant. Electron., 1 (2016)]). Finite difference methods are easy to implement numerically in the scalar apporoximation (see our tutorial), but the 2D discretization of the field leads to high memory requirements and computation times, and could also lead to inaccuracies for high order modes. Because we consider axiosymmetric index \(n(r)\) profiles, we want to simplify the system to solve a 1D problem that only depends on the radial coordinate \(r\), allowing us to increase the accuracy and decrease the computation time.
Moreover, solving directly the 2D case rises the problem of degeneracy, where there is an infinite number of possibilities for any spatially degenerate sub-space. For graded-index fiber, the size of the group of (quasi-)degenerate modes increases for high order modes, as it corresponds to modes for which the value of \(l+m\) is the same. Technically, all possible combinations of these groups of degenerate modes are modes of the fiber, so all these basis are equivalent. However, from a practical point of view, it is convenient to work in a particular basis. Typically, for graded-index fibers, Linearly Polarized (LP) modes, Orbital Angular Momentum (OAM) modes, and Laguerre-Gauss are different valid representations. So ideally, we want to generate specifically the same representation for all the modes, not a random combination, which is what a 2D (or 3D) solver will output. In the following, we will focus on the OAM modes.
The 2D scalar Helmholtz equation for a propagating mode can be written in the cylindrical coordinate system as
\begin{equation}
\partial_r^2 \psi(r,\phi) + \frac{1}{r}\partial_r \psi(r,\phi)+ \frac{1}{r^2}\partial_\phi^2 \psi(r,\phi)
+ \left[ k_0^2 n^2(r) - \beta^2\right]\psi(r,\phi) = 0\,\,, \tag{1}
\label{eq:hel_radial}
\end{equation}
where \(\psi\) is the optical field, \(\phi\) is the azimuthal coordinate, \(\beta\) is the propagation constant, and \(k_0 = 2\pi/\lambda\) with \(\lambda\) the wavelength.
Because the refractive index only depends on the radial coordinate \(r\) for a perfect graded-index fiber, we can separate the variables \(r\) and \(\phi\). We are then looking for the orbital angular momentum modes of the form:
\begin{equation}
\psi_{ml} = f_l(r) e^{i m \phi}\,\,. \tag{2}
\end{equation}
with \(l\) the radial order and \(m\) the azimuthal order, which also corresponds to the orbital angular momentum. Injecting this expression in equation \ref{eq:hel_radial} leads to the 1D equation
\begin{equation}
d_r^2 f_l(r) + \frac{1}{r}d_r f_l(r) + \left[ k_0^2 n^2(r) - \beta^2 - \frac{m^2}{r^2}\right]f_l(r) = 0 \tag{3}
\label{eq:hel_oam}
\end{equation}
The singularity at \(r = 0\) arising from the \(\frac{1}{r}\) term makes direct finite difference methods unstable. We can use the transformation:
\begin{equation}
g_l(r) = \frac{1}{f_l(r)} d_r f_l(r), \tag{4}
\label{eq:transformation}
\end{equation}
and rewrite equation \ref{eq:hel_oam} as a quadratic Ricatti equation:
\begin{equation}
d_r g_l(r) + P(r) + Q(r)g_l(r) + g^2_l(r) = 0, \tag{5}
\end{equation}
where
\begin{align}
Q(r) &= \frac{1}{r}, \tag{6}\\
P(r) &= k_0^2 n^2(r) - \beta^2 - \frac{m^2}{r^2}. \tag{7}
\end{align}
A finite-difference approximation of such an equation leads to the recursive expression:
\begin{equation}
1+h g_l^{n+1} = \frac{1}{1+hQ_n/2}\left(h^2 P_n-2+\frac{1-h Q_n/2}{1+h g_l^n}\right), \tag{8}
\end{equation}
where \(g_l^{n} = g_l(r_n)\), \(Q_{n} = Q(r_n)\), \(P_{n} = P(r_n)\), and \(h = r_{n+1} - r_n\) is the step size. The expression is formulated this way in [L.S. Tamil, Appl. Opt., 9 (1991)], a general solution for the Riccati's equation can be found here [G. File and T. Aga, Egypt. j. basic appl. sci., 3 (2016)].
The expression \ref{eq:transformation} can then be discretized as:
\begin{equation}
f_l^{n+1} = f_l^{n} \left( 1+h g_l^{n}\right). \tag{9}
\end{equation}
To find the first steps to initialize the iteration, we need to consider the boundary conditions at the center of the fiber core:
\begin{align}
\left.\frac{df}{dr}\right|_{r=0} = 0 & \quad \text{for } m = 0, \tag{10}\\
f(r=0) = 0 & \quad \text{for } m \neq 0. \tag{11}
\end{align}
For \(m=0\), we discretize the functions at \(r_n = n h-1/2\), and initialize the functions with \(f_l^0 = 1\) and \(g_l^0 = 0\). For \(m\neq0\), we discretize the the functions at \(r_n = n h\), and initialize the functions with \(f_l^0 = 0\), \(f_l^1 = h\) and \(g_l^1 = (1-h^2P_1)/h\). For a given value of \(m\), the propagation constants \(\beta_{ml}\) that satisfy the Helmholtz equation, corresponding to the propagating modes, are the ones for which the field vanishes at large values of \(r\).
The steps to find the modes of the fiber are the following:
- We start with \(m=0\), and perform a coarse scan of the propagation constant values between \(\beta_{min} = k_0 n_{min}\) and \(\beta_{max} = k_0 n_{max}\).
- We choose \(r_N > a\) large enough to assume that the field at this point, and thus \(f_N\), should be vanishingly small. The number of times \(f_N(\beta)\) changes sign between \(\beta_{min} = k_0 n_{min}\) and \(\beta_{max} = k_0 n_{max}\) gives us the number of propagation modes for the current value of \(m\). It corresponds to the maximal radial number \(l\) admissible for the azimuthal number \(m\).
- We then use a binary search algorithm to find, at a minimum computational cost, the accurate admissible values of \(\beta\) for each \(l\), i.e. the values that minimize \(f_N\) under a given tolerance value.
- We then increment the value of \(m\), and repeat the procedure.
- We stop when no solution is found for the current value of \(m\).
We implemented this technique in our module Python pyMMF. In the next tutorial, we will compare the numerical results to full vectorial simulations and to the 2D FDTD approach.