Full-Vectorial Fiber Mode Solver Based on a Discrete Hankel Transform

: It is crucial to be time and resource-efﬁcient when enabling and optimizing novel applications and functionalities of optical ﬁbers, as well as accurate computation of the vectorial ﬁeld components and the corresponding propagation constants of the guided modes in optical ﬁbers. To address these needs, a novel full-vectorial ﬁber mode solver based on a discrete Hankel transform is introduced and validated here for the ﬁrst time for rotationally symmetric ﬁber designs. It is shown that the effective refractive indices of the guided modes are computed with an absolute error of less than 10 − 4 with respect to analytical solutions of step-index and graded-index ﬁber designs. Computational speeds in the order of a few seconds allow to efﬁciently compute the relevant parameters, e.g., propagation constants and corresponding dispersion proﬁles, and to optimize ﬁber designs.


Introduction
Optical glass fibers can support arbitrary numbers of guided modes depending on the actual waveguide design. In combination with the linear and nonlinear material dispersion, the propagation constants of these modes, i.e., the waveguide dispersion, determines the overall dispersion (of whatever order) of the fiber, which is relevant for various applications in pulsed or CW operation. For example and in combination with the split-step approach [1] or the frequency-domain Runge-Kutta approach [2], the nonlinear pulse propagation can be analyzed numerically. Such studies allow us to understand and tailor effects such as the pulse broadening for supercontinuum generation either in the fundamental mode [2] or even in higher-order modes (HOMs) [3]. Similarly, the possibility to tailor and design the HOM dispersion is attractive for enhancing the performance of telecommunication links, see for example [4]. Even novel types of fiber-based quantum sources have been realized by exploiting the properties of HOMs [5]. In CW operation, intra-modal dispersion of specialty fiber designs allows establishing phase-matching for efficient nonlinear frequency generation via four-wave mixing, e.g., in the short-wavelength UV regime [6]. Another crucial example are few-mode or multimode fibers for spatial multiplexing in next-generation telecommunication networks [7] or lens-less endoscopy [8] as the differential modal dispersion is relevant (besides other effects) for the propagation robustness of the individual modes [9]. Recently, machine learning algorithms were used for the first time to optimize a fiber design with respect to the propagation robustness of the guided modes [10]. These examples demonstrate that there is a large need for time and resource-efficient concepts to compute and optimize the guided modes (the propagation constants) of arbitrary fiber designs. Indeed, a lot of fibers have a rotationally symmetric waveguide profile, in particular as fabrication concepts of fiber preforms such as the widely used modified chemical vapor deposition (MCVD) [11] can fabricate symmetric profiles with great flexibility. Thus, the work presented here is focused on fibers with rotationally symmetric waveguide designs. In addition, the focus here is on computing the vector modes as scalar approximations of the corresponding wave equations are only valid under certain boundary conditions such as weak guiding (small refractive index variations). In addition, modes with orbital angular momentum (OAM), which became very prominent in recent years for modal multiplexing in telecommunication network [12] or quantum cryptography through fiber links [13], are conveniently described as a superposition (change of basis) of vector modes.
To numerically compute the guided fiber modes of a given refractive index profile, a transversal discretization concept has to be chosen [14]. This can be a local strategy such as a finite element [15] or a finite difference [16][17][18] approach for which a discretization i.e., mesh (either uniform or non-uniform) is used to linearize the differential operators in the corresponding wave equation(s). Afterward, a corresponding set of equations has to be solved. As an opposite approach, global methods are based on expanding the mode profiles as finite series of corresponding basis functions. Afterward, Galerkin or collocation methods are used to obtain algebraic sets of equations whose unknowns are the series expansion coefficients, see e.g., [19][20][21]. In combination with the powerful tool of the fast Fourier transform (FFT), a 2D Fourier expansion has been used in [22] to compute the vectorial modes of arbitrary 2D waveguide profiles with acceptable computational resources. Although computational times were rather slow (around a few minutes), the inherent flexibility of the series expansion approach yielded small errors around 10 −3 regarding the computed propagation constants for various fiber designs.
In the tradition of the global series expansion concepts, it will be demonstrated here for the first time, that a 1D Fourier-Bessel series expansion in combination with a discrete Hankel transform (DHT) enables fast (<10 s) and accurate (errors as low as 10 −4 ) computation of the guided vector modes and their propagation constants in rotationally symmetric fibers.

