Multi-channel computational electromagnetics

We need full-wave solutions of Maxwell’s equations to make quantitative predictions about wave propagation and scattering in multi-channel systems such as disordered media, aperiodic metasurfaces, diffractive optical elements, channel [de]multiplexers, and photonic circuits. However, these simulations are often prohibitively slow given the multiple length scales and the large number of input channels involved. Each input channel requires the solution of one distinct scattering problem. The Optics in Complex Systems group aims to overcome this computational bottleneck.

Conventional methods to solve Maxwell’s equations numerically discretize the system on a basis and then solve for all the unknowns in this basis set. However, such a basis may contain millions to billions of elements, while we only need to extract a relatively small number of data points from it (e.g., the amplitude of light scattered to the directions or waveguide modes of interest). An additional loop over the many input channels makes the inefficiency of a full-basis solution even more severe, as illustrated in the figure below.

We develop a more efficient strategy that computes only the quantities of interest. Any set of linear responses in the frequency domain can be written as a generalized scattering matrix CA^-1B where A discretizes the Maxwell differential operator, the columns of matrix B contain different input source profiles of interest, and the rows of matrix C contain different output projection profiles of interest. Instead of evaluating the full-basis solution X = A^-1B, we directly compute the orders-of-magnitude-smaller projected matrix CA^-1B through a partial LU factorization of an augmented matrix K = [[A, B]; [C, 0]] as shown in the figure below. Doing so is equivalent to a partial Gaussian elimination that projects the Green’s function A^-1 onto the inputs B and outputs C, bypassing the forward and backward substitutions that solve for the many unknown variables in A^-1B. This approach, which we call “augmented partial factorization” (APF) [1], is ideal for multi-channel problems since it jointly handles all inputs without a loop over them. It is exact and applicable to any discretization scheme in any dimension.

APF can reduce the computing time by three to seven orders of magnitude while lowering memory usage (see figure below) without compromising accuracy [1]. This drastic speed-up enabled us to perform full-wave simulations of scattering matrix tomography and other imaging methods [2] and demonstrate the enhanced correlations in the backscattering direction for entangled photon pairs scattered from disorder [3].

We also developed multi-channel gradient computation using APF and used it for the inverse design of aperiodic metasurfaces [4], as shown below for a broad-angle beam splitter.

We developed an open-source APF solver called MESTI (Maxwell’s Equations Solver with Thousands of Inputs) in Julia [5] and MATLAB [6]. It uses finite-difference discretization on the Yee grid. MESTI supports any relative permittivity profile, any set of input source profiles, any set of output projection profiles, and comes with PML and the common boundary conditions. Multithreading and MPI are supported via the MUMPS linear-solver package.

Related publications

1. Fast multi-source nanophotonic simulations using augmented partial factorization, Ho-Chun Lin, Zeyu Wang, and Chia Wei Hsu. Nature Computational Science 2, 815–822 (2022).
2. Full-wave simulations of tomographic optical imaging inside scattering media, Zeyu Wang, Yiwen Zhang, and Chia Wei Hsu. arXiv:2308.07244.
3. Coherent backscattering of entangled photon pairs, Mamoon Safadi, Ohad Lib, Ho-Chun Lin, Chia Wei Hsu, Arthur Goetschy, and Yaron Bromberg. Nature Physics 19, 562–568 (2023).
4. Fast multichannel inverse design through augmented partial factorization, Shiyu Li, Ho-Chun Lin, and Chia Wei Hsu. ACS Photonics 11, 378–384 (2024).
5. github.com/complexphoton/MESTI.jl
6. github.com/complexphoton/MESTI.m