{'fraction_non_alphanumeric': 0.04313857425890621, 'fraction_numerical': 0.019214234981869836, 'mean_word_length': 4.241379310344827, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An adaptive boundary element method for the transmission problem with hyperbolic metamaterials

In this work we present an adaptive boundary element method for computing the electromagnetic response of wave interactions in hyperbolic metamaterials. One unique feature of hyperbolic metamaterial is the strongly directional wave in its propagating cone, which induces sharp transition for the solution of the integral equation across the cone boundary when wave starts to decay or grow exponentially. In order to avoid a global refined mesh over the whole boundary, we employ a two-level a posteriori error estimator and an adaptive mesh refinement procedure to resolve the singularity locally for the solution of the integral equation. Such an adaptive procedure allows for the reduction of the degree of freedom significantly for the integral equation solver while achieving desired accuracy for the solution. In addition, to resolve the fast transition of the fundamental solution and its derivatives accurately across the propagation cone boundary, adaptive numerical quadrature rules are applied to evaluate the integrals for the stiff matrices. Finally, in order to formulate the integral equations over the boundary, we also derive the limits of layer potentials and their derivatives in the hyperbolic media when the target points approach the boundary.


Background
Hyperbolic metamaterials are a class of anisotropic electromagnetic materials for which one of the principal components of their relative permittivity or permeability tensors attain opposite sign of the other two principal components: In the above, the subscripts and ⊥ denotes the component parallel and perpendicular to the anisotropy axis respectively, and there holds Reε ⊥ · Reε < 0 or Reµ ⊥ · Reµ < 0.
Hyperbolic metamaterials can be realized, for instance, by metal-dielectric structures or by embedding arrays of metallic wires in a dielectric matrix by restricting free-electron motion to certain directions [19,34,37]. More recently, hexagonal boron nitride (hBN), α-phase molybdenum trioxide (α-MoO 3 ), α-phase vanadium pentoxide (α-V 2 O 5 ) and a few others emerge as natural hyperbolic materials that attain opposite signs for the in-plane and outof-plane components of the dielectric tensor [4,14,15,32,38]. Assume that the media is nonmagentic such that the permeability reduces to the unit tensor, then the dispersion relation for the time-harmonic (with e −iωt dependence) Maxwell's equations ∇ × E = iωµ 0 µH, ∇ × H = −iωε 0 εE is given by where ε 0 and µ 0 are the free-space permittivity or permeability, k 0 is the free-space wavenumber, and k 1 , k 2 and k 3 are the x 1 , x 2 and x 3 components of the wave vector respectively in the Cartesian coordinate. The first term in (1.1) corresponds a spherical isofrequency surface for the transverse electric (TE) polarized waves, while the second term gives rise to a hyperboloidal isofrequency surface for the transverse magnetic (TM) polarized waves when Reε ⊥ ·Reε < 0. It is seen that the TM waves remain propagating with arbitrarily large wave vectors in the hyperbolic medium, as opposed to evanescent in an isotropic medium. This unique property leads to many interesting applications of hyperbolic metamaterials ranging from sub-wavelength light manipulation and imaging to spontaneous and thermal emission modification [12,24,26,35,36].

Problem formulation
In this paper, we investigate the computation of the electromagnetic response from the wave interactions in the hyperbolic metamaterials. We focus on the two-dimensional problem when the medium is invariant along the x 3 direction and the wave is TM-polarized with the magnetic field H = (0, 0, u) T . The TE-polarized case is less interesting as it leads to an isotropic problem with the dispersion relation given by the first term of (1.1), and various existing computational methods can be applied to solve the problem. The Maxwell's equations in the hyperbolic medium for the TM-polarized polarization reduce to the scalar wave equation ∇ · (A∇u) + k 2 0 u = 0, (1.2) where the coefficient matrix .

Computational challenges
One unique feature of a hyperbolic medium is the strongly directional wave propagation inside a cone with the half cone angle given by arctan − Re ε 1 Re ε 2 . This induces sharp transition of the solution across the cone boundary when wave starts decaying/growing exponentially. The domain discretization methods such as finite element or finite difference are computationally very expensive to resolve such singular behaviors of the solution. Here we propose a boundary element method for solving for the transmission problem (1.4). Integral equation solvers have played an increasing role in the computational electromagnetics in the past several decades due to its powerfulness in solving large-scale problems by discretization over the boundaries of the objects only; see, for instance, the monographs [1,9,10,11,25,30] and the references therein. The application of the integral equation method for the hyperbolic media requires us to address several new computational challenges as described below. First, the fundamental solution in the hyperbolic medium is strongly directional in the propagating cone and it decays exponentially across the cone boundary (see Figure 1). As a result, when computing the stiff matrices in the boundary element method, one needs to apply adaptive numerical quadrature rules to evaluate the integrals with the fundamental solution or its derivatives as kernels in order to achieve sufficient accuracy for discretization. Second, the solution of the integral equation formulation along the interface Γ attains sharp transitions when wave front reaches the boundary, especially for hyperbolic media with small loss (see examples in Section 4). Here we employ a two-level a posteriori error estimator and an adaptive mesh refinement procedure to resolve the singularity of the integral equation solution in an accurate and efficient manner. The theory and computation for the adaptive boundary element methods (BEM) are mature in solving elliptic boundary value problems [18]. By using a posteriori error estimator, the adaptive procedure chooses a sequence of meshes such that the numerical error decays in an optimal manner with increasing dimension of the approximation spaces. There exist a variety of error estimators for elliptic boundary value problems, including residual type estimators, space enrichment type estimators, averaging estimators, etc [5,6,7,8,16,17,20,21,39]. The two-level a posteriori error estimator proposed here for the hyperbolic transmission problem belongs to the family of the space enrichment type estimators. The principal idea is to use improved approximation of solutions u h/2 and ∂νu h/2 obtained over a uniform refined mesh with mesh size h/2 to replace the exact solution u and ∂νu in the numerical error u − u h and ∂νu − ∂νu h . The Dörfler strategy is then applied to mark and refine the mesh where local errors u h/2 − u h + ∂νu h/2 − ∂νu h are large. Such an adaptive procedure allows for the reduction of the degree of freedom significantly while achieving desired accuracy for the solution, as demonstrated by the numerical examples in Section 4. The goal of our work in this paper is to demonstrate the efficacy and accuracy of the adaptive algorithm for the two dimensional problems. Its application in three dimensions will be investigated in the forthcoming work. The rest of the paper is organized as follows. In Section 2 we introduce layer potentials and derive their limits as the target points approach the boundary. The limiting formulas recover the formulas in the isotropic medium when ε 1 = ε 2 . The boundary integral equation for the transmission problem (1.4) is then formulated in Section 2. The adaptive Galerkin boundary element method is desribed in Section 3, where we introduce the adaptive numerical quadrature and the two-level a posteriori error estimator. Several numerical examples are given in Section 4 to illustrate the accuracy and efficiency of the adaptive algorithm. The paper is concluded with brief remark about the proposed computational approach and the future work along this direction.
2 Layer potentials and boundary integral equations for the transmission problem

Layer potentials and integral operators
Here and henceforth, for a hyperbolic material with ε ⊥ (x) ≡ ε 1 and ε (x) ≡ ε 2 , we let be the A-deformed distance between x and y, where A is given in (1.3) and the function √ z is understood as an analytic function defined in the domain C\{−it : t ≥ 0} such that √ z = |z| 1 2 e 1 2 i arg z . Note that when Imε j > 0 (j = 1, 2),r(x, y) lies in the first quadrant of the complex plane. Letν = Aν be the A-deformed normal vector over the interface Γ. Correspondingly, the derivative of a given function ϕ along the directionν is defined as 0 (kr(x, y)) be the fundamental solution, which satisfies (1.2) when x = y and is outgoing when |x−y| → ∞. Here H It is well-known that the single layer potential v(x) is continuous throughout R 2 . In what follows, we derive the limits of the double layer potential and the derivatives of two layer potentials as x approaches Γ. The limiting formulas recover the classical limiting formulas when ε 1 = ε 2 = 1 [28].

Lemma 2.1
The double-layer potential w(x) with the continuous density ϕ can be continuously extended from Ω toΩ and R 2 \Ω to R 2 \Ω respectively with the limit w(x ± hν(x)).
√ ε 1 ε 2 ln(r(x, y)) be the fundamental solution of (1.2) when k 0 = 0, and w 0 be the corresponding double layer potential: Note that the difference of two double layer potentials w(x) and w 0 (x) is continuous in R 2 , thus it suffices to verify (2.2) for w 0 (x). The proof can be further reduced to the special case when the density function ϕ ≡ 1, this is because for an arbitrary density function ϕ, one can write the double layer potential as and the latter is continuous throughout R 2 when ϕ is continuous. Next we verify the assertion by assuming that ϕ(x) ≡ 1 and showing that be the small disk with radius δ centered at x and ∂B δ (x) be its boundary. It follows from the Green's formula that where ν denotes the unit normal exterior to B δ (x). A direct calculation yields By evaluating the above integral explicitly and noting that lim

Proposition 2.2
The derivative of the single-layer potential v(x) with the continuous density ϕ can be continuously extended from Ω toΩ and R 2 \Ω to R 2 \Ω respectively with the limit ∇v(x ± hν(x)) ·ν(x).
Proof. This can be observed from the formula wherein we have used the relation ∇ x Φ(x, y) = −∇ y Φ(x, y). Note that the second term is continuous in R 2 , thus an application of Lemma 2.1 leads to (2.4).
Lemma 2.4 Let w(x) be the double layer potential defined in (2.1) and ϕ ∈ C 1 (Γ), then there holds ∇w(x ± hν(x)) ·ν(x), and d ds denotes the derivative with respect to the arc length.
By using the equation (1.2) and the relation ∇ Hence, applying the integration by parts leads to (2.8) Now for fixed z ∈ Γ, we set x = z + hν(z). By virtue of (2.8) and Lemma 2.3, it follows that A straightforward calculation leads to Similarly, by setting x = z − hν(z), one can obtain the same formula for ∂w − (z) ∂ν .

Boundary integral equation formulation for the transmission problem
. Letν j = A j ν be the A j -deformed normal vector over the interface Γ. We define the integral operators S j , K j , K j and N j for x ∈ Γ as follows: Note that Φ j (x, y) in the lossy hyperbolic medium attains the same singularity as the fundamental solution with ε 1 = ε 2 = 1, thus from the standard theory of the boundary integral operators (cf. [28]), we have the following lemma for the above integral operators.
be the difference of two integral operators with associated kernels. The volume integral operators P j and Q j are defined as Applying the Green's formula in Ω 1 and using the formula (2.3), we obtain for x ∈ Ω 1 that Similarly, for x ∈ Ω 2 there holds (2.14) Recall that f 1 (x) and f 2 (x) are localized sources with compact support, the volume integrals above only need to be evaluated over their support regions. By taking the limit of (2.13) and (2.14) when x approaches the interface Γ and applying Lemma 2.1, we achieve the integral equations on Γ: (2.13) and (2.14) respectively, and taking the limit when h → 0 yields 1 2  18) respectively, and applying the continuity condition (1.5) for the wave field across the interface, we obtain the following system of integral equations: (1) g (2) .

Galerkin boundary element method
be the basis of the space U h and V h respectively. By , the Galerkin approximation (3.1) leads to the following linear system: (3.2) In the above, c (1) and c (2) represent the unknown vectors that take the following form: 2 , · · · , c 1 , c 2 , · · · , c The (m, n)-th entry of the matrices N h , K h , S h , and I h are given by

4)
I h and K h are the transposes of I h and K h respectively. In light of (2.12), we evaluate N h (m, n) using the formula The m-th element for the vectors g h and g (2) h are given by h (m) = g (2) , φ (2) m .

