Mathematical and Numerical Modeling of On-Threshold Modes of 2-D Microcavity Lasers with Piercing Holes

Abstract

This study considers the mathematical analysis framework aimed at the adequate description of the modes of lasers on the threshold of non-attenuated in time light emission. The lasers are viewed as open dielectric resonators equipped with active regions, filled in with gain material. We introduce a generalized complex-frequency eigenvalue problem for such cavities and prove important properties of the spectrum of its eigensolutions. This involves reduction of the problem to the set of the Muller boundary integral equations and their discretization with the Nystrom technique. Embedded into this general framework is the application-oriented lasing eigenvalue problem, where the real emission frequencies and the threshold gain values together form two-component eigenvalues. As an example of on-threshold mode study, we present numerical results related to the two-dimensional laser shaped as an active equilateral triangle with a round piercing hole. It is demonstrated that the threshold of lasing and the directivity of light emission, for each mode, can be efficiently manipulated with the aid of the size and, especially, the placement of the piercing hole, while the frequency of emission remains largely intact.

1. Introduction

Lasers are complicated photonic sources of light, and their understanding and design have been always closely tied to efficient and reliable modeling. Here, from the onset of research into lasers, it became clear that purely electromagnetic analysis of laser cavities is very useful. Therefore, it was also applied to the microcavity lasers, introduced in the 1990s, where the light is confined inside a transparent dielectric cavity of micrometer size. For arbitrary shapes, the mode analysis was performed usually with the aid of geometrical optics (see review papers [1,2,3,4,5] and references therein) and, for simple circular-disk shapes, with Maxwell equations (see reviews [6,7,8]). On the one hand, such analysis enabled one to explain the fact that the lasers emitted light on discrete frequencies, via the concept of natural modes of laser cavities, i.e., discrete eigenstates of the electromagnetic field as solutions to eigenvalue problems. This made the frequencies of lasing predictable, although the eigenfrequencies of open cavities are only complex-valued while lasers emit light that does not decay in time.

On the other hand, such analysis was still unable to explain another fundamental property of laser, namely, that each mode started lasing only above a certain threshold. This term reflects the fact that laser cavities differ from more conventional microwave cavities by the presence of so-called active regions. The latter is filled in with gain material, i.e., a material, which is able, at the microscopic level, to emit the light due to certain quantum-mechanical mechanisms. To enable such mechanisms to work, one has to provide an external influx of power called pumping, and the term threshold relates to the intensity of the pumping.

Therefore, to explain the threshold of lasing, the active region has to be introduced into the electromagnetic model of the laser. Then, immediately, it can be found that the eigenfrequencies of the open cavities equipped with active regions can obtain purely real values. Moreover, this happens when the gain index of the active region material takes a value, specific for each mode of the cavity. After that, it becomes evident that the classical electromagnetic eigenvalue problem should be modified to address the threshold of gain (proportional to pump power), together with the real-valued frequency, as two components of an eigenvalue. This idea nicely agrees with the experimentally observed fact that each mode of the laser cavity has a different threshold pump power, and its value is closely tied to the mode field structure.

Such observations lied in the core of the lasing eigenvalue problem (LEP) approach, suggested in [9,10,11] and applied to the on-threshold analysis of lasing modes of 2-D circular [10,11,12,13], and non-circular cavities: Limacon [14], ellipse [15], kite [16], square and other regular polygons [17]. As is easy to see, a presence of both active regions with gain material and lossy regions with absorptive material can be taken into account, in LEP, without any difficulty. Therefore, recently on-threshold mode analyses were published for the plasmonic nanolasers based on a silver strip in a quantum wire [18] and a silver tube in such a wire [19], assumed to be made of gain material. In today’s laser engineering, the direction of research associated with periodic arrays of metal or dielectric nanoparticles and nanowires, placed inside or on top of the quantum well (active layer) is very interesting and much promising. These configurations can be also treated with LEP approach. This was demonstrated by the analysis of modes of an infinite grating of circular quantum wires [20] and a binary grating made of pairs of silver and quantum wires [21]. In the latter case, it was found that the localized surface-plasmon modes have higher thresholds than the so-called lattice or grating modes, appearing thanks to periodicity. Still, from our point of view, the most impressive demonstration of the predictive power of LEP in the laser physics is the recent explanation of the mystery: Why elliptic or similar 2-D lasers emit light not on whispering-gallery modes but on bow-tie modes [15]. As we showed in that paper, the mystery is solved if one introduces a centrally located partial active region, mimicking realistic injection electrode, which better overlaps with the bow-tie mode fields than with the whispering-gallery ones.

Note that other LEP-like formulations exist [22,23,24,25,26]. They differ from what is presented by us in the way of introduction of the gain, for instance, as the imaginary part of not the refractive index but of the dielectric permittivity in active region, or a product of the imaginary part of the refractive index and the frequency.

