Numerical Estimation of Multimode Fiber Modes and Propagation Constants: 

Part 2: Bent Fibers 

We saw in the first part of the tutorial that the profiles and the propagation constants of the propagation modes of a straight multimode fiber can easily be avulated for an arbitrary index profile by inverting a large but sparse matrix. Under some approximations [1], a portion of fiber with a fixed radius of curvature satisfies a similar problem that can be solved with the same numerical tools, as we illustrate with the PyMMF Python module [2]. Moreover, when the modes are known for the straight fiber, the modes for a fixed radius can be approximate by inverting a square matrix of size the number of propagating modes [1]. It allows fast computation of the modes for different radii of curvature.  

Effect of bending

Introduction

We base the following on the results of Ref. [1]. In this study, the authors state that even when the bending radius is very large, the changes on the effective index of the fiber are not that small, leading a perturbation theory to provide inaccurate results.

Let's consider a short portion of fiber centered at the origin of the coordinate system with a constant bent radius \(\rho\) (see Figure 1). The center of curvature is localized at the position \((\rho,0,0)\).

Figure 1. Geometry of a bent segment of a multimode fiber.

The curvature will affect propagation through three effects:

1. Geometrical effect, the bending forces guided modes to follow the curvature to keep the light inside the fiber,
2. Compression effect, the matter compression, and dilatation will modify the local refractive index,
3. Modification of the cross-section, the deformation of the fiber will modify the circular shape of the fiber.

We will neglect the latter term as the cross-section is unchanged to the first order in \(x/\rho\) [1].

Geometrical effect of bending

