On the robustness of bound states in the continuum in waveguides with lateral leakage channels

Bound states in the continuum (BICs) are trapped or guided modes with frequencies in radiation continua. They are associated with high-quality-factor resonances that give rise to strong local field enhancement and rapid variations in scattering spectra, and have found many valuable applications. A guided mode of an optical waveguide can also be a BIC, if there is a lateral structure supporting compatible waves propagating in the lateral direction, i.e., there is a channel for lateral leakage. A BIC is typically destroyed (becomes a resonant or a leaky mode) if the structure is slightly perturbed, but some BICs are robust with respect to a large family of perturbations. In this paper, we show (analytically and numerically) that a typical BIC in optical waveguides with a left-right mirror symmetry and a single lateral leakage channel is robust with respect to any structural perturbation that preserves the left-right mirror symmetry. Our study improves the theoretical understanding on BICs and can be useful when applications of BICs in optical waveguides are explored.


I. INTRODUCTION
A bound state in the continuum (BIC) is a trapped or guided mode that decays to zero in a spatial direction along which radiation modes with the same frequency and wavevector (when appropriate) can propagate to or from infinity [1][2][3]. In photonics, BICs have been found on various structures including waveguides with a local distortion [4], periodic structures sandwiched between two homogeneous media [5][6][7][8], waveguides with lateral leakage channels [9][10][11][12][13][14][15][16][17][18][19], rotationally symmetric periodic structures surrounded by a homogeneous medium [20,21], and anisotropic multilayer structures [22]. Applications of the BICs are mostly related to resonances with arbitrarily large quality factors (Q factors) that appear when the structure or the wavevector are properly perturbed [23][24][25][26][27]. High-Q resonances give rise to strong local field enhancement and abrupt variations in transmission/reflection spectra that are useful for lasing, sensing, switching, and nonlinear optics applications.
An important theoretical question is about the robustness (i.e., the continual existence) of the BICs with respect to small structural perturbations. If a BIC is protected by a symmetry, i.e., there is a symmetry mismatch between the BIC and the radiation modes, it is typically robust with respect to symmetry-preserving perturbations. For periodic structures with the in-plane inversion symmetry and the up-down mirror symmetry, some BICs unprotected by symmetry are also robust with respect to perturbations preserving these two symmetries [28][29][30][31][32]. An optical waveguide may have a lateral leakage * Corresponding author: ljyuan@ctbu.edu.cn channel, if the structure away from the waveguide core has guided modes that can propagate in the lateral direction. A BIC on such a structure may or may not be symmetry-protected. In a recent work [18], Bykov et al. studied a waveguide with a rectangular core in a background with a slab, showed that even the BICs unprotected by symmetry are robust with respect to variations in the structural parameters.
In this paper, we analyze the robustness of BICs in optical waveguides with respect to arbitrary structural perturbations. Based on an all-order perturbation method first developed in [30], we show analytically that a generic BIC in a waveguide with a lateral mirror symmetry and a single radiation channel is robust with respect to any perturbation that preserves the lateral mirror symmetry. The rest of this paper is organized as follows. In Sect. II, we present the necessary background concerning the waveguide structure and the BICs. In Sect. III, we formulate a related scattering problem and present some useful properties of the scattering solutions. Section IV contains the main result about robustness, including the precise conditions and the proof. In Sect. V, we present numerical examples to illustrate the continuous existence of BICs when structural parameters are varied. The paper is concluded with a brief discussion in Sect. VI.

