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

Abstract

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 find geometrical characteristics of microring lasers that provide a significant increase in the directivity of lasing emission, while maintaining low lasing thresholds.

1. 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 complex-valued 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 integral equations [13] on the boundary of the microcavity laser. However, there is no full equivalence between GCFEP and the eigenvalue problem for the system of Muller BIEs [14]. Namely, it was proven in [15] that for each eigenfunction of GCFEP there is a corresponding eigenvector of the system of Muller BIEs. Still, the assertion in the opposite direction is not true: there is one more problem that is reduced to the Muller BIEs, called “turned inside out GCFEP” [15]. If GCFEP and the turned inside out GCFEP together have only the trivial solutions, then the system of Muller BIEs has only the trivial solution [15], and the resolvent set of the corresponding operator-valued function is not empty. This result is important for the theoretical investigation of the spectrum of the eigenvalue problem. Using it and the fundamental results of the theory of projection methods for holomorphic Fredholm operator-valued functions [16,17], the convergence of a Nystrom method was proven in [8].

Recently, for numerical simulation of more complicated 2D microcavity lasers, namely, active cavities with piercing holes [18], a modified version of the Muller BIEs, together with a trigonometric Galerkin discretization technique, was proposed [19,20]. Mathematically, this means that there is an additional region (the hole) inside the cavity domain, and hence, an additional boundary in the integral-equation formulation. This makes the theoretical analysis more difficult compared with [8,14], as well as [15], where the problems with one boundary were investigated, as it was done originally by Muller [9]. In [21], the authors generalized results of [15] and clarified the connection between GCFEP and the eigenvalue problem for the system of Muller BIEs in this more complicated situation.

In [19,20], the authors investigated the directivities, spectra, and thresholds of the on-threshold modes of eccentric microring lasers. For such circular microcavity lasers with non-concentric circular air holes, explicit expressions for the matrix elements were obtained in [19,20]. Together with an account of the symmetry, this made the calculations much faster and more stable. Additionally, the analysis of the numerical experiments in [19,20] demonstrated the exponential convergence of the Galerkin method.

The main idea of the present work is to provide, using results of [21], a rigorous proof of convergence of the Galerkin method proposed previously in [19,20] for the numerical modeling of lasers with piercing holes, and to derive the accuracy estimates for the approximate eigenvalues and eigenfunctions. Our consideration, similar to [8], is based on the fundamental results of the theory of holomorphic operator-valued functions. Using the Galerkin method, we build a sequence of finite-dimensional holomorphic operator-valued functions that regularly approximate the original holomorphic Fredholm operator-valued function. This enables us to apply the results of the general theory to the numerical analysis of the proposed method. Particularly, we prove that if the generalized eigenfunctions are 2π-periodic and analytic in a strip of the complex plane about the real axis, then the approximate eigenvalues and eigenfunctions exponentially converge to the exact ones as the number of the basis functions increases (see Section 3 of the paper, Theorem 4, estimate ii). Thus, the numerical results of [19,20], where the exponential convergence of the Galerkin method was observed, now obtain firm mathematical ground, as for the circular boundaries the generalized eigenfunctions are infinitely smooth.

Preliminarily, in Section 2, we follow [18] and briefly recall the main steps of reducing of the original problem to the nonlinear eigenvalue problem for the set of Muller BIEs. Section 3 presents the equations related to the trigonometric Galerkin discretization of the mentioned BIEs.

To demonstrate the practical efficiency of the algorithm, in Section 4, we show that an air hole in a circular active cavity located at a certain place and with a suitable radius can lead to a notable growth in the directivity of the lasing emission together with the preservation of the low thresholds. This agrees with the physical experiments described in [22].

