PaperlayerView exact revision on arXiv

Source-linked AI summary

Fast Nanophotonic Inverse Design using Precomputed Numerical Green Functions

Thi Phuong Thao Nguyen, James Wang, Constantine Sideris

arXiv:2610.02821v1physics.opticscs.NAmath.NAphysics.comp-ph

TL;DR

Nanophotonic inverse design relies on repeated full-wave simulations, making optimization computationally expensive. The paper extends PNGF to nanophotonic design, using precomputed numerical Green functions and low-rank updates to evaluate designs and gradients. In a broadband directional-coupler example, PNGF-enabled DBS achieved a speedup of more than 2400× over FDTD evaluation.

  • Problem

    Repeated full-wave simulations for objective evaluations and gradients make inverse-design optimization costly, particularly for large or 3D structures.

    Source: §1 · p. 2

  • Method

    The paper extends PNGF to nanophotonic inverse design, treating design changes as scatterers and evaluating their effects with precomputed quantities and low-rank updates.

    Source: §1 · p. 2, §2 · p. 3, §4 · p. 9

  • Results

    More than 2400× speedup over Lumerical FDTD was achieved for broadband directional-coupler optimization using PNGF-accelerated DBS.

    Source: §4.2 · p. 11

  • Takeaways & Limitations

    The method reproduces full-wave solutions without approximation and can be applied to a wide range of photonic and other electromagnetic devices.

    Source: §6 · p. 17

    The method has not yet been extended to fully vectorial three-dimensional problems, and APF precomputation for large 3D problems may face significant factorization-memory demands.

    Source: §6 · p. 17

Abstract

from arXiv · show

Inverse design has been transformative in nanophotonics, providing an automated means of realizing high-performance, non-intuitive devices. However, the primary bottleneck of existing inverse design approaches is their reliance on time-consuming and computationally expensive full-wave electromagnetic simulations. In this work, we extend the Precomputed Numerical Green Function (PNGF) method for the first time to nanophotonic inverse design, accelerating evaluation of the forward problem by multiple orders of magnitude with no loss in accuracy compared to the full-wave solution. After a single, fully parallelizable precomputation step, both the objective function and its gradient can be evaluated through simple calculations whose cost scales linearly with the size of the design region. A low-rank matrix update technique further reduces the cost of objective function evaluation, yielding sub-millisecond computation times per iteration. Design examples using both gradient-based level-set optimization and tile-flipping direct binary search often achieve more than three orders of magnitude speedup in design time, resulting in an ultrafast photonic inverse design platform that can design new devices from scratch in seconds.

1. Introduction

The paper addresses the cost of repeated full-wave simulations in nanophotonic inverse design by extending PNGF to evaluate modified structures and gradients using precomputed quantities. The method reduces computation to the design region and reports design-time speedups of three to four orders of magnitude in its examples.

1. Introduction

  • Repeated full-wave simulations for objective and gradient evaluations are a major computational cost in inverse design, especially for large or 3D structures.Machine-learning alternatives also require many simulations for training and may be inaccurate in underrepresented regions.Source: §1 · p. 2
  • Three to four orders of magnitude faster design times than conventional finite-difference methods are reported for the examples, reducing hours to days of computation to seconds.The examples include a 2D power splitter and a broadband directional coupler.Source: §1 · p. 3

2. Fast Structure Evaluation using Precomputed Numerical Green Functions

  • PNGF treats each design-region modification as a scatterer and adds its equivalent-source field to the original field to obtain the modified structure's field exactly.The approach exploits the fact that much of the simulation domain remains static during optimization.Source: §2 · p. 3, Figure 1 · p. 4

2.1. Evaluation of Modified Structures via Low-Rank Matrix Updates

  • Precomputed numerical Green functions let low-rank updates evaluate candidate modifications exactly, with cost scaling linearly in the design-region size.The matrices encode design-region responses, objective outputs, and incident fields.Source: §2 · p. 3, §2.1 · p. 5

2.2. Precomputation of the Numerical Green Functions

  • The Green-function matrices can be computed from independent unit-source simulations or obtained together through augmented partial factorization, then reused across optimization runs.For the paper's examples, precomputation took under 2 minutes for the power splitter and under 20 minutes for the broadband coupler.Source: §2.2 · p. 7

2.3. Gradient Calculation

  • The gradient for all design variables costs O(Ndes) operations per excitation–objective pair, avoiding the additional adjoint simulation used in conventional gradient-based design.It is obtained from the current PNGF matrices at negligible additional cost.Source: §2.3 · pp. 7–8

3. Runtime Scaling of PNGF