Although the usefulness of LEP is quite obvious, its comprehensive mathematical theory has not been provided so far, however, certain efforts were presented in [27]. Development of such theory is actually the aim of our work. Here, the LEP grounding meets some difficulties. These difficulties are in the fact that the theory of operator-valued functions of a two-component vector parameter is not yet sufficiently developed. Therefore, to achieve our goal, we introduce a generalized complex-frequency eigenvalue problem (GCFEP) for the modes of open cavity with active region.

We prove the theorems on the properties of its spectrum, i.e., the set of complex eigenfrequencies, each of which depends continuously on the additional real parameter, the gain/loss index in the active/absorptive region. Then, we observe that the LEP is embedded in GCFEP as a particular case. Further, we discuss the reduction of GCFEP and LEP for 2-D lasers to the search of the spectrum of the coupled Muller boundary integral equations and discretization of the latter equations using the Nystrom algorithm. The convergence of that algorithm is proved. We illustrate the presented approach by computing the frequencies, the threshold gains, and the modal fields of the 2-D laser shaped as an active equilateral triangle with a round piercing hole. Note that the effect of piercing holes on the directionality of emission of laser was studied in [28], however, only with the aid of geometrical optics, which is sufficiently accurate only for huge-size cavities.

2. Analytical Regularization of the Generalized Complex-Frequency Eigenvalue Problem

In the beginning of this section, we consider the statement of GCFEP for 2-D optical resonators with piercing holes. The problem statement combines two physical models of mode emission from 2-D microcavities: Complex-frequency eigenvalue problem (CFEP) and LEP (see, e.g., [16,29], respectively). The geometry of investigated microcavity lasers is shown in Figure 1. Here, the domain ${\Omega }_1$ is a hole in the cavity, the main body of the resonator is ${\Omega }_2$, and the unbounded domain ${\Omega }_e$ is the environment of the resonator. These regions are separated by the boundaries ${\Gamma }_1$ and ${\Gamma }_2$. We assume that they are twice continuously differentiable curves, and n is the outer normal unit vector ether to ${\Gamma }_1$ or ${\Gamma }_2$, depending on the context.

We also assume that the refractive index ${\nu }_e$ of the environment of the resonator is positive and coincides with the refractive index ${\nu }_1$ of the piercing hole. This value is given, we write it as ${\nu }_1={\nu }_e={\alpha }_e>0$. The refractive index in the domain ${\Omega }_2$ is complex-valued, ${\nu }_2={\alpha }_a-i\gamma$. Here, ${\alpha }_a>0$ is the given real part of ${\nu }_2$, and $\gamma \in ℝ$ is the real-valued parameter of GCFEP (the loss/gain index). If the cavity is passive and lossless, then $\gamma =0$. For lossy cavities, $\gamma <0$. If the region ${\Omega }_2$ is filled in with a gain material, then $\gamma >0$.

Besides, we assume that the electromagnetic field does not depend on the variable $x_3$ and depends on time as $~exp(-ikct)$. Here, as usual, c denotes the speed of light in vacuum. We suppose that the wavenumber k is complex-valued and unknown (it is the eigenvalue of GCFEP). Following [30], we are looking for k on the Riemann surface L of the function $\ln k$. Since the electromagnetic field does not vary along the $x_3$ axis, all E- and H-components are represented in terms of a scalar eigenfunction of GCFEP, $u\in U\backslash \left\{0\right\},$ which is either the $E_3$ or $H_3$ component for E- and H-polarization, respectively. By the symbol U we denote the space of all complex-valued continuous on ${\overline{\Omega }}_1$ and ${\overline{\Omega }}_2$ and twice continuously differentiable on ${\Omega }_1$ and ${\Omega }_2$ functions.

For a given $\gamma \in ℝ$, an eigenfunction of GCFEP $u\in U\backslash \left\{0\right\}$ together with a corresponding eigenvalue $k\in L$ have to satisfy the Helmholtz equations,

$$\Delta u+k_e^{2}u=0,\hspace{1em}x\in {\Omega }_{1,e},$$

(1)

$$\Delta u+k_a^{2}u=0,\hspace{1em}x\in {\Omega }_2,$$

(2)

the transmission conditions,

$$u^{-}=u^{+},\hspace{1em}{\eta }_e\frac{\partial u^{-}}{\partial n}={\eta }_a\frac{\partial u^{+}}{\partial n},\hspace{1em}x\in {\Gamma }_1,$$

(3)

$$u^{-}=u^{+},\hspace{1em}{\eta }_a\frac{\partial u^{-}}{\partial n}={\eta }_e\frac{\partial u^{+}}{\partial n},\hspace{1em}x\in {\Gamma }_2,$$

(4)

and the outgoing Reichardt radiation condition [30,31],

$$u(r,\phi )={\displaystyle \sum _{l=-\infty }^{\infty }{a}_l}{H}_l^{(1)}\left(k_{e}r\right)e^{il\phi },\hspace{1em}r\ge R_0.$$

(5)