Due to the bending, the longitudinal wavenumber \(k_z(x)\), i.e. the projection of the wave-vector on the propagation axis \(z\), is no longer constant across the fiber cross-section. Let's call \(\beta'\) the value of the wavenumber at the center of the fiber, i.e. for \(x=0\), \(\beta'=k_z(0)\).

For the light to stay confined inside the fiber core, for any given propagating mode, one should have the optical field to be in phase at any plane orthogonal to the axis of the fiber. This gives:

\begin{equation}
k_z(x)\, \theta \, \rho(x) = cst = \beta \, \theta \, \rho \tag{1}
\end{equation}

As the curvature can be expressed as \(\rho(x) = \rho-x\), we can then express the longitudinal wavenumber as

\begin{equation}
k_z(x) = \frac{\beta'}{1-x/\rho} \approx \beta'(1+x/\rho)\tag{2}
\end{equation}

The latter approximation holds for a curvature radius large compared to the fiber core radius.

In the Helmholtz equation, we can use

\begin{equation}
k_z^2(x) \approx \beta'^2\left(1+2\frac{x}{\rho}\right)
\label{eq:ksq_bending}\tag{3}
\end{equation}

Compression effect of bending

Bending introduces an inhomogeneous deformation of the fiber core that undergoes compression for \(x>0\) and dilatation for \(x<0\) along the \(y\)-axis. Longitudinal deformation will in turn induces a transverse deformation linked by the Poisson coefficient \(\sigma\) of the material. The components of the deformation tensor reads

\begin{align}
\epsilon_{xx} &= \epsilon_{zz} = -\sigma \epsilon_{yy} \tag{4}\\
\epsilon_{yy} &= -\frac{x}{\rho}
\end{align}

The relative change of density can be written:

\begin{equation}
\frac{\Delta \eta}{\eta} = -\frac{\Delta V}{V} = - \epsilon_{yy} (1-2\sigma) \tag{5}
\end{equation}

with \(V\) the volume.

The Gladstone–Dale relation approximate \((n-1)/\eta\) to a constant:

\begin{equation}
\frac{n'(\mathbf{r})-1}{\eta+\Delta\eta} =  \frac{n(\mathbf{r})-1}{\eta} \tag{6}
\end{equation}

We then find

\begin{equation}
n'(\mathbf{r})= n(\mathbf{r})+\frac{x}{\rho}(n(\mathbf{r})-1)(1-2\sigma) \tag{7}
\end{equation}

We will use the following approximation:

\begin{equation}
n'^2(\mathbf{r}) \approx n^2(\mathbf{r})\left(1+2\frac{x}{\rho}\,\frac{n(\mathbf{r})-1}{n(\mathbf{r})}(1-2\sigma)\right) \tag{8}
\end{equation}

In the Helmholtz equation \(n'\) only appears in \(k_0^2n'^2\). As \(\beta'\) in the curved fiber is close to \(k_0 n\) to the first order, we can write:

\begin{equation}
k_0n'^2(\mathbf{r}) \approx k_0n^2(\mathbf{r})+ 2\beta' x\frac{n(\mathbf{r})-1}{\rho n(\mathbf{r})}(1-2\sigma)
\label{ksq_comp} \tag{9}
\end{equation}

We note that the geometrical effect of bending in equation \ref{eq:ksq_bending} and the effect of compression in equation \ref{ksq_comp} both contribute with a similar form but with an opposite sign. The effect of compression can be included as renormalizing the geometrical effect by a factor \(\xi\) with

\begin{equation}
\xi(\mathbf{r}) = 1-\frac{n(\mathbf{r})-1}{n(\mathbf{r})}(1-2\sigma) \tag{10}
\end{equation}

Modified Helmholtz equation

We then rewrite the scalar Helmholtz equation as

\begin{equation}
\left[\Delta_{\perp} +n^2(\mathbf{r})k_0^2-\beta'^2\left(1+\frac{2x}{\rho}\xi(\mathbf{r})\right)\right] \ket{\psi'} = 0 \tag{11}\label{eq:modified_HE}
\end{equation}

with \(\ket{\psi'}\) the field in the curved fiber.

We now have different options to solve this equation numerically:

Solution 1
Similarly to the case of the straight fiber, we can rewrite the problem as a eigenvalue problem:

\begin{gather}
\mathbf{\hat{A}'}\ket{\psi'} = \beta'^2 \ket{\psi'} \label{eq:EVP_curv} \tag{12}\\
\text{with } \mathbf{\hat{A}'} = \left(1-\frac{2x}{\rho}\xi(\mathbf{r})\right)\left[\Delta_{\perp} +n^2(\mathbf{r})k_0^2\right]
\end{gather}

The elements of the matrix \(\mathbf{A'}\) representing the operator \(\mathbf{\hat{A}'}\) projected in a discretized space basis read

\begin{equation}
\mathbf{A'}_{ij} = \left(1-\frac{2x_i}{\rho}\xi(x_i,y_i)\right)\mathbf{A}_{ij} \tag{13}
\end{equation}

with \(\mathbf{A}\) the matrix representing the operator for the straight fiber.

We can then solve this equation numerically as we did for the unperturbed fiber. This equation is very similar to the eigenvalue problem in the case of the straight fiber. The only difference is a factor \(1/\left(1+x\xi(\mathbf{r})/\rho\right)^2\) approximated by \(\left(1-2x\xi(\mathbf{r})/\rho\right)\) in equation \ref{eq:EVP_curv}. This factor does not change the number of non-zero terms in \(\mathbf{A'}\) compared to the case of the straight fiber, leading to a similar computational complexity of the system for numerical resolutions. This approach is the one adopted in the PyMMF Python module [2]. 

Note that this solution requires, for a given waveguide, to recalculate the eigenvalues and eigenvectors of a large matrix for each value of the radius. If one wants to calculate different radii for the same fiber, a new approach is required to decrease the computational time.

Solution 2

Another approach is to express the field \(\ket{\psi'}\) in the disordered waveguide as a linear combination of the straight waveguide modes and then solve an eigenvalue problem to find these coefficients, knowing the propagation constants of the straight fiber modes. As it requires to calculate the straight fiber modes first, this approach is efficient only when one wants to find the modes for different radii of the same fiber.

Let's express \(\ket{\psi'}\) in the basis of the straight waveguide modes

\begin{equation}
\ket{\psi'}=\sum_ i^M'{c_i\ket{\psi_i} }
\label{eq:decomp} \tag{14}
\end{equation}

where \(M'\) is the number of modes of the straight fiber we want to express our new modes into. 

Important note: The set of guided modes of the straight fiber is not a complete basis of the space. Rigorously, to obtain exact results, one should express the new modes as linear combinations of all the modes of the straight fiber system, including evanescent modes and non-guided propagating modes (modes that propagate both in the core and the cladding). If we do so, the current method would have no advantage compared to the previous one in terms of calculation speed. However, due to the fact that the system is weakly perturbed, modes will couple preferentially to neighboring modes. If one takes \(M' = M\), the low order modes will be correctly represented, as neighboring modes are propagating modes, but modes close to the cut-off which will overlap with non-guided modes of the straight fiber will give inaccurate propagation constants. The correct number of modes to take into account depends on the system, the correct way to proceed is to add modes until the solution converges.

By definition, the eigenmodes \(\ket{\psi_i}\), \(i \in [1,N]\) of the straight fiber and their corresponding propagation constants \(\beta_i\) satisfy

\begin{equation}
\left[\Delta_{\perp} +n^2(\mathbf{r})k_0^2-\beta_i^2\right] \ket{\psi_i} = 0
\label{eq:beta_straight} \tag{15}
\end{equation}

By injecting \ref{eq:decomp} in \ref{eq:modified_HE} and using \ref{eq:beta_straight}, we obtain

\begin{equation}
\sum_i^{M'}c_i\left[\beta_i^2-\beta'^2\left(1+\frac{2x}{\rho}\xi(\mathbf{r})\right)\right]\ket{\psi_i} = 0 \tag{16}
\end{equation}

We project the previous relation on the \(j^\text{th}\) mode \(\bra{\psi_j}\) of the straight waveguide,

\begin{equation}
\sum_ic_i\beta_i^2\braket{\psi_j}{\psi_i}-\sum_ic_i\beta'^2\braket{\psi_j}{\psi_i}-\frac{2}{\rho}\sum_ic_i\beta'^2\bra{\psi_j}x\xi(\mathbf{r})\ket{\psi_i} = 0
\quad \forall j \in [1,N]
\label{eq:int_psiprime} \tag{17}
\end{equation}

We use the orthogonality relation of the normalized guided modes

\begin{equation}
\braket{\psi_j}{\psi_i}=\iint_S \psi_j^*(x,y)\psi_i(x,y)dx dy = \delta_{ij} \quad \forall j \in [1,N] \tag{18}
\end{equation}

We simplfy the relation \ref{eq:int_psiprime}:

\begin{equation}
c_j\beta_j^2-c_j\beta'^2-\frac{2}{\rho}\beta'^2\sum_ic_i \bra{\psi_j}x\xi(\mathbf{r})\ket{\psi_i}=0 \quad \forall j \in [1,N] \tag{19}
\end{equation}

\begin{gather}
\frac{1}{\beta_j^2}\left[c_j+\frac{2}{\rho}\sum_i c_j \Gamma_{ji}\right] = \frac{1}{\beta'^2}c_j \quad \forall j \in [1,N] \label{eq:gamma} \tag{20}\\
\text{with}\quad \Gamma_{ij}=\bra{\psi_i}x\xi(\mathbf{r})\ket{\psi_j}=\iint_S \psi_j^*(x,y)\psi_i(x,y)\xi(x,y)xdx dy \tag{21}
\end{gather}

This results has been demonstrated in this form originally in [1]. Note that this relation is an eigenvalue problem for \(1/\beta'^2\) which can be directly solved.

To further simplify the resolution of the system, we use the fact that the variation of the propagation constants is small in the weakly guided approximation (\(\Delta n/n_{max} = 1-n_{min}/n_{max} \ll 1\)) as \(\beta' \in [n_{min}k_0,n_{max}k_0]\).
We can then write

\begin{gather}
\beta'=\beta_{min}+\delta\beta' \quad \delta \beta' \ll \beta_{min} \tag{22}\\
\beta_j=\beta_{min}+\delta\beta_j \quad \delta \beta_j \ll \beta_{min} \tag{23}
\end{gather}

We inject into \ref{eq:gamma}:

\begin{equation}
\frac{1}{\beta_{min}^2} \left[1-2\frac{\delta\beta_j}{\beta_{min}}\right]\left[c_j+\frac{2}{\rho}\sum_i c_j \Gamma_{ji}\right] = 
\frac{1}{\beta_{min}^2} \left[1-2\frac{\delta\beta'}{\beta_{min}}\right]c_j \tag{24}
\end{equation}

It follows, by neglecting the terms in \(\Gamma_{ij}\delta\beta_j\):

\begin{equation}
\left[\delta\beta_j-\delta\beta'\right]c_j-\frac{\beta_{min}}{\rho}\sum_i c_j \Gamma_{ji} = 0 \tag{25}
\end{equation}

Using \(\delta\beta_j-\delta\beta' = \beta_j-\beta'\), we finally obtain:

\begin{equation}
\beta_jc_j - \frac{\beta_{min}}{\rho}\sum_i c_j \Gamma_{ji} = \beta'c_j \tag{26}
\end{equation}

This relation can be written as an eigenvalue problem:

\begin{gather}
\mathbf{B}\,\mathbf{c} = \beta' \,\mathbf{c} \label{eq:EVP_B} \tag{27}\\
\text{with } \mathbf{B}_{ij} = \beta_i \delta_{ij}-\frac{\beta_{min}}{\rho}\bra{\psi_i}\xi(\mathbf{r})x\ket{\psi_j} \tag{28}
\end{gather}

The relation \ref{eq:EVP_B} can then be solved numerically for any value of the curvature. \(\mathbf{B}\) is a \(M\) by \(M\) matrix, as the number of modes \(M\) is typically much smaller than the number of points \(N^2\) in discretized transverse plane, solving equation \ref{eq:EVP_B} is much faster than solving equation \ref{eq:EVP_curv}.

Acknowledgment

I want to thank T. Cizmar and T. Tyc for their kind replies to my questions about their papers.

Bibliography

[1] M. Plöschner, T. Tyc, and T. Čižmár, Nature Photonics Seeing through chaos in multimode fibres 9, 529–535 (2015).[2] S. M. Popoff, github.com: pyMMF (2018),  zenodo.1419006.

[2] S. M. Popoff, github.com: pyMMF (2018),  zenodo.1419006.