Electromagnetic mode matching in a Wilson basis: optical fiber connections with a gap

On the basis of two optical fibers with an optional lateral and longitudinal displacement in a homogeneous background medium, we describe a general full-vectorial Wilson-basis discretized mode-matching method that evaluates the converged electromagnetic fields after all resonances in the possible gap cavity have settled. Wilson basis functions feature strong localization in both the spatial and the spectral domain, which allows for efficient modal electromagnetic field expansions, adequate truncation of field propagation operators, and sparse translation operators, which in turn allow to make ad hoc electromagnetic field-matching operators. For physical contact connections between single-mode fibers with a mode-field diameter mismatch, we obtain attenuation curves that are right between those obtained from approximation methods that either effectively match the electric field or the magnetic field under the assumption of a vanishing reflection. For fibers separated by a growing gap, constructive and destructive interference patterns in the cavity are computed by the successive application of Love’s equivalence principle and the propagation operator. By leveraging the physical width of the Wilson basis functions and the stepsize of the propagation operator, the initial operators may be reused in solving the interface problem for other wavelengths. For multi-mode fiber connections, we provide attenuation curves on a modal electromagnetic field level, as well as for overfilled and core-confined target encircled-flux compliant launches. A comparison to geometrical-optics based approaches shows attenuation differences in the order of several hundredth of a dB, and although that is small, it is significant for modern connection attenuation specifications. In the final example of connections between regular and trench-assisted multi-mode fibers, we notice that the relative change in the cumulative near-field power distribution can be significant, despite a marginally small attenuation. The influence of a core diameter and/or numerical aperture mismatch can be examined with a deliberate lateral misalignment.


