Exponentially Convergent Galerkin Method for Numerical Modeling of Lasing in Microcavities with Piercing Holes

: The paper investigates an algorithm for the numerical solution of a parametric eigenvalue problem for the Helmholtz equation on the plane specially tailored for the accurate mathematical modeling of lasing modes of microring lasers. The original problem is reduced to a nonlinear eigenvalue problem for a system of Muller boundary integral equations. For the numerical solution of the obtained problem, we use a trigonometric Galerkin method, prove its convergence, and derive error estimates in the eigenvalue and eigenfunction approximation. Previous numerical experiments have shown that the method converges exponentially. In the current paper, we prove that if the generalized eigenfunctions are analytic, then the approximate eigenvalues and eigenfunctions exponentially converge to the exact ones as the number of basis functions increases. To demonstrate the practical effectiveness of the algorithm, we ﬁnd geometrical characteristics of microring lasers that provide a signiﬁcant increase in the directivity of lasing emission, while maintaining low lasing thresholds.


Introduction
Various two-dimensional (2D) models of microdisk and microring lasers (see, e.g., [1,2]) can be investigated with the aid of a specific electromagnetic eigenvalue problem adapted to calculate the threshold values of gain, in addition to the emission frequencies, which is called the lasing eigenvalue problem (LEP) [3][4][5][6][7]. For 2D microcavity lasers with uniform gain, LEP was reduced in [8] to a nonlinear eigenvalue problem for the system of the Muller boundary integral equations (BIEs). This system, obtained by Muller in [9], is widely used in the analysis of electromagnetic-wave scattering from 2D and 3D homogeneous dielectric objects with smooth boundaries [10,11]. This is because Muller BIEs are the Fredholm second-kind equations, which guarantee the convergence of their numerical solutions. By the same reasons, the eigenmodes of fully active [6,8] and passive [12] microcavities can be calculated using Muller BIEs. Many authors, as in [12], have used a physical model called the complex-frequency eigenvalue problem (CFEP). It is based on the search for complexvalued natural frequencies of open passive resonators. To be able to build a general theory for both LEP and CFEP models, a generalized model was proposed in [8]. It obtained the following name: generalized complex-frequency eigenvalue problem (GCFEP) [8]. The reason for reducing GCFEP to the Muller BIEs was to get a system of weakly singular

