Mathematical derivations and approximation ledger

This document derives the equations implemented by the repository. It is written so that a new investigator can reconstruct the method without the original conversation. Symbols follow the code wherever possible.

1. Conventions and normalized units

The implementation uses

\[ \lambda_0=c=\mu_0=1, \qquad k_0=\frac{\omega}{c}=2\pi f, \]

where frequency is the dimensionless frequency \(f\). A physical vacuum wavelength of 780 nm therefore corresponds to \(f=1\). Lengths in the input configuration are measured in \(\lambda_0\), and a dimensionless complex pole \(f_p\) maps to

\[ \nu_p=\frac{c_\mathrm{SI}}{780\,\mathrm{nm}}f_p, \qquad \omega_p=2\pi\nu_p. \]

We use the \(e^{-i\omega t}\) convention. Passive loss and the retarded-sheet selector have positive imaginary permittivity; decaying temporal poles lie in the lower half frequency plane. For

\[ f_p=f_r+i f_i,\qquad f_i<0, \]

the field-amplitude decay rate is \(-2\pi f_i\) in normalized angular units and

\[ Q=-\frac{f_r}{2f_i}. \]

Materials are currently linear, isotropic, nondispersive, and nonmagnetic in the physical region. The design interpolation is

\[ \epsilon(\mathbf r)=\epsilon_\mathrm{clad} +\rho(\mathbf r)\Delta\epsilon, \qquad \Delta\epsilon=n_\mathrm{core}^2-n_\mathrm{clad}^2. \]

The default indices are 2.0 and 1.0.

2. Continuous frequency-domain Maxwell problem

Eliminating \(\mathbf H\) from Maxwell's equations gives the driven electric field problem

\[ \mathcal A(\omega,\rho)\mathbf E =\left[\nabla\times\mu^{-1}\nabla\times -\omega^2\epsilon(\rho)\right]\mathbf E =\mathbf d. \]

The source vector \(\mathbf d\) represents an even, z-oriented, centered electric dipole. Overall source normalization does not affect pole location or Q; it does affect the numerical scale of the response.

The semi-infinite mirrors and transverse radiation condition make the effective operator frequency dependent and non-Hermitian. Reciprocity is retained, so after consistent discretization

\[ A^T=A, \]

but generally \(A^\dagger\ne A\). That distinction is why the adjoint bilinear form uses transpose, not complex conjugation.

3. Yee discretization and discrete de Rham structure

3.1 Unknown locations

Cell-centered material samples are stored on an \(n_x\times n_y\times n_z\) Cartesian grid. The three electric components live on staggered Yee edges. The discrete curl maps these electric edges to the corresponding dual faces.

Let \(D_x,D_y,D_z\) be open-boundary forward-difference matrices. With Kronecker products implied, the discrete curl is

\[ C= \begin{bmatrix} 0&-D_z&D_y\\ D_z&0&-D_x\\ -D_y&D_x&0 \end{bmatrix}, \]

and the edge gradient is

\[ G=\begin{bmatrix}D_x\\D_y\\D_z\end{bmatrix}. \]

Because the tensor-product differences commute,

\[ CG=0 \]

to machine precision (exactly zero in the current sparse construction). The solver therefore does not need the grad-div penalty used in an earlier collocated prototype.

3.2 Cell-to-edge material interpolation

Let \(W_x,W_y,W_z\) average cell-centered permittivity to the appropriate staggered electric edges, including cladding completion at exterior half cells. Stack these maps as

\[ W=\begin{bmatrix}W_x\\W_y\\W_z\end{bmatrix}. \]

For cell permittivity vector \(\boldsymbol\epsilon_c\), the edge samples are

\[ \boldsymbol\epsilon_e=W\boldsymbol\epsilon_c. \]

The topology derivative must return through \(W^T\), not through an ad hoc cell average. This exact transpose is tested by directional finite differences.

3.3 Discrete operator

Let \(T_e\) be the diagonal electric transformation-optics PML tensor on Yee edges and \(T_{\mu^{-1}}\) the magnetic-inverse tensor on curl faces. The assembled operator is

\[ A_\mathrm{bare}(f,\rho) =C^T\operatorname{diag}(T_{\mu^{-1}})C -(2\pi f)^2 \operatorname{diag}\!\left(T_e\odot W\boldsymbol\epsilon_c(\rho)\right). \]

