Derivation of the radiation-aware 2.5D optimizer

1. Scope, normalization, and time convention

All lengths are normalized by the free-space target wavelength $\lambda_0$. The target ordinary frequency is $f_0=1$, angular frequency is

\[ \omega=2\pi f, \]

and $c=1$. Fields have the time convention $e^{-i\omega t}$. A passive quasinormal-mode (QNM) pole therefore lies in the lower half plane,

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

with

\[ Q=-\frac{\omega_r}{2\omega_i}. \]

The physical dielectric is one patterned slab of thickness $h$, with an in-plane density $\rho(x,y)$ that is constant through the slab:

\[ \epsilon(x,y)=\epsilon_c+\Delta\epsilon\rho(x,y), \qquad \Delta\epsilon=n_\mathrm{core}^2-n_\mathrm{clad}^2. \]

The homogeneous half spaces above and below have $\epsilon_c=n_\mathrm{clad}^2$. The present implementation is a scalar polarization reduction. Its purpose is rapid, radiation-aware topology discovery. The full-vector 3D Maxwell code remains authoritative.

2. Scalar open-slab equation

Inside the patterned layer, the scalar field obeys

\[ -\left(\nabla_\parallel^2+\partial_z^2\right)\psi =\lambda\epsilon(x,y)\psi, \qquad \lambda=\omega^2. \tag{1} \]

This equation can represent one decoupled polarization exactly only in special geometries. Here it is a controlled reduced model: it retains diffraction, vertical standing-wave structure, all retained light-cone channels, and their outgoing radiation loss, while omitting vector polarization mixing.

The x/y domain is a periodic supercell. It contains the finite free-form defect, the immutable air strip, and several explicit fixed mirror periods on both sides. A sufficiently long mirror makes interaction between periodic cavity copies exponentially small. Mirror-length convergence must still be checked for a device candidate.

3. In-plane Fourier expansion

At zero supercell Bloch wavevector, use orthonormal plane waves

\[ \phi_{\mathbf G}(x,y)=\frac{1}{\sqrt A} e^{i\mathbf G\cdot\mathbf r_\parallel}, \]

where

\[ \mathbf G=\left(\frac{2\pi m_x}{L_x}, \frac{2\pi m_y}{L_y}\right). \]

The retained orders are $-M_x\le m_x\le M_x$ and $-M_y\le m_y\le M_y$. The dielectric multiplication matrix is

\[ E_{\mathbf G\mathbf G'} =\frac{1}{A}\int_A \epsilon(x,y)e^{-i(\mathbf G-\mathbf G')\cdot\mathbf r_\parallel} d^2r. \tag{2} \]

For a cell-centred density grid with $N_p$ pixels, the code evaluates the Galerkin projection directly:

\[ E=\frac{1}{N_p}\Phi^\dagger \operatorname{diag}(\boldsymbol\epsilon)\Phi, \tag{3} \]

where $\Phi_{p\mathbf G}=e^{i\mathbf G\cdot\mathbf r_p}$. This construction is explicitly symmetrized against roundoff and is Hermitian for real density. It also makes the exact pixel derivative transparent.

4. Vertical finite-element basis

Only the physical slab interval $-h/2\le z\le h/2$ is meshed. Let $N_n(z)$ be continuous piecewise-linear nodal functions. Expand

\[ \psi(x,y,z)=\sum_{\mathbf G,n} c_{\mathbf G n}\phi_{\mathbf G}(x,y)N_n(z). \tag{4} \]

On a vertical element of length $\Delta z$, the stiffness and mass matrices are

\[ K_z^{(e)}=\frac{1}{\Delta z} \begin{bmatrix}1&-1\\-1&1\end{bmatrix}, \qquad M_z^{(e)}=\frac{\Delta z}{6} \begin{bmatrix}2&1\\1&2\end{bmatrix}. \tag{5} \]

After assembly, the closed-slab volume matrices are

\[ K_{\mathbf G n,\mathbf G'm} =\delta_{\mathbf G\mathbf G'} \left[(K_z)_{nm}+|\mathbf G|^2(M_z)_{nm}\right], \tag{6} \]

and

\[ M_{\mathbf G n,\mathbf G'm} =E_{\mathbf G\mathbf G'}(M_z)_{nm}. \tag{7} \]