Adaptive numerical integration
The entries of the local stiff matrices (3.3) -(3.5) boil down to the evaluation of integral in the form of e mn =

Em En
Θ(x, y)φ n (y) φ m (x) ds y ds x , in which Θ(x, y) represents the kernel Φ j (x, y) or ∂Φ j (x, y) ∂ν j (y) , and φ m (x) represents the basis functions in U h or V h . As pointed out in Section 1.3, the fundamental solution in the hyperbolic medium is strongly directional in the propagating cone C = {x ; x T (ReA) −1 x > 0}, and it attains sharp transition near the cone boundary ∂C. To capture the variation of the kernels accurately, we employ the adaptive numerical quadrature to compute e mn . In more details, for a given small real number τ > 0, let us introduce the domain that includes the cone boundary. τ is chosen so that Ω τ contains the region where Θ(x, y) attains very large derivatives. Let E mn = {x − y ; x ∈ E m , y ∈ E n } be the set of relative locations between the source and target points when evaluating e mn . If E mn ∩ Ω τ = ∅, then Θ(x, y) changes smoothly over the domain (x, y) ∈ E m × E n , thus there is no need for adaptivity and one can still apply fast algorithms, such as the fast multipole methods, to evaluate e mn for fixed m and all n satisfying E mn ∩ Ω τ = ∅ in an efficient manner [9,23]. Otherwise, if E mn ∩ Ω τ = ∅, we compute e mn adaptively to resolve the kernels accurately as described below. Note that when m = n, the kernel Θ(x, y) is weakly singular and the singular part of e mn can be evaluated analytically. Hence we only need to consider the case when Θ(x, y) is nonsingular. By a change of variable, e mn is expressed as in the parameter space. The integral in (3.7) is computed via the adaptive Lobatto quadrature rule [22]. Let R , and then calculating the integral using the Lobatto quadrature rule over these four subregion regions. The recursive procedure stops when the relative difference of the two approximations at level and + 1 is smaller than the prescribed tolerance. We refer the readers to [2] for more details of the recursive procedure. Since the region Ω τ is usually thin with small τ (τ = 0.1 is chosen in variety of numerical experiments demonstrated in Section 4), the cardinal number of the set {(m, n)|E mn ∩ Ω τ = ∅} M 2 . In addition, the recursive adaptive Lobatto quadrature rule for the integral inside Ω τ convergences fast, thus the adaptive integration only accounts for a small percentage of the overall cost in assembling the stiff matrices.

