Differentiable Waveguide Method

• Differentiable Waveguide Method is a set of formulations that make waveguide analysis differentiable with respect to design variables in optimization and analog computing contexts.
• Techniques include eigenmode solvers for modal inverse design, differentiable diffraction methods for EUV lithography, and waveguide circuits that approximate temporal differentiation.
• These methods enhance design efficiency by leveraging automatic differentiation and structured operator mappings to reduce computational cost and improve accuracy.

“Differentiable Waveguide Method” denotes a family of waveguide-centered formulations in which either the electromagnetic waveguide solver itself is made differentiable with respect to design variables, or the waveguide structure is engineered so that its scattering response realizes a differentiation operator. In contemporary arXiv literature, the phrase is used in at least three technically distinct senses: differentiable eigenmode-based inverse design of optical waveguides, end-to-end differentiable waveguide diffraction solvers for inverse lithography, and waveguide-junction architectures whose transfer function approximates temporal derivatives (Gray et al., 2024, Es'kin et al., 24 Jun 2026, MacDonald et al., 2023). A broader operator-theoretic lineage also includes scattering-matrix, Dirichlet-to-Neumann, Fourier, and nonlinear-eigenvalue formulations that are not themselves presented as autodiff frameworks but expose structured matrix maps that are amenable to sensitivity analysis (Roddick, 2017, Jarlebring et al., 2015, Goodwill et al., 25 Sep 2025).

1. Terminological scope and conceptual variants

The term does not refer to a single standardized algorithm. In one usage, it describes a differentiable electromagnetic waveguide mode solver for inverse design of modal dispersion, where the full chain

$$\text{design parameters} \rightarrow \text{geometry} \rightarrow \varepsilon(\mathbf r,\omega) \rightarrow \text{eigenmodes} \rightarrow \text{dispersion metrics} \rightarrow \text{objective}$$

is made compatible with automatic differentiation (Gray et al., 2024). In another, it denotes a differentiable implementation of the classical waveguide diffraction method for EUV masks, where Fourier expansion, layerwise eigenproblems, global mode coupling, reflected-field synthesis, and wafer-plane loss evaluation are treated as a single differentiable computation graph (Es'kin et al., 24 Jun 2026). In a third, quite different sense, it denotes a waveguide-based analogue temporal differentiator in which interconnected stubs and junctions are designed so that the transmission coefficient reproduces the local frequency dependence of a temporal derivative operator (MacDonald et al., 2023).

Usage	Core object	Representative paper
Differentiable mode solver	Eigenmodes, dispersion, adjoint pullbacks	(Gray et al., 2024)
Differentiable waveguide diffraction	Layered-mask modal solver with reverse-mode AD	(Es'kin et al., 24 Jun 2026)
Waveguide temporal differentiator	Scattering synthesis of $(j\omega)^n$-type response	(MacDonald et al., 2023)

This multiplicity matters because the phrase combines two distinct notions of differentiability. In the inverse-design literature, “differentiable” means differentiable with respect to parameters inside an optimization graph. In the analogue-computing literature, “differentiable” refers to a device that performs mathematical differentiation on a signal. A plausible implication is that any encyclopedic treatment must separate solver differentiability from operator realization.

2. Differentiable eigenmode solvers for modal inverse design

A central modern meaning of the differentiable waveguide method is the differentiable eigenmode framework developed for inverse design of waveguide dispersion (Gray et al., 2024). The method is motivated by the claim that existing inverse-design tools based on finite-difference-time-domain models are poorly suited for optimizing waveguide modes for adiabatic transformation or perturbative coupling, especially in dispersion engineering of optical waveguides involving broad bandwidths and anisotropic, frequency-dependent dielectric response. The key technical contribution is custom back-propagation through an electromagnetic eigenmode solver so that gradients of objectives depending on mode indices, fields, group index, and GVD can be computed efficiently with adjoint sensitivities (Gray et al., 2024).

The forward model begins from parameterized waveguide geometry and general tensor dielectric models. In the demonstrated example, the waveguide is a rib in thin-film lithium niobate parameterized by top/core width, core thickness, and normalized etch fraction, with sidewall angle fixed by fabrication assumptions. The solver explicitly supports anisotropic dielectric tensors and frequency-dependent dispersion, and it uses analytically precomputed

$$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$

because group index and GVD depend on those derivatives (Gray et al., 2024).

A critical requirement is smooth geometry-to-grid dependence. To achieve this, the method applies anisotropic subpixel dielectric smoothing based on Kottke et al. Near an interface with unit normal $\hat n$ between materials $\varepsilon_1$ and $\varepsilon_2$, the smoothed tensor $\tilde\varepsilon$ is constructed through a differentiable map involving local rotation into the interface-aligned frame and matrix-valued transforms $\tau,\tau^{-1}$. The framework then computes

$$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$

using precompiled Jacobian and Hessian information for the smoothing map (Gray et al., 2024). This smoothing stage is described as essential to make the geometry-to-material map differentiable and stable under optimization.

The eigenproblem is adapted from the plane-wave expansion formulation used in MIT Photonic Bands. The fixed-$k$ problem is

$(j\omega)^n$0

and the physically relevant fixed-$(j\omega)^n$1 problem becomes nonlinear because $(j\omega)^n$2 depends on frequency (Gray et al., 2024). The solver handles the fixed-$(j\omega)^n$3 case by iteratively solving fixed-$(j\omega)^n$4 problems and updating $(j\omega)^n$5 with Newton’s method.

A notable feature is that group index and GVD are computed from a single mode solution rather than by finite differences in frequency. The group index is expressed through mode-energy terms, the $(j\omega)^n$6-derivative of the operator, and explicit material-dispersion terms involving $(j\omega)^n$7; GVD is then obtained via differentiation of the group-index observable and $(j\omega)^n$8 (Gray et al., 2024). This avoids finite-difference noise and reduces the number of required eigen/adjoint solves.

The differentiable core lies in custom AD pullbacks for the fixed-$(j\omega)^n$9 and fixed-$$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$0 eigenproblems. For the fixed-$$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$1 Hermitian eigensystem, the method derives sensitivities with respect to both $$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$2 and $$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$3, including explicit dependence on upstream sensitivities to eigenvalue and eigenvector. For the nonlinear fixed-$$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$4 problem, it introduces the nonlinearity via an equality constraint and differentiates implicitly (Gray et al., 2024). The expensive part of the gradient computation is then one adjoint solve per relevant mode/frequency rather than one solve per parameter, so the cost scales with objective-dependent mode solves, not with the number of geometric parameters.

The paper states both theoretically and empirically that the back-propagation cost is effectively independent of parameter count, labeling the method as $$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$5 mode solves per iteration and reporting a roughly constant gradient/primal cost ratio as the number of parameters increases from $$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$6 to $$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$7 (Gray et al., 2024). In the demonstrated TFLN SHG problem, optimization of three design parameters converged in eight steps, reduced computational cost by about $$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$8 relative to exhaustive search, and identified broadband optical frequency-doubling designs in the $$\frac{d\varepsilon}{d\omega}, \qquad \frac{d^2\varepsilon}{d\omega^2},$$9–$\hat n$0m range (Gray et al., 2024). For center wavelengths $\hat n$1 nm, $\hat n$2 nm, $\hat n$3 nm, and $\hat n$4 nm, the reported quasi-phase-matching FWHM bandwidths were approximately $\hat n$5–$\hat n$6 nm, or $\hat n$7–$\hat n$8 THz, for $\hat n$9 cm interaction lengths (Gray et al., 2024).

3. Differentiable waveguide diffraction for inverse lithography

A second major meaning of the term is the differentiable waveguide method for EUV mask inverse lithography (Es'kin et al., 24 Jun 2026). Here the “waveguide method” is the classical rigorous diffraction solver for layered periodic masks, but reformulated as an end-to-end differentiable forward model. The target problem is inverse lithography for reflective masks at $\varepsilon_1$0 nm, with the mask consisting of a single absorber layer above a multilayer mirror and the design variable taken as the absorber permittivity distribution $\varepsilon_1$1 in 2D, later generalized to $\varepsilon_1$2 in 3D (Es'kin et al., 24 Jun 2026).

The forward model is governed by layerwise Helmholtz equations with periodicity in $\varepsilon_1$3, TE-polarized oblique incidence, and interface continuity of tangential fields. Each patterned layer is represented in a Fourier basis,

$\varepsilon_1$4

and yields a generalized eigenvalue problem

$\varepsilon_1$5

Mode coupling across layers produces a global linear system

$\varepsilon_1$6

from which reflected free-space diffraction orders and wafer-plane intensity are synthesized (Es'kin et al., 24 Jun 2026).

The differentiable aspect is explicit: the paper states that the map from optimization variables to reflected field is completely determined by a sequence of algebraic operations and is therefore differentiable with respect to the design variables. The full chain includes absorber-density parameterization, construction of $\varepsilon_1$7, Fourier coefficients, assembly of the layer eigenproblem, eigendecomposition, assembly of the global linear system, solution for amplitudes, reflected-field synthesis, intensity formation, and loss computation. Gradients are obtained not by finite differences and not by a hand-derived adjoint, but by reverse-mode automatic differentiation through the full solver, implemented in PyTorch (Es'kin et al., 24 Jun 2026).

Two parameterizations are used. The pixel-wise density formulation writes

$\varepsilon_1$8

with $\varepsilon_1$9 pixels and a loss

$\varepsilon_2$0

A Fourier-parameterized projection instead represents a band-limited latent function and uses a spectral regularizer $\varepsilon_2$1 (Es'kin et al., 24 Jun 2026). This second parameterization reduces dimensionality and promotes smoother walls.

The experiments use realistic absorber materials TaBN, La, and U, with a multilayer mirror of 30 Ru/Be/Sr periods, absorber thickness $\varepsilon_2$2 nm, incidence angle $\varepsilon_2$3, and mask period $\varepsilon_2$4 nm (Es'kin et al., 24 Jun 2026). The target condition is intensity matching on a remote plane $\varepsilon_2$5 at $\varepsilon_2$6, with $\varepsilon_2$7 m in the experiments.

Reported results include preservation of target optical behavior after binarization of the relaxed density, suppression of side lobes, and material-dependent differences in image quality. For $\varepsilon_2$8 harmonics and 5000 epochs, the reported maximum reflected-to-incident intensity ratios are $\varepsilon_2$9 for TaBN, $\tilde\varepsilon$0 for La, and $\tilde\varepsilon$1 for U; qualitatively, Lanthanum gives the best central maximum, whereas Uranium gives the closest match to the desired field (Es'kin et al., 24 Jun 2026). For a three-strip target with $\tilde\varepsilon$2 harmonics, the Fourier-parameterized method reproduces essentially the same target as pixel-wise optimization with a wall-clock time reduction from $\tilde\varepsilon$3 s to $\tilde\varepsilon$4 s, reported as a $\tilde\varepsilon$5-fold speed-up (Es'kin et al., 24 Jun 2026). The framework is also extended to 3D with $\tilde\varepsilon$6, TaBN absorber, and 500 optimization epochs (Es'kin et al., 24 Jun 2026).

The same paper introduces the waveguide neural operator, which replaces the linear-system solve $\tilde\varepsilon$7 by a multilayer perceptron trained via the physics residual

$\tilde\varepsilon$8

The WGNO remains end-to-end differentiable and is described as mesh-independent and physics-informed, but the paper reports that it does not yet provide a speed advantage because simultaneous training outweighs the faster approximate solve (Es'kin et al., 24 Jun 2026). This suggests that, in this branch of the literature, the differentiable waveguide method is less a new approximation than a differentiable reformulation of an established rigorous solver.

4. Waveguide structures as differentiation operators

A third usage appears in electromagnetic analogue computing, where the waveguide itself implements differentiation of a signal envelope (MacDonald et al., 2023). The basic architecture is a set of parallel-plate waveguides treated as transmission lines supporting TEM propagation, with one input guide, one output guide, and $\tilde\varepsilon$9 side stubs connected at a junction, so the total number of interconnected waveguides is

$\tau,\tau^{-1}$0

At a parallel junction with equal characteristic impedances and subwavelength cross-section, scattering is described by

$\tau,\tau^{-1}$1

Closed-ended or open-ended stubs then return delayed copies of the signal with known phase conditions (MacDonald et al., 2023).

For a closed-ended stub of length $\tau,\tau^{-1}$2, the round-trip delay is

$\tau,\tau^{-1}$3

and the reflected wave acquires a sign inversion. In the four-waveguide example with two closed stubs, the output becomes

$\tau,\tau^{-1}$4

which is directly related to the finite-difference definition

$\tau,\tau^{-1}$5

Thus, for small $\tau,\tau^{-1}$6, the transmission behaves as a temporal differentiator on the envelope of a sinusoidally modulated input (MacDonald et al., 2023).

The central synthesis formula for $\tau,\tau^{-1}$7 identical stubs of equal length $\tau,\tau^{-1}$8 and reflection coefficient $\tau,\tau^{-1}$9 is

$$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$0

By choosing stub length and termination so that the transmission zero lies at the carrier frequency $$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$1, the device approximates the normalized narrowband differentiator response around $$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$2. For air-filled guides, the paper states that the required phase condition is met for closed stubs at $$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$3 and open stubs at $$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$4 (MacDonald et al., 2023).

The method extends to higher-order and fractional differentiation by cascading blocks and matching the overall transfer function to

$$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$5

Because each stage is reflective near its notch, exact cascade analysis uses the Redheffer star product rather than weak-coupling approximations (MacDonald et al., 2023). The paper further invokes the Riemann–Liouville derivative to interpret the fractional case in the time domain, while noting that the device actually approximates the corresponding frequency response over a finite bandwidth.

Numerical examples include an 8 GHz first-order differentiator, transmission-mode differentiation around 8 GHz and reflection-mode differentiation around 4 GHz in the same 2-stub structure, a second-order differentiator with working frequency range approximately $$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$6, and a fractional differentiator with order $$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$7 and usable bandwidth about $$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$8 (MacDonald et al., 2023). The paper is explicit that the method acts on the envelope of a narrowband sinusoidally modulated signal, not on arbitrary unrestricted broadband waveforms. In encyclopedic terms, this literature uses “differentiable waveguide method” to mean waveguide-based realization of a differentiation operator, not gradient-based optimization.

5. Operator-theoretic and numerical antecedents

Several papers do not use the phrase in the modern autodiff sense but supply structured formulations that are highly relevant to differentiable waveguide analysis. One line concerns scattering matrices. A method for Euclidean waveguides with cylindrical ends computes the stationary scattering matrix $$\frac{d\tilde{\varepsilon}}{d\omega}, \qquad \frac{d^2\tilde{\varepsilon}}{d\omega^2}$$9 and its arbitrarily high derivatives with respect to the spectral parameter $k$0 by gluing interior and exterior Neumann-to-Dirichlet maps and differentiating the resulting operator system recursively (Roddick, 2017). The central formula

$k$1

is then supplemented by an explicit recursive formula for $k$2 (Roddick, 2017). This is differentiability with respect to spectral parameter rather than geometry, but it establishes a rigorous framework for derivative-aware waveguide scattering.

A related domain-decomposition approach for metallic waveguide circuits constructs local impedance-to-impedance maps on subdomains and merges them through Schur-complement coupling, producing the global scattering matrix as an explicit matrix composition (Goodwill et al., 25 Sep 2025). The paper does not present an adjoint or autodiff method, but its divide-and-conquer structure exposes the global response as a composition of modal data, local BIE solves, and sparse block inverses. The paper explicitly notes mode cutoff and adaptive basis truncation as obstacles to smooth differentiability (Goodwill et al., 25 Sep 2025). This suggests that differentiable waveguide solvers benefit from fixed modal bases and parameter regimes away from thresholds.

A second antecedent is the nonlinear eigenvalue formulation of periodic waveguides. The tensor infinite Arnoldi method and its waveguide specialization WTIAR write the waveguide problem as

$k$3

with polynomial FEM blocks and DtN boundary blocks that depend nonlinearly on $k$4 through square-root terms (Jarlebring et al., 2015). The solver repeatedly uses exact derivatives of $k$5 with respect to the eigenvalue, together with recurrence formulas for the derivatives of the DtN coefficients. The paper is not an inverse-design work, but it exposes a highly structured derivative-aware computational graph for eigenvalue-dependent waveguide analysis (Jarlebring et al., 2015).

A third lineage appears in Fourier and DtN formulations. Fourier methods for harmonic scalar waves in general waveguides map smooth blocks conformally to strips, encode shape through

$k$6

and derive stable Riccati equations for reflection and transmission matrices (Andersson et al., 2013). Stable DtN-map recursion for piecewise-uniform 2D optical waveguides similarly propagates local modified DtN maps segment by segment through operator recurrences (Wang et al., 2012). Neither paper claims autodiff, but both formulate the waveguide problem through explicit operator compositions that are plausible sensitivity-analysis backbones.

At a different geometric scale, the immersed interface method for cylindrical optical waveguides modifies only the finite-difference rows whose stencils cross refractive-index interfaces, using explicit jump matrices derived from Maxwell continuity conditions (Horikis, 2011). Because the resulting sparse eigenproblem depends algebraically on one-sided refractive indices and interface locations, this discretization is structurally favorable for differentiation, although the paper does not derive gradients and the interface-row assignment is only piecewise smooth (Horikis, 2011).

6. Applications, limitations, and present directions

The application range of differentiable waveguide methods is broad but stratified by physical model. In eigenmode inverse design, the demonstrated application is broadband SHG phase matching in etched TFLN waveguides, including quasi-TE$k$7 fundamental and second-harmonic modes, with $k$8 cm interaction length and broadband optimization over multiple frequencies (Gray et al., 2024). In diffraction-based lithography, the application is gradient-based inverse lithography of EUV reflective masks with realistic absorber materials and multilayer mirrors (Es'kin et al., 24 Jun 2026). In analogue computing, the application is first-order, higher-order, and fractional temporal differentiation of modulated signals using passive waveguide junctions (MacDonald et al., 2023). At system level, end-to-end differentiable design of geometric waveguide displays combines non-sequential Monte Carlo polarization ray tracing with a differentiable transfer-matrix thin-film solver to optimize more than one thousand coating thicknesses and billions of non-sequential ray-surface intersections, improving throughput from $k$9 to $(j\omega)^n$00 and improving eyebox and FoV uniformity by about $(j\omega)^n$01 and $(j\omega)^n$02, respectively (Yang et al., 7 Jan 2026). Although that work is not labeled a “differentiable waveguide method” in exactly the same sense, it belongs to the same differentiable-waveguide design ecosystem.

The main limitations recur across subfields. Differentiability often requires smoothing or fixed operator structure. In the eigenmode-dispersion framework, subpixel smoothing is required because otherwise mode frequencies and fields vary non-smoothly with geometry (Gray et al., 2024). In scattering and modal methods, thresholds, degeneracies, and changing numbers of propagating modes create nonsmooth parameter dependence and can alter even the dimension of the scattering operator (Roddick, 2017, Goodwill et al., 25 Sep 2025). In inverse lithography, the authors identify differentiation through eigendecomposition, complex-valued linear solves, and post hoc hard binarization as sensitive steps (Es'kin et al., 24 Jun 2026). In waveguide analogue differentiation, the transfer-function approximation is inherently narrowband and accompanied by reflection and insertion loss because the device is passive and operates near a transmission zero (MacDonald et al., 2023).

A second limitation is model scope. The differentiable mode-solver work is demonstrated on 2D cross-sections, although the authors state that the formulation and adjoint method apply to 3D eigenmode models as well (Gray et al., 2024). The differentiable EUV waveguide method is built around periodic layered-mask structure, TE-polarized incidence, and target matching on a remote plane without a full projection-optics or resist model (Es'kin et al., 24 Jun 2026). The waveguide-display optimization keeps geometry fixed and optimizes mainly coating thicknesses, noting that mirror tilts and rotations remain future work (Yang et al., 7 Jan 2026).

The present trajectory of the field points toward broader end-to-end wavewave design graphs. One direction is clear expansion from hand-derived adjoint pullbacks to general differentiable physics engines; another is replacement of the most expensive linear algebra components by learned surrogates such as WGNO, albeit without current wall-clock advantage (Es'kin et al., 24 Jun 2026). A plausible implication is that future differentiable waveguide methods will combine rigorous operator structure with learned accelerators, while retaining explicit control of modal branches, boundary operators, and physically meaningful observables.

In summary, the differentiable waveguide method is best understood not as a single algorithm but as a research program: the recasting of waveguide analysis and design into differentiable operator pipelines, together with a parallel usage in which waveguide networks are designed to realize differentiation itself. Across these meanings, the unifying theme is the replacement of black-box waveguide analysis by structured mappings—modal, scattering, transfer, or diffraction mappings—whose internal algebra is exposed well enough to support either gradient propagation or direct synthesis of a differentiating transfer function (Gray et al., 2024, Es'kin et al., 24 Jun 2026, MacDonald et al., 2023).