PNGF evaluation time is independent of the surrounding simulation-domain size when the design region is fixed, while conventional FDTD and FDFD runtimes grow with the domain. With a fixed simulation domain, PNGF evaluation scales linearly with the number of design parameters on one core; apparent sublinearity on 128 cores reflects underused parallelism on smaller problems.

3. Runtime Scaling of PNGF

  • With a fixed design region, PNGF evaluation time stays constant as the surrounding domain grows, while FDTD and FDFD runtimes increase.This makes PNGF especially advantageous for structures embedded in large simulation environments.Source: Figure 3 · p. 8, §3 · p. 8
  • On one core, PNGF evaluation time scales linearly with design-parameter count, reaching millisecond-range or lower times per evaluation.The tested design region grows from 10 × 10 to 50 × 50 tiles.Source: Figure 3 · p. 8, §3 · p. 9
  • The 128-core trend appears sublinear for smaller problems because they do not fully use the available parallelism.Source: Figure 3 · p. 8, §3 · p. 9

4. Accelerated DBS using PNGF

PNGF-accelerated Direct Binary Search (DBS) evaluates tile changes efficiently, supporting inverse design of a 1550 nm power splitter and a broadband directional coupler. The coupler achieves near-even splitting with high transmission across 1490–1610 nm, while its optimization is reported to be over 2400 times faster than using FDTD.

4. Accelerated DBS using PNGF

  • DBS flips discrete design tiles and retains a change when it improves the objective; PNGF updates are performed only for accepted flips.The more expensive PNGF matrix update is avoided for rejected flips.Source: §4 · p. 9

4.1. 1550 nm 50:50 Power Splitter with Improved Connectivity

  • The 1550 nm splitter uses a 2.4 µm × 2.4 µm design region tiled into 30 × 30 silicon-or-air elements, with mirror symmetry enforced.Its PNGF precomputation took under two minutes.Source: §4.1 · p. 9
  • For the 1550 nm splitter, standard DBS reached >99.4% total transmission in 12.1 seconds, while connectivity-enhanced DBS reached >99.1% in 12.2 seconds.The modified design has far fewer isolated single-tile features and better connectivity; the figure compares structures and their field distributions.Source: §4.1 · pp. 9–10, Figure 4 · p. 10

4.2. Broadband Directional Coupler using PNGF-Accelerated DBS

  • The broadband coupler uses a 4 µm coupling region and a 924-parameter design, optimizing transmission across sampled wavelengths while penalizing uneven splitting and low total transmission.Seven wavelengths were sampled; its precomputation took less than 11 minutes total.Source: §4.2 · pp. 10–11
  • Over 1490–1610 nm, the directional coupler averages a 50.61:49.39 splitting ratio and 98.03% total transmission; PNGF optimization is over 2400× faster than FDTD.The PNGF run took 82.9 seconds; the same optimization using Lumerical FDTD on the same machine was estimated at approximately 55 hours.Source: §4.2 · p. 11, Figure 7 · p. 12

5. Level-Set Method

The level-set approach represents device boundaries with a continuous function whose zero contour defines the structure, and uses localized knot updates to retain PNGF’s sparse-update structure. It produces a 50:50 splitter and a broadband coupler, with the latter combining high transmission and balanced splitting with a smooth design and at least 35× speedup over conventional full-wave solvers.

5. Level-Set Method

  • The level-set function’s zero contour defines the device boundary, and its permittivity profile is projected onto the level-set grid and volumetrically averaged onto the simulation grid.Radial basis functions on a coarse grid parameterize the function; cutoff radii localize each knot update and preserve PNGF’s sparse-update structure.Source: §5 · p. 13
  • The level-set boundary evolves incrementally along its normal direction according to the Hamilton–Jacobi equation.The zero contour separates the two permittivity regions, while the function is constructed from radial basis functions centered at coarse-grid control knots.Source: §5 · p. 13
  • The level-set optimization can start from a random guess, but may encounter near-nondifferentiability and local minima, especially from an intuitive but suboptimal initialization.Greedy gradient descent with momentum is used for the splitter to improve convergence and reduce the tendency to become trapped in local minima.Source: §5 · p. 13, §5.1 · p. 13

5.1. 50:50 Power Splitter

  • The level-set power splitter reached 99.1% total transmission, or 49.5% per port, in 314 iterations and 27.8 seconds.The optimization used greedy gradient descent with momentum; Fig. 8 shows the final structure and field distribution.Source: §5.1 · p. 13, Figure 8 · p. 14, Table 3 · p. 14

5.2. Broadband Directional Coupler

  • 98.25% average total transmission and a 50.4 : 49.6 average splitting ratio were achieved across the broadband coupler’s operating band, with at least 35× speedup over conventional full-wave solvers.The design converged in 367 iterations in under 5 minutes and has a smooth, fabrication-friendly coupling region.Source: §5.2 · p. 14, Figure 9 · p. 15, Table 4 · p. 15