II. GUIDED MODES AND BICS
We consider a three-dimensional (3D) z-invariant lossless dielectric waveguide with a core embedded in a layered background. The dielectric function ε of the waveguide depends only on the two transverse variables x and y, that is, ε = ε(r) for r = (x, y). We assume the structure is symmetric in x, the core is contained in the region given by |x| < W/2 for some W > 0, and the background is a layered medium with a dielectric function ε s that depends only on y. Thus In addition, we assume the background itself is a 2D waveguide, typically consisting of a slab, a substrate and a cladding. We denote the dielectric constants of the cladding and the substrate by ε 0 and ε 2 , respectively, and assume ε 2 ≥ ε 0 ≥ 1. As an example, we show a ridge waveguide in Fig. 1 trapezoid. For such a 3D z-invariant waveguide, a guided mode decays to zero as r = x 2 + y 2 → ∞, and depends on z as exp(iγz), where γ is a real constant (the propagating constant). The electric field of the guided mode can be written as Re[E e i(γz−ωt) ], where ω is the angular frequency, and E is a vector function that depends only on r. From the frequency-domain Maxwell's equations, it is easy to show that E satisfies where k = ω/c is the freespace wavenumber, c is the speed of light in vacuum, and e 3 = (0, 0, 1) is the unit vector in the z direction. For simplicity, we call E the electric field and normalize the guided mode such that The layered background given by the dielectric function ε s (y) is a 2D waveguide with transverse-electronic (TE) and transverse-magnetic (TM) modes. To distinguish these modes from the eigenmodes of the original 3D waveguide, we call them slab modes. A TE slab mode is characterized by a zero y component of the electric field and a scalar function u(y) satisfying the eigenvalue problem: (11) with the condition u(y) → 0 as y → ±∞, where η te is the real positive effective index of the slab mode (kη te is the propagation constant). Since the slab mode must decay exponentially in the substrate and the cladding, we have η te > √ ε 2 . In addition, we assume the TE slab mode is normalized such that 1 L +∞ −∞ |u(y)| 2 dy = 1, for the same L in Eq. (5). In the region given by x > W/2, if this TE slab mode also depends on z as exp(iγz) and is outgoing or tends to zero as x → +∞, then its electric field is given by where Notice that if γ < kη te , then α te > 0 and the field given in Eq. (12) is outgoing as x → +∞. If γ > kη te , then α te is pure imaginary with a positive imaginary part, and the field decays to zero as x → +∞. For x < −W/2, the TE slab mode with the same z dependence has an electric field given by where and it is either outgoing or tends to zero as x → −∞.
The case for the TM slab mode is similar. It has a zero y component in the magnetic field and a scalar function v(y) satisfying the eigenvalue problem d dy with the condition v(y) → 0 as y → ±∞, where η tm > √ ε 2 > 0 is the effective index of the mode. Similar to the TE modes, we normalize v(y) such that |v(y)| 2 dy = 1.
Assuming the same dependence on z, the electric field of the TM slab mode can be written down as where and v ′ (y) = dv/dy. The TM slab mode with its electric field given in Eqs. (18) and (19) is either outgoing (if γ < kη tm ) or evanescent (if γ > kη tm ) as x → ±∞. At a fixed frequency, the layered background has a finite number of slab modes. We order the TE and TM slab modes separately according to their effective indices, and denote them as u j , η te j for j = 0, 1, . . ., N − 1, and v j , η tm j for j = 0, 1, . . ., M −1, where N and M are the numbers of TE and TM modes, respectively. It is known that the TE and TM modes are interlaced, namely, their effective indices satisfy In the slab region given by |x| > W/2, a general timeharmonic field is a superposition of the finite number of slab modes and a continuum of radiation modes [33]. If the field depends on z as e iγz and is outgoing as x → ±∞, then the electric field can be written as where the "+" and "−" subscripts are chosen for x > W/2 and x < −W/2, respectively, and In the above, α te j , α tm j , f te +,j , f te −,j , f tm +,j , f tm −,j are defined for the jth TE or TM mode following Eqs. (13), (20), (14), (16), (21) and (22); a ± j and b ± j are the coefficients of the slab modes; and the terms with the subscript "⋄" are radiation modes that depend continuously on a real wavenumber β for the y direction. The wavenumber for the x direction is where ε 2 is the dielectric constant of the substrate. For each β, f te ±,⋄ and f tm ±,⋄ are the TE and TM radiation modes [33], and a ± ⋄ and b ± ⋄ are the corresponding coefficients. If E is the electric field of a guided mode (of the 3D waveguide) with a frequency ω and a propagation constant γ, then it must satisfy Eq. (23) for properly cho- A regular guided mode satisfies the condition kη te 0 < γ, thus all α te j , α tm j and α ⋄ (β) are pure imaginary with positive imaginary parts, and E → 0 as x → ±∞. A BIC is a special guided mode with a positive propagation constant satisfying This means that α te 0 and probably other α tm j or α te j are positive, thus the layered background has at least one slab mode that can radiate power to x = ±∞ and propagate in the z direction as the BIC. Since the BIC must decay to zero as x → ±∞, it is clear that a ± 0 and probably the coefficients of other slab modes (that radiate out power in the x direction) must be zero. We are particularly interested in BICs satisfying In that case, α te 0 is real, all other α te j , α tm j , and α ⋄ (β) are pure imaginary. Therefore, only the fundamental TE slab mode can radiate power to x = ±∞, and the coefficient a ± 0 of the BIC must be zero. Due to the reflection symmetry of the structure in the x direction, we have either a + 0 = a − 0 or a + 0 = −a − 0 , if the BIC satisfies (9) or (10), respectively.