Coupled Set of 1D Wave Equations
In dielectric materials such as fused silica (no free charges and currents), Maxwell's equations are with absolute permittivity = r 0 ( r is the relative permittivity i.e., the dielectric constant) and absolute permeability µ = µ r µ 0 (µ r is the relative permeability). In dielectric materials, it is known that µ r ≈ 1 and r = n 2 with refractive index n. In the following, it is assumed that the refractive index has a certain spatial profile, i.e., ( r) = const. Using the common Ansatz E( r, t) = E( r)exp(−iωt) and multiplication of the third Maxwell equation with ∇× as well as using the fourth Maxwell equation yields For a homogeneous medium with ( r) = const, Equation (3) would yield the wellknown Helmholtz equation. To find solutions for ( r) = const, the following strategy can be used. Applying where the first term vanishes due to Maxwell's equations, Equation For ( r) with cylindrical symmetry ( ( r) − → (r)), it is advantageous to work in cylindrical coordinates, i.e., which yields a set of coupled equations for the E ρ and E φ components At this point, it is crucial to understand that if solutions for E ρ and E φ are known, the E z component as well as the components of the magnetic H field would follow from Maxwell's equations. To find the solutions for E ρ and E φ , the common Ansatz with integer m ≥ 0 yields For each m, Equations (12) and (13) constitute a coupled set of eigenvalue equations for A ρ (ρ) and A φ (ρ) with eigenvalue β 2 . In [22] the very same set of equations has been used in Cartesian coordinates and solutions were obtained by the aforementioned global expansion of the field amplitudes A x (x, y) and A y (x, y) as 2D finite Fourier series. Using the fact that the differential operators transform well under a 2D FFT, i.e., they correspond to products in the frequency space, solutions were obtained by finding the eigenvectors and eigenvalues of the corresponding matrix via algebraic routines. Indeed, it appears to be attractive to expand A ρ (ρ) and A φ (ρ) as 1D finite Fourier series as well and to follow the same strategy as in [22]. However, this is not possible for various technical reasons. In particular, it would be required that the corresponding sampling includes data points at ρ = 0, which is incompatible with the 1/ρ and 1/ρ 2 operations in Equations (12) and (13). As a solution to this issue and as mentioned, a finite Fourier-Bessel series expansion in combination with a DHT will be used in the following. Before dealing with the details of the actual numerical implementation and its validation in Section 5, the DHT and the transformation rules for the individual differential operators are introduced in the following sections.

The Discrete Hankel Transform
In [23], it has been discussed that any space-limited function f (ρ) with f (ρ > R) = 0 can be expanded as a finite Fourier-Bessel series where J n (ρ) is the n-th order Bessel function and j ni its i-th root. The Fourier-Bessel coefficients f i are given by where F(ν) is the n-th order continuous Hankel transform of f (ρ), which is defined as Although F(ν) cannot be band-limited if f (ρ) is space-limited, F(ν) can still be expanded analogously as a finite Fourier-Bessel series if it is effectively zero for ν > V: with coefficients where the n-th order continuous inverse Hankel transform has been used. Sampling f (ρ) at discrete points ρ nk = j nk R/j nN and sampling F(ν) at discrete points ν nm = j nm V/j nN yields Now, if for a given R (or V), V (or R) is chosen such that V = j nN /R (or R = j nN /R), Equations (20) and (21) define a DHT that links the discrete data sets F(ν nm ) and f (ρ nk ). Through the remainder of this report, the following slightly different DHT between scaled data sets F(ν nm )/J n+1 (j nm ) and f (ρ nk )/J n+1 (j nk ) will be used, which has been introduced in [23,24]: The corresponding inverse transform is The transformation kernel obeys T nN m,k = T nN k,m . An important advantage of this Kernel (which is not unique, see [23]) is, that it satisfies Parseval's theorem, i.e., energy is conserved during the actual transformation. Note that the DHT can readily deal with complex data sets because the real and the imaginary part transform independently, i.e.,

Transformation of the Differential Operators
Now that the basics of the DHT are established, it remains to be investigated how the differential operators in Equations (12) and (13) transform under the DHT. There are two types of operators that are relevant: The Bessel operator B n = ∂ 2 /∂ρ 2 + ρ −1 ∂/∂ρ − n/ρ 2 of order n, which occurs in Equations (12) and (13).