Equation (7) is the source of topology coupling between Fourier harmonics.

Exact centered-emitter parity sector

For the enforced reflection-symmetric density and a centered scalar emitter, the coupled field is even under $x\mapsto-x$, $y\mapsto-y$, and $z\mapsto-z$. Let $P_x$, $P_y$, and $P_z$ contain normalized even reflection pairs (and the unpaired zero harmonic or middle FEM node). The production atom-coupling basis can therefore use

\[ P=P_x\otimes P_y\otimes P_z, \qquad K_r=P^\dagger K P, \qquad M_r=P^\dagger M P, \qquad \Sigma_r=P^\dagger\Sigma P. \tag{S} \]

This is an invariant-subspace restriction, not a symmetry penalty or a field approximation. The source is projected with the same orthonormal map and fields are lifted through $P$ before evaluating real-space energies and topology derivatives. On the regression basis, the reduced and full atom responses agree to 2.09e-15 relative error, their gradients agree to 2.05e-15, and their field correlation is unity to roundoff. The current 10 x 4 x 9 Fourier/FEM launch basis falls from 1,701 to 275 unknowns (6.185x); the ratio approaches eight as all three orders grow.

5. Exact vertical outgoing boundary

Outside the slab, every in-plane Fourier harmonic propagates independently. For channel $\mathbf G$, define

\[ k_{z,\mathbf G}(\lambda) =\sqrt{\epsilon_c\lambda-|\mathbf G|^2}. \tag{8} \]

An outgoing upper-half-space field is proportional to $e^{+ik_z(z-h/2)}$. The outgoing lower-half-space field is proportional to $e^{-ik_z(z+h/2)}$. In both cases the outward-normal derivative is

\[ \partial_n\psi=i k_z\psi. \tag{9} \]

Integrating Equation (1) by parts gives the weak boundary contribution

\[ -\int_{\partial\Omega}v^*\partial_n\psi,dS =-i k_z v^*\psi. \]

Thus the two surfaces contribute the channel-diagonal self-energy

\[ \Sigma_{\mathbf G}(\lambda) =-ik_{z,\mathbf G}(\lambda) \left(e_0e_0^T+e_{N_z-1}e_{N_z-1}^T\right). \tag{10} \]

For a propagating channel, $k_z$ is primarily real and Equation (10) has a negative imaginary part: energy leaves the slab. For an evanescent channel, $k_z=i\alpha$ with $\alpha>0$, so $-ik_z=+\alpha$ is a real reactive boundary impedance.

No PML thickness, strength, or onset is present. Equation (10) is exact for the scalar equation, homogeneous cladding, and retained Fourier basis.

6. Riemann sheet

The square root in Equation (8) has a branch point at every channel threshold. A QNM calculation must remain on one outgoing analytic sheet. At the target frequency $\lambda_t=(2\pi f_t)^2$, classify a channel as open when

\[ |\mathbf G|^2<\epsilon_c\lambda_t. \]

During nonlinear continuation:

open channels retain the root with positive real part;
closed channels retain the root with positive imaginary part.

This classification is frozen during the pole solve. Selecting a new sign from each intermediate $\lambda$ can jump sheets near a threshold and was observed to prevent basis convergence.

The derivative on the chosen sheet is

\[ \frac{d k_{z,\mathbf G}}{d\lambda} =\frac{\epsilon_c}{2k_{z,\mathbf G}}, \]

\[ \Sigma'_{\mathbf G}(\lambda) =-i\frac{\epsilon_c}{2k_{z,\mathbf G}} \left(e_0e_0^T+e_{N_z-1}e_{N_z-1}^T\right). \tag{11} \]

The implemented derivative agrees with a direct complex finite difference to relative error $3.46\times10^{-10}$.

7. Nonlinear eigenvalue problem

Combining the volume and boundary terms gives

\[ T(\lambda,\rho)c =\left[K+\Sigma(\lambda)-\lambda M(\rho)\right]c=0. \tag{12} \]

This is nonlinear in $\lambda$ because the radiation boundary depends on frequency. Treating $\Sigma$ as constant would give the wrong pole and the wrong topology gradient.

For a reference value $\lambda_n$, solve the frozen generalized problem

\[ \left[K+\Sigma(\lambda_n)\right]c =\mu(\lambda_n)M c \tag{13} \]

near the tracked pole. If $l$ is the corresponding left eigenvector,

\[ \mu'(\lambda_n) =\frac{l^\dagger\Sigma'(\lambda_n)c} {l^\dagger M c}. \tag{14} \]

