## Finding the modes

We consider the scalar Helmholtz equation:

$$\left[\Delta_{\perp} +n^2(x,y)k_0^2-\beta^2\right] \ket{\psi} = 0$$

Which can be directly expressed as an eigenvalue problem:

$$\mathbf{A} \ket{\psi} = \beta^2 \ket{\psi} \label{eq:EVP}\\ \text{with } \mathbf{A} = \Delta_{\perp} +n^2(x,y)k_0^2$$

The solutions are the eigenvalues $$\beta_p^2$$ and the corresponding eigenvectors $$\ket{\psi_p}$$ are the transverse mode profiles.

Let's consider a discretization of the $$(x,y)$$ plane with a square mesh so that $$x(m) = m\, dh$$ and $$y(n) = n\, dh$$. For $$dh$$ small enough compared to the wavelength, the transverse part of the Laplacian in the Cartesian coordinate system reads:

\begin{gather}
\Delta_{\perp} = \partial_x^2 + \partial_y^2\\
\text{with} \left[\partial_x^2\ket{\psi}\right]_m \approx \frac{1}{dh^2}\left[\ket{\psi}_{m+1}-2\ket{\psi}_m+\ket{\psi}_{m-1}\right]\\
\text{and} \left[\partial_y^2\ket{\psi}\right]_n \approx \frac{1}{dh^2}\left[\ket{\psi}_{n+1}-2\ket{\psi}_n+\ket{\psi}_{n-1}\right]
\end{gather}

The matrix representation of the operator $$\mathbf{A}$$ for a square area of size $$L = N\, dh$$ is a $$N^2$$ by $$N^2$$ sparse matrix with only 5 non-zero diagonals. The main diagonal represents the central part of Laplacian plus the dielectric constant term:

\mathbf{A}_{ii} = -\frac{2}{dh^2} + k_0^2n^2(\mathbf{r}_i)

The unique parameter $$i \in [0,N^2]$$ indexes the transverse position. The other non-zero coefficients represent the lateral terms of the Laplacian in the $$x$$ and $$y$$ direction:

\begin{align}
\mathbf{A}_{i,i+1} &= \frac{1}{dh^2},\quad \text{for } i \not\equiv N-1 \pmod N\\
\mathbf{A}_{i+1,i} &= \frac{1}{dh^2},\quad \text{for } i \not\equiv N-1 \pmod N\\
\mathbf{A}_{i,i+N} &= \frac{1}{dh^2},\quad \text{for } i \in [0,N(N-1)]\\
\mathbf{A}_{i+N,i} &= \frac{1}{dh^2},\quad \text{for } i \in [0,N(N-1)]
\end{align}

These four contributions has exactly $$N(N-1)$$ non-zero elements, corresponding to close boundary conditions that force the points on the edge of the square window to have a missing neighbor. Close boundary conditions may create artificial unwanted reflections. However, as we are only interested here in the guided modes, which exponentially decay as a function of $$r$$ outside the fiber core, boundary conditions at the edge of the window would not affect the results as long as the window size $$L$$ is sufficiently large compared to the core size.

The matrix $$\mathbf{A}$$ can be very large but it is also very sparse and Hermitian, allowing numerical tools to efficiently find the eigenvalues and the eigenvectors with reasonable computing resources. In Python, we use the function scipy.sparse.eigh from the scipy module [3]. To select only the propagating mode we only keep the solution for which the following condition is met:

\beta > \beta_{min}

with $$\beta_{min} = k_0 n_{min}$$ and $$n_{min}$$ being the optical index in the cladding. When this condition is not satisfied, the mode propagates both in the cladding and in the core and is thus not confined in the core. As we are interested here only in the propagating modes, these modes are rejected. If we wanted to accurately account for them, more attention should be paid to the boundary conditions as their transverse mode profile would extend away from the core and would be sensitive to the edges of the observation window.

## Bibliography

[1] K. Okamoto, Fundamentals of optical waveguides, Academic press (2006).
[2] D. Marcuse, Light transmission optics, Van Nostrand Reinhold New York (1972).
[3] S. M. Popoff, github.com: pyMMF (2018),  zenodo.1419006.