Wilson basis expansions of electromagnetic wavefields: a suitable framework for fiber optics

To detach large vectorial full-wave electromagnetic scattering problems from the specific, dissimilar bases that are typically associated with the optical waveguiding structures comprising an optical interface, we propose and demonstrate the expansion of wavefields in a single universal Wilson basis. We construct a Wilson basis by following the derivation in the paper by Daubechies, Jaffard and Journé. This basis features exponentially decaying basis functions in both the spatial and spectral domain. The strong localization in phase-space, with basis functions that strongly resemble wavefields, allow for efficient expansions of high-frequency electromagnetic fields. In a Wilson basis, the interface scattering problem is effectively separated from the physical configuration. For the evaluation of multiple, laterally displaced interface configurations, one may reuse electromagnetic fields in the Wilson basis, because the translation operator is sparse and diagonally dominant. We consider actual reflection-transmission problems comprising optical fibers and homogeneous media in the Wilson basis framework in a companion paper. There, the localization in the spectral domain aids the conversion of the numerical scheme to generate electromagnetic fields in homogeneous media due to Wilson-basis discretized electromagnetic sources. In this paper, we review the Wilson basis construction, demonstrate the expansion of modal electromagnetic fields in an optical fiber, and complex-source beams that are tilted with respect to the optical axis. Even for largely tilted beams (up to 60°), despite being highly oscillatory in the cross-sectional plane, the fields are well represented by a finite number of higher-order Wilson basis functions.


