Comparative modal analysis in micro–nano-optical fiber tapers using spectral parameter power series method and exact modes method

This work presents a comparative theoretical analysis of spatial modal evolution in micro/nano-optical fiber (MNF) tapers. The study proposes the use of the Spectral Parameter Power Series (SPPS) Method and compares its performance with results from the so-called Exact Modes Method (EMM) and the Finite Element Method (FEM) (the method employed by the COMSOL© software in which the computations were implemented). By using these techniques, the modal analysis and intensity evolution are discussed along different sections of the optical fiber taper. Furthermore, the data are compared considering experimental values from a real micro/nano-optical fiber taper sample. The SPPS method offers a competitive accuracy and versatility to deal with graded index profiles, its computational costs are low, and its implementation is relatively easy. The results from the SPPS method fit to those of the EM method, which sometimes involves intricated models, and those of the FEM, which may require more computational time. The SPPS method offers an average relative error of less than 5% with respect to the exact method with less computational cost compared to the FEM method for radii bigger than 2 μm at 1550 nm.


Introduction
Optical fibers which are adiabatically tapered via fused biconical process, lead to micro/nano-optical fiber structures. The interest in such structures comes from a range of applications that include probing nonlinear phenomena and atomic physics, sensing, implementing atomic mirrors and light sources with a supercontinuum spectrum, atom trapping and guiding, and optical observations of singlemolecule dynamics (see [1][2][3][4][5][6] and references therein).
One way to analyze these devices is by determining their propagated modes while computing the spatial evolution of the electric/magnetic field intensity patterns, as the taper gets thinner. This changing geometry imposes different conditions for the modes propagating through it. In the thinnest section of the optical fiber taper, the core can be assumed to vanish while it is replaced by the clad. Air or vacuum will act as the new clad, now with a refractive index contrast for which the weakly guiding condition will not hold [2]. The problem becomes more complex when multilayered fibers are considered.
The correct identification of the guided and radiated modes as the taper gets thinner allows to compute the energy that each of these modes carries. Typical analysis includes considering different methods for different sections: a weakly guiding approximation is used for the widest region of the taper; a double cladding structure for the middle section; and an exact mode eigenvalue equation for the thinnest region of the taper [7][8][9]. The study is carried out at several points along the taper to get information about the mode evolution, with certain accuracy.
Despite great advances in numerical analysis and computer sciences, the numerical methods which are often used in the solution of second-order linear differential equations face some limitations. Some of them may not be able to deal with the complex discrete spectrum of a problem. This condition may prevent the use of methods such as the shooting method, but the problem still can be solved with finite differences. However, if some problem has boundary conditions including the spectral parameter, finite differences might not be applicable at all. For some methods, if the spectral parameter takes different values, then new computations are required for each one. The Spectral Parameter Power Series (SPPS) method can be used advantageously in all those situations (see [10][11][12]). The WKB method offers an analytic asymptotic approach to solve the singular differential equation considered in this work, but it is conditioned to the spectral parameter being sufficiently large. This limitation is not present in the SPPS method [13]. Numerically, the Finite Element Method offers a solid performance when applied to solve the same problem. The SPPS method delivers a comparable precision with a shorter computational time [13,Sect. 6], [14,15]. The SPPS method both gives an analytical approach and offers a practical numerical implementation.
Its key feature used in this work is the possibility it offers to generate characteristic equations in an explicit form, turning the spectral problem into a search of zeros of the obtained approximate analytic function.
Here we propose the use of the SPPS method for the study of the different sections of the tapers. This method has shown a good performance in its application to the modal analysis of step index, multilayered and graded index optical fibers [13]. Its performance is analyzed comparing its results with those of the Exact Modes Method and with results from the COMSOL© software, based on the FEM. In contrast to this method and to the Finite Cloud Method (FCM) that allow to consider fibers with an intricate geometric structure [16], in this case the SPPS method depends on a radial symmetry. However, it can study cylindrical waveguides with almost arbitrary graded (in a wide sense) index profiles, and it can be used to compute the fundamental mode, low-order modes and high-order modes. In some tests with multimode waveguides, hundreds of modes have been correctly identified [13,Sect. 6].
The paper is organized as follows: in Sect. 2, the main aspects required to implement the SPPS Method are presented. In Sect. 3, the Exact Modes Method is introduced. In Sect. 4, the process to generate a taper is detailed, and finally in Sect. 5, the results of the computation of the modes are discussed.