6. Outlook

The paper concludes that PNGF reduces the computational problem from the full simulation domain to the design region and can speed up optimization without sacrificing accuracy. It reports compatibility with multiple parameterizations and numerical methods, while identifying fully vectorial 3D extension as future work.

6. Outlook

  • PNGF achieves speedups of up to several orders of magnitude over conventional full-wave simulation without compromising accuracy.The method is derived directly from Maxwell’s equations with no approximations.Source: §6 · p. 17
  • PNGF is compatible with incremental-update optimization and parameterizations including DBS, level sets, rotatable cells, grating couplers, and spline-based boundaries.The formalism also applies to finite-element and integral-equation discretizations, in addition to the demonstrated finite-difference schemes.Source: §6 · p. 17
  • Extending PNGF to fully vectorial 3D problems is future work, where APF precomputation may be memory-impractical for large problems.The paper proposes iterative methods such as GMRES, domain decomposition, and advanced preconditioning to mitigate the factorization-memory overhead.Source: §6 · p. 17

S1. Equivalence Theorem

The supplementary derivation replaces scatterers with equivalent polarization currents and expresses their fields through free-space Green functions. Its discretized form identifies the design-region and modified-pixel selections used to form the equivalent-source system.

S1. Equivalence Theorem

  • A scatterer’s effect is represented by an equivalent polarization current, separating the total field into incident and scattered components.The scattered field is expressed using free-space Green functions, and an integral equation solves for polarization current collocated with the scatterers.Source: S1. Equivalence Theorem · p. 20

S2. Discretized Equivalence Theorem

  • In the discretized nonmagnetic formulation, the total electric field satisfies AE = B, with A the wave operator and B the source profile.The operator is formed from a background permittivity profile; changes from that profile are represented as permittivity modifications.Source: S2. Discretized Equivalence Theorem · pp. 20–21

S3. Selection Matrices for Encoding Modifications

  • Selection matrices G and U encode the design region and modified pixels, allowing the equivalent-source system to be projected onto the modified pixels.The projected system is the equivalent-source system used in the main text.Source: S3. Selection Matrices for Encoding… · p. 21

S4. Woodbury Matrix Identity

The Woodbury identity computes an inverse after a rank-k matrix update through a smaller inverse, without approximation. Applied to PNGFs, rank updates allow the precomputation to be performed once and reused as the physical structure changes.

S4. Woodbury Matrix Identity

  • The Woodbury identity replaces inversion of an updated n × n matrix with inversion of a smaller matrix and matrix multiplications, giving the same exact result under its stated rank and invertibility conditions.The required conditions include rank-k factors and a nonzero determinant for the updated matrix.Source: S4. Woodbury Matrix Identity · p. 21

S5. Matrix Rank Update of the Precomputed Numerical Green Functions

  • For structure changes during optimization, rank updates refresh the precomputed PNGFs rather than rerunning the expensive precomputation stage.The updated inverse is substituted into the PNGF definitions, so recurring full-inverse terms reduce to the precomputed matrices.Source: S5. Matrix Rank Update of the… · p. 22

S6. Computational Complexity Analysis

Tile-flip evaluation scales linearly with design-region size under the stated size assumptions, while updating the system matrices scales quadratically. The mode-overlap projection maps fields to modal amplitudes, whose squared magnitudes give coupled optical power.

S6. Computational Complexity Analysis

  • Under Nobj, Ninc, M ≪ Ndes, tile-flip evaluation has O(Ndes) complexity, while updating the system matrices scales as O(Ndes^2).The evaluation bound follows from grouping matrix operations; the matrix-update cost is dominated by a dense multiplication.Source: S6. Computational Complexity Analysis · p. 23

S7. Mode Overlap Integral

  • The modal fields are sampled or interpolated onto the simulation Yee grid, and the target mode profile is obtained from a waveguide cross-section eigenmode problem.For a propagation-invariant structure, the field dependence is written as exp(−iβ_m z).Source: S7. Mode Overlap Integral · p. 23
  • The mode-overlap integral extracts a desired forward output-mode amplitude from the simulated field using the target mode profile.The discretized overlap is represented by matrix S; the magnitude squared of the complex overlap coefficient gives power coupled into the target mode.Source: S7. Mode Overlap Integral · p. 24

S8. Modified Direct Binary Search (DBS) for Improved Connectivity

The modified DBS combines connectivity-aware pixel selection with a cleanup stage to reduce isolated features while preserving device performance within a prescribed objective tolerance.

S8. Modified Direct Binary Search (DBS) for Improved Connectivity

  • Modified DBS combines connectivity-aware search with post-processing cleanup to address isolated silicon islands and air holes in standard DBS designs.The method prioritizes pixels adjacent to existing air holes, encouraging connected void growth before considering background pixels.Source: S8. Modified Direct Binary Search… · p. 24
  • After convergence, cleanup flips isolated pixels only when the objective degradation stays within a prescribed tolerance.Each pixel is checked against its neighbors, and the objective is reevaluated for a proposed flip.Source: S8. Modified Direct Binary Search… · pp. 24, 26
  • The search accepts improving pixel flips, restarts the scan after an accepted change, and stops when no single-pixel flip improves the objective.Pixels adjacent to air holes are considered before other pixels; standard DBS instead selects unchecked pixels without that stated preference.Source: S8. Modified Direct Binary Search… · pp. 25–26

S9. Level-set parameterization

The level-set design uses Gaussian radial basis functions centered on a coarse knot grid, then limits each control-knot update and its gradient effects to a local cutoff region.

S9.1. Radial Basis Functions (RBF) for level-set Function Construction

  • The level-set function is built as a sum of Gaussian radial basis functions centered at knots on a coarse grid.The spread parameter is set to 2 or 3 times the knot-grid resolution to provide overlap between neighboring knots and avoid ill-conditioned systems or overly smoothed designs.Source: S9.1. Radial Basis Functions (RBF) for… · pp. 26–27

S9.2. Single-Knot Update with Cutoff Radius

  • Updating one control coefficient affects only pixels within a cutoff-radius ball centered on its knot.The cutoff radius is 3 or 4 times the knot-grid resolution, so the entire design region need not be updated.Source: S9.2. Single-Knot Update with Cutoff… · p. 27

S9.3. Gradient with Respect to the Parameterization

  • The gradient with respect to control knots is updated element-wise over the local cutoff subsets rather than recomputed across the full grid.Each level-set point belongs to only a small number of these subsets.Source: S9.3. Gradient with Respect to the… · p. 27

S9.4. Volumetric Averaging

The supplement describes volumetric averaging for mapping level-set permittivity to simulation cells and gives inverse-design examples for bends, arbitrary-ratio splitters, and a mode converter using rotatable DBS.

S9.4. Volumetric Averaging

  • Simulation-cell permittivity is obtained by volume- or area-averaging the level-set permittivity, with a matrix D mapping level-set values to simulation-cell values.The derivative with respect to the control coefficients is obtained by downsampling the level-set-grid derivative.Source: S9.4. Volumetric Averaging · p. 27

S10.1. Waveguide Bend

  • 99.3% transmission efficiency was achieved for the waveguide bend with standard DBS; modified DBS reached approximately 99.5%.The standard-DBS design took 17 seconds; the modified-DBS optimization took 24.7 seconds.Source: S10.1. Waveguide Bend · p. 28, Figure S1 · p. 28

S10.2. Arbitrary Power Splitter

  • Standard DBS examples target 30:70, 20:80, and 10:90 power-splitting ratios, while the modified-DBS examples incorporate connectivity considerations.The standard-DBS structures and their normalized field distributions are shown in Fig. S2; modified-DBS structures and distributions are shown in Fig. S3.Source: S10.2. Arbitrary Power Splitter · p. 29, Figure S2 · p. 30, Figure S3 · p. 32

S10.3. Mode Converter using Rotatable DBS

  • 97.5% conversion efficiency was achieved for the TM0-to-TM1 mode converter after 15 sweeps over 190 cells, in 36.1 seconds.The input waveguide is positioned off-center to introduce asymmetry between the symmetric TM0 input mode and antisymmetric TM1 output mode.Source: S10.3. Mode Converter using Rotatable… · p. 33, Figure S4 · p. 33

S11. Additional Inverse Design Examples using the Level-Set Method

Level-set optimization produced a high-throughput waveguide bend and a compact, efficient mode converter; Y-branch initialization also reduced the iterations and runtime needed to refine power splitters.

S11.1. Waveguide Bend

  • The waveguide bend delivered almost 99.5% of input power to the output port.The 2.4 × 2.4 μm² junction reached this performance after 370 iterations in 33.3 seconds.Source: S11.1. Waveguide Bend · pp. 33–34

S11.2. Mode Converter

  • The mode converter reached 98.6% conversion efficiency in a 1.4 × 2.4 μm² region after 730 iterations in 28.8 seconds.Efficiency was calculated at the output using mode overlap.Source: S11.2. Mode Converter · p. 34

S11.3. Arbitrary-Ratio Power Splitters with a Y-Junction Initial Guess

Related papers

Loading 2610.02821v1…