Constructing the scattering matrix for optical microcavities as a nonlocal boundary value problem

# Constructing the scattering matrix for optical microcavities as a nonlocal boundary value problem

## Abstract

We develop a numerical scheme to construct the scattering () matrix for optical microcavities, including the special cases with parity-time and other non-Hermitian symmetries. This scheme incorporates the explicit form of a nonlocal boundary condition, with the incident light represented by an inhomogeneous term. This approach resolves the artifact of a discontinuous normal derivative typically found in the -matrix method. In addition, we show that by excluding the aforementioned inhomogeneous term, the non-Hermitian Hamiltonian in our approach also determines the Periels-Kapur states, and it constitutes an alternative approach to derive the standard -matrix result in this basis. Therefore, our scheme provides a convenient framework to explore the benefits of both approaches. We illustrate this boundary value problem using one-dimensional and two-dimensional scalar Helmholtz equations. The eigenvalues and poles of the matrix calculated using our approach show good agreement with results obtained by other means.

OCIS codes: (140.3945) Microcavities; (290.5825) Scattering theory; (080.6755) Systems with special symmetry.

## I Introduction

Driven by advances in nanofabrication capabilities and their applications to integrated optics, understanding resonances and wave transport in optical microcavities (1); (2) has been one of the most energized subjects in modern optics. These compact optical structures also offer a unique opportunity to study non-Hermitian phenomena (3) and wave chaos (4); (5) in a well-controlled manner. To probe these properties of optical microcavities, one approach resorts to the scattering matrix formalism (6), which was an essential tool in the understanding of resonances in nuclear physics (7); (8), particle physics (9) and quantum field theory (10), which also played a crucial role in the study of wave transport in various fields including condensed matter systems (11), optics (12) and microwave networks (13).

In essence, the scattering matrix, denoted by an energy- or frequency-dependent , connects a set of incoming channels to their corresponding outgoing channels , both defined outside the scattering potential. Therefore, it takes the openness of the system into account, and the conservation of optical flux in the absence of gain and loss is manifested by the unitarity of (i.e., ). When is analytically continued into the complex- plane, its poles (i.e., where its eigenvalues approach infinity) correspond to the resonances of the system, whose wave functions only connect to the outgoing channels, now also evaluated at complex frequencies (14).

As already pointed out by Wigner and Eisenbud’s early work in nuclear physics (8), the calculation of the matrix can be understood as a nonlocal boundary value problem (BVP), which was derived using an orthogonal basis of the system and contains explicitly the real-valued frequencies of this basis. More specifically, this orthogonal basis was defined with vanishing normal derivatives at the boundary of the system, and as a result, the expansion of an arbitrary state using a finite number of these basis functions has a discontinuous normal derivative in general as an artifact (8); (15). Alternative, the expansion can be carried out using quasinormal modes (16) (i.e., the Gamow states (14)) or the Periels-Kapur states (17); (18), both defined with purely outgoing boundary conditions. These approaches, however, do not remove the artifact in the normal derivative, due to the lack of incoming flux that is inherent in the scattering process. We note that in literature the modal expansion approach, regardless of the specific basis, is referred to as the -matrix method (15) in general.

Although the consequence of the aforementioned artifact may be insignificant with the introduction of the Bloch operator (19) and a large number of basis functions (20), having an approach at one’s disposal that addresses this conceptual problem may prove valuable in some cases. In addition, it is often favorable in optical systems to adopt a self-contained approach, i.e., without requiring a priori knowledge or assigning phenomenological parameters to a large set of basis functions. These two goals can be met in principle by using either time-dependent numerical simulations of the Maxwell’s equations (21); (22) or boundary integral equations (23); (24). However, the former do not offer much physical insight on the properties of the matrix, and the latter require that the Green’s function is known inside the optical microcavity and hence do not apply, for example, when the refractive index varies smoothly due to a localized thermal source (25) or a modulated optical gain and loss profile (26).

To overcome these limitations, we propose a finite-difference approach in this article that solves for the steady state of the Maxwell’s equations inside a last scattering surface (LSS), outside which the scattered flux does not reenter other parts of the system. As we will show in Sec. I, accomplishing the two aforementioned goals in one dimension (1D) is effortless, because the boundary conditions are local and involve only the wave function and its spatial derivative at either end of the system. In higher dimensions however, a local boundary condition in the form that holds for all ’s on the LSS does not exist in general. Instead, we use a nonlocal boundary condition in its finite difference form, resulting in (1) a self-contained scheme to construct the matrix without using a modal expansion that (2) resolves the artifact in the normal derivative at the LSS. In addition, this approach via a nonlocal BVP produces the same non-Hermitian Hamiltonian that determines the Periels-Kapur states (17) (or constant-flux (CF) states (18) for the Maxwell’s equations specifically), which constitutes an alternative approach to derive the standard -matrix result in this basis. Hence our scheme provides a convenient framework to explore the benefits of both approaches when constructing and comprehending the matrix.

This paper is organized as follows: In Sec. II we discuss the BVP and -matrix approaches in 1D using the scalar Helmholtz equation. Despite the simplicity of the boundary conditions, the discussion already reveals the fundamental connection between the matrix and the CF states in our approach. We also show explicitly that the (normal) derivative of an arbitrary scattering state cannot be captured accurately by the -matrix approach with a finite number of basis functions. In Sec. III we exemplify the nonlocal BVP approach for the scattering of transverse-magnetic (TM) waves in two dimensions (2D), and the treatment of transverse-electric (TE) waves is very similar. We then apply this scheme to parity-time () (27); (28); (29); (30); (31); (32); (26); (33); (35); (36); (37); (38); (39); (34); (40); (39); (42); (43); (41); (46); (44); (45) and rotation-time () symmetric (39); (47) optical microcavities, focusing on the spontaneous symmetry breaking of the -matrix eigenvalues. Finally, we give some concluding remarks in Sec. IV.

## Ii BVP in 1D

We begin by considering the 1D Helmholtz equation

 [∂2x+ε(x)k2]Ψ(x,k)=0, (1)

where is the electric field, is the dielectric constant for an optical microcavity placed between and , and is the free-space wave vector. For an incoming wave from the left (denoted by , we write the formal solution of as

 Ψ(x,k)={Ψ−L+rLΨ+L,(x<−L/2)tLΨ+R,(x>−L/2) (2)

where are the outgoing wave functions on the left and right of the system and are the reflection and transmission coefficients. By requiring that both and are continuous at (denoted by ), the boundary conditions are then simply

 ∂xΨ|L =ik(2−Ψ|L), (3) ∂xΨ|R =ikΨ|R (4)

by eliminating and , which are local without involving both or in a single expression.

Before we embark on our quest to higher dimensions, we note an important feature of Eq. (3): without the constant term , the boundary conditions become the same as those imposed by the CF states (18), i.e., with purely outgoing waves. Similarly, an incoming wave from the right simply adds an additional constant term in Eq. (4). Therefore, starting from the non-Hermitian Hamiltonian that determines the (outgoing) CF states in the interior of the system, one can obtain the wave function in the scattering problem by turning an eigenvalue problem to an inhomogenous equation. Below we give the specific forms of this non-Hermitian Hamiltonian and the inhomogenous term using the finite-difference method.

To start, we discretize the 1D space into equally spaced points, with the left (right) boundary of the optical microcavity placed at the middle of the 0th and 1st (th and th) points. The separation of two neighboring grid points is then given by , and the Helmholtz equation takes the following form

 1Δ2[Ψi+1−2Ψi+Ψi−1]+εik2Ψi=0, (5)

where are the values of the wave function and the dielectric constant at these points. The boundary conditions (3) and (4) can be rewritten as

 Ψ0=2+iq2−iqΨ1+η, (6) ΨN+1=2+iq2−iqΨN, (7)

where is dimenionless. The constant term in Eq. (6) is due to the incoming wave as mentioned previously, and by dropping it we recover the generalized eigenvalue problem that determines the CF states with purely outgoing boundary condition (48):

 Hψm=−q2mεψm. (8)

is a column vector containing the values in the th CF state and is a diagonal matrix with elements . The CF frequencies (similar to the resonances or the poles of the matrix) is given by . Note that the real-valued free-space wave vector , instead of the complex-valued CF frequencies ’s, appears in the non-Hermitian Hamiltonian :

 Hij=[−2+2+iq2−iq(δi,1+δi,N)]δij+(δi+1,j+δi−1,j), (9)

which is tri-diagonal. We note that it is not possible to write the corresponding equation for resonances as such a generalized eigenvalue problem, because would contain the complex-valued resonance frequencies yet to be determined.

Now with the constant term in the boundary condition (6), the wave function in the scattering problem is given by

 HΨ+F=−q2εΨ, (10)

and the inhomogeneous term is a column vector with a single non-zero element, i.e., when light is incident from the left. Similarly, for the scattering of a right-incident wave, the only non-zero element of is . The equation above can be put in a more explicit form to obtain :

 Ψ=−[H+q2ε]−1F. (11)

The transmission and reflection coefficients can then be calculated using

 rL=2Ψ1−(2+iq)2−iq,tL=2ΨN2−iq (12)

for left incidence (and similarly for right incidence), and the matrix is given by

 S=(rLtLtRrR) (13)

with when the system has Lorentz reciprocity (49); (50); (51).

For comparison, here we also discuss the modal expansion approach to construct . By inserting to Eq. (10) and utilizing Eq. (8), we find

 Ψ=−ΔL∑mψmψTmq2−q2mF, (14)

where we have used the “self-orthogonality” of the CF states (52) in the following form:

 ψTmεψn=LΔδm,n. (15)

For left incidence and taking , we find

 ψ(x)=−ηΔL∑mψm(x)ψm(0)k2−k2m=2ikL∑mψm(x)ψm(0)k2−k2m. (16)

This expression is identical to that used in the standard derivation of the -matrix method in the CF basis, which we outline below.

As mentioned in the introduction, the expansion of in a finite number of basis functions introduces an artifact to the normal derivative at the LSS. Therefore, the standard derivation of the -matrix method resorts to the Green’s theorem instead to take into consideration the boundary condition of , resulting in

 Missing or unrecognized delimiter for \left (17)

We note that the scattered waves in , as well as the CF states, produce () at () after taking the derivative, and the corresponding boundary terms above are all canceled. Hence we find , with which we immediately recover Eq. (16). Once ’s and ’s are known, the matrix can be constructed again by applying Eq. (12) and the corresponding expressions for .

In Fig. 1(a) we show the total wave function inside a half-gain-half-loss optical microcavity with symmetry (27); (28); (29) and a left incident wave, whose refractive index satisfies (30); (31); (32); (26); (33); (35); (36); (38); (39); (34); (40); (42); (43); (41); (46); (44); (45); (47). Good agreement between ’s given by Eqs. (11) and (14) are obtained using 50 CF states. Nevertheless, the artifact of at the boundary of the microcavity in the -matrix approach can be readily seen in Fig. 1(b), where we plot the optical flux given by (up to a prefactor). The BVP approach, on the other hand, gives a good agreement with the analytical result (53)

 Im[Ψ∗∂xΨ]L=k(1−|rL|2),rL=G+iFD, (18)

where , , , and , , , , , .

## Iii Nonlocal BVP in 2D

In this section we elucidate how the matrix is constructed in our scheme as a nonlocal BVP in 2D. Similar to the 1D case discussed in the previous section, we show that the non-Hermitian Hamiltonian in our approach is also the one determines the CF states with purely outgoing boundary condition.

### iii.1 Construction of the S matrix

Before we proceed, we note that in the calculation of a CF state or a resonance, one basically assumes a homogenous source residing inside an optical microcavity, reflected by the imaginary part of its complex frequency. In a scattering problem, the source is an external one instead and one needs to find a way to distinguish in the boundary condition the known incident wave and the unknown scattered waves. In the 1D case shown in the previous section, the incident wave adds an inhomogeneous term (i.e., ) to the non-Hermitian eigenvalue problem that determines the outgoing CF states [c.f. Eqs.(8) and (10)]. This separation of incident and outgoing waves holds even when the boundary condition becomes nonlocal in 2D as we show below.

To illustrate this property, we again consider the scalar Helmholtz equation and use the polar coordinates. A circular LSS of radius encloses the optical microcavity with dielectric constant (see Fig. 2), outside which we assume and adopt

 Ψ−m=H−m(nekr)H−m(nekR)eimθ,Ψ+m=H+m(nek∗r)H+m(nek∗R)eimθ (19)

as our incoming and outgoing channels. Here are the second and first Hankel functions and is the angular momentum number, which also serves as the channel index. Note that we do not restrict the free-space wave vector to be real, which enables the calculation of the resonances as the complex-valued poles of the matrix.

Suppose the incident wave impinges on the LSS in the th channel . The total field and its radial derivative outside the LSS can then be written as

 Ψ>=Ψ−m0+∑mSm,m0Ψ+m, (20) ∂Ψ>∂r=V−m0(r)eim0θ+∑mSm,m0V+∗m(r)eimθ, (21)

where . Next we discretize the 2D space on a polar grid, with the circular LSS placed at the middle of the th and th rings (see Fig. 2). On each ring there are equally spaced grid points with spacing . We then write the total field on the th and th rings as

 ΨNr,ν=∑mbmeimθν, (22) ΨNr+1,ν=∑mbm(1+cmΔr)eimθν, (23)

where is the index for the grid points in the azimuthal direction, is the coresponding azimuthal angle, and is the uniform spacing between two consecutive grid points in the radial direction.

To derive the nonlocal boundary condition and the matrix, we eliminate the two coefficients in the example of TM polarization (with the electric field perpendicular to the 2D cavity plane), and the case of TE polarization can be treated in a very similar fashion. Using the continuity of both and its radial derivative, the left hand sides of Eqs. (20) and (21) on the LSS can be approximated by

 ΨNr+1,ν+ΨNr,ν2 =∑mbm(1+cmΔr2)eimθν, (24) ΨNr+1,ν−ΨNr,νΔr =∑mbmcmeimθν, (25)

 Sm,m0 =bmcm−V−m0δm,m0V+∗m, (26) Sm,m0 =bm(1+cmΔr2)−δm,m0, (27)

and we have dropped the argument in . The product can be eliminated to derive a more concise form of the matrix:

 Sm,m0=bm−(1−V−m0Δr2)δm,m01−V+∗mΔr2. (28)

This is the expression we will use in our numerical examples, which requires obtaining using the Fourier transform of ’s:

 bm=∑νΔθ2πe−imθνΨNr,ν. (29)

To find , we eliminate in Eqs. (26) and (27) and derive an expression for , which when substituted into Eq. (23) gives our nonlocal scattering boundary condition:

 ΨNr+1,ν=∑ν′Oν,ν′ΨNr,ν′+f(m0)ν, (30)

where

 Oν,ν′ ≡ Δθ2π∑mH+m(nek∗R+)H+m(nek∗R−)eim(θν−θν′), (31) f(m0)ν = (V−m0−V+∗m0)Δr1−V+∗m0Δr2eim0θν, (32)

and is the radius of the th ring. Note that the term in Eq. (30) is the manifestation of the incoming wave, and by dropping it we again recover the (nonlocal) outgoing boundary condition that defines the CF states (54), similar to Eq. (6) in 1D.

Eq. (30) can then be inserted into the discretized Helmholtz equation on the th ring to eliminate (48), which gives rise to the following matrix equation:

 (H+k2ε)ψ+F(m0)=0. (33)

The column vector represents the total wave function, i.e.,

 ψ=(ψ1,1…ψ1,Nθψ2,1…ψNr,1…ψNr,Nθ)T. (34)

and is a column vector of zeros except for the last elements, which are given by . has the same form as its 1D counterpart, i.e., a diagonal matrix with the values of the dielectric constant on the discretized grid. has rows and columns; it is the same non-Hermitian Hamiltonian that determines the CF states (54); (48):

 Hψn=−k2nεψn, (35)

and it consists a banded matrix and a block in the lower right corner. is symmetric with nonzero elements on the 0, , and diagonals

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩H(μ−1)Nθ+ν,(μ−1)Nθ+ν=−2(Δr)2−2(rμΔθ)2,H(μ−1)Nθ+ν,(μ−1)Nθ+ν+1=1(rμΔθ)2,H(μ−1)Nθ+ν,μNθ+ν=rμ+12(Δr)2√rμrμ+1,

and comes from the nonlocal boundary condition (30):

 H′ν,ν′=1(Δr)2RR−Oν,ν′. (36)

can be checked to be also symmetric using and the definition of in Eq. (30). Once is obtained for each incoming channel by solving Eq. (33), i.e.,

 ψ=−(H+k2ε)−1F(m0) (37)

as in the 1D case, we know immediately the Fourier coefficients from Eq. (29), with which the construction of the matrix is completed using Eq. (28).

In the limit of a fine grid (), Eq. (28) becomes

 Sm,m0Δr→0−−−→bm−δm,m0. (38)

We note that the second term in this expression does not depend on the dielectric constant inside the LSS, and hence it can be regarded as the result of a “direct scattering” process (6). The first term then corresponds to the “resonance-assisted” scattering process, and to understand its determining factors, we resort to the modal expansion of in the CF basis, which takes the following form in 2D:

 ψ=−ΔrΔθπR2∑nψnψTnk2−k2nF(m0). (39)

Note that we have used the following normalization of the CF basis

 ψTnεψn′=πR2ΔrΔθδn,n′ (40)

such that is dimensionless and its value does not scale with the discretization, i.e., the expression above becomes

 ∫systemεΨnΨn′rdrdϕ=πR2δn,n′ (41)

in the continuous limit.

given by Eq. (39) has the typical resonant denominator with ’s being the CF frequencies. If we project each CF state at the LSS onto the outgoing channel function (which is equivalent to a Fourier transform), i.e.,

 Ψn|r=R=∑mz(n)mΨ+m|r=R=∑mz(n)meimθ, (42)

the inner product singles out the th Fourier coefficient :

 ψTnF(m0)=2RΔrΔθ(V−m0−V+∗m0)z(n)m0, (43)

and consequently,

 Sm,m0Δr→0−−−→2R(V+∗m0−V−m0)Rm,m0−δm,m0, (44)

where

 Rm,m0=∑nz(n)mz(n)m0k2−k2n (45)

is the matrix (15). This expression indicates clearly that the contribution of a particular CF state to the scattering process is not only dependent on the resonant denominator; it also depends on the spatial overlaps between this CF state and the incoming and outgoing channels at the LSS, represented by and respectively. In the simplest example where the system is isotropic and the angular momentum is a good quantum number (e.g., a circular microdisk cavity), only the CF states with the same angular momentum as the incoming (and outgoing) channel contribute to the scattering process.

To show that given by Eq. (44) is consistent with the standard -matrix result in the CF basis, we apply the Green’s theorem to the interior of the LSS, which gives us

 an =1πR2∫systemεΨnΨdx =1πR∫LSS[Ψ∂rΨn−Ψn∂rΨ]Rdθk2−k2n. (46)

Similar to the 1D case, the outgoing channels are cancelled in the boundary integral, which can be proved rigorously by applying the Green’s theorem again to the exterior of the system, where both and satisfy

 [∇2+n2ek2]X=0,X=Ψ,Ψn. (47)

Note that the same wave vector for the total field and the CF states is crucial to remove the outgoing waves from the boundary integral in Eq. (46), and we end up with

 an=1πR∫LSS[Ψ−m0∂rΨn−Ψn∂rΨ−m0]Rdθk2−k2n (48)

for incoming wave in the th channel. Unlike the matrix in other basis (15), here the matrix does not appear in the expansion coefficient . By realizing that the expansion (42) holds not just on the LSS but also in the exterior region, we find

 ∂rΨn|r=R=∑mbmV+∗meimθ (49)

and the total wave function in the interior of the system is given by

 Ψ(x)=2R(V+∗m0−V−m0)∑nz(n)m0Ψn(x)k2−k2n. (50)

Once we substitute the on the LSS by Eq. (20) and project both sides of the equation above onto the outgoing channels, we immediately recover the matrix given by Eq. (44).

### iii.2 Results

To test our approach based on a nonlocal BVP in 2D, we first calculate the poles of the matrix for a circular microdisk cavity with a uniform index. As mentioned in the introduction, the poles of the matrix correspond to complex-valued resonances of the optical microcavity. This connection is due to the diverging eigenvalues of the matrix at its poles, meaning that an infinitesimal incoming amplitude leads to finite outgoing waves. For a circular microdisk cavity with a uniform index, the angular momentum number is a conserved quantity and the LSS is chosen as the disk boundary. The TM resonances can be found by solving the following analytical expression:

 H+′m(kR)J+′m(nkR)H+m(kR)J′m(nkR)=n. (51)

In Fig. 3 we compare the poles of the matrix calculated by the BVP approach and this analytical expression, and good agreement is found for different poles with .

Next we inspect the matrix constructed by the BVP approach from a different perspective, i.e., the symmetry property of its eigenvalues in the presence of and symmetries. It was found that in a -symmetric system ’s undergo spontaneous symmetry breaking as a function of frequency or system size (36): symmetry warrants . When in this expression are the same, the corresponding eigenstate of the matrix is in the -symmetric phase with , i.e., the optical flux in the corresponding scattering eigenstate is conserved, even though the system is non-Hermitian in the presence of gain and loss. In the -broken phase , and we find instead (or equivalently, ), which represent a pair of amplified and attenuated scattering eigenstates.

In Fig. 4 we show the eigenvalues of a microdisk cavity with refractive index , which satisfies not only symmetry (here changes to ) but also symmetry (39); (47), i.e., . The eigenvalues of the matrix for an -symmetric structure display similar spontaneous breaking as those in their -symmetric counterparts, and these properties are manifested nicely by the matrix given by Eqs. (28) using the BVP approach [Fig. 4(a)].

It can be easily seen that when a system has both and symmetries, their symmetric (broken) phases for the matrix coincide due to the unimodular property of . Therefore, the scattering eigenstates in the symmetric phase should also possess simultaneously the properties of both and symmetries. To understand and differentiate these properties, we turn to the eigenvectors of the matrix, which are the projection coefficients of the scattering eigenstates onto the incoming and outgoing channels. The cylindrical channels specified by Eq. (19) are -symmetric, i.e., , where the minus sign introduced by the parity operator (again ) in the exponent is canceled by performing a time reversal (i.e., becomes in the exponent; it also changes to ). Now if is the incident wave in an eigenstate of with eigenvalue , then by performing a combined operation on its scattered wave (i.e., ) the new incoming state should also be an eigenstate of the matrix due to symmetry. It then follows that should hold for all ’s in the -symmetric phase, and the proportional constant can be set as 1 by choosing a proper global phase of . In other words, all ’s can be made real in the -symmetric phase. If we apply a similar analysis of symmetry, we find that the channel functions transform according to , which leads to . Therefore, we find

 v−m=±(−)mvm(vm∈R) (52)

as a result of these two symmetries, which is captured nicely by the result of the BVP approach [Fig. 4(b)].

We also note that a 2D structure with both and symmetries also has mirror symmetry about the axis (i.e., the axis perpendicular to the parity axis in symmetry) (39), which imposes the following property:

 v−m=±(−)mvm(vm∈C). (53)

The overall sign corresponds to scattering eigenstates that are even and odd functions about the axis, respectively. Again the mirror symmetry about the axis is observed nicely in the BVP approach, both in the symmetry-broken phase [Fig. 4(c)] and symmetric phase [Fig. 4(d)]. In the former and hence ; In the latter we find is symmetric about both the horizontal and vertical axes instead.

The corresponding resonant modes with resonant frequencies are shown in Fig. 5, which are calculated as the poles of the matrix constructed using the BVP method. These poles differ by less than from the results of an iterative method we outline below, in both their real and imaginary parts. As briefly mentioned in Sec. II, the difference between a CF state and a resonance lies in the wave vector in the exterior of an optical microcavity: A CF state features a real-valued while a resonance has the same complex-valued resonant frequency as in the interior of the microcavity. Therefore, if we replace by a CF frequency in the eigenvalue problem (35) that determines the CF states and repeat this procedure until converges, we end up with the same complex-valued frequency in both the interior and exterior of the microcavity, which is a resonance.

It is important to note that while the resonant modes of the system have mirror symmetry about the axis here, they do not possess a - and -symmetric phase: performing a combined or operation does not leave a resonant mode unchanged. Its outgoing waves outside the microcavity are now turned into incoming waves, which by definition give a zero of the matrix, whose complex-valued frequency is the complex conjugate of the original resonance (36). Nevertheless, since the resonant modes are the scattering eigenstates at the poles of the matrix, they bare resemblance to the latter when is evaluated at a real-valued frequency, as can be seen by comparing Figs. 4(c),(d) and Figs. 5(a),(b).

## Iv Conclusion and discussion

In summary, we have presented a fintie-difference scheme to construct the matrix for optical microcavities as a BVP. The boundary condition for the total field is simple in 1D but becomes nonlocal in 2D, which appears as an inhomogeneous term and also in the non-Hermitian Hamiltonian that determines the CF states. We have verified that our approach is consistent with the -matrix method in the basis of CF states, and it addresses the artifact in the normal derivative of the total field typically found in the -matrix approach.

For applications such as enhancing light-matter interactions and sensing, often it requires accurate knowledge of wave function inside and on the boundary of optical microcavities. In such cases, the BVP method proposed here provides an economic alternative to the modal expansion approach, as the latter requires a large number of basis functions to provide the same level of accuracy. For example, at least 500 basis functions and five times more computational time are needed to capture the symmetry properties of the scattering eigenvalues shown in Fig. 4(a), whether CF states or the orthogonal states with a vanishing radial derivative at the LSS are used (55); (36).

We also note that there are several other efficient numerical approaches to construct the matrix, such as finite-different-time-domain methods (21); (22) already mentioned in the introduction and the method of auxiliary sources (56); (57). The advantages of our approach are that it provides a conceptually clear construction and a numerically straightforward implementation, and it can be applied to three dimensional structures using techniques similar to those developed for binary gratings (58). Our approach can also be applied to a network of optical microcavities (59), and it can treat continuous variations of the refractive index both inside and between these cavities, all enclosed by the LSS.

Finally, we note that while the poles of the matrix are independent of the choice of the incoming and outgoing channels, the eigenvalues of do depend on such choices in general. Only when the incoming and outgoing channels are transformed in the same way do the eigenvalues of stay unchanged, because then merely experiences a similar transformation. In our discussion of 2D TM waves, the specific forms of the channels have been chosen to simplify the notations in the derivation of the nonlocal boundary condition and the matrix. When different channels are used, for example, by changing the angular dependence of to (36), we effectively perform a permutation on the incoming channels, which is not a similar transformation with unchanged outgoing channels. Therefore, the -matrix eigenvalues and their symmetric (symmetry-broken) phases change as a result in general. Exploring this freedom of choosing the channel functions may lead to a close resemblance between the spontaneous symmetry breaking of the matrix and the corresponding close-cavity modes in - and -symmetric systems, similar to the finding in 1D heterostructures (37).

Funding. National Science Foundation (NSF) (DMR-1506987).

Acknowledgement. The author acknowledges support from NSF.

### References

1. R. K. Chang and A. J. Campillo, Optical processes in microcavities (World Scientific, 1996).
2. K. J. Vahala, Optical Microcavities (World Scientific, 2004).
3. H. Cao and J. Wiersig, “Dielectric microcavities: Model systems for wave chaos and non-Hermitian physics,” Rev. Mod. Phys. 87, 61 (2015).
4. J. U. Nöckel and A. D. Stone. “Ray and wave chaos in asymmetric resonant optical cavities,” Nature (London) 385, 45 (1997).
5. C. Gmachl, F. Capasso, E. E. Narimanov, J. U. Nöckel, A. D. Stone, J. Faist, D. L. Sivco, and A. Y. Cho. “High-power directional emission from microlasers with chaotic resonators,” Science 280, 1556 (1998).
6. W. Suh, Z. Wang, and S. Fan, “Temporal coupled-mode theory and the presence of non-orthogonal modes in lossless multimode cavities,” IEEE J. Quant. Elec. 40, 1511 (2004).
7. J. A. Wheeler, “On the mathematical description of light nuclei by the method of resonating group structure,” Phys. Rev. 52, 1107 (1937).
8. E. P. Wigner and L. Eisenbud, Phys. Rev. 72, 29 (1947).
9. Y. Nagashima. Scattering Matrix, in Elementary Particle Physics: Quantum Field Theory and Particles, Vol. 1 (Wiley-VCH, Weinheim, Germany, 2010).
10. M. E. Peskin and D. V. Schroeder, An Introduction to Quantum Field Theory (Westview, US, 1995).
11. C. W. J. Beenakker, “Random-matrix theory of quantum transport,” Rev. Mod. Phys. 69, 731 (1997).
12. R. G. Newton, Scattering Theory of Wave and Particles (Dover Publications, Mineola, New York, 2002).
13. R. H. Dicke, “A computational method applicable to microwave ntetworks,” J. Appl. Phys. 18, 873 (1947).
14. G. Gamow, “Zur quantentheorie des atomkernes,” Z. Phys. 51, 204 (1928).
15. A. M. Lane and R. G. Thomas, “R-matrix theory of nuclear reactions,” Rev. Mod. Phys. 30, 257 (1958).
16. G. Breit, “Scattering matrix of radioactive states,” Phys. Rev. 58, 1068 (1940).
17. P. L. Kapur and R. Peierls, “The dispersion formula for nuclear reactions,” Proc. Roy. Soc. A 166, 277 (1937).
18. H. E. Türeci, A. D. Stone and B. Collier, “Self-consistent multimode lasing theory for complex or random lasing media,” Phys. Rev. A 74, 043822 (2006).
19. C. Bloch, “Une formulation unifièe de la théorie des réactions nucléaires,” Nucl. Phys. 4, 503 (1957).
20. R. Szmytkowski and J. Hinze, “Convergence of the non-relativistic and relativistic R-matrix expansions at the reaction volume boundary,” J. Phys. B: At. Mol. Opt. Phys. 29, 761 (1996).
21. P. Yang and K. N. Liou, “Finite-difference time domain method for light scattering by small ice crystals in three-dimensional space,” J. Opt. Soc. Am. A 13, 2072 (1996).
22. T. Shibata and T. Itoh, “Generalized-scattering-matrix modeling of waveguide circuits using FDTD field simulations,” IEEE Trans. Microwave Theory Tech. 46, 1742 (1998).
23. J. Wiersig, “Boundary element method for resonances in dielectric microcavities,” J. Opt. A: Pure Appl. Opt. 5, 53 (2003).
24. http://homerreid.github.io/scuff-em-documentation/ (accessed July 31, 2017).
25. S. F. Liew, L. Ge, B. Redding, G. S. Solomon, and H. Cao, “Controlling a microdisk laser by local refractive index perturbation,” Appl. Phys. Lett. 108, 051105 (2016).
26. K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, “Beam dynamics in symmetric optical lattices,” Phys. Rev. Lett. 100, 103904 (2008).
27. C. M. Bender and S. Boettcher, “Real spectra in non-Hermitian Hamiltonians naving symmetry,” Phys. Rev. Lett. 80, 5243 (1998).
28. C. M. Bender, S. Boettcher, and P. N. Meisinger, “-symmetric quantum mechanics,” J. Math. Phys. 40, 2201 (1999).
29. C. M. Bender, D. C. Brody, and H. F. Jones, “Complex extension of quantum mechanics,” Phys. Rev. Lett. 89, 270401 (2002).
30. R. El-Ganainy, K. G. Makris, D. N. Christodoulides, and Z. H. Musslimani, “Theory of coupled optical -symmetric structures,” Opt. Lett. 32, 2632 (2007).
31. S. Klaiman, U. Gunther, and N. Moiseyev, “Visualization of branch points in -symmetric waveguides,” Phys. Rev. Lett. 101, 080402 (2008).
32. Z. H. Musslimani, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, “Optical solitons in periodic potentials,” Phys. Rev. Lett. 100, 030402 (2008).
33. A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides, “Observation of PT-symmetry breaking in complex optical potentials,” Phys. Rev. Lett. 103, 093902 (2009).
34. C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, “Observation of parity-time symmetry in optics,” Nature Phys. 6, 192 (2010).
35. S. Longhi, “-symmetric laser absorber,” Phys. Rev. A 82, 031801(R) (2010).
36. Y. D. Chong, L. Ge, and A. D. Stone, “-symmetry breaking and laser-absorber modes in optical scattering systems,” Phys. Rev. Lett. 106, 093902 (2011).
37. L. Ge, Y. D. Chong, and A. D. Stone, “Conservation relations and anisotropic transmission resonances in one-dimensional -symmetric photonic heterostructures,” Phys. Rev. A 85, 023802 (2012).
38. P. Ambichl, K. G. Makris, L.Ge, Y. Chong, A. D. Stone, and S. Rotter, “Breaking of symmetry in bounded and unbounded scattering systems,” Phys. Rev. X 3, 041030 (2013).
39. L. Ge and A. D. Stone, “Parity-time symmetry breaking beyond one dimension: the role of degeneracy,” Phys. Rev. X 4, 031011 (2014).
40. Z. Lin, H. Ramezani, T. Eichelkraut, T. Kottos, H. Cao, and D. N. Christodoulides, “Unidirectional invisibility induced by PT-symmetric periodic structures,” Phys. Rev. Lett. 106, 213901 (2011).
41. A. Regensburger, C. Bersch, M. A. Miri, G. Onishchukov, D. N. Christodoulides, and U. Peschel, “Parity-time synthetic photonic lattices,” Nature (London) 488, 167 (2012).
42. L. Feng, Y.-L. Xu, W. S. Fegadolli, M.-H. Lu, J. E. B. Oliveira, V. R. Almeida, Y.-F. Chen, and A. Scherer, “Experimental demonstration of a unidirectional reflectionless parity-time metamaterial at optical frequencies,” Nature Mater. 12, 108 (2013).
43. L. Feng, Z. J.Wong, R.-M.Ma, Y.Wang, and X. Zhang, “Singlemode laser by parity-time symmetry breaking,” Science 346, 972 (2014).
44. B. Peng, S. K. Özdemir, F. Lei, F. Monifi, M. Gianfreda, G. L. Long, S. Fan, F. Nori, C. M. Bender, and L. Yang, “Parity-time-symmetric whispering-gallery microcavities,” Nature Phys. 10, 394 (2014).
45. L. Chang, X. Jiang, S. Hua, C. Yang, J. Wen, L. Jiang, G. Li, G. Wang, and M. Xiao, “Parity-time symmetry and variable optical isolation in active-passive-coupled microresonators,” Nature Photon. 8, 524 (2014).
46. H. Hodaei, M. A. Miri, M. Heinrich, D. N. Christodoulides, and M. Khajavikhan, “Parity-timeâsymmetric microring lasers,” Science 346, 975 (2014).
47. L. Ge, K. G. Makris, D. N. Christodoulides, and L. Feng, “Scattering in - and -symmetric multimode waveguides: Generalized conservation laws and spontaneous symmetry breaking beyond one dimension,” Phys. Rev. A 92, 062135 (2015).
48. L. Ge, Ph. D. thesis, Yale University (2010).
49. H. A. Haus, Waves and Fields in Optoelectronics (Prentice-Hall, Englewood Cliffs, NJ, 1984), pp. 56â-61.
50. R. E. Collin, Field Theory of Guided Waves (McGraw-Hill, New York, 1960).
51. L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media (Pergamon Press, Oxford, 1960).
52. L. Ge, Y. D. Chong, and A. D. Stone, “Steady-state ab initio laser theory: generalizations and analytic results,” Phys. Rev. A 82, 063824 (2010).
53. L. Ge and L. Feng, “Optical-reciprocity-induced symmetry in photonic heterostructures and its manifestation in scattering -symmetry breaking,â Phys. Rev. A 94, 043836 (2016).
54. H. E. Türeci, L. Ge, S. Rotter and A. D. Stone, “Strong interactions in multimode random lasers,” Science 320 643 (2008).
55. Y. D. Chong, L. Ge, H. Cao, and A. D. Stone, “Coherent perfect absorbers: time-reversed lasers,” Phys. Rev. Lett. 105, 053901 (2010).
56. P. C. Waterman, “Matrix formulation of electromagnetic scattering,” Proc. IEEE 53, 805 (1965).
57. D. Maystre, S. Enoch and G. Tayeb, Scattering Matrix Method Applied to Photonic Crystals, edited by K. Yasumoto (CRC Press, 2010).
58. M. G. Moharam, Eric B. Grann, and Drew A. Pommet, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” J. Opt. Soc. Am. A 12, 1068 (1995).
59. L. Ge, “Symmetry-protected zero-mode laser with a tunable spatial profile,” Phys. Rev. A 95, 023812 (2017).
72375