Introduction
Passive optical networks for indoor applications typically rely on multi-mode optical fibers to facilitate high-bandwidth communication over distances up to several hundred meters (IEEE 2010). For flexibility in connectivity, the link should be able to support multiple connections, without sacrificing the power budget or transmission distance. The relatively long fiber length is a strong motivation to make an expansion of the physical field (due to an excitation) in terms of the set of guided modes, i.e., the set of source-free solutions to Maxwell's equations that may propagate unattenuated along the fiber. Although the set of guided modes are generally convenient vector fields to work with, for instance to analyze propagation effects such as constructive-and destructive interference patterns (Tsekrekos et al. 2007), or differential mode delay (Smink 2009), they are inherently specific to the waveguide under consideration and the wavelength of choice. In larger electromagnetic problems, for instance reflection-transmission problems at an interface that consists of two dissimilar waveguides at various offsets, it becomes impractical to work with two different bases on either side of the interface. The use of a single universal spatially localized basis would be more convenient, particularly if translations entail sparse operations.
The concept of combining global and local formulations are not unusual. For instance, in the eigencurrent approach to the linear embedding via Greens operators (LEGO) technique (Lancellotti et al. 2009) the computational domain is subdivided into smaller domains call bricks. Each brick encloses a complex structure with equivalent sources expanded in local Rao-Wilton-Glisson (RWG) functions on its boundary. To reduce the number of unknowns of the overall system, the set of eigenfunctions for each brick are spanned by RWG functions. Many of the equivalent-source eigenfunctions do not radiate, resulting in a quite sparse overall linear system of equation containing the interaction of the bricks.
On our agenda are large vectorial full-wave optical problems that consist of dissimilar optical fiber connectors with finite alignment precision and possibily separated by a slab of homogeneous material. For flexibility, a framework employing global-local basis functions may be more suitable to cope with connections of highly dissimilar waveguiding structures, for example consisting of optical fibers and integrated photonic circuits. The ability to tailor the basis functions to the physical size of the optical wavefields and match to the spatial frequencies allows for a high degree of flexibility.
Another example of combining global and local formulations is a beam-based phasespace source representation (Arnold 2002). It offers the localized behavior and flexibility of geometrical optics in the spatial domain, without sacrificing the uniform properties of spectral techniques in the spectral domain. A popular approach to discretize phase-space is the use of Gabor frames, named after Gabor for his idea in 1946 to describe one-dimensional audio signals in terms of two-dimensional elementary functions that contain both the local time and frequency information (Gabor 1946). The expansion functions are obtained by translations and modulations of a generating window function g, which essentially represent a discretization of the windowed spatial Fourier transform, mathematically formulated as g qr ðxÞ ¼ e j2paqx gðx À brÞ; q; r 2 Z; a; b; x 2 R: Throughout this paper, we denote the unit imaginary number by j. The window g, also called Gabor atom, can be moved in phase-space to an arbitrary coordinate on a grid spanned by ðaq; brÞ (Daubechies et al. 1991).
In order to expand an arbitrary function f ðxÞ 2 L 2 R ð Þ through coefficientsc qr by f ðxÞ ¼ X q;r c qr g qr ðxÞ; requires that the collection of functions g qr in Eq. (1) span L 2 ðRÞ. To describe f with a minimum number of coefficients, Gabor proposed efficient sampling by selecting the product ab ¼ 1, and proposed a Gaussian window function g proportional to exp Àpx 2 ð Þ because it achieves good localization in both the spatial and spectral domain (Casazza 2001). Unfortunately, that Gabor scheme is numerically unstable (Daubechies et al. 1991). Actually, the Balian-Low theorem, developed much later, states that it is impossible to achieve good localization in both the spatial and spectral domain simultaneously for these choices (Benedetto et al. 1994). In particular, the Balian-Low theorem states that the function g(x) in Eq. (1) or its Fourier transformĝðnÞ have an infinitely long tail, i.e., either R x 2 jgðxÞj 2 dx ¼ 1 or R n 2 jĝðnÞj 2 dn ¼ 1. The expansion in Eq.
(2) is possible if one chooses ab\1, by using frames that were defined first by Duffin and Schaeffer (1952) in a paper on a deep mathematical problem involving nonharmonic Fourier series (Christensen 2014). By their definition, a set of functions g qr are a frame if there are nonvanishing finite frame bounds A and B, that satisfy Akf k 2 X q;r for all functions f 2 L 2 ðRÞ. The inner products in Eq.
(3) are denoted by h:; :i, and are linear in the second argument. It was due to fundamental mathematical research by Daubechies et al. (1986) and Janssen (1981) as well as others in the eighties that the connection was made between frames, wavelets and Gabor systems (Christensen 2014). In abstract frame theory studies, Gabor frames are often called Weyl-Heisenberg frames (Casazza 2001). The choice ab\1 leads to an oversampling of phase-space, rendering the frame, that is the collection of functions in Eq. (1), linearly dependent (Christensen 2002). In case ab 1, the associated frames are called tight when the frame bounds satisfy A ¼ B (Daubechies et al. 1991). A tight Gabor frame may be redundant, but does not contain redundant functions when ab ¼ 1. For ab ¼ 1, the Balian-Low theorem prevents that the resulting basis has good phase-space localization, apart from the reported numerical instabilities (Daubechies et al. 1991). For a Gabor frame, the expansion in Eq. (2) is not unique (Jing 1999), however an expansion may still be obtained in a practical way with the construction of a canonical dual Gabor frame (Christensen 2014).
Although the redundancy may seem undesirable, Christensen (2014) has argued that the less restrictive frame condition make frames more flexible than orthonormal bases. Frames are not only appealing for theoretical studies, but have also been employed in many applications, including electromagnetics (Dilz and van Beurden 2017).
Fortunately, for large-scale problems where redundancy is undesirable, it turns out to be possible to construct an orthonormal basis that achieves good localization in both the spatial and spectral domain. In the work of Wilson (1987), the construction of bases with strong localization in phase-space entails basis functions with two peaks in the spectral domain. He called these functions generalized Wannier functions. Because of the two peaks, the spectral localization becomes duo-localization, to which the Balian-Low theorem does not apply. Daubechies et al. (1991) were inspired by this idea, and were able to construct an orthonormal basis starting from the same analytical constituents as Gabor frames. This leads to a Wilson basis with trigonometric bases constituents.
There are quite a number of research activities that involve Wilson bases, but these are predominantly of mathematical nature. For instance, studies that focus on its connection with Gabor frames (Bownik et al. 2017;Kaushik and Panwar 2013), construction of Wilson bases with arbitrary shapes (Bittner 1999;Chui and Shi 2000), or studies of discretized Wilson bases (Bolcskei et al. 1996;Kutyniok and Strohmer 2005). So far, the use of Wilson bases in electromagnetics has only been advocated by Arnold (2002).
In this paper we want to promote and demonstrate electromagnetic field expansions in a Wilson basis. The strong localization in phase-space with constituents that have a strong resemblence to wavefields allows for efficient field representations in this basis (Floris and de Hon 2012). Because of the orthonormality, the expansion may be obtained by straightforward projection integrals. In Sect. 3, we show that spatial field translations may be achieved directly in the Wilson basis with a sparse, diagonally dominant operator. This is an appealing feature, and useful for the evaluation of electromagnetically large reflection-transmission problems consisting of optical fibers with different but finite lateral alignment configurations. It means that computational effort to compute the field expansion is needed only once, and reflection-transmission problems may be generalized in terms of Wilson basis functions and coefficients. Due to that separation in two stages, we will dedicate this paper to Wilson basis expansion of electromagnetic waves. In a companion paper (Floris and de Hon 2018), we address reflection-transmssion problems in a Wilson basis. In that paper, the convergence of spectal convolution integrals of dyadic Green's functions for electromagnetic field evaluation due to Wilson-basis discretized source distributions benefit from the spectral localization of the Wilson basis. This allows to extend to interfaces consisting of homogeneous half-spaces or homogeneous slabs. The effective localization in the spectral domain greatly improves the efficiency of the scheme.
In Sect. 2, we shall first give a brief summary of two ways to construct a Wilson basis given by Daubechies et al. (1991), where this subject is treated in great mathematical detail with abundant proofs. We show the construction through evaluation of an inverse Zak transform, as well as through Daubechies' algorithm derived using the connection of Wilson bases and tight frames. Then, in Sect. 4, we construct Cartesian two-dimensional Wilson basis functions, and describe two methods to expand electromagnetic fields on planes in these 2D Wilson basis functions. First, we discuss a method that was proposed by Sullivan et al. (1989), that makes use of a Taylor-series expansion of the function to be expanded around the spatial coordinate where the basis function is also conveniently localized. It turns out that the resulting integrals are readily evaluated in the spectral domain. Although the method is particularly convenient for generalized Wannier functions, we show that the same approach is still valid for Wilson bases as proposed by Daubechies et al. (1991). For practical reasons, we also review the construction of numerical cubatures to evaluate inner products of two-dimensional Cartesian Wilson basis functions and fields that are represented in cylindrical polar coordinates.
In Sect. 5, we briefly summarize the construction of complex-source beam (CSB) wavefields, which are electromagnetic fields that satisfies Maxwell's equations everywhere (except at a branch cut disc) and is constructed by inserting a complex source coordinate in the dyadic counterpart of the scalar Green's function. A CSB is easily evaluated, convenient to test interface problems consisting of homogeneous media, and is a Maxwellcompliant model of a laser illuminating an optical fiber in the companion paper (Floris and de Hon 2018).
In Sects. 6 and 7, we give examples of complex-source beams and modal electromagnetic fields in the Wilson basis. First, we consider three CSBs illuminating a crosssectional plane at three launch angles h z 2 ð0 ; 40 ; 60 Þ at a distance equal to the radius of the source disc, so that the CSBs are diverging. The examples demonstrate that the highfrequency oscillations due to large angles of incidence are still well represented using a limited number of higher-order Wilson-basis functions.
The Wilson-basis expansion of modal electromagnetic fields involve multiple cubature rings to account for the discontinuities in the refractive index profile. By means of an example, we demonstrate the power distribution of an arbirary guided modal field in Wilson-basis coefficient space. It shows that including only a few levels of higher-order basis functions are already sufficient to accurately represent modal electromagnetic fields in the Wilson basis.