2.
The common derivative operator ∂/∂ρ, which occurs as ∂/∂ρ ln( ) in Equations (12) and (13) and additionally in Equation (12) as a term ∂/∂ρ A ρ ∂/∂ρ ln( ) . In the following, it will be demonstrated that the Bessel operator B n transforms well for any DHT order n, i.e., it corresponds to a simple product in the frequency space. It will also be shown that a proper transformation rule of the operator ∂/∂ρ can be established for a 0-th order DHT. Thus, for actually solving the coupled equation set given by Equations (12) and (13), a 0-th order DHT will be used (see Section 5). However, this restriction to a 0-th order DHT establishes a crucial advantage: because the sampling of the DHT is performed at the zeros of the corresponding Bessel functions and as J 0 (ρ) does not have a zero at ρ = 0, there are no issues with the remaining 1/ρ and 1/ρ 2 operations in Equations (12) and (13).

Transformation of the Bessel Operator
For the continuous Hankel transform of a function f (ρ) as defined in Equation (16) it is known that g(ρ) = B n f (ρ) transforms as In the following, it will be shown that a corresponding transformation rule exists for the DHT as well. Using the finite Fourier-Bessel series defined in Equation (14), it follows that Substitution of the well-known identities yields Discrete sampling of g(ρ) at ρ nk = j nk R/j nN and applying the DHT yields Sorting the sums and substitution of the DHT kernel (Equation (24)) leads to By using the orthogonality relation which is not a closed analytical result (see the discussions in [25]) but is numerically precise for sufficiently large N (for N > 30 the error is less than 10 −7 [25]), it follows that Now, using the definition of the Fourier-Bessel coefficients (see Equation (15)) and recalling that ν nm = j nm /R, it finally follows that which is the sought after discrete transformation rule of the Bessel operator B n .

Transformation of the Derivative Operator
For the continuous Hankel transform, it is known that Such modification (±1) of the order of the Hankel transform seems to be incompatible with the discussed discrete implementation as the sampling in the spatial space (or the frequency space) is determined by the zeros of the n-th order Bessel function. However, in the following a solution for the case n = 0 (0-th order DHT) will be presented, which relies on the workaround to not analyze the transformation behavior of the operator ∂/∂ρ but of the operator ρ −1 ∂/∂ρ. Indeed, this operator can easily be introduced in differential equations by substituting any ∂/∂ρ operation with ρρ −1 ∂/∂ρ. Similar to the last section, using the finite Fourier-Bessel series as defined in Equation (14), it follows for any DHT order n that where the identity has been used. Discrete sampling at ρ nk = j nk R/j nN and applying the DHT yields At this point the discussion is restricted to the case n = 0 because using the identity J −1 (ρ) = −J 1 (ρ), substitution of the DHT kernel (Equation (24)), and sorting of the sums yields Now the identity can be used, which has been introduced in [25] for any order n. Similar to the aforementioned orthogonality relation of Equation (33), the identity of Equation (41) cannot be deduced analytically but is numerically precise for sufficiently large N. We found that the absolute error of Equation (41) is less than 10 −5 (10 −4 ) for N > 100 (N > 20) for n = 0. Using the identity of Equation (41) yields and using the definition of the Fourier-Bessel coefficients (see Equation (15)) yields This result can also be written in a compact form as with the matrix Thus, the operator ρ −1 ∂/∂ρ transforms such that it corresponds to a scaling with 1/R 2 and a matrix multiplication in the frequency space, which can be implemented numerically quite easily. Note that the definition of the matrix χ m,i has been chosen such (without the factor 1/R 2 ) that it can be computed independently from the problem-specific choice of R (or V). This is important because numerically computing the integral in Equation (45) is feasible but requires quite some time for large N (see discussions and results presented in Section 5.2). Thus, for repeated applications of the DHT it is certainly beneficial to compute the matrix χ m,i once in advance and to simply load it during run-time. For validation of Equation (44), a simple Gaussian function f (ρ) = exp(−(aρ) 2 /2) on the spatial domain [0, R = 30a] was considered, whose derivative can easily be computed analytically. In Figure 1a, a comparison of the (sampled) analytical derivative and a 0-th order DHT (for details of the actual implementation see next section), followed by the matrix multiplication as given by Equation (44) and followed by the corresponding inverse DHT is presented for N = 100. In addition, Figure 1b presents the mean absolute deviation of the DHT-based result with respect to the analytical result for an increasing number of spatial mesh points i.e., Fourier-Bessel coefficients. As the identity of Equation (41) becomes more and more precise for increasing N, the same trend is observed in Figure 1b and for N > 100 the mean absolute deviation is already less than 10 −5 for the considered Gaussian function. Note that the convergence follows a 1/N 2 trend, see the dashed line in Figure 1b

