The method in sequence

This page is the shortest complete account of what the solver actually does. It is meant to be read before any result: every number elsewhere on this site is an output of the chain below, and the chain is what should be evaluated first.

The single most important thing to understand is what the analytical quality factor is not. There is no empirical formula for \(Q\), no fitted regression, and no estimate from a mirror reflectivity. We construct an approximate open-boundary Maxwell operator for the cavity, solve it for a complex quasinormal-mode frequency \(\widetilde\omega\), and read \(Q\) off that pole. Everything else is machinery in service of making the very small imaginary part of \(\widetilde\omega\) meaningful.

1. What \(Q\) means for an open resonator

Use the time convention \(E(t)\propto e^{-i\widetilde\omega t}\). An outgoing quasinormal mode has

\[ \widetilde\omega=\omega_r+i\omega_i,\qquad \omega_i<0, \]

so the field amplitude decays as \(e^{\omega_i t}\) and the energy, being quadratic in the field, as \(e^{2\omega_i t}\). Writing the energy decay as \(U(t)=U_0e^{-\kappa t}\) gives \(\kappa=-2\omega_i\), and with \(Q=\omega_r/\kappa\),

\[ \boxed{\;Q=-\frac{\Re\widetilde\omega}{2\,\Im\widetilde\omega}\;} \]

That conversion is trivial. The entire theoretical problem is the step before it: find the complex resonance frequency of an open nanobeam accurately enough that \(\Im\widetilde\omega\), which is five orders of magnitude below \(\Re\widetilde\omega\), carries information.

2. Choose a background that is already most of the answer

The cavity is a silicon nanobeam with a finite array of etched holes. Rather than treating the whole structure as a perturbation of free space, we make the unperforated rectangular waveguide the background:

\[ \epsilon(\mathbf r)=\epsilon_{\rm wg}(y,z)+\Delta\epsilon_{\rm holes}(\mathbf r), \qquad \Delta\epsilon=\epsilon_{\rm air}-\epsilon_{\rm Si}<0 \ \text{inside a hole}. \]

The complicated object is now a known uniform waveguide plus a finite set of localized perturbations. This is what makes the reduction possible: the unknown only has to live inside the holes.

3. The outgoing integral equation

Source-free frequency-domain Maxwell,

\[ \nabla\times\nabla\times\mathbf E-\frac{\omega^2}{c^2}\epsilon\,\mathbf E=0, \]

split on the background dielectric and inverted with the outgoing Green operator \(G^{+}_{\rm wg}\) of the bare waveguide, becomes the homogeneous Lippmann–Schwinger equation

\[ \mathbf E=\frac{\omega^2}{c^2}\,G^{+}_{\rm wg}(\omega)\,\Delta\epsilon\,\mathbf E . \]

This is already an eigenvalue problem: a resonance occurs where \(K(\omega)=k_0^2G^{+}_{\rm wg}\Delta\epsilon\) has eigenvalue one. In the implementation we locate zeros of the finite Galerkin operator rather than claiming a convergent Fredholm determinant, because the electric Maxwell Green function carries a non-compact longitudinal contribution. That distinction is recorded deliberately; see observable conventions.

4. Solve only for the polarization inside the holes

Define \(\mathbf P=\Delta\epsilon\,\mathbf E\), which vanishes everywhere outside the hole volumes. Then

\[ \boxed{\;\bigl[I-k_0^2\,\Delta\epsilon\,G^{+}_{\rm wg}(\omega)\bigr]\mathbf P=0\;} \]

There is no voxelized air domain and no simulation box. The unknown exists only where the geometry actually differs from the background.

5. The bare-waveguide Green tensor

Because the unperforated waveguide is translationally invariant along the beam axis \(x\), transform along it:

