version of the boundary-integral method described by Heider [@HeiderBIM] in which we directly search for the quasibound states satisfying outgoing-wave boundary conditions at infinity [@NockelMcBook1]. The second approach isolates the high-Q modes by harmonic inversion [@Mandelshtam:2001fk] of temporal field variation collected in a finite-difference time domain simulation with perfectly matched layer boundary condition (MEEP with post-processing by Harminv) [@MEEP-citation]. In the direct quasibound-state calculation, we obtain complex wavenumbers $k_{\text{QB}}=k-i\,\kappa$ to determine Q using Eq. (\[eq:QfactorDef\]).

Quasibound states\[subsec:Quasibound-states\]
---------------------------------------------

In the boundary-integral approach, we leverage the fact that the dielectric defining the billiard is uniform so that a Green-function description of the interior and exterior fields purely in terms of the boundary is possible. At the dielectric interface, the electric field is assumed to satisfy the boundary conditions for TM polarization (electric field perpendicular to the plane). After discretization along the interface, this leads to a nonlinear eigenvalue problem $A(k_{\text{QB}})\,\boldsymbol{u}=\boldsymbol{0}$ where $\boldsymbol{u}$ contains the source values of the electric field and its normal derivative, and $A$ is a matrix obtained from the field matching equations. A non-trivial solution requires searching for $k_{\text{QB}}$ in the complex plane, which we do using a predictor-corrector method [@HeiderBIM]. The method requires discretization of the boundary 