2. 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 ${\mathsf{\Omega }}_1$, the main body of the resonator is denoted as ${\mathsf{\Omega }}_2$, and the environment of the resonator is ${\mathsf{\Omega }}_3$. The boundaries ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ separate these regions. We suppose that the boundaries ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ are twice continuously differentiable, and $n_1$ and $n_2$ 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 ${\nu }_o$ of the hole ${\mathsf{\Omega }}_1$ and the environment around the resonator ${\mathsf{\Omega }}_3$ are given. The complex-valued refractive index of the domain ${\mathsf{\Omega }}_2$ is ${\nu }_i={\alpha }_i-i\gamma$. We denote the given real part of ${\nu }_i$ by ${\alpha }_i>0$ and the imaginary part, which is the real-valued parameter of GCFEP, by $\gamma \in ℝ$. The case of $\gamma =0$ corresponds to the passive cavity (i.e., without material losses), $\gamma <0$ is for the cavity with lossy material, and if the region ${\mathsf{\Omega }}_2$ is filled in with a gain material, then $\gamma >0$. In the latter case, the imaginary part of ${\nu }_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 $\sim \exp \left(-ikct\right)$. Herein, the speed of light in a vacuum was denoted by $c$. We were looking for complex values of $k$ on the Riemann surface $\mathbb{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\in U\backslash \left\{0\right\}$, 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 $\overline{{\mathsf{\Omega }}_1}$, $\overline{{\mathsf{\Omega }}_2}$, and $\overline{{\mathsf{\Omega }}_3}$ and twice continuously differentiable on ${\mathsf{\Omega }}_1$, ${\mathsf{\Omega }}_2$, and ${\mathsf{\Omega }}_3$.

For each $\gamma \in ℝ$, the eigenvalues $k\in \mathbb{L}$ and the eigenfunctions $u\in U\backslash \left\{0\right\}$ of GCFEP have to satisfy the Helmholtz equations,

$$\mathsf{\Delta }u+k_o^{2}u=0, x\in {\mathsf{\Omega }}_1,$$

(1)

$$\mathsf{\Delta }u+k_i^{2}u=0, x\in {\mathsf{\Omega }}_2,$$

(2)

$$\mathsf{\Delta }u+k_o^{2}u=0, x\in {\mathsf{\Omega }}_3,$$

(3)

the transmission conditions,

$$u^{-}=u^{+}, {\eta }_o\frac{\partial u^{-}}{\partial n_1}={\eta }_i\frac{\partial u^{+}}{\partial n_1}, x\in {\mathsf{\Gamma }}_1,$$

(4)

$$u^{-}=u^{+}, {\eta }_i\frac{\partial u^{-}}{\partial n_2}={\eta }_o\frac{\partial u^{+}}{\partial n_2}, x\in {\mathsf{\Gamma }}_2,$$

(5)

and the outgoing Reichardt radiation condition [23],

$$u\left(\rho ,\phi \right)={\displaystyle \sum }_{l=-\infty }^{\infty }{a}_l{H}_l^{\left(1\right)}\left(k_o\rho \right)\exp \left(il\phi \right), \rho \ge R_0.$$

(6)

Here, the polar coordinates of point $x$ are denoted by $\left(\rho ,\phi \right),$ $k_o=k{\nu }_o$, $k_i=k{\nu }_i$. In Equations (4) and (5), we have the dependence of the coefficients on the polarization; namely, ${\eta }_{o,i}={\nu }_{o,i}^{-2}$ and ${\eta }_{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_l^{\left(1\right)}\left(z\right).$ The functions $u\in U$ in (4) and (5), which are related to the boundary conditions, have the following limit values (see, e.g., [24], p. 68):

$$\frac{\partial u^{\pm }}{\partial n_i}\left(x\right)=\underset{h\to +0}{\lim }(n_i\left(x\right), \mathrm{grad} u\left(x\pm hn\left(x\right)\right), x\in {\mathsf{\Gamma }}_i, i=1,2,$$

(7)

which are expected to exist uniformly on ${\mathsf{\Gamma }}_{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 $\mathbb{L}$ by ${\mathbb{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\in {\mathbb{L}}_0$. Equation (6) is interchangeable to the common Sommerfeld radiation condition in the case of $\mathrm{Im} k=0$,

$$\left(\frac{\partial }{\partial \rho }-i{k}_o\right)u=o\left(\frac{1}{\sqrt{\rho }}\right), \rho \to \infty .$$

(8)

The case of $\mathrm{Im} k>0$ corresponds to the situation when $u$ exponentially decays as $\rho \to \infty$. The alternative case, $\mathrm{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\in \mathbb{L},\gamma \in ℝ$, and $u$, which satisfies (3) and (6):

$${{\displaystyle \int }}_{{\mathsf{\Gamma }}_R}{u}^{-}\left(y\right)\frac{\partial G_o\left(x,y\right)}{\partial n\left(y\right)}\mathrm{d}l\left(y\right)-{{\displaystyle \int }}_{{\mathsf{\Gamma }}_R}{G}_o\left(x,y\right)\frac{\partial u^{-}\left(y\right)}{\partial n\left(y\right)}\mathrm{d}l\left(y\right)=0.$$

(9)

Here, $x\in {\mathsf{\Omega }}_3$, $G_o=\left(i/4\right)H_0^{\left(1\right)}\left(k_o\left|x-y\right|\right)$. We denote ${\mathsf{\Gamma }}_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\in {\mathbb{L}}_0$ on $\gamma \in ℝ$ [8]. In the case of the passive cavity, where $\gamma \le 0$, without losses or with them, the GCFEP statement conforms with the usual statement of CFEP [12]. At this point, $\mathrm{Im} k<0$ for all the eigenvalues $k\in {\mathbb{L}}_0$. The alternative case is the active cavity, where $\gamma >0$, and the imaginary part of $k\in {\mathbb{L}}_0$ can be equal to or greater than zero. The pair $\left(k,\gamma \right)$, where $\gamma$ 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 ${\mathsf{\Omega }}_1,{\mathsf{\Omega }}_2$, and ${\mathsf{\Omega }}_3$, respectively,

$$u\left(x\right)=-{{\displaystyle \int }}_{{\mathsf{\Gamma }}_1}\frac{\partial G_o\left(x,y\right)}{\partial n_1\left(y\right)}{u}^{-}\left(y\right)\mathrm{d}l\left(y\right)+{{\displaystyle \int }}_{{\mathsf{\Gamma }}_1}{G}_o\left(x,y\right)\frac{\partial u^{-}\left(y\right)}{\partial n_1\left(y\right)}\mathrm{d}l\left(y\right), x\in {\mathsf{\Omega }}_1,$$

(10)

$$\begin{array}{c}u\left(x\right)={{\displaystyle \int }}_{{\mathsf{\Gamma }}_1}\frac{\partial G_o\left(x,y\right)}{\partial n_1\left(y\right)}{u}^{+}\left(y\right)\mathrm{d}l\left(y\right)-{{\displaystyle \int }}_{{\mathsf{\Gamma }}_1}{G}_o\left(x,y\right)\frac{\partial u^{+}\left(y\right)}{\partial n_1\left(y\right)}\mathrm{d}l\left(y\right)-\\ -{{\displaystyle \int }}_{{\mathsf{\Gamma }}_2}\frac{\partial G_i\left(x,y\right)}{\partial n_2\left(y\right)}{u}^{-}\left(y\right)\mathrm{d}l\left(y\right)+{{\displaystyle \int }}_{{\mathsf{\Gamma }}_2}{G}_i\left(x,y\right)\frac{\partial u^{-}\left(y\right)}{\partial n_2\left(y\right)}\mathrm{d}l\left(y\right), x\in {\mathsf{\Omega }}_2,\end{array}$$

(11)

$$u\left(x\right)={{\displaystyle \int }}_{{\mathsf{\Gamma }}_2}\frac{\partial G_o\left(x,y\right)}{\partial n_2\left(y\right)}{u}^{+}\left(y\right)\mathrm{d}l\left(y\right)-{{\displaystyle \int }}_{{\mathsf{\Gamma }}_2}{G}_o\left(x,y\right)\frac{\partial u^{+}\left(y\right)}{\partial n_2\left(y\right)}\mathrm{d}l\left(y\right), x\in {\mathsf{\Omega }}_3,$$

(12)

where $G_i=\left(i/4\right)H_0^{\left(1\right)}\left(k_i\left|x-y\right|\right)$. 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\in \mathbb{L}$ and $\gamma \in ℝ$ (see [8,18]). Now, we introduce the notations,

$$u_j\left(x\right)=u^{+}\left(x\right)=u^{-}\left(x\right), x\in {\mathsf{\Gamma }}_j, j=1,2,$$

(13)

$$v_1=\frac{{\eta }_i+{\eta }_o}{2{\eta }_o}\frac{\partial u^{+}}{\partial n_1}=\frac{{\eta }_i+{\eta }_o}{2{\eta }_i}\frac{\partial u^{-}}{\partial n_1}, x\in {\mathsf{\Gamma }}_1,$$

(14)

$$v_2=\frac{{\eta }_i+{\eta }_o}{2{\eta }_i}\frac{\partial u^{+}}{\partial n_2}=\frac{{\eta }_i+{\eta }_o}{2{\eta }_o}\frac{\partial u^{-}}{\partial n_2}, x\in {\mathsf{\Gamma }}_2,$$

(15)

and denote the space of continuous on ${\mathsf{\Gamma }}_j, j=1,2$, functions with the maximum norm by $C_j=C\left({\mathsf{\Gamma }}_j\right),$ $j=1,2$, $C=C_1\times C_2$, and $W=C\times 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]:

$$A\left(k,\gamma \right)w=\left(I+B\left(k,\gamma \right)\right)w=0,$$

(16)

$$B=\left(\begin{array}{cccc}{B}_1^{\left(1,1\right)}& B_1^{\left(1,2\right)}& B_1^{\left(1,3\right)}& B_1^{\left(1,4\right)}\\ B_1^{\left(2,1\right)}& B_1^{\left(2,2\right)}& B_1^{\left(2,3\right)}& B_1^{\left(2,4\right)}\\ B_2^{\left(3,1\right)}& B_2^{\left(3,2\right)}& B_2^{\left(3,3\right)}& B_2^{\left(3,4\right)}\\ B_2^{\left(4,1\right)}& B_2^{\left(4,2\right)}& B_2^{\left(4,3\right)}& B_2^{\left(4,4\right)}\end{array}\right), w=\left(\begin{array}{c}{u}_1\\ v_1\\ u_2\\ v_2\end{array}\right)$$

$$\left(B_j^{l,m}\left(k,\gamma \right)g\right)\left(x\right)={{\displaystyle \int }}_{{\mathsf{\Gamma }}_j}{K}_j^{\left(l,m\right)}\left(k,\gamma ;x,y\right)g\left(y\right)\mathrm{d}l\left(y\right).$$

Here, we denote $u_j$ or $v_j$, $j=1,2$ by the function $g$. The kernels have the following forms [18]:

$$K_j^{\left(1,1\right)}=-K_j^{\left(3,3\right)}=\frac{\partial G_o\left(x, y\right)-\partial G_i\left(x,y\right)}{\partial n_j\left(y\right)}, x\in {\mathsf{\Gamma }}_j, y\in {\mathsf{\Gamma }}_j, j=1,2,$$

$$K_j^{\left(1,2\right)}=-K_j^{\left(3,4\right)}=\frac{2\left({\eta }_o{G}_i\left(x,y\right)-{\eta }_i{G}_o\left(x,y\right)\right)}{{\eta }_i+{\eta }_o}, x\in {\mathsf{\Gamma }}_j, y\in {\mathsf{\Gamma }}_j, j=1,2,$$

$$K_1^{\left(1,3\right)}=\frac{\partial G_i\left(x,y\right)}{\partial n_2\left(y\right)}, K_1^{\left(1,4\right)}=-\frac{2{\eta }_o\partial G_i\left(x,y\right)}{{\eta }_o+{\eta }_i}, x\in {\mathsf{\Gamma }}_1, y\in {\mathsf{\Gamma }}_2,$$

$$K_j^{\left(2,1\right)}=-K_j^{\left(4,3\right)}=\frac{{\partial }^2{G}_o\left(x, y\right)}{\partial n_j\left(x\right)\partial n_j\left(y\right)}-\frac{{\partial }^2{G}_i\left(x,y\right)}{\partial n_j\left(x\right)\partial n_j\left(y\right)}, x\in {\mathsf{\Gamma }}_j, y\in {\mathsf{\Gamma }}_j, j=1,2,$$

$$K_j^{\left(2,2\right)}=-K_j^{\left(4,4\right)}=\frac{2{\eta }_o}{{\eta }_o+{\eta }_i}\frac{\partial G_i\left(x,y\right)}{\partial n_j\left(y\right)}-\frac{2{\eta }_i}{{\eta }_o+{\eta }_i}\frac{\partial G_o\left(x,y\right)}{\partial n_j\left(y\right)}, x\in {\mathsf{\Gamma }}_j, y\in {\mathsf{\Gamma }}_j, j=1,2,$$

$$K_1^{\left(2,3\right)}=\frac{{\partial }^2{G}_i\left(x,y\right)}{\partial n_1\left(x\right)\partial n_2\left(y\right)}, K_1^{\left(2,4\right)}=-\frac{2{\eta }_o}{{\eta }_o+{\eta }_i}\frac{\partial G_i\left(x,y\right)}{\partial n_1\left(y\right)}, x\in {\mathsf{\Gamma }}_1, y\in {\mathsf{\Gamma }}_2,$$

$$K_2^{\left(3,1\right)}=-\frac{\partial G_i\left(x,y\right)}{\partial n_1\left(y\right)}, K_2^{\left(3,2\right)}=\frac{2{\eta }_o{G}_i\left(x,y\right)}{{\eta }_o+{\eta }_i}, x\in {\mathsf{\Gamma }}_2, y\in {\mathsf{\Gamma }}_1,$$

$$K_2^{\left(4,1\right)}=-\frac{{\partial }^2{G}_i\left(x, y\right)}{\partial n_1\left(y\right)\partial n_2\left(x\right)}, K_2^{\left(4,2\right)}=-\frac{2{\eta }_o}{{\eta }_o+{\eta }_i}\frac{\partial G_i\left(x,y\right)}{\partial n_2\left(y\right)}, x\in {\mathsf{\Gamma }}_2, y\in {\mathsf{\Gamma }}_1.$$

Some of the kernels $K_j^{\left(q,s\right)}$ have logarithmic singularities and the others are continuous [13]. Consequently, the operator $B\left(k,\gamma \right):W\to W$ is compact, and the operator $A\left(k,\gamma \right):W\to W$ is Fredholm with index zero for every $k\in \mathbb{L}$ and $\gamma \in ℝ$ [13].

If $u\in U$ is an eigenfunction of problem (1)–(6) corresponding to an eigenvalue $k\in \mathbb{L}$ for a value of the parameter $\gamma \in ℝ$, 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\in W$ of (16) with the same values of $k$ and $\gamma$. 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 ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ term by term. However, the following result holds true (see Theorem 4 [21]). For each $\gamma \in ℝ$ and $k\in {\mathbb{I}}_{+}$ problem (16) has only the trivial solution $w=0, w\in W$. Here, ${\mathbb{I}}_{+}$ denotes the strictly positive imaginary semi-axis of ${\mathbb{L}}_0$.

3. Galerkin Method

In the current section, we present a trigonometric Galerkin method for the numerical solution of problem (16). Assume that each contour ${\mathsf{\Gamma }}_j$ has a parameterization ${\rho }_j\left(t\right)=\left({\rho }_j^1\left(t\right),{\rho }_j^2\left(t\right)\right)$, where ${\rho }_j^1\left(t\right)=f_j\left(t\right)\cos t, {\rho }_j^2\left(t\right)=f_j\left(t\right)\sin t, t\in \left[0,2\pi \right], j=1,2$. Then, for any given $\gamma \in ℝ$, we have

$$\left(B^{\left(l,m\right)}\left(k\right)w^{\left(m\right)}\right)\left(t\right)=\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }{K}^{\left(l,m\right)}\left(k;t,\tau \right)w^{\left(m\right)}\left(\tau \right)\mathrm{d}\tau .$$

Here, $l, m=1,2,3,4, y=y\left(\tau \right)\in {\mathsf{\Gamma }}_j, j=1,2,$

$$K^{\left(l,m\right)}\left(k;t,\tau \right)=2\pi K_j^{\left(l,m\right)}\left(k;x,y\right), w^{\left(m\right)}\left(\tau \right)=w^{\left(m\right)}\left(y\right)\left|{\rho }_j^{\prime }\left(\tau \right)\right|.$$

For the construction and investigation of the Galerkin method, it is convenient to consider the problem (16) in the Hilbert space $H={\left(L_2\right)}^4$, where $L_2$ denotes the space of square integrable functions with the inner product

$$\left(u,v\right)=\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }u\left(\tau \right)\overline{v\left(\tau \right)}\mathrm{d}\tau , u,v\in L_2.$$

By $T_n\subset L_2$, we denote the subspace of all trigonometric polynomials of the order no greater than $n$ with complex coefficients. Then, $H_n={\left(T_n\right)}^4\subset H$ is the subspace with the elements of the form,

$${\mathrm{w}}_n=\left(\begin{array}{c}{w}_n^{\left(1\right)}\\ w_n^{\left(2\right)}\\ w_n^{\left(3\right)}\\ w_n^{\left(4\right)}\end{array}\right), w_n^{\left(1\right)},w_n^{\left(2\right)},w_n^{\left(3\right)},w_n^{\left(4\right)}\in T_n.$$

By $p_n:H\to H_n$, we define the following projection operator:

$$p_{n}w=\left(\begin{array}{c}{\mathsf{\Phi }}_n{w}^{\left(1\right)}\\ {\mathsf{\Phi }}_n{w}^{\left(2\right)}\\ {\mathsf{\Phi }}_n{w}^{\left(3\right)}\\ {\mathsf{\Phi }}_n{w}^{\left(4\right)}\end{array}\right), w^{\left(1\right)},w^{\left(2\right)},w^{\left(3\right)},w^{\left(4\right)}\in L_2.$$

Here, ${\mathsf{\Phi }}_n:L_2\to T_n$ is the Fourier operator,

$$\left({\mathsf{\Phi }}_n{w}^{\left(m\right)}\right)\left(t\right)={\displaystyle \sum }_{q=-n}^n{c}_q\left(w^{\left(m\right)}\right){\phi }_q\left(t\right), m=1,2,3,4.$$

For $q=-n,\dots ,n$, the vectors ${\phi }_q\left(t\right)=\exp \left(iqt\right)$ form the orthonormal basis in the space $T_n$. We rewrite Equation (16) as follows

$$w^{\left(l\right)}+{\displaystyle \sum }_{m=1}^4{B}^{\left(l,m\right)}\left(k\right)w^{\left(m\right)}=0,l=1,2,3,4.$$

(17)

We look for approximate solutions $w_n^{\left(1\right)}, w_n^{\left(2\right)}, w_n^{\left(3\right)}, w_n^{\left(4\right)}\in T_n$ of the system of Equation (17) in the form

$$w_n^{\left(m\right)}\left(t\right)={\displaystyle \sum }_{q=-n}^n{\alpha }_q^{\left(m\right)}{\phi }_q\left(t\right), n\in N, m=1,2,3,4.$$

Therefore, we have

$$w_n^{\left(l\right)}+{\displaystyle \sum }_{m=1}^4{B}^{\left(l,m\right)}\left(k\right)w_n^{\left(m\right)}=0, l=1,2,3,4, n\in N.$$

We calculate the unknowns ${\alpha }_q^{\left(m\right)}$ using the Galerkin method,

$$\left(w_n^{\left(l\right)},{\phi }_p\right)+{\displaystyle \sum }_{m=1}^4(B^{\left(l,m\right)}\left(k\right)w_n^{\left(m\right)},{\phi }_p)=0, p=-n,\dots ,n,$$

(18)

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:

$${\alpha }_p^{\left(l\right)}+{\displaystyle \sum }_{m=1}^4{\displaystyle \sum }_{q=-n}^n{h}_{pq}^{\left(l,m\right)}\left(k\right){\alpha }_q^{\left(m\right)}=0, p=-n,\dots ,n,$$

(19)

where $l=1,2,3,4$,

$$h_{pq}^{\left(l,m\right)}\left(k\right)=\frac{1}{4{\pi }^2}{{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^{2\pi }{K}^{\left(l,m\right)}\left(k;t,\tau \right)\exp \left(-ipt\right)\exp \left(iq\tau \right)\mathrm{d}t\mathrm{d}\tau .$$

The system of linear algebraic Equation (19) is equivalent to the finite-dimensional linear operator equation

$$A_n\left(k\right)w_n\equiv p_{n}A\left(k\right)w_n\equiv \left(I+p_{n}B\left(k\right)\right)w_n\equiv \left(I+B_n\left(k\right)\right)w_n=0.$$

(20)

Here, $k\in L, A_n:H_n\to H_n, I$ is the unitary operator in the space $H_n$. As usual, we denote by $\rho \left(A_n\right)$ and by $\sigma \left(A_n\right)$, the regular and the characteristic sets of the operator-valued function $A_n\left(k\right)$, respectively. Let also $N^{\prime }, N^{″}, N^{‴},\dots$ be infinite sequences of the set of all natural numbers $N$.

Theorem 1.

For any given$\gamma \in ℝ$, the following statements are true:

• If$k_0$is an eigenvalue of$A\left(k\right)$, then for each$n\in N$there exists an eigenvalue$k_n$of$A_n(k$) such that$k_n\to k_0$ $\left(n\in N\right)$.
• If for each$n\in N$there exists an eigenvalue$k_n$of$A_n(k$), such that$k_n\to k_0\in \mathbb{L}$ $\left(n\in N\right)$, and$w_n$is a normalized eigenfunction of$A_n(k_n$), then

  (i)

  $k_0$is an eigenvalue of$A(k$),

  (ii)

  ${\left\{w_n\right\}}_{n\in N}$is a discretely compact sequence and its cluster points are normalized eigenfunctions of$A\left(k_0\right)$.

• For every compact$L_0\subset \rho \left(A\right)$, the sequence${\left\{A_n\left(k\right)\right\}}_{n\in N}$is stable on$L_0$, i.e., there exist$n\left(L_0\right)$and$c\left(L_0\right)$, such that$L_0\subset \rho \left(A_n\right)$,$A_n{\left(k\right)}^{-1}\le c\left(L_0\right)$for all$k\in L_0$and$n\ge n\left(L_0\right)$.

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 ${\left\{w_n\right\}}_{n\in N^{\prime }}$ of vectors from the space $H_n$ discretely converges to the limit $w\in H$ if $\Vert w_n-p_{n}w\Vert \to 0$ as $n\to \infty , n\in N^{\prime }$. Discrete convergence of the vectors will be denoted as $w_n\to w\left(n\in N^{\prime }\right)$. The sequence of elements ${\left\{w_n\right\}}_{n\in N^{\prime }}$ is called discretely compact if, for each subsequence ${\left\{w_n\right\}}_{n\in N^{″}}$, $N^{″}\subseteq N^{\prime }$, there exists a subset $N^{‴}\subseteq N^{″}$ and a vector $w\in H$, such that $w_n\to w\left(n\in N^{‴}\right)$.

Consider a bounded linear operator $A\in ℒ\left(H, H\right)$ and a sequence of finite-dimensional bounded linear operators ${\left\{A_n\right\}}_{n\in N}$. It is said that the sequence of operators ${\left\{A_n\right\}}_{n\in N}$ approximate the operator $A$, if for any vector $w\in H$ we have

$$\Vert A_n{p}_{n}w-p_{n}Aw\Vert \to 0\left(n\to \infty \right)$$

If the discrete convergence of vectors $w_n\to w\left(n\in N\right)$ implies the discrete convergence of their images, $A_n{w}_n\to Aw\left(n\in N\right)$, then the sequence of operators ${\left\{A_n\right\}}_{n\in N}$ is said to converge discretely to $A$.

The sequence of operators ${\left\{A_n\right\}}_{n\in N}$ is regular if from the boundedness of the sequence of vectors ${\left\{w_n\right\}}_{n\in N}$ (thanks to the estimate $\Vert w_n\Vert \le \mathrm{const}\left(n\in N\right))$ and from the discrete compactness of the sequence of their operator images ${\left\{A_n{w}_n\right\}}_{n\in N}$ follows the discrete compactness of the sequence of the vectors ${\left\{w_n\right\}}_{n\in N}$. If a sequence of operators ${\left\{A_n\right\}}_{n\in 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 ${\left\{A_n\left(k\right)\right\}}_{n\in N}$ regularly converges on $\mathbb{L}$ to an operator-valued function $A\left(k\right)$, if for each converging numerical sequence $k_n\to k_0\in \mathbb{L}\left(n\to \infty \right)$, the operator sequence ${\left\{A_n\left(k_n\right)\right\}}_{n\in N}$ regularly converges to the operator $A\left(k_0\right)$.

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\left(k\right):H\to H$ is holomorphic and Fredholm on $\mathbb{L}$, and its regular set is not empty. The holomorphicity and the Fredholm property of $A\left(k\right):W\to W$ were proved in Theorem 2 in [18]. For $A\left(k\right):H\to 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\in {\mathbb{I}}_{+},$ Equation (16) has only a trivial solution in the space $W$. The operator $B\left(k\right)$ is weakly singular, therefore, any solution of Equation (16) from $H$ must belong to the space $W$ and can only be trivial for $k\in {\mathbb{I}}_{+}$. The operator-valued function $A\left(k\right):H\to H$ is a Fredholm one, so for it we have ${\mathbb{I}}_{+}\subset \rho \left(A\right)$.

(b2) For any $n\in N,$ the operator-valued function $A_n\left(k\right):H_n\to H_n$ is holomorphic and Fredholm on $\mathbb{L}$. Indeed, $A\left(k\right):H\to H$ is holomorphic on $\mathbb{L}$. As the operator $p_n$ is linear and bounded, $A_n\left(k\right)=p_{n}A\left(k\right):H_n\to H_n$ has the same property. The Fredholm property of the operator-valued function $A_n\left(k\right)$ is obvious because of its finite-dimensionality.

(b3) On each compact set $L_0\subset \mathbb{L},$ the norms $\Vert A_n\left(k\right)\Vert$ are bounded uniformly in the parameters $n\in N$ and $k\in {\mathbb{L}}_0$. Indeed, from the definition of the operator $A_n\left(k\right)$ and the equality

$$\Vert p_n\Vert =1,$$

(21)

It follows that

$$\Vert A_n\left(k\right)\Vert \le \Vert A\left(k\right)\Vert , n\in N, k\in \mathbb{L},$$

However, because of

$$\Vert B^{\left(l,m\right)}{\left(k\right)\Vert }^2\le \frac{1}{4{\pi }^2}{{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^{2\pi }{\left|K^{\left(l,m\right)}\left(k;t,\tau \right)\right|}^2\mathrm{d}t\mathrm{d}\tau , l, m=1,2,3,4,$$

The following estimate is correct:

$$\Vert A\left(k\right)\Vert \le c\left(k\right), k\in \mathbb{L},$$

(22)

where $c\left(k\right)$ is a continuous function on $\mathbb{L}$:

$$c\left(k\right)=1+\frac{1}{2\pi }{\displaystyle \sum }_{l,m=1}^4{\left({{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^{2\pi }{\left|K^{\left(l,m\right)}\left(k;t,\tau \right)\right|}^2\mathrm{d}t\mathrm{d}\tau \right)}^{1/2}.$$

It is easy to see that to complete the verification of the required property, it suffices to compute the maximum of the function $c\left(k\right)$ on the given compact set $L_0\subset \mathbb{L}$.

(b4) For each fixed value $k\in \mathbb{L}$, the operator sequence ${\left\{A_n\left(k\right)\right\}}_{n\in N}$ approximates the operator $A\left(k\right)$. Indeed, by the definition of the operators $A\left(k\right):H\to H$ and $p_n:H\to H_n$, for any vector $w\in H$ we have

$$\Vert A_n\left(k\right)p_{n}w-p_{n}A\left(k\right)w\Vert =\Vert p_{n}A\left(k\right)p_{n}w-p_{n}A\left(k\right)w\Vert \le \Vert p_n\Vert \Vert A\left(k\right)\Vert \Vert p_{n}w-w\Vert \to 0\left(n\in N\right).$$

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\in \mathbb{L},$ the operator sequence ${\left\{A_n\left(k\right)\right\}}_{n\in N}$ is regular. Indeed, the discrete compactness of the sequence of vectors ${\left\{A_n\left(k\right)w_n\right\}}_{n\in N}$ means that for any $N^{\prime }\subseteq N$, there exist , such that the sequence $\left\{A_n\left(k\right)w_n=w_n+B_n\left(k\right)w_n\right\}, n\in N^{″}$ converges discretely to some $z\in H$. If the sequence ${\left\{w_n\right\}}_{n\in N^{″}}$ is bounded, then there is a weakly convergent subsequence ${\left\{w_n\right\}}_{n\in N^{‴}}$, $N^{‴}\subset N^{″}$. As it is known, the compact operator $B\left(k\right)$, takes it to a strongly converging one to some vector $u\in H$:

$$\Vert B\left(k\right)w_n-u\Vert \to 0, n\in N^{‴}$$

Hence, by virtue of the inequality

$$\Vert B_n\left(k\right)w_n-p_{n}u\Vert \le \Vert p_n\Vert \Vert B\left(k\right)w_n-u\Vert$$

and equality (21), it follows that the sequence ${\left\{B_n{w}_n\right\}}_{n\in N^{‴}}$ converges discretely to $u\in H$. Thus, ${\left\{w_n\right\}}_{n\in N^{‴}}$ converges discretely to the vector $w=z-u\in H$, and the definition of the regularity of the sequence ${\left\{A_n\left(k\right)\right\}}_{n\in 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\left(k\right)$. We denote by $G\left(A,k_0\right)$ the generalized eigenspace, i.e., the closed linear hull of all the generalized eigenfunctions of $A\left(k\right)$ corresponding to $k_0$. As the operator $p_n$ is linear, the next theorem follows from [17].

Theorem 2.

Assume that $\gamma \in ℝ$ is given, $k_0$ is an eigenvalue of $A\left(k\right)$, and $L_0\subset \mathbb{L}$ is a compact set with the boundary ${\Gamma }_0\subset \rho \left(A\right)$ so that $L_0\cap \sigma \left(A\right)=\left\{k_0\right\}$. For each $n\in N$, we denote by ${\epsilon }_n$ the maximum of the approximation error over $k\in {\Gamma }_0$ and $w\in G\left(A,k_0\right)$,

$${\epsilon }_n=\mathit{sup}\Vert \left\{A_n\left(k\right)p_{n}w-p_{n}A\left(k\right)w\Vert : k\in {\Gamma }_0, w\in G\left(A,k_0\right), \Vert w\Vert =1\right\}.$$

(23)

Then, ${\epsilon }_n\to 0$ $\left(n\in N\right)$ and the following estimations hold for almost all $n\in N$:

(i)

$\left|k_n-k_0\right|\le c{\epsilon }_n^{1/\kappa }$for all$k_n\in \sigma \left(A_n\right)\cap L_0$, where$\kappa =\kappa \left(k_0,A\right)$is the multiplicity of the pole$k_0$of the operator-valued function$A^{-1}\left(k\right)$;

(ii)

$\left|{\overline{k}}_n-k_0\right|\le c{\epsilon }_n$, where${\overline{k}}_n$is the weighted (proportionally to their algebraic multiplicities) mean of all the eigenvalues of$A_n\left(k\right)$in$L_0$,${\overline{k}}_n={\sum }_{k\in \sigma \left(A_n\right)\cap L_0}{\mu }_k\cdot k$,${\mu }_k=\nu \left(k,A_n\right)/\nu \left(k,A\right)$, where$\nu \left(\cdot ,\cdot \right)$is the algebraic multiplicity of the corresponding eigenvalue$k;$

(iii)

$max\{\left|k_n-k_0\right|:k_n\in \sigma \left(A_n\right)\cap L_0\}\le c{\epsilon }_n^{1/l_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, ${\epsilon }_n$ is defined in (23), $G_0\left(A,k_0\right)$ is the eigenspace of $A\left(k\right)$ corresponding to the eigenvalue $k_0\in \mathbb{L},$ ${\left\{k_n\right\}}_{n\in N}$ and ${\left\{w_n\right\}}_{n\in 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\to k_0$ $\left(n\in N\right)$, and ${\delta }_n$ is defined by the equality

$${\delta }_n=\mathit{max}\Vert \left\{A_n\left(k_0\right)p_n{w}_0-p_{n}A\left(k_0\right)w_0\Vert : w_0\in G_0\left(A,k_0\right), \Vert w_0\Vert =1\right\}.$$

(24)

Then, for each eigenfunction $w_n$ there exists an eigenfunction$w_0=w_0\left(w_n\right)\in G_0\left(A,k_0\right)$, such that the following error estimate holds for almost all$n\in N$:

$$\Vert w_n-w_0\Vert \le c\left({\epsilon }_n^{1/\kappa }+{\delta }_n\right).$$

Using the results from [27,28], pp. 270, 271, we derive the following approximation error estimates.

Theorem 4.

Suppose that the conditions of Theorem 2 are fulfilled, ${\epsilon }_n$ is defined in (23), ${\delta }_n$ is defined in (24), and $G_{\eta ,m}\left(0,2\pi \right)$ is the Gevrey space [28], p. 271. Then, the following error estimates are valid:

(i)

$\mathit{max}\left\{{\epsilon }_n,{\delta }_n\right\}\le c{n}^{-m}$, when${\mathrm{d}}^{i}w/\mathrm{d}{x}^i\in L_2,$ $i=0,1,\dots ,m,$for any generalized eigenfunction$w\in G\left(A,k_0\right);$

(ii)

$\mathit{max}\left\{{\epsilon }_n,{\delta }_n\right\}\le c{n}^{-m}{e}^{-\eta n}$, when$w\in G_{\eta ,m}\left(0,2\pi \right),$ $\eta >0,$ $m\ge 0,$for any generalized eigenfunction$w\in G\left(A,k_0\right)$.

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.

4. 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 $\gamma$ 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 ${\mathsf{\Omega }}_3$ and in the hole is equal ${\nu }_e=1$, the real part of the refractive index in the active region is ${\alpha }_i=2.63$, the dimensionless quantities are $\kappa =k{a}_2$, $d=\left|O_1-O_2\right|/a_2$, $r=a_1/a_2$. Here, the domains ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_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 $\left(r,d\right)$ satisfy the inequality $d+0.8r\le 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 $\left(r,d\right)$ 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).

In addition to obtaining a high value of for directivity $D$, it was necessary to maintain low values of the threshold $\gamma$. It is more convenient to search for low values of $\gamma$ by maximizing the values of the function $\mathrm{T}=-{\log }_{10}\gamma$. Figure 3 and Figure 4 show the dependences of the normalized wavenumber $\kappa$ and the threshold gain index $\gamma$ (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 $\kappa$ ranges from 5.85 to 6.3. Mostly, $\kappa$ takes on values close to 5.85 and increases only when the cavity contours are close to each other. The values of $\gamma$ remain low in the region under the straight line $d=-0.7r+0.5$ and in a small rectangle with vertices A = (0.01, 0.5), B = (0.01, 0.8), C = (0.02, 0.5), and D = (0.02, 0.8). Furthermore, above the straight line $d=-0.7r+0.5$, the values of $\gamma$ increase.

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 $\gamma$ remained low enough. For this, we considered the target functions of the following forms:

$$F=TD,$$

(25)

$$F=T+D,$$

(26)

$$F=T+10D,$$

(27)

where $T=-{\log }_{10}\gamma .$ It was assumed that among the points of the local maxima of the target functions (25)–(27), we would find such pairs of values $\left(r,d\right)$ for which a high value of the directivity $D$ and a low value of the threshold $\gamma$ 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 ${\alpha }_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. Next, by intersecting the sets of the local maxima of the considered target functions, for further research, we chose the following pairs of values $\left(r,d\right)$: (0.02 and 0.62), (0.5 and 0.15), and (0.17 and 0.5). For the function (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.

Table 1. Local maxima of the function $F=TD$.

No.	$\mathit{r}$	$\mathit{d}$	${\mathit{\beta }}_{\mathbf{r}\mathbf{a}\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{n}}$	${\mathit{\beta }}_{\mathbf{d}\mathbf{e}\mathbf{g}\mathbf{r}\mathbf{e}\mathbf{e}\mathbf{s}}$	$\mathit{T}$	$\mathit{D}$	$\mathit{F} = \mathit{T}\mathit{D}$
1	0.0247	0.6296	5.8 × 10−4	0.0334	4.6665	8.3574	39.0007
2	0.5053	0.1498	3.1410	179.9666	4.4007	5.4207	23.8555
3	0.4759	0.1065	1.7217	98.6445	5.0208	3.7432	18.7944
4	0.1766	0.4928	3.1410	179.9666	3.4664	5.3417	18.5169
5	0.0955	0.4045	2.2916	131.3014	4.9353	3.5744	17.6409
6	0.5618	0.1019	1.7217	98.6445	4.8159	3.4912	16.8138
7	0.0763	0.4486	2.8281	162.0373	4.7391	3.4036	16.1303
8	0.2616	0.2518	0.8500	48.6986	4.9647	3.2298	16.0356
9	0.4894	0.2286	3.1410	179.9666	3.4134	4.2434	14.4847

,

Table 2. Local maxima of the function $F=T+D$.

No.	$\mathit{r}$	$\mathit{d}$	${\mathit{\beta }}_{\mathbf{r}\mathbf{a}\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{n}}$	${\mathit{\beta }}_{\mathbf{d}\mathbf{e}\mathbf{g}\mathbf{r}\mathbf{e}\mathbf{e}\mathbf{s}}$	$\mathit{T}$	$\mathit{D}$	$\mathit{F} = \mathit{T} + \mathit{D}$
1	0.0264	0.6286	5.83 × 10−4	0.0334	4.5927	8.4594	13.0521
2	0.5052	0.1526	3.1410	179.9666	4.3621	5.4639	9.8260
3	0.1737	0.5020	3.1410	179.9666	3.4109	5.4120	8.8229
4	0.4765	0.1054	1.7217	98.6445	5.0289	3.7363	8.7652
5	0.0956	0.4028	2.2916	131.3014	4.9470	3.5644	8.5114
6	0.1326	0.4436	2.7163	155.6340	4.0975	4.1075	8.2050
7	0.2620	0.2471	0.8500	48.6986	4.9940	3.2060	8.1999
8	0.0607	0.4515	2.2916	131.3014	4.9077	3.2713	8.1790

and

Table 3. Local maxima of the function $F=T+10D$.

No.	$\mathit{r}$	$\mathit{d}$	${\mathit{\beta }}_{\mathbf{r}\mathbf{a}\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{n}}$	${\mathit{\beta }}_{\mathbf{d}\mathbf{e}\mathbf{g}\mathbf{r}\mathbf{e}\mathbf{e}\mathbf{s}}$	$\mathit{T}$	$\mathit{D}$	$\mathit{F} = \mathit{T} + 10\mathit{D}$
1	0.0290	0.6271	5.83 × 10−4	0.0334	4.4824	8.5175	89.6577
2	0.5042	0.1642	3.1410	179.9666	4.1982	5.5524	59.7217
3	0.1689	0.5192	3.1410	179.9666	3.2977	11.4688	58.0176
4	0.4065	0.3040	2.3587	135.1434	3.1912	4.2054	45.2453

show the values of the numerical characteristics of the cavity, and the lasing modes corresponding to the points of the local maximum. Here, $\beta$ 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.

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].

5. 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.