The same difference operators occur on both sides of the magnetic term, and the material/PML mass is diagonal, so the bare matrix is complex symmetric.

4. Transformation-optics PML

No x-directed PML is used. The semi-infinite periodic mirrors supply the longitudinal radiation/band-gap condition. In y and z the coordinate stretches are

\[ s_y(y)=1+i\sigma_y(y),\qquad s_z(z)=1+i\sigma_z(z), \]

with

\[ \sigma_u(u)=S \left[\frac{\max(|u|-u_\mathrm{start},0)} {u_\mathrm{outer}-u_\mathrm{start}}\right]^m. \]

There is no x stretch, so \(s_x=1\). Transformation optics gives the electric tensor factors

\[ T_e=\operatorname{diag} \left(s_y s_z,\frac{s_z}{s_y},\frac{s_y}{s_z}\right) \]

and the inverse-magnetic factors

\[ T_{\mu^{-1}}=\operatorname{diag} \left(\frac{1}{s_y s_z},\frac{s_y}{s_z},\frac{s_z}{s_y}\right). \]

They are sampled on their correct edge/face sublattices. The dual samples are shifted a half cell outward, with reflection-safe absolute-coordinate logic.

The design envelope ends before pml_start_y and pml_start_z. That ordinary cladding buffer is essential: allowing design voxels into an absorber lets an optimizer create PML-dependent numerical resonances.

5. Topology parameterization and immutable atom channel

Let \(p_j\in[0,1]\) be raw design variables on permitted cell centers. A normalized compact-support density filter is

\[ \widetilde p_i=\sum_j F_{ij}p_j, \qquad F_{ij}=\frac{\max(R-|\mathbf r_i-\mathbf r_j|,0)} {\sum_k\max(R-|\mathbf r_i-\mathbf r_k|,0)}. \]

The smoothed Heaviside projection is

\[ \rho_i= \frac{\tanh(\beta\eta)+\tanh[\beta(\widetilde p_i-\eta)]} {\tanh(\beta\eta)+\tanh[\beta(1-\eta)]}, \]

with derivative

\[ \frac{\partial\rho_i}{\partial\widetilde p_i} =\frac{\beta\,\operatorname{sech}^2[\beta(\widetilde p_i-\eta)]} {\tanh(\beta\eta)+\tanh[\beta(1-\eta)]}. \]

After projection, terminal mirror-buffer cells are set to one and every cell with \(|x|<0.5\lambda_0\) is set to zero. Neither region is a design variable. For a full-density gradient \(g_\rho\), the raw-variable gradient is

\[ g_p=F^T\left(g_\rho\odot \frac{\partial\rho}{\partial\widetilde p}\right). \]

Optional reflection projection averages all eight x/y/z reflected copies of a design vector or gradient. This constrains topology but does not reduce the electromagnetic solve to a parity subspace in the legacy open-transverse backend.

5.1 Exact y/z electromagnetic parity sector

The opt-in reflection-closed backend uses periodic y/z Yee differences at the outer edge of the transverse PML. Cell centers reflect as \(i\mapsto n-1-i\), while lower staggered edges reflect as \(i\mapsto-i\pmod n\). These closed index sets make the reflection sector an exact invariant subspace of the discrete operator. For the centered z-directed source, orthonormal component bases select y parities \((+,-,+)\) and z-sector parities \((-,-,+)\) for \((E_x,E_y,E_z)\). Writing their direct sum as \(Q\), the reduced solve is

\[ (Q^T A Q)e_r=Q^T d,\qquad e=Qe_r. \]

This is a transpose projection, consistent with complex symmetry. On the fine profile it reduces 68,640 field unknowns to 16,960. It also projects the periodic lead cell and surface plane before dense condensation/decimation. The transverse closure is behind the PML and therefore still requires the normal absorber-convergence study. It is a boundary-discretization change, so topology may be transferred but poles and Q may not. Exactness evidence, benchmarks, and the unresolved x-reflection issue are recorded in SYMMETRY_ACCELERATION.md.

6. Turning an infinite periodic mirror into a boundary condition

This section is the central 2.5D/3D construction. The mirror is periodic only outside the finite defect: for example, the right half-space begins at the right defect plane and repeats for all larger x. That is exactly the setting in which a surface Green function is a boundary condition. Global periodicity is neither assumed nor desired.