Implementation and Validation of the DHT-Based Mode Solver
Based on the discussions and results presented so far, the coupled set of eigenwert Equations (12) and (13) was solved numerically by expanding A ρ (ρ) and A φ (ρ) as finite Fourier-Bessel series (with N − 1 coefficients) and by finding the eigensolutions, i.e., eigenvalues (squared propagation constants) and eigenvectors (Fourier-Bessel expansion coefficients of A ρ (ρ) and A φ (ρ)). The actual implementation of the DHT-based solver was written in Python 3 (64 bit) and all results presented in the following were computed on a standard computer (AMD Ryzen 52600, 6 cores at 3.4 GHz). For the computation of the eigensolutions, a dedicated solver for large (sparse) algebraic eigenvalue problems based on the implicitly restarted Lanczos method (IRLM) was used [26]. As the number of points (N − 1) is equally large in the spatial and in the frequency space, the solving could be performed in either space. The transformation between these spaces is given by the 0th-order DHT in combination with the transformation rules of the differential operators. In the used implementation, the solving was performed in the spatial domain but tests with solving in the frequency space did not reveal any significant differences with respect to the required computational time. Note that a zero-padding could be implemented in the future to enhance the resolution in the spatial space while keeping the number of Fourier-Bessel coefficients in the frequency space reasonably low, also see discussions in Section 7. Under such conditions, the solving should be performed in the lower dimensional frequency space and might enable significant reductions of the required computational times.

Step-Index Fiber Design
As a first validation case for the DHT-based mode solver, a standard step-index fiber design was evaluated. Note that for such a design, the terms that scale with ∂/∂ρ ln( ) in Equations (12) and (13) are effectively zero and only the transformation rule of the Bessel operator has to be considered. The core radius was 3 µm, the refractive index of the cladding was 1.42, the light wavelength was 1064 nm, and the core was given a 9 × 10 −3 higher refractive index than the cladding. The V-parameter of this configuration is around 2.84 (NA of 0.16) and one would expect this fiber to guide the HE 11 , the TE 01 , the TM 01 , and the HE 21 mode. Indeed, for m = 0 two guided modes (TE 01 and TM 01 ) where found by the DHT-based mode solver, for m = 1 a single-mode (HE 11 ) was found, and for m = 2 another mode (HE 21 ) was found. The computational window was set to be 20× the core radius to allow all (guided) solutions i.e., modes to be effectively space-limited. Note that no investigations regarding the optimal computational window have been carried out but it is obvious that it scales with the wavelength because evanescent fields are more pronounced at longer wavelengths. As expected, the TE 01 , TM 01 , and HE 21 modes were nearly degenerated with propagation constants that only differed at the fifth decimal place.
In Figure 2a the computed effective refractive indices, i.e., the propagation constants divided by k 0 , of the HE 11 and the TE 01 mode are shown for an increasing number N of spatial mesh points, i.e., Fourier-Bessel coefficients. In addition, Figure 2b shows the deviations of the computed effective refractive indices from their analytical solutions, which can easily be computed for step-index profiles [1]. As can be seen, already for N ≥ 50 the deviation from the analytical solutions is less than 10 −3 and becomes less than 10 −4 for N ≥ 750. The fact that the deviation is larger for the TE 01 than for the HE 11 mode is certainly explained by the more pronounced evanescent field of the TE 01 mode in combination with the Gibbs phenomenon, which causes ringing artifacts due to the idealized perfect refractive index step at the core-to-cladding boundary. This might also explain the rather slow convergence in comparison to a 1/N 2 convergence as indicated by the dashed line in Figure 2b. A typical computational time required by the DHT-based mode solver to compute the HE 11 mode is presented in Figure 2c. This is the time that is required by the corresponding Python code to construct the DHT Kernel (see Equation (24)) and to compute the largest eigenvalue for m = 1 with the IRLM solver. Two different modes of operation are compared in Figure 2c: (i) For each N the IRLM solver uses random noise as a guess or (ii) the solver uses the previous solution as a guess. Note that the latter method could only be implemented properly as the selected DHT Kernel obeys Parseval's theorem (otherwise, normalization could not be performed properly) and the speed advantage corresponds roughly to a factor of two. In general, 10% of the total computational time is required to construct the DHT Kernel, which could be avoided in the future by computing the Kernel for the various N in advance and later only load it during run-time.
Indeed, the proposed method requires sufficiently fine meshes, i.e., sufficiently large numbers of Fourier-Bessel coefficients to adequately compute the HOMs of the large-modearea (LMA) step-index designs as they have complex radial profiles. As an example, an LMA fiber with a core radius of 15 µm (NA of 0.16) was considered, which supports nine modes for m = 1. It has been found that 500 mesh points, i.e., Fourier-Bessel coefficients, are sufficient to compute the effective refractive index of the lowest-order HE 11 mode with an absolute error of less than 6 × 10 −5 . However, with the same mesh, the highest-order mode (HE 15 ) is computed with an absolute error of 2 × 10 −4 .

Graded-Index Fiber Design
As a validation case of the DHT-based mode solver for a design with non-vanishing ∂/∂ρ ln( ) terms, for which the transformation rule of the derivative operator has to be taken into account as well, a graded-index fiber design [8] with refractive index n 0 = 1.47 at ρ = 0 (the center of the fiber) was considered. The scaling parameter b controls the steepness of the graded profile and was set to 11.6 µm, which corresponds to n(ρ = 3 µm) = 1.45. In real graded-index fibers, n(ρ) is constant beyond the core radius (in the cladding). For the comparative purpose of this study, it is sufficient to deal with the non-truncated profile as it is well-known that only the highest-order modes, which are not considered in the following discussion, are affected by the truncation [27]. Thus, the lower-order modes are well space-limited and no caution regarding a proper selection of the computational window is required. The propagation constants of the linearly polarized modes of a graded-index design, which are solutions to a corresponding (approximated) scalar wave equation, can be computed analytically [27]: where the index k = 0, 1, 2 . . . corresponds to the number of radial nodes. Furthermore, the index l = · · · − 1, 0, 1 . . . corresponds to the orbital angular momentum. The linearly polarized mode with l = k = 0 corresponds to the HE 11 vector mode and the approximation actually used to obtain the aforementioned scalar wave equation (commonly for the E z field component) involves neglecting the terms that scale with ∂/∂ρ ln( ) in the wave equation. Thus, one expects that the DHT-based mode solver computes propagation constants equal to the analytical solutions of Equation (47) if the ∂/∂ρ ln( ) terms are set to zero. This is confirmed in Figure 3a: without the ∂/∂ρ ln( ) terms, the computed effective refractive index of the HE 11 mode deviates only less than 10 −5 from the analytical result for N > 100 and less than 10 −6 for N > 300. Here, the convergence follows a 1/N 2 trend, see the dashed line in Figure 3a, which can be explained by the absence of the Gibbs phenomenon. If the ∂/∂ρ ln( ) terms are included in the DHT-based mode solver, the deviation is larger but still small, i.e., less than 10 −4 for N > 100. For the investigated design, i.e., the choice of the scaling parameter b, these results indicate that the ∂/∂ρ ln( ) terms do not play a major role and can be neglected for sufficiently fine meshing (N > 100) without significant reduction of accuracy.
Besides the comparison to the analytically computed propagation constant of the linearly polarized mode, a second comparison is presented in Figure 3b. Here, an implementation of the full-vectorial 2D FFT-based mode solver presented in [22] was used to independently compute the propagation constant of the HE 11 mode. As can been from Figure 3b, the effective refractive indices computed with the DHT-based and the FFT-based solvers deviate only by less than 10 −5 for N > 100 if the ∂/∂ρ ln( ) terms are considered in the DHT-based mode solver and, as before, the convergence follows a 1/N 2 trend. Without these terms, the deviation is larger but still less than 10 −4 in all cases. This result is consistent as the FFT-based mode solver accounts for the spatially varying refractive index profile. Similar to Figure 3a, the results indicate that the ∂/∂ρ ln( ) terms only play a minor role for a sufficiently fine mesh. Note that the overall smaller deviation in comparison with the discussed step-index design can be explained by the smoothness of the graded-index profile, i.e., the absence of the Gibbs phenomenon. In Figure 3c, typical computational times required by the DHT-based solver to compute the HE 11 mode of the graded-index design are shown. If the ∂/∂ρ ln( ) terms are included, the computational time is dominated by the computation of the matrix of Equation (45) that is required to perform any ∂/∂ρ operations. If the matrix does not need to be computed, the computational speed is roughly two orders of magnitude faster and does not exceed 10 s for N < 400 (using at each step the previous solution as a guess). It must be emphasized that the applied code implementation computes the matrix of Equation (45) during run-time, which is not required as discussed. Thus, if the computation of the matrix of Equation (45) and the kernel of Equation (24) is performed off-line, it is expected that a computational time less than 10 s for N < 500 can be achieved.