III. SCATTERING SOLUTIONS
For ω and γ satisfying condition (30), the fundamental TE slab mode can propagate in the slab re- This implies that we can consider scattering problems using the fundamental TE slab mode as incident waves. The solutions of the scattering problems are needed in Sect. IV for proving the robustness of BICs. If (with the assumed z dependence e iγz ) is the incident wave in the left slab region (x < −W/2), then the total field contains a reflected wave Ru − (r), for in the left slab region, and a transmitted wave T u + (r) in the right slab region (x > W/2), where R and T are the reflection and transmission coefficients. Due to the reflection symmetry in x, if we specify u − (r) as the incident wave in the right slab region, the total field contains a reflected wave and a transmitted wave in the right and left slab regions, respectively, and the reflection and transmission coefficients are exactly the same. In addition, if condition (30) is satisfied, the fundamental TE slab mode provides the only radiation channel, thus, the 2 × 2 scattering matrix R T T R is unitary. This implies that To construct a scattering solution satisfying the symmetry condition (9) or (10), we need to specify incident waves in both left and right slab regions. For condition (9), we let the incident waves be Cu + for x < −W/2 and Cu − for x > W/2, where C is a nonzero constant, then any solution of this scattering problem satisfies Although the far field given in the right hand side above already satisfies condition (9), the total field E(r) may not, because the solution is not unique if a BIC exists at the same ω and γ. However, we can define a vector field E(r) as in Eq. (8), thenÊ and (E +Ê)/2 solve the same scattering problem. We can replace our original solution by (E +Ê)/2, then the new solution, still denoted as E(r), satisfies condition (9). By choosing a proper constant C, we can ensure that condition (7) is also satisfied. Since |R + T | = 1, we have R + T = e iϕ for a real phase ϕ. If we let C = e −iϕ/2 , then Eq. (32) becomes Due to the possible non-uniqueness, the above solution may not satisfy condition (7). However, we can definẽ E as in Eq. (6), thenẼ and (E +Ẽ)/2 solve the same scattering problem. If we replace the original solution by (E +Ẽ)/2, then the new E(r) satisfies condition (7).
Since we are concerned with the case where a BIC exists at the same ω and γ, we denote the electric fields of the BIC and the scattering solution by E 0 and E (s) , respectively. Assuming both E 0 and E (s) satisfy conditions (7) and (9), we can replace E (s) by E (s) + C 1 E 0 for a real constant C 1 (if necessary), such that the following orthogonality condition is satisfied: In summary, we have constructed a scattering solution E (s) that satisfies (7), (9) and (34). Similarly, we can construct another scattering solution that satisfies conditions (7), (10) and (34).