6.1 Block-tridiagonal periodic half-space

Choose one geometric mirror period and collect its electromagnetic unknowns into \(q_n\), where \(n=0,1,2,\ldots\) counts cells away from the defect. Yee locality makes the half-space block tridiagonal:

\[ V^Tq_{n-1}+H(f)q_n+Vq_{n+1}=0. \]

The period can contain arbitrary etched cross-sectional detail. Complexity changes the entries and size of \(H,V\); it does not change the reduction. The only structural requirements are eventual periodicity and finite-range coupling across the chosen period boundary.

The left and right leads are constructed independently because their surface orientation and mirror phase are reversed. The period adjacent to the defect must have exactly the same Yee-averaged material phase as an explicitly meshed continuation.

6.2 Exact condensation of internal slices

Iterating a dense Green function over every slice in a period is wasteful. Partition one period into a retained transverse surface plane \(s_n\) and internal slices \(i_n\). In the code's orientation, define within-period blocks \(H_{00},H_{0i},H_{i0},H_{ii}\), plus outward cross-period blocks \(F_{i0}\) and \(B_{0i}\). Eliminating the internal variables with \(H_{ii}^{-1}\) gives

\[ \begin{aligned} D_s &=H_{00}-H_{0i}H_{ii}^{-1}H_{i0},\\ D &=D_s-B_{0i}H_{ii}^{-1}F_{i0},\\ \alpha&=-H_{0i}H_{ii}^{-1}F_{i0},\\ \beta &=-B_{0i}H_{ii}^{-1}H_{i0}. \end{aligned} \]

Here \(D_s\) is the exposed terminal-plane block, \(D\) is the repeated bulk plane block, and \(\alpha,\beta\) couple adjacent condensed planes. Reciprocal assembly gives \(\beta=\alpha^T\) to numerical precision. This is an exact Schur complement, not a modal truncation: all internal period degrees of freedom remain represented through frequency-dependent dense blocks.

For the fine profile this reduces one six-slice lead cell from 10,296 unknowns to a 1,716-unknown surface plane before recursion.

6.3 Retarded surface Green function by cyclic reduction

Let \(g_s\) be the inverse seen at the exposed plane of the semi-infinite chain. Formally it obeys a Dyson fixed point. With the implementation's orientation,

\[ g_s^{-1}=D_s-\alpha g_c\beta, \]

where \(g_c\) is the continuation half-space Green function. Direct fixed point iteration can converge slowly near a band edge. Lopez-Sancho cyclic reduction eliminates every other plane at each iteration, so the effective couplings span twice the distance each time.

Starting with surface block \(S=D_s\), bulk block \(D\), and couplings \(\alpha,\beta\), the code applies

\[ \begin{aligned} S'&=S-\alpha D^{-1}\beta,\\ D'&=D-\alpha D^{-1}\beta-\beta D^{-1}\alpha,\\ \alpha'&=\alpha D^{-1}\alpha,\\ \beta'&=\beta D^{-1}\beta. \end{aligned} \]

The process stops when

\[ \frac{\|\alpha\|+\|\beta\|}{\|D\|}<\text{tolerance}, \]

and then \(g_s=S^{-1}\). A vanishing positive lead loss selects the retarded real-axis sheet. A computed Dyson residual checks the returned inverse.

6.4 Schur-complement attachment to the finite defect

Write the finite-defect and semi-infinite-lead equations as

\[ \begin{bmatrix} A_D&B\\ B^T&A_L \end{bmatrix} \begin{bmatrix}e\\q\end{bmatrix} =\begin{bmatrix}d\\0\end{bmatrix}. \]

The lead equation gives

\[ q=-A_L^{-1}B^T e. \]

Substitution produces the exact defect-only system

\[ \left(A_D-BA_L^{-1}B^T\right)e=d. \]

Only the lead surface is coupled to the defect, so the needed block of \(A_L^{-1}\) is \(g_s\). Therefore

\[ \Sigma(f)=B g_s(f)B^T, \]

and with independent left/right mirrors

\[ A_\mathrm{eff}(f,\rho) =A_D(f,\rho)-\Sigma_L(f)-\Sigma_R(f). \]

