## Abstract

We present the use of the Douglas-Gunn Alternating Direction Implicit finite difference method for computationally efficient simulation of the electric field propagation through a wide variety of optical fiber geometries. The method can accommodate refractive index profiles of arbitrary shape and is implemented in a tool called BPM-Matlab. We validate BPM-Matlab by comparing it to published experimental, numerical, and theoretical data and to commercially available state-of-the-art software. It is user-friendly, fast, and is available open-source. BPM-Matlab has a broad scope of applications in modeling a variety of optical fibers for diverse fields such as imaging, communication, material processing, and remote sensing.

Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

## 1. Introduction

Tailoring light delivery through optical fibers provides the intriguing possibility to dispense with optical elements at the fiber output, opening up avenues towards miniaturized light delivery in fields such as endoscopy, communication, material processing, remote sensing, and beyond [1–11]. Endoscopy has been proven to be attractive in the context of minimally invasive *in vivo* medical imaging, helping keep morbidity and mortality rates low. The field of optical communication thrives by utilizing optical fibers for higher and faster data transmission capacity by exploiting various properties of light such as polarization, wavelength, electric field amplitude, and phase [5,12]. An accurate understanding of the propagation of a shaped beam through a given fiber is required to realise the full potential of such methods.

Optical propagation through fibers can be determined through the experimental acquisition of the Transmission Matrix (TM) of a fiber at hand [13,14], but this experimental TM is only valid as long as the fiber conformation is unchanged; any appreciable change in fiber bending or twisting may alter the TM which would have to be redetermined. Theoretical methods based on perturbation theory exist [15] to estimate the TM for a well-defined bending. The correlation between the random speckle patterns generated at the distal end of a stationary fiber and one that is bent within a range of bending angles, termed the memory effect [16], aids in dynamically redetermining the TM. Recent demonstrations of this have utilized fast binary deformable mirror devices (DMD) for real time wavefront shaping in multimode fibers [17] or digital phase conjugation for multicore fibers [18]. Nevertheless, numerical tools must supplement these theoretical methods to accurately simulate beam propagation and to calculate the memory effect, especially for non-trivial fiber structures and for cases in which perturbation theory may no longer be valid.

Multiple toolboxes such as BeamLab [19] and software packages such as Lumerical [20], based on the Beam Propagation Method (BPM) and/or the Finite-Difference Eigenmode solver and Eigenmode Expansion propagator, are widely used by the scientific community for optical fiber modeling. Although these simulation tools can model a multitude of applications in the regimes of free-space, waveguide, and non-linear optics, they are neither free as in beer (gratis) nor free as in speech (libre). Moreover, many of these tools make symmetry assumptions pertaining to the fiber/component geometry by simulating only a fraction (one half or one quarter) of the required geometry, thereby restricting complex arbitrary structure definitions. To the best of our knowledge, they are also all restricted to CPU operation and do not support GPU acceleration despite the computation times in some cases being prohibitively slow.

In this article, we present a numerical simulation tool called BPM-Matlab in which the Douglas-Gunn Alternating Direction Implicit (DG-ADI) method is used to efficiently model the electric field propagation in a wide variety of optical fiber geometries with arbitrary refractive index profiles. Our approach incorporates the important practical use cases of tapering, twisting, and bending deformations and may also offer solutions for modes for arbitrary refractive index profiles even in the presence of bending.

The key aspect of this work lies in the use of the fast DG-ADI method to model field propagation through non-symmetric refractive index profiles of arbitrary shape combined with the software’s availability for free (gratis) and as open-source (libre). As we will explain in Sec. 2., this method lends itself much better to high computational efficiency than the more common Crank-Nicolson method. Although ADI BPM was originally proposed for modeling beam propagation through waveguides by Hadley *et al.* [21], the lack of open-source software harnessing the advantages of the above method has until now restricted its wide applicability in the optics and photonics research community.

The potential impact is significant due to the code’s 1) high computational efficiency due to the DG-ADI method, 2) optional support for even faster GPU-accelerated execution, 3) support for fiber bending, twisting, and tapering, 4) inclusion of a mode solver for refractive index profiles of arbitrary shape, and 5) user-friendliness.

We give a brief description of the theoretical background governing the numerical simulations in Sec. 2., explaining the employed DG-ADI method, the incorporation of fiber bending and LP mode analysis into the framework, and the corresponding implementations in MATLAB. In Sec. 3., we discuss the applicability of BPM-Matlab and provide experimental validations for the numerical results. Finally, we conclude with comparisons to other state-of-the-art simulation software in Sec. 4. The applicability and versatility of BPM-Matlab is demonstrated with examples of single-mode, multi-mode, and multicore optical fibers, for conditions where the absence of any symmetry assumptions in geometry is necessary. The simulation tool is gratis, open-source, fast, and supports optional CUDA acceleration. It is available at https://gitlab.gbar.dtu.dk/biophotonics/BPM-Matlab.

## 2. Theory

In this section, we will present the Finite Difference Beam Propagation Method (FD-BPM) for monochromatic light, motivate the use of the Douglas-Gunn Alternating Direction Implicit (DG-ADI) finite difference method for numerical calculations, and then present its implementation in MATLAB.

The beam propagation method (BPM) solves the paraxial Helmholtz equation for the spatial distribution of the complex electric field $E(x,y,z)$ at any transverse plane along $z$ for the case when the electric field distribution at the input plane $z=0$, $E(x,y,z=0)$, and the three dimensional refractive index distribution of the optical fiber, $n(x,y,z)$ are all known. The paraxial approximation is well satisfied as long as the angles between the local $k$-vectors and the $z$-axis are small, which is the case for weakly guiding structures such as fibers with small difference in refractive index between the core and the cladding. Coupled Mode Theory (CMT) could also be used to solve Helmholtz equation for modeling beam propagation, but the solution would be based on the eigenmodes of the waveguide [22,23]. The need for knowledge of guided modes of waveguides hinders beam propagation simulations through non-trivial refractive index profiles using CMT. In this respect, BPM is advantageous since it propagates the field in the space domain without decomposing it into modes. Another popular numerical method for beam propagation is the finite-difference time-domain (FDTD) method which solves the time-dependent Maxwell’s curl equations [24]. FDTD is relatively computationally inefficient as the simulation of beam propagation in the space- and time-domains requires the entire three-dimensional computational domain to be gridded and the field data to be stored over the whole volume for each time step, while the BPM, modelling a steady-state solution for electric field propagation only in the spatial or frequency domains, only needs to store the two-dimensional field data over one transverse cross-section at a time [22]. The two main methods that could be followed to implement a BPM algorithm are the Fast Fourier Transform Beam Propagation Method (FFT-BPM) and the Finite Difference Beam Propagation Method. FD-BPM is generally favoured over FFT-BPM for simulating light propagation in integrated optical elements due to its ability to simulate waveguiding structures with discontinuities in the refractive index, variability in grid point spacing, improvement in accuracy and computational speed, and operational comfort [25–31]. The explicit FD-BPM [32] algorithm gives rise to unstable solutions such as checkerboard patterns if the step size $\Delta z$ along the longitudinal direction is not sufficiently small [29,33]. The unconditional stability of the Crank-Nicolson scheme in finite difference approximation provides more accurate solutions compared to the explicit method [29]. As we explain later in Secs. 2.1 and 2.2, the DG-ADI further surpasses the Crank-Nicolson scheme by localizing the system couplings in each step to only one transverse dimension at a time, thereby providing vastly better computational efficiency without any significant reduction in accuracy.

#### 2.1 FD-BPM

For FD-BPM [25,34], with the three-dimensional electric field and refractive index of the optical fiber denoted as $E$ and $n(x,y,z)$, the scalar paraxial Helmholtz equation can first be represented in the continuous form as

This equation is not valid for strongly guiding structures that simultaneously contain refractive indices very different to $n_0$ and guiding structures of size less than the wavelength of the light. This limitation will be quantified in Sec. 2.4.

Performing a discretization of Eq. (1) into its finite difference approximated form on an equidistant rectangular grid in $x$ and $y$ leads to

Representing $E$ as a column vector, this system of linear equations can be represented in matrix format as

where, neglecting boundary conditions,The right-hand side of Eq. (4) is an explicit step that can be evaluated by a simple matrix multiplication whereafter $\textbf{E}^{z+\Delta z}$ can be solved for in a subsequent implicit step using various solvers such as the “matrix left division" operator in MATLAB.

#### 2.2 Douglas-Gunn alternating direction implicit FD-BPM

The Crank-Nicolson model of finite difference procedure leads to Eq. (4) where there are five non-zero diagonals in each of the matrices $\textbf{B}$ and $\textbf{C}$. Despite the conceptual simplicity of the couplings that this system describes, being local in both the $x$ and $y$ directions, standard methods of solving such systems such as with the linear algebra package LAPACK are, in fact, orders of magnitude less efficient than for a simple tridiagonal system, such as representing coupling only in the $x$ or the $y$ directions [35]. However, the alternating direction implicit (ADI) finite-difference method [26] splits Eq. (3) into two tridiagonal linear equations of $\tfrac {\Delta z}{2}$ step sizes each. Equation (2) can be rearranged as

We can now write Eq. 5 in a form that can be treated by the DG-ADI method discussed in Refs. [36–39] where the propagation through a distance $\Delta z$ is represented in two steps as

The above matrices describe free-space beam propagation in a medium of uniform refractive index of $n_0$, and all matrix elements can be described by just three numbers; $a_x$, $a_y$, and $d$. If we had included in these matrices the additional phase term containing the actual spatially dependent refractive index, each matrix element might have been unique and the computation would require vastly more memory accesses and would therefore be much slower. Instead, we incorporate the additional phase term after the DG-ADI propagation steps described above and loop over all elements to calculate the final propagated electric field as $\textbf{E}^{z+\Delta z}_{prop}= \textbf{E}^{z+\Delta z}\,\frac {k_0}{2jn_0}\,\left [n_{p,q}^{2}(z)-n_0^{2}\right ]$.

