BPM-Matlab: an open-source optical propagation simulation tool in MATLAB

: 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.


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][2][3][4][5][6][7][8][9][10][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.

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][26][27][28][29][30][31].The explicit FD-BPM [32] algorithm gives rise to unstable solutions such as checkerboard patterns if the step size ∆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.

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 where E p,q is the electric field at (p∆x, q∆y, z), p ∈ 0, 1, 2, . . ., N x − 1, q ∈ 0, 1, 2, . . ., N y − 1, and N x , N y are the number of pixels in the x and y directions.Integrating Eq. ( 2) in the interval [z, z + ∆z] and application of the trapezoidal rule of integration leads to the Crank-Nicolson finite difference representation of the relation between electric fields at transverse planes z and z + ∆z as where .
Representing E as a column vector, this system of linear equations can be represented in matrix format as where, neglecting boundary conditions, ⎦ are sparse block-tridiagonal matrices of order N × N, where N = N x N y .Note that in these matrices, the displaced diagonals start after N x rows or columns, and correspond to partial differentiation in the y direction.To take the Dirichlet boundary conditions into account, every N th x entry in the diagonals containing a x or −a x elements must be set to 0. The right-hand side of Eq. ( 4) is an explicit step that can be evaluated by a simple matrix multiplication whereafter E z+∆z can be solved for in a subsequent implicit step using various solvers such as the "matrix left division" operator in MATLAB.

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 B and 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 ∆z 2 step sizes each.Equation ( 2) can be rearranged as and ]︁ E p,q can be considered as an additional phase term which we will include later.So, Eq. 3 still stands true for Eq. 5, but with the modified coefficients We can now write Eq. 5 in a form that can be treated by the DG-ADI method discussed in Refs.[36][37][38][39] where the propagation through a distance ∆z is represented in two steps as with I being the identity matrix, P 1 being a tridiagonal matrix, , and P 2 and Q being block tridiagonal matrices: , where d = 1 − 2a x − 2a y .Here, P 1 and P 2 each have non-zeros only along three diagonals and Q is similar in structure to C, having non-zeros along five diagonals.The Dirichlet boundary condition can be taken into account here by modifying the matrices in the same way as before.
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 ]︁ .

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 where n(x, y) is the refractive index distribution in the absence of bending that transforms into n ′ (x, y) when the fiber is bent through a bending radius of curvature R measured from the center of the fiber located at x 0 .p e comprises the Poisson ratio and photo-elastic tensor components.
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) i z ∆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 Ti z ∆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 ∆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.

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 C matrix described above.The reason that it is not necessary to consider the whole of Eq. 4 is that the B matrix on the left-hand side describes the same propagation as C, just as an implicit step rather than an explicit step.When solving for eigenvectors, this difference is irrelevant and matrix 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 η is defined as the squared norm of the inner product integral of the present field with the field distribution of the mode [45] where E(x, y) is the electric field distribution and E l,m,p (x, y) is the electric field of the mode with azimuthal mode index l, radial mode index m, and parity p (even or odd).Equation ( 9) represents the specific mode's normalized fractional power of the total propagating E-field.The remaining power is either in other guided modes or lost through unguided, radiative modes.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 nm, the BPM-Matlab estimation of the mode field diameter of an air-clad cylindrical waveguide of refractive index 1.46 and diameter 400 nm still agreed with that of Lumerical to within 1%, while the same comparison for a waveguide of diameter 300 nm was underestimated by 6%, with the deviation further increasing for decreasing waveguide diameter.

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][49][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.

Straight multimode fiber propagation
In our first validation model, we simulate launching a Gaussian beam with an arbitrary wavelength λ = 800 nm and beam waist w 0 = 2.5 µm into a straight 1 mm long MMF of core radius a = 10 µm, n core = 1.4833, n 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 µm along the y-axis and enters the MMF at a 5 • angle to the facet normal in air, corresponding to 3.37 • 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 ∆x = ∆y = 0.067 µm and ∆z = 0.4 µ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× speedup.

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 µm, cladding refractive index n clad = 1.52, and NA = 0.1, where 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 µm for two different wavelengths, the results of which are provided in Table 1.For 1320 nm, n clad = 1.447,NA = 0.117, and V = 2.28 while for 1550 nm, n clad = 1.440,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 η 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, LP 01 and LP 03 .The fiber has n clad = 1.45, n core = 1.4666, core diameter a = 6 µ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 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 LP 01 .

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 µm to 3 µm over a distance of 10 mm.The fiber has n clad = 1.45, n core = 1.46, and the wavelength is 1000 nm.At the initial core diameter the fiber supports eight modes (LP 01 , LP 11e , LP 11o , LP 21e , LP 21o , LP 02 , LP 31e , and LP 31o ), whereas at the final core diameter, only the fundamental 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.

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][54][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.

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.

Fig. 1 .
Fig. 1.Comparison of BPM-Matlab and BeamLab simulation results.(a) Intensity and (b) phase of the electric field at the proximal end of the MMF.Electric field intensity distribution at the distal end of the MMF simulated using (c) BPM-Matlab and (d) BeamLab.

Fig. 2 .
Fig. 2. Mode intensity profiles of the five lowest order fiber modes supported by (a) a straight and (b) a bent MMF simulated using BPM-Matlab.The BPM-Matlab results are in good agreement with the mode intensities represented in Ref. [41], Fig. 4.

Fig. 3 .
Fig. 3. BPM-Matlab numerical simulations of different LP modes' normalized fractional power in a bent MMF excited with the (a) LP 01 and (b) LP 03 modes.The LP 03 can be seen to be more resilient to bending-induced mode coupling.See the text for details.

Fig. 4 .
Fig. 4. Simulation of tapering a fiber from the multimode regime to the single mode regime.The fiber tapers from a core diameter of 10 µm to 3 µm over a distance of 10 mm, followed by an untapered segment of length 5 mm.(a): The intensity in the xz-plane.The dotted lines show the core width.(b): The power within the simulation window.The dotted lines show at which z-positions each mode group is no longer guided and starts to experience loss.The z-axis of (b) also applies to (a).

Fig. 5 .
Fig.5.Simulated intensity distributions at a distance of z = 5 mm after the fiber distal end after phase calibration and applying a quadratic phase at the proximal end for MCFs.Left column: Non-twisted fiber.Right column: Twisted fiber (helical period 8 mm).Top row: Straight fiber.Bottom row: Bent fiber (RoC = 6 cm).The pattern formed by the non-twisted fiber is seen to have changed in both the positions and amplitudes of the spots when the fiber is bent.A dashed magenta circle centered at x = 0, y = 0 is added to all subfigures to clearly depict the shift in pattern in (c).