## Abstract

We develop a numerical scheme to construct the scattering ($S$) 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 $\mathcal{R}$-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 $\mathcal{R}$-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 1D and 2D scalar Helmholtz equations. The eigenvalues and poles of the $S$ matrix calculated using our approach show good agreement with results obtained by other means.

© 2017 Chinese Laser Press

## 1. 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 (*S*) 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 *S* matrix, denoted by an energy- or frequency-dependent $S(\omega )$, connects a set of incoming channels ${\mathrm{\Psi}}^{-}$ to their corresponding outgoing channels ${\mathrm{\Psi}}^{+}$, 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 $S(\omega )$ [i.e., $S(\omega ){S}^{\u2020}(\omega )=1$]. When $S(\omega )$ is analytically continued into the complex-$\omega $ 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 $S$ matrix can be understood as a nonlocal boundary value problem (BVP), which was derived using an orthogonal basis of the system and explicitly contains 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 $\mathrm{\Psi}$ using a finite number of these basis functions has a discontinuous normal derivative in general as an artifact [8,15]. Alternatively, 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 $\mathcal{R}$-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–25]. However, the former do not offer much physical insight on the properties of the $S$ 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 [26] or a modulated optical gain and loss profile [27]. We also note that, while boundary integral equations can be regarded as a nonlocal boundary condition (and related to the Ewald–Oseen extinction theorem [23]), they typically involve solving for the wave function and its normal derivative *simultaneously* at the boundary [24].

To overcome these limitations, we propose a finite-difference approach in this paper that solves for the steady state of the Maxwell’s equations inside a last scattering surface (LSS), outside of which the scattered flux does not reenter other parts of the system. As we will show in Section 1, accomplishing the two aforementioned goals in 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 $f[\mathrm{\Psi}(\mathit{x},k),\nabla \mathrm{\Psi}(\mathit{x},k)]=0$ that holds for all $\mathit{x}$’s on the LSS does not exist in general. Instead, we use a nonlocal boundary condition $f[\mathrm{\Psi}(\mathit{x},k),{\int}_{\mathrm{LSS}}\mathrm{d}{\mathit{x}}^{\prime}O(\mathit{x},{\mathit{x}}^{\prime})\mathrm{\Psi}({\mathit{x}}^{\prime},k)]=0$ in its finite difference form, resulting in (1) a self-contained scheme to construct the $S$ 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 $\mathcal{R}$-matrix result in this basis. Hence our scheme provides a convenient framework to explore the benefits of both approaches when constructing and comprehending the $S$ matrix.

This paper is organized as follows: In Section 2 we discuss the BVP and $\mathcal{R}$-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 $S$ matrix and the CF states in our approach. We also explicitly show that the (normal) derivative of an arbitrary scattering state cannot be accurately captured by the $\mathcal{R}$-matrix approach with a finite number of basis functions. In Section 3 we exemplify the nonlocal BVP approach for the scattering of TM waves in 2D, and the treatment of TE waves is similar. We then apply this scheme to parity-time ($\mathcal{PT}$) [27–47] and rotation-time ($\mathcal{RT}$) symmetric [40,48] optical microcavities, focusing on the spontaneous symmetry breaking of the $S$-matrix eigenvalues. Finally, we give concluding remarks in Section 4.

## 2. BVP IN 1D

We begin by considering the 1D Helmholtz equation

where $\mathrm{\Psi}(x,k)$ is the electric field, $\u03f5(x)$ is the dielectric constant for an optical microcavity placed between $-L/2$ and $L/2$, and $k$ is the free-space wave vector. For an incoming wave from the left {denoted by ${\mathrm{\Psi}}_{L}^{-}(x,k)\equiv \mathrm{exp}[ik(x+L/2)]$}, we write the formal solution of $\mathrm{\Psi}(x,k)$ as*local*without involving both ${\partial}_{x}\mathrm{\Psi}{|}_{L,R}$ or $\mathrm{\Psi}{|}_{L,R}$ 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 $2ik$, the boundary conditions become the same as those imposed by the CF states [18], i.e., with purely outgoing waves. Similarly, an incoming wave ${\mathrm{\Psi}}_{R}^{-}\equiv \mathrm{exp}[-ik(x-L/2)]$ from the right simply adds an additional constant term $-2ik$ 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 inhomogeneous equation. Below we give the specific forms of this non-Hermitian Hamiltonian $H$ and the inhomogeneous term $F$ using the finite-difference method.

To start, we discretize the 1D space into $N+2$ equally spaced points, with the left (right) boundary of the optical microcavity placed at the middle of the 0th and 1st [$N$th and $(N+1)$th] points. The separation of two neighboring grid points is then given by $\mathrm{\Delta}=L/N$, and the Helmholtz equation takes the following form:

Now with the constant term $\eta $ in the boundary condition in Eq. (6), the wave function in the scattering problem is given by

and the inhomogeneous term $\mathit{F}$ is a column vector with a single non-zero element, i.e., ${F}_{1}=\eta $ when light is incident from the left. Similarly, for the scattering of a right-incident wave, the only non-zero element of $\mathit{F}$ is ${F}_{N}=\eta $. The equation above can be put into a more explicit form to obtain $\mathbf{\Psi}$:The transmission and reflection coefficients can then be calculated using

For comparison, here we also discuss the modal expansion approach to construct $S$. By inserting $\mathbf{\Psi}=\sum _{m}{a}_{m}{\mathit{\psi}}_{m}$ to Eq. (10) and utilizing Eq. (8), we find

For left incidence and taking $\mathrm{\Delta}\to 0$, we find

This expression is identical to that used in the standard derivation of the $\mathcal{R}$-matrix method in the CF basis, which we outline below.

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

We note that the scattered waves in $\mathrm{\Psi}$, as well as the CF states, produce $-ik$ ($ik$) at $x=0$ ($L$) after taking the derivative, and the corresponding boundary terms above are all canceled. Hence, we find ${[\mathrm{\Psi}{\partial}_{x}{\psi}_{m}-{\psi}_{m}{\partial}_{x}\mathrm{\Psi}]}_{0}^{L}=-[{\mathrm{\Psi}}_{L}^{-}(0){\partial}_{x}{\psi}_{m}(0)-{\psi}_{m}(0){\partial}_{x}{\mathrm{\Psi}}_{L}^{-}(0)]=2ik{\psi}_{m}(0)$, with which we immediately recover Eq. (16). Once ${\mathit{\psi}}_{m}$’s and ${q}_{m}$’s are known, the $S$ matrix can be constructed again by applying Eq. (12) and the corresponding expressions for ${r}_{R}$, ${t}_{R}$.

In Fig. 1(a) we show the total wave function inside a half-gain-half-loss optical microcavity with $\mathcal{PT}$ symmetry [28–30] and a left incident wave, whose refractive index satisfies $n(-x)={n}^{*}(x)$ [27,31–48]. Good agreement between $\mathrm{\Psi}$’s given by Eqs. (11) and (14) are obtained using 50 CF states. Nevertheless, the artifact of ${\partial}_{x}\mathrm{\Psi}$ at the boundary of the microcavity in the $\mathcal{R}$-matrix approach can be readily seen in Fig. 1(b), where we plot the optical flux given by $\mathrm{Im}[{\mathrm{\Psi}}^{*}{\partial}_{x}\mathrm{\Psi}]$ (up to a prefactor). The BVP approach, on the other hand, gives a good agreement with the analytical result [54]

## 3. NONLOCAL BVP IN 2D

In this section, we elucidate how the $S$ 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 $H$ in our approach is also the one that determines the CF states with a purely outgoing boundary condition.

#### A. Construction of the $\mathsf{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., $\mathit{F}$) 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 $R$ encloses the optical microcavity with dielectric constant $\u03f5(\mathit{x})$ (see Fig. 2), outside of which we assume $\u03f5(\mathit{x})={n}_{e}^{2}>1$ and adopt

Suppose the incident wave impinges on the LSS in the ${m}_{0}$th channel ${\mathrm{\Psi}}_{{m}_{0}}^{-}$. The total field and its radial derivative outside the LSS can then be written as

To derive the nonlocal boundary condition and the $S$ matrix, we eliminate the two coefficients ${b}_{m}$, ${c}_{m}$ 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 similar fashion. Using the continuity of both $\mathrm{\Psi}$ and its radial derivative, the left-hand sides of Eqs. (20) and (21) on the LSS can be approximated by

This is the expression we will use in our numerical examples, which requires obtaining ${b}_{m}$ using the Fourier transform of ${\mathrm{\Psi}}_{{N}_{r},\nu}$’s:

To find ${\mathrm{\Psi}}_{{N}_{r},\nu}$, we eliminate ${S}_{m,{m}_{0}}$ in Eqs. (26) and (27) and derive an expression for ${b}_{m}{c}_{m}$, which when substituted into Eq. (23) gives our nonlocal scattering boundary condition:

Equation (30) can then be inserted into the discretized Helmholtz equation on the $N$th ring to eliminate ${\mathrm{\Psi}}_{{N}_{r}+1,\nu}$ [49], which gives rise to the following matrix equation:

The column vector $\mathit{\psi}\equiv \sqrt{r}\mathbf{\Psi}$ represents the total wave function, i.e.,

In the limit of a fine grid (${\mathrm{\Delta}}_{r}\to 0$), Eq. (28) becomes

We note that the second term in this expression does not depend on the dielectric constant inside the LSS; 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 $\mathit{\psi}$ in the CF basis, which takes the following form in 2D:$\mathit{\psi}$ given by Eq. (40) has the typical resonant denominator with ${k}_{n}$’s being the CF frequencies. If we project each CF state at the LSS onto the outgoing channel function ${\mathrm{\Psi}}_{m}^{+}$ (which is equivalent to a Fourier transform), i.e.,

To show that ${S}_{m,{m}_{0}}$ given by Eq. (45) is consistent with the standard $\mathcal{R}$-matrix result in the CF basis, we apply the Green’s theorem to the interior of the LSS, which gives us

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 $\mathrm{\Psi}$ and ${\psi}_{n}$ satisfy

Note that the same wave vector $k$ for the total field $\mathrm{\Psi}$ and the CF states ${\psi}_{n}$ is crucial to remove the outgoing waves from the boundary integral in Eq. (47), and we end up with

Once we substitute the $\mathrm{\Psi}(\mathit{x})$ on the LSS by Eq. (20) and project both sides of the equation above onto the outgoing channels, we immediately recover the $S$ matrix given by Eq. (45).

#### B. Results

To test our approach based on a nonlocal BVP in 2D, we first calculate the poles of the $S$ matrix for a circular microdisk cavity with a uniform index. As mentioned in the introduction, the poles of the $S$ matrix correspond to complex-valued resonances of the optical microcavity. This connection is due to the diverging eigenvalues of the $S$ 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 $m$ 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:

In Fig. 3 we compare the poles of the $S$ matrix calculated by the BVP approach and this analytical expression, and good agreement is found for different poles with $m\in [\mathrm{0,12}]$.

Next, we inspect the $S$ matrix constructed by the BVP approach from a different perspective, i.e., the symmetry property of its eigenvalues in the presence of $\mathcal{PT}$ and $\mathcal{RT}$ symmetries. It was found that, in a $\mathcal{PT}$-symmetric system, ${s}_{m}$’s undergo spontaneous symmetry breaking as a function of frequency or system size [37]: $\mathcal{PT}$ symmetry warrants ${s}_{m}^{*}={s}_{{m}^{\prime}}^{-1}$. When $m,{m}^{\prime}$ in this expression are the same, the corresponding eigenstate of the $S$ matrix is in the $\mathcal{PT}$-symmetric phase with $|{s}_{m}|=1$, 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 $\mathcal{PT}$-broken phase $m\ne {m}^{\prime}$, and we find $|{s}_{m}|={|{s}_{{m}^{\prime}}|}^{-1}$ instead (or equivalently, ${\mathrm{log}}_{10}{|{s}_{m}|}^{2}=-{\mathrm{log}}_{10}{|{s}_{{m}^{\prime}}|}^{2}$), which represent a pair of amplified and attenuated scattering eigenstates.

In Fig. 4 we show the eigenvalues of a microdisk cavity with refractive index $n(\mathit{x})=1.5+0.4\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\theta $, which satisfies not only $\mathcal{PT}$ symmetry (here $\mathcal{P}$ changes $\theta $ to $-\theta $) but also $\mathcal{RT}$ symmetry [40,48], i.e., $n(r,\theta )={n}^{*}(r,\theta +\pi )$. The eigenvalues of the $S$ matrix for an $\mathcal{RT}$-symmetric structure display similar spontaneous breaking as those in their $\mathcal{PT}$-symmetric counterparts, and these properties are nicely manifested by the $S$ matrix given by Eq. (28) using the BVP approach [Fig. 4(a)].

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

as a result of these two symmetries, which is nicely captured by the result of the BVP approach [Fig. 4(b)].We also note that a 2D structure with both $\mathcal{PT}$ and $\mathcal{RT}$ symmetries also has mirror symmetry about the $\theta =\pm \pi /2$ axis (i.e., the axis perpendicular to the parity axis in $\mathcal{PT}$ symmetry) [40], which imposes the following property:

The overall $\pm $ sign corresponds to scattering eigenstates that are even and odd functions about the $\theta =\pm \pi /2$ axis, respectively. Again, the mirror symmetry about the $\theta =\pm \pi /2$ axis is nicely observed in the BVP approach, both in the symmetry-broken phase [Fig. 4(c)] and symmetric phase [Fig. 4(d)]. In the former $\mathcal{PT}\mathrm{\Psi}\ne \mathrm{\Psi}$, $\mathcal{RT}\mathrm{\Psi}\ne \mathrm{\Psi}$ and hence $|\mathrm{\Psi}(r,\theta )|\ne |\mathrm{\Psi}(r,-\theta )|,|\mathrm{\Psi}(r,\theta )|\ne |\mathrm{\Psi}(r,\theta +\pi )|$; in the latter, we find $|\mathrm{\Psi}|$ is symmetric about both the horizontal and vertical axes instead.

The corresponding resonant modes with resonant frequencies ${k}_{m}R=3.5042-0.2492i,4.5989-0.6445i$ are shown in Fig. 5, which are calculated as the poles of the $S$ matrix constructed using the BVP method. These poles differ by less than 0.01% from the results of an iterative method we outline below, in both their real and imaginary parts. As briefly mentioned in Section 2, 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 $k$, while a resonance has the same complex-valued resonant frequency ${k}_{m}$ as in the interior of the microcavity. Therefore, if we replace $k$ by a CF frequency ${k}_{m}$ in the eigenvalue problem in Eq. (35), which determines the CF states, and repeat this procedure until ${k}_{m}$ converges, we end up with the same complex-valued frequency ${k}_{m}$ 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 $\theta =\pm \pi /2$ axis, here, they do not possess a $\mathcal{PT}$- and $\mathcal{RT}$-symmetric phase: performing a combined $\mathcal{PT}$ or $\mathcal{RT}$ 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 $S$ matrix, whose complex-valued frequency is the complex conjugate of the original resonance [37]. Nevertheless, because the resonant modes are the scattering eigenstates at the poles of the $S$ matrix, they bear resemblance to the latter when $S$ is evaluated at a real-valued frequency, as can be seen by comparing Figs. 4(c) and 4(d) and Figs. 5(a) and 5(b).

## 4. CONCLUSION AND DISCUSSION

In summary, we have presented a finite-difference scheme to construct the $S$ 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. Although the uniqueness and existence of the solutions for a nonlocal BVP are not guaranteed in general, here we can rest assured as the underlying physical process (i.e., scattering) is deterministic. We have verified that our approach is consistent with the $\mathcal{R}$-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 $\mathcal{R}$-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 [37,57].

We also note that there are several other efficient numerical approaches to construct the $S$ matrix, such as finite-different-time-domain methods [21,22] already mentioned in the introduction and the method of auxiliary sources [58,59]. The advantages of our approach are that it provides a conceptually clear construction and a numerically straightforward implementation, and it can be applied to 3D structures using techniques similar to those developed for binary gratings [60]. Our approach can also be applied to a network of optical microcavities [61], 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 $S$ matrix are independent of the choice of the incoming and outgoing channels, the eigenvalues of $S$ do depend on such choices in general. Only when the incoming and outgoing channels are transformed in the same way do the eigenvalues of $S$ stay unchanged because then $S$ 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 $S$ matrix. When different channels are used, for example, by changing the angular dependence of ${\mathrm{\Psi}}_{m}^{-}$ to ${e}^{-im\theta}$ [37], we effectively perform a permutation on the incoming channels, which is not a similar transformation with unchanged outgoing channels. Therefore, the $S$-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 $S$ matrix and the corresponding close-cavity modes in $\mathcal{PT}$- and $\mathcal{RT}$-symmetric systems, similar to the finding in 1D heterostructures [38].

## Funding

Directorate for Mathematical and Physical Sciences (MPS) (DMR-1506987); National Science Foundation (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–111 (2015). [CrossRef]

**4. **J. U. Nöckel and A. D. Stone, “Ray and wave chaos in asymmetric resonant optical cavities,” Nature (London) **385**, 45–47 (1997). [CrossRef]

**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–1564 (1998). [CrossRef]

**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. Quantum Electron. **40**, 1511–1518 (2004). [CrossRef]

**7. **J. A. Wheeler, “On the mathematical description of light nuclei by the method of resonating group structure,” Phys. Rev. **52**, 1107–1122 (1937). [CrossRef]

**8. **E. P. Wigner and L. Eisenbud, “Higher angular momenta and long range interaction in resonance reactions,” Phys. Rev. **72**, 29–41 (1947). [CrossRef]

**9. **Y. Nagashima, *Scattering Matrix, in Elementary Particle Physics: Quantum Field Theory and Particles* (Wiley-VCH, 2010), Vol. 1.

**10. **M. E. Peskin and D. V. Schroeder, *An Introduction to Quantum Field Theory* (Westview, 1995).

**11. **C. W. J. Beenakker, “Random-matrix theory of quantum transport,” Rev. Mod. Phys. **69**, 731–808 (1997). [CrossRef]

**12. **R. G. Newton, *Scattering Theory of Wave and Particles* (Dover, 2002).

**13. **R. H. Dicke, “A computational method applicable to microwave networks,” J. Appl. Phys. **18**, 873–878 (1947). [CrossRef]

**14. **G. Gamow, “Zur quantentheorie des atomkernes,” Z. Phys. **51**, 204–212 (1928). [CrossRef]

**15. **A. M. Lane and R. G. Thomas, “R-matrix theory of nuclear reactions,” Rev. Mod. Phys. **30**, 257–353 (1958). [CrossRef]

**16. **G. Breit, “Scattering matrix of radioactive states,” Phys. Rev. **58**, 1068–1074 (1940). [CrossRef]

**17. **P. L. Kapur and R. Peierls, “The dispersion formula for nuclear reactions,” Proc. R. Soc. A **166**, 277–295 (1938). [CrossRef]

**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). [CrossRef]

**19. **C. Bloch, “Une formulation unifièe de la théorie des réactions nucléaires,” Nucl. Phys. **4**, 503–528 (1957). [CrossRef]

**20. **R. Szmytkowski and J. Hinze, “Convergence of the non-relativistic and relativistic R-matrix expansions at the reaction volume boundary,” J. Phys. B **29**, 761–777 (1996). [CrossRef]

**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–2085 (1996). [CrossRef]

**22. **T. Shibata and T. Itoh, “Generalized-scattering-matrix modeling of waveguide circuits using FDTD field simulations,” IEEE Trans. Microwave Theory Tech. **46**, 1742–1751 (1998). [CrossRef]

**23. **D. N. Pattanayak and E. Wolf, “Scattering states and bound states as solutions of the Schrodinger equation with nonlocal boundary conditions,” Phys. Rev. D **13**, 913–923 (1976). [CrossRef]

**24. **J. Wiersig, “Boundary element method for resonances in dielectric microcavities,” J. Opt. A **5**, 53–60 (2003). [CrossRef]

**25. **http://homerreid.github.io/scuff-em-documentation/ (accessed July 31, 2017).

**26. **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). [CrossRef]

**27. **K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, “Beam dynamics in PT symmetric optical lattices,” Phys. Rev. Lett. **100**, 103904 (2008). [CrossRef]

**28. **C. M. Bender and S. Boettcher, “Real spectra in non-Hermitian Hamiltonians naving PT symmetry,” Phys. Rev. Lett. **80**, 5243–5246 (1998). [CrossRef]

**29. **C. M. Bender, S. Boettcher, and P. N. Meisinger, “PT-symmetric quantum mechanics,” J. Math. Phys. **40**, 2201–2229 (1999). [CrossRef]

**30. **C. M. Bender, D. C. Brody, and H. F. Jones, “Complex extension of quantum mechanics,” Phys. Rev. Lett. **89**, 270401 (2002). [CrossRef]

**31. **R. El-Ganainy, K. G. Makris, D. N. Christodoulides, and Z. H. Musslimani, “Theory of coupled optical PT-symmetric structures,” Opt. Lett. **32**, 2632–2634 (2007). [CrossRef]

**32. **S. Klaiman, U. Gunther, and N. Moiseyev, “Visualization of branch points in PT-symmetric waveguides,” Phys. Rev. Lett. **101**, 080402 (2008). [CrossRef]

**33. **Z. H. Musslimani, K. G. Makris, R. El-Ganainy, and D. N. Christodoulides, “Optical solitons in PT periodic potentials,” Phys. Rev. Lett. **100**, 030402 (2008). [CrossRef]

**34. **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). [CrossRef]

**35. **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,” Nat. Phys. **6**, 192–195 (2010). [CrossRef]

**36. **S. Longhi, “PT-symmetric laser absorber,” Phys. Rev. A **82**, 031801 (2010). [CrossRef]

**37. **Y. D. Chong, L. Ge, and A. D. Stone, “PT-symmetry breaking and laser-absorber modes in optical scattering systems,” Phys. Rev. Lett. **106**, 093902 (2011). [CrossRef]

**38. **L. Ge, Y. D. Chong, and A. D. Stone, “Conservation relations and anisotropic transmission resonances in one-dimensional PT-symmetric photonic heterostructures,” Phys. Rev. A **85**, 023802 (2012). [CrossRef]

**39. **P. Ambichl, K. G. Makris, L. Ge, Y. Chong, A. D. Stone, and S. Rotter, “Breaking of PT symmetry in bounded and unbounded scattering systems,” Phys. Rev. X **3**, 041030 (2013). [CrossRef]

**40. **L. Ge and A. D. Stone, “Parity-time symmetry breaking beyond one dimension: the role of degeneracy,” Phys. Rev. X **4**, 031011 (2014). [CrossRef]

**41. **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). [CrossRef]

**42. **A. Regensburger, C. Bersch, M. A. Miri, G. Onishchukov, D. N. Christodoulides, and U. Peschel, “Parity-time synthetic photonic lattices,” Nature (London) **488**, 167–171 (2012). [CrossRef]

**43. **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,” Nat. Mater. **12**, 108–113 (2013). [CrossRef]

**44. **L. Feng, Z. J. Wong, R.-M. Ma, Y. Wang, and X. Zhang, “Singlemode laser by parity-time symmetry breaking,” Science **346**, 972–975 (2014). [CrossRef]

**45. **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,” Nat. Phys. **10**, 394–398 (2014). [CrossRef]

**46. **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,” Nat. Photonics **8**, 524–529 (2014). [CrossRef]

**47. **H. Hodaei, M. A. Miri, M. Heinrich, D. N. Christodoulides, and M. Khajavikhan, “Parity-time symmetric microring lasers,” Science **346**, 975–978 (2014). [CrossRef]

**48. **L. Ge, K. G. Makris, D. N. Christodoulides, and L. Feng, “Scattering in PT- and RT-symmetric multimode waveguides: Generalized conservation laws and spontaneous symmetry breaking beyond one dimension,” Phys. Rev. A **92**, 062135 (2015). [CrossRef]

**49. **L. Ge, “Steady-state ab initio laser theory and its applications in random and complex media,” Ph.D. thesis (Yale University, 2010).

**50. **H. A. Haus, *Waves and Fields in Optoelectronics* (Prentice-Hall, 1984), pp. 56–61.

**51. **R. E. Collin, *Field Theory of Guided Waves* (McGraw-Hill, 1960).

**52. **L. D. Landau and E. M. Lifshitz, *Electrodynamics of Continuous Media* (Pergamon, 1960).

**53. **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). [CrossRef]

**54. **L. Ge and L. Feng, “Optical-reciprocity-induced symmetry in photonic heterostructures and its manifestation in scattering PT-symmetry breaking,” Phys. Rev. A **94**, 043836 (2016). [CrossRef]

**55. **S. Rotter, J.-Z. Tang, L. Wirtz, J. Trost, and J. Burgdörfer, “Modular recursive Greens function method for ballistic quantum transport,” Phys. Rev. B **62**, 1950–1960 (2000). [CrossRef]

**56. **H. E. Türeci, L. Ge, S. Rotter, and A. D. Stone, “Strong interactions in multimode random lasers,” Science **320**, 643–646 (2008). [CrossRef]

**57. **Y. D. Chong, L. Ge, H. Cao, and A. D. Stone, “Coherent perfect absorbers: time-reversed lasers,” Phys. Rev. Lett. **105**, 053901 (2010). [CrossRef]

**58. **P. C. Waterman, “Matrix formulation of electromagnetic scattering,” Proc. IEEE **53**, 805–812 (1965). [CrossRef]

**59. **D. Maystre, S. Enoch, and G. Tayeb, *Scattering Matrix Method Applied to Photonic Crystals*, K. Yasumoto, ed. (CRC Press, 2010).

**60. **M. G. Moharam, E. B. Grann, and D. A. Pommet, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” J. Opt. Soc. Am. A **12**, 1068–1076 (1995). [CrossRef]

**61. **L. Ge, “Symmetry-protected zero-mode laser with a tunable spatial profile,” Phys. Rev. A **95**, 023812 (2017). [CrossRef]