The self-energy is dense only within its supported boundary rows. It must not be expanded into a dense matrix over the entire defect. Doing that in an early implementation caused a 33 GiB validation run even though the physical surface plane had only 420 unknowns.

6.5 Why boundary translation is the decisive test

If \(\Sigma\) exactly represents the periodic half-space, inserting any integer number of explicit mirror periods between the defect and \(\Sigma\) cannot change the field inside the original defect. The repository compares the direct boundary field with a domain extended by two periods:

\[ \varepsilon_\mathrm{trans} =\frac{\|e_\mathrm{extended}|_D-e_\mathrm{direct}\|} {\|e_\mathrm{direct}\|}. \]

This simultaneously tests period phase, left/right orientation, Yee material averaging at the interface, condensation, surface recursion, and self-energy embedding. It is stronger than checking the surface Dyson equation alone.

7. Driven response, LDOS, and reciprocal adjoint

Solve

\[ A(f,p)e=d \]

at a real frequency. The complex colocated dipole response is

\[ h(f,p)=d^T e=d^T A^{-1}d. \]

Up to fixed physical constants and source normalization, the projected LDOS is

\[ \mathrm{LDOS}=\Im h. \]

Differentiate the state equation:

\[ A,de+(dA)e=0, \qquad de=-A^{-1}(dA)e. \]

Then

\[ dh=-d^T A^{-1}(dA)e. \]

Because \(A^T=A\) and \(Ae=d\),

\[ d^T A^{-1}=e^T, \qquad dh=-e^T(dA)e. \]

Thus a colocated reciprocal dipole needs no additional adjoint solve.

For one cell-density perturbation, the material part of the operator gives

\[ dA=-(2\pi f)^2\Delta\epsilon\, \operatorname{diag}(T_e)W\,d\rho. \]

Therefore

\[ \frac{\partial h}{\partial\rho} =(2\pi f)^2\Delta\epsilon\, W^T\left(T_e\odot e^2\right). \]

There is no complex conjugate on \(e^2\). The fixed-frequency log-LDOS gradient is

\[ \frac{\partial\log(\Im h)}{\partial\rho} =\frac{(2\pi f)^2\Delta\epsilon} {\Im h} \Im\left[W^T(T_e\odot e^2)\right], \]

followed by projection/filter backpropagation.

LDOS is effective for nucleating a source-coupled resonance from scratch, but a single fixed frequency eventually rewards frequency alignment as well as linewidth. The later stages therefore recenter the response or optimize a pole-Q surrogate while retaining independent real-axis validation.

8. Resonance poles and real-axis local continuation

8.1 Nonlinear resonance equation

Quasinormal modes satisfy

\[ A_\mathrm{eff}(f_p,p)v=0. \]

This is nonlinear in frequency because both the Maxwell mass term and mirror self-energies depend on \(f\).

Evaluating the surface recursion directly at \(\Im f<0\) is dangerous. An outgoing quasinormal field grows spatially in the semi-infinite lead, and a naive decimation can move from the retarded to the advanced sheet. The local continuation therefore samples only the retarded real axis.

8.2 First-order pole estimate

Around a real center \(f_c\),

\[ A(f_c+\delta)\approx A_0+\delta A_1, \qquad A_1=\left.\frac{dA_\mathrm{eff}}{df}\right|_{f_c}. \]

The derivative is evaluated by real-axis differences and includes \(-d\Sigma_L/df-d\Sigma_R/df\). Solve

\[ A_0^{-1}A_1v=\mu v. \]

Since \((A_0+\delta A_1)v=0\),

\[ \delta=-\frac{1}{\mu}, \qquad f_p=f_c+\delta. \]

The real center is iterated toward \(\Re f_p\), while vector overlap with the previous checkpoint selects the same source-coupled mode. Freezing the mirror self-energy derivative is invalid: on one development case it changed an apparent Q=556 pole to Q=3.7.

8.3 General implicit pole derivative

Let \(w^TA=0\) be the left null vector. Differentiating \(Av=0\) with respect to a real design parameter \(p_j\) and left-multiplying by \(w^T\) gives

\[ \frac{df_p}{dp_j} =-\frac{w^T A_{p_j}v}{w^T A_f v}. \]

For a simple pole of a complex-symmetric matrix, left and right vectors can be chosen with \(w=v\) under the transpose bilinear form.

