Boundary integral equation method for resonances in gradient index cavities designed by conformal transformation optics

In the case of two-dimensional gradient index cavities designed by the conformal transformation optics, we propose a boundary integral equation method for the calculation of resonant mode functions by employing a fictitious space which is reciprocally equivalent to the physical space. Using the Green’s function of the interior region of the uniform index cavity in the fictitious space, resonant mode functions and their far-field distributions in the physical space can be obtained. As a verification, resonant modes in limaçon-shaped transformation cavities were calculated and mode patterns and far-field intensity distributions were compared with those of the same modes obtained from the finite element method.

geometries from the circular shape and for the calculation of highly-excited resonant states (short wavelength regime) since it uses only the field information at the discretized cavity boundary, while FEM which is based on the 2D discretized space causes a heavy computational load 20,21 . In addition, because an outgoing-wave boundary condition at infinity is naturally built into the Green's function in BEM, one does not need to artificially construct an absorbing boundary environment at the outermost computational region such as the perfectly matched layer (PML), which is required in FEM for the truncation of the unlimited outside region 22 . Thus, BEM has been extensively used in calculating resonant modes in various 2D uniform index dielectric cavities, e.g., deformed cavities 23,24 , coupled cavities 25 , and annular cavities 26 .
In this paper, we show that conventional BEM can still be used to calculate resonant modes of TCs with a spatially-varying refractive index profile determined by an optical conformal mapping. To this end, we introduce a reciprocal virtual (RV) space which is related to the physical space by the inverse conformal mapping. The cavity in the RV space has a uniform refractive index and there one can build a boundary integral equation (BIE) for the interior region of the cavity.
The paper is organized as follows. The wave equations in a wholly transformed (WT) space are derived from Maxwell's equations in an original virtual (OV) space. And then we formulate BEM for the resonant modes of TCs by using the RV space. Finally, we demonstrate the validity of our BEM with examples of limaçon-shaped TCs.