The two-level a posteriori error estimator and mesh refinement
Let Γ h = {E 1 , E 2 , · · · , E M } be a mesh over Γ with the mesh size h, andΓ h = {Ê 1 ,Ê 2 , · · · ,Ê 2M } be a uniform refinement of Γ h with the mesh size h/2. The corresponding Galerkin solution in the finite-dimensional space respectively. We use a two-level a posteriori error estimator where the exact solution in the numerical error is replaced by the solutionφ h obtained over the uniformly refined meshΓ h . To this end, we define the first h − h/2 based estimators as follows: and h ]. The estimators η (1) and η (2) defined above require the computation of the solutions φ h and φ h at both discretization levels, which could be computationally burdensome. Furthermore, φ h is expected to be more accurate than φ h , thus the latter becomes a temporary result that is not useful once η (1) and η (2) are calculated. In order to reduce the computational cost and avoid such redundancy, afterφ h is computed, we useφ h to approximate the solution φ h by projecting the refined solutionφ h over the finite-dimensional space H h [18]. In addition, we localize the estimators by using h 1/2 -weighted H 1 and L 2 seminorm for φ (1) h and φ (2) h , respectively. The localization yields error indicator over each element that is computable and can be used to design the mesh refinement strategy. More precisely, we define the second error estimatorsη (1) andη (2) as follows: where the error indicators over each element E m are given by , (3.9) and Π (3.10) Recall thatφ (1) h andφ (2) h are the Galerkin approximations of the solution for the transmission problem and its normal derivative, we see thatη provides a h 1/2 -weighted H 1 -seminorm estimator for the solution u over the interface Γ. We point out thatη (j) (j = 1, 2) alone has been used as a posteriori error estimator for solving elliptic boundary value problems with Dirichlet or Neumann boundary conditions, and it has been shown that the estimator is efficient and reliable in the sense that the true error for the numerical solution is bounded below and above by the estimatorη (j) [18]. Here we combine the two together in (3.10) as the error estimator for the transmission problem in the hyperbolic media. A variety of numerical examples in Section 4 demonstrate the effectiveness of the proposed error estimator. The theoretical investigation on its robustness remains to be investigated in the future.
The Dörfler strategy is employed to mark and refine the mesh Γ h using the error estimator η. For given 0 < γ < 1, we find the minimal set E ⊂ Γ h such that (3.11) and each marked element in E is then divided into two sub-elements with equal size. The complete adaptive strategy starts from an initial discretization of the interface Γ. The calculation of the estimator (3.10) and the refinement procedure (3.11) is repeated until η < σ for certain prescribed tolerance σ. This is summarized in the following algorithm:

Algorithm: The adaptive boundary element method
Given tolerance σ and the initial mesh size h, generate the initial Γ   (1) and φ (2) , respectively. The corresponding numerical solutions returned from the adaptive algorithm are denoted byφ (1) h andφ (2) h , which are obtained over the refined meshΓ ( ) h . We also define the relative errors by lettinĝ where φ (1) and φ (2) are reference solutions obtained with high-order accuracy.
Example 1 Assume that the interface Γ is an ellipse with semi-axes a = 2 and b = 1 respectively. The interior domain is a hyperbolic medium while the exterior domain is vacuum. The permittivity values in Ω 1 and Ω 2 are given by ε (1) 2 = −2 + 0.02i and ε (2) 1 = ε (2) 2 = 1, respectively. The wavenumber k 0 = 1, and the source is located in Ω 1 such that f 1 = −δ(x) and f 2 = 0. Due to the smoothness of the interface, we compute the reference solutions φ (1) and φ (2) with high-oder accuracy by the Nyström scheme, which is a spectral method using the trigonometry interpolant over Γ [28]. By parameterizing the interface Γ with x 1 = a cos θ and x 2 = b sin θ, the initial mesh Γ (0) h is generated with grid points (a cos θ m , b sin θ m ) Table 1 shows the mesh sizes and the corresponding numerical errors for different levels of refinements when the adaptive procedure is applied, in which M ( ) denotes the number of grid points being used in the mesh Γ ( ) h . As expected, the numerical errorê (1) for φ (1) is relatively small over the initial mesh Γ (0) h , whileê (2) is large due to the singular behavior of φ (2) . It is observed that the two-level a posteriori error estimatorη is effective in identifying the solution singularities and a local mesh refinement near the singularities reduces the numerical error significantly after each refinement. The adaptive procedure terminates with = 4, and there holds M (4) = 253 for the mesh Γ (4) h ; see Figure 2 (left) for a plot of the mesh Γ (4) h . The numerical solutionsφ (1) h andφ (2) h in the final stage of the adaptive procedure are plotted in Figure 3, which are solved over the meshΓ (4) h . The corresponding numerical errors areê (1) = 0.0044 andê (2) = 0.0745, respectively. As a comparison, if the mesh is quasi-uniform, one needs a much more refined meshΓ h with a total number of 2M grid points to achieve a comparable accuracy, in whichM = 700. This is illustrated in the last column of Table 1. In Figure 4 we also plot the wave field in the domain Ω 1 and Ω 2 , which are computed via the formulas (2.13) and (2.14) using the solutionsφ (1) h andφ (2) h obtained from the adaptive algorithm. The strong directional propagation of the wave inside the hyperbolic medium and the multiple reflections by the interface Γ is clearly seen.  (1) and e (2) for the adaptive BEM and non-adaptive BEM.  h andφ (2) h obtained from the adaptive algorithm. θ represents the polar angle over the ellipse.  The adaptive procedure also terminates with = 4 and the mesh Γ (4) h attains a total of 249 grids points as shown in Figure 2 (right). The final numerical solutionsφ (1) h and φ (2) h are plotted in Figure 5, and the corresponding numerical errors areê (1) = 0.0034 and e (2) = 0.0616, respectively. For completeness, we collect all the mesh sizes and the corresponding numerical errors for different levels of refinements in Table 2. The computation with a quasi-uniform mesh to achieve a comparable accuracy is also shown. We see that the number of degrees of freedom is reduced by about 4 times when the adaptive procedure is applied. Figure 6 demonstrates the wave field in the domain Ω 1 and Ω 2 . The wave is strongly directional while penetrating through the interior hyperbolic medium and being reflected at the interface.  (1) and e (2) for the adaptive BEM and non-adaptive BEM.  h andφ (2) h obtained from the adaptive algorithm. θ represents the polar angle over the ellipse.
1i. The exterior domain Ω 2 is occupied by vacuum. When the source is located at (0.3, 0.1) and the wavenumber k 0 = 1, the mesh sizes and the corresponding numerical errors for the adaptive algorithm are shown in Table  3. Here we obtain the reference solutions φ (1) and φ (2) with high-oder accuracy by using a uniform mesh with large degree-of-freedom. The final mesh Γ (3) h attains a total of 238 grids points as shown in Figure 7 (top). The numerical errors for the solutions obtained fromΓ (3) h areê (1) = 0.0011 andê (2) = 0.0722, respectively. In contrast, the same level of accuracy is obtained by a quasi-uniform meshΓ h with a total number of 2M grid points, in which M = 672. Figure 8 plots the wave field in the domain Ω 1 and Ω 2 . It is observed that multiple reflections by the interface induce strong wave interactions inside Ω 1 .  (1) and e (2) for the adaptive BEM and non-adaptive BEM.   (2) 2 = 1 + 0.05i for the interior and exterior domain, respectively. The wavenumber k 0 = 2π, and the source is located in exterior domain with f 1 = 0 and f 2 = −δ(x − x 0 ), where the source location is x 0 = (0.5, 0.3). The adaptive boundary element still successfully leads to a reduction of numerical errors to the desired accuracy when the mesh refinement is performed locally, as demonstrated by Table 4. Due to the singularity of the solution, one needs a quasi-uniform mesh with the mesh size h max = 0.0013 to achieve an accuracy that is obtained by an adaptive mesh with h max = 0.02. The number of degrees of freedom for the former is more than 5 times higher than the latter. The final mesh in the adaptive procedure and the corresponding wave field in the domain are shown in Figures 7 (right) and 9, respectively.   Example 5 In the final example, the interior domain Ω 1 attains a wedge shape and it is placed in the vacuum. Such geometry has been be used for super-focusing of the electromagnetic waves near the sharp tip [31]. The permittivity values for the hyperbolic medium are ε (1) 1 = 2, ε (1) 2 = −3 + 0.03i. The source is located in Ω 1 so that f 1 = −δ(x − x 0 ) with x 0 = (0.1, 0.1) and f 2 = 0, and the wavenumber k 0 = 1. The adaptive algorithm terminates after three mesh refinements, which yield the numerical solution with an accuracy ofê (1) = 0.0011 andê (2) = 0.0722 (see Table 5). As demonstrated in Figure 11, the wave propagates toward the tip while being reflected by the boundary of the wedge.