IV. ROBUSTNESS OF BICS
In this section, we study the robustness of BICs in optical waveguides based on an all-order perturbation method developed in our previous works [30,32]. The robustness refers to the continual existence of a BIC under small structural perturbations. The original unperturbed waveguide is given by a dielectric function ε 0 (r) satisfying conditions (1) and (2). It is assumed that the unperturbed waveguide has a non-degenerate BIC with an electric field E 0 (r), a frequency ω 0 and a propagation constant γ 0 . Moreover, ω 0 and γ 0 must satisfy condition (30), where k should be replaced by k 0 = ω 0 /c, and η te 0 and η tm 0 are effective indices of the fundamental TE and TM slab modes at frequency ω 0 . Without loss of generality, we assume E 0 (r) satisfies conditions (7) and (9). The dielectric function of the perturbed waveguide is given by where δ is a small real number, s(r) is a real O(1) function satisfying condition (1) and s(r) = 0 if |x| > W/2 and if |y| is sufficiently large. In the following, we show that under a generic condition, the perturbed waveguide has a BIC with a frequency ω near ω 0 , a propagation constant γ near γ 0 , and an electric field E(r) near E 0 (r). Following the procedure developed in [30,32], we construct the BIC in the perturbed waveguide by expanding E(r), γ and k = ω/c in power series of δ: Note that for any real γ near γ 0 , the governing equations (3) and (4) have a solution satisfying an outgoing radiation condition as x → ±∞, but the frequency ω is complex in general. Similarly, for any real ω near ω 0 , there is an outgoing solution with a complex γ. These outgoing solutions with a complex ω or a complex γ are the resonant and leaky modes, respectively. Our objective is to determine a BIC that decays to zero as r → ∞, and it only exists for a particular pair (ω, γ) near (ω 0 , γ 0 ). Therefore, both ω and γ must be determined together with the field. Inserting expansions (36)-(38) into Eqs. (3) and (4), and comparing the coefficients of δ j for j ≥ 1, we obtain the following equation for E j (r): where and F j , for j > 1, are given in Appendix A. Both L and B are differential operators independent of j. L is the differential operator associated with Eq. (3) for ω 0 , γ 0 and the unperturbed waveguide. Since E 0 satisfies condition (7), we can verify that BE 0 also satisfies (7). The vector function F j depends on γ m , k m and E m for 0 ≤ m < j. For each j ≥ 1, we need to determine a real k j and a real γ j , show that E j can be solved from Eq. (39), and it decays to zero exponentially as r → ∞ and satisfies conditions (7) and (9). We establish the result recursively. For the jth step, it is assumed that for each m satisfying 0 ≤ m < j, we already have a real k m , a real γ m , and a vector function E m satisfying (7) and (9). If such an E j exists, we can show (see Appendix A) that where E (s) is a scattering solution constructed in Sec. III (for ω 0 , γ 0 , and the unperturbed waveguide) satisfying conditions (7), (9) and (34). Replacing L E j by the right hand side of Eq. (39), we have The above can be written as where Since BE 0 satisfies condition (7), a 11 and a 21 are real. In addition, a 12 > 0, and a 22 = 0 due to the orthogonality condition (34). Therefore, all entries of A are real. Meanwhile, since F j depends on γ m , k m and E m for 0 ≤ m < j, γ m and k m are real, and E m satisfy condition (7), F j also satisfies condition (7), and thus b 1j and b 2j are real. Therefore, when A is invertible, γ j and k j can be solved, and they are real. Although Eqs. (41) and (42) are derived assuming Eq. (39) has a solution, the linear system (43) does not depend on E j . In the following, we show that if γ j and k j satisfy (43) and A is invertible, then Eq. (39) indeed has a solution E j that decays to zero exponentially as r → ∞ and satisfies conditions (9) and (7).
Since F j depends on E m for 0 ≤ m < j and E m → 0 exponentially as r → ∞, the right hand side of Eq. (39) tends to zero as r → ∞ and can be regarded as a source distributed around the waveguide core. Therefore, we require the solutions of Eq. (39), if they exist, to satisfy an outgoing radiation condition. Since only the fundamental TE slab mode can radiate out power (in the x direction), the outgoing radiation condition gives rise to for some coefficients d ± j . Since L E 0 = 0, Eq. (39) is a singular equation, and it has solutions if and only if the right hand side is orthogonal with the kernel of L . Since the BIC is non-degenerate, the kernel of L is the onedimensional vector space spanned by E 0 . Therefore, if condition (41) is satisfied, Eq. (39) has solutions satisfying (50). The solutions are not unique, but all solutions of Eq. (39) must have the same asymptotic coefficients d ± j . This is so, because the difference of two solutions of Eq. (39) satisfies a homogeneous equation and does not radiate out power to infinity. Furthermore, since all E m , for 0 ≤ m < j, satisfy condition (9) and the perturbation profile s(r) satisfies condition (1), it is easy to verify that E j also satisfies condition (9). Therefore, d + j = d − j . The second condition, Eq. (42), can be used to show that d ± j = 0. For any positive h, let Ω h be the domain given by |x| < h and |y| < ∞, then Eqs. (39) and (42) clearly imply that In Appendix B, we show that where C = 0 is the constant given in Eq. (33). This implies that d ± j = 0, and thus E j → 0 exponentially as r → ∞.
It can be easily verified that the right hand side of Eq. (39) satisfies condition (7). Thus, for any solution E j , we can construct a vector fieldẼ j as in Eq. (6), theñ E j and (E j +Ẽ j )/2 also satisfy Eq. (39) and they decay to zero exponentially as r → ∞. Therefore, we can replace E j by (E j +Ẽ j )/2 (if necessary), and assume E j satisfies condition (7).
The above proof is for a BIC satisfying the symmetry condition (9). If the BIC satisfies condition (10), we need to use the scattering solution that also satisfies (10), then d + j = −d − j , and d + j − d − j should appear in the right hand side of Eq. (52). If s(r) does not satisfy condition (1), the above proof fails, because we no longer have d + j = ±d − j , then d ± j = 0 in general, and E j does not decay to zero as r → ∞. To ensure that matrix A is invertible, the BIC of the unperturbed waveguide is required to satisfy where E (s) satisfies the orthogonality condition (34). In summary, if the original waveguide is symmetric in x and has a non-degenerate BIC satisfying conditions (30) and (53), and the perturbation is also symmetric in x, then the perturbed waveguide has a BIC with a slightly different frequency and a slightly different propagation constant.