GRin cavities Designed by conformal transformation optics
The theories of TO have been independently proposed by Ulf Leonhardt (conformal TO with Helmholtz equation in complex plane) 1 and Sir John Pendry (general TO with Maxwell's equations) 2 as design methodologies of their optical invisibility cloaks. The general TO is based on the form invariance of Maxwell's equations under a general coordinate transformation; as the consequences of the form invariance, electric and magnetic fields are renormalized and the constitutive parameters (i.e., the electric permittivity ε and the magnetic permeability μ) are changed to anisotropic tensor quantities 2,27 .
Here, we will derive the conformal TO from the general TO and then obtain the wave equation for infinitely-long cylindrical dielectric cavity in a WT space. We start from an OV space described by Cartesian coordinates (u, v, w), where the Maxwell's equations in frequency domain (fields ~ e −iωt ) without sources or currents in linear isotropic dielectric media are written as 13 In Eq. (1.1), E and H are time-independent complex electric and magnetic fields, respectively; the magnetic permeability μ is given by μ 0 μ r , where μ 0 and μ r are the vacuum permeability and the relative permeability, respectively; and the electric permittivity ε is given by ε ε r 0 , where ε 0 and ε r are the vacuum permittivity and the relative permittivity, respectively.
Under a general coordinate transformation from the OV space to a WT space described by coordinates (x, y, z), Maxwell's equations that keep their forms invariant in the WT space transforms as The relation of the constitutive parameters between the OV space and the WT space are given by where Λ is the Jacobian matrix of the transformation from the OV space to the WT space, and is given by (1 4) According to the so-called material interpretation 28 or manipulation-of-material viewpoint 29 , one can think that the effect of coordinate transformation is encoded in the redefined fields and anisotropic constitutive parameters in unchanged flat Cartesian space. In conformal TO, the coordinate transformation is given by an analytic function ζ = f(η) = x(u, v) + iy(u, v) of a complex variable η = u + iv with z = w, so the Jacobian matrix Λ is reduced to a block diagonal form, www.nature.com/scientificreports www.nature.com/scientificreports/ . The analytic function ζ = f(η) always produces a conformal coordinate transformation, which satisfies the Cauchy-Riemann equations, Due to these relations and the orientation-preserving property of conformal mappings (i.e., Λ > det 0), the constitutive parameters in the WT space turn out to have a simple diagonal form with the scale factor of a conformal mapping squared as the third diagonal element as follows 30 , . For this case of conformal TO, let us consider an infinitely-long cylindrical dielectric cavity with a circular cross section of radius R 0 = 1 in the OV space ( Fig. 1(a)) from which cylindrical cavity with a deformed cross section in a WT space ( Fig. 1(b)) can be obtained by a conformal mapping, taking z-axis and w-axis along the length of the cylinder 14,19 . Hereafter, we use the dimensionless coordinates (x, y) rescaled as (x/R 0 , y/R 0 ). For the transverse magnetic (TM) polarized mode of the cylindrical cavity, only the z-component, E z ′ among the components of electric field is non-zero, i.e.,  where ∆′ ≡ ∂ + ∂ x y 2 2 , and the spatially-varying refractive index n′(r) in the WT space is given by where n 0 is the refractive index of the cavity in an OV space which is equal to ε μ r r ; the free space wave number is k = ω/c, where c is the speed of light in vacuum given by ε μ 1/ 0 0 ; Ω in(ex) denotes the interior (exterior) region of the cavity. Equation (1.9) can be also obtained starting from the 2D Helmholtz equation under conformal mapping, which is a well-known result in waveguide theory 31 . At the dielectric interface, both are continuous across the cavity boundary 32 , The normal derivative is defined as ∂ ⊥ ≡ p(r) ⋅ ∇| r where p(r) is the outward normal unit vector at point r on the boundary curve Γ of the interior (exterior) region of the cavity. The TM modes of a cylindrical cavity in the WT space are supported by the refractive index n′ originating from the anisotropic electric permittivity, ε′, not from the anisotropic magnetic permeability, μ′ as can be seen from Eq. (1.8). It is noted that the electric field wave function ′ E x y ( , ) z in the WT space and its counterpart function E w (u, v) in the OV space are the same, as can be shown by Eqs. (1.2b) and (1.5).
For the transverse electric (TE) polarized mode of a cylindrical cavity in a WT space, among the components of the magnetic field, only the z-component, From these equations, the following wave equation for the TE mode can be obtained by using another Maxwell equation, ∇′ ⋅ (μ′H′) = 0, TE modes of a cylindrical cavity in the WT space are supported with the refractive index n′ originating from the anisotropic magnetic permeability, μ′, not from the anisotropic electric permittivity, ε′, as one can see from Eq. (1.12). It is also noted that the magnetic field wave functions

Boundary Element Method for Resonant Modes in Transformation Cavities
The GRIN cavity in WT space obtained from the homogeneous disk cavity in the OV space by an optical conformal mapping is surrounded by inhomogeneous material that fills infinitely extended outer region, which is an unphysical situation. Thus, we replace the outside inhomogeneous index with 1, i.e., the refractive index of air, as shown in Fig. 1(c) and then call the same cavity a transformation cavity (TC) in the physical space, which has gradually varying refractive index profile, n(x, y) in the interior region of the cavity according to Eq. (1.10) 14 . In other words, conformal mapping is only applied to the interior region (Ω ∼ in ) of the cavity in the OV space to obtain a TC in the physical space. Hereafter, primed quantities and operators in the WT space used in the previous section will be replaced with unprimed ones in the physical space for convenience. Now let us formulate a BEM for optical resonant modes in a TC with inhomogeneous refractive index in physical space. We have to solve the following scalar wave equation For asymptotically large r, the resonant mode function ψ(r) must satisfy the outgoing-wave boundary condition, where h k (θ) is the far-field angular distribution of the radiation emission. The outgoing-wave boundary condition leads to solutions exponentially decaying in time with discrete complex eigenvalues k with Im(k) < 0. So their angular frequency ω becomes complex. The lifetime τ of these so-called 'resonant modes' or 'resonances' is given by the imaginary part of the angular frequency as τ = −1/[2Im(ω)]; the quality factor Q of each resonance is with the oscillation period of light wave T = 2π/Re(ω). The complex wave function ψ equals E z and H z for TM and TE polarizations, respectively. The continuity relations at the cavity-air interface are given by where ψ in(ex) and ∂ ⊥ ψ in(ex) are wave functions and their normal derivatives from the interior (exterior) of the cavity, respectively, and n in (=n 0 |dζ/dη| −1 ) is the inhomogeneous refractive index from the interior region of the cavity evaluated at the boundary 32 .
In BEM, the above scalar wave equation is replaced by corresponding BIEs by using the Green's function and Green's second identity, and then the cavity boundary Γ is discretized. In our case, the Green's function for the interior region in physical space is not known since the refractive index is spatially varying there, so we cannot obtain the BIE from the interior region of the cavity, while for the exterior region of the uniform refractive index the BIE of the conventional BEM can be used. The Green's function in the exterior region is defined as the solution of where δ(r − r′) is the Dirac δ-function, and r, r′ are arbitrary points in the exterior region of the cavity. The Green's function G ex (r, r′; k) for the exterior region is given by the zeroth order Hankel function of the first kind, By multiplying Eq. (2.1) by the Green's function, subtracting the resulting equation from the product of Eq. (2.5) and ψ(r), and applying Green's second identity, we obtain the resulting boundary integral equation (BIE) as follows 17 where ∂ ⊥ is the normal derivative at point r on the boundary curve Γ of region Ω ex and s ≡ s(r), which is the arc length along Γ at r in the physical space; PV means Cauchy principal value integration. In order to build the BIE for the interior region of cavity, we introduce the RV space which is obtained from the physical space by the inverse conformal mapping, η = f −1 (ζ), as shown in Fig. 1(c, d). In general, the inverse conformal mapping is not a one-to-one mapping and therefore the RV space corresponding to the OV space should be selected. In the RV space, functions, vectors, differential operators, and other relevant symbols will be expressed with tildes. Under the inverse conformal mapping, the scalar wave Eq. (2.1) transforms to www.nature.com/scientificreports www.nature.com/scientificreports/ We note that the cavity in the RV space has a constant refractive index profile as shown in Fig. 1(d). Thus, for the interior region of the cavity in the RV space, we can build the BIE, in on the boundary of the cavity in the RV space is given by the zeroth order Hankel function of the first kind, 14) can be discretized as follows, In the above equation, we used the relation, ψ ψ = ∼ m m , which comes from the fact that electric field wave functions in the physical space and their counterpart functions in the RV space are the same as mentioned above and the normal derivative relation, is evaluated at the midpoint  s m ; details of derivation of these relations can be found in ref. 33 . Thus, for the case of TM resonant modes, the 2N × 2N matrix M(k) is defined by Finally, the wave functions in the interior and exterior regions of cavities can be obtained from following BIEs (see Eq. (2.7a)) as follows~~~~~~~~~~~~~ψ The wave function ψ in (r′) in the interior region of the cavity in the physical space can be obtained , where r′ is related to ′  r by the conformal mapping.

numerical calculation in Limaçon-Shaped transformation cavities
The conformal transformation, which maps the unit circle to the limaçon 14,34 , is given by where n 0 is the constant refractive index of a circular cavity in the OV space. The refractive index  n u v ( , ) derived from Eq. (2.9) in the RV space is given bỹ , n 0 = 2.0, and β = 1.0 in the physical space and the RV space, respectively. We note that the refractive index  n u v ( , ) given by the inverse conformal transformation of Eq. (3.1) is a double-valued function. The proper branch of refractive index  n u v ( , ) must be chosen so that the circular cavity in the RV space is same as the circular cavity in the OV space.
Using the BEM formulated in the previous section, we obtained TM resonances when  = . 0 15, n 0 = 1.8, and β = 0.769. The positive scaling factor β is introduced to obtain anisotropic WGMs supported by TIR which is called conformal Whispering Gallery Mode (cWGM) 14 . The condition for TIR in generic TCs is given by |dζ/dη| −1 ≥ 1; in this case, the condition becomes  β β ≤ = + 1/(1 2 ) max and we use β = β max which is the minimal condition for TIR. The dimensionless resonant wave numbers k res R 0 can be obtained from = M kR det[ ( )] 0 0 through the root-finding algorithm such as the Newton-Raphson method. All of resonances in the range of 9.7 < Re[kR 0 ] < 10.7 and 0 < Im[kR 0 ] < −0.2 are marked in the complex k-space as shown in Fig. 2. Since each mode has correspondence with the resonances in a uniform index circular cavity in the OV space which are labeled by the azimuthal mode number m and the radial mode number l, we simply denote the modes by (m, l). Additionally, our BEM can generate the unphysical modes with Im[kR 0 ] = 0, so-called 'spurious solutions' , as in the conventional BEM 17 . We can sort out them by comparing the Q-factor and the corresponding near-field patterns.
In order to confirm the validity of our BEM, we compared the results of BEM with that of the corresponding simulation performed by COMSOL Multiphysics v.5.3, a commercial FEM-based electromagnetic solver. In Fig. 3, we depict the normalized values of |∂ ⊥ ψ| 2 and |ψ| 2 along the arc length s normalized with total arc length L for a resonant mode, (14,1) with the complex wave number k res R 0 = 9.785 − i 0.00158 and one can see that the two results match almost exactly. The (14,1) cWGM has a high-Q factor (Q ≈ 3097) because of the evanescent leakage by TIR.
We plotted the near-field wave patterns and the far-field intensity distributions of the cWGM obtained by the two methods in Fig. 4. The far-field distribution exhibits bidirectionality due to the tunneling emission at the rightmost position of the limaçon shaped TC where the refractive index is lowest. In Fig. 4(a), the top, bottom-left, and bottom-right patterns are the intensity, the real part, and the imaginary part of the complex wave function obtained from BEM, respectively. Figure 4(b) is the results for the same mode obtained from FEM. As shown in Fig. 4(c), the far-field intensity distribution obtained with BEM also agree well with the result obtained with FEM. Incidentally, the BEM based on the Green's function can obtain the wave function value at any spatial point using only ∂ ⊥ ψ and ψ at the cavity boundary and it is one of the notable advantages of this method. www.nature.com/scientificreports www.nature.com/scientificreports/ In our case, all resonances except those with m = 0 are exist in nearly degenerate pairs due to the breaking of rotational symmetry. The mode pairs are formed with even-and odd-parity with respect to the mirror symmetry axis (x-axis) and are very close to each other in the complex k-space as shown in Fig. 2. In Fig. 5, we depicted   www.nature.com/scientificreports www.nature.com/scientificreports/ the mode intensity patterns and the corresponding far-field intensity distributions of a nearly degenerate (14,1) cWGM pair. The even-parity ( Fig. 5(a)) and the odd-parity mode (Fig. 5(b)) with respect to x-axis have k res R 0 = 9.785240667 − i 0.0015797513 and k res R 0 = 9.785240670 − i 0.0015797508, respectively. The far-field intensity distributions of the cWGM pair show out of phase interference patterns but all exhibit bidirectional emission features as shown in Fig. 5(c).
Using the BEM, we also obtained another resonant mode in a limaçon-shaped TC with = . 0 24  , n 0 = 2.0, and β = 1.0. The complex wave number of the mode is k res R 0 = 11.913 − i 0.107 and the intensity pattern of the resonant mode is shown in Fig. 6(a). This resonant mode has low-Q factor (Q ≈ 56) and the corresponding far-field intensity distribution is shown in Fig. 6(c). Contrary to the bidirectional far-field distribution of the above (14,1) cWGM, this low-Q mode has a unidirectional far-field intensity distribution. From these results, one can note that the Q-factors and emission directionalities of the resonant modes in TCs depend on the scaling factor β as well as the deformation parameter  35 . We also obtained the resonant mode under the same parameters by FEM and the resultant wave number of the mode is k res R 0 = 11.913 − i 0.108, which agrees well with the BEM result. Also, the intensity pattern of the corresponding mode and the far-field intensity distribution obtained by FEM are shown in Fig. 6(b,c), respectively, which nearly match the results of BEM as well.

conclusion
The cWGMs in TCs can simultaneously possess useful properties, such as an ultrahigh-Q factor and directional light emission, which are seemingly incompatible. In order to obtain the resonant mode functions in TCs designed by conformal TO, we have developed a pure BEM based on the fact that the Green's functions of interior and exterior regions are known in the RV and the physical spaces, respectively, which are connected by a conformal mapping. For the verification of our BEM, we have respectively, which are connected by a conformal mapping. For the verification of our BEM, we have calculated resonant modes in limaçon-shaped TCs and compared those with the corresponding results from FEM. Complex wave numbers, mode patterns, and far-field intensity distributions of the resonant modes obtained by the BEM almost match those obtained by FEM. Like the conventional BEM for homogeneous dielectric cavities, the newly proposed BEM has advantages in computing resonances in the TCs, e.g. a relatively simple formalism and efficiency, especially in finding highly-excited states. It can also be used for a measure of reliability for FEM calculations in non-piecewise-constant media. We expect that our method will be extended to the calculation of resonant modes in TCs with more complex geometries.