Conclusions
For the first time, a full-vectorial fiber mode solver based on a DHT has been developed and discussed. For rotationally symmetric refractive index profiles, the electrical field amplitudes A ρ (ρ) and A φ (ρ) can globally be expanded as finite Fourier-Bessel series. Subsequently, the corresponding coupled wave equations yield an algebraic eigenwert problem, which can be solved numerically to obtain the expansion coefficients and the propagation constants, i.e., the effective refractive indices, of the guided modes. The DHT allows to conveniently compute the action of the differential operators of the coupled wave equations in the frequency space. The corresponding transformation rules for the Bessel operator and the common derivative operator have been developed and validated here for the first time. The actual implementation of the DHT-based solver was written in Python 3 and was successfully proof-tested for a step-index and a graded-index design as well. For both designs, it has been demonstrated that the effective refractive indices of the guided modes are computed with an absolute error of less than 10 −4 with respect to analytical solutions at a reasonably low number of Fourier-Bessel coefficients. For the graded-index design, a convergence following a 1/N 2 trend has been found. In addition, it has been shown that computational times of a few seconds are feasible, in particular, if transformation kernels and further matrices, which are required to compute the action of differential operators, are computed offline. In the future, the DHT-based solver will be a time and resource-efficient tool to compute the vector modes of arbitrary (rotationally symmetric) fiber designs. For example, it might be used for fast and accurate computation of dispersion profiles. Furthermore, the introduced transformation rules of the differential operators under the DHT might become valuable for other problems that require solving differential equations as well.

Outlook
As an outlook for future work, there are potential improvements but also obstacles of the introduced method, which need to be analyzed critically. For example and beyond common fiber designs, the Gibbs phenomenon will become quite pronounced for high-index contrast designs such as silicon on insulator (SOI) waveguides [28]. Indeed, compensation via an increase in spatial mesh points, i.e., Fourier-Bessel coefficients, might be a viable counter-strategy but will increase the computational time. Similar problems will certainly arise for thin but pronounced refractive index wells and dips for example in Bragg fibers [29] (or other Bragg-type waveguides [30]) or quantum well waveguides [31]. With respect to such designs, another disadvantage of the DHT-based solving approach is the inflexible meshing because the sampling is determined by the zeros of the corresponding Bessel function and cannot be chosen freely. However, a promising improvement to address both issues (Gibbs phenomenon and inflexible meshing) might be a zero-padding as it is common in discrete Fourier transforms. On one hand, such zero-padding could maintain reasonable high mesh resolutions in the space domain. On the other hand, it would keep the number of Fourier-Bessel coefficients reasonably low to maintain acceptable computational times. Another major improvement of the demonstrated DHT-based mode solving approach might be its combination with a common Fourier row expansion (in combinations with discrete/fast Fourier transforms) in the angular coordinate for 2D refractive index designs with higher-order rotational symmetries. Such an approach would overcome the present limitation to rotationally symmetric designs, i.e., specialty fiber designs and concepts such as photonic crystal fibers (PCFs) could also be studied. In that context, a recently developed two-dimensional Fourier transform in polar coordinates [32], which is built upon the introduced DHT, might be helpful. However, some waveguide concepts will remain to be not addressable with the presented approach, in particular, if it is a priori known that evanescent fields do not comply with the DHT-related requirement regarding being (effectively) space-limited as for example for leaky modes in (Bragg-type, anti-resonant, etc.) hollow-core waveguides [33,34].