For the linearized problem \(A(f_c)+\delta A_f(f_c)\), its topology derivative carries

\[ A_p(f_c)+\delta A_{fp}(f_c). \]

With nondispersive topology material this yields the code coefficient

\[ \frac{4\pi^2\Delta\epsilon \left(f_c^2+2f_c\delta\right)}{v^T A_fv} \]

multiplying \(W^T(T_e\odot v^2)\).

The locally predicted pole is a function \(g(f_c,p)\), and continuation imposes the fixed point \(f_c=\Re g\). A second real-axis derivative supplies

\[ g_c\approx-\delta\frac{v^TA_{ff}v}{v^TA_fv}. \]

For the fixed-point motion,

\[ \frac{df_c}{dp} =\frac{\Re(g_p)}{1-\Re(g_c)}, \qquad \frac{df_p}{dp}=g_p+g_c\frac{df_c}{dp}. \]

Finally,

\[ \log Q=\log(\Re f_p)-\log[-2\Im f_p], \]

\[ \frac{d\log Q}{dp} =\frac{\Re(df_p/dp)}{\Re f_p} -\frac{\Im(df_p/dp)}{\Im f_p}. \]

This direct pole gradient is useful locally but inherits the first-order frequency approximation. Full real-axis rational continuation is the acceptance authority.

9. Real-axis rational continuation and linewidth validation

9.1 Meromorphic response

Near isolated resonances, the retarded response is meromorphic:

\[ h(f)=b(f)+\sum_j\frac{r_j}{f-p_j}, \]

where \(b(f)\) is analytic background, \(p_j\) are lower-half-plane poles, and \(r_j\) are complex residues. Sampling \(h\) on the real axis avoids direct complex-frequency lead-sheet ambiguity.

9.2 AAA representation

AAA selects support samples \((z_k,h_k)\) and forms a barycentric rational approximant

\[ r(z)= \frac{\displaystyle\sum_{k=1}^m\frac{w_k h_k}{z-z_k}} {\displaystyle\sum_{k=1}^m\frac{w_k}{z-z_k}}. \]

At each greedy iteration the weights minimize residuals at nonsupport samples. Poles are zeros of the barycentric denominator. Candidate filters require a lower-half-plane pole within the window, sufficient decay relative to sample spacing, and a non-negligible residue. Interleaved frequency folds test both held-out complex-response error and pole stability.

AAA is not differentiated because support-point selection is discrete and can jump under an arbitrarily small topology perturbation.

9.3 Local Fano-linewidth cross-check

A complex residue and analytic background produce an asymmetric Fano-like real-axis LDOS rather than a pure Lorentzian. The local validator therefore fits a smooth background plus the real/imaginary contribution of a complex simple pole. Its linewidth gives an independent Q estimate. Agreement with AAA is necessary but not sufficient; the rejected PML mode also had apparent pole/linewidth agreement.

9.4 Periodic-lead branch points and the sheet-local protocol

The semi-infinite periodic lead is not globally meromorphic in frequency. Opening or closing a propagating Bloch channel creates a square-root branch point in its surface Green function. A model of the form

\[ h(f)=b(f)+\sum_j\frac{r_j}{f-p_j} \]

requires \(b(f)\) to be analytic on the fitted interval and therefore must not be forced across a lead band edge. The wide seven-probe sweep still covers the complete local cluster and both edge neighborhoods. For pole continuation, sheet-local-v2 selects, before evaluating confirmatory data, the dense local samples within plus/minus two halfwidths of the explicit tracking pole. A geometry gate requires at least 1.75 halfwidths of data on both sides. The rational order ensemble, Fano fit, and SIMO Loewner check use that one-sheet interval. A separately solved staggered grid drawn only from the same interval remains the external response-prediction holdout.

On the exact y/z-reduced fine grid, the legacy merged fit crossed sharp lead features near 0.99844 and 1.01451. Its confirmatory NRMSE was 1.49e-4 and its linewidth disagreement was 10.87%. With the full cluster retained for identity but the pole fit restricted to the predeclared sheet-local interval, 209 acquisition samples give interleaved NRMSE 2.68e-10; an untouched 28-point confirmation gives 4.35e-12, and linewidth disagreement is 1.36%. No response threshold was loosened.