\[ G_{\rm wg}(x,\mathbf r_\perp;x',\mathbf r'_\perp;\omega) =\int_{\mathcal C_\beta}\frac{d\beta}{2\pi}\, e^{i\beta(x-x')}\,\widetilde G_{\rm wg}(\beta;\mathbf r_\perp,\mathbf r'_\perp;\omega). \]

The 3-D waveguide problem becomes a continuum of 2-D transverse problems labelled by longitudinal momentum \(\beta\). The contour \(\mathcal C_\beta\) is not a detail: it is what makes the solution outgoing rather than incoming.

Cladding first. For homogeneous cladding with \(k_c=n_c\omega/c\) and longitudinal/transverse projectors \(P_L=\mathbf k\mathbf k^T/k^2\), \(P_T=I-P_L\),

\[ \widetilde G_c(\mathbf k,\omega)=\frac{P_T}{k^2-k_c^2}-\frac{P_L}{k_c^2}. \]

The second term is the longitudinal/self part. Simplified formulations tend to mishandle it because in real space it is tied to the singular local structure of the electric dyadic; we keep it explicitly in momentum space.

Then the core, nonperturbatively. With \(\chi_b=\epsilon_{\rm Si}-\epsilon_{\rm air}\), the cross-sectional Dyson equation \(G_{\rm wg}=G_c+G_ck_0^2\chi_bG_{\rm wg}\) is solved, not iterated:

\[ G_{\rm wg}=\bigl[I-k_0^2\chi_bG_c\bigr]^{-1}G_c . \]

Silicon is never treated as a small perturbation. Guided modes then emerge as poles of \(A_{\rm wg}(\beta,\omega)=I-k_0^2\chi_bG_c\), i.e. from \(\det A_{\rm wg}(\beta_g,\omega)=0\), with the residue at a simple pole fixing the guided-state normalization. Guided and radiating fields are not two pieces of physics glued together afterwards; they are different analytic structures of one waveguide resolvent.

The transverse problem is represented in normalized Legendre functions whose Fourier transforms are known in closed form — over beam width \(W\), \(\phi^{(W)}_m(y)=\sqrt{(2m+1)/W}\,P_m(2y/W)\) with \(\Phi^{(W)}_m(q_y)=\sqrt{W(2m+1)}\,(-i)^m j_m(q_yW/2)\) — so the matrix elements are explicit Bessel/rational integrals rather than a transverse grid.

6. Where the imaginary part physically comes from

This is the part worth reading twice. The transverse denominator contains \(\beta^2+q_y^2+q_z^2-k_c^2\), so at fixed \(\beta\) with \(\kappa^2=k_c^2-\beta^2\) the radial integral meets a real singularity on the light shell. Under the outgoing prescription, Sokhotski–Plemelj gives

\[ \int_0^Q\frac{q\,H(q)}{q^2-\kappa^2-i0}\,dq =\int_0^Q\frac{q\bigl[H(q)-H(\kappa)\bigr]}{q^2-\kappa^2}\,dq +\frac{H(\kappa)}{2}\ln\frac{Q^2-\kappa^2}{\kappa^2} +\boxed{\frac{i\pi}{2}H(\kappa)} . \]

The boxed term is the exact on-shell radiation contribution. We never assign a phenomenological linewidth to radiation modes or an arbitrary imaginary part to free-space channels. The continuum boundary condition itself generates the imaginary part of the Green function — the classical-wave analogue of a discrete state acquiring a width by coupling to a continuum.

7. Fields inside each hole

For hole \(j\), map the ellipse with \(x-x_j=r_{xj}\rho\cos\theta\), \(y=r_{yj}\rho\sin\theta\), \(0\le\rho\le1\), and expand the polarization in real Zernike functions across the ellipse, Legendre polynomials through the thickness, and Cartesian polarization:

\[ \boldsymbol h_{j,nm\sigma\ell s} =Z^{\sigma}_{nm}(\rho,\theta)\,P_\ell(2z/T)\,\mathbf e_s . \]

A hole is a handful of smooth functions shaped like a hole, not thousands of cubic cells. In the production Quan model the retained hole basis is P2-z1 against a y3-z1 waveguide basis.

8. Assemble the hole-to-hole propagator

Projecting each hole basis onto the waveguide Green function at fixed \(\beta\) gives \(C(\beta)\), where the thickness integral is exactly orthogonal by construction and the elliptical-disk integral is smooth and quadrature-friendly. Concatenating all holes,

\[ G^{+}_H(\omega)=\int_{\mathcal C_\beta}\frac{d\beta}{2\pi}\, C(-\beta)^{T}\,G_{\rm wg}(\beta,\omega)\,C(\beta). \]

Note \(C(-\beta)^T\) and not \(C^\dagger\). At complex frequency the operator must remain an analytic function of \(\omega\); conjugation would destroy that. Reciprocity gives a complex-symmetric structure instead. This matrix carries hole \(i\) → all guided and radiation channels → hole \(j\), for every pair at once, which is what allows radiation from different parts of the cavity to interfere. High \(Q\) is that interference; omitting even a weak radiation subspace can leave the field profile right and the linewidth meaningless.

9. The cavity as a nonlinear eigenproblem

With the geometry vector \(\mathbf p\) (longitudinal spacings and both hole radii — 45 free coordinates for the Quan device),

\[ \boxed{\;A(\omega,\mathbf p)=I-k_0^2\,\Delta\epsilon\,G^{+}_H(\omega,\mathbf p), \qquad A(\widetilde\omega,\mathbf p)\,x=0.\;} \]

We search complex \(\omega\) until the operator acquires a null vector. Because many poles exist, the branch is tracked by right-null-vector overlap with the previous mode, together with the target symmetry sector and polarization content, so the optimizer cannot silently hop onto a neighbouring resonance with a more flattering \(Q\).

Then, and only then, \(Q=-\omega_r/(2\omega_i)\).

10. The gradient is what makes it an optimizer

With left and right null vectors \(y,x\), differentiating \(A(\widetilde\omega,\mathbf p)x=0\) and left-multiplying by \(y^\dagger\) annihilates the unknown eigenvector derivative:

\[ \boxed{\;\frac{d\widetilde\omega}{dp_j} =-\frac{y^\dagger (\partial_{p_j}A)\,x}{y^\dagger (\partial_\omega A)\,x}\;} \]

Nothing is differentiated through Newton's method, through an eigensolver history, or through a time-domain ringdown. The denominator is also the QNM residue normalization. The full contract, including the boundary (Hadamard) form for arbitrary interfaces, is on the QNM boundary gradient page. Converting to the observable,

\[ \frac{dQ}{dp}=-\frac{1}{2\omega_i}\frac{d\omega_r}{dp} +\frac{\omega_r}{2\omega_i^{2}}\frac{d\omega_i}{dp}, \]

which is also the reason high-\(Q\) optimization is delicate: the second coefficient scales as \(\omega_i^{-2}\), so a model can have the resonance wavelength right and a field overlap above 0.999 and still produce a useless \(dQ/dp\).

11. One step of the loop

Each accepted step is: evaluate the pole and gradient; project the proposed move into the frequency-neutral subspace,

\[ \Re\bigl(\nabla_{\mathbf p}\widetilde\omega\bigr)\cdot\Delta\mathbf p=0, \]

so that improvement comes from reducing loss rather than from sliding the resonance; apply fabrication constraints (positive radii, minimum ligaments, no overlap, fixed topology, symmetry); take a small trust-region step; then re-solve the pole from scratch on the tracked branch. The gradient proposes and the nonlinear operator disposes.

Q = -Re w / 2 Im w                 what we want
        ^
outgoing pole   A(w,p) x = 0       nonlinear eigenproblem
        ^
hole-to-hole propagator  G_H+      coherent, all pairs at once
        ^
Zernike x Legendre hole basis      unknown lives only in the holes
        ^
guided poles + radiation continuum + evanescent
        ^
G_wg = (I - k0^2 chi G_c)^-1 G_c   Dyson-resummed silicon core
        ^
G_c = P_T/(k^2-kc^2) - P_L/kc^2    exact cladding Green tensor

then: implicit dw/dp  ->  dQ/dp  ->  project out Re dw
      ->  constrained trust-region step  ->  re-solve the pole  ->  repeat

12. What this establishes, and what it does not

The construction is internally exact in the sense that both the real and the imaginary part of \(\widetilde\omega\) come from solving Maxwell's pole equation in a retained subspace. It is not a claim that the retained subspace reproduces the physical linewidth. It is a reduced Galerkin approximation, and its absolute \(Q\) carries a systematic bias that is itself geometry-dependent:

• Quan–Lončar, matched endpoints. Analytical \(7455.16\to18037.3\) against independently fitted FDTD \(15449.47\to59337.16\), a \(3.84\times\) FDTD gain. A further eight-step continuation of the same ascent reached analytical \(30307.26\) but was never evaluated in FDTD; it is retained under a separate name and supports nothing.
• Design23, matched endpoints. Analytical \(8958.74\to717003\) — an \(80\times\) analytical gain — against FDTD \(18991.94\to58240.27\), a \(3.07\times\) FDTD gain. The surrogate bias drifts from \(0.472\) to \(12.31\) along the very trajectory it is driving.
• What transfers. Rank, not magnitude: Spearman \(0.911\) within the fixed 12 ps protocol segment, \(0.774\) over all qualified observations. Per-step sign agreement is explicitly not claimed.
• How it fails. The failure mode is a direction the retained operator cannot see. The pre-repair operator once scored a genuinely flat direction as a \(3.6\times\) gain, with near-unity field overlap throughout. The model can drive its retained radiation operator toward a null while the device keeps radiating into channels the basis omits.

Both device gains remain provisional pending the matched convergence panel; the current promotion gates and the ringdown-adequacy conditions are on the status and claims page. The honest summary of the method is that it is a fast, structured, differentiable surrogate for the physical loss functional — useful for its gradient, not for its absolute linewidth.