V. NUMERICAL EXAMPLES
In this section, we present some numerical examples to demonstrate the robustness of BICs in optical waveguides. We start with a silicon ridge waveguide with a rectangular ridge, a SiO 2 substrate, and an air cladding.
The dielectric constants of the cladding, the ridge and the slab, and the substrate are ε 0 = 1, ε 1 = 12.25, and ε 2 = 2.1025, respectively. The geometric parameters, as shown in Fig. 1, are W = 0.68µm, h r = 0.07µm, h s = 0.15µm, and θ = 0. This waveguide has a BIC with freespace wavenumber k 0 ≈ 0.7266(2π/µm) (i.e. freespace wavelength λ 0 ≈ 1.3762µm) and propagation constant γ 0 ≈ 1.4794(2π/µm). In Fig. 2(a), we show the dispersion curves of the fundamental TE and TM slab modes (solid blue curves) and the dispersion curves of some resonant modes (solid red lines, for real part of k = ω/c only). The BIC, marked as " * " in Fig. 2(a), satisfies condition (30) and is a special point on a band of resonant modes. In Fig. 2(b), we show the Q factor of this particular band for γ near γ 0 . The BIC has been scaled to satisfy conditions (7) and (5) with L = 1µm. The z-components of the electric field and a scaled magnetic field of the BIC are shown in Fig. 2(c) and (d), respectively. The scaled magnetic field is obtained by multiplying the free space impedance to the original magnetic field. Due to the normalization, Eq. (5), all field components are dimensionless. From Fig. 2(c), it is clear that the BIC satisfies symmetry condition (10). Following Sect. III, we calculate a scattering solution satisfying conditions (7), (10)  It is easy to verify that the BIC satisfies condition (53). Therefore, according to the theory developed in Sect. IV, the BIC should be robust with respect to any perturbation that preserves condition (1), i.e., the left-right mirror symmetry. To validate the theory, we consider a perturbation profile for |x| < W/2 and h s < y < h s + h r , and s(r) = 0, otherwise, and study the waveguide numerically for δ ∈ [0, 1].
The numerical results confirm that the BIC exists continuously as δ is increased from 0 to 1. The propagation constant γ and freespace wavenumber k of the BIC are shown as functions of δ (solid blue curves) in Fig. 3(c) and (d), respectively. For δ = 0.2, the BIC is obtained with γ ≈ 1.4947(2π/µm) and k ≈ 0.7288(2π/µm). In Fig. 3(a) and (b), we show the Q factor and the dispersion curve, respectively, for a band of resonant modes containing the BIC [shown as the " * " in Fig. 3(b)]. In Fig. 3(c) and (d), we also show first order approximations, γ ≈ γ 0 + γ 1 δ and k ≈ k 0 + k 1 δ, as the red dashed lines, where γ 1 and k 1 are first order terms used in Sect. III. Based on the numerical solutions of the BIC and the scattering solution of the unperturbed waveguide, we obtain γ 1 ≈ 0.5240(1/µm) and k 1 ≈ 0.0781(1/µm). As a second test for the theory, we consider a waveguide with an isosceles trapezoidal ridge. For the same set of parameters, ε 0 , ε 1 , ε 2 , W , h r , and h s , we study the waveguide allowing θ, shown in Fig. 1, to increase from 0 • to 10.2 • . The numerical results reveal a BIC that depends on θ continuously. The propagation constant γ and freespace wavenumber k of the BIC are shown as functions of θ in Fig. 4(c)  the BIC, as usual, is a special point on a band of resonant modes. The Q factor and the dispersion curve of this band are shown in Fig. 4(a) and (b), respectively. The BIC is obtained with (γ, k) ≈ (1.6082, 0.7544)(2π/µm) and shown as the " * " in Fig. 4(b). Additional numerical results indicate that the BIC ceases to exist for θ > 10.2 • . It appears that the BIC in the waveguide with θ ≈ 10.2 • fails to satisfy condition (53), i.e., a 21 = 0 and matrix A is not invertible. In Fig. 4(e), we show a 21 as a function of θ. It indicates clearly that a 12 can become zero when θ is further increased.