10. Differentiable Loewner pole surrogate

Given scalar response samples \((f_i,h_i)\), split alternating samples into left data \((\mu_i,v_i)\) and right data \((\lambda_j,w_j)\). The Loewner and shifted-Loewner matrices are

\[ L_{ij}=\frac{v_i-w_j}{\mu_i-\lambda_j}, \qquad S_{ij}=\frac{\mu_i v_i-\lambda_j w_j}{\mu_i-\lambda_j}. \]

SVDs of \([L\;S]\) and \([L;S]\) give left/right compression bases \(Y,X\). The reduced pencil is

\[ S_r x=f_p L_r x, \qquad L_r=Y^HLX,quad S_r=Y^HSX. \]

During one local topology step, \(X,Y\) are frozen. For left generalized eigenvector \(y\),

\[ df_p= \frac{y^H(dS_r-f_p,dL_r)x}{y^H L_r x}. \]

Every entry of \(L,S\) is linear in the response samples, so this produces \(\partial f_p/\partial h_i\). Chain it with the reciprocal Maxwell response gradient from Section 7:

\[ \frac{df_p}{dp_j} =\sum_i\frac{\partial f_p}{\partial h_i} \frac{\partial h_i}{\partial p_j}. \]

The log-Q formula then gives a real topology gradient. Because compression subspaces and pole identities can change over a finite step, the Loewner model is rebuilt after every accepted update and a fresh AAA validator decides acceptance.

11. Overlapping modes: shared-pole multi-output model

11.1 Why one scalar trace is inadequate

When poles overlap, analytic background and neighboring residues can trade off against rational model order. Different scalar observables see the same physical poles but different residues. Requiring common poles across several outputs adds the missing physical constraint.

Let \(h_m(f)\) be output \(m\), with seven outputs in the current solver. The shared-pole model is

\[ h_m(f)=\sum_{k=0}^{d}c_{mk}x(f)^k +\sum_{j=1}^{J}\frac{r_{mj}}{f-p_j}, \]

where \(x\) is a normalized frequency, \(p_j\) are common to every output, and \(c_{mk},r_{mj}\) are output specific.

11.2 Variable projection

For trial poles \(p_j\), form a basis

\[ \Phi(f;p)= \left[1,x,\ldots,x^d, \frac{1}{f-p_1},\ldots,\frac{1}{f-p_J}\right]. \]

The output coefficients are linear and are eliminated exactly by complex least squares:

\[ C(p)=\arg\min_C\|\Phi(p)C-H\|_F^2. \]

Only pole real parts and positive decay rates \(p_j=a_j-i\gamma_j\), \(\gamma_j>0\), remain in the nonlinear bounded least squares. Interleaved holdouts refit without each test subset and predict the withheld complex responses.

11.3 SIMO Loewner check

For a single-input/multiple-output system, build one scalar Loewner block for each output and stack the blocks vertically. A common reduced pencil then estimates shared poles without the nonlinear variable-projection parameterization. Rank perturbations test model-order dependence.

11.4 Authoritative combined gate

The validator uses the wide dataset to reacquire the complete source-coupled cluster and compares its tracked pole with the dense local result. Under the production sheet-local-v2 protocol, response prediction and continuation use the predeclared one-sheet target subset described in Section 9.4. Each of seven probes is fit by high-order AAA with four interleaved holdouts at two independently capped model orders, 29 and 32. The accepted pole is the component-wise median of all 56 order/probe/holdout poles, not a maximum-order full fit. The response and pole gates use the worst result across both orders. The legacy merged-window protocol remains available only to reproduce older artifacts.

The gate requires:

\[ \begin{array}{ll} \text{independent held-out model-order caps} &=\{29,32\},\\ \text{worst holdout response NRMSE} &\le 10^{-5},\\ \text{worst holdout pole displacement} &\le0.15\ |\Im p|,\\ \text{relative holdout Q range} &\le0.15,\\ \text{wide/local Q difference} &\le0.05,\\ \text{wide/local pole displacement} &\le0.15\ |\Im p|,\\ \text{linewidth Q difference} &\le0.10,\\ \text{near-consensus SIMO Q difference} &\le0.10, \end{array} \]

with at least four near-consensus SIMO ranks. These are empirical numerical acceptance tolerances, not universal physical constants.