SPPS Method
The SPPS method allows to solve analytically and numerically Sturm-Liouville equations. The SPPS representation has been obtained for Bessel-type singular Sturm-Liouville equations, higher-order Sturm-Liouville equations, Dirac-type systems of equations and applied for the study of different models in mathematics, physics and engineering (see [17,18], and references therein). Here we are interested in obtaining the solution u of the perturbed Bessel equation where p is a real number, p ≥ − 1/2, q is a complex valued continuous function in (0, d] satisfying a growing limit |q(x)|≤ Cx α for α > − 2, r 0 is a complex valued function, and is the (complex) spectral parameter [19]. The solution u can be associated with the electromagnetic field components of the modes propagating inside a cylindrical optical waveguide.
For the sake of simplicity, in the study of this waveguide the radius is normalized (d = 1) for the core during the construction of the mathematical model. This normalization [13,Sect. 2] improves the accuracy of the integrals that will be used in the computation of the formal powers. This in no way means a loss in the relevance or the implications of the size of the radius since its participation is relocated in the other coefficients. Moreover, inhomogeneous refractive index profiles, as for example from graded index fibers, can be studied and the SPPS method is able to distinguish different propagation constants conforming clusters [13, Sect. 6.1].
As a first requirement of the SPPS method, a non-vanishing solution u 0 of the homogeneous equation (r 0 (x) = 0) is needed. In the case that it is not available, the SPPS method allows its construction.
It has been demonstrated that if q(x) ≥ 0, x ∈ (0, d], then such solution u 0 exists and has the form [19] where the functions Y (j) are defined through the system of recursive integrals starting with Y (0) ≡ 1 and for j = 1, 2, … As a second step, the solution u 0 is used to construct a unique solution for (1) for any value of . This solution has the form [19] where the functions X (j) are defined through the system of recursive integrals with X (0) ≡ 1, and for j = 1, 2, … Here we will consider, after a shift of the spectral parameter as described in [13], the following coefficients for (1) which relate this equation to the scalar Helmholtz equation in cylindrical coordinates which models the wave propagation in a typical optical fiber: is the function describing the almost arbitrary refractive index profile, is the propagation coefficient, l is known as the order of the linearly polarized modes LP lm , with m representing the m-th solution of (7) for a single value of l, (r) represents the electric field intensity pattern of the conventional modes (see [20]), a is the radius of the fiber's core, and is the wavelength.
It is possible to fix the value of k and then find the values of for which the equation holds. The procedure consists in the solution of (1) separately for the nucleus ( x ≤ 1) and the cladding ( x > 1 ). Let us distinguish these solutions calling them u core and u clad , respectively. Both solutions must satisfy the following boundary conditions In many works considering fibers with multiple layers, only the most inner section of the fiber conforms the nucleus. This opens for the possibility of having an inhomogeneous cladding. For us, all the layers but the last one (the one we call the cladding, which can be air) are included in the nucleus. As for the mathematical considerations, the SPPS method can handle finite numbers of step discontinuities if the functions q and r are piecewise continuous on the whole segment of interest [13,Sect. 4]. So only one boundary occurs separating an inhomogeneous (multilayered or graded) core and a homogeneous cladding. However, care must be taken when adding the physical boundary conditions 2 of the problem (as the continuity of the field components or the situations that change the weakly guiding conditions). A physically meaningful solution for the cladding is well known [13] where K l is the modified Bessel function of the second kind of order l, n 2 is the (constant) refractive index of the cladding, and C is an arbitrary constant. From the boundary conditions (8), the characteristic equation of the problem can be obtained.
Substituting (9) in (10), one gets: where = a √ 2 − k 2 n 2 2 and n 2 is the (constant) refractive index of the cladding. The pairs of values ( , k ) satisfying simultaneously (11) and the propagation condition n 2 k < < n 1 k , where n 1 is the maximum value of the refractive index in the core region, correspond to the guided modes in the optical taper.

Computation of the field intensity patterns
With the use of the SPPS method, it is possible to find the values of ∈ 2 n 1 , 2 n 2 satisfying (11) for a fixed k (or vice versa). Then the solution u core (x) is found. Next, the value of C from (9) is found from (8), in the following way The solution then can be completed concatenating u core (x) and u clad (x). Finally, one can go back to the original variable: u(ax) → (r). The whole process is repeated for each considered value of l.

Exact modes method
The eigenvalue equation arising from the Exact Modes Method can be obtained considering the nano-optical fiber as consisting of just the refractive index of the cladding n 2 , followed by the refractive index of the surrounding medium (air) n 3 = 1. The propagation coefficient β, and therefore the modal parameters U and W, of the fundamental mode HE 11 are found using the exact mode equation [8,9] where J l is the Bessel function of the first kind of order l. The derivatives are (HE 1m and EH 1m ): Now, the electric field components of the fundamental mode HE 11 are expressed inside the core ( 0 < r < a ) as follows [21]: And outside the core ( a ≤ r < ∞ ) they are given by and i s the complex identity number.
In order to keep the analysis as simple as possible, in this case we are not using the full vector description of (13) the optical field as we consider that the radial component is enough to show the evolution of the optical field and mode intensity profile distribution along the taper. Such component was denoted by (r) for the SPPS method. The same computations can be done with the respective magnetic component, and with the use of these all the field components can be obtained following standard procedures (see, e.g., [22]).

Experimental setup to fabricate optical fiber tapers
The experimental setup employed to fabricate micro-and nano-optical fibers consists of two computer-controlled horizontal translation stages and a third stage to hold a ceramic micro-heater (CMH-7019) [23]. A laser diode at an operating wavelength of 1550 nm and an appropriate photodiode were used to measure the losses during the pulling process, which allowed us to control the shape of the micro-/nanofiber (MNF). The experimental setup diagram is shown in Fig. 1. The CMH has a 2-mm aperture where the optical fiber is inserted from the top. For this, the translation stage that controls the CMH first has a vertical movement; once the fiber is in the keyhole, the CMH starts to move horizontally, along the stripped fiber. The advantage of using the CMH is that we can reach controlled temperatures.

Experimental optical fiber taper fabrication results
Let us consider the use of a Single Mode Fiber SMF-28 to manufacture micro-and nano-optical fibers (MNF). The main parameters to obtain such fibers were: 46 mm of total pulling distance, 10 s of pre-heating time and temperature of 1200 °C of the CMH. Figure 2a shows the taper shape obtained experimentally, and Fig. 2b shows the minimum waist diameter close to 400 nm obtained for the total length of the fiber of 65 mm. Table 1 shows the comparison of the fundamental mode solved with different methods at 1550 nm.
Considering the micro-nano-fibers results experimentally and with the aim to explore the mode field evolution, we used three different methods. The methods SPPS, COMSOL (FEM) and exact mode solution in MATLAB were compared for the different structures of fibers as shown in Fig. 3. Figure 3 shows the evolution of the fields for different values of r, where it is clear that for values of r larger than 2 microns, both methods are closer to the exact method. For example, for Fig. 3c) with r = 2.04 mμ, the SPPS method shows an average relative error of 5.0% while COMSOL shows an average relative error of about 5.3% for a wavelength of 1550 nm. For values r < 2 mμ, the FEM tends to be closer to the exact method. Figure 4a shows results for a wavelength of 1300 nm. The results are similar to the 1550 nm case; however, for small radii, the SPPS method tends to be closer to the exact  method, for this wavelength. The average error for the SPPS method is as high as 43% in the core region, but it is better than the 98% error of the FEM method. Note that the normalization magnitude occurs at a radius near 500 nm (the core-cladding interface). This is the peak magnitude for all but the FEM's curve, in which it is located at 0 nm. Figure 4b) shows an experimental result using a BP209IR/M laser profilometer, where the field distribution in the taper can be observed. A comparison with the numerical techniques shows that they have an acceptable approximation. As a complement, Fig. 5 shows the measured field distribution in Fig. 5a), and an approximated representation of  Fig. 5b), both for a wavelength of 1300 nm.
Although both methods are a good approximation to the exact and experimental method, the SPPS method has a computational cost saving compared to COMSOL (of at least 75%) [14,15].

Conclusion
The results show that the SPPS method offers a good approximation similar to the COMSOL method for radii larger than 2 mμ at 1550 nm. For radii smaller than 2 μm, both methods deviate slightly from the exact method, with the COMSOL method showing a better approximation. However, for a wavelength of 1300 nm the SPPS method shows a better approximation for small radii. Furthermore, it can be seen by experimentation that all three methods have a good approximation.
It is difficult to say that this is an error in the computation, since the other methods also have to deal with more complicated structures than the step index variation with weakly guiding conditions. From this and from the previous results in the implementation of the SPPS method, we conclude that its application to the analysis of tapers obtained from graded index fibers (not only quadratic but with arbitrary profiles) would yield a good accuracy. The study of nonlinear effects was out of the scope of this work (powers below the lowest threshold for their occurrence can be supposed).
However, the SPPS method can be used to solve the nonlinear Schrödinger equation, from which solitons can be studied [24]. This opens the possibility to include nonlinear effects in future works. Data availability We can make it available under confidentiality for review, but it will be made available in the public domain upon publication.

Conflict of interest
The authors declare that they have no conflict of interest whatsoever.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.