#### 2.3 Implementation in MATLAB

We used MATLAB mex functions written in C for implementing the above DG-ADI method. Implementing the aforementioned computations as MATLAB mex functions enables combining the speed of subroutines written in C with the versatility and user-friendliness of MATLAB. In solving Eqs. (6) and 7, the tridiagonal matrix algorithm (Thomas algorithm) [40] was used to solve for the $E$ column vectors on the left hand sides of the equations, corresponding to a MATLAB “matrix left division".

Fiber bending is taken into account by modifying the refractive index as a function of its transverse coordinates [41] as

The effect of a taper on the fiber profile is determined by gradually scaling the features of refractive index profile $n(x,y)$ along the radial direction. Defining the taper scaling $S$ as the ratio of the fiber’s final diameter to its initial diameter and denoting the total z propagation length as $L_z$, the initially specified $x$ and $y$ positions as well as the radii $r$ of all the geometric features that make up the refractive index distribution $n(x,y)$ are multiplied by $1-(1-S)\frac {i_z \Delta z}{L_z}$ at step number $i_z$. Twisting of the fiber profile is implemented by gradually rotating the $x$ and $y$ coordinates of the geometric features an angle $T i_z \Delta z$ around the origin, where $T$ is the twisting rate of the fiber in radians per unit length of propagation. Both tapering and twisting functionalities of BPM-Matlab impose further restrictions on step size since a large $\Delta z$ could result in discretization artefacts. The size of the $z$-steps is user-specified independently for each segment of the fiber and it is up to the user to verify that the step size is small enough to obtain converged results.

There are no symmetry conditions explicitly or implicitly assumed anywhere while modeling the fiber in BPM-Matlab. Hence, the code is highly versatile providing the capability to simulate refractive index profiles of arbitrary shape and the important practical use cases of bending, twisting, and tapering deformations.

#### 2.4 Fiber modes and modal power distribution

BPM-Matlab propagates the field in the space domain without decomposing the beam into modes. However, the mode picture can be very useful and is a common way to visualise the field distribution during propagation.

The $E$-field envelope of linearly polarized (LP) modes for straight step-index fibers is described by a piece-wise combination of Bessel functions, the analytical expression of which can be found in many other works [42,43]. However, to be able to perform this visualization for more complicated fiber geometries and bending, one needs to be able to find the set of modes for refractive index profiles of arbitrary shape. BPM-Matlab includes a mode solver for that purpose, solving for the eigenvectors of the $\textbf{C}$ matrix described above. The reason that it is not necessary to consider the whole of Eq. 4 is that the $\textbf{B}$ matrix on the left-hand side describes the same propagation as $\textbf{C}$, just as an implicit step rather than an explicit step. When solving for eigenvectors, this difference is irrelevant and matrix $\textbf{C}$ will have the same eigenvectors as the equation as a whole.

The refractive index profile used here is $n'(x,y)$ as described in Eq. 8, so that bending is also taken into account. The eigenvectors are found using the built-in Matlab eigs function, which in turn employs the Krylov-Schur method [44], finding a user-specified number of modes with eigenvalues closest to 1 in the complex plane.

The LP modes’ normalized fractional power $\eta$ is defined as the squared norm of the inner product integral of the present field with the field distribution of the mode [45]

During beam propagation, BPM-Matlab can calculate and plot the normalized fractional power in any choice of set of LP modes, independent of fiber used or level of bending present. For example, in Sec. 3., we will calculate and show the coupling between straight fiber LP modes as a field distribution propagates through a fiber bent in the x and y directions.

We can also use the mode solver to study the limitations of the scalar paraxial Helmholtz equation. We calculated the mode field diameter of the fundamental mode of a small strongly guiding waveguide and compared it to the results of simulations in the commercial software Lumerical, which solves for waveguide modes in a fully vectorial, non-paraxial manner. For light with a wavelength of $1000\,\textrm {nm}$, the BPM-Matlab estimation of the mode field diameter of an air-clad cylindrical waveguide of refractive index $1.46$ and diameter $400\,\textrm {nm}$ still agreed with that of Lumerical to within 1%, while the same comparison for a waveguide of diameter $300\,\textrm {nm}$ was underestimated by 6%, with the deviation further increasing for decreasing waveguide diameter.

## 3. Validation

In this section, we present the BPM-Matlab results of simulations with various geometries such as single-mode fibers (SMFs), multimode fibers (MMFs), and multicore fibers (MCFs). These results are validated by comparing them to published experimental, numerical, and theoretical data and results of the commercially available BeamLab toolbox.

MMFs are widely studied for their applicability in endoscopic illumination [14,46,47] and signal collection [48–50]. By shaping the light appropriately into the proximal end of a MMF, one can achieve a variety of beam shapes at the distal end, even forming a focal point at a short distance in front of the fiber facet. A larger number of fiber modes provide more degrees of freedom and correspondingly better control of the output field pattern, which increase the complexity and variety of the fiber’s beam shaping abilities.