GCFEP and Nonlinear Eigenvalue Problem for the Set of Muller BIEs
The formulation of GCFEP for 2D microcavity lasers with piercing holes is given in [18]. A generic geometry of the analyzed microcavities is shown in Figure 1. The air hole is domain Ω , the main body of the resonator is denoted as Ω , and the environment of the resonator is Ω . The boundaries Г and Г separate these regions. We suppose that the boundaries Г and Г are twice continuously differentiable, and and are the outer normal unit vectors to them, respectively. Material properties of the laser cavity can be characterized using either the dielectric permittivity or the refractive index. For non-magnetic materials, these two options are equivalent to each other. We used the latter choice because this is customary in optics and photonics research.
Thus, we assumed that the positive refractive index of the hole Ω and the environment around the resonator Ω are given. The complex-valued refractive index of the domain Ω is = − . We denote the given real part of by > 0 and the imaginary part, which is the real-valued parameter of GCFEP, by ∈ ℝ. The case of = 0 corresponds to the passive cavity (i.e., without material losses), < 0 is for the cavity with lossy material, and if the region Ω is filled in with a gain material, then > 0. In the latter case, the imaginary part of is called the gain index. We assumed that the electromagnetic field does not depend on the variable and depends on the time, as ∼ exp(− ). Herein, the speed of light in a vacuum was denoted by . We were looking for complex values of on the Riemann surface of the function ln . Because of the independence of the electromagnetic field on the variable, we are dealing with the scalar eigenfunctions of GCFEP ∈ ∖ {0}, each of which is the third element of the density vector E or H for the E-and H-polarization, respectively. We use the notation for the space of functions, which are complex-valued and continuous on Ω , Ω , and Ω and twice continuously differentiable on Ω , Ω , and Ω .
For each ∈ ℝ , the eigenvalues ∈ and the eigenfunctions ∈ ∖ {0} of GCFEP have to satisfy the Helmholtz equations, the transmission conditions, Material properties of the laser cavity can be characterized using either the dielectric permittivity or the refractive index. For non-magnetic materials, these two options are equivalent to each other. We used the latter choice because this is customary in optics and photonics research.
Thus, we assumed that the positive refractive index ν o of the hole Ω 1 and the environment around the resonator Ω 3 are given. The complex-valued refractive index of the domain Ω 2 is ν i = α i − iγ. We denote the given real part of ν i by α i > 0 and the imaginary part, which is the real-valued parameter of GCFEP, by γ ∈ R. The case of γ = 0 corresponds to the passive cavity (i.e., without material losses), γ < 0 is for the cavity with lossy material, and if the region Ω 2 is filled in with a gain material, then γ > 0. In the latter case, the imaginary part of ν i is called the gain index.
We assumed that the electromagnetic field does not depend on the variable x 3 and depends on the time, as ∼ exp(−ikct). Herein, the speed of light in a vacuum was denoted by c. We were looking for complex values of k on the Riemann surface L of the function ln k. Because of the independence of the electromagnetic field on the x 3 variable, we are dealing with the scalar eigenfunctions of GCFEP u ∈ U\{0}, each of which is the third element of the density vector E or H for the E-and H-polarization, respectively. We use the notation U for the space of functions, which are complex-valued and continuous on Ω 1 , Ω 2 , and Ω 3 and twice continuously differentiable on Ω 1 , Ω 2 , and Ω 3 .
For each γ ∈ R, the eigenvalues k ∈ L and the eigenfunctions u ∈ U\{0} of GCFEP have to satisfy the Helmholtz equations, the transmission conditions, and the outgoing Reichardt radiation condition [23], Here, the polar coordinates of point x are denoted by (ρ, ϕ), k o = kν o , k i = kν i . In Equations (4) and (5), we have the dependence of the coefficients on the polarization; namely, η o,i = ν −2 o,i and η o,i = 1 for the H-and E-polarization, respectively. The Hankel function of the first kind with the index l is denoted by H (4) and (5), which are related to the boundary conditions, have the following limit values (see, e.g., [24], p. 68): which are expected to exist uniformly on Γ 1,2 . The series in (6) converges uniformly and absolutely for any eigenfunction of GCFEP; besides, it is important to note that it is an infinitely term wise differentiable [8].
We denote the main sheet of L by L 0 and suppose that it is branch-cut along the negative imaginary semi-axis. At this point, we note that three types of GCFEP eigenfunctions exist, depending on the location of the eigenvalue k ∈ L 0 . Equation (6) is interchangeable to the common Sommerfeld radiation condition in the case of Im k = 0, The case of Im k > 0 corresponds to the situation when u exponentially decays as ρ → ∞ . The alternative case, Im k < 0 entails the eigenfunction u growing exponentially at infinity. An important note for our consideration is that the following property is true [8,18,23] for any k ∈ L, γ ∈ R, and u, which satisfies (3) and (6): Here, x ∈ Ω 3 , G o = (i/4)H (1) 0 (k o |x − y|). We denote Γ R as the circle with a big enough radius R, which center is located at x. This fact helps us explore all the eigenfunction types within the same framework.
We need to remember about the dependence of the imaginary part of k ∈ L 0 on γ ∈ R [8]. In the case of the passive cavity, where γ ≤ 0, without losses or with them, the GCFEP statement conforms with the usual statement of CFEP [12]. At this point, Im k < 0 for all the eigenvalues k ∈ L 0 . The alternative case is the active cavity, where γ > 0, and the imaginary part of k ∈ L 0 can be equal to or greater than zero. The pair (k, γ), where γ and k are positive, and the corresponding eigenfunction u satisfy all the conditions of LEP [6]. Particularly, condition (8) holds true.
Following [18], we use the integral representations of the eigenfunctions of the problem (1)-(6) in the domains Ω 1 , Ω 2 , and Ω 3 , respectively, where G i = (i/4)H (1) 0 (k i |x − y|). Equations (10) and (11) are well known (see, e.g., [24], p. 68). Equation (12) also holds true as we have Equation (9) for each value of parameters k ∈ L and γ ∈ R (see [8,18]). Now, we introduce the notations, and denote the space of continuous on Γ j , j = 1, 2, functions with the maximum norm by C j = C Γ j , j = 1, 2, C = C 1 × C 2 , and W = C × C. Furthermore, we denote the identical operator in the space W by I. Then, any solution of GCFEP (1)-(6) in terms (13)-(15) satisfies the following nonlinear eigenvalue problem for the set of Muller BIEs [18]: Here, we denote u j or v j , j = 1, 2 by the function g. The kernels have the following forms [18]: Some of the kernels K (q,s) j have logarithmic singularities and the others are continuous [13]. Consequently, the operator B(k, γ) : W → W is compact, and the operator A(k, γ) : W → W is Fredholm with index zero for every k ∈ L and γ ∈ R [13].
If u ∈ U is an eigenfunction of problem (1)-(6) corresponding to an eigenvalue k ∈ L for a value of the parameter γ ∈ R, then, defined by (13)- (15), functions u j and v j belong to the Banach spaces C j , j = 1, 2, respectively, and form a nontrivial solution w ∈ W of (16) with the same values of k and γ. This was proved in Theorem 3 of [21]. The assertion in the opposite direction relative to the statement of this theorem is not true, as, as in [18], we did not substitute representations (10)-(12) into (4) and (5), but added the limit values of them and their normal derivatives from both sides of the boundaries Γ 1 and Γ 2 term by term. However, the following result holds true (see Theorem 4 [21]). For each γ ∈ R and k ∈ I + problem (16) has only the trivial solution w = 0, w ∈ W. Here, I + denotes the strictly positive imaginary semi-axis of L 0 .