Here, $k_{a,e}=k{\nu }_{a,e}$ is the wavenumber in the corresponding region. The coefficients in (3) and (4) depend on polarization, namely, ${\eta }_{a,e}={\nu }_{a,e}^{-2}$ and b ${\eta }_{a,e}=1$ for the H-polarized and E-polarized fields, respectively. As usual, we denote the Hankel function of the first kind and index l by $H_l^{(1)}(z)$, and the polar coordinates of the point x by r and $\phi$. The limit values of the function $u\in U$ in boundary conditions (3) and (4) have the following definitions (see, e.g., [32], p. 68): The limits

$$\frac{\partial u^{\pm }}{\partial n}(x)=\underset{h\to +0}{\lim }(n(x),\mathrm{grad}u(x\pm hn(x)),\hspace{1em}x\in {\Gamma }_{1,2},$$

(6)

are supposed to exist uniformly on ${\Gamma }_{1,2}$. It is important to note that for any eigenfunction of GCFEP, the series in (5) converges uniformly and absolutely and is infinitely termwise differentiable [30].

Denote the principal sheet of L by $L_0$ and assume that it has a branch cut along the negative imaginary semi-axis. There are three types of eigenfunctions of GCFEP depending on the location of the corresponding eigenvalue $k\in L_0$. If $\mathrm{Im}k=0$, then (5) is equivalent to the usual Sommerfeld radiation condition,

$$\left(\frac{\partial }{\partial r}-i{k}_e\right)u=o\left(\frac{1}{\sqrt{r}}\right),\hspace{1em}r\to \infty .$$

(7)

If $\mathrm{Im}k>0$, then $u$ decays exponentially as $r\to \infty$, while, if $\mathrm{Im}k<0,$ then the eigenfunction $u$ exponentially grows at infinity. For our consideration it is important that for any $k\in L$, $\gamma \in ℝ$, and u satisfying (5), the following equality is true [30,31]:

$${\displaystyle \underset{{\Gamma }_R}{\int }{u}^{-}}(y)\frac{\partial G_e(x,y)}{\partial n(y)}dl(y)-{\displaystyle \underset{{\Gamma }_R}{\int }{G}_e}(x,y)\frac{\partial u^{-}(y)}{\partial n(y)}dl(y)=0,\hspace{1em}x\in {\Omega }_e.$$

(8)

Here, $G_e=\frac{i}{4}{H}_0^{(1)}(k_e|x-y|)$ and ${\Gamma }_R$ is a circle of sufficiently large radius R with the center at x. This fact permits us to investigate all the types of eigenfunctions in one framework.

The imaginary parts of the eigenvalues $k\in L_0$ depend on the parameter $\gamma \in ℝ$. If $\gamma \le 0$, then the cavity is passive (lossless or lossy, as it was mentioned above), and the statement of GCFEP coincides exactly with the statement of conventional CFEP [29]. In this case, as follows from the complex Poynting theorem, $\mathrm{Im}k<0$ for all the eigenvalues $k\in L_0$—see Equation (10) of [11]. If $\gamma >0$, then the cavity is active, and some of the eigenvalues $k$ can have $\mathrm{Im}k<K$, where $K>0,$ i.e., the imaginary part of some $k\in L_0$ can be equal to or greater than zero. If there exits $\gamma >0$, such that the eigenvalue $k$ of GCFEP is positive, then the pair $\left(k,\gamma \right)$ and the corresponding eigenfunction $u$ satisfy all the conditions of LEP [11,16,27]. Such values of $\gamma$ are different for different eigenvalues $k$. They are the values of the threshold material gain that is needed to compensate for the radiation losses and provide not attenuating in the time field function for the corresponding $k$.

Theorem 1. 

For each$\gamma \in ℝ$the positive imaginary semi-axis${\mathrm{Im}}_{+}$of$L_0$is free of the eigenvalues$k$of GCFEP.

Proof of Theorem 1. 

We prove this theorem by direct calculations in the same way as the uniqueness theorem for the transmission problem for the Helmholtz equations is proved (see [32], p. 100). Namely, we assume that a triple $\gamma \in ℝ$, $k\in {\mathrm{Im}}_{+}$, $u\in U$ satisfy (1)–(5) and apply Green’s theorem (see, e.g., [32], p. 68) to the function u and its complex conjugate in all the domains of the problem. Further, we use the transmission conditions (3), (4) and the condition at infinity (5). Here we note that for any $k\in {\mathrm{Im}}_{+}$ the function $u$ and all its derivatives exponentially decay at infinity, therefore, all the integrals on the domain ${\Omega }_e$ exist. Finally, equating to zero the imaginary parts of both sides of the constructed equality and analyzing they signs, we conclude that the function u is zero on the plane. Thus, u is not an eigenfunction of GCFEP. □

To apply the method of analytical regularization [33,34] for GCFEP, we use the following integral representations of the eigenfunctions $u$:

$$u(x)=-{\displaystyle \underset{{\Gamma }_1}{\int }{u}^{-}}(y)\frac{\partial G_e(x,y)}{\partial n(y)}dl(y)+{\displaystyle \underset{{\Gamma }_1}{\int }{G}_e}(x,y)\frac{\partial u^{-}(y)}{\partial n(y)}dl(y),\hspace{1em}x\in {\Omega }_1,$$

(9)

$$\begin{array}{l}u(x)={\displaystyle \underset{{\Gamma }_1}{\int }{u}^{+}}(y)\frac{\partial G_e(x,y)}{\partial n(y)}dl(y)-{\displaystyle \underset{{\Gamma }_1}{\int }{G}_e}(x,y)\frac{\partial u^{+}(y)}{\partial n(y)}dl(y)\\ \hspace{1em}-{\displaystyle \underset{{\Gamma }_2}{\int }{u}^{-}}(y)\frac{\partial G_a(x,y)}{\partial n(y)}dl(y)+{\displaystyle \underset{{\Gamma }_2}{\int }{G}_a}(x,y)\frac{\partial u^{-}(y)}{\partial n(y)}dl(y),\hspace{1em}x\in {\Omega }_2,\hspace{0.33em}\end{array}$$

(10)

$$u(x)={\displaystyle \underset{{\Gamma }_2}{\int }{u}^{+}}(y)\frac{\partial G_e(x,y)}{\partial n(y)}dl(y)-{\displaystyle \underset{{\Gamma }_2}{\int }{G}_e}(x,y)\frac{\partial u^{+}(y)}{\partial n(y)}dl(y),\hspace{1em}x\in {\Omega }_e.$$

(11)

where $G_a=\frac{i}{4}{H}_0^{(1)}(k_a|x-y|)$. Equations (9) and (10) are well known (see, e.g., [32], p. 68). Equation (11) is also true since we have (8) for any $k\in L$ and $\gamma \in ℝ$.

Let

$$u_j(x)=u^{+}(x)=u^{-}(x),\hspace{1em}x\in {\Gamma }_j,\hspace{0.33em}\hspace{0.33em}j=1,2.$$

(12)

$$v_1=\frac{{\eta }_a+{\eta }_e}{2{\eta }_e}\frac{\partial u^{+}}{\partial n}=\frac{{\eta }_a+{\eta }_e}{2{\eta }_a}\frac{\partial u^{-}}{\partial n},\hspace{1em}x\in {\Gamma }_1,\hspace{0.33em}{v}_2=\frac{{\eta }_e+{\eta }_a}{2{\eta }_a}\frac{\partial u^{+}}{\partial n}=\frac{{\eta }_e+{\eta }_a}{2{\eta }_e}\frac{\partial u^{-}}{\partial n}, \hspace{0.33em}x\in {\Gamma }_2,$$

(13)

and let $C_j=C({\Gamma }_j)$ be the Banach space of continuous functions on ${\Gamma }_j$, $j=1,2$, supplied with the maximum norm, $C=C_1\times C_2$, and $W=C\times C$. Denote by $I$ the identical operator in the space W. Then GCFEP (1)–(5) is equivalent to the following nonlinear eigenvalue problem [27] (see also [35]):

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

(14)

Here,

$$B=\left(\begin{array}{cccc}{B}_1^{1,1}& B_1^{1,2}& B_1^{1,3}& B_1^{1,4}\\ B_1^{2,1}& B_1^{2,2}& B_1^{2,3}& B_1^{2,4}\\ B_2^{3,1}& B_2^{3,2}& B_2^{3,3}& B_2^{3,4}\\ B_2^{4,1}& B_2^{4,2}& B_2^{4,3}& B_2^{4,4}\end{array}\right),w=\left(\begin{array}{c}{u}_1\\ v_1\\ u_2\\ v_2\end{array}\right),\left(B_j^{q,s}g\right)(x)={\displaystyle \underset{{\Gamma }_j}{\int }{K}_j^{q,s}}(x,y)g(y)dl(y).$$

(15)

The function $g$ denotes either $u_j$ or $v_j$, $j=1,2$. The kernels have the form [27],

$$K_1^{1,1}=-K_2^{3,3}=\frac{\partial (G_e(x,y)-G_a(x,y))}{\partial n(y)},K_1^{1,2}=-K_2^{3,4}=\frac{2({\eta }_e{G}_a(x,y)-{\eta }_a{G}_e(x,y))}{{\eta }_a+{\eta }_e},$$

(16)

$$K_1^{1,3}=-K_2^{3,1}=\frac{\partial G_a(x,y)}{\partial n(y)},\hspace{0.33em}{K}_1^{1,4}=-\frac{2{\eta }_e{G}_a(x,y)}{{\eta }_e+{\eta }_a},K_1^{2,1}=-K_2^{4,3}=\frac{{\partial }^2(G_e(x,y)-G_a(x,y))}{\partial n(x)\partial n(y)},$$

(17)

$$K_1^{2,2}=-K_2^{4,4}=\frac{2{\eta }_e}{{\eta }_a+{\eta }_e}\frac{G_a(x,y)}{\partial n(x)}-\frac{2{\eta }_a}{{\eta }_a+{\eta }_e}\frac{\partial G_e(x,y)}{\partial n(x)},K_1^{2,3}=-K_2^{4,1}=\frac{{\partial }^2{G}_a(x,y)}{\partial n(x)\partial n(y)},$$

(18)

$$K_1^{2,4}=-K_2^{4,2}=-\frac{2{\eta }_e}{{\eta }_e+{\eta }_a}\frac{\partial G_a(x,y)}{\partial n(x)},K_2^{3,2}=\frac{2{\eta }_e{G}_a(x,y)}{{\eta }_a+{\eta }_e},$$

(19)

As proved in [27], some of the kernels $K_j^{q,s}$ have the logarithmic singularities and all the others are continuous. Therefore, for each $k\in L$ and $\gamma \in ℝ$ the operator $B(k,\gamma )$ is compact [27].

Theorem 2. 

The operator$A(k,\gamma )$has a bounded inverse operator for each$\gamma \in ℝ$and$k\in I{m}_{+}$. For any given$\gamma \in ℝ$the set of all the eigenvalues$k$of the operator-valued function$A(k)$can be only a set of discrete points on$L$having finite algebraic multiplicities. Each eigenvalue$k$depends continuously on the parameter$\gamma \in ℝ$and can appear and disappear only at zero and at infinity on$L$.

Proof of Theorem 2. 

The first assertion of the theorem is derived directly from the Fredholm alternative (see, e.g. [36], p. 47), the compactness of the operator $B(k,\gamma )$, and Theorem 1. Arguing as in the proof of Lemmas 1–4 in [27] and following [37], we see that for any given $\gamma \in ℝ$ the operator-valued function $A(k)$ is holomorphic in $k\in L$. Therefore, the second assertion of the theorem follows from Proposition A. 8.4, p. 422, [38]. Following Proposition 6.1, p. 1148, [39], we prove that the operator-valued function $A(k,\gamma )$ is jointly continuous at any point $(k,\gamma )$ in $L\times ℝ$. Thus, the last assertion of the theorem follows immediately from Theorem 4.3, [40]. □

If $\gamma$ is less than or equal to zero, then the assertions of Theorem 2 correspond to CFEP, while the following corollary from Theorem 2 describes the set of all the eigenvalues of LEP.

Corollary 1. 

Assume that for some positive$\gamma$the intersection of the set of all the eigenvalues$k$of the operator-valued function$A(k)$and the positive real semi-axis of the principal shit of L is not empty. Then this intersection can be only a set of discrete eigenvalues of$A(k)$having finite algebraic multiplicities.

3. Nyström Method

Now, following [41], p. 69, we present the Nyström method for numerical solution of nonlinear eigenvalue problem (14). Note that this method was applied, in the simplest form, in [16] and then sophisticated in [17], to take full account of possible symmetry lines of a 2-D cavity. Assume that each contour ${\Gamma }_j$ has a parameterization $r_j(t)=(r_j^1(t),r_j^2(t))$, where $r_j^1(t)=f_j(t)\cos t$, $r_j^2(t)=f_j(t)\sin t$, and $t\in [0,2\pi ]$. We write

$$L(t,\tau )=Q(t,\tau )\ln \left(4{\sin }^2\frac{t-\tau }{2}\right)+P(t,\tau ),$$

(20)

where $L=K_j^{q,s}$, $P=P_j^{q,s}$, $Q=Q_j^{q,s}$, j = 1, 2, $q,s=1,2,3,4$. It is easy to see that all of the following functions are continuous:

$$Q_1^{1,1}=-Q_2^{3,3}=-\frac{\left((x-y)\cdot n(y)\right)}{|x-y|}\times \frac{k_e{J}_1(k_e|x-y|)-k_a{J}_1(k_a|x-y|)}{4\pi },$$

(21)

$$Q_1^{1,2}=-Q_2^{3,4}=\frac{{\eta }_e{J}_0(k_a|x-y|)-{\eta }_a{J}_0(k_e|x-y|)}{2\pi ({\eta }_a+{\eta }_e)},Q_1^{1,3}=-Q_2^{3,1}=-\frac{k_a{J}_1(k_a|x-y|)\left((x-y)\cdot n(y)\right)}{4\pi |x-y|},$$

(22)

$$Q_1^{1,4}=-\frac{{\eta }_e{J}_0(k_a|x-y|)}{2\pi ({\eta }_e+{\eta }_a)},Q_2^{3,2}=\frac{2{\eta }_e{J}_0(k_a|x-y|)}{4\pi ({\eta }_a+{\eta }_e)},$$

(23)

$$\begin{array}{l}{Q}_1^{2,1}=-Q_2^{4,3}=\frac{((x-y)\cdot n(y))((x-y)\cdot n(x))}{|x-y{|}^2}\times \frac{\left(k_e^2{J}_2(k_e|x-y|)-k_a^2{J}_2(k_a|x-y|)\right)}{4\pi }\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{\left(k_e{J}_1(k_e|x-y|)-k_a{J}_1(k_a|x-y|)\right)}{4\pi }\times \frac{(n(x)\cdot n(y))}{|x-y|},\end{array}$$

(24)

$$Q_1^{2,2}=-Q_1^{4,4}=\frac{\left((x-y)\cdot n(x)\right)}{|x-y|}\times \frac{{\eta }_e{k}_a{J}_1(k_a|x-y|)-{\eta }_a{k}_e{J}_1(k_e|x-y|)}{2\pi ({\eta }_a+{\eta }_e)},$$

(25)

$$Q_1^{2,3}=-Q_2^{4,1}=\frac{k_a{J}_1(k_a|x-y|)(n(x)\cdot n(y))}{4\pi |x-y|}-\frac{k_a^2{J}_2(k_a|x-y|)((x-y)\cdot n(y))((x-y)\cdot n(x))}{4\pi |x-y{|}^2},$$

(26)

$$Q_1^{2,4}=-Q_2^{4,2}=-\frac{{\eta }_e{k}_a{J}_1(k_a|x-y|)\left((x-y)\cdot n(x)\right)}{2\pi ({\eta }_a+{\eta }_e)|x-y|},$$

(27)

$$P_j^{q,s}(t,\tau )=K_j^{q,s}(t,\tau )-Q_j^{(q,s)}(t,\tau )\ln \left(4{\sin }^2\frac{t-\tau }{2}\right).$$

(28)

Let ${\mathrm{\Xi }}_h={\left\{t_{l,h}\right\}}_{l=0}^{2n-1}$ be a uniform grid on $[0,2\pi ]$ with the mesh size $h=\pi /n$, i.e., $t_l=t_{l,h}=lh$. We approximate the integrals with continuous kernels using the quadrature formula of trapezoidal rule,

$${\displaystyle {\int }_0^{2\pi }g}(\tau )d\tau \approx h{\displaystyle \sum _{l=0}^{2n-1}g}(t_l).$$

(29)

For the logarithmic parts of kernels, we use the quadrature formula for approximation of integrands by trigonometric polynomials,

$${\displaystyle {\int }_0^{2\pi }l}n\left(4{\sin }^2\frac{t-\tau }{2}\right)g(\tau )d\tau \approx {\displaystyle \sum _{l=0}^{2n-1}{R}_l^{(h)}}(t)g(t_l),$$

(30)

where

$$R_l^{(h)}(t)=-2h{\displaystyle \sum _{m=1}^{n-1}\frac{1}{m}}\cos m(t-t_l)-\frac{h}{n}\cos n(t-t_l).$$

(31)

Let

$${\mathcal{K}}_{j; l}^{(q,s)}(t)=\left(R_l^{(h)}(t)Q_j^{(q,s)}(t,t_{l,h})+h{P}_j^{(q,s)}(t,t_{l,h})\right).$$

(32)

Then for the integral operators in (14) we have

$$\left(K_j^{q,s}{g}_j\right)(x)={\displaystyle \underset{{\Gamma }_j}{\int }{K}_j^{(q,s)}}(x,y)g_j(y)dl(y)\approx {\displaystyle \sum _{l=0}^{2n-1}{\mathcal{K}}_{j; l}^{(q,s)}}(t)g_j^{(h)}(t_{l,h})|r_{j^{\prime }}(t_{l,h})|,$$

(33)

where $q,s=1,2,3,4,$ and $l=0,\dots ,2n-1.$ Applying approximations (33) and equating both sides of the functional equality following from (14) on the mesh points ${\mathrm{\Xi }}_h$, we reduce (14) to the following finite-dimensional nonlinear eigenvalue problem:

$$A_h(k,\gamma )w_h=(I+B_h(k,\gamma ))w_h=0.$$

(34)

Here, $w_h$ is the vector with the entries $w_{h,j;i}=w_j^{(h)}(t_{i,h})$. Denote by $\sigma (A)$ and $\sigma (A_h)$ the spectrum of $A(k)$ and $A_h(k)$, respectively.

Theorem 3. 

Let$\gamma \in ℝ$be given. Assume that an eigenvalue$k_0$belongs to the spectrum$\sigma (A)$. Then there exists a sequence of eigenvalues${\left\{k_h\right\}}_{h\in (0,\overline{h})}$, $k_h\in \sigma (A_h)$, such that$k_h\to k_0$as$h\to 0$. On the other hand, if${\left\{k_h\right\}}_{h\in (0,\overline{h})}$, $k_h\in \sigma (A_h)$, is a converging sequence of eigenvalues such that$k_h\to k_0\in L$as$h\to 0$, then$k_0\in \sigma (A)$. Moreover, if$k_h\in \sigma (A_h)$, $A_h(k_h)w_h=0$, where$\Vert w_h\Vert =1$, and$k_h\to k_0\in L$, $w_h\to w_0$as$h\to 0$, then$k_0\in \sigma (A)$and$A(k_0)w_0=0$, $\Vert w_0\Vert =1$.

Proof of Theorem 3. 

All the assertions of the theorem are derived directly from Theorem 2, [42]. Indeed, as it follows from Theorem 2, for any given $\gamma \in ℝ$ the operator-valued functions $A(k)$ and $A_h(k)$ are holomorphic in $k\in L.$ The first operator is Fredholm with zero index and has a bounded inverse operator for each $\gamma \in ℝ$ and $k\in {\mathrm{Im}}_{+}$, the second operator is finite-dimensional. Therefore, to conclude the proof, we have to check that the sequence ${\left\{A_h(k)\right\}}_{h\in (0,\overline{h})}$ regularly approximates $A(k)$ on L in the sense of [42]. This is easy to verify using Theorems 12.8, 12.13, pp. 202, 209, [36]. □

4. Numerical Results for LEP for Pierced Equilateral Triangle Laser

In this section, we apply the formulation presented above to the analysis of the eigenvalues of electromagnetic-field problem associated with a 2-D dielectric resonator shaped as equilateral triangle with arbitrarily located round piercing hole. Triangular-shaped 2-D lasers attract attention as possible candidates for more directive emission. The classical CFEP for triangular dielectric cavities were studied in [43,44,45] using the geometrical optics and Muller integral equations. Discretization of the latter was performed with the aid of a different technique that entailed much larger orders of matrix truncation needed to achieve comparable accuracy.

Unlike [43,44,45], we approximate the smooth “triangular” boundary ${\Gamma }_2$ using the parametric representation, originally developed by Bickley [46],

$$r_m(t)=\frac{a\sqrt{3}}{4}\left[f_m(t)-{\displaystyle \sum _{k=1}^K{\alpha }_k}{f}_m\left((3k-1)t\right)\right],\hspace{1em}m=1,\hspace{0.33em}2,{\alpha }_k={\left[3^k(3k-1)k!\right]}^{-1}{\displaystyle \prod _{j=1}^k(5-3j)},$$

(35)

where $t\in \left[0,2\pi \right],f_1(z)=\cos z,$ and $f_2(z)=\sin z.$ In computations, we take K = 5 and assume that the side of the equilateral triangle is $a=1$.

Being interested in the laser applications, we consider only the LEP for such a cavity, which is, at first, assumed uniformly active. That is, we are looking for the LEP eigenvalue pairs, $(ka,\gamma )$ and numerically solve the obtained nonlinear spectral algebraic problem by the residual inverse iteration method. In computations, we assume that the microcavity material has refractive index ${\alpha }_a=2.63$, the environment is air with ${\alpha }_e=1$, and consider the H-polarized modes.

To validate our code, we have compared its results with those of [16] for a smooth “kite” contour, where a different discretization of the same integral equation was used. The coincidence of results was observed within an arbitrary number of digits, controlled by the order of the interpolation polynomials. Note that results of [16] were supported by the experimental measurements.

Besides of the frequency and threshold gain, another important characteristic of the mode emission is the far-field directivity, which is defined as (see [16,17]),

$$D=2\pi P^{-1}{\left|\mathsf{\Phi }({\phi }_{\max })\right|}^2,$$

(36)

where

$$P={\displaystyle {\int }_0^{2\pi }{\left|\mathsf{\Phi }(\phi )\right|}^2}d\phi ,$$

(37)

$$\mathsf{\Phi }(\phi )={\displaystyle \underset{{\Gamma }_2}{\int }\left(i{k}_e\left[n(y)\cdot (\cos \phi ,\sin \phi )\right]u_2(y)+\frac{2{\eta }_2}{{\eta }_e+{\eta }_2}{v}_2(y)\right)}{e}^{-i{k}_e\left(y\cdot (\cos \phi ,\sin \phi )\right)}dl(y),\hspace{1em}\phi \in [0,2\pi ],$$

(38)

and ${\phi }_{\max }$ is the direction of the maximum intensity of radiation.

Note that for the modes with $m=0$ of circular microcavities the directivity of emission is D = 1, and for all modes with $m>0$, D = 2 [16]. Effects of the small circular holes on the directivity of emission from microcavity lasers were previously considered in [28] for stadium-like and circular cavities, however, only on the basis of the geometrical optics technique.

Figure 2a shows near and far fields (i.e., the absolute value of the magnetic field function) of the (9,9,10,o) mode normalized by the maximum value. Here, we use the notation (n₁, n₂, n₃, e) for the even with respect to the $x_1$-axis mode. The indices n₁, n₂, and n₃ correspond to the numbers of maxima of the function |u/max(u)| along the upper left-hand side, the lower left-hand side, and the right-hand side of the equilateral triangle, respectively. This mode is one of the two modes, which have minimum thresholds and the normalized frequency of lasing, $ka$, laying in between 23.5 and 26.5.

Consider now the same mode of the microcavity of the same shape but with a piercing hole bounded by the circle ${\Gamma }_1$. At first, we fix the relative radius at r/a = 0.021 and vary the position of the center $(ox,0)$ of the hole in finite interval, $ox\in \left[-0.5,\hspace{0.17em}0.25\right]$. The plots of dependences of the frequency of lasing, and the threshold gain index, and the directivity of emission on $ox$ are presented in Figure 3. These plots enable making elementary optimization of the hole position, for the given mode, with either the threshold gain index or the directivity as a target function.

Here, we indicate the values corresponding to the maximum directivity by the hollow squares. This happens at ox = 0.17. Note nearly periodic variations of the studied quantities, with the “period” close to the half-wavelength in the cavity material. Here, the minima (maxima) of the mode threshold gain correspond to the hole position at the nodes (hot spots) of the mode field without the hole. This is because the hole either does not affect or spoils, respectively, the overlap coefficient between the mode electric field and the active region [11].

Figure 4 shows the near and far fields of the (9,9,10,o) mode for the same as above triangular microcavity with holes at four different positions on the x-axis.

Consider now the same active triangle microcavity with a round hole, the center of which is fixed at the point with the Cartesian coordinates $(0,0.17)$ and suppose that the relative radius of the hole $r/a$ varies between 0.005 and 0.037. The curves in Figure 5 show that for a relatively small hole the emission frequency grows up monotonically (i.e., redshifts) with hole’s radius, in comparison with the frequency of the same mode of the cavity without a hole. The threshold gain index at first goes down slightly and reaches minimum value at $r/a=0.015.$ However, if the hole gets larger, so that its diameter approaches half-lambda in the cavity material, then the threshold grows up catastrophically. The directivity D of emission of the analyzed mode varies significantly and reaches a maximum at $r/a=0.021.$ In Figure 5, we indicate the values corresponding to this relative radius, by hollow squares. The near and far field patterns in Figure 2b correspond to the optimum hole radius $r/a=0.021$ as well to the optimum hole position ox = 0.17. As known, the lasing modes with higher directivities can be useful in many applications.

Finally, for better understanding of the mode characteristics, Figure 6 presents near and far fields of the (9,9,10,o) mode for the triangular microcavity with holes of different radiuses.

As one can see, the hole radius can also be used as a tool of the manipulation of the mode threshold gain and the directivity of emission. Still, it is less efficient than the hole location. An explanation of this fact can be seen in the mode overlap coefficient dependence on these parameters. It is apparently more important to pierce a laser cavity at the proper location (either in the mode electric field minimum or maximum) than to tune it to a certain size. In general, the issue of optimization of mode emission characteristics should be addressed separately for each mode in each specific laser.

5. Conclusions

We have presented the general framework of GCFEP and the main results of its analysis with the aid of the theory of operator-valued functions depending on parameter applying the methods of analytical regularization. We have also linked it to LEP as a full-wave classical electromagnetic theory of laser modes on threshold. Additionally, we highlighted the main steps in the conversion of GFEP and LEP for a 2-D laser with a partial active region to a set of four coupled boundary integral equations of the Muller type. We have further explained the discretization of these equations with the Nystrom quadrature method and proved its convergence. Finally, we have precisely calculated the on-threshold characteristics of a typical lasing mode of a microcavity shaped as an equilateral triangle having a circular piercing hole. In the numerical experiments, we have varied the position of the piercing hole on the x₁ axis in the cavity and the radius of the hole and measured the changes in the lasing frequencies, directionalities, and thresholds. Our numerical investigation has shown that a hole of a suitable radius and located at a certain place can lead to a notable growth of the directivity of lasing mode with the conservation of its low threshold. Hence, a small piercing hole’s radius and position in the 2-D equilateral triangular dielectric microcavity laser can be used as an engineering tool to control efficiently the directivity of emission.

We would like to emphasize that the laser model presented here possesses full mathematical rigor for the analysis of the on-threshold laser modes. This is its strength but also its limitation as it is not able to reveal time dynamics of lasing and does not account for the non-linear effects. A lesser limitation can be seen in the assumption that the material gain index in active region does not depend on the wavelength. Note that the way to overcome this limitation has been proposed in [13]. It is necessary to introduce a normalized spectral function of the gain and look for its amplitude coefficient as eigenvalue. This can be seen as a possible direction of future work. Besides, the cavity model can be refined. For instance, it is known that shaping the cavity with the aid of the molecular beam epitaxy results in the “bleaching” [15]. This can be accounted for by introducing a narrow passive region along the cavity rim and the hole rim.

Another direction of research should be the development of a numerical optimization code, based on extremely efficient code for numerical analysis, presented here. Here, a target function should be either the mode threshold gain or the directivity of its emission or a suitable combination of both. For optimization of the laser cavity shape and the active region shape, one can follow, for instance, the guidelines presented in the works [47,48].