#### 3.1 Straight multimode fiber propagation

In our first validation model, we simulate launching a Gaussian beam with an arbitrary wavelength $\lambda =800\,\textrm {nm}$ and beam waist $w_0=2.5\,\mathrm{\mu} \textrm {m}$ into a straight 1 mm long MMF of core radius $a=10\,\mathrm{\mu} \textrm {m}$, $n_{\mathrm {core}}=1.4833$, $n_{\mathrm {cl}}=1.4533$, and V-number $V = 23.3$. Figure 1 shows the normalized intensity profiles at the proximal and distal ends. The input beam is offset by $5\,\mathrm{\mu} \textrm {m}$ along the y-axis and enters the MMF at a $5^{\circ }$ angle to the facet normal in air, corresponding to $3.37^{\circ }$ inside the fiber due to refraction. For this specific example, executed on a Dell Latitude E6530 laptop PC with an Intel i7-3630QM CPU, BeamLab took 15.8 minutes and BPM-Matlab took 2.2 minutes for evaluation with a resolution of $\Delta x = \Delta y = 0.067\,\mathrm{\mu} \textrm {m}$ and $\Delta z = 0.4\,\mathrm{\mu} \textrm {m}$. The resulting intensity profiles, shown in Fig. 1(c,d) agree, indicating that BPM-Matlab produces the same results as state-of-the-art software. Execution of BPM-Matlab on an Nvidia GeForce GTX 970 GPU offers a further $9\times$ speedup.

#### 3.2 Bent fiber propagation in multi- and single-mode fibers

MMFs are sensitive to bending [15,41], which is why most experiments demonstrating beam shaping through MMFs are performed with a fixed fiber. This is not transferable to applied endoscopy, and hence studying the sensitivity of fibers to bending is highly important. It has been shown that some modes are more resilient to bending than others [45], which is relevant for practical implementation of MMF beam shaping.

Figure 2 shows the simulated mode field intensity profiles of the five lowest order modes of a MMF in a straight and a bent configuration, calculated using BPM-Matlab’s mode solver. The fiber has core radius $a = 12.5\,\mathrm{\mu} \textrm {m}$, cladding refractive index $n_{\textrm {clad}} = 1.52$, and $\mathrm {NA} = 0.1$, where $NA = \sqrt {n_{\textrm {core}}^{2}-n_{\textrm {clad}}^{2}}$. The bending radius was set to 1.24 cm and the wavelength to 1064 nm. The mode profiles of the straight MMF shown in Fig. 2(a) are consistent with the profiles dictated by analytically calculated profiles based on the well known formulas involving Bessel and modified Bessel functions [42,43]. Moreover, the BPM-Matlab results for both bent and straight fibers agree with Ref. [41], Fig. 4, demonstrating that bending is modeled correctly.

Furthermore, to illustrate that bending loss is correctly estimated in BPM-Matlab, we performed a simulation study of fiber bending loss in a Corning SMF-28 single-mode fiber with core diameter $a = 4.1\,\mathrm{\mu} \textrm {m}$ for two different wavelengths, the results of which are provided in Table 1. For 1320 nm, $n_{\textrm {clad}} = 1.447$, $\mathrm {NA} = 0.117$, and V = 2.28 while for 1550 nm, $n_{\textrm {clad}} = 1.440$, $\mathrm {NA} = 0.117$, and V = 1.94. The results are in excellent agreement with the data provided in Ref. [41], Fig. 7.

Numerical and experimental works have suggested that the higher order modes of a MMF exhibit relative insensitivity towards bend-induced distortions in spite of their large mode areas [51,52]. We reproduce this finding in Fig. 3, in which we evaluate the normalized fractional power $\eta$ in each guided mode according to Eq. (9) of the total propagating electric field as a function of propagation distance and for two different injected modes, $\mathrm {LP_{01}}$ and $\mathrm {LP_{03}}$. The fiber has $n_{\textrm {clad}} = 1.45$, $n_{\textrm {core}} = 1.4666$, core diameter $a = 6\,\mathrm{\mu} \textrm {m}$ and the wavelength is 800 nm. The first 0.5 mm is straight, the next 1 mm is bent in the x direction, followed by a 1 mm bend in the y direction and finally another 0.5 mm straight section. The bending radii of curvature are 10 mm. Although the mode overlap to all 30 guided modes are calculated, only the 6 most excited modes are plotted in the figure.

Several observations can be made. In the first two segments, due to symmetry, no modes with odd parity are excited at all, but this changes once the bend in the $y$-direction starts. We can also see mode coupling as back and forth oscillations with different periodicity. Finally, we can see that the fractional power of the higher order mode $\mathrm {LP_{03}}$ of the bent MMF stays relatively constant under propagation through bends, illustrating the relative stability reported in Ref. [52] of the higher order mode against bending compared to the lower order mode $\mathrm {LP_{01}}$.

#### 3.3 Tapering a multimode fiber to a single-mode fiber