Galerkin Method
In the current section, we present a trigonometric Galerkin method for the numerical solution of problem (16). Assume that each contour Γ j has a parameterization ρ j (t) = Then, for any given γ ∈ R, we have Here, l, m = 1, 2, 3, 4, y = y(τ) ∈ Γ j , j = 1, 2, For the construction and investigation of the Galerkin method, it is convenient to consider the problem (16) in the Hilbert space H = (L 2 ) 4 , where L 2 denotes the space of square integrable functions with the inner product By T n ⊂ L 2 , we denote the subspace of all trigonometric polynomials of the order no greater than n with complex coefficients. Then, H n = (T n ) 4 ⊂ H is the subspace with the elements of the form, n ∈ T n . By p n : H → H n , we define the following projection operator: Here, Φ n : L 2 → T n is the Fourier operator, For q = −n, . . . , n, the vectors ϕ q (t) = exp(iqt) form the orthonormal basis in the space T n . We rewrite Equation (16) as follows We look for approximate solutions w n , w (4) n ∈ T n of the system of Equation (17) in the form w (m) Therefore, we have We calculate the unknowns α (m) q using the Galerkin method, where l = 1, 2, 3, 4. As the trigonometric functions are orthonormal, we can rewrite Equation (18) in the form of the following system of linear algebraic equations: where l = 1, 2, 3, 4, The system of linear algebraic Equation (19) is equivalent to the finite-dimensional linear operator equation Here, k ∈ L, A n : H n → H n , I is the unitary operator in the space H n . As usual, we denote by ρ(A n ) and by σ(A n ), the regular and the characteristic sets of the operator-valued function A n (k), respectively. Let also N , N , N , . . . be infinite sequences of the set of all natural numbers N. Theorem 1. For any given γ ∈ R, the following statements are true: 1.
If k 0 is an eigenvalue of A(k), then for each n ∈ N there exists an eigenvalue k n of A n (k) such that k n → k 0 (n ∈ N).

2.
If for each n ∈ N there exists an eigenvalue k n of A n (k), such that k n → k 0 ∈ L (n ∈ N), and w n is a normalized eigenfunction of A n (k n ), then (i) k 0 is an eigenvalue of A(k), (ii) {w n } n∈N is a discretely compact sequence and its cluster points are normalized eigenfunctions of A(k 0 ).

3.
For every compact L 0 ⊂ ρ(A), the sequence {A n (k)} n∈N is stable on L 0 , i.e., there exist n(L 0 ) and c(L 0 ), such that L 0 ⊂ ρ(A n ), A n (k) −1 ≤ c(L 0 ) for all k ∈ L 0 and n ≥ n(L 0 ). The proof of this theorem is based on the general results of the discrete convergence theory [25] applied for the investigation of approximate methods in the eigenvalue problem, where the parameter appears non-linearly [16]. Therefore, let us preface it with some definitions from [16].
As it is said, the sequence {w n } n∈N of vectors from the space H n discretely converges to the limit w ∈ H if w n − p n w → 0 as n → ∞, n ∈ N . Discrete convergence of the vectors will be denoted as w n → w(n ∈ N ) . The sequence of elements {w n } n∈N is called discretely compact if, for each subsequence {w n } n∈N , N ⊆ N , there exists a subset N ⊆ N and a vector w ∈ H, such that w n → w(n ∈ N ) .
Consider a bounded linear operator A ∈ L(H, H) and a sequence of finite-dimensional bounded linear operators {A n } n∈N . It is said that the sequence of operators {A n } n∈N approximate the operator A, if for any vector w ∈ H we have If the discrete convergence of vectors w n → w(n ∈ N) implies the discrete convergence of their images, A n w n → Aw(n ∈ N) , then the sequence of operators {A n } n∈N is said to converge discretely to A.
The sequence of operators {A n } n∈N is regular if from the boundedness of the sequence of vectors {w n } n∈N (thanks to the estimate w n ≤ const(n ∈ N)) and from the discrete compactness of the sequence of their operator images {A n w n } n∈N follows the discrete compactness of the sequence of the vectors {w n } n∈N . If a sequence of operators {A n } n∈N is regular and wherein approximates the operator A, then it is said that it regularly approximates A. The regular convergence of a sequence of operators is defined in similar way.
It is said that a sequence of operator-valued functions {A n (k)} n∈N regularly converges on L to an operator-valued function A(k), if for each converging numerical sequence k n → k 0 ∈ L(n → ∞) , the operator sequence {A n (k n )} n∈N regularly converges to the operator A(k 0 ).
Proof. Let us verify, that in the case under consideration, all conditions (b1)-(b5) of Theorem 2 of [16] are satisfied. Then, all the assertions of Theorem 1 hold true.
(b1) The operator-valued function A(k) : H → H is holomorphic and Fredholm on L, and its regular set is not empty. The holomorphicity and the Fredholm property of A(k) : W → W were proved in Theorem 2 in [18]. For A(k) : H → H , these properties are established similarly with the replacement of estimates for all norms in the space of continuous functions W on the corresponding estimates in the space of functions integrable with the square H.
In Theorem 4 from [21], it was established that, for any k ∈ I + , Equation (16) has only a trivial solution in the space W. The operator B(k) is weakly singular, therefore, any solution of Equation (16) from H must belong to the space W and can only be trivial for k ∈ I + . The operator-valued function A(k) : H → H is a Fredholm one, so for it we have I + ⊂ ρ(A).
(b2) For any n ∈ N, the operator-valued function A n (k) : H n → H n is holomorphic and Fredholm on L. Indeed, A(k) : H → H is holomorphic on L. As the operator p n is linear and bounded, A n (k) = p n A(k) : H n → H n has the same property. The Fredholm property of the operator-valued function A n (k) is obvious because of its finite-dimensionality.
(b3) On each compact set L 0 ⊂ L, the norms A n (k) are bounded uniformly in the parameters n ∈ N and k ∈ L 0 . Indeed, from the definition of the operator A n (k) and the equality p n = 1, It follows that A n (k) ≤ A(k) , n ∈ N, k ∈ L, However, because of The following estimate is correct: where c(k) is a continuous function on L: It is easy to see that to complete the verification of the required property, it suffices to compute the maximum of the function c(k) on the given compact set L 0 ⊂ L.
(b4) For each fixed value k ∈ L, the operator sequence {A n (k)} n∈N approximates the operator A(k). Indeed, by the definition of the operators A(k) : H → H and p n : H → H n , for any vector w ∈ H we have A n (k)p n w − p n A(k)w = p n A(k)p n w − p n A(k)w ≤ p n A(k) p n w − w → 0(n ∈ N).
The tendency to zero is a consequence of the tendency to zero of the norm of the remainder term of the segment of the Fourier series for any function from L 2 , estimate (22), and equality (21).
(b5) For each fixed value k ∈ L, the operator sequence {A n (k)} n∈N is regular. Indeed, the discrete compactness of the sequence of vectors {A n (k)w n } n∈N means that for any N ⊆ N, there exist N ⊆ N , such that the sequence {A n (k)w n = w n + B n (k)w n }, n ∈ N converges discretely to some z ∈ H. If the sequence {w n } n∈N is bounded, then there is a weakly convergent subsequence {w n } n∈N , N ⊂ N . As it is known, the compact operator B(k), takes it to a strongly converging one to some vector u ∈ H: B(k)w n − u → 0, n ∈ N Hence, by virtue of the inequality B n (k)w n − p n u ≤ p n B(k)w n − u and equality (21), it follows that the sequence {B n w n } n∈N converges discretely to u ∈ H. Thus, {w n } n∈N converges discretely to the vector w = z − u ∈ H, and the definition of the regularity of the sequence {A n (k)} n∈N is satisfied.
As usual, we denote various positive constants that do not depend on n by the same letter c. Let k 0 be an eigenvalue of A(k). We denote by G(A, k 0 ) the generalized eigenspace, i.e., the closed linear hull of all the generalized eigenfunctions of A(k) corresponding to k 0 . As the operator p n is linear, the next theorem follows from [17].