Wilson basis construction
In order to efficiently approach large scale problems, Sullivan et al. (1989) argue that calculations should involve orthogonal, complete, generic, local bases with useful analytic properties. The locality should apply to both position and spatial frequency (momentum in their case), to minimize the basis set size. The generalized Wannier functions introduced by Wilson (1987) satisfy these requirements. These basis functions are of the form We regard x as a spatial variable. The associated spectral functionsf q ðnÞ are described bŷ and have two peaks along the spectral coordinate n, located at the spatial frequencies Àq=2 and q=2, with/ þ q and/ À q centered around zero. This amounts to an effective duo-localization, to which the Balian-Low theorem does not apply, because that is tied to the specific structure of g qr in Eq. (1).
The formulation of Wilson bases by Daubechies et al. (1991) is of a similar form as Eq. (4), however is more elegant in the sense that it entails the construction of only one single symmetric, real-valued function h, rather than many different functions f q for every q 2 N. Together with translations and modulations by means of trigonometric functions, the Wilson basis functions read w 'n ðxÞ ¼ hðx À n=2Þ; ' ¼ 0; n even; 0; ' ¼ 0; n odd; ffiffi ffi 2 p hðx À n=2Þ cosð2p'xÞ; '[ 0; ' þ n even; ffiffi ffi 2 p hðx À n=2Þ sinð2p'xÞ; '[ 0; ' þ n odd: In the original formulation in Daubechies et al. (1991), the ' ¼ 0 case was defined by hðx À nÞ for all n. However, to achieve consistent displacement of the function h due to the parameter n for all values of ', we have decided to change the indexing accordingly. The Wilson basis is orthonormal because the inner product of two basis functions satisfies Z where Ã denotes complex conjugation. With the aid of the (forward) Fourier transformation defined byŵ we obtain the representation of the spectral counterpart of the Wilson basis functions in Eq. (6), which features the spectral duo-localization bŷ w 'n ðnÞ ¼ e jpnnĥ ðnÞ; ' ¼ 0; n even; Now we only need to numerically construct the function h, and its Fourier transform counterpartĥ. To achieve that, we adopted the procedure from the paper by Daubechies et al. (1991), and give a short summary here. In that paper, it is first shown, in great detail, that there is an equivalence between demanding orthonormality for all generalized Wannier functions w qr in Eq. (4), and the frame condition in Eq. (3). More specifically, to achieve orthonormality among all functions w qr with q 2 Z [ 0 ; r 2 Z leads to infinitely many requirements on the shape of the single functionĥ. These conveniently reduce to only one condition when transformed to the Zak domain. That unitary transformation to the Zak domain is achieved by the Zak transform U Z , which is defined in Daubechies et al. (1991) as where in our caseĜðt; sÞ ¼ ðU Zĝ Þ is the Zak-transform of a Gaussian window function. In particular, we define a Gaussian window functionĝ m ðnÞ in the spectral domain that may be scaled along the spectral coordinate by a parameter m 2 R aŝ The peculiar amplitude ð2mÞ 1=4 ensures that the inverse Fourier transform g m ðxÞ is conveniently related toĝ m ðnÞ through g m ðxÞ ¼ĝ1 m ðxÞ. The functionĥ in Eq. (9) may be obtained through the inverse Zak transformation of a functionĤ, defined by Daubechies et al. (1991) The specific fraction in the definition ofĤðt; sÞ ensures that the Zak-equivalent condition jĤðt; sÞj 2 þ jĤðt; s þ 1 2 Þj 2 ¼ 2 associated with the orthonormality requirement (Daubechies et al. 1991), is satisfied. Indeed, this is easily verified through direct substitution ofĤðt; sÞ from Eq. (12) into this condition, and the use of Eq. (10) to rewriteĜðt; s þ 1Þ ¼ e j2ptĜ ðt; sÞ by letting l ! l À 1.
For a window function g m ðnÞ with m ¼ 1= ffiffi ffi 2 p , we evaluatedĥðnÞ numerically using an N-point trapezoidal rule, and tested the accuracy using an adaptive Gauss-Krondrod rule to construct a reference functionĥ ref ðnÞ for 1001 coordinates in n 2 0; 8 ½ . We configured the Gauss-Kronrod quadrature to achieve accuracy down to the order 10 À14 . The shape ofĥðnÞ is shown in the lower two graphs of Fig. 1. The absolute difference between the two integration methods are shown in both graphs of Fig. 2.
The graphs on the right show that it suffices to limit the sum over l in the Zak-transform in Eq. (10) to bse À l 2 Àl max ; l max ½ with l max ¼ 2 to achieve machine precision for n 2 À8; 8 ½ .
Here, bse denotes the integer nearest to s. The graph on the left shows that the absolute difference reduces two orders of magnitude each time the trapezoidal scheme employs two additional samples.  6) depend on translations and modulations of the function hðxÞ, that is shown in the top two graphs. Its spectral counterpartĥðnÞ is shown in the bottom two graphs. The graphs on the right show that these functions have exponential decay in both the spatial and spectral domain This example for the construction ofĥðnÞ demonstrates fast convergence for a short and simple numerical code. Efficient evaluation of spectral Wilson basis functions is highly valuable, for instance to achieve fast convergence of convolution integrals of spectral dyadic Green's functions for the evaluation of electromagnetic fields due to source distributions conveniently discretized in localized Wilson basis functions, especially when short distances are concerned. We do however also require that the inverse Fourier transformed counterpart hðxÞ is available to expand electromagnetic fields in the spatial domain.
Although the inverse Fourier transform hðxÞ may also be computed in a similar fashion as Eq. (12), it requires more work in the sense that it comprises the evaluation of four inverse Zak-transforms of scaled and modulated versions of the functionĤ (Daubechies et al. 1991).
To avoid that, we continue with an alternative construction of the functions hðxÞ and hðnÞ, using the algorithm developed by Daubechies et al. (1991). They observed that the specific construction of a Wilson basis is linked to tight frames (A ¼ B in Eq. 3), in the sense that skillfully choosing specific linear combinations of frame functions removes the redundancy entirely. To paraphrase Daubechies et al. (1991): ''it is a small miracle that the result is an orthonormal basis!''. Although it is claimed that the algorithm is short and straightforwardly programmed (Daubechies et al. 1991), we dedicate the appendix to give a short and streamlined summary of the derivation of the algorithm. The appeal of using their algorithm is that the evaluation of the functions hðxÞ andĥðnÞ at their respective arguments entail the evaluation of a weighted sum of frame functions with coefficients a qr that have Gaussian decay. The functionĥðnÞ is constructed througĥ a qrĝqr ðnÞ; withĝ qr ðnÞ ¼e jpqnĝ m ðn À rÞ; a qr ¼ The graph on the right shows that it suffices to limit l max to two to achieve machine precision for n 2 À8; 8 ½ . The left graph shows rapid decay of the absolute difference for increasing N and the inverse Fourier transformed function hðxÞ has a similar construction, i.e., qr a qr g qr ðxÞ; with g qr ðxÞ¼ðÀ1Þ qr e Àj2prxĝ 1 m ðx À q=2Þ: The coefficients a qr are obtained with the aid of recurrence relations for f k and b k qr that are respectively given in Eqs. (40) and (42) in the appendix. The coefficients a qr need to be determined only once. Due to the similarity of construction, the evaluation of spatial and spectral basis functions require approximately an equal number of numerical operations to achieve an equal level of accuracy. Moreover, due to the closed form representation, certain mathematical operations such as inner products between displaced Wilson basis functions may be evaluated analytically, as we will show in Sect. 3. Now we are left to determine the coefficients A and B, which depend on the choice of the parameter m in the Gaussian function. With the aid of the Zak transform in Eq. (10), the frame bounds A and B are then found by calculating In view of the Gaussian decay of the windowĝ m , and the small domain of s and t, the infinite sum in Eq. (10) can be reduced to incorporate only a few terms around l ¼ 0. We have selected m ¼ 1= ffiffi ffi 2 p for the Gaussian window, and obtained the bounds A ¼ 1:529887182 and B ¼ 2:491627873. In Fig. 1, the function hðxÞ (top) as well as its Fourier transformĥðnÞ (bottom) are shown, both on a linear and logarithmic scale. The figures on the right clearly visualize the exponential decay in both domains.
We have truncated the sum over q and r to AE31 in Eq. (45), which is justified by the Gaussian decrease in q and r as seen in Eq. (44) in the appendix. To visualize the decay of The absolute value of the coefficients ja qr j demonstrating exponential decay in q. The value for r associated with each line is given in the legend. The coefficients a qr are defined in Eq. (45) and are used to construct h andĥ in Fig. 1 the coefficients a qr along the q and r index, we have plotted ja qr j as function of q for increasing indices r in Fig. 3. The evaluation of the coefficients a qr depend on the sum over k in Eq. (13) for a qr . The maximum absolute error of the functions h andĥ due to the truncation in parameter k decays as expðÀ1:6kÞ, so that summation up to k ¼ 16 achieves an accuracy of the order 10 À12 . We have subsequently truncated the set a qr and only take the 120 dominant coefficients into account for the construction of h andĥ.
These graphs for hðxÞ andĥðnÞ for m ¼ 1= ffiffi ffi 2 p in Fig. 1 are different from those presented by Daubechies et al. (1991), although the algorithm and the coefficients for A and B are the same. We have strong reasons to suspect that the graphs presented in the otherwise excellent paper by Daubechies et al. (1991) are wrong (Arnold 2011).
In view of reflection-transmission problems, we want to expand electromagnetic fields in terms of Wilson basis functions at two-dimensional interface planes that are oriented perpendicular to the z-axis.