Introduction
Detachable low-attenuation optical fiber connections are indispensable to configure optical networks in a flexible manner. In optical systems that employ the space-division multiplexing technique to further increase the bandwidth of a single fiber (Richardson 2010), sub-micrometer fiber alignment is needed to preserve the specific linear combinations of the modal electromagnetic fields that span spatio-multiple transmission channels (Nagase 2014). In structured indoor cabling systems that rely on cascaded multi-mode fiber technology with a 50 m core diameter, every connection should make and maintain physical contact between the accurately aligned fibers in order to keep the penalty on the power budget below several tenths of a dB. A further reduction of the connection attenuation would allow for more inline connections without sacrificing the total link length, and without exceeding the tight optical power budgets that are just shy of 2 dB for modern IEEE802.3 Ethernet networks (http://www.ieee802.org/3 2010) and will be tightened even further as data rates increase (Spurgeon 2000). Other (potentially) detachable optical interfaces that bring fiber optics closer to the transceivers or compute engines include free-space optics lenses inside a transmitter and receiver optical subassembly (Lam 2011), 45 • mirrors that are integrated in board level optics (Fettweis et al. 2019), or direct butt-coupled fibers to rectangular waveguides on chip (Hauffe et al. 2001).
The examples mentioned above involve electromagnetic connections across planar interfaces. However, the actual material configurations in the half spaces on either side of those interfaces may be slightly, or very different. Even if the electromagnetic fields on either side of an interface may efficiently be described in terms of modal expansions, those modal expansions are generally mutually unrelated. Every time one would encounter such electromagnetic interface problems, one would have to expand one set of modes (including the radiation field) in terms of the other: a tedious task indeed. This can be circumvented by expressing every modal expansion in terms of the same complete set of orthogonal basis functions, i.e., a universal interface. We have discussed such a universal interface in terms of Wilson basis functions in (Floris and de Hon 2018a, b). The Wilson basis functions are duo-localized, i.e., they decay exponentially both in the spatial and in the spectral domain (Wilson 1987;Daubechies et al. 1991). In the Wilson basis, high-frequency electromagnetic wavefields can be represented through a parsimonious multidimensional set of coefficients (Arnold 2002;Floris and de Hon 2018a). In (Floris and de Hon 2018a, b), we introduced the Wilson basis universal interface approach, applied it to the electromagnetic scattering of incident modes at the end-face of a multi-mode fiber, and reciprocally, the coupling of light into a fiber. The major advantage of our Wilson-basis approach from a computational mode-matching perspective, is that it becomes irrelevant whether the coefficients represent fields that are associated with a round fiber, a rectangular dielectric waveguide or any other physical medium for which a complete set of (modal) transverse electromagnetic field solutions can be computed.
Below, we argue that the full benefits of the universal interface approach for computational electromagnetics really become manifest when two (or more) structures are connected along their planar interfaces. The main contribution of this paper is to demonstrate this through examples of two, not necessarily aligned nor identical, optical fibers that may be longitudinally separated by a homogeneous medium. In particular, we have devised a Wilson-basis discretized mode-matching method that relies on a recurrence scheme to accelerate the convergence of the electromagnetic fields that resonate in the cavity formed by two successive fiber interfaces with a possibly vanishing gap in between.

Page 3 of 36 730
The effective duo-localization of the Wilson basis ensures the rapid convergence of the spectral Green's function integrals, needed for the evaluation of the electromagnetic wavefields that emanate from Wilson-basis discretized electromagnetic source distributions in a homogeneous medium. From the Wilson-basis coefficients that describe a given transverse electromagnetic field at the end-face of the physical medium, it is possible to directly configure two electromagnetic source distributions that produce one-way forward and one-way backward propagating electromagnetic fields in a homogeneous medium in such a way that boundary conditions are satisfied at the interface (Floris and de Hon 2018b). The underlying Green's function integrals constitute a propagation operator that synthesizes the electromagnetic field at a chosen distance from the source plane. The successive application of Love's equivalence principle and the propagation operator allows for the evaluation of the electromagnetic fields at integer multiples of the selected distance almost effortlessly . It is convenient that the coupling matrix between Wilson basis functions and laterally displaced Wilson basis functions is sparse (Floris and de Hon 2018a), so that any adjustment of the lateral alignment of any of the the physical media (for instance a fiber) with respect to the interface is easily evaluated in the Wilson basis prior to executing the mode-matching procedure.
Our first aim is to demonstrate the effectiveness of the Wilson-basis universal interface computational framework in solving multiple interface problems. Our second aim is to discuss several interesting realistic fiber connection applications, outlined in Sect. 2. The main ingredients of our Wilson-basis interface framework: the Wilson basis, the Wilson-basis expansion of electromagnetic fields, and the construction of Wilson-basis discretized oneway propagating fields in homogeneous space are summarized in Sect. 3.

Optical fiber interconnect problems
In this manuscript, we aim to serve a diverse readership. If your primary interest is computational electromagnetic modeling of realistic interface problems, then you can find the theoretical underpinning in Sects. 4-6. Conversely, if you are mainly looking for a detailed discussion of modern optical fiber interconnect problems, you can find extensive practical simulation results in Sects. 7-11. Below, we present a roadmap to these sections to help navigate the paper.
In our Wilson-basis interface framework, all reflection-transmission problems that involve excitations and reflections that are spanned by the guided modal electromagnetic fields of a single optical fiber may be solved through the procedure described in Sect. 4. To cope with all other excitations, and to avoid having to discretize the continuous spectrum that constitutes the radiation field, we exploit the small transverse contrast in the refractive index profile to temporarily replace the fiber by a homogeneous background medium and solve the reflection-transmission radiation-field problem for two homogeneous half-spaces in Sect. 5.
In Sect. 6, we consider a cavity formed by two fiber interfaces separated by a slab of homogeneous material (e.g. air), and describe a loop operator that evaluates the evolution of electromagnetic fields that make a loop inside the cavity, i.e., that arrive back at the first interface after two successive reflections of respectively the first and the second slab-tofiber interface. Particularly for configurations with highly reflective interfaces, we would like to avoid having to evaluate a Neumannn series through a direct term-by-term application of the loop operator in order to obtain the converged electromagnetic fields in the cavity, together with the reflection in the first fiber, transmission in the second fiber and the radiation fields. Instead, we describe an iterative method that decomposes the looped electromagnetic field in terms of the field in the previous loop and a residual field. The underlying Neumann-series for the converged electromagnetic fields reduces to evaluating geometric progressions of known vector fields, and evaluating the loop operator for a sequence of quickly decaying residual fields.
To test the approach, we consider single-mode fiber connections with a vanishingly small gap in Sect. 7 and compare the connection attenuation to the approximated case that spans the incident field in the first fiber by the guided modes of the second fiber by using projection (overlap) integrals under the assumption that the reflected field vanishes (Neumann 1988). There are two possible choices for the projection operation, one effectively matching the electric field, the other the magnetic field. As the mode-field diameter mismatch between the two fibers is increased, the effective impedance mismatch leads to two attenuation curves for the two choices of the approximate case, and interestingly we find that the curve produced by our mode-matching approach is right in the middle of the two. For the relatively small mode-field diameter mismatches that we considered here, the differences between the full mode-matching method and the approximation are only a few hundredths of a dB. However, in case of dissimilar waveguides or a significant separation, the assumption of a vanishing reflected field is no longer valid. We demonstrate examples of this for a fixed lateral displacement and a range of longitudinal separations in Sect. 8.
In Sect. 9, we apply the full mode-matching method to physical contact connections between two multi-mode fibers and give the attenuation curves for lateral misalignments up to 6 m for each of the 342 guided modes of a regular multi-mode fiber. In order to make judicious, reproducible and repeatable attenuation measurements in practice, international standards organizations such as IEC (IEC 61280-4-1 2021) and TIA (TIA-455-203-A 0000) demand an encircled-flux compliant launch, which means that the cumulative near-field power distribution measured at the fiber end-face of the plug under test shall keep within tight bounds that are imposed on a handful of target levels. Even though there are commercially available mode conditioner devices that aid in achieving an encircledflux compliant launch in practice, from a theoretical perspective the five prescribed targets are hardly restrictive in view of the 342 guided modal electromagnetic fields that one can mathematically choose from. To be able to study the impact of fiber alignment and fiber geometry parameters such as the fiber core diameter and numerical aperture on the attenuation, we follow a mode-continuum approximation (Daido et al. 1979) to narrow down to a sensible combination that employs all modes in such a way that the modal amplitudes are sampled from a smooth modal power distribution function.
In an earlier publication, we have shown that it is also possible to generate encircledflux compliant ray-based launches in multi-mode fibers (Floris et al. 2019), and to determine design curves that balance core diameter mismatch, numerical aperture mismatch and lateral misalignment for reference-grade connectors specifications that achieve less than a tenth of a dB of attenuation when it makes a physical contact connection to another reference-grade connector (Floris et al. 2012). Here, we make a comparison between the vectorial full-wave mode-matching technique and geometrical optics, which reveals that both approaches are in agreement for lateral misalignments up to about 4 m , after which the attenuation curve produced by the geometrical optics method is a bit higher.
We evaluated the attenuation curves for each of the speckle patterns that together produce a spatio-temporal target encircled-flux compliant near-field pattern. The variation in the attenuation is a measure for modal noise, and in (Floris et al. 2019), we have shown that this coherent wavefield phenomenon can alternatively be studied in terms of incoherent rays with the aid of the cumulative near-field patterns of the individual speckle patterns. This is also discussed in Sect. 9.
In Sect. 10, we examine the attenuation and return loss for two longitudinally separated multi-mode fibers and a lateral offset. Such connections were also studied with ray-based methods (Di Vita and Rossi 1978;van Etten et al. 1985;Bude and Van Etten 1987) albeit for very large separations and for an overfilled launch. Here, we consider every mode as an individual excitation. Via weighted combinations of modal excitations, we construct realistic encircled-flux compliant launches. This approach enables us to study the response of the diverging fields in the cavity on the attenuation in the receiving fiber for the whole transition from small to very large separations.
We conclude in Sect. 12 with examples of connections between regular and trenchassisted multi-mode fibers. Even though additional core diameter and numerical aperture mismatches might not lead to significant changes in the attenuation, it does show up as a relative change in the cumulative near-field pattern in the sense that either the entire curve shifts upwards (narrower distribution) or downwards (wider distribution). The observation that a deliberate misalignment alters the shape of the curve suggests that relative encircledflux measurements can reveal very subtle differences in the refractive index profile across the fiber connection interface. It may be possible to exploit this sensitivity for profile measurement purposes.

Decomposition of fiber modes to one-way propagating fields in homogeneous space
The Wilson basis consists of basis functions w n (x) that are trigonometric modulated and translated versions of a single window function (x) . The spatial translation is denoted by index n and the modulation by index . We adopted the algorithm of Daubechies et al. (Daubechies et al. 1991) for the construction of (x) and gave a streamlined summary in (Floris and de Hon 2018a). With the aid of the forward Fourier transformation the spectral counterpart of the Wilson basis functions are described by The spectral basis functions have two peaks due to the trigonometric modulations of the Wilson basis functions. Further, the function ̂= F T ( ) decays exponentially for a Gaussian window function (x) (Floris and de Hon 2018a). (1) We proceed by giving a brief summary on the construction of one-way propagating fields in homogeneous space in the Wilson basis from (Floris and de Hon 2018b). In particular, we emphasize the construction of specific one-way forward propagating field ̄ + (blue wiggles) and one-way backward propagating field ̄ − (red wiggles) that match to the waveguide-specific guided mode ̄ (green wiggle) as shown schematically in Fig. 1.
Let ( t , z) = E x , E y , H x , H y T describe the transverse electromagnetic field components of an arbitrary wavefield along the transverse coordinate t = x x + y y at a longitudinal coordinate z. In the remainder of this section, we assume that represents an arbitrary guided mode in the fiber. We expand the vectorfield in the Wilson basis by writing where ̄ (z) is a four-vector containing Wilson basis coefficients and w ( t ) = w x n x (x)w y n y (y) is a 2D Wilson basis function. The multi-index = x , n x , y , n y is used to compress notation (and we shall also use the multi-indices and for this purpose), and the parameter d scales the Wilson basis function to the size of the modal electromagnetic fields. The Wilson basis coefficients ̄ (z) are computed numerically using fixed-point cubatures as detailed in (Floris and de Hon 2018a).
We proceed with the construction of a one-way forward propagating field ̄ + (z) and a one-way backward propagating field ̄ − (z) , in such a way that the continuity of the transverse electromagnetic field ̄ (z 0 ) =̄ + (z 0 ) +̄ − (z 0 ) is satisfied across the interface z = z 0 . The field propagation in the positive and negative z direction are relative to the assumed exp(j t) dependence on time t, where denotes the angular frequency.
In the first step, we invoke Love's equivalence principle (Love 1901) and construct equivalent magnetic sources for the electric field of ̄ (z) at z = z 0 . This leads to a source distribution = J x , J y , K x , K y T in the Wilson basis where the superscript K indicates that only magnetic sources are present.
In the next step, we configure a second electric source distribution by invoking Love's equivalence principle again, but now on the magnetic field that was generated by the sources ̄ K in Eq. (4), i.e., Fig. 1 For a fiber mode ̄ incident on a fiber end-face (single green wiggle in the fiber), a decomposition into a one-way forward ̄ + (blue wiggles) and one-way backward ̄ − (red wiggles) propagating fields in the homogeneous medium satisfy the boundary conditions. (Color figure online) where the superscript J, K emphasizes that ̄ contains a distribution of electric equivalent sources of the (magnetic) field that in turn was generated by magnetic sources. The source distribution ̄ J,K in Eq. (5) generates the same electromagnetic field as the source distribution ̄ K in Eq. (4) for z ≥ z 0 . This becomes apparent upon inspection of the four-by-four matrix ̄ , ( z) with z = z − z 0 , which evaluates the double integral (Floris and de Hon 2018b) In Eq. (6), ̃ T ( t , z) = 1 2 e −jk z | z|M is the spectral dyadic Green's function with, in which the impedance Z = √ 0 ∕ = √ 0 ∕( 0 n 2 ) relates to the refractive index n, permittivity 0 and permeability 0 . In view of the branch cut associated with the square root and with k t = | t | , we define Im(k z ) < 0 , and Re(k z ) ≥ 0 if Im(k z ) = 0 , associated with the physical Riemann surface. Furthermore, M in Eq. (7) depends on For z → 0 , the magnetic ̄ K sources in Eq. (4) generate an electric field that changes sign across the z = z 0 interface, whereas the electric ̄ J,K sources generate a magnetic field that changes sign across z = z 0 . Hence, the electromagnetic field with the + subscript is one-way propagating in the positive z-direction and matches the electric field of ̄ at the positive side of the interface at z ↓ z 0 . On the opposite side of the interface for z < z 0 , the electromagnetic fields vanish due to destructive interference of the fields generated by the sum of the sources ̄ K and ̄ J,K . The − subscript denotes a one-way propagating electromagnetic field in the negative z-direction for z < z 0 that matches the electric field of ̄ at the interface z ↑ z 0 , while the field vanishes for z > z 0 .
In a similar fashion, one can deploy equivalent electric sources supplemented with magnetic sources to achieve a one-way propagating electromagnetic field that has the same transverse magnetic field as ̄ at the interface. Upon combining Eqs (9) and (12), the fiber mode ̄ is matched to two one-way propagating fields ̄ + and ̄ − in homogeneous space by choosing (Floris and de Hon 2018b) The sum ̄ + (z) +̄ − (z) at z = z 0 matches ̄ , because the sum of the fields ̄ E+ (z) +̄ E− (z) describes the desired electric field of ̄ (albeit with a factor 2), and has a zero magnetic field, whereas the sum of fields ̄ H+ (z) +̄ H− (z) describes the desired magnetic field of ̄ (albeit with a factor 2), and has a zero electric field.
To evaluate these fields at an arbitrary longitudinal z coordinate, one only needs to evaluate the matrix ̄ , (z − z 0 ) for a compact set of source Wilson source distributions (Floris and de Hon 2018b). As an added bonus, the fields may be evaluated at subsequent longitudinal coordinates z = m z , with m = 1, … , M almost effortlessly by the repetitive application of Love's equivalence principle and re-use of the propagation operation ̄ , ( z) as shown in . Moreover, by ensuring that the quantities d∕ and ( z)∕d are left unchanged, one can re-use the pre-computed coefficients ̄ , ( z) for other wavelengths too.
We proceed with a brief notation-consistent summary of the solution strategy to reflection-transmission problems at the end-face of an optical fiber from (Floris and de Hon 2018b).

Reflection-transmission problems for fiber end-faces
For fiber end-face reflection-transmission problems as sketched in Fig. 2, we now apply the field decomposition method of Sect. 3. For the forward propagating fiber mode ̄ + ,m that is incident from the left on the interface at z = z 0 in the top sketch of Fig. 2, the decomposition results in a one-way forward field ̄ + ,m and a one-way backward field ̄ − ,m in the homogeneous medium. For brevity, we denote the mode number m for the combined radial and angular indices of the pertaining fiber mode.
In a similar fashion for the backward propagating m th fiber mode ̄ − ,m in the fiber (second sketch in Fig. 2), the one-way forward field ̄ + ,m and one-way backward field ̄ − ,m in the homogeneous medium satisfy the boundary conditions. As we use the homogeneous-space wavefields ̄ + ,m , ̄ − ,m , ̄ + ,m , and ̄ − ,m as building blocks, we deliberately chose to assign two separate symbols ̄ and ̄ to emphasize that the associated fiber modes are propagating forward and backward in the fiber, respectively.

Consider an excitation
,m written as a linear combination of forward propagating guided fiber modes ̄ + ,m weighted with coefficients c I m (green wiggle in Fig. 3). In the homogeneous medium, the transmitted field ̄ T,guided = ,m is associated with the forward-propagating fiber modes (blue wiggles).
We aim to cancel the unphysical one-way backward field ∑ m c I m̄ − ,m from the homogeneous medium by expanding it in terms of one-way backward fields ̄ − ,m , for which we already determined the scattered fields (viz. second schematic in Fig. 2). This gives rise to a reflection in the fiber ̄ R,guided = ∑ m c R m̄ − ,m (black wiggles) and the second contribution to the transmission ̄ −R,guided = ∑ m c R m̄ + ,m (dark red wiggles). All that remains is a (typically small (Floris and de Hon 2018b) incident from the right, which in turn excites the fields ̄ R,rad in the fiber and ̄ +,rad in the homogeneous medium (purple wiggles) that are associated with the radiation spectrum of the fiber. The solution strategy to the evaluation of these residual fields is discussed in Sect. 5. The vector R that contains the coefficients c R m of the guided reflected modal electromagnetic fields is obtained by solving the system of equations

Fig. 2
Top: a one-way forward propagating fiber mode ̄ + ,m is decomposed into a one-way forward propagating field ̄ + ,m and a one-way backward propagating field ̄ − ,m in homogeneous space. Bottom: a backward propagating fiber mode ̄ − ,m is decomposed into a one-way forward propagating field ̄ + ,m and a one-way backward propagating field ̄ − ,m in the homogeneous medium. Then, the residual field ,m excites the radiation fields in the fiber and the homogeneous medium.
Previously (Floris and de Hon 2018b), we considered situations where the excitation of the residual fields in the radiation spectrum of the fiber was negligibly small. However, for large excitations from the homogeneous medium, for example due to multiply scattered resonating fields in a gap between two fibers, the power carried by the residual fields can be substantial.
To avoid the burden of discretizing the radiation spectrum of the fiber modes just for the computation of the residual fields, we exploit the small contrast in the (transverse) refractive index profile by replacing the fiber by an effective homogeneous background medium, and solve the residual-field reflection-transmission for this effective medium. We adopt this approach below to evaluate the reflected field in the Wilson basis for an interface between two homogeneous half-spaces.

Reflection-transmission problems of two homogeneous half-spaces
For an interface between two homogeneous, isotropic, lossless, dielectric media, the spectral representation of the reflection is directly related to the incident electromagnetic field through Fresnel reflection coefficients for the s-polarized and p-polarized field components. The well-established approach has recently allowed the identification of conditions for incident evanescent plane-wave excitation to experience total internal reflection (Lakhtakia and Mackay 2008).
The incident Wilson-basis discretized transverse electromagnetic field formulation in Eq.
(3) has a spectral representation with the four-vector ̄ I = [̄ I,E ,̄ I,H ] T containing two transverse electric Wilson-basis coefficients ̄ I,E and two transverse magnetic Wilson-basis coefficients ̄ I,H oriented along Cartesian unit vectors. Without loss of generality, we assume that the interface between two homogeneous media is located at z = 0.
To distinguish between the s-and p-polarisation in Eq. (15), we introduce polar coordinates k x = k t cos and k y = k t sin , and rewrite the unit vectors x and y as where k t = t ∕|k t | = cos x + sin y is the unit vector in the t direction. By substituting Eq. (16) in Eq. (15), we recognize the s-and p-polarizations of the electric field with amplitudes a s, and a p, in the Wilson basis Here, the numbers 1 and 2 between the square brackets denote the Wilson-basis discretized Cartesian E x and E y electric field components. In a similar fashion, we write the transverse magnetic field as with Wilson basis coefficients Here, the numbers 3 and 4 between the square brackets denote the Wilson-basis discretized Cartesian H x and H y magnetic field components.
The spectral representation of the reflected field is similar to Eq. (15) and reads (Lakhtakia and Mackay 2008) and propagates in the negative z-direction. The reflection operator R produces a four-vector with components and applies the Fresnel reflection coefficients that depend on t . In view of the branch cut associated with the square root The wavenumber k 1 of the incident and reflected fields in the first medium n 1 are related to the free-space The transmitted field in the second medium n 2 has wavenumber k 2 = k 0 n 2 .
By applying the transformation in Eq. (16) to Eq. (22), and inserting the amplitudes in Eqs (18) and (20), the reflection operator R ̄ I produces a spectral field representation with components oriented along Cartesian unit vectors.
We wish to evaluate the integral in Eq. (21) and expand the field R ( t , 0) again in the Wilson basis through the inner-product By inserting Eq. (21) into Eq. (25) and changing the order of integration, we arrive at the spectral integral representation for the Wilson-basis coefficients for the reflected field, By bringing the coefficients ̄ I outside the integral, we recognize from Eq. (24) that Eq. (26) can be written as a linear combination of six spectral integrals with m ∈ (1, … , 6) , and where the reflections R m are given by, In (Floris and de Hon 2018b), we used fixed-point cubatures on a cylindrical polar grid of abscissae to evaluate the spectral Green's function integrals in Eq. (6). Because the integrals in Eq. (27) have a similar form, we re-use the same approach to evaluate them. In particular, we make use of one, two or four non-overlapping cubatures, each cubature centered about a different peak of the Wilson basis source function w in the k x k y -plane. For all Wilson-basis test functions w under consideration in Eq. (21), we re-use the same cubature grids.
With the aid of the transformations k x = k t cos( ) and k y = k t sin( ) with the Jacobian d t = k t dk t d , the k −2 t singularities in Eq. (28) are conveniently removed. In view of the branch cuts at k t = k 1 and k t = k 2 , we subdivide the integration over k t in three regions. Without loss of generality, we assume n 1 < n 2 . In the opposite case, only the coefficients r p and r s in Eq. (23) change sign as n 1 → n 2 and n 2 → n 1 . In the first region with k t ∈ 0, k 1 , we deployed 300 abscissae along the radial coordinate, and 1200 along the angular coordinate and applied the substitution u = √ k 2 1 − k 2 t to aid the convergence of the six integrals. In the second region with k t ∈ k 1 , k 2 , we also deployed 300 radial and 1200 angular abscissae, however, we did not apply any substitutions. In the third region with k t ∈ k 2 , ∞ , we applied the substitution u = √ k 2 t − k 2 2 and again deployed 300 radial and 1200 angular abscissae.
Fortunately, we only need to evaluate the integrals for the non-displaced Wilsonbasis source functions w , with (l, n) = (0, 0) and n ∈ (0, 1) for 1 ≤ ≤ l max applied to the pairs ( x , n x ) and ( y , n y ) (Floris and de Hon 2018b). The integrals with spatially displaced Wilson-basis source functions w + with = (0, s x , 0, s y ) are the same as those for the non-displaced cases when the test functions are displaced according to w →w − (Floris and de Hon 2018b).
Consider an interface between air with n 1 = 1 and a glass material with n 2 = 1.452 that resembles a fiber cladding at a wavelength of 850 nm . To be able to adequately truncate the set of test function w for the expansion of an arbitrary reflected field ̄ R in Eq. (26), we visualize the behavior of the absolute maximum of the integrals I m of (28) Eq. (27) in Fig. 4 for all source functions w with ≤ 4 and all test functions w with ≤ 6. As the distance | x | + | y | between the Wilson source function w and the test function w increases, the coupling decreases quite rapidly. The behavior for the integrals I m associated with m ∈ {2, 5} (red curves with • markers) decay a bit differently due to the product k x k y in the integrand in Eq. (28), compared to the integrals with m ∈ {1, 3, 4, 6} that either depend on k 2 x or k 2 y (black curves with + markers). For an increasing spatial separation | n x | + | n y | , the same integrals seem to decay very slowly as shown in Fig. 5, which indicates that the reflected field due to a confined incident field can potentially be very wide.
Although it may appear that an a priori truncation of | n x | + | n y | ≤ 12 is inadequate, we want to emphasize that the actual reflected field depends on the specific linear combination of these integrals. For electromagnetic wavefields that resonate in a gap between two fibers, we estimate that we achieve an accuracy up to 5 digits judging from the power balance of the transmitted and reflected fields using the aforementioned truncations of the Wilson basis test functions.
We proceed in Sect. 6 with a procedure to evaluate the resonances in a gap between two fibers.

Fig. 5
For an air-glass interface, the absolute maximum of the integrals in I m in Eq. (27) decay relatively slowly as the spatial distance between | n x | + | n y | between the Wilson source and test functions increases Fig. 6 Sketch of a reflection-transmission problem comprising two laterally and axially displaced optical fibers in a homogeneous medium and an excitation by an arbitrary combination of guided fiber modes (green wiggle in the first fiber). The triple wiggles denote electromagnetic fields that are readily evaluated through one-way field decompositions and propagation (Sect. 3). The small purple wiggles represent approximated radiating fiber modes (Sect. 5). We wish to minimize the mismatch across the dashed interface at z + 1 = z 1 + , with ↓ 0 , solving the multiple-scattering problem Page 15 of 36 730

A resonating cavity between two fibers
Consider two optical fibers with flat end-faces that are displaced both laterally and axially with respect to each other in a homogeneous medium (see Fig. 6). The Wilson-basis discretized incident field ̄ I (green wiggle) in the first fiber propagates to the first interface at z = z 1 . For readability, we henceforth suppress the multi-index in the subscript of all Wilson-basis fields. Associated with the incident field ̄ I , the first contribution to the reflection ̄ R,guided 0 (red wiggles) and the transmission ̄ 1+ 0 (blue wiggles) in the homogeneous medium satisfy the boundary conditions at the first interface z = z 1 as detailed in Sect. 4. The subscript 0 indicates that the reflection does not include contributions due to multiply scattered fields resonating in the gap. For convenience, we denote the propagation of a Wilson-basis discretized vector field ̄ over a gap distance z = z 2 − z 1 as propagation operator (̄ , z) , which abbreviates the construction of equivalent sources in Eqs (4) and (5) to one-way field excitation in Eq (13) as discussed in Sect. 3. The field ̄ 2+ 0 = (̄ 1+ 0 , z) (cyan wiggles) impinges on the second interface at z = z 2 and excites the reflected field ̄ 2− 0 = (̄ 2+ 0 ) (gold wiggles), where is an abbreviation of the reflection operator (discussed in Sect. 4). In the second fiber, the guided transmitted field ̄ T,guided 0 (black wiggles) and the radiation field ̄ T,rad 0 (purple wiggles) are excited (as discussed in Sect. 4). By applying the propagation operator to the reflected field ̄ 2− 0 , the field ̄ 1− 0 = (̄ 2− 0 , − z) denotes the first arrival on the right-hand side of the dashed interface at z = z 1 + with ↓ 0 (orange wiggles). To improve readability, we let ̄ 0 =̄ 1− 0 denote the fields at the dashed interface, and suppress the superscript 1− in the remainder of this section.
We temporarily remove the excitation ̄ I from the first fiber and place equivalent sources for the field ̄ 0 (dark blue straight arrows) at z = z + 1 where z + 1 = z 1 + with ↓ 0 as shown in Fig. 7. The evolution of the field ̄ 0 along the loop contour to ̄ 1 on the right-hand side of the dashed interface may be evaluated by a loop operator that we define as The field ̄ 0 (dark blue wiggles) excites a backward-propagating guided field ̄ R,guided 1 (black wiggle) and a radiating field ̄ R,rad 1 (purple wiggles) in the first fiber and excites a forwardpropagating reflected field ̄ 1+ 1 = (̄ 0 ) (dark red wiggles) that propagates to the second interface and arrives as ̄ 2+ 1 = (̄ 1+ 1 ) (red wiggles). At the second interface, it excites a forward-propagating guided field ̄ T,guided 1 (black wiggles) and a radiating field ̄ T,rad 1 (purple (29) 1 = (̄ 0 ) = ( ( ( (̄ 0 ), z)), − z). Fig. 7 Sketch of the reflected field in the first fiber and the transmitted field in the second fiber due to a specific excitation ̄ 0 at the dashed interface that was generated in Fig. 6. Along the trajectory of the first loop, the excitation ̄ 0 reflects off the first fiber (̄ 1+ 0 ) , propagates to the second interface (̄ 2+ 0 ) , reflects (̄ 2− 0 ) , and propagates back to the dashed interface ̄ 1 wiggles). The reflection ̄ 2− 1 = (̄ 1 2+) (blue wiggles) propagates back to the dashed interface and arrives as ̄ 1 = (̄ 2− 1 ) (cyan wiggles) and completes the loop operator. The application of the loop operator may be continued indefinitely to obtain ̄ l with l ∈ 2, 3, … , and allows to accumulate all contributions to the scattered fields ̄ R,guided l , ̄ R,rad l , ̄ T,guided l , and ̄ R,rad l in the process. The repeated application of the loop operator leads to a Neumann series ̄ = ∞ ∑ n=0 n (̄ 0 ) that describes the converged field we wish to evaluate. Depending on the gap configuration and the material properties, convergence of the Neumann series may be slow. Moreover, the finite numerical precision of the electromagnetic fields and the loop operator may interfere with reaching a stabilized situation in which all field resonances in the cavity have settled.
Our strategy is to rewrite the loop operator in an iterative form, where every step leads to a new geometric progression of a known field, and a residual field. By applying the Gram-Schmidt process, we separate the one-time looped field ̄ 1 in terms of the original field ̄ 0 and a residual field ̄ 1 , i.e., By applying the loop operator on Eq. (30), the second loop ̄ 2 = (̄ 1 ) can be written as As the loop operator continues to iterate, the n-times looped field n (̄ 0 ) reads The total electromagnetic field at the dashed interface ̄ = ∑ ∞ n=0̄ n comprises the sum of geometric series for the vector fields ̄ 0 , ̄ 1 , (̄ 1 ) , 2 (̄ 1 ) , … , all with a factor s 0 . Because |s 0 | < 1 , the total electromagnetic field converges to We want to avoid the numerical computation of a large number of loops of the first residual field n (̄ 1 ) for n = 1, 2, … in Eq. (33). Therefore, we also make a decomposition of the loop of the first residual field by writing which yields a second residual field ̄ 2 that generally carries less power than the first residual field ̄ 1 , because in the loop ̄ 1 → (̄ 1 ) , power is coupled into the two waveguides. Analogous to Eqs (30)-(32), we express the n th -loop of the first residual field ̄ 1 , (32) n =̄ 0 s n 0 +̄ 1 s n−1 0 + (̄ 1 )s n−2 0 + … + n−1 (̄ 1 ).
which comprises geometric series with a factor s 1 for the first residual field ̄ 1 and looped versions of the second residual field ̄ 2 . Because the factor |s 1 | < 1 , the infinite sum in Eq. (33) converges to By substituting Eq. (36) in Eq. (33), the total electromagnetic field ̄ becomes a weighted sum of ̄ 0 , ̄ 1 and an infinite sum over looped versions of ̄ 2 . For every subsequent loop, one can make a new decomposition according to Eq. (34) with (̄ 1 , s 1 ,̄ 2 ) → (̄ l , s l ,̄ l+1 ) , where l denotes the loop number.
An iterative scheme for the evaluation of the Neumann series for the residual field is given by Eq. (36) by letting (̄ 1 , s 1 ,̄ 2 ) → (̄ l , s l ,̄ l+1 ) . The converged total electromagnetic field is evaluated by The sum over l should converge rapidly, because with every loop, the newly constructed residual field ̄ l carries less power than the previous residual field ̄ l−1 . For optical fiber connections, it suffices to evaluate the first three loops with ∈ (1, 2, 3) in the sum in Eq. (37).
With the aid of Eq. (37), and by re-using the coefficients s l , the converged electromagnetic fields ̄ 1+ , ̄ 2+ and ̄ 2− are easily evaluated as well. The converged-state electromagnetic field ̄ +̄ 1+ is the reflected field in the first fiber, whereas ̄ 2+ +̄ 2− is the converged transmitted field in the second fiber. Also, we have kept track of the contributions to the modal amplitudes of the guided modal electromagnetic fields in both fibers as the iterative scheme progressed.
Upon combining the excitation situation in Fig. 6, and the converged response due to the equivalent source at the dashed interface in Fig. 7, the multiple scattering problem is properly solved.
To give a feel for the speed of the convergence of the outlined approach, consider two identical standard single-mode fibers that are perfectly aligned along the transverse coordinates, and only separated by a vanishingly small gap of air with z ↓ 0 . Upon expansion of the fundamental mode in the Wilson basis ( ̄ I in Fig. 6) at a wavelength of 1310 nm , we truncated the set of Wilson basis coefficients when unit power was achieved with an absolute accuracy of 10 −6 (further details on the expansion in Sect. 7). At the fiber to air interface at z = z 1 , the first reflected field carried by ̄ R,guided 0 +̄ R,rad 0 is a factor r = 3.39 × 10 −2 of the incident power. The power carried by the field ̄ 2+ 0 that is incident on the air to fiber interface at z = z 2 is (1 − r ) = 0.9661 . The associated reflected power is a factor r lower, rendering the power carried by the field ̄ 0 effectively (1 − r ) r = 3.27 × 10 −2 . After evaluation of the first intra-cavity loop ̄ 0 →̄ 1 , the power carried by ̄ 1 is as low as (1 − r ) 3 r = 3.76 × 10 −5 . Due to the geometric progression in Eq. (33) with s 0 = 3.39 × 10 −2 (obtained by means of Eq. (30)), all the fields have converged with an absolute accuracy of 10 −6 after a single loop correction! The power carried by the first residual field ̄ 1 is as low as 6.67 × 10 −8 . Without the Neumann-series geometric progression, i.e. for s 0 = 0 and a first residual field ̄ 1 carrying 3.76 × 10 −5 of power, it would take at least three loop evaluations to achieve the same level of accuracy. We apply the mode-matching technique to single-mode fiber connections with a modefield diameter mismatch in Sect. 7 and apply an axial misalignment in Sect. 8. In Sects. 9 and 10, we consider multi-mode fiber connections with lateral and axial misalignment respectively for encircled-flux compliant and overfilled launch conditions. As a closing example in Sect. 12, we evaluate the change in the near-field pattern when coupling a target encircled-flux compliant launch to a regular and a trench-assisted multi-mode fiber.

Single-mode fiber connections with a vanishingly small gap
To model a physical contact connection, we apply the mode-matching approach as outlined in Sect. 6 here for a vanishingly small gap of z ↓ 0 . We assume a perfect lateral alignment and an excitation by the fundamental mode at a wavelength of 1310 nm in the first fiber. For both the transmitting and receiving fibers, we assume a step-index profile comprising a pure silica cladding and a germanium-doped core with a core radius a = 4.06 m.
We define the mode-field diameter (MFD) of the fundamental mode as the diameter of the contour where the normalized power density distribution has decayed to exp(−2) . We model eleven different receiving fibers by choosing the germanium mole fraction at equidistant samples between 0.03 and 0.05, resulting in mode-field diameters between 10 and Fig. 8 Single-mode fiber connection attenuation as function of the mode-field diameter of the receiving fiber ( MFD2 ) for two transmitting fibers with MFD1 respectively 8.7 m (ascending curves) and 10 m (descending curves). The solid red curves are obtained by the mode-matching technique, and are sandwiched between two approximation methods that effectively match either the electric field (blue + markers) or the magnetic field (green • markers) with overlap integrals. The red dashed curve is a Gaussian approximation by Marcuse (Marcuse 1977), and the thick black dashed curves that almost overlap our solid red curves are obtained with FIMMWAVE. (Color figure online) 8.7 m . With the aid of Sellmeier's equation for optical glass (Bach and Neuroth 1998), we evaluate the refractive index at the wavelength of choice and compute the modal electromagnetic fields for the pertaining optical waveguide (Smink 2009).
The attenuation curves as a function of the MFD of the receiving fiber ( MFD2 ) in Fig. 8 are shown for a launch fiber with MFD1 = 8.7 m (ascending curves) and for a launch fiber with MFD1 = 10 m (descending curves). The attenuation is evaluated as −10 log 10 (P 2 ∕P 1 ) , where P 1 and P 2 denote the total power carried by the guided modes in the first and second fibers, respectively.
The red solid attenuation curve (without markers) is obtained after solving through the procedure as outlined in Sect. 6 for a Wilson-basis discretized launch ̄ I =̄ 1 fun that excites the fundamental mode ( f un in subscript) in the first fiber (1 in superscript).
We have pre-computed a set of coefficients ̄ , ( z) using Eq. (6) for z ↓ 0 at a wavelength of 850 nm for d = 3 m for the decomposition of the modal electromagnetic fields into their equivalent one-way propagating fields in homogeneous-space (cf. Eq. (13)). To be able to re-use the coefficients ̄ , (0) at the wavelength = 1310 nm , we preserved the quantity ∕d by choosing the scaling parameter d ≈ 4.6 m for the expansion of the electromagnetic fields. Further, we included all the Wilson basis coefficients between the multi-indices min = (0, −22, 0, −22) and max = (3, 22, 3, 22) for all the field expansions, including the converged fields ̄ 1+ , ̄ 2+ , ̄ 2− and ̄ (as shown in Fig. 7), which in turn include all stabilized resonances due to the first three residual fields ̄ 1 , ̄ 2 and ̄ 3 .
For validation, we have chosen to compare our vectorial full-wave Wilson-basis discretized mode-matching results for this axisymmetric (aligned) full-contact case involving single-mode transmitting and receiving fibers with results obtained using a fully-licensed PhotonD suite (http://www.photond.com 2022). In particular, we have solved the pertaining HE 11 modes using the FIMMWAVE implementation of the EigenMode Expansion (EME) method (Gallagher and Felici 2003) with the full-vectorial General Fiber Solver (GFS), and have simulated the aligned full contact fiber connection through the Simple Joint operation. In FIMMWAVE one has to take care to set the (radial) cladding thickness judiciously. In this simple example, the difference in FIMMWAVE between setting it in the range 15 − 20 m and a larger value of 50 m is very small (about 0.01 dB). However, the results produced by FIMMWAVE are much more sensitive to the choice of cladding thickness for the case of two misaligned single-mode fibers with a gap to be discussed in Sect. 8.
The red solid line in Fig. 8 is sandwiched between two solid curves that are both obtained using approximation methods. When assuming a negligibly small reflected field ̄ R , the transverse electromagnetic fields are continuous across the fiber interface when , rendering the transmitted field ̄ T directly proportional to the excitation ̄ I . The transmitted field in the second fiber may be expressed as where the fundamental mode ̄ 2 fun and radiating fields ̄ 2,rad are orthonormal with respect to the scalar product defined in Eq. (14). The coupling coefficient = ⟨̄ 1 fun ,̄ 2 fun ⟩ is the outcome of an overlap integral and effectively matches the electric field of the excitation ̄ 1 fun in the first fiber to the electric field in the second fiber in a minimum-norm sense (Neumann 1988), and results in the solid curve with blue + markers in Fig. 8. The overlap integrals are straightforwardly evaluated in the Wilson basis through a summation of scalar products of vectorfield coefficients in Eq. (14). One might just as well have chosen for an effective magnetic field matching by letting → ⟨̄ 2 fun ,̄ 1 fun ⟩ * , where * denotes complex conjugation. This choice results in the solid curve with green • markers in Fig. 8.
The differences between the full mode-matching approach and the approximation methods are relatively small, but increase as the mode-field diameter mismatch between the two fibers becomes larger. The effective electric (or magnetic) field matching with overlap integrals does not account for a difference in the impedance of the electromagnetic fields across the connection interface. The full mode-matching approach does achieve continuous transverse electric and magnetic fields across the connection interface. In case an approximation method using only overlap integrals is preferred, one would do best to combine both choices for the coupling coefficient by letting 2 = ⟨̄ 1 fun ,̄ 2 fun ⟩ + ⟨̄ 2 fun ,̄ 1 fun ⟩ * . The red dashed curve in Fig. 8 is evaluated by Marcuse's scalar approximation that uses a Gaussian electric field distribution (Marcuse 1977). To make a fair comparison, we configured the width of the Gaussian electric fields according to Marcuse's approximation equation given in (Marcuse 1977) that depends only on the normalized frequency With this careful choice for the mode-field diameter, the attenuation due to Gaussian approximation is close to the full mode-matching model, albeit a bit on the high side.

Two misaligned single-mode fibers with a gap
As in Sect. 7, we start with two single-mode fibers with a MFD of 8.7 m at a wavelength of 1310 nm . Here, we consider a fixed lateral misalignment of 3 m and break the physical contact by gradually increasing the distance between the two fibers Fig. 9 The attenuation as function of the gap size between two single-mode fibers with a 3 m lateral offset for = 1310 nm (blue • marker), = 1550 nm (black × marker) and = 1625 nm (red + marker). The black solid curve is the = 1310 nm case simulated with FIMMWAVE. The attenuation due to power coupling into the fundamental mode of the receiving power is shown above the vertical scale brake, whereas the curves below the scale brake are due to the total received power including radiating fields. (Color figure online) through a set of propagation operators (̄ , z) → sweep (̄ , m z) with m = 2, … , M that makes use of the repeated application of Love's equivalence principle and re-use of the propagation operation ̄ , ( z) for the first step (de Hon and Floris 2018). Specifically, we evaluated ̄ , ( z) in Eq. (6) with z = 35.1 nm for Wilson basis functions that are scaled with d = 3 m , and included all Wilson-basis test functions w that satisfy | x | + | y | ≤ 6 and | n x | + | n y | ≤ 6 with respect to all the Wilson-basis source functions w . Re-use of ̄ , ( z) is possible by preserving the quantities d∕ and ( z)∕d , which results in step sizes z = (54.1, 64.0, 67.1) nm for the respective telecom wavelengths = (1310, 1550, 1625) nm and electromagnetic field expansions using Wilson basis scalings d = (4.6, 5.5, 5.7) m.
The solid attenuation curves that are shown above the vertical scale brake in Fig. 9 are for = 1310 nm (blue • marker), = 1550 nm (black × marker) and = 1625 nm (red + marker) for four gap regions. The step from one region to the next is achieved by iterating an additional step operator step (̄ , z) → step (̄ , 2 z) that was initialized by which became available after completing the sweep in the first region. For the 1310 nm wavelength, laterally misaligned, variable gap case, the comparison of our vectorial full-wave Wilson-basis discretized mode-matching results (blue • marker) with those generated using the PhotonD suite (solid black curve) is not as straightforward. In FIMMWAVE, the Simple Joint connection has to be replaced by a Free Space Joint based on a plane-wave expansion. It is probably due to the symmetry breaking offset that the FIMMWAVE results are now much more sensitive to the radial cladding index. FIM-MWAVE produces results that match well with ours for a radial cladding thickness of 15 m . However, some caution is advised, since for cladding thickness values of 30 m to 50 m , the attenuation curves deviate from about 0.3 dB to 6 dB.
We have demonstrated that for single-mode fiber interfaces both our method and FIM-MWAVE/FIMMPROP produces rigorous full-wave solutions, where the GFS solver in the latter relies on a more traditional Bessel-function basis (or rather cylindrical harmonics in 2-D). We would like to point out that we expand the transverse fields in terms of a 2-D Wilson basis with the appealing effective localisation in both the spatial and spectral domains and the transverse shift properties. This basis serves as a universal interface. From Sect. 9 onwards, we shall discuss multi-mode applications of our method.
The three dashed curves below the scale brake represent the attenuation for the total power coupled into the receiving fiber. By including the radiating fields, this simulation effectively resembles the constructive and destructive interference in a gap between two glass materials. As the gap size increases, the spatial coherence slowly diminishes, whereas the attenuation to the guided mode attenuation increases significantly due to the fixed lateral misalignment of 3 m combined with the widening of the modal electromagnetic fields.
We checked with a separate propagation operator step that was evaluated with Eq. (6) directly for z = 70 m that the results in the fourth region were not impaired by the large number of concatenation operations that were applied to the relatively short-distance propagation operators. Furthermore, by keeping track of the total power carried by the fields in the cavity at every evaluated cross-section, we verified that the employed propagation operators properly conserve energy, confirming that sufficient interactions between Wilson-basis source and test functions were included.
In Fig. 9, the attenuation is highest for the = 1310 nm case due to the relatively small mode-field diameter. The fixed 3 m lateral offset is much less impactful on the coupling efficiency for the spatially wider fields associated with the longer wavelengths.
The total power that does not contribute to the transmission in the second fiber is reflected power in the first fiber. The part of the reflected power that is carried by the backward propagating fundamental mode is shown as the return loss in Fig. 10.
We have interpolated the reflected power curve with a cubic spline prior to evaluating the logarithm, in order to preserve the deep dips in the return loss curves that are associated with the near-perfect transmission that occurs when the cavity length approaches integer multiples of the half wavelength. The return loss for the wavelengths = 1310 nm (blue • marker), = 1550 nm (black × marker) and = 1625 nm (red + marker) respectively settle to −14.70 dB , −14.75 dB , and −14.77 dB as the gap size tends to infinity.
The solid black and green curves without markers respectively indicate an encircledflux compliant and overfilled launch for a gap between two regular multi-mode fibers with a 3 m offset that are discussed in Sect. 10.

Multi-mode fiber connections with a vanishingly small gap
We repeated the example in Sect. 7, but now for two identical cylindrically symmetrical graded-index multi-mode fibers with a parabolically-shaped germanium-doped core profile with a diameter of 50 m and a numerical aperture NA = √ n 2 co − n 2 cl = 0.2 between the core refractive index n co and the homogeneous pure-silica cladding index n cl . At a wavelength of 850 nm , each fiber has 180 guided one-way forward propagating modal We have expanded the transverse electromagnetic fields of the 180 fiber modes in a Wilson basis using d = 3 m for the Wilson basis function scaling. For the expansion, we used all the Wilson basis coefficients between the multi-indices min = (0, −28, 0, −28) and max = (3, 28, 3, 28) . We populated the dyad ̄ , ( z) (Eq. (6)) in the propagation operator (̄ , z) with all Wilson-basis test functions w that satisfy | x | + | y | ≤ 6 and | n x | + | n y | ≤ 6 with respect to all the Wilson-basis source functions w .
For every of the 180 guided modes in the first fiber as excitation, we have evaluated the total power coupled into the guided modes of the second fiber for lateral misalignments up to 6 m . Fig. 11 shows the attenuation (red • symbols), and piece wise linear interpolations (black curves). The attenuation is defined as −10 log 10 (P 2 ∕P 1 ) , where P 1 and P 2 denote the total power carried by the guided modes of the first and the second fiber, respectively.
The curves are clustered in groups, and the left-most group consist of excitations by the high-order modes from the 18 th mode group. Every subsequent group is associated a lower mode group, which contain the more core-confined modes. The identification of the mode groups is discussed in Sect. 10 (and displayed graphically by the blue • symbols in Fig. 15) in the context of an increasing gap between the fibers and a fixed 3 m lateral offset.
It is interesting to note that some curves have strong oscillations and are not monotonically increasing. We verified that these dips are not artifacts due to the truncation in the Wilson basis expansion of the electromagnetic fields. They are not the result of the number of iterations for the residual fields in the cavity being too low either. The oscillations are due to the high spatial frequency of the modal electromagnetic fields that locally couple better to the relatively small number of guided modes in the receiving fiber for certain lateral offsets. For the lower-order excitations, the number of excited modes in the receiving fiber are much larger, effectively smoothing the pertaining curves. Below, we proceed to study overfilled and encircled-flux compliant launches with the aid of the 342 modal electromagnetic fields, in order to arrive at practical attenuation curves.

Overfilled launch
A spatially and spectrally wide launch generated by a light-emitting diode (LED) has applications in legacy networks as well as in measurement instruments (Refi 1986), and optimum coupling may be achieved with a lens (Hasegawa et al. 1978;Park et al. 1999). We assume that such an overfilled launch is spatio-temporally incoherent and that the spatial interference of the 342 modal electromagnetic fields lead to speckle patterns that blend to a cylindrically symmetrical intensity pattern due to the large spectral width of the LED source. To simulate an overfilled launch incident on an interface between two fibers, we adopt the common assumption that all of the 342 modes are excited with equal power. Instead of adding the fields in a coherent manner and solving the reflection-transmission problem for a large number of speckle patterns, we treat the excitation of each fiber mode separately and then coherently add the results. The average of all the curves in Fig. 11 represents the attenuation curve of an overfilled launch which is shown with red + symbols in Fig. 12. For a lateral offset of about 2 m , the attenuation matches well with a simple approximation method (thin black solid line) given by Neumann and Weidhaas (Neumann and Weidhaas 1976).
Overfilled launches are no longer considered suitable for fiber connection measurements, because vertical cavity surface emitting laser (VCSEL) sources (deployed in large Fig. 12 The attenuation due to lateral misalignment of two multi-mode optical fibers with a vanishingly small 1 nm gap. For an overfilled launch (red + marks), the results match an approximation method (thin black solid line) when the offset is larger than 2 m . An incoherent target-EF compliant launch (blue × symbols), and a coherent target-EF compliant launch (red solid line) produce identical attenuation curves. A geometrical-optics approach (blue dashed curve) leads to a slightly steeper curve. The purple solid vertical attenuation regions correspond to launches that touch the EF upper-and lower-bounds as shown in Fig. 13. (Color figure online) datacenter networks) generate strongly core-confined launches, leading to much higher optical power coupling efficiencies.

Encircled-flux compliant launch
Every launch may be distinguished by its time-averaged near-field pattern at the end-face of the multi-mode fiber and the associated normalized cumulative radial intensity distribution is known as encircled flux (EF). Specifically for connection attenuation measurements, international standards (cf. (IEC 61280-4-1 2009(IEC 61280-4-1 , 2021) demand that the EF curve of the launch lies within tight bounds (between red keep-out zones) associated with four target power levels (green dots) that are defined at four prescribed radial coordinates as shown in the leftmost graph in Fig. 13.
We want to avoid the burden of simulating a light source (either a VCSEL or a LED) that is connected to a specific (proprietary) mode scrambling device. Instead, we aim to describe an encircled-flux compliant incident field ̄ I = ∑ m c m̄ + m by assigning appropriate modal amplitudes c m to the modal electromagnetic fields ̄ + m , that individually propagate through a factor exp(−jk z,m z) . Of course, the 342 modal amplitudes c m cannot be uniquely determined by just five constraints that are imposed on the intensity distribution. However, by demanding that the launch is spatio-spectrally incoherent (Floris et al. 2019), it is possible to obtain a unique relation between the desired radial intensity distribution I( ) and the model power distribution MPD( ) which is given by (Daido et al. 1979) Fig. 13 The cumulative near-field distributions (shown left) of the launches are encircled-flux compliant when they do not cross the vertical red dashed lines. An overfilled launch (black dashed curve) exceeds the three lower bounds. The four encircled-flux targets are vertically centered between the respective upper-and lower bounds and are highlighted with four green • marks. An incoherent launch configured by sampling a modal power distribution (blue solid curve) meets the targets. Also, the average (red solid curve) of 60 coherent launches (thin grey curves) as observed by a slow detector meets the targets. A geometrical-optics based launch (blue dashed curve) has a similar shape. Two extreme launches (purple solid curves) are configured to touch the EF upper-and lower bounds. The right figure shows the deviations of the EF compared to the geometrical-optics based launch. (Color figure online) where is the radial coordinate normalized with respect to the fiber core radius. The coefficient ∈ (0, 1) can be thought of as a continuous normalized mode group number that samples the modal power distribution function. The relation in Eq. (39) was derived by following the mode-continuum approximation, which assumes that all modes with comparable propagation coefficient carry equal power (Daido et al. 1979).
To make the connection between the modal amplitudes c m and the modal power distribution function, we first write the propagation coefficients as where k 0 is the free-space wavenumber, and f = (n 2 co − n 2 cl )∕(2n 2 co ) is the contrast in the refractive index profile. The parameter m is explained below.
In a related scalar problem with an infinite parabolic refractive index profile (Ghatak and Thyagarajan 1998), the propagation coefficients k z,m of the linearly-polarized LP m modal electromagnetic fields are of the same form as given in Eq. (40). As m → (l, m) with l and m respectively denoting the azimuthal and radial indices, the parameter m → M∕M ub directly connects the propagation coefficients k z,m → k z, (l,m) to the mode group number M = 2m + l − 1 (Ghatak and Thyagarajan 1998). The upper bound M ub = ⌊k 0 n co R f � f ∕2⌋ equals 18 for a nominal multi-mode fiber that has a core radius R f = 25 m and numerical aperture of NA = 0.2 at a wavelength of = 850 nm.
In the vectorial full-wave (VFW) case, the modes have comparable propagation coefficients that are clustered in 18 distinct groups, which we consider the vectorial full-wave counterpart of mode groups. We count the number of modes in each group as N mg (m) and sample the modal power distribution in Eq. (39) for each mode with the VFW variable m that is connected to the VFW propagation coefficients k z,m in Eq. (40). With these ingredients, we select the modal amplitudes as where exp(j m ) is a phase term with m ∈ (0, 2 ) that will be discussed below with respect to three different launches.
First, we assume a launch that is spatio-temporally incoherent in the sense that the contribution of every mode to the near-field pattern (and the reflection-transmission problem) is considered to be independent of the contributions of all other modes. The EF curve (blue solid curve in Fig. 13) that is obtained by summation of all the separate 342 cumulative near-field patterns that are weighted by |c m | 2 (rendering the phase m of every mode irrelevant) nicely checks all four EF targets. The associated attenuation curve is shown in Fig. 12 with blue × symbols. For reference, we have added the EF curve that is obtained by following a geometrical-optics approach (blue dashed curve) in Fig. 13 that also samples the modal power distribution in Eq. (39) as outlined in (Floris et al. 2019). The associated attenuation curve (blue dashed curve in Fig. 12) follows the vectorial full-wave approach very well, albeit that it develops a visible mismatch from lateral offsets larger than 2 m . Beyond an offset of 4 m , the mismatch is about 0.05 dB . Despite a small mismatch in the shape of the encircled-flux distribution, the (real-ray) geometrical optics model yields a slightly pessimistic attenuation prediction because the guided rays are confined to the core region only, and can not represent the presence of a tail in the cladding region. Second, we have configured a temporally incoherent launch that consists of a sequence of spatially coherent fields, where every field has has coefficients c m according to Eq. (41) however with m ∈ U(0, 2 ) sampled arbitrarily from a uniformly distribution for all 342 modes. Every field represents a non-symmetric speckle pattern with its own radial intensity distribution as shown by the 60 thin grey curves in Fig. 13. The average intensity distribution as observed by a slow detector camera is cylindrically symmetric and also satisfies the EF targets (solid curve with red + markers in Fig. 13). The associated mean attenuation distribution (red curve overlapping the one with the blue × markers), as well as the ensemble of individual field contributions (grey curves) are shown in Fig. 12. The variation in the attenuation among the thin grey curves for a fixed lateral offset connection can be observed as modal noise, and can deteriorate the eye diagram of a modulated signal generated by a spectrally narrow laser source (Epworth 1979;Efimov 2019;Koonen 1986).
As a third approach, we applied the Nelder-Mead minimization algorithm to the 342 unknown amplitudes in an attempt to minimize the absolute distance between the EF curve and the four EF targets. Although not shown here, the approach to solve the underdetermined system of equations leads to a solution (de Hon et al. 2015). Serving merely as an example, we applied this approach also to the upper and lower bounds of the target EF values and highlighted the results through solid purple curves in Fig. 13. The variation in attenuation associated with these two extreme EF launches are shown by vertical purple solid lines in Fig. 12. It is by coincidence that this range seems to align with the variations obtained by 342 spatially coherent speckle-pattern launches (from the second approach), because the EF bounds were originally determined to limit the variation in attenuation by no more than 0.1 dB for relatively low-loss connections.
We proceed with the introduction of an increasing gap between the two fibers.

Two multi-mode fibers with a gap
To break the physical contact between the two multimode fibers, we repeat the procedure as outlined in Sect. 9 but now for gap sizes that are integer multiples of z = 70.2 nm . In this section, we consider a fixed 3 m lateral offset between the transmitting and receiving fiber. By (again) using all the Wilson basis coefficients between the multiindices min = (0, −28, 0, −28) and max = (3, 28, 3, 28) , the total transmitted and reflected power (including radiating fields) deviates less than 0.1 dB compared to the unit launch power (for every launch mode), we could increase the gap size to 70 m . In this case, the diffracted field inside the open cavity formed by the gap extends to a region of about 70 * (28 * 3) 2 ∕(0.850) 3 ≈ 0.8 million cubic wavelengths. For each of the 180 guided modes in the launch fiber, the attenuation as a function of the separation z are highlighted by red • symbols in Fig. 14.
The solid black curves are sinusoidal interpolations, and the blue − symbols indicate the attenuation when considering the total received power (including radiating fields) with respect to the incident power.
It is interesting to observe that the attenuation curves are clustered in distinct groups. The highest cluster of attenuation curves are due to 19 modal excitations from the 18 th mode group. Due to the lateral offset of 3 m , these specific high-order modal launches couple about half of the power into the guided modes of the receiving fiber (attenuation about 3 dB at z = 0 ), whereas the other half couples into the (approximated) radiating fields (blue horizontal − symbols at 0 dB). Despite the 3 m lateral offset, the excitation of radiating fields is insignificant for modal launches that are associated with the first 14 mode groups.
To better visualize the behavior of all 180 modes as the gap size increases, the minimum and the range of the attenuation curves are shown with blue • symbols in Fig. 15 as a function of the propagation coefficient k z,m of the respective modes. The number of modes in Fig. 15 The top graph for 180 individual modal excitations (labeled by k z,m ) the minimum attenuation into the guided modes of the receiving fiber for the three gap ranges listed in the legend. The lower graph shows the variation in attenuation. The number of modes that form a mode group are given between brackets each mode group is denoted between parentheses above the bottom subfigure in Fig. 15. The smaller the propagation coefficient k z,m for the cluster of modes, the higher the associated mode group number. Note that the fundamental mode, labeled m = 1 , has the highest propagation coefficient k z,1 and is the only mode in the first mode group. The fundamental mode launch diverges slowly, and excites only guided modes in the receiving fiber, even if the gap size is increased to 71 m.
The variation in attenuation is about 0.6 dB for all modes as shown by blue • symbols in the bottom subfigure of Fig. 15. This 0.6 dB variation was also seen in the single-mode example in Fig. 9. It is noteworthy that the distribution of curves at z ↓ 0 in Fig. 14 are identical to the distribution of curves in Fig. 11 along the vertical axis for a 3 m lateral offset.
For a second region with a separation of z = 17.5 m , we employed a spatially larger propagation operator step (̄ , z) including interactions up to | n x | + | n y | ≤ 14 . This operator was evaluated directly using Eq. (6). The attenuation curves shown in Fig. 16 were made with an additional sweep operator sweep (̄ , z) that was configured for z = 35.1 nm and was chosen spatially smaller, including interactions up to | n x | + | n y | ≤ 6.
The attenuation curves are still somewhat clustered for the mode groups 15 to 18, as can also be seen by the red × markers in the top graph of Fig. 15. However, within these mode groups, the attenuation curves are very different and include curves that vary about 0.5 dB and include curves that are almost horizontal as seen by the spread of red × markers in the lower subfigure of Fig. 15. Also, the sinusoidal interference patterns in Fig. 16 no longer synchronized.
More interesting are the peculiar near-constant attenuation curves for several mainly high-order modal launches. One might expect that the specific combination of a relatively large gap and a fixed 3 m lateral misalignment counteracts the interference in the cavity. However, the total received power in the receiving fiber is also near-constant for these specific launches. We are confident that we included sufficient Wilson-basis source and test functions in the expansions that rely on Eq. (6), because the power carried by the electromagnetic fields remains unchanged across the gap. Also, the sum of the transmitted and reflected power remains equal to the incident power for every step in the sweep. These power-balance checks are sensitive to incautious truncations of the field expansions and propagation operators.
As the gap-size is increased to about z = 35 m , the near-constant attenuation behavior is now also seen for relatively lower-order mode groups, as seen by the black + markers in Fig. 15. We suspect that the relatively diverging fields in the cavity cause the gap-size insensitivity for specific regions of z.
Although we also evaluated the curves for the region near z = 70 m , we do not want to clog Fig. 15 with more markers. For this distance, the diverging beam fields that are excited by the high-order modes in the first fiber start to noticeably exceed the edge of the truncated spatial domain.
To ascertain that the computational domain is large enough, we ran three simulations including all Wilson coefficients for which |n x |, |n y | ≤ 28 , or 30, or 32, and observed that the power loss of the residual fields (corrections to the first loop) for the most divergent m = 18 modal field excitation reduces from 12% , to 4% , and 0.24% respectively. Here it should be noted that this first residual field only carries about 10 −3 of the input power.
We have combined all computations in Fig. 17 for a spatio-temporal incoherent EF launch (blue curve with blue × markers) and for an overfilled launch (black curve with red + markers).
The overfilled launch excites all modes equally strong (cf. Sect. 9.1), whereas the target EF launch is more core confined (cf. Sect. 9.2).
It is noteworthy that the attenuation curves at z = 0 for both launches agree with the curves in Fig. 12 for the physical contact case with 3 m lateral offset. Although hardly noticeable from the figure, the envelope of the overfilled launch case decays slightly faster than the encircled-flux launch case, due the the stronger presence of the high-order modes. We proceed with a few examples that show the evolution of the EF distribution when coupling a target launch from a regular fiber to a trench-assisted multi-mode fiber.

Regular multi-mode fibers connected to trench-assisted multi-mode fibers
Nowadays, trench-assisted multi-mode fibers are commonly deployed in datacenter networks and there they can get connected to regular multi-mode fibers that are already part of existing installed cabling infrastructure. To ensure that an optical link still meets the powerbudget requirements when fiber types are mixed, the attenuation at each interface should be sufficiently low, preferably below several tenths of a decibel. To make a fair assessment on the connector performance of trench-assisted multi-mode fiber cables, a global debate has started on the applicability and functionality of an encircled-flux target launch (and the respective bounds) in trench-assisted multi-mode fibers. Because trench designs are different among fiber suppliers, we did a preliminary study on the attenuation for various trench dimensions and trench positions in the cladding region (de Hon et al. 2015). In that assessment, we assumed that both the regular and trench-assisted fibers have a parabolic core refractive-index profile with a diameter of 50 m and a numerical aperture of 0.2. Even though the presence of the trench effectively narrows to modal electromagnetic fields in the trench-assisted fiber, it hardly impacted the Fig. 18 Both subfigures show the relative change in EF when a target EF launch in a regular multi-mode fiber couples into a trench-assisted multi-mode fiber. Left: both fibers have nominal dimensions and are perfectly aligned (blue curve without marker), with a 3 m lateral offset (black curve with + marker) and with a 6 m lateral offset (red curve with × marker). Right: the trench-assisted multi-mode fiber has a 49 m core diameter and has no lateral misalignment (blue solid curve without marker), or has a 3 m lateral offset (black curve with + marker) or has a 6 m lateral offset (red curve with × marker). The blue dashed curve represents the EF change for a perfectly aligned trench-assisted multi-mode fiber with a 51 m core diameter. The dashed curve with ◊ markers represents a receiving fiber with a 51.75 m core, a NA of 0.215 and a 10 m lateral offset. (Color figure online) coupling matrix for connections to a regular fiber, resulting in attenuation values that were well below a tenth of a decibel.
Here, we redo some of those simulations with the aid of the mode-matching approach that was described in Sect. 6, and look at the impact on the cumulative near-field distribution as the target encircled-flux launch couples from a regular fiber to a trench-assisted fiber. We assume that a typical trench has a width of about 3 m and a constant trench refractive index n tr that satisfies the ratio (n cl − n tr ) f = 0.0045(n co − n cl ) relative to the core refractive index contrast f = (n 2 co − n 2 cl )∕(2n 2 co ) . The trench position is assumed to be 1 m away from the parabolic core in the cladding region.
As shown by the blue solid curve (without markers) in the leftmost subfigure of Fig. 18, the tail of the near-field pattern is a bit more core confined (the curve is no longer horizontal). More importantly, it stays well within the bounds for a standards compliant launch, which means that its application can in principle be extended to connection attenuation measurements with modern multi-mode fibers. Even for an additional large lateral misalignment of 3 m (black curve with + marker) and 6 m (red × marker), the launch remains EF compliant.
In the rightmost subfigure of Fig. 18, the parabolic core of the trench-assisted multimode fiber is reduced to 49 m , and the three solid EF curves effectively move further upwards and touch the upper bounds. When the trench-assisted multi-mode receiving fiber has a 51 m core diameter, the near-field pattern widens and the EF curve shifts downwards (blue dashed curve in Fig. 18), and gets close to the lower bounds.
A numerical aperture mismatch affects the width of the near-field distribution in the opposite way as a core diameter mismatch. The near-field pattern widens when the receiving fiber has a smaller numerical aperture than the launch fiber, whereas the near-field pattern becomes more core confined when the numerical aperture is larger.
The specific change in the cumulative near-field pattern due to a simultaneous mismatch in the core diameter and the numerical aperture might be revealed through deliberate lateral misalignments. Consider for instance the cumulative near-field pattern in a trenchassisted receiving fiber that has a core diameter of 51.75 m , a NA of 0.215, and a lateral misalignment of 10 m with respect to the launch fiber. The larger receiving fiber CD (for equal NA) stretches the support of the EF towards larger radii, and the associated EF curve downwards, whereas the larger receiving fiber NA (for equal CD) shrinks the support of the EF towards smaller radii, and the associated EF curve upwards, albeit not in equal measure. That is to say, a receiving fiber CD increase is not fully cancelled by a simultaneous appropriate NA decrease. This is revealed by the asymmetric change in the EF (purple dashed curve in Fig. 18 for a deliberate lateral offset. Hence, it might be possible to unravel the two parameters by monitoring the EF changes across a connection interface, and thus provide an indirect measurement of such small-scale differences.

Discussion and conclusions
We have described a Wilson-basis discretized mode-matching technique that evaluates the converged electromagnetic fields for configurations that have two consecutive interfaces forming a resonating cavity. We gave examples of two optical fibers that are laterally and longitudinally misaligned in a homogeneous air medium. To avoid countless evaluations of fields that propagate forward and backward between the two interfaces (which we call looped fields), we devised an iterative method that yields a Neumann series of looped fields with every step.
The sum of these Neumann series converges very rapidly for fiber-air interfaces. For the limiting case of two identical and perfectly aligned fibers, the evaluation of the first Neumann series already yields the same level of accuracy as evaluating three consecutive loops of the initial looped field. To make sure all resonances have decayed, we included three Neumann-series (converged loops) for all simulations.
For single-mode fiber connections with a vanishingly small gap, we have evaluated the attenuation curves due to an impedance mismatch caused by a numerical aperture mismatch. The curves are right in the middle of those that are obtained with a common approximation technique that assumes a vanishing reflection and projects the transverse electric (or magnetic) field of the excitation to the transmitted field using overlap integrals. For the mismatches in the mode-field diameter that can occur for typical single-mode fiber applications (up to 2.5 m ), the attenuation difference between the approximation method and our approach is only a few hundredth of a dB.
However, the approximation methods are no longer applicable when the gap size increases. We have demonstrated the capability of the full mode-matching method by expanding the gap up to 70 m and a 3 m lateral misalignment for three telecom wavelengths. The consecutive application of Love's equivalence principle and a propagation operator for just 32.1 nm allowed us to bridge the large gap of 70 m , sparing us the time-consuming effort to evaluate the propagation operator in Eq. (6) for every step. Even though the energy distribution among the Wilson-basis source and test functions in the propagation operator suggests that (spatial) truncations cannot be made (cf. Figs 4 and 5), we could truncate the operator due to the spatio-spectral confinement of the electromagnetic fields that are associated with the optical fibers. The conservation of the power carried by the fields that propagate across the gap are a valuable indicator for field degradations due to truncations. We were able to preserve the power balance to at least five digits. The attenuation and return-loss curves show strong interference effects for the short gap, but also decay of the resonances and increase in the overall attenuation as the gap size increases.
For multi-mode fiber connections, we have determined the attenuation curves for lateral and longitudinal misalignment for the excitation of each of the 180 guided modes (not counting complex conjugate counterparts). Particularly the high-order mode excitations reveal interesting attenuation dynamics, such as variations in the coupling efficiency as the lateral misalignment increases (cf. Fig. 11) or damping interference effects for the specific gap size near 17 m (cf. Fig. 16).
A multi-mode fiber launch in practice comprises a specific distribution of modal electromagnetic fields. We first evaluated a naive implementation of an overfilled launch that excites all modes equally strong and in a spatio-temporal independent manner. The resulting attenuation curve matches well to approximations by Neumann and Weidhaas. To construct a target EF launch, which is standardized for attenuation measurements, we sampled a modal power distribution to end up with a specific linear combination of modal amplitudes. By distributing the phases of every mode randomly, we generated a large number of speckle patterns, each with its own cumulative near-field pattern and attenuation. The average cumulative near-field pattern approaches the target encircled-flux levels, while the variation in the attenuation due to the speckle patterns is a measure for the modal noise (cf. Fig. 12). Similar to the single-mode fiber case, the repetitive application of Love's equivalence principle and a short propagation operator of just 70.2 nm , allowed us to cover a distance of 70 m . In this case, the diffracted field inside the open cavity formed by the gap extends to a region of about 1 million cubic wavelengths.
As a concluding example that couples a launch from a regular nominal multi-mode fiber to a modern trench-assisted multi-mode fiber, we have considered the impact on the cumulative near-field pattern rather than the attenuation. Even though the presence of the trench in the receiving fiber effectively narrows the intensity pattern, the EF curve stays within the standardized bounds for encircled-flux compliant launches. This means that existing modal conditioners can be extended with trench-assisted multi-mode fibers to generate a standards compliant launch that is in principle suitable for connection attenuation measurements with modern-day multi-mode fibers. If the receiving fiber has a smaller core diameter and/or a larger numerical aperture, the transmitted field shrinks further. Vice versa, if the receiving fiber with a larger core diameter and/or a smaller numerical aperture, the transmitted field widens. We have shown that a simultaneous core diameter and numerical aperture mismatch across a fiber connection can manifest in asymmetrical EF curves when applying a deliberate lateral misalignment. A model-based reconstruction is a promising means to determine the actual fiber geometry parameters of the receiving fiber.

Funding Information
The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
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/.