Theorem 2.
Assume that γ ∈ R is given, k 0 is an eigenvalue of A(k), and L 0 ⊂ L is a compact set with the boundary Γ 0 ⊂ ρ(A) so that L 0 ∩ σ(A) = {k 0 }. For each n ∈ N, we denote by ε n the maximum of the approximation error over k ∈ Γ 0 and w ∈ G(A, k 0 ), Then, ε n → 0 (n ∈ N) and the following estimations hold for almost all n ∈ N: (i) |k n − k 0 | ≤ cε 1/κ n for all k n ∈ σ(A n ) ∩ L 0 , where κ = κ(k 0 , A) is the multiplicity of the pole k 0 of the operator-valued function A −1 (k); (ii) k n − k 0 ≤ cε n , where k n is the weighted (proportionally to their algebraic multiplicities) mean of all the eigenvalues of A n (k) in L 0 , k n = ∑ k∈σ(A n )∩L 0 µ k · k, µ k = ν(k, A n )/ν(k, A), where ν(·, ·) is the algebraic multiplicity of the corresponding eigenvalue k; (iii) max{|k n − k 0 | : k n ∈ σ(A n ) ∩ L 0 } ≤ cε 1/l n n , where l n is the number of the different eigenvalues of A n in L 0 .
The next theorem follows from [26].
Theorem 3. Suppose that the conditions of Theorem 2 are fulfilled, ε n is defined in (23), G 0 (A, k 0 ) is the eigenspace of A(k) corresponding to the eigenvalue k 0 ∈ L, {k n } n∈N and {w n } n∈N are some sequences of eigenvalues k n of A n (k) and normalized eigenfunctions w n of A n (k), such that k n → k 0 (n ∈ N), and δ n is defined by the equality Then, for each eigenfunction w n there exists an eigenfunction w 0 = w 0 (w n ) ∈ G 0 (A, k 0 ), such that the following error estimate holds for almost alln ∈ N: Using the results from [27,28], pp. 270, 271, we derive the following approximation error estimates.
We solved the nonlinear eigenvalue problem (20) using the residual inverse iteration algorithm [29]. If the boundaries of the active cavity and the piercing hole were nonconcentric circles, then the entries of the Galerkin's matrix had the explicit expressions calculated carefully in [19,20]. We used them in the next section.