12. PML, grid, and mode-identity gates

A reliable resonance must be invariant under numerical choices that do not change the physical open system.

12.1 Transverse PML gate

For several absorber starts, strengths, and thicknesses:

recompute both wide and local multi-probe datasets;
identify the same pole using the combined validator;
reject any case with an upper-half-plane pole, failed holdouts, or mode switch; and
bound pole and Q changes relative to the baseline.

PML invariance is evaluated on the shared pole, not on the maximum LDOS peak.

12.2 Grid/domain gate

Transfer projected geometry to a controlled new grid, but reacquire the spectrum from a broad window. Compare dimensionless pole, Q, localization, and later port/atom observables. Refine one controlled aspect at a time when possible so a large shift can be attributed to x spacing, transverse spacing, slab resolution, defect length, or outer-domain size.

12.3 Mode tracking

At local continuation stages, use a bilinear/vector-overlap score with the previous pole mode. At rational stages, require agreement across frequency windows, probes, holdouts, linewidth, and Loewner ranks. Frequency proximity alone is insufficient in a dense cluster.

13. Field localization, mode volume, Purcell factor, and beta factor

The present diagnostic integrates a driven-field electric-energy proxy over selected cell masks, schematically

\[ \eta_R=\frac{\sum_{\mathbf r\in R}|E(\mathbf r)|^2} {\sum_\mathbf r|E(\mathbf r)|^2}. \]

It shows whether the driven response is concentrated in the protected air strip or inner defect. It is not a rigorously normalized quasinormal mode and cannot by itself give mode volume.

For a nondispersive, closed, lossless cavity one often writes

\[ V_\mathrm{eff}(\mathbf r_a) =\frac{\int\epsilon(\mathbf r)|E(\mathbf r)|^2d^3r} {\epsilon(\mathbf r_a)|\hat{\mathbf d}\cdot E(\mathbf r_a)|^2} \]

and

\[ F_P=\frac{3}{4\pi^2} \left(\frac{\lambda}{n}\right)^3\frac{Q}{V_\mathrm{eff}}. \]

For this open, radiating, eventually dispersive structure, a QNM normalization or a directly calibrated Green-tensor LDOS is safer. The atom objective should eventually use the projected Green tensor at the actual atom position and dipole orientation, possibly averaged over trap uncertainty.

The implemented first atom-aware objective uses that direct Green response. For the centered z-oriented source at the atomic frequency,

\[ F_a(p)=\frac{\Im[d_a^T A(f_a,p)^{-1}d_a]} {\Im[d_a^T A_{\rm hom}(f_a)^{-1}d_a]}. \]

The homogeneous calculation retains the same grid, source, PML, and longitudinal boundary but sets n_core == n_clad. Its denominator is fixed during topology optimization, so the gradient is exactly the fixed-frequency log-LDOS gradient from Section 7. In the weak-coupling one-mode limit this optimizes a quantity proportional to \(g^2/\kappa\), hence Purcell enhancement or cooperativity for a fixed atomic linewidth. It does not isolate bare \(g\); a pole-residue or QNM-normalized objective is required for that claim. The complete estimator and independent acceptance contract are recorded in ATOM_COUPLING_OBJECTIVE.md.

Total decay should be separated into the desired guided-waveguide channel and all other channels:

\[ \Gamma_\mathrm{tot}=\Gamma_\mathrm{wg}+\Gamma_\mathrm{rad} +\Gamma_\mathrm{abs}, \qquad \beta_\mathrm{wg}=\frac{\Gamma_\mathrm{wg}}{\Gamma_\mathrm{tot}}. \]

Equivalently for a cavity,

\[ \kappa=\kappa_\mathrm{wg}+\kappa_\mathrm{parasitic}, \qquad Q=\frac{\omega_r}{\kappa}. \]

Maximizing total Q alone can suppress useful waveguide outcoupling. A later objective should combine atom-position coupling \(g\) or Green LDOS, parasitic Q, and a target external coupling/critical-coupling condition. This requires guided-mode flux projections or port self-energies not yet present.

14. Fabrication continuation

The current density filter controls numerical length scale but does not prove fabricability. A credible fabrication stage should add:

beta continuation toward a binary design;
explicit minimum solid and minimum void radii;
erosion/nominal/dilation scenarios sharing one design;
slab-thickness and refractive-index uncertainty;
sidewall angle, corner rounding, and etch bias;
connectivity and mechanically isolated-feature checks; and
robust objectives evaluated at all manufacturing scenarios.

The protected air strip remains a hard mask in every scenario.

15. Approximation and trust ledger

Component	Mathematical status	Current trust condition
Yee curl-curl discretization	controlled finite-difference approximation	`curl(grad)=0`, reciprocity, gradient tests, grid convergence
transformation PML	finite absorbing-layer approximation	material clearance and transverse-PML invariance
periodic x boundary	exact for the discretized eventually periodic lead	condensation/full-period agreement, Dyson residual, two-period translation
infinitesimal lead loss	retarded-sheet selector	convergence as selector tends to zero in a band gap
topology filter/projection	chosen parameterization	exact transpose gradient and hard masks
fixed-frequency LDOS gradient	analytic for the discrete problem	directional finite difference
local pole continuation	first-order real-axis approximation	center convergence, overlap tracking, rational cross-check
implicit pole-Q gradient	analytic for the local continued problem	directional finite difference and independent step validation
AAA pole	data-driven rational continuation	holdouts, residues, model order, window, linewidth
Loewner gradient	differentiable frozen-subspace surrogate	rebuild after step; AAA decides acceptance
shared-pole fit	nonlinear system-identification model	probes, folds, polynomial/pole order, SIMO cross-check
energy fractions	driven localization proxy	never interpret as QNM mode volume

16. Reproducibility checklist

For any claimed improvement record:

complete config.json and source checkpoint;
code version or archive SHA-256;
Python/package versions and BLAS thread settings;
grid, unknown counts, mirror period/phase, and PML clearances;
air-strip maximum density and reflection error;
boundary/Dyson/translation and gradient checks when relevant;
wide/local frequency arrays and all seven complex responses;
validator metrics, held-out poles, and acceptance booleans;
PML and grid comparison tables;
runtime and maximum resident set size; and
the exact distinction between spectral, PML, grid, and device claims.

The code and reports in this package provide worked examples of each stage.

17. Shifted cooperativity and constrained moving-asymptote updates

For post-nucleation optimization, define \(\Omega(p)=\operatorname{Re}\widetilde f_p(p)\), where \(\widetilde f_p\) is the differentiable Loewner continuation from Section

The proposal objective and its total topology derivative are

\[ J_s(p)=\log\frac{\operatorname{Im}h(\Omega(p),p)} {\operatorname{Im}h_{\rm hom}(\Omega(p))}, \]

\[ \frac{dJ_s}{dp}= \left.\frac{\partial\log\operatorname{Im}h}{\partial p}\right|_\Omega +\left[ \frac{\partial\log\operatorname{Im}h}{\partial\Omega} -\frac{\partial\log\operatorname{Im}h_{\rm hom}}{\partial\Omega} \right]\operatorname{Re}\frac{d\widetilde f_p}{dp}. \]

The frequency derivatives are centered differences using only retarded real-axis solves. This adaptation never replaces the exact frequency-dependent mirror self-energy by a conventional complex-frequency linear eigenproblem and never directly decimates the lead below the real axis.

For \(\gamma=-\operatorname{Im}\widetilde f_p\), the signed atom detuning and its derivative are

\[ \delta_a=\frac{\operatorname{Re}\widetilde f_p-f_a}{\gamma}, \qquad \frac{d\delta_a}{dp}=\frac{\operatorname{Re}(d\widetilde f_p/dp)}{\gamma} +\frac{(\operatorname{Re}\widetilde f_p-f_a) \operatorname{Im}(d\widetilde f_p/dp)}{\gamma^2}. \]

A conservative separable moving-asymptote subproblem maximizes the shifted objective subject to nondecreasing fixed-atom LDOS, nondecreasing pole Q, the detuning band, physical density bounds, and an outer trust box. Several provisional real-axis gradient refreshes may occur inside one box, but only a fresh wide/local/confirmatory rational validation can accept the cumulative topology. Accepted steps additionally require strict independently evaluated improvements in fixed-atom LDOS, shifted LDOS, and shared-pole Q, followed by the normal checkpoint-bound PML gate. The complete estimator and trust contract are in SHIFTED_COOPERATIVITY_ARCHITECTURE.md.