VI. CONCLUSION
A 3D z-invariant optical waveguide may have a lateral leakage channel if the waveguide core is embedded in a background which itself is a 2D waveguide (typically a slab waveguide). A guided mode of the 3D waveguide is a BIC if its propagation constant γ is smaller than the propagation constant of the fundamental TE slab mode. In this paper, we studied BICs satisfying condition (30), so that only the fundamental TE slab mode can have the same z dependence as the BIC and still propagate in the lateral direction. For waveguides with the left-right mirror symmetry, i.e., condition (1), we showed that any BIC satisfying conditions (30) and (53) is robust with respect to structural perturbations that preserve condition (1). The left-right mirror symmetry is imposed so that the solutions satisfy either condition (9) or condition (10), and the radiation channels to the left and right (via the fundamental TE slab mode) are essential the same. Condition (53) is given assuming the BIC and the scattering solution are properly scaled, have the same lateral symmetry, and are orthogonal to each other, i.e., they are required to satisfy conditions (7), (34), and (9) or (10).
It is straightforward to extend our result to waveguides without the left-right mirror symmetry, but with leakage channels in only one lateral direction, for example, the positive x direction. A guided mode of such a waveguide is a BIC, if its propagation constant γ is smaller than the propagation constant of the right fundamental TE slab mode (of the layered structure for x > W/2). If γ is larger than the propagation constants of all other right slab modes and all left slab modes (of the layered structure for x < −W/2, if they exist), then the BIC exists in a continuum with only one radiation channel. Such a BIC, if it is non-degenerate and satisfies the generic condition (53), is robust with respect to any structure perturbation. The proof is slightly simpler than the one given in Sect. III, because we no longer need to worry about the symmetry conditions (9) and (10). It is only necessary to ensure the BIC and the scattering solution satisfy conditions (7) and (34).
The BICs studied in this paper are quite different from those in periodic structures sandwiched between two homogeneous media. They have a different number of direction along which the field is confined, and they exist in radiation continua provided by different waves. The left-right mirror symmetry assumed in this paper and the up-down mirror symmetry for biperiodic structures with propagating BICs [32], serve the same purpose, that is, to lock together the radiation channels in the opposite directions. For biperiodic structures, besides the up-down mirror symmetry, the in-plane inversion symmetry is also needed to ensure the robustness of propagating BICs [32]. In contrast, for waveguides with lateral leakage channels, robustness of BICs does not require any additional symmetry. This difference is related to our assumption that the waveguide is z-invariant and the perturbation is zindependent. Chongqing Statistics Postgraduate Supervisor Team Since L E 0 = 0, we have we can simplify the above as where Integrating both sides of Eq. (55) on a disk of radius a (denoted as D a ) and using Gauss's Law, we obtain where Γ a is the circle of radius a. Since E 0 → 0 exponentially as r → ∞, the line integral in the right hand side above tends to zero as a → ∞. Therefore R 2 E 0 · L E j dr = 0. Since we assume E j → 0 exponentially as r → ∞, E (s) · L E j is also integrable on R 2 and the integral is also zero.

APPENDIX B
Unlike the case considered in Appendix A, E j is outgoing as x → ±∞, and thus bounded at infinity, Since L E (s) = 0, we have Since E j and E (s) decay to zero exponentially as y → ±∞, we can integrate E (s) · L E j on Ω h , apply Gauss' Law, and obtain where e 1 = (1, 0, 0) is the unit vector along the x axis. Based in the asymptotic formulae (33) and (50), it is easy to verify that That is Eq. (52).