Numerical Results
Optimization of geometrical parameters of the microring resonator is aimed at finding such microlaser configurations for which a high value of the directivity D (see exact definition in [19,20]) and a low value of the threshold γ will be obtained. The modes that meet these requirements are of primary interest during the design of microlasers.
The studies were carried out for the H-polarization, because, for thin flat 3D laser cavities, which can be approximated with 2D models, the values of the thresholds for the H-polarized modes are lower than of the E-polarized modes [5]. This is because such a reduction dimensionality entails the replacement of the bulk refractive index with its effective value, which depends on the polarization. The results were obtained, as in [19,20], for the following parameter values: the refractive index in the domain Ω 3 and in the hole is equal ν e = 1, the real part of the refractive index in the active region is α i = 2.63, the dimensionless quantities are κ = ka 2 , d = |O 1 − O 2 |/a 2 , r = a 1 /a 2 . Here, the domains Ω 1 and Ω 2 are circles with centers O 1 and O 2 and radii a 1 and a 2 , respectively.
In our analysis, we used the same mode classification as in [19,20], with index e and o denoting the even and odd eigenfunction symmetry, respectively, with respect to the line of symmetry, which is the x 1 -axis. We note that, in the ideally circular cavity, the modes with small radial indices had their fields compressed to the cavity rim. This feature is reflected by their specific name as whispering-gallery modes.
Let us first investigate the dependence of the directivity D on the relative distance d between the center of the cavity and the center of the hole, and the relative radius of the hole r for the modes (11, 1, e/o). Figure 2 shows the dependence of D on d and r in the region, where the cavity contours do not cross each other. The value of D practically does not increase if the values of the parameters (r, d) satisfy the inequality d + 0.8r ≤ 0.4. Above the straight line described by the equation d = −0.8r + 0.4, an increase in the directivity coefficient D was observed for both modes. The regions of values (r, d) in which D is maximal are clearly distinguishable. For the odd mode, this is a vicinity of the point 0.32 and 0.49, for the even mode, this is a vicinity of the point 0.02 and 0.63. We see that the value of D for the even modes is higher than for the odd modes, in addition, obtaining quasi-unidirectional emission is impossible for odd modes [19,20]. Therefore, we carried out further studies for the mode (11, 1, e). the hole is equal = 1, the real part of the refractive index in the active region is = 2.63, the dimensionless quantities are = , = | − |/ , = / . Here, the domains Ω and Ω are circles with centers and and radii and , respectively.
In our analysis, we used the same mode classification as in [19,20], with index e and o denoting the even and odd eigenfunction symmetry, respectively, with respect to the line of symmetry, which is the -axis. We note that, in the ideally circular cavity, the modes with small radial indices had their fields compressed to the cavity rim. This feature is reflected by their specific name as whispering-gallery modes.
Let us first investigate the dependence of the directivity on the relative distance between the center of the cavity and the center of the hole, and the relative radius of the hole for the modes (11, 1, / ). Figure 2  is maximal are clearly distinguishable. For the odd mode, this is a vicinity of the point 0.32 and 0.49, for the even mode, this is a vicinity of the point 0.02 and 0.63. We see that the value of for the even modes is higher than for the odd modes, in addition, obtaining quasi-unidirectional emission is impossible for odd modes [19,20]. Therefore, we carried out further studies for the mode (11, 1, ). In addition to obtaining a high value of for directivity , it was necessary to maintain low values of the threshold . It is more convenient to search for low values of by maximizing the values of the function T = − log . Figures 3 and 4 show the dependences of the normalized wavenumber and the threshold gain index (in short, the threshold) for the mode (11, 1, e) on the relative distance between the centers of the cavity and the hole and the relative radius of the hole . We see that the normalized wavenumber ranges from 5.85 to 6.  In addition to obtaining a high value of for directivity D, it was necessary to maintain low values of the threshold γ. It is more convenient to search for low values of γ by maximizing the values of the function T = − log 10 γ. Figures 3 and 4 show the dependences of the normalized wavenumber κ and the threshold gain index γ (in short, the threshold) for the mode (11, 1, e) on the relative distance d between the centers of the cavity and the hole and the relative radius of the hole r. We see that the normalized wavenumber κ ranges from 5.85 to 6.   We searched for such pairs of values of and , for which took a high value, while at the same time the value of remained low enough. For this, we considered the target functions of the following forms: = + , where = − log . It was assumed that among the points of the local maxima of the target functions (25)-(27), we would find such pairs of values ( , ) for which a high value of the directivity and a low value of the threshold would be obtained. For the found points of the local maxima, a check should be performed with a control value. The value of the target function must be no less than the value obtained if the problem for the microdisk resonator (without an air hole) is considered. Solving the problem for the microdisk resonator with the same radius , and the same refractive index , we have = 5.2074, = 2. This means that, for the target function (25), the control value is = 10.4147, for the target function (26), the control value is = 7.2074, and for the target function (27), the control value is = 25.2074. Figure 5 shows the points of local maxima of the considered target functions. For function (25), there are nine points of the local maxima; for function (26), there are eight points of the local maxima; and for target function (27), there are four points of the local maxima. Observing the results, we see that some points are presented on all three panels of Figure 5. Next, by intersecting the sets of the local maxima of the considered target functions, for further research, we chose the following pairs of values ( , ): (0.02 and   We searched for such pairs of values of and , for which took a high value, while at the same time the value of remained low enough. For this, we considered the target functions of the following forms: = + , where = − log . It was assumed that among the points of the local maxima of the target functions (25)- (27), we would find such pairs of values ( , ) for which a high value of the directivity and a low value of the threshold would be obtained. For the found points of the local maxima, a check should be performed with a control value. The value of the target function must be no less than the value obtained if the problem for the microdisk resonator (without an air hole) is considered. Solving the problem for the microdisk resonator with the same radius , and the same refractive index , we have = 5.2074, = 2. This means that, for the target function (25), the control value is = 10.4147, for the target function (26), the control value is = 7.2074, and for the target function (27), the control value is = 25.2074. Figure 5 shows the points of local maxima of the considered target functions. For function (25), there are nine points of the local maxima; for function (26), there are eight points of the local maxima; and for target function (27), there are four points of the local maxima. Observing the results, we see that some points are presented on all three panels of Figure 5. Next, by intersecting the sets of the local maxima of the considered target functions, for further research, we chose the following pairs of values ( , ): (0.02 and We searched for such pairs of values of d and r, for which D took a high value, while at the same time the value of γ remained low enough. For this, we considered the target functions of the following forms: where T = − log 10 γ. It was assumed that among the points of the local maxima of the target functions (25)- (27), we would find such pairs of values (r, d) for which a high value of the directivity D and a low value of the threshold γ would be obtained. For the found points of the local maxima, a check should be performed with a control value. The value of the target function must be no less than the value obtained if the problem for the microdisk resonator (without an air hole) is considered. Solving the problem for the microdisk resonator with the same radius a 2 , and the same refractive index α i , we have T = 5.2074, D = 2. This means that, for the target function (25), the control value is F = 10.4147, for the target function (26), the control value is F = 7.2074, and for the target function (27), the control value is F = 25.2074. Figure 5 shows the points of local maxima of the considered target functions. For function (25), there are nine points of the local maxima; for function (26), there are eight points of the local maxima; and for target function (27), there are four points of the local maxima. Observing the results, we see that some points are presented on all three panels of Figure 5 (25) in Figure 5a, the points under consideration are numbered as 1, 2, and 4, respectively. For functions (26) and (27) in panels (b) and (c) of Figure 5, the points under consideration are 1, 2, and 3, respectively.
(a) (b) (c)  Tables 1-3 show the values of the numerical characteristics of the cavity, and the lasing modes corresponding to the points of the local maximum. Here, is the angle showing the emission direction, i.e., the target of the main beam in the far-field patterns (see Figure 6). The points in the tables are numbered in descending order of the value of the corresponding target function.   Tables 1-3 show the values of the numerical characteristics of the cavity, and the lasing modes corresponding to the points of the local maximum. Here, β is the angle showing the emission direction, i.e., the target of the main beam in the far-field patterns (see Figure 6). The points in the tables are numbered in descending order of the value of the corresponding target function.   Thus, in the course of the numerical experiments, we found that a quasiunidirectional emission could occur both at a small hole radius and at a relatively large hole radius. In the case of a small hole radius, the main beam was directed oppositely to the direction of the hole shift, while in the case of a large hole, the main beam was in the same direction as the hole shift. The maximum directivity was obtained with a small  Figure 6 shows the near-and far-field patterns for the points selected by intersecting the sets of the local maxima of the considered target functions (25)- (27), among which the pair 0.02 and 0.63 is of the greatest practical interest, because the directivity D is maximal, and the value of the threshold is the smallest among all three points. As it is known [20], in the case where the center of the cavity and the center of the hole coincide, the directivity factor is D = 2.00 and T = 5.2074. This means that if choosing any of the considered pairs of values (r,d), the directivity D becomes at least 2.5 times higher, while the threshold γ does not increase significantly.
Thus, in the course of the numerical experiments, we found that a quasi-unidirectional emission could occur both at a small hole radius and at a relatively large hole radius. In the case of a small hole radius, the main beam was directed oppositely to the direction of the hole shift, while in the case of a large hole, the main beam was in the same direction as the hole shift. The maximum directivity was obtained with a small relative radius of the piercing hole. These phenomena were studied by physical experiments in [22].

Conclusions
We presented the main steps in the reduction of GCFEP for a 2D laser with a piercing hole to a set of four coupled boundary integral equations of the Muller type. We explained the discretization of these equations with the Galerkin method and proved its convergence. We obtained the error estimates for the approximate eigenvalues and eigenfunctions dependent on the smoothness of the generalized eigenfunctions.
Finally, we calculated the on-threshold characteristics of the lasing modes of a circular microcavity with a shifted hole. In the numerical experiments, we varied the position of the piercing hole and the radius of the hole, and computed the changes in the lasing frequencies, directionalities, and thresholds. Our numerical investigation showed that a hole with a suitable radius located at a certain place could lead to notable growth of the directivity of the perturbed whispering-gallery mode emission, together with the preservation of its low threshold. Hence, a piercing hole radius and position in the 2D eccentric microcavity laser can be used as an engineering tool to efficiently control the directivity of emission.