To demonstrate the tapering functionality of the toolbox, we show here a simulation in which a fiber is tapered from a core diameter of $10\,\mathrm{\mu} \textrm {m}$ to $3\,\mathrm{\mu} \textrm {m}$ over a distance of $10\,\textrm {mm}$. The fiber has $n_{\textrm {clad}} = 1.45$, $n_{\textrm {core}} = 1.46$, and the wavelength is $1000\,\textrm {nm}$. At the initial core diameter the fiber supports eight modes ($\mathrm {LP_{01}}$, $\mathrm {LP_{11e}}$, $\mathrm {LP_{11o}}$, $\mathrm {LP_{21e}}$, $\mathrm {LP_{21o}}$, $\mathrm {LP_{02}}$, $\mathrm {LP_{31e}}$, and $\mathrm {LP_{31o}}$), whereas at the final core diameter, only the fundamental $\mathrm {LP_{01}}$ mode is guided.

Using the mode solver, we solve for the above modes and excite the fiber with a superposition of them all with equal powers and random phases. Thus, one-eighth (12.5%) of the power is in the fundamental mode. In Fig. 4(a), we show the intensity in the $xz$-plane, showing the transition from a region in which we see interference between the eight injected modes to a region in which only the fundamental mode is left. Figure 4(b) shows the power within the simulation window decreasing as the higher-order modes become unguided as the fiber narrows. The dotted lines show the $z$-positions at which each group of higher-order modes is no longer guided and starts to leak out of the fiber. The final power converges to 12.5%, consistent with adiabatic projection of the initial fundamental mode into the final fundamental mode, while the power in all the other modes is lost.

#### 3.4 Multicore fiber propagation with distal end focusing

Fibers with multiple single-mode cores, referred to here as multicore fibers (MCFs), are promising for the realization of lensless endoscopes, especially in the nonlinear imaging regime, owing to their minimal modal dispersion, weak intercore coupling, and better bending resilience with respect to spatial, temporal, and polarization distortions [53–55]. Through beam shaping, using for instance a spatial light modulator, a laser beam can be divided into beamlets and coupled into the individual cores of the MCF. The amplitude and phase difference between the beamlets results in an interference pattern in the field propagating from the fiber distal end into a hypothetical sample medium, for example forming a single spot focus [56,57].

It has been shown that a MCF with a periodic core pattern, such as the traditional hexagonal pattern, gives rise to significant side lobes in the distal field image [58]. However, the aperiodicity of patterns such as the Fermat’s golden spiral core pattern [58] and pseudo-random or random patterns [54] have all been shown to reduce the side-lobes. Recent theoretical and experimental studies have shown that a conformationally invariant (bending resilient) MCF can be fabricated with twisted cores, thereby providing promising fibers for lensless endoscopes [59].

The simulated intensity profiles of the focal spot pattern made in free space near the distal end of a 41 cm long MCF demonstrated in Ref. [59] using BPM-Matlab are presented in Fig. 5 in the presence and absence of fiber bending and fiber twisting. Although non-twisted and twisted MCFs provide similar focal patterns when kept straight as seen in Fig. 5(a,b), Fig. 5(c,d) indicates that the focal pattern formed by a non-twisted MCF is sensitive to bending, whereas the one corresponding to a twisted MCF stays invariant with respect to bending. These results are congruous with the earlier observations stated in Ref. [59], demonstrating that BPM-Matlab models twisted fibers correctly.

## 4. Conclusion

We have presented the application of the Douglas-Gunn Alternating Direction Implicit method to the problem of finite-difference beam propagation in a numerical simulation tool called BPM-Matlab. The tool is capable of modeling light propagation in a wide variety of optical fiber geometries with refractive index profiles of arbitrary shape incorporating fiber bending, twisting, and tapering deformations and solve for modes for arbitrary refractive index profiles even in the presence of bending. It is provided for open-source and gratis download at https://gitlab.gbar.dtu.dk/biophotonics/BPM-Matlab. Due to the use of the DG-ADI method, BPM-Matlab provides fast results with optional CUDA acceleration on GPUs. The simulation tool is compared against state-of-the-art simulation software BeamLab and other relevant experimental and theoretical research outcomes to establish the accuracy and versatility of the code-base. Initial results based on this work can be found in Refs. [60,61].

The absence of axial symmetry assumptions in BPM-Matlab enables the user to model any arbitrary refractive index geometries. A few carefully selected cases demonstrating the versatility of the code are shown, and other geometries may easily be simulated by straightforward adaptation of the code. BPM-Matlab is beneficial in a wide variety of fields owing to its ability to model arbitrary refractive indices in both the transverse and longitudinal directions. In imaging modalities, BPM-Matlab could be used in designing experiments based on light delivery through optical fibers such as beam shaping and scanning at fiber distal end for lensless endoscopes utilizing confocal point-scanning microscopy or light-sheet microscopy, fiber-based spatially-offset optical coherence tomography (SO-OCT), etc. BPM-Matlab can also aid in the field of optical communications in designing, for instance, the optimal taper profile for minimal power loss and maximum mode conversion in a photonic lantern, or in the field of material processing by designing MCF-based laser ablation for nanoparticle synthesis, and many other experimental modalities. We expect the optics and photonics research community to adapt the open-source code-base in respect of their specific experimental parameters and thus to benefit from the versatility of our software suite.