Field displacement in a Wilson basis
To address electromagnetic reflection-transmission problems for an interface that consist of two spatially displaced waveguides in a universal Wilson basis, we want to be able to displace electromagnetic fields in a Wilson basis as well. To achieve that, we need a translation operator, that maps all Wilson basis function to the set of displaced Wilson basis functions.
For example, it may readily be recognized from Eq. (6) that the operation n ! n À 2 for all basis functions essentially leads to w 'n ðxÞ ! w 'n ðx þ 1Þ which facilitates a displacement of the expanded field to x ! x À 1 along the spatial coordinate. So by shifting Wilson basis coefficients in coefficient space already achieves integer spatial displacements.
To achieve an arbitrary displacements so that x ! x À x D , we perform an expansion of displaced Wilson basis function in terms of non-displaced Wilson basis functions. Fortunately, because of the localization of the basis functions, these operations are sparse, in the sense that only a few inner products need to be evaluated. It is convenient to compute these inner products in the spectral domain, Depending on the parity of ' þ n and ' 0 þ n 0 associated the Wilson basis functionsw 'n ðnÞ andw ' 0 n 0 ðnÞ defined in Eq. (9), the integral in Eq. (16) is equal to the weighted sum of one, two or four integrals of the form For instance, with ' þ n even and ' 0 þ n 0 odd, Eq. (16) equals Upon replacingĥðn þ 'Þ and its complex conjugate by the expressions in Eq. (13) and rearranging terms, we infer that the integral I may be written as a Ã qr a q 0 r 0e jpð' 0 q 0 À'qÞ hĝ qþn;rÀ' ;ĝ q 0 þn 0 À2xd;r 0 À' 0i; where the inner products of shifted and modulated frame functionsĝ are readily evaluated in Eq. (44).
We have evaluated the integrals in Eq. (16), and visualize the power density distributions due to the projections of displaced Wilson basis functions w ';n ðx À x d Þ onto the nondisplaced counterparts w ' 0 ;n 0 ðxÞ in Figs. 4 and 5. From left to right, the subfigures consider the expansion of w ';n ðx À x d Þ for increasing levels of ' from zero to two, and displacements x d 2 ð0; 1Þ. For convenience, we display the expansions of basis functions with n ¼ 0, however these results generalize all cases of n even, although then n 0 ! n 0 þ n. In Fig. 4, the power density is grouped by ' 0 , and each of the six levels for ' 0 is labeled and highlighted by a different color.
By orthonormality, at x d ¼ 0, the basis functions couple only to the same non-displaced counterparts. For displacements 0 \ x d \ 1, the coupling is strongest among basis functions with a level ' 0 equal to '. The power coupling rapidly decreases as the difference j' À ' 0 j becomes larger. For instance, in the middle and rightmost subfigures, the coupling is dominant among functions with ' 0 ¼ ', followed by basis functions of the two Fig. 4 From left to right, the Wilson basis functions w '0 ðx À x d Þ for ' ¼ 0; 1; 2 are displaced by x d and projected onto the non-displaced counterparts w ' 0 ;n 0 ðxÞ. Grouped by ' 0 , the strongest coupling occurs among basis functions that satisfy ' 0 ¼ ', and have lower coupling when the mismatch in j' 0 À 'j becomes larger. (Color figure online) neighboring levels ' 0 ¼ ' À 1 and ' 0 ¼ ' þ 1. As x d ! 1, the power coupling occurs only to basis functions with ' 0 ¼ '.
In Fig. 5, we visualize the same information, but now grouped by n 0 in each subfigure. Again, by orthonormality, at x d ¼ 0, the basis functions couple only to the same nondisplaced counterparts as highlighted through the (overall descending) red curve. For 0 \ x d \ 1, the coupling is dominant to only a few neighboring functions, and are weaker to basis functions that are spatially farther away. Particularly, we highlight the direct spatial neighbor n 0 ¼ 1 through the (symmetric) black curve. In the left case, for ' ¼ 0, the contributions to the n 0 ¼ 1 functions are those that satisfy ' 0 [ 0, because there are no basis functions tied to the ' 0 ¼ 0 and n 0 odd case. As a result, the contributions to the n 0 ¼ 1 functions is more dominant in x d in the subsequent graphs with ' [ 0. As expected, for x d ¼ 1, the basis function w 'n ðx À 1Þ couples entirely to w ';nþ2 ðxÞ as highlighted through the (overall ascending) green curves.
This demonstrates that an arbitrary spatial displacement of an electromagnetic field in a Wilson basis may be achieved through an operator that is sparse and dominant along the diagonal.
The numerical accuracy of the evaluation of inner products among Wilson basis functions depends on the number of coefficients a q;r and a q 0 ;r 0 included in Eq. (19). To supplement the absolute accuracy of the construction of Wilon basis functions (cf. Fig. 3), the truncation effect on the inner products of basis function w 'n and w ' 0 n 0 in Eq. (16) for x d ¼ 0 are shown in Fig. 6.
The solid green (most rapidly descending) line visualizes the rapid decay of j1 À hw 0n ðxÞ; w 0n ðxÞij upon including an increasing number of coefficients a qr and a q 0 r 0 . It reaches 10 À12 upon incorporating at least 60 coefficients in both sums in Eq. (19). The dashed orange line represents the inner product of two neighboring basis functions hw 0n ðxÞ; w 0;nþ2 ðxÞi and experiences a similar decay. The (thin) grey curves show the decay of all inner products hw 0n ðxÞ; w 0;nþ2p ðxÞi with p 2 Z n f0g. The larger the value of p, the slower the decay in the number of coefficients, although its maximum decreases rapidly as well. As an example, highlighted through the black dashed line, the inner product hw 0n ðxÞ; w 0;12 ðxÞi always keeps below 10 À6 , but has a long tail that slowly decays. The envelope of all the curves is highlighted through the solid red line, and essentially quantifies the number of coefficients that should be included in the construction of the Wilson basis to achieve a desired level of accuracy in the inner products. Inner products of combinations of Wilson basis functions with ' [ 0 and ' 0 [ 0 do not change the envelope, and for that reason we did not include them in the figure.
Numerical evaluation of integrals with Wilson basis functions are relevant for the electromagnetic field expansions that we discuss in the next section, and demands finite integration bounds. To assess the impact of finite integration bounds on the accuracy of Wilson basis inner products hw 'n ðxÞ; w ' 0 n 0 ðxÞi, we have employed an adaptive numerical integration scheme with an absolute tolerances set to 10 À10 . Furthermore, we chose x 2 ½Àx max ; x max .
The graph with black cross marks (Â) shows that the accuracy of Wilson basis inner products is bound to 10 À4 because of the finite integration bound set to x max ¼ 4. The graph consists of four overlapping lines, and each line represents the worst offending deviation from orthonormality for basis functions with ' ¼ 0, and ' 0 2 ð0; 3Þ. Before the lines saturate when more than 40 coefficients are included in the construction of the Wilson basis functions, the decay nicely follows the red envelope curve. Upon extending the integration bounds to x max ¼ 8, denoted through the lines with blue plus marks (þ), the accuracy of 10 À8 is achieved upon including upward of a 100 coefficients. x max = 8 Fig. 6 The accuracy of the inner products \ w 0;n ; w 0;nþ2p [ evaluated in the spectral domain for p 2 N unequal to zero are plotted in grey as function of the number of coefficients a qr included in Eq. (19). The solid green curve represents j1 À\ w 0;n ; w 0;n [j, and the orange dashed curve for \ w 0;0 ; w 0;2 [ has similar decay, whereas \ w 0;0 ; w 0;12 [ , highlighted by the black dashed curve is lower but has a longer tail. The envelope is highlighted in red. For inner products evaluated numerically in the spatial domain with finite integration bounds x max ¼ 4, shown with black crosses, the accuracy is bound to 10 À4 , whereas for x max ¼ 8, marked by blue plusses, 10 À8 may be achieved. (Color figure online)

Field expansion in a Wilson basis
In this section, we show two possible ways to expand a function in a Wilson basis. Because we work with basis functions that are orthonormal with respect to the norm in Eq. (7), the coefficients c 'n that describe a function hðxÞ ¼ P ';n c 'n w 'n ðxÞ follow from In fiber optics, we are primarily interested in fields on 2-dimensional cross-sectional planes in cylindrical polar coordinates. Before we construct a numerical integration procedure to evaluate the integrals for these problems, we first review the approach as proposed by Sullivan et al. (1989) for a similar expansion as Eq. (20) but in terms of generalized Wannier functions. They recognized that integrals of the product of monomials and basis functions are readily evaluated in terms of derivatives of the spectral basis function at the zero coordinate (Sullivan et al. 1989). A Taylor-series expansion of the function h in Eq. (20) about the coordinate where the basis function peaks, leads to a sum of integrals of this form. Particularly for generalized Wannier basis function, many of these integrals subsequently vanish because of symmetry arguments as well as due to the localization property of the basis functions (Sullivan et al. 1989).
To assess whether the approach in Sullivan et al. (1989) is useful in for our case, we follow a very similar procedure as outlined in the appendix of Sullivan et al. (1989) but now for Wilson basis functions. The Wilson basis functions in Eq. (6) have a shape that depends on the parity of the sum ' þ n, whereas the generalized Wannier functions satisfy w mn ðxÞ ¼ w m0 ðx À nÞ consisting of a set of (infinitely many) functions f m as shown in Eq. (4). Hence, the equations are different, but the approach is the same.
First, we replace the function h in Eq. (20) by its Taylor-series expansion about the spatial coordinate x ¼ n=2 (where the Wilson basis functions w 'n are also localized), and obtain The integral over a polynomial times a Wilson basis function is readily evaluated in terms of derivatives of the spectral basis functions at the origin. Specifically, consider the k thderivative of an expðÀjpnnÞ phase-shifted spectral Wilson basis function in terms of the Fourier integral of Eq. (8). We obtain which on the right-hand side includes the integral we wanted to evaluate in Eq. (21). The phase-term expðÀjpnnÞ on the left-hand side conveniently cancels the phase-term that is present inŵ 'n in Eq. (9). Combining Eqs. (21) and (22), we infer that the coefficients c 'n are readily evaluated in the spectral domain by In this form, only the functionĥ depends on n. The derivatives of this function may subsequently be evaluated by the sum in Eq. (45) albeit upon replacing g qr by its k th derivative with respect to the argument. In contrast to the case of generalized Wannier functions in Sullivan et al. (1989), we can not identify terms that conveniently vanish for most k in Eq. (23). Moreover, this method demands that the k th derivative of the function h with respect to its argument can be determined accurately, and preferably efficiently.
With the numerical scheme for the Wilson basis available, and given that electromagnetic fields are generally obtained numerically as well, we resort to numerical integration to determine the Wilson basis coefficients. We proceed by constructing a numerical counterpart for the evaluation of Eq. (20) using numerical quadratures but for two-dimensional functions. In view of interfaces that consist of two dissimilar waveguides at each side of an interface, displacements with respect to the center of the interface are conveniently specificied in Cartesian coordinates. In a fashion similar to Eq. (20), we expand a two-dimensional function h(x, y) in two-dimensional basis functions that are simply the product of two basis functions associated with the x-and y-coordinates. As a result, each two-dimensional basis function has four indices, which we abbreviate for convenience in terms of the multi-index i ¼ ð' x ; n x ; ' y ; n y Þ, so that w i ðx; yÞ ¼ w 'xnx ðxÞw 'yny ðyÞ. The coefficients due to the expansion of h(x, y) are obtained by the integral where r T ¼ xu x þ yu y and dr T is short for dxdy. In order to match the width of the basis functions more appropriately to the electromagnetic fields, we introduce a scaling parameter d. To preserve the orthonormality property, the scaling also modifies the amplitude of the basis functions, i.e., There are a large number of guided modal electromagnetic fields associated with circularly symmetric optical fibers that we want to expand in the Wilson basis. For example, a standard multi-mode fiber at a wavelength of 850 nm supports about 180 guided modal fields (Floris et al. 2014). The computation of the modal fields in cylindrical coordinates comprises the numerical computation of fields along the radial coordinate, and the evaluation of an analytical Fourier expansion along the angular coordinate (Smink 2009). Depending on the shape of the refractive index profile of the fibers, electromagnetic field components (or their derivatives) may be discontinuous, for instance at the core-cladding interface. Therefore, we require that the numerical integration along the radial coordinate can be subdivided into appropriate subintervals.
To expand electromagnetic fields that are conveniently described in cylindrical polar coordinates in terms of Cartesian Wilson basis functions, one can of course always employ two nested adaptive quadratures, and make use of substitutions. However, in view of the large number of Wilson basis coefficients we wish to evaluate, that would become very time consuming.
We use cubatures, i.e., quadratures for more than one dimension, that are specifically constructed to numerically integrate over an annular region using a fixed cylindrical grid of abscissae of the form Here, the vector hðq; uÞ ¼ fðq; uÞw i ðq cos u; q sin uÞ contains the electromagnetic field vector f, and the Wilson basis function w i . The weights D uv and abscissae q u and u v are known in literature, or can be constructed in such a way that Eq. (26) is solved exactly for functions h that are polynomial in x and y, with x ¼ q cos u and y ¼ q sin u, up to a desired total degree (Peirce 1957). With those substitutions, the integration of polynomials in x and y over an annular region may be separated into two quadratures. Along the angular coordinate, a quadrature with equidistant abscissae with constant weights (achieved through a Chebyshev-Gauss quadrature, or a trapezoidal rule) suffices. Along the radial coordinate, a quadrature that accounts for the Jacobian q and integration bounds q 0 and q 1 may be constructed using, for example, a Golub-Welsch algorithm (Gil et al. 2007). An overview of known cubatures over a unit disc (i.e., for q 0 ¼ 0 and q 1 ¼ 1) with weights and abscissae to achieve any set degree of polynomials with as few coefficients as possible are listed in Cools and Kim (2000). For our purpose, it suffices to employ a Gauss-Legendre quadrature along the radial coordinate, and incorporate the Jacobian q directly. Neither electromagnetic fields nor the Wilson basis functions are polynomials, hence we need to ensure we take sufficient coefficients in both quadratures.

Complex source beam
In this section, we briefly describe the construction of electromagnetic complex source beam wavefields in homogeneous space. These fields may be mathematically generated by inserting a complex-valued spatial source coordinate in the dyadic counterpart of the scalar Green's function Felsen 1986, 1989). The method was first introduced in the scalar case by Deschamps (1971) in the context of complex rays, by the realization that ray formulations with complex amplitude-and phase functions still satisfy the eikonal equation. The physical interpretation of these complex rays correspond to Gaussian beams (Deschamps 1972).
The electromagnetic fields due to source distributions may be formalized in terms of a volume integral over a dyadic Green's function operating on a source vector distribution. In case of a combination of electric and magnetic dipole sources, the volume integral reduces to The vectors J u and K u denote the electric and magnetic dipole unit vectors respectively, and J d and K d the respective amplitudes. G denotes the dyadic Green's function comprising the electromagnetic field response (e and m subscripts) due to electric and magnetic dipoles (e and m subscripts). We assume an exp jxt ð Þ dependence on time t, where x denotes the angular frequency. For a homogeneous medium with permittivity e and permeability l equal to the free-space permeability l 0 , the refractive index n satisfies n 2 ¼ e=e 0 , with e 0 the free-space permittivity. For such media, the four dyads in Eq. (27) may be derived from Maxwell's equations (e.g., see Felsen and Marcuvitz 1994). To briefly show that, consider the inhomogeneous dyadic wave equation, where the wavenumber k ¼ k 0 n 2 contains the free-space wavenumber k 0 ¼ 2p=k related to the free-space wavelength k. I is the identity dyad, and r and r 0 are the observation and source coordinates. Eq. (28) may be rewritten to the dyadic equivalent of Helmholtz equation, with the aid of the dyadic identity r Â r Â G ee ¼ rr Á G ee À r 2 G ee . Hence the term rr Á G ee may be obtained from the dyadic counterpart of the Maxwell-Ampère equation, and satisfies rr Á G ee ¼ À jxe ð Þ À1 rr Á Idðr À r 0 Þ. Combining these last two equations allows to rewrite Eq. (28) as, Using the identity Àdðr À r 0 Þ ¼ r 2 þ k 2 ð Þ Gðr; r 0 Þ, with Gðr; r 0 Þ representing the scalar Green's function G ¼ ð4pjr À r 0 jÞ À1 expðÀjkjr À r 0 jÞ; it is easily seen that the dyadic Green's function, generated by electric dipoles of arbitrary orientation, is related to G by The dyadic representation of the associated magnetic field is obtained through the Maxwell-Faraday equation and reads Similarly, a magnetic dipole source leads to G mm ¼ Àjxe I þ k À2 rr ð Þ G and G em ¼ ÀrG Â I.
Insertion of a complex source coordinate r 0 ¼ jb in the argument of the Green's function, with b ¼ bu b ; b 2 R, leads to a complex source beam (CSB) in the beam direction u b . Care must be taken in the evaluation of the associated electromagnetic fields, because of the branch cut associated with r À r 0 j j in Eq. (30). We keep to the so-called source region for all the observation coordinates r by ensuring that Re r À r 0 f g[ 0 (Heyman and Felsen 1989). We are free to choose the electromagnetic dipole amplitude and orientations. We follow Cullen and Yu (1979), and let the electric and magnetic dipole vectors be oriented perpendicular with respect to each other, and with respect to the beam direction u b . Furthermore, we choose the amplitude of the electric and magnetic source amplitude vectors in such a way that J d ¼ ÀZ 0 K d , where Z 0 ¼ ffiffiffiffiffiffiffiffiffiffiffi l 0 =e 0 p denotes the freespace impedance. This achieves almost equal on-axis contribution from the J and K dipole sources. The beam has a Gaussian-like shape. We now proceed with the expansion of this field in a Wilson basis.

Expansion of a complex-source beam in a Wilson basis
We want to expand complex-source beams (CSB) in a Wilson basis for several reasons. The Gaussian beam-shaped electromagnetic field is a very convenient representation for an optical light source, that may be used to illuminate an optical fiber. This reflectiontransmission problem can be solved upon constructing scattered electromagnetic fields that have equivalent electromagnetic source distributions located at the interface, in such a way that the sum of the incident and scattered electromagnetic fields match the transmitted field in the fiber at the interface. CSBs are particularly useful to test our implementation for the electromagnetic field evaluation at an interface z ¼ z 1 due to equivalent sources at another interface z ¼ z 0 . The reason is that CSBs are easily evaluated at any interface including z ¼ z 1 , as long as we keep away from the branch cut. We demonstrate the approach for solving electromagnetic reflection-transmission problems in a Wilson basis in a companion paper (Floris and de Hon 2018).
We have expanded three complex source beams in a Wilson basis with the aid of numerical cubatures of Eq. (26). We constructed one complex source beam with a beam vector b ¼ b z u z using b z ¼ 70 lm at a wavelength of k ¼ 1:55 lm in free-space, and subsequented tilted about the x-axis in the positive y-direction by h z ¼ 0 , h z ¼ 40 and h z ¼ 60 .
For the expansion in the plane z ¼ 70 lm, we used a cubature disc with bounds q 0 ¼ 0 and q 1 ¼ 100 lm, using 400 samples along the radial coordinate and 1600 samples along the angular coordinate u 2 0; 2p ½ . At this distance, the electromagnetic field is diverging, yet not very much in that the total power of the electromagnetic field captured within the cubature disc is equal up to at least 6 digits, for all three tilted beams. We did however displace the cubature in the y-direction to the coordinate y 0 ¼ 70 tanðh y Þ lm, to account for the beam displacement. We have expanded the field in basis functions with l x ; l y 2 ð0; 6Þ and n x ; n y 2 ðÀ40; 40Þ. This leads to a grand total of 321489 Wilson basis coefficients. We are only interested in the coefficients that significantly contribute to the expansion. By orthonormality of the Wilson basis, the power through the cross-sectional plane evaluated by the integral over the z-directed Poyning vector may be determined directly in the Wilson basis. Fortunately, we only need to consider a small subset of coefficients as shown in Fig. 7. For the three tilded CSBs, the deviation from the total power of the CSB that is incident on the cubature disc is shown as function of increasing number of expansion coefficients.
For normal incidence, a power accuracy of 10 À5 is achieved through less than 500 coefficients. The larger the angle of incidence, the more coefficients are necessary because those CSBs are more spread out and contain higher spatial frequencies.
For instance, in blue, at the top of Fig. 8, we have plotted the real part of the highly oscillatory x-component of the electric field of a CSB that is tilted along the y-axis with h z ¼ 60 , along the line x ¼ 0, z ¼ 70 lm. It shows that the beam field is asymmetric and diverging, whereas in red, for an observation at a cross-sectional plane at z ¼ 5 lm, the beam is highly collimated. However, at this distance, the tilted source disc crosses the observation plane, causing a sudden cut-off of the electromagnetic fields. This is to demonstrate that one must remember to account for the source disc of a CSB. At the bottom in Fig. 8, we have visualized the size and shape of the pertaining Wilson basis  Fig. 8 In blue, in the top graph, the highly oscillatory electric field for a beam tilted by h z ¼ 60 is plotted along the line x ¼ 0, z ¼ 70 lm. In red, along the line x ¼ 0, z ¼ 5 lm, the field is highly collimated, however is suddenly cut-off because the observation plane crosses the tilted source disc. In the lower graph, functions that are scaled to be commensurate with the size of the CSB through Eq. (25) with d ¼ 5:47 lm.
Judging from the frequency of the oscillations, it may be noticed that it suffices to include basis functions up to the level ' ¼ 3. As evidence for that, we have plotted the power distribution for the entire CSB in the Wilson basis for different combinations of ' x and ' y in Fig. 9. Without paying attention to the specific lines in each of the four figures yet, it may already be noticed that the electromagnetic fields do not heavily oscillate along the x-coordinate, because the power levels are the highest in the leftmost figure for ' x ¼ 0 and decrease consistenly for increasing levels of ' x . Actually, in case h z ¼ 0, denoted by the blue line, the decay of the power distribution is the same for ' x as it is for ' y , due to the comparable oscillations of the electromagnetic field along the x-and y coordinates. As the angle h z is increased to 40°and 60°, the power distribution changes to higher levels in ' y .
We have demonstrated the expansion of complex source beams in a Wilson basis. In view of reflection-transmission problems for interfaces involving optical fibers, we proceed with an example of the Wilson basis expansion of guided modal electromagnetic fields.

Expansion of modal electromagnetic fields in a Wilson basis
We mentioned before that circularly symmetric optical fibers may have discontinuities in the refractive index profile along the radial coordinate. For instance, step-index optical fibers, or trench-assisted multi-mode optical fibers have discontinuities, or a step in the derivative, for instance at the core-cladding interface. To cope with those, we employ multiple cubatures that integrate over annular regions.
As an example, we show an electric field component of a linearly polarized LP 3;1 mode in Fig. 10 that is associated with a trench-assisted space-division multiplexing (SDM) fiber with a profile that is shown in Sillard et al. (2014). In this example, we have employed a scalar approximation for the electromagnetic fields, which is common in the field of spacedivision multiplexing, and facilitates a labeling of mode groups. Each mode group consists of a finite number of guided modes that have similar optical properties, such as phase velocities. The number of mode groups are a measure for the number of independent transmission channels that may be formed in the SDM fiber. Fig. 9 Grouped by ' x the power distribution of three entire complex-source beams are plotted as function of ' y . For normal incidence (h z ¼ 0), the power distribution in ' x and ' y are similar, whereas for larger angles of incidence, the power distributions of the y-directed electromagnetic fields change to higher levels in ' y to represent the higher spatial frequencies observed at the cross-sectional plane In Fig. 10, we have plotted the real part of the x-component of the electric field of a relatively high-frequency LP 3;1 modal electromagnetic field, that is associated with the fourth mode group.
On the right, we have plotted the power distribution in Wilson-basis coefficient space. We limit the expansions of the modal electromagnetic fields to ' x ; ' y 2 ð0; 3Þ. By grouping the power distribution to integer levels of ' x þ ' y , we observe that low-order Wilson basis functions dominantly contribute to the field expansion, especially for n x small. At each coordinate n x , the sum of the power over all n y is included. As the level ' x þ ' y increases, the density of accumulated power decreases. The combinations for ' x þ ' y [ 3 stay below 10 À6 and are not added to the graph. As shown in Fig. 7, upon incorporating about 600 dominant coefficients, unit power is achieved to an accuracy of 10 À5 .

Discussion
We have demonstrated Wilson-basis expansions of electromagnetic fields, in preparation to solving electromagnetically large reflection-transmission problems in a more general setting in a companion paper (Floris and de Hon 2018). The characterization of electromagnetic fields through coefficients in a Wilson basis, which in turn optimally samples phase-space with constituents that share similarity to Gaussian wavefields (Arnold 2002), effectively detaches the electromagnetic scattering problem from the specific optical structures or devices present at the interface. The global-local Wilson basis functions thus allow for efficient expansions of high-frequency electromagnetic fields (cf. Fig. 7), with coefficients that are obtained through straightforward projection integrals. Furthermore, a spatial translation of an electromagnetic field in the Wilson basis, for instance due to a reconfiguration of an interface problem, comprises a sparse, diagonally dominant translation operator (cf. Figs. 4 and 5). Although the expansion needs to be computed first to a desired degree of accuracy, the flexibility to reconfigure an interface by reusing l x + l y = 0 l x + l y = 1 l x + l y = 2 l x + l y = 3 Fig. 10 The power density in the Wilson-basis coefficients is visualized through the power accumulation for each n x -coordinate for integer levels of ' x þ ' y (shown right) that is associated with the Re E x ðx; yÞ f g field of the LP 3;1 mode (shown left) computational effort, and the generalization of scattering problems in the Wilson basis that come in return are rewarding.
The construction of the Wilson basis in the paper by Daubechies et al. (1991) has received much attention, but primarily in a mathematical context. Following their approach, we have demonstrated that Wilson basis functions may indeed be constructed via an inverse Zak transformation, or via Daubechies' algorithm that was derived through the connection to theory of frames. Although both approaches have their appealing aspects, our prime objective is to promote and demonstrate the use of the basis functions for electromagnetic field problems.
We have demonstrated the expansion of complex-source beam (CSB) wavefields in a Wilson basis. CSBs are Maxwell-compliant wavefields that may be used as sources to illuminate an optical interface. We have demonstrated the illumination of an interface by a divergent CSB for three angles of incidence between 0°and 60°, and visualized the power distribution in coefficient space of the Wilson basis (cf. Fig. 9). In a similar fashion, we have expanded modal electromagnetic fields in a Wilson basis. In a companion paper (Floris and de Hon 2018), we consider interface problems consisting of optical fibers and homogeneous media, to which modal fields and CSBs in a Wilson basis are key ingredients. Moreover, the strong localization of Wilson basis functions in both the spatial and the spectral domain aid in the convergence of Green's function integrals in Floris and de Hon (2018).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Appendix: A summary of a numerical recipe for hðxÞ andĥðnÞ In this appendix, we give a short, streamlined and notation consistent summary of the construction of the numerical agorithm that was given in the mathematical paper of Daubechies et al. (1991). With the aid of the spectral Gaussian window functionĝ m that was defined in Eq. (11), one can make a frame in the spectral domain that is of a similar form as shown in Eq. (1), but with x ! n. Upon choosing a ¼ 0:5 and b ¼ 1, we obtain the spectral frame functionŝ g qr ðnÞ ¼ e jpqnĝ m ðn À rÞ; q; r; 2 Z; n 2 R: In Daubechies et al. (1991), the window in Eq. (11) is defined along the spatial coordinate, but later in that paper, the spatial and spectral arguments n and x are interchanged. We want to stick to consistent notation for spatial and spectral functions, and deliberately defined both the windowĝ m in Eq. (11) and the frameĝ qr in Eq. (33) in the spectral domain. By doing so, we deviate from the procedure in Daubechies et al. (1991) by effectively interchanging forward and inverse Fourier transform. However, because the functions h andĥ both are real-valued even functions, there are no essential differences. The construction of the functionĥ involves applying an operatorP on g, in such a way that it satisfieŝ where the operatorP is defined by P ¼ X q;rP qr ;P qrĝm ¼ hĝ qr ;ĝ m iĝ qr : In order to obtainĝ qr such that it constitutes a frame, we require that it satisfies the frame bounds in Eq.
(3), so that the operatorP is bound as well by AI P BI; 0\A B\1: Here, I is the identity operator. As a first estimate forP, we chooseP est ¼ AþB 2 I. With the aid of Eq. (36), we recognize that the norm is bounded, so that we know that the operatorP À1=2 exists and is bounded as well (Daubechies et al. 1991). The approach we adopt to determineP À1=2 is to rewrite it first as, In this form, we recognize that a Taylor series expansion ofP À1=2 converges, because of Eq. (37). In combination with Eq. (34), the Taylor series expansion reads, This expression can be cast into a recurrence scheme, making it suitable for numerical evaluation. The coefficient in the sum can be calculated efficiently using the hypergeometric recurrence relation Upon inserting the expansion of operatorP of Eq. (35) into Eq. (39), we may subsequently write the k repeated applications of I À 2P=ðA þ BÞ on the window functionĝ m as where b k qr satisfies the following recurrence relation, The initial condition b 0 q 0 r 0 is obtained from Eq. (33), by observing that ĝ m ðnÞ ¼ĝ 00 ðnÞ ¼ X q;rĝ qr ðnÞd q0 d r0 : The inner products hĝ qr ;ĝ q 0 r 0 i in Eq. (42) are linear in the second argument, and are known in closed form, hĝ qr ;ĝ q 0 r 0 i ¼e À mp 2 ðrÀr 0 Þ 2 e À p 8m ðqÀq 0 Þ 2 e ÀjpðqÀq 0 Þðrþr 0 Þ=2 : Finally, upon combining Eqs. (39), (40) and (41), and changing the order of summation, we infer that the functionĥ may numerically be constructed bŷ The advantage of changing the order or summation is that the coefficients a qr need to be computed only once. The inverse Fourier transformed window h has a similar construction, i.e., hðxÞ ¼ 2 ffiffiffiffiffiffiffiffiffiffiffiffi A þ B p X qr a qr g qr ðxÞ; with g qr ðxÞ ¼ ðÀ1Þ qr e Àj2prxĝ 1 m ðx À q=2Þ: