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

Li Ge

doi:10.1364/PRJ.5.000B20

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

Keywords

(080.6755) Systems with special symmetry.(290.5825) Scattering theory (140.3945) Microcavities

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 (ω)$ , 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 $S (ω)$ [i.e., $S (ω) S^{†} (ω) = 1$ ]. When $S (ω)$ 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 $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 $Ψ$ 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 $R$ -matrix method ^[15] in general.

Sign up for Photonics Research TOC. Get the latest issue of Photonics Research delivered right to you！Sign up now

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 [Ψ (x, k), \nabla Ψ (x, k)] = 0$ that holds for all $x$ ’s on the LSS does not exist in general. Instead, we use a nonlocal boundary condition $f [Ψ (x, k), \int_{LSS} d x^{'} O (x, x^{'}) Ψ (x^{'}, 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 $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 $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 $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 ( $PT$ ) ^[27–47] and rotation-time ( $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 $Ψ (x, k)$ is the electric field, $ϵ (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 $Ψ_{L}^{-} (x, k) \equiv \exp [i k (x + L / 2)]$ }, we write the formal solution of $Ψ (x, k)$ as where $Ψ_{L, R}^{+} \equiv \exp [\mp i k (x \pm L / 2)]$ are the outgoing wave functions on the left and right of the system and $r$ , $t$ are the reflection and transmission coefficients. By requiring that both $Ψ (x, k)$ and $\partial_{x} Ψ (x, k)$ are continuous at $x = \mp L / 2$ (denoted by $\dots |_{L, R}$ ), the boundary conditions are then simply by eliminating $r_{L}$ and $t_{L}$ , which are local without involving both $\partial_{x} Ψ |_{L, R}$ or $Ψ |_{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 $2 i k$ , the boundary conditions become the same as those imposed by the CF states ^[18], i.e., with purely outgoing waves. Similarly, an incoming wave $Ψ_{R}^{-} \equiv \exp [- i k (x - L / 2)]$ from the right simply adds an additional constant term $- 2 i k$ 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 $Δ = L / N$ , and the Helmholtz equation takes the following form: where $Ψ_{i}, ϵ_{i} (i = 1,2, \dots, N)$ are the values of the wave function and the dielectric constant at these points. The boundary conditions in Eqs. (3) and (4) can be rewritten as where $q \equiv k Δ$ is dimensionless. The constant term $η \equiv - 4 i q / (2 - i q)$ in Eq. (6) is due to the incoming wave $Ψ_{L}^{-}$ , as mentioned previously, and by dropping it we recover the generalized eigenvalue problem that determines the CF states with purely outgoing boundary condition ^[49]: $ψ_{m}$ is a column vector containing the values $Ψ_{i}$ in the $m$ th CF state, and $ϵ$ is a diagonal matrix with elements $ϵ_{i}$ . The CF frequencies (similar to the resonances or the poles of the $S$ matrix) is given by $k_{m} = q_{m} / Δ$ . Note that the real-valued free-space wave vector $k$ , instead of the complex-valued CF frequencies $k_{m}$ ’s, appears in the $N \times N$ non-Hermitian Hamiltonian $H$ : 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 $H$ would contain the complex-valued resonance frequencies yet to be determined.

Now with the constant term $η$ in the boundary condition in Eq. (6), the wave function in the scattering problem is given by and the inhomogeneous term $F$ is a column vector with a single non-zero element, i.e., $F_{1} = η$ when light is incident from the left. Similarly, for the scattering of a right-incident wave, the only non-zero element of $F$ is $F_{N} = η$ . The equation above can be put into a more explicit form to obtain $Ψ$ :

The transmission and reflection coefficients can then be calculated using for left incidence (and similarly for right incidence), and the $S$ matrix is given by with $t_{L} = t_{R}$ when the system has Lorentz reciprocity ^[50–52].

For comparison, here we also discuss the modal expansion approach to construct $S$ . By inserting $Ψ = \sum_{m} a_{m} ψ_{m}$ to Eq. (10) and utilizing Eq. (8), we find where we have used the “self-orthogonality” of the CF states ^[53] in the following form:

For left incidence and taking $Δ \to 0$ , we find

This expression is identical to that used in the standard derivation of the $R$ -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 $R$ -matrix method resorts to the Green’s theorem instead to take into consideration the boundary condition of $ψ$ , resulting in

We note that the scattered waves in $Ψ$ , as well as the CF states, produce $- i k$ ( $i k$ ) at $x = 0$ ( $L$ ) after taking the derivative, and the corresponding boundary terms above are all canceled. Hence, we find ${[Ψ \partial_{x} ψ_{m} - ψ_{m} \partial_{x} Ψ]}_{0}^{L} = - [Ψ_{L}^{-} (0) \partial_{x} ψ_{m} (0) - ψ_{m} (0) \partial_{x} Ψ_{L}^{-} (0)] = 2 i k ψ_{m} (0)$ , with which we immediately recover Eq. (16). Once $ψ_{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 $PT$ symmetry ^[28–30] and a left incident wave, whose refractive index satisfies $n (- x) = n^{*} (x)$ ^[27,31–48]. Good agreement between $Ψ$ ’s given by Eqs. (11) and (14) are obtained using 50 CF states. Nevertheless, the artifact of $\partial_{x} Ψ$ at the boundary of the microcavity in the $R$ -matrix approach can be readily seen in Fig. 1(b), where we plot the optical flux given by $Im [Ψ^{*} \partial_{x} Ψ]$ (up to a prefactor). The BVP approach, on the other hand, gives a good agreement with the analytical result ^[54] where $G \equiv q s_{1} s_{2}$ , $F \equiv u_{1} s_{1} c_{2} + u_{2} s_{2} c_{1}$ , $D \equiv c_{1} c_{2} - g s_{1} s_{2} - i (h_{1} s_{1} c_{2} + h_{2} s_{2} c_{1})$ , and $q = (n_{1} / n_{2} - n_{2} / n_{1}) / 2$ , $s_{j} = \sin (n_{j} k L / 2)$ , $c_{j} = \cos (n_{j} k L / 2)$ , $u_{j} = (n_{j} - 1 / n_{j}) / 2$ , $g = (n_{1} / n_{2} + n_{2} / n_{1}) / 2$ , $h_{j} = (n_{j} + 1 / n_{j}) / 2 (j = 1,2)$ .

$(a) Total wave function and (b) its flux depicted by black thick lines for a half-gain-half-loss microcavity with light incident from the left. The wave vector k=12/L and refractive indices n1=n2*=2−0.2i are used. The expansion in Eq. (14) with 50 CF states is plotted by the red thin lines as a comparison, which can barely be distinguished from the black line in (a) but shows a significant deviation near the left boundary in (b). The black dot in (b) shows the analytical result at x=−L/2 given by Eq. (18).$

Figure 1.(a) Total wave function and (b) its flux depicted by black thick lines for a half-gain-half-loss microcavity with light incident from the left. The wave vector $k = 12 / L$ and refractive indices $n_{1} = n_{2}^{*} = 2 - 0.2 i$ are used. The expansion in Eq. (14) with 50 CF states is plotted by the red thin lines as a comparison, which can barely be distinguished from the black line in (a) but shows a significant deviation near the left boundary in (b). The black dot in (b) shows the analytical result at $x = - L / 2$ given by Eq. (18).

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 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., $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 $ϵ (x)$ (see Fig. 2), outside of which we assume $ϵ (x) = n_{e}^{2} > 1$ and adopt $Ψ_{m}^{-} = \frac{H_{m}^{-} (n_{e} k r)}{H_{m}^{-} (n_{e} k R)} e^{i m θ}, Ψ_{m}^{+} = \frac{H_{m}^{+} (n_{e} k^{*} r)}{H_{m}^{+} (n_{e} k^{*} R)} e^{i m θ},$ (19)as our incoming and outgoing channels. Here $H_{m}^{\mp}$ are the second and first Hankel functions and $m$ is the angular momentum number, which also serves as the channel index. Note that we do not restrict the free-space wave vector $k$ to be real, which enables the calculation of the resonances as the complex-valued poles of the $S$ matrix.

Figure 2.Schematic of an optical microcavity (shaded area) and the circular LSS (solid line) in 2D. The finite-difference grid is indicated by the dots and dashed lines.

Suppose the incident wave impinges on the LSS in the $m_{0}$ th channel $Ψ_{m_{0}}^{-}$ . The total field and its radial derivative outside the LSS can then be written as $Ψ_{>} = Ψ_{m_{0}}^{-} + \sum_{m} S_{m, m_{0}} Ψ_{m}^{+},$ (20) $\frac{\partial Ψ_{>}}{\partial r} = V_{m_{0}}^{-} (r) e^{i m_{0} θ} + \sum_{m} S_{m, m_{0}} V_{m}^{+ *} (r) e^{i m θ},$ (21)where $V_{m}^{\pm} (r) \equiv n_{e} k H_{m}^{\pm'} (n_{e} k r) / H_{m}^{\pm} (n_{e} k R)$ . Next we discretize the 2D space on a polar grid ^[55], with the circular LSS placed at the middle of the $N_{r}$ th and $(N_{r} + 1)$ th rings (see Fig. 2). On each ring, there are $N_{θ}$ equally spaced grid points with spacing $Δ_{θ} = 2 π / N_{θ}$ . We then write the total field on the $N_{r}$ th and $(N_{r} + 1)$ th rings as $Ψ_{N_{r}, ν} = \sum_{m} b_{m} e^{i m θ_{ν}},$ (22) $Ψ_{N_{r} + 1, ν} = \sum_{m} b_{m} (1 + c_{m} Δ_{r}) e^{i m θ_{ν}},$ (23)where $ν$ is the index for the grid points in the azimuthal direction, $θ_{ν}$ is the corresponding azimuthal angle, and $Δ_{r}$ is the uniform spacing between two consecutive grid points in the radial direction.

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 $Ψ$ and its radial derivative, the left-hand sides of Eqs. (20) and (21) on the LSS can be approximated by $\frac{Ψ_{N_{r} + 1, ν} + Ψ_{N_{r}, ν}}{2} = \sum_{m} b_{m} (1 + \frac{c_{m} Δ_{r}}{2}) e^{i m θ_{ν}},$ (24) $\frac{Ψ_{N_{r} + 1, ν} - Ψ_{N_{r}, ν}}{Δ_{r}} = \sum_{m} b_{m} c_{m} e^{i m θ_{ν}},$ (25)which lead to $S_{m, m_{0}} = \frac{b_{m} c_{m} - V_{m_{0}}^{-} δ_{m, m_{0}}}{V_{m}^{+ *}},$ (26) $S_{m, m_{0}} = b_{m} (1 + c_{m} \frac{Δ_{r}}{2}) - δ_{m, m_{0}},$ (27)and we have dropped the argument $R$ in $V_{m}^{\pm}$ . The product $b_{m} c_{m}$ can be eliminated to derive a more concise form of the $S$ matrix: $S_{m, m_{0}} = \frac{b_{m} - (1 - V_{m_{0}}^{-} \frac{Δ_{r}}{2}) δ_{m, m_{0}}}{1 - V_{m}^{+ *} \frac{Δ_{r}}{2}} .$ (28)

This is the expression we will use in our numerical examples, which requires obtaining $b_{m}$ using the Fourier transform of $Ψ_{N_{r}, ν}$ ’s: $b_{m} = \sum_{ν} \frac{Δ_{θ}}{2 π} e^{- i m θ_{ν}} Ψ_{N_{r}, ν} .$ (29)

To find $Ψ_{N_{r}, ν}$ , 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: $Ψ_{N_{r} + 1, ν} = \sum_{ν^{'}} O_{ν, ν^{'}} Ψ_{N_{r}, ν^{'}} + f_{ν}^{(m_{0})},$ (30) $O_{ν, ν^{'}} \equiv \frac{Δ_{θ}}{2 π} \sum_{m} \frac{H_{m}^{+} (n_{e} k^{*} R^{+})}{H_{m}^{+} (n_{e} k^{*} R^{-})} e^{i m (θ_{ν} - θ_{ν^{'}})},$ (31) $f_{ν}^{(m_{0})} = \frac{(V_{m_{0}}^{-} - V_{m_{0}}^{+ *}) Δ_{r}}{1 - V_{m_{0}}^{+ *} \frac{Δ_{r}}{2}} e^{i m_{0} θ_{ν}},$ (32)where $R^{-}$ is the radius of the $N_{r}$ th ring. Note that the term $f_{ν}^{(m_{0})}$ 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 ^[56], similar to Eq. (6) in 1D.

Equation (30) can then be inserted into the discretized Helmholtz equation on the $N$ th ring to eliminate $Ψ_{N_{r} + 1, ν}$ ^[49], which gives rise to the following matrix equation: $(H + k^{2} ϵ) ψ + F^{(m_{0})} = 0 .$ (33)

The column vector $ψ \equiv \sqrt{r} Ψ$ represents the total wave function, i.e., $ψ = {(ψ_{1,1} \dots ψ_{1, N_{θ}} ψ_{2,1} \dots ψ_{N_{r}, 1} \dots ψ_{N_{r}, N_{θ}})}^{T}$ (34)and $F^{(m_{0})}$ is a column vector of zeros except for the last $N_{θ}$ elements, which are given by $R f_{ν}^{(m_{0})} / (\sqrt{R^{-}} Δ_{r}^{2}) (ν = 1,2, \dots, N_{θ})$ . $ϵ$ has the same form as its 1D counterpart, i.e., a diagonal matrix with the values of the dielectric constant on the discretized grid. $H$ has $(N_{r} \times N_{θ})$ rows and columns; it is the same non-Hermitian Hamiltonian that determines the CF states ^[49,56]: $H ψ_{n} = - k_{n}^{2} ϵ ψ_{n},$ (35)and it consists a banded matrix $H_{0}$ and a $(N_{θ} \times N_{θ})$ block $H^{'}$ in the lower right corner. $H_{0}$ is symmetric with nonzero elements on the 0, $\pm 1$ , and $\pm N_{θ}$ diagonals ${\begin{cases} H_{(μ - 1) N_{θ} + ν, (μ - 1) N_{θ} + ν} & = - \frac{2}{{(Δ_{r})}^{2}} - \frac{2}{{(r_{μ} Δ_{θ})}^{2}} \\ H_{(μ - 1) N_{θ} + ν, (μ - 1) N_{θ} + ν + 1} & = \frac{1}{{(r_{μ} Δ_{θ})}^{2}} \\ H_{(μ - 1) N_{θ} + ν, μ N_{θ} + ν} & = \frac{r_{μ + \frac{1}{2}}}{{(Δ_{r})}^{2} \sqrt{r_{μ} r_{μ + 1}}} \end{cases},$ (36)and $H^{'}$ comes from the nonlocal boundary condition in Eq. (30): $H_{ν, ν^{'}}^{'} = \frac{1}{{(Δ_{r})}^{2}} \frac{R}{R^{-}} O_{ν, ν^{'}} .$ (37) $H^{'}$ can be checked to be also symmetric using $H_{- m}^{+} (z) = {(- 1)}^{m} H_{m}^{+} (z)$ and the definition of $O_{ν, ν^{'}}$ in Eq. (30). Once $ψ$ is obtained for each incoming channel by solving Eq. (33), i.e., $ψ = - {(H + k^{2} ϵ)}^{- 1} F^{(m_{0})},$ (38)as in the 1D case, we immediately know the Fourier coefficients ${b_{m}}$ from Eq. (29), with which the construction of the $S$ matrix is completed using Eq. (28).

In the limit of a fine grid ( $Δ_{r} \to 0$ ), Eq. (28) becomes $S_{m, m_{0}} \overset{Δ_{r} \to 0}{\to} b_{m} - δ_{m, m_{0}} .$ (39)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 $ψ$ in the CF basis, which takes the following form in 2D: $ψ = - \frac{Δ_{r} Δ_{θ}}{π R^{2}} \sum_{n} \frac{ψ_{n} ψ_{n}^{T}}{k^{2} - k_{n}^{2}} F^{(m_{0})} .$ (40)Note that we have used the following normalization of the CF basis $ψ_{n}^{T} ϵ ψ_{n^{'}} = \frac{π R^{2}}{Δ_{r} Δ_{θ}} δ_{n, n^{'}},$ (41)such that $Ψ_{n} \equiv ψ_{n} / \sqrt{r}$ is dimensionless, and its value does not scale with the discretization, i.e., the expression above becomes $\int_{system} ϵ Ψ_{n} Ψ_{n^{'}} r d r d ϕ = π R^{2} δ_{n, n^{'}}$ (42)in the continuous limit.

$ψ$ 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 $Ψ_{m}^{+}$ (which is equivalent to a Fourier transform), i.e., $Ψ_{n} |_{r = R} = \sum_{m} z_{m}^{(n)} Ψ_{m}^{+} |_{r = R} = \sum_{m} z_{m}^{(n)} e^{i m θ},$ (43)the inner product $ψ_{n}^{T} F^{(m_{0})}$ singles out the $m_{0}$ th Fourier coefficient $z_{m_{0}}$ : $ψ_{n}^{T} F^{(m_{0})} = \frac{2 R}{Δ_{r} Δ_{θ}} (V_{m_{0}}^{-} - V_{m_{0}}^{+ *}) z_{m_{0}}^{(n)},$ (44)and, consequently, $S_{m, m_{0}} \overset{Δ_{r} \to 0}{\to} \frac{2}{R} (V_{m_{0}}^{+ *} - V_{m_{0}}^{-}) R_{m, m_{0}} - δ_{m, m_{0}},$ (45)where $R_{m, m_{0}} = \sum_{n} \frac{z_{m}^{(n)} z_{m_{0}}^{(n)}}{k^{2} - k_{n}^{2}}$ (46)is the $R$ matrix ^[15]. This expression clearly indicates 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 $z_{m_{0}}^{(n)}$ and $z_{m}^{(n)}$ , 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 $S_{m, m_{0}}$ given by Eq. (45) is consistent with the standard $R$ -matrix result in the CF basis, we apply the Green’s theorem to the interior of the LSS, which gives us $a_{n} = \frac{1}{π R^{2}} \int_{system} ϵ Ψ_{n} Ψ d x = \frac{1}{π R} \frac{\int_{LSS} {[Ψ \partial_{r} Ψ_{n} - Ψ_{n} \partial_{r} Ψ]}_{R} d θ}{k^{2} - k_{n}^{2}} .$ (47)

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 $ψ_{n}$ satisfy $[\nabla^{2} + n_{e}^{2} k^{2}] X = 0, X = Ψ, Ψ_{n} .$ (48)

Note that the same wave vector $k$ for the total field $Ψ$ and the CF states $ψ_{n}$ is crucial to remove the outgoing waves from the boundary integral in Eq. (47), and we end up with $a_{n} = \frac{1}{π R} \frac{\int_{LSS} {[Ψ_{m_{0}}^{-} \partial_{r} Ψ_{n} - Ψ_{n} \partial_{r} Ψ_{m_{0}}^{-}]}_{R} d θ}{k^{2} - k_{n}^{2}}$ (49)for the incoming wave in the $m_{0}$ th channel. Unlike the $R$ matrix in other bases ^[15], here the $S$ matrix does not appear in the expansion coefficient $a_{n}$ . By realizing that the expansion in Eq. (43) holds not just on the LSS but also in the exterior region, we find $\partial_{r} Ψ_{n} |_{r = R} = \sum_{m} z_{m}^{(n)} V_{m}^{+ *} e^{i m θ},$ (50)and the total wave function in the interior of the system is given by $Ψ (x) = \frac{2}{R} (V_{m_{0}}^{+ *} - V_{m_{0}}^{-}) \sum_{n} \frac{z_{m_{0}}^{(n)} Ψ_{n} (x)}{k^{2} - k_{n}^{2}} .$ (51)

Once we substitute the $Ψ (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: $\frac{H_{m}^{+^{'}} (k R) J_{m}^{+^{'}} (n k R)}{H_{m}^{+} (k R) J_{m}^{'} (n k R)} = n .$ (52)

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 [0,12]$ .

Figure 3.Resonances of a microdisk cavity with a uniform index $n = 1.5$ . The crosses are the analytical results given by Eq. (52), and the circles are the poles of the $S$ matrix constructed using Eq. (28).

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 $PT$ and $RT$ symmetries. It was found that, in a $PT$ -symmetric system, $s_{m}$ ’s undergo spontaneous symmetry breaking as a function of frequency or system size ^[37]: $PT$ symmetry warrants $s_{m}^{*} = s_{m^{'}}^{- 1}$ . When $m, m^{'}$ in this expression are the same, the corresponding eigenstate of the $S$ matrix is in the $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 $PT$ -broken phase $m \neq m^{'}$ , and we find $| s_{m} | = {| s_{m^{'}} |}^{- 1}$ instead (or equivalently, $\log_{10} {| s_{m} |}^{2} = - \log_{10} {| s_{m^{'}} |}^{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 (x) = 1.5 + 0.4 \sin θ$ , which satisfies not only $PT$ symmetry (here $P$ changes $θ$ to $- θ$ ) but also $RT$ symmetry ^[40,48], i.e., $n (r, θ) = n^{*} (r, θ + π)$ . The eigenvalues of the $S$ matrix for an $RT$ -symmetric structure display similar spontaneous breaking as those in their $PT$ -symmetric counterparts, and these properties are nicely manifested by the $S$ matrix given by Eq. (28) using the BVP approach [Fig. 4(a)].

$(a) Spontaneous symmetry breaking of S-matrix eigenvalues sn in a microdisk cavity with PT and RT symmetries. Its refractive index is given by n(x)=1.5+0.4 sin θ, the imaginary part of which is shown schematically by the inset in (b). (b) A real-valued eigenvector of S at kR=4 in the PT- and RT-symmetric phase. The blue (pink) bars show symmetric and antisymmetric components with opposite m’s. The corresponding wave function is shown in (d), where the cavity boundary is marked by the white circle. The wave function of a scattering eigenstate in the broken-symmetry phase at kR=4 is shown in (c) as a comparison.$

Figure 4.(a) Spontaneous symmetry breaking of $S$ -matrix eigenvalues $s_{n}$ in a microdisk cavity with $PT$ and $RT$ symmetries. Its refractive index is given by $n (x) = 1.5 + 0.4 \sin θ$ , the imaginary part of which is shown schematically by the inset in (b). (b) A real-valued eigenvector of $S$ at $k R = 4$ in the $PT$ - and $RT$ -symmetric phase. The blue (pink) bars show symmetric and antisymmetric components with opposite $m$ ’s. The corresponding wave function is shown in (d), where the cavity boundary is marked by the white circle. The wave function of a scattering eigenstate in the broken-symmetry phase at $k R = 4$ is shown in (c) as a comparison.

It can be easily seen that, when a system has both $PT$ and $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 $PT$ and $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 $PT$ -symmetric, i.e., $PT Ψ_{m}^{-} = Ψ_{m}^{+}$ , where the minus sign introduced by the parity operator (again $θ \to - θ$ ) 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 $Ψ^{-} (r, θ) = \sum_{m} v_{m} Ψ_{m}^{-} (r, θ)$ is the incident wave in an eigenstate of $S$ with eigenvalue $s_{n}$ , then by performing a combined $PT$ operation on its scattered wave [i.e., $Ψ^{+} (r, θ) = s_{n} \sum_{m} v_{m} Ψ_{m}^{+} (r, θ)$ ] the new incoming state ${\bar{Ψ}}^{-} (r, θ) \equiv PT Ψ^{+} (r, θ) = s_{n}^{*} \sum_{m} v_{m}^{*} Ψ_{m}^{-} (r, θ)$ should also be an eigenstate of the $S$ matrix due to $PT$ symmetry. It then follows that $v_{m} = g v_{m}^{*}$ should hold for all $m$ ’s in the $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 $PT$ -symmetric phase. If we apply a similar analysis of $RT$ symmetry, we find that the channel functions transform according to $RT Ψ_{m}^{-} = {(-)}^{m} Ψ_{- m}^{+}$ , which leads to $v_{- m}^{*} = \pm {(-)}^{m} v_{m}$ . Therefore, we find $v_{- m} = \pm {(-)}^{m} v_{m} (v_{m} \in R),$ (53)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 $PT$ and $RT$ symmetries also has mirror symmetry about the $θ = \pm π / 2$ axis (i.e., the axis perpendicular to the parity axis in $PT$ symmetry) ^[40], which imposes the following property: $v_{- m} = \pm {(-)}^{m} v_{m} (v_{m} \in C) .$ (54)

The overall $\pm$ sign corresponds to scattering eigenstates that are even and odd functions about the $θ = \pm π / 2$ axis, respectively. Again, the mirror symmetry about the $θ = \pm π / 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 $PT Ψ \neq Ψ$ , $RT Ψ \neq Ψ$ and hence $| Ψ (r, θ) | \neq | Ψ (r, - θ) |, | Ψ (r, θ) | \neq | Ψ (r, θ + π) |$ ; in the latter, we find $| Ψ |$ is symmetric about both the horizontal and vertical axes instead.

The corresponding resonant modes with resonant frequencies $k_{m} R = 3.5042 - 0.2492 i, 4.5989 - 0.6445 i$ 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.

Figure 5.Resonant modes corresponding to the scattering eigenstates in Figs. 4(c) and 4(d).

It is important to note that, while the resonant modes of the system have mirror symmetry about the $θ = \pm π / 2$ axis, here, they do not possess a $PT$ - and $RT$ -symmetric phase: performing a combined $PT$ or $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 $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 $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 $Ψ_{m}^{-}$ to $e^{- i m θ}$ ^[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 $PT$ - and $RT$ -symmetric systems, similar to the finding in 1D heterostructures ^[38].