Newton's method for $F(\lambda)=\mu(\lambda)-\lambda=0$ is

\[ \lambda_{n+1}=\lambda_n -\alpha\frac{\mu(\lambda_n)-\lambda_n} {\mu'(\lambda_n)-1}, \tag{15} \]

with damping $0<\alpha\le1$ and a bounded correction. The complex field overlap selects the same eigenpair from each frozen solve. The final reported residual is

\[ r_\mathrm{NEP} =\frac{\|T(\lambda)c\|_2} {\|[K+\Sigma(\lambda)]c\|_2}. \tag{16} \]

The accepted production and refined results have residuals $1.30\times10^{-8}$ and $4.50\times10^{-9}$.

8. Nonlinear eigenvalue topology derivative

Let $l$ satisfy

\[ l^\dagger T(\lambda,\rho)=0. \]

Differentiate Equation (12) for one real density perturbation:

\[ \left[\Sigma'(\lambda)-M\right]c\,d\lambda -\lambda(dM)c+T\,dc=0. \]

Premultiplying by $l^\dagger$ removes the field derivative because $l^\dagger T=0$:

\[ l^\dagger\left[\Sigma'-M\right]c\,d\lambda -\lambda l^\dagger(dM)c=0. \]

Therefore

\[ \boxed{ d\lambda =\lambda\frac{l^\dagger(dM)c} {l^\dagger[\Sigma'(\lambda)-M]c} }. \tag{17} \]

Equation (17) is the central adjoint identity. The $\Sigma'$ term is not optional; it is the normalization correction for a dispersive open boundary.

9. Pixel sensitivity without one matrix per pixel

Let the reconstructed right and left vertical coefficient vectors at pixel $p$ be

\[ R_n(p)=\sum_{\mathbf G}e^{i\mathbf G\cdot\mathbf r_p}c_{\mathbf G n}, \qquad L_n(p)=\sum_{\mathbf G}e^{i\mathbf G\cdot\mathbf r_p}l_{\mathbf G n}. \]

From Equation (3), the derivative of the dielectric Fourier matrix with respect to one density pixel gives

\[ l^\dagger\frac{\partial M}{\partial\rho_p}c =\frac{\Delta\epsilon}{N_p} \sum_{n,m}L_n^*(p)(M_z)_{nm}R_m(p). \tag{18} \]

All pixel sensitivities are therefore obtained from two Fourier reconstructions and a small vertical contraction. No per-pixel Maxwell solve or matrix assembly is required.

10. Q and frequency objective derivatives

Because $\omega=\sqrt\lambda$,

\[ d\omega=\frac{d\lambda}{2\omega}. \tag{19} \]

For $Q=-\Re\omega/(2\Im\omega)$,

\[ d\log Q =\frac{\Re(d\omega)}{\Re\omega} -\frac{\Im(d\omega)}{\Im\omega}. \tag{20} \]

The differentiable score is

\[ S=\log Q -w_f\left(\frac{f-f_t}{f_t}\right)^2 -w_g\frac{1}{N_d}\sum_{p\in D}\rho_p(1-\rho_p), \tag{21} \]

where the last term weakly discourages gray density. The complete raw-variable directional derivative, including the topology map below, agrees with finite differences to relative error $2.70\times10^{-5}$.

Pure Q maximization can reduce atom coupling. The current implementation therefore enforces protected-gap and defect energy fractions as hard step gates. A future version should differentiate a calibrated atom-position Green-tensor or QNM-volume objective instead.

11. Density filter, projection, and hard geometry

Raw design variables $p_j\in[0,1]$ exist only in the allowed defect region. The normalized compact filter is

\[ \widetilde p_i=\sum_jF_{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)}. \tag{22} \]

The smooth Heaviside projection is

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

with

\[ \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)]}. \tag{24} \]

For a full density gradient $g_\rho$,

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

The exact transpose is regression-tested. Reflection projection averages the x, y, and xy reflected copies. Because that projection is self-adjoint, the same operation is applied to gradients.

After projection:

the complete region $|x|<0.5\lambda_0$ is reset to zero;
terminal mirror buffers and explicit periodic mirrors are reset to their fixed binary geometry.

Neither hard region is a design variable. The protected-strip maximum is exactly zero in every included checkpoint.

12. Mode metrics and identity

At each physical pixel, the energy proxy is

\[ U_p=\epsilon_p\Re\left[R^\dagger(p)M_zR(p)\right]. \tag{26} \]

It defines defect and protected-gap fractions used for gating. This is a localization diagnostic, not a rigorous dispersive QNM energy or mode volume.

On one basis, candidate identity uses the generalized energy overlap

\[ C_M=\frac{|c_a^\dagger M c_b|} {\sqrt{|c_a^\dagger M c_a|\,|c_b^\dagger M c_b|}}. \tag{27} \]

Across bases, reconstruct the complex midplane fields on the same physical density grid and use

\[ C_E=\frac{|E_a^\dagger E_b|} {\|E_a\|_2\|E_b\|_2}. \tag{28} \]

Equation (28) is invariant to arbitrary complex phase and different numbers of basis coefficients. Frequency is only a secondary hint. This gate was added after a smoke-basis Q=136.9 candidate proved orthogonal to the refined mode.

13. Trust region and acceptance

The raw-variable score gradient is normalized by its root-mean-square value. A bounded proposal is

\[ p_\mathrm{trial}=\Pi_{[0,1]} \left(p+\alpha\frac{\nabla S}{\operatorname{RMS}(\nabla S)}\right), \tag{29} \]

followed by reflection projection. The step is halved until all gates pass:

score increases;
Q is finite and positive;
same-basis overlap exceeds its threshold;
defect and protected-gap fractions exceed their thresholds;
frequency stays in the continuation window;
nonlinear residual is below threshold;
hard topology invariants remain exact.

The independent cross-basis validator does not differentiate or reuse the production eigenvector coefficients. It reconstructs field hints and solves the richer nonlinear eigenproblems from scratch.

14. Analytic slab validation

For a homogeneous slab and a fixed in-plane propagation constant $\beta$, the fundamental even scalar/TE mode satisfies

\[ q\tan\left(\frac{qh}{2}\right)=\alpha, \tag{30} \]

where

\[ q=\sqrt{\epsilon_\mathrm{core}\lambda-\beta^2}, \qquad \alpha=\sqrt{\beta^2-\epsilon_c\lambda}. \]

At $\beta=8$, the analytic normalized frequency is 0.7688154870. The vertical-FEM results are:

vertical nodes	frequency	relative error
5	0.7706122140	2.3370e-3
9	0.7692604928	5.7882e-4
17	0.7689264829	1.4437e-4

The approximately quadratic error reduction is consistent with linear finite elements and independently validates the DtN sign for evanescent channels.

15. What the validated optimization establishes

The production basis has 1,701 field unknowns. Four accepted steps give

\[ Q:15.3533\rightarrow17.3124\rightarrow20.3130 \rightarrow25.0269\rightarrow33.1158. \]

On the independently richer 3,025-unknown basis:

\[ Q:16.2539\rightarrow32.8077. \]

The gains are 2.1569 and 2.0184. Initial and final field correlations are 0.98387 and 0.98490. Therefore the same reduced-model mode improves on both bases.

This establishes that the 2.5D optimizer is possible and working. It does not establish an absolute physical Q because the scalar reduction omits vector polarization and the z-oriented atomic Green tensor.

16. Route to a vector 2.5D formulation

The next reduced-model upgrade should replace the scalar basis with several vector slab modes. Schematically,

\[ \mathbf E(x,y,z)=\sum_{\nu,\mathbf G} c_{\nu\mathbf G}\, \mathbf e_{\nu\mathbf G}(z)e^{i\mathbf G\cdot\mathbf r_\parallel}. \]

The homogeneous cladding boundary becomes a dyadic TE/TM channel admittance. The nonlinear-eigenvalue derivative retains the same abstract form as Equation (17): only $K$, $M$, and $\Sigma$ become block-vector matrices. The topology filter, trust region, checkpointing, physical-grid field correlation, and cross-basis acceptance architecture can be reused.

Until that vector model is independently validated, every promising scalar topology must be exported to the full-vector 3D Maxwell solver and reacquired without carrying its Q.