## Funding

Teknologi og Produktion, Det Frie Forskningsråd (7017-00021); Otto Mønsteds Fond (20-81-0035); Engineering and Physical Sciences Research Council (EP/P030017/1).

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are available in Ref. [62]

(https://doi.org/10.11583/DTU.13622081)

## References

**1. **G. Zouridakis, * Biomedical Technology and Devices Handbook* (Boca Raton: CRC, 2003).

**2. **S. Mukherjee, J. S. Wysock, C. K. Ng, M. Akhtar, S. Perner, M. M. Lee, M. A. Rubin, F. R. Maxfield, W. W. Webb, and D. S. Scherr, “Human bladder cancer diagnosis using multiphoton microscopy,” Proc. SPIE **7161**, 716117 (2009). [CrossRef]

**3. **A. Ohkawa, H. Miwa, A. Namihisa, O. Kobayashi, N. Nakaniwa, T. Ohkusa, T. Ogihara, and N. Sato, “Diagnostic Performance of Light-Induced Fluorescence Endoscopy for Gastric Neoplasms,” Endoscopy **36**(6), 515–521 (2004). [CrossRef]

**4. **E. R. Andresen, S. Sivankutty, V. Tsvirkun, G. Bouwmans, and H. Rigneault, “Ultrathin endoscopes based on multicore fibers and adaptive optics: a status review and perspectives,” J. Biomed. Opt. **21**(12), 121506 (2016). [CrossRef]

**5. **P. J. Winzer, D. T. Neilson, and A. R. Chraplyvy, “Fiber-optic transmission and networking: the previous 20 and the next 20 years [invited],” Opt. Express **26**(18), 24190–24239 (2018). [CrossRef]

**6. **D. Noordegraaf, P. M. W. Skovgaard, M. D. Nielsen, and J. Bland-Hawthorn, “Efficient multi-mode to single-mode coupling in a photonic lantern,” Opt. Express **17**(3), 1988–1994 (2009). [CrossRef]

**7. **J. D. Love, W. M. Henry, W. J. Stewart, R. J. Black, S. Lacroix, and F. Gonthier, “Tapered single-mode fibres and devices. i. adiabaticity criteria,” IEE Proc.-J: Optoelectron. **138**(5), 343–354 (1991). [CrossRef]

**8. **D. J. Richardson, J. M. Fini, and L. E. Nelson, “Space-division multiplexing in optical fibres,” Nat. Photonics **7**(5), 354–362 (2013). [CrossRef]

**9. **W. R. Seitz, “Chemical sensors based on fiber optics,” Anal. Chem. **56**(1), 16A–34A (1984). [CrossRef]

**10. **M. A. Farahani and T. Gogolla, “Spontaneous raman scattering in optical fibers with modulated probe light for distributed temperature raman remote sensing,” J. Lightwave Technol. **17**(8), 1379–1391 (1999). [CrossRef]

**11. **N. Van de Giesen, S. C. Steele-Dunne, J. Jansen, O. Hoes, M. B. Hausner, S. Tyler, and J. Selker, “Double-ended calibration of fiber-optic raman spectra distributed temperature sensing data,” Sensors **12**(5), 5471–5485 (2012). [CrossRef]

**12. **A. E. Willner, H. Huang, Y. Yan, Y. Ren, N. Ahmed, G. Xie, C. Bao, L. Li, Y. Cao, Z. Zhao, J. Wang, M. P. J. Lavery, M. Tur, S. Ramachandran, A. F. Molisch, N. Ashrafi, and S. Ashrafi, “Optical communications using orbital angular momentum beams,” Adv. Opt. Photonics **7**(1), 66–106 (2015). [CrossRef]

**13. **S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media,” Phys. Rev. Lett. **104**(10), 100601 (2010). [CrossRef]

**14. **M. Plöschner, V. Kollárová, Z. Dostál, J. Nylk, T. Barton-Owen, D. E. K. Ferrier, R. Chmelík, K. Dholakia, and T. Čižmár, “Multimode fibre: Light-sheet microscopy at the tip of a needle,” Sci. Rep. **5**(1), 18050 (2015). [CrossRef]

**15. **M. Plöschner, T. Tyc, and T. Čižmár, “Seeing through chaos in multimode fibres,” Nat. Photonics **9**(8), 529–535 (2015). [CrossRef]

**16. **I. Freund, M. Rosenbluh, and S. Feng, “Memory effects in propagation of optical waves through disordered media,” Phys. Rev. Lett. **61**(20), 2328–2331 (1988). [CrossRef]

**17. **O. Tzang, E. Niv, S. Singh, S. Labouesse, G. Myatt, and R. Piestun, “Wavefront shaping in complex media with a 350 kHz modulator via a 1D-to-2D transform,” Nat. Photonics **13**(11), 788–793 (2019). [CrossRef]

**18. **N. Stasio, D. B. Conkey, C. Moser, and D. Psaltis, “Light control in a multicore fiber using the memory effect,” Opt. Express **23**(23), 30532–30544 (2015). [CrossRef]

**19. **“Beamlab,” https://www.codeseeder.com/.

**20. **“Lumerical solutions, inc.,” http://www.lumerical.com/tcad-products/mode/.

**21. **G. R. Hadley and R. E. Smith, “Full-vector waveguide modeling using an iterative finite-difference method with transparent boundary conditions,” J. Lightwave Technol. **13**(3), 465–469 (1995). [CrossRef]

**22. **W.-P. Huang, C. Xu, and B. Little, * Computer-Aided Analysis and Design in Guided-Wave Optoelectronics* (Springer, 1995), pp. 423–428.

**23. **S. E. Miller, “Coupled wave theory and waveguide applications,” The Bell Syst. Tech. J. **33**(3), 661–719 (1954). [CrossRef]

**24. **A. Taflove and K. R. Umashankar, “Review of FD-TD numerical modeling of electromagnetic wave scattering and radar cross section,” Proc. IEEE **77**(5), 682–699 (1989). [CrossRef]

**25. **Y. Chung and N. Dagli, “An assessment of finite difference beam propagation method,” IEEE J. Quantum Electron. **26**(8), 1335–1339 (1990). [CrossRef]

**26. **J. Yamauchi, T. Ando, and H. Nakano, “Propagating beam analysis by alternating-direction implicit finite-difference method,” Electron. Commun. Jpn. Part II: Electron. **75**(9), 54–62 (1992). [CrossRef]

**27. **J. Yamauchi, T. Ando, and H. Nakano, “Beam-propagation analysis of optical fibres by alternating direction implicit method,” Electron. Lett. **27**(18), 1663–1665 (1991). [CrossRef]

**28. **D. Yevick, “A guide to electric field propagation techniques for guided-wave optics,” Opt. Quantum Electron. **26**(3), S185–S197 (1994). [CrossRef]

**29. **G. L. Pedrola, * Beam Propagation Method for Design of Optical Waveguide Devices* (John Wiley & Sons, Ltd, 2015).

**30. **R. Accornero, M. Artiglia, G. Coppa, P. Vita, G. Lapenta, M. Potenza, and P. Ravetto, “Finite difference methods for the analysis of integrated optical waveguides,” Electron. Lett. **26**(23), 1959–1960 (1990). [CrossRef]

**31. **R. Scarmozzino and R. M. Osgood, “Comparison of finite-difference and fourier-transform solutions of the parabolic wave equation with emphasis on integrated-optics applications,” J. Opt. Soc. Am. A **8**(5), 724–731 (1991). [CrossRef]

**32. **Y. Chung and N. Dagli, “Analysis of z-invariant and z-variant semiconductor rib waveguides by explicit finite difference beam propagation method with nonuniform mesh configuration,” IEEE J. Quantum Electron. **27**(10), 2296–2305 (1991). [CrossRef]

**33. **H. M. Masoudi and J. M. Arnold, * Parallel-Processing Finite-Difference Beam Propagation Methods* (Springer, 1995), pp. 429–434.

**34. **D. Yevick and B. Hermansson, “Efficient beam propagation techniques,” IEEE J. Quantum Electron. **26**(1), 109–112 (1990). [CrossRef]

**35. **E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. J. Dongarra, J. Du Croz, S. Hammarling, A. Greenbaum, A. McKenney, and D. Sorensen, * LAPACK Users’ Guide (Third Ed.)* (Society for Industrial and Applied Mathematics, 1999).

**36. **J. M. McDonough, * Lectures on computational numerical analysis of partial differential equations* (University of Kentucky, 2008).

**37. **J. Douglas and H. H. Rachford, “On the numerical solution of heat conduction problems in two and three space variables,” Trans. Am. Math. Soc. **82**(2), 421 (1956). [CrossRef]

**38. **N. N. Yanenko, * The Method of Fractional Steps* (Springer, 1971).

**39. **J. Douglas and J. E. Gunn, “A general formulation of alternating direction methods,” Numer. Math. **6**(1), 428–453 (1964). [CrossRef]

**40. **B. N. Datta, * Numerical Linear Algebra and Applications, Second Edition* (Society for Industrial and Applied Mathematics, 2010), 2nd ed.

**41. **R. T. Schermer and J. H. Cole, “Improved bend loss formula verified for optical fiber by simulation and experiment,” IEEE J. Quantum Electron. **43**(10), 899–909 (2007). [CrossRef]

**42. **C.-L. Chen, * Foundations for Guided-Wave Optics* (John Wiley & Sons, Inc., 2006).

**43. **A. W. Snyder and J. D. Love, * Optical Waveguide Theory* (Springer, 1984).

**44. **G. W. Stewart, “A Krylov-Schur Algorithm for Large Eigenproblems,” SIAM J. Matrix Anal. & Appl. **23**(3), 601–614 (2002). [CrossRef]

**45. **J. Demas, L. Rishøj, and S. Ramachandran, “Free-space beam shaping for precise control and conversion of modes in optical fiber,” Opt. Express **23**(22), 28531–28545 (2015). [CrossRef]

**46. **T. Čižmár and K. Dholakia, “Shaping the light transmission through a multimode optical fibre: complex transformation analysis and applications in biophotonics,” Opt. Express **19**(20), 18871–18884 (2011). [CrossRef]

**47. **M. Plöschner and T. Čižmár, “Compact multimode fiber beam-shaping system based on gpu accelerated digital holography,” Opt. Lett. **40**(2), 197–200 (2015). [CrossRef]

**48. **S. Bianchi and R. Di Leonardo, “A multi-mode fiber probe for holographic micromanipulation and microscopy,” Lab Chip **12**(3), 635–639 (2012). [CrossRef]

**49. **T. Čižmár and K. Dholakia, “Exploiting multimode waveguides for pure fibre-based imaging,” Nat. Commun. **3**(1), 1027 (2012). [CrossRef]

**50. **Y. Choi, C. Yoon, M. Kim, T. D. Yang, C. Fang-Yen, R. R. Dasari, K. J. Lee, and W. Choi, “Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber,” Phys. Rev. Lett. **109**(20), 203901 (2012). [CrossRef]

**51. **J. M. Fini and S. Ramachandran, “Resistance of higher order modes to bend-induced mode coupling and distortion,” in Conference on Lasers and Electro-Optics/Quantum Electronics and Laser Science Conference and Photonic Applications Systems Technologies, (Optical Society of America, 2007), pp. 1–2.

**52. **J. M. Fini and S. Ramachandran, “Natural bend-distortion immunity of higher-order-mode large-mode-area fibers,” Opt. Lett. **32**(7), 748–750 (2007). [CrossRef]

**53. **Y. Kim, S. C. Warren, J. M. Stone, J. C. Knight, M. A. A. Neil, C. Paterson, C. W. Dunsby, and P. M. W. French, “Adaptive multiphoton endomicroscope incorporating a polarization-maintaining multicore optical fibre,” IEEE J. Sel. Top. Quantum Electron. **22**(3), 171–178 (2016). [CrossRef]

**54. **S. Sivankutty, V. Tsvirkun, G. Bouwmans, D. Kogan, D. Oron, E. R. Andresen, and H. Rigneault, “Extended field-of-view in a lensless endoscope using an aperiodic multicore fiber,” Opt. Lett. **41**(15), 3531–3534 (2016). [CrossRef]

**55. **S. Sivankutty, E. R. Andresen, G. Bouwmans, T. G. Brown, M. A. Alonso, and H. Rigneault, “Single-shot polarimetry imaging of multicore fiber,” Opt. Lett. **41**(9), 2105–2108 (2016). [CrossRef]

**56. **E. R. Andresen, G. Bouwmans, S. Monneret, and H. Rigneault, “Toward endoscopes with no distal optics: video-rate scanning microscopy through a fiber bundle,” Opt. Lett. **38**(5), 609–611 (2013). [CrossRef]

**57. **E. R. Andresen, G. Bouwmans, S. Monneret, and H. Rigneault, “Two-photon lensless endoscope,” Opt. Express **21**(18), 20713–20721 (2013). [CrossRef]

**58. **S. Sivankutty, V. Tsvirkun, O. Vanvincq, G. Bouwmans, E. R. Andresen, and H. Rigneault, “Nonlinear imaging through a fermat’s golden spiral multicore fiber,” Opt. Lett. **43**(15), 3638–3641 (2018). [CrossRef]

**59. **V. Tsvirkun, S. Sivankutty, K. Baudelle, R. Habert, G. Bouwmans, O. Vanvincq, E. R. Andresen, and H. Rigneault, “Flexible lensless endoscope with a conformationally invariant multi-core fiber,” Optica **6**(9), 1185–1189 (2019). [CrossRef]

**60. **M. Veettikazhy, A. K. Hansen, D. Marti, K. Dholakia, and P. E. Andersen, “Numerical comparison of robustness of shaped beam delivery through multimode and multicore fibre against fibre bending,” in * Adaptive Optics and Wavefront Control for Biological Systems VI*, vol. 11248T. G. Bifano, S. Gigan, and N. Ji, eds., International Society for Optics and Photonics (SPIE, 2020), pp. 43–46.

**61. **M. Veettikazhy, A. K. Hansen, D. Marti, K. Dholakia, and P. E. Andersen, “Numerical comparison of robustness of multimode and multicore fibre sensitivity against fibre bending,” in 2019 Conference on Lasers and Electro-Optics Europe European Quantum Electronics Conference (CLEO/Europe-EQEC), (2019), pp. 1.

**62. **M. Veettikazhy, A. K. Hansen, D. Marti, S. M. Jensen, A. L. Borre, E. R. Andresen, K. Dholakia, and P. E. Andersen, “Data underpinning: BPM-Matlab - An open-source optical propagation simulation tool in MATLAB,” (2021).