Discussions
An adaptive boundary element method was presented in this paper for solving the transmission problem with hyperbolic metamaterials. Compared to the discretization with a quasi-uniform mesh, the adaptive approach is able to resolve the singular behavior of the solution with local mesh refinement, which reduces the number of the degrees of freedom and the overall computational cost significantly. There are several theoretical and computational issues to be explored along this direction. Although numerical examples show the efficacy of the error estimator and the accuracy of the adaptive procedure, the robust analysis of the estimator and the convergence analysis of the algorithm have not been carried out yet. As pointed out previously, this work is mainly on the demonstration of the adaptive algorithm for the two dimensional problems, its application in three dimensions has not been explored. It is expected that the adaptive algorithm would offer even larger reduction on the computational cost for 3D simulations. However, the computation becomes more challenging for mesh refinement and numerical integration in 3D. Note that for the full Maxwell's equations, the Dyadic Green's functions used in the integral formulation is highly anisotropic with coexistence of the cone-like pattern due to emission of the extraordinary TM-polarized waves and elliptical pattern due to emission of ordinary TE-polarized waves [33]. Finally, efficient integral equation methods for hyperbolic metamaterials with unbounded domains (e.g. layered media) and other settings of practical interest remain to investigated. The propagation nature of waves with arbitrarily large wave vectors inside the propagating cone requires new treatement in developing the computational algorithm.