MEENT: DIFFERENTIABLE ELECTROMAGNETIC SIMULATOR FOR MACHINE LEARNING

## Abstract

Abstract Electromagnetic (EM) simulation plays a crucial role in analyzing and designing devices with sub-wavelength scale structures such as solar cells, semiconductor devices, image sensors, future displays and integrated photonic devices.
Specifically, optics problems such as estimating semiconductor device structures and designing nanophotonic devices provide intriguing research topics with far-reaching real world impact.
Traditional algorithms for such tasks require iteratively refining parameters through simulations, which often yield sub-optimal results due to the high computational cost of both the algorithms and EM simulations.
Machine learning (ML) emerged as a promising candidate to mitigate these challenges, and optics research community has increasingly adopted ML algorithms to obtain results surpassing classical methods across various tasks.
To foster a synergistic collaboration between the optics and ML communities, it is essential to have an EM simulation software that is user-friendly for both research communities.
To this end, we present meent , an EM simulation software that employs rigorous coupled-wave analysis (RCWA). Developed in Python and equipped with automatic differentiation (AD) capabilities, meent serves as a versatile platform for integrating ML into optics research and vice versa.
To demonstrate its utility as a research platform, we present three applications of meent : 1) generating a dataset for training neural operator, 2) serving as an environment for the reinforcement learning of nanophotonic device optimization, and 3) providing a solution for inverse problems with gradient-based optimizers.
These applications highlight meent ’s potential to advance both EM simulation and ML methodologies.
The code is available at https://github.com/kc-ml2/meent with the MIT license to promote the cross-polinations of ideas among academic researchers and industry practitioners.

## 1 Introduction

Figure: Figure 1: Summary. Simulation Algorithm depicts the process flow of electromagnetic simulation algorithm, namely RCWA, in meent. Applications present the most representative problems that meent can be utilized.
Refer to caption: x1.png

Harnessing light-matter interaction to design or analyze a device with sub-wavelength scale structure has a wide range of applications, including high-efficiency solar cells [, ], ultra-thin metalenses and displays [, ], optical metrology for semiconductor fabrication [, ], X-ray diffraction for material analysis [, ], optical computation [, ], and so on. Their implication to the real world is far-reaching, leading to improved renewable energy production, enhanced user experience, and next-generation computation.
Electromagnetic (EM) simulation plays a crucial role in such applications, which also poses a challenging problem due to its time-consuming nature for precise calculation [] and iterative characteristic of optimization.

Machine learning (ML) is a promising candidate to solve such a problem, and with the advent of deep learning, optics community has been made successful efforts [, , ] to leverage modern machine learning techniques to find both better optimization algorithms [, , , , ] and faster electromagnetic simulators [, , ].
These researches show significant potential of ML for uncovering new insights and expediting scientific discoveries. On the other hand, from the viewpoint of ML, designing or analyzing such devices are ideal environments for developing new ML algorithms, as they provide both ample simulation data to train ML models and well-defined optimization goals, often called figure of merit (FoM), for the real world applications.

However, the integration of ML into computational optics presents several challenges. Traditional EM simulation software are often written in languages like C or MATLAB and therefore are not easily compatible with ML frameworks which are predominantly developed in Python and require the compatibility with automatic differentiation (AD).
Furthermore, the scarcity of public data in certain domains, such as semiconductor fabrication industry, compounded by stringent intellectual property regulations, poses significant obstacles, especially to ML researchers lacking domain expertise.
Overcoming these barriers requires innovative approaches in generating and sharing data to enable ML researchers to explore new frontiers in computational optics.

In response to these challenges, we introduce meent, a Python-native differentiable EM simulator. meent is based on rigorous coupled-wave analysis (RCWA) [, , , , , , , ], a high-throughput, deterministic EM simulation algorithm that is widely adopted in optics across academia and industries.
And meent is a full-wave simulator, that approaches the exact solution with more Fourier modes included.

Key features of meent include its compatibility with automatic differentiation (AD) [, ] for modeling and optimizing devices in a continuous space.
AD compatibility in ML toolchain is pivotal, and while some existing tools support vector modeling and others support AD, none offer both functionalities simultaneously.
For developer ergonomics, meent is developed to be compatible with three different backend frameworks: NumPy [], JAX [], and PyTorch [].
By supporting multiple backends, meent facilitates easy adoption among researchers with varying levels of domain expertise and different backend preferences.

We showcase the utility of meent with various applications of ML to optics. First we present how to use meent to analyze and design a metasurface, whose sub-wavelength scale structure is carefully designed to achieve unprecedented control of light [, , ]. We also illustrate how to use meent in optical metrology [], one of the most successful industrial applications within the semiconductor fabrication, that serves to estimate the dimensions of device structures between process steps, thereby effectively monitoring excursions and maximizing yield due to its non-destructive nature and high throughput capabilities.

By enabling each user to generate datasets tailored to specific research needs, meent can democratize the access to EM simulation data. We hope that meent will facilitate collaboration between ML and optics researchers and thereby accelerate scientific discovery in computational optics.

Our contributions are summarized as follows:

- •
Development of meent, a Python-native EM simulation software under MIT license supporting automatic differentiation and continuous space operation in ML frameworks.
- •
Demonstration of meent’s versatility with
examples of ML algorithms, including Fourier neural operator, model-based RL, and gradient-based optimizers.
- •
Documentation of meent with detailed explanations and instructions.

## 2 Related Work

##### EM simulation algorithms.

There exist several methods for full-wave
EM simulation, each offering distinct advantages.

Finite difference time domain (FDTD) operates within the real space and time domain, employing the finite difference method. It discretizes space into grids and iteratively solves the function at these grid points over successive time steps
[].
Notable open-source software packages include
[Meep](https://github.com/NanoComp/meep) [],
[gprMAX](https://github.com/gprMax/gprMax) [],
[OpenEMS](https://github.com/thliebig/openEMS) [],
[ceviche](https://github.com/fancompute/ceviche) [] and
[FDTD++](http://www.fdtdxx.com).

The Finite element method (FEM) is a general technique used to solve differential equations across diverse fields.
Like FDTD, it functions within real space but supports both frequency and time domain. The key advantage of FEM is its capability to handle complex geometries through the utilization of unstructured grids, called meshes [].
Two prominent commercial tools are Ansys HFSS and COMSOL Multiphysics.
For open-source software, there are
[FEniCS](https://github.com/FEniCS) [],
[Elmer FEM](https://github.com/ElmerCSC/elmerfem) [],
[deal.II](https://github.com/dealii/dealii) [] and
[FreeFEM](https://github.com/FreeFem/FreeFem-sources) [].

In the realm of RCWA simulators, which operates in Fourier space and frequency domain, several open-source software options are available, each with its unique characteristics and advantages. [Reticolo](https://arxiv.org/abs/2101.00901) [] and [S4](https://github.com/victorliu/S4) [] stand out as classical open-source software, built upon a strong foundation of physics principles. These tools, implemented in MATLAB and C++ respectively, have earned recognition and have been extensively utilized in numerous research endeavors, establishing themselves as prominent representatives among open-source RCWA tools.
With the emergence of ML, the significance of Python-native code has grown substantially, prompting optics researchers to familiarize themselves with Python and its associated technologies. In recent years, there has been a proliferation of Python-based RCWA codes. [MAXIM](https://www.sciencedirect.com/science/article/abs/pii/S0010465521000163) [], for instance, offers a graphical user interface tailored to assist novice users in electromagnetic simulations. Additionally, [rcwa-tf](https://github.com/scolburn54/rcwa_tf) [], [gRCWA](https://github.com/weiliangjinca/grcwa) [], and [TORCWA](https://github.com/kch3782/torcwa) [] are notable for their support of automatic differentiation, a feature particularly valuable for gradient-based optimization tasks.

##### ML applications in optics.

Assisted with physical simulators, ML is being actively embraced across scientific domains to substitute heavy simulations with deep models that serve as surrogate solvers, offering increased robustness to hidden noise.
Seminal works such as physics-informed neural network (PINN) [, , ] and neural operators [, , , , , , ] showed their potential as surrogate EM solvers [, ]. Reinforcement learning (RL) also showed its efficacy in the scientific domains, such as magnetic control of tokamak plasmas [] and classical mechanics [, ].

Representative examples of surrogate EM solver include MaxwellNet [], an instance of PINN. Fourier neural operator was used in [], where optimization of nanophotonic device [] is also conducted. Deep generative model was used in [] to reduce computational cost compared to traditional optimization algorithm, and the feasibility of using model-free RL was demonstrated in [, , ]. Our example explores the possibility of applying model-based RL to device optimization, rooted on RNN-based world model [, , , , ].

## 3 Meent: electromagnetic simulation framework

##### Electromagnetic simulation algorithm.

meent uses rigorous coupled-wave analysis (RCWA) [, , , , , , , ], which is based on Faraday’s law and Ampére’s law of Maxwell’s equations [],

$$ $\nabla\times\mathbf{E}=-j\omega\mu_{0}\mathbf{H},\qquad\nabla\times\mathbf{H}= j\omega\varepsilon_{0}\varepsilon_{r}\mathbf{E},$ (1) $$

where $\mathbf{E}$ and $\mathbf{H}$ are electric and magnetic field in real space, $j$ denotes the imaginary unit number, i.e. $j^{2}=-1$, $\omega$ denotes the angular frequency, $\mu_{0}$ is vacuum permeability, $\varepsilon_{0}$ is vacuum permittivity, and $\varepsilon_{r}$ is relative permittivity.
As illustrated in Figure 1, it is a technique used to solve PDEs in Fourier space, aiming to estimate optical properties such as diffraction efficiency or field distribution.
We reserve the detail of RCWA for Appendices A and G for readers interested in delving into the fundamentals of RCWA.

##### Geometry modeling.

In order to simulate real-world phenomena, it is important to faithfully replicate them within a computational environment. In this context, geometries, or device structures, must be transferred and modeled to facilitate the execution of virtual experiments through simulation.

Figure: (a) Representation
Refer to caption: extracted/5658627/sections/meent_framework/images/modeling-geometry.png

meent offers support for two modeling types: raster and vector. Analogous to the image file format, raster-type represents data as an array, while vector-type utilizes a set of objects, with each object comprising vertices and edges, as shown in Figure 2(a).
Due to their distinct formats, each method employs different algorithms for space transformation, resulting in different types of geometry derivatives, including topological and shape derivatives, as depicted in Figure 2(b).
The topological derivative yields the gradient with respect to the permittivity changes of every cell in the grid, while the shape derivative provides the gradient with respect to the deformations of a shape.

These two modeling methods have distinct advantages and are suited to different applications - raster for freeform metasurface design and vector for optical critical dimension (OCD) metrology.
Freeform metasurface design is the realm of device optimization that typically using raster-type representation consisting of two materials.
The process involves iteratively altering the material of individual cells to find the optimal structure, where gradients may help effectively determining which cell to modify.
OCD involves estimating the dimensions of specific design parameters, such as width or height of a rectangle,
aligning naturally with vector-type and shape derivative.
Using raster in this scenario severely limits the resolution of the parameters due to its discretized representation.
In mathematical terms, raster operates in a discrete space, whereas vector operates in a continuous space, as illustrated in Figure 2(c).

##### Program sequence.

We provide a detailed explanation of the functions in meent and discuss the simulation sequence with examples in Appendix H. It explains the function codes and technical details that help users understand the inner workings and follow step by step.

##### Computation performance.

meent supports three numerical computation libraries as backends: NumPy, JAX, and PyTorch.
The computational performance evaluated
by backend, by architecture (64-bit and 32-bit), and by computing device (CPU and GPU) are provided in Appendix I as detailed benchmark tests. We emphasize that there is no single setup that is universally optimal, therefore it is recommended to conduct preliminary benchmark tests to find the best configuration before proceeding to the main experiments.

##### Applications.

Optics-focused applications are introduced in Appendix J, detailing the inverse design of nanophotonic devices aimed at optimizing device performance through EM simulation and automatic differentiation.

## 4 Meent in action: machine learning algorithms applied to optics problems

Here, we present how meent can be used in applying machine learning (ML) to optics problems.
First, we explore neural PDE solvers for Maxwell’s equations, using meent as a data generator.
Then we delve into device optimization through reinforcement learning (RL), utilizing meent as an RL environment.
Lastly, we address inverse problems within the semiconductor metrology domain, leveraging meent as a comprehensive solution for both simulation and optimization.

Figure: Figure 3: Metagrating and its abstract representations. (a) An example of metagrating sized 16 cells, where the grating layer is bounded by air and glass layers. The grating pattern periodically repeats along x axis and all layers remains constant along y axis. (b) Abstract matrix representation $u\in{\{1,-1\}}^{k\times k}$, the front view of metagrating. (c) Array representation of the grating pattern, $g\in{\{1,-1\}}^{k}$.
Refer to caption: extracted/5658627/sections/applications/images/metasurface2.png

Throughout Sections 4.1 and 4.2, analysis and design of metagrating beam deflector are performed.
A metagrating is a specific type of metasurface that is arranged in a periodic pattern and is primarily used to direct light into specific angle $\theta$ as shown in Figure 3. At the grating layer, a material is placed on $k$ uniform cells, and has the constraint of minimum feature size (MFS). MFS refers to the smallest contiguous grating cells the device can have.
Our figure of merits (FoMs) from the beam deflector include deflection efficiency $\eta\in[0,1]$ and $x$-component of electric field $\mathbf{E}$.

### 4.1 Fourier neural operator: prediction of electric field distribution

We provide two representative baselines of neural PDE solvers: (a) Image-to-image model, UNet [] and operator learning model, Fourier neural operator (FNO) []. The solvers learn to predict electric field, given a matrix representation of a metagrating.

##### Problem setup.

Our governing PDE that describes electric field distribution can be found by substituting left-hand side into right-hand side of Equation 1,

$$ $\begin{split}\nabla\times\nabla\times\mathbf{E}&=\omega^{2}\mu_{0}\varepsilon_ {0}\varepsilon_{r}\mathbf{E}.\\ \end{split}$ (2) $$

The aim of the neural PDE solver is to predict the electric field from the Maxwell’s equation in transverse magnetic (TM) polarization case.
Here we consider only $x$-component of the electric field, denoted as $v$.

Firstly, under the constraint of MFS, device patterns are sampled from uniform distribution, $g=[e_{1},...,e_{k}]\sim\mathrm{Unif}(\{-1,1\})$, and laid on matrix $u$ in Figure 3.
The corresponding ground truth $v$ is then calculated with meent by solving Equation 2.
For each physical condition, we generated 10,000 pairs of $(u,v)$ as a training dataset.
$u$ and $v$ have the same dimension, $x\in\mathbb{R}^{k\times k}$, and both are discretized on regular grids.
MFS was chosen as 4 which is more granular than 8 in the previous work [], which demonstrates that users can easily generate tailor-made datasets with meent.

Let $\mathcal{O}$ be an operator that maps a function $u(x)$ to a function $v(x)$ describing electric field, s.t. $v(x)=\mathcal{O}(u)(x)$.
An approximator of $\mathcal{O}$ represented by a neural network is updated using various losses, and since the solution space of a PDE is highly dependent on physical conditions, we assessed the robustness of baseline models across nine conditions and collectively report in Appendix D.2.

Figure: (a) Electric field
Refer to caption: x2.png

##### Fourier neural operator.

The effectiveness of FNO for solving Maxwell’s equation in our metagrating beam deflector is exhibited in Figure 4(a).
We follow techniques from [], in which original FNO is adapted to light scattering problem by applying batch normalization [], adding zero-padding to the input and adopting Gaussian linear error unit (GELU) activation [].
We further improved FNO’s parameter efficiency by applying Tucker factorization [], where a model’s weight matrices are decomposed into smaller matrices for low-rank approximation.
In addition to field prediction capability, we also show zero-shot super-resolution (trained in lower resolution, tested on higher resolution) capability in Figure 4(b), which is claimed to be a major contribution of FNO [].

Remarkably, FNO outperformed UNet by a significant margin (76% lower mean error) with only 1/10 parameters of UNet when both models were trained with $L_{2}$ loss, see D.2 for more detail.
The moderate performance of UNet in other PDE solvers [, ] contrasts with its poor performance in our task, which we attribute to its inability to capture detailed structures around the grating area. More information on model training is provided in Appendix D.1.
Additionally, in terms of wall time per electric field calculation, FNO required 0.23 seconds for inference, whereas physical simulation took 0.84 seconds, implying the utility of neural operator as a surrogate EM solver.

### 4.2 Model-based reinforcement learning: metasurface design

Here we demonstrate that meent can be used as an environment to train a model-based reinforcement learning (RL) agent, whose model learns how EM field evolves according to the change of the metagrating structure.
The objective of an RL agent here is to find the metagrating structure that yields high deflection efficiency, by flipping the material of a cell between silicon and air.
Here, MFS is set to 1, i.e. there is no MFS constraint.

##### Problem setup.

Here an RL agent undergoes fully-observable Markov decision process described as sequence of tuples $\{s_{t},a_{t},r_{t},s_{t+1}\}_{t=1}^{T}$, where the next state $s_{t+1}$ is determined by the dynamics model $p(s_{t+1}\ |\ s_{t},a_{t})$ and the action is the index of cells to flip, $a_{t}\in\{0,1,...,k-1\}$.
The state $s_{t}$ and the reward $r_{t}$ depends on which RL algorithm is used and will be defined shortly after.
Throughout the decision process, the agent learns to flip cells that maximizes the deflection efficiency $\eta_{t}$.

For training purpose, we implemented simple Gymnasium [] environment called
deflector-gym(^1^11[https://github.com/kc-ml2/deflector-gym](https://github.com/kc-ml2/deflector-gym)), which is built on top of of meent.
Given an input action $a_{t}$, the environment modifies current structure $g_{t}$ and outputs FoMs such as deflection efficiency ${\eta}_{t}$ or electric field $v_{t}$ in deterministic manner.

##### Model-based RL vs Model-free RL.

One way to categorize RL algorithms is whether it has an explicit dynamics model $p\approx p_{\theta}(s_{t+1}\ |\ s_{t},a_{t})$, where $p_{\theta}$ is some neural network.
Model-based RL (MBRL) agent utilizes the model $p_{\theta}$ to produce simulated experiences, from which the policy is improved [, ].
Therefore, MBRL agent is considered to be more sample efficient than model-free agent, i.e.  requires less interactions with actual environment to train. Sample efficiency can be crucial in computational science when the simulation cost is high.

For an MBRL algorithm, we chose DreamerV3 [], the first algorithm that solved ObtainDiamond task of MineRL [].
DreamerV3 is based on recurrent state-space model (RSSM) [] for modeling dynamics in the latent space. Developed to be robust at varying scale of observations and rewards across different tasks, it solved numerous tasks with a single set of hyperparameters, most of which were reused here as well.
We compare this with Deep Q Network (DQN) [], a model-free algorithm, adapted from [].

Figure: Figure 5: Learning curve across 3 random seeds. Historical maximum efficiencies during training phase are plotted.
Refer to caption: extracted/5658627/sections/applications/images/dv3-training-curve.png

Same as in [], our DQN agent receives the reward $r_{t}=\eta_{t}-\eta_{t-1}$ and observes the grating pattern as the state $s_{t}=g_{t}$. On the other hand, DreamerV3 agent receives reward $r_{t}=\eta_{t}$ and observes the grating pattern plus the electric field, $s_{t}=(g_{t},v_{t})$, to enable the dynamics model $p_{\theta}$ to learn underlying physics of the transitions of electric fields. For further training details and the motivation behind the above reward engineering, we refer to Appendix E.1

As was the case in other tasks. DreamerV3 agent showed improved sample efficiency in our task as seen in Figure 5 when compared to DQN.
Not only does it optimizes the structure faster, but also achieves higher deflection efficiency.
As a side remark, we mention that the training of RL agents can be massively accelerated with the parallelization of meent with Ray/RLlib [].
Simple comparison of training speed between single worker and multiple workers are shown in Appendix E.1.

Electric field prediction of the dynamics model.
In addition to the sample efficiency, another advantage of using MBRL in this task is that the dynamics model $p_{\theta}(s_{t+1}\ |\ s_{t},a_{t})$ can be used to predict the electric field of the next state, which not only makes the dynamics model interpretable but also suggest another way of developing an EM surrogate solver in addition to neural operators.

Figure: Figure 6: Comparison between ground truth and prediction. Rollout trajectory with the MBRL dynamics model (below). Given previous electric fields from step 1 to 5, the model predicts 25 future electric fields with actions the agent has actually taken, but without any access to ground truth electric field. All of the images are real part of the field.
Refer to caption: x4.png

Figure 6 shows the prediction of the model compared to the ground truth calculation from meent. One interesting observations is that, even when the prediction at a certain time step deviates from the ground truth, the model does not compound the error but is able to correct itself to converge to the correct estimation. The robustness of the prediction of the dynamics model is also illustrated in Appendix E.2, where the dynamics model was able to reproduce the correct electric field configuration even for a difficult problem that a neural operator fails to estimate correctly.

### 4.3 Inverse problem: OCD metrology

Figure: Figure 7: Schematic Diagram of Spectrum Fitting. This diagram illustrates the key components of the approach, including the vector of design parameters (P) from the real device, the set of spectra (S) derived from these parameters, the synthesized spectra at iteration i ($\mathbf{\hat{\textit{S}}_{\textit{i}}}$), and the estimated design parameters at the corresponding iteration ($\mathbf{\hat{\textit{P}}_{\textit{i}}}$). The methodology involves assessing the distance between $\mathbf{\hat{\textit{S}}_{\textit{i}}}$ and S using a loss function, with subsequent updates to $\mathbf{\hat{\textit{P}}_{\textit{i}}}$ to minimize this distance during each iteration.
Refer to caption: extracted/5658627/sections/applications/images/ocd-spectrum_fitting.png

A semiconductor device is a three-dimensional stack of layers, rendering direct measurement of parameters beneath the surface unfeasible without causing damage.
OCD offers a solution to this challenge by redirecting the observation target from the dimensions of physical device to its spectral characteristics (spectrum).
Consequently, OCD becomes an inverse problem: we deduce the dimensions of the structure in real space, which are the causal factors of spectrum shape, from observations.
The solution involves a probabilistic and iterative process known as spectrum fitting. This necessitates optimization in continuous space, which can be achieved using meent with vector modeling.

##### Spectrum fitting.

Figure 7 shows the process of finding solution using spectrum fitting.
The goal is to estimate P using S, as direct observation of P is impossible.
To achieve this, we initially create a virtual structure with limited prior knowledge provided by domain experts, and sample the initial parameters $\mathbf{\hat{\textit{P}}_{\textit{0}}}$ from a suitable distribution.
Subsequently, we generate $\mathbf{\hat{\textit{S}}_{\textit{0}}}$ from $\mathbf{\hat{\textit{P}}_{\textit{0}}}$ through simulation.
We then employ a distance metric as a loss function to quantify the discrepancy between S and $\mathbf{\hat{\textit{S}}_{\textit{0}}}$ to compare P and $\mathbf{\hat{\textit{P}}_{\textit{0}}}$.
The process is followed by backpropagation, which computes gradients of the distance with respect to each element of $\mathbf{\hat{\textit{P}}_{\textit{0}}}$. Following that, $\mathbf{\hat{\textit{P}}_{\textit{0}}}$ is updated to $\mathbf{\hat{\textit{P}}_{\textit{1}}}$.
This iterative process can be generalized using $\mathbf{\hat{\textit{P}}_{\textit{i}}}$ and $\mathbf{\hat{\textit{S}}_{\textit{i}}}$, where i denotes the iteration number.
As iterations progress, $\mathbf{\hat{\textit{S}}_{\textit{i}}}$ gradually converges towards S, and it is expected that $\mathbf{\hat{\textit{P}}_{\textit{i}}}$ will similarly approach the parameter set we seek to obtain, P.

##### Demonstration.

We now pivot our approach to meent. Rather than seeking P, we utilize meent to observe the behavior of optimization algorithms, a subject of keen interest for ML researchers. As an example, we introduce a case study involving eight design parameters with the details provided in Appendix F.
Employing spectrum fitting, we present the optimization curve of five distinct gradient-based algorithms: Momentum, Adagrad [], RMSProp [], Adam [] and RAdam [].
All algorithms share identical
$\mathbf{\hat{\textit{P}}_{\textit{0}}}$ to ensure consistency,
and evaluated repeatedly with 10 different initial conditions to mitigate the randomness associated with initial point location, a critical factor in local optimization algorithms.

Figure: Figure 8: Loss Curves of Various Gradient-based Algorithms for Spectrum Fitting. Each plot illustrates the change of loss over iterations, with the y axis represented in logarithmic scale.
Refer to caption: x5.png

Our purpose in this section is to demonstrate the capability of meent rather than introducing novel algorithms or achieving precise predictions. Accordingly, we provide a concise overview of the results.
Observing the optimization performance depicted in Figure 8, it is notable that Adagrad and RMSprop exhibit a more rapid initial decrease in loss. However, as the process progresses, the slope of the curves plateaus, and they finish at a relatively higher value (near 1E-4).
The Adam-family algorithms, namely Adam and RAdam, exhibit similarities in the early stages of optimization.
Furthermore, these
maintain the initial slope
during iterations, indicating that better fitting can be expected by increasing the number of iterations.
Adam achieves the optimum result overall, while RAdam demonstrates a more consistent performance across different cases, indicating superior robustness.

## 5 Conclusion

In our work, we introduce meent, a full-wave, differentiable electromagnetic simulation framework.
Through its capability for vector-type modeling and automatic differentiation, meent operates seamlessly within a continuous space, overcoming the limitations inherent in raster-type modeling for geometry representation.
We demonstrate with examples of applications how to use meent as a valuable tool to generate data for ML as well as a comprehensive solver for inverse problems, expanding its role beyond that of a simple electromagnetic simulator.
This versatility makes meent an invaluable framework to both machine learning and optics.

##### Limitations.

RCWA is restricted to time-harmonic cases, making it incapable of simulating scenarios where the fields do not exhibit steady-state sinusoidal behavior over time. However, this inherent algorithmic characteristic, which cannot be overcome, enhances RCWA’s speed, making it invaluable for cases that adhere to these constraints.

## Acknowledgements

We express our gratitude to Jean-Paul Hugonin, an external collaborator at Laboratoire Charles Fabry in Palaiseau, France, for his invaluable assistance in enhancing our understanding of RCWA and reticolo [],
Yannick Augenshtein at Flexcompute, for insightful discussion on computational optics,
Hoyun Choi at department of physics and astronomy, Seoul National University, for through overview of operator learning,
SheepRL [] team of Orobix, especially Federico Belotti (University of Milano-Bicocca) and Michele Milesi for discussion on the DreamerV3 code and the algorithm,
and Softmax team at KC for sharing their computational resources.
S.K., K.O., Chaejin P., J.S., S.N., Chanhyung P., J.P., S.H., and M.S.J. acknowledges the support by the MOTIE (Ministry of Trade, Industry & Energy, 1415180303) and KSRC (Korea Semiconductor Research Consortium, 20019357) support program for the development of the future semiconductor device.

## Appendix

## Appendix A A brief introduction to RCWA

RCWA is based on Faraday’s law and Ampére’s law of Maxwell’s equations [],

$$ $\nabla\times\mathbf{E}=-j\omega\mu_{0}\mathbf{H},\qquad\nabla\times\mathbf{H}= j\omega\varepsilon_{0}\varepsilon_{r}\mathbf{E},$ (3) $$

where $\mathbf{E}$ and $\mathbf{H}$ are electric and magnetic field in real space, $j$ denotes the imaginary unit number, i.e. $j^{2}=-1$, $\omega$ denotes the angular frequency, $\mu_{0}$ is vacuum permeability, $\varepsilon_{0}$ is vacuum permittivity, and $\varepsilon_{r}$ is relative permittivity.
After Fourier transform, $\mathbf{E}$ and $\mathbf{H}$ in real space become $\boldsymbol{\mathfrak{S}}$ and $\boldsymbol{\mathfrak{U}}$ in Fourier space, respectively, and Equation 3 then turns into

$$ $\displaystyle(-j\tilde{k}_{z})\boldsymbol{\mathfrak{S}}=\boldsymbol{{\Omega}}_ {L}\boldsymbol{\mathfrak{U}},\qquad(-j\tilde{k}_{z})\boldsymbol{\mathfrak{U}}= \boldsymbol{{\Omega}}_{R}\boldsymbol{\mathfrak{S}},$ (4) $$

where the matrices $\boldsymbol{\Omega}_{L}$ and $\boldsymbol{\Omega}_{R}$ are composed of wavevector matrices and convolution matrices that lie in Fourier space, $\tilde{k}_{z}$ is a normalized wavevector in the $z$-direction.
We can find $\boldsymbol{\mathfrak{S}}$ by merging Equation 4 in a single matrix equation as

$$ $\displaystyle(-j\tilde{k}_{z})^{2}\boldsymbol{\mathfrak{S}}=\boldsymbol{{ \Omega}}_{LR}^{2}\boldsymbol{\mathfrak{S}},$ (5) $$

where the matrix $\boldsymbol{\Omega}_{LR}^{2}$ is a matrix product of $\boldsymbol{\Omega}_{L}$ and $\boldsymbol{\Omega}_{R}$.
As the form implies, this equation can be solved by eigendecomposition of $\boldsymbol{\Omega}_{LR}^{2}$ to obtain the eigenvectors $\boldsymbol{\mathfrak{S}}$ and the eigenvalues $(-j\tilde{k}_{z})$. Then by substituting the eigenvectors for $\boldsymbol{\mathfrak{S}}$ in Equation 4, the corresponding solution of $\boldsymbol{\mathfrak{U}}$ can be obtained.

$\boldsymbol{\mathfrak{S}}$ and $\boldsymbol{\mathfrak{U}}$, which we have just computed, represent the set of electromagnetic modes (field representation in Fourier space) within a given medium.
To understand their interaction with incoming and outgoing light, we employ boundary conditions to ascertain their respective weights of the modes, or in other words, the coefficients in linear combination of the modes. These coefficients describe the extent of each mode’s influence on the overall field distribution.
Notably, coefficients at the input and output interfaces are designated as diffraction efficiencies, also called the reflectance and transmittance, serving as the primary purpose of RCWA.
Subsequently, the inverse transformation from Fourier space to real space enables the reconstruction of the field distribution.

## Appendix B Computing Resources

**Table 1: Hardware specification**
|  | CPU | clock | # threads | GPU |
| --- | --- | --- | --- | --- |
| Alpha | Intel Xeon Gold 6138 | 2.00GHz | 80 | TITAN RTX |
| Beta | Intel Xeon E5-2650 v4 | 2.20GHz | 48 | GeForce RTX 2080ti |
| Gamma | Intel Xeon Gold 6226R | 2.90GHz | 64 | GeForce RTX 3090 |
| Softmax | Intel i9-13900K | 3.00GHz | 32 | GeForce RTX 4090 |

## Appendix C Code snippets for FoMs

Parameters of Code 1, 2

- •
pattern_input: The grating pattern $g$ or $u$
- •
wavelength: The wavelength of light $\lambda$
- •
fourier_order: The Fourier truncation order of RCWA
- •
deflected_angle: The desired deflection angle $\theta$
- •
field_res: The resolution of the field

Please refer to Appendix G for more physical conditions in meent.

Figure: Code 1: Method for caculating electric field $v$

Figure: Code 2: Method for caculating deflection efficiency $\eta$

## Appendix D Neural PDE solver

### D.1 Training details

##### Dataset

We split 10,000 pairs of $(u,v)$ into 8000 training pairs and 2000 test pairs for each of nine physical conditions. An instance of $u$ is sized 1$\times$256$\times$256, each element indicating whether it is filled or empty. $v$ is sized 2$\times$256$\times$256, one channel for real part and another for imaginary part of electric field, each element expressing the intensity of electric field.
The set of nine physical conditions are shown in the first column of Table 5, and the set is followed from [].
Fourier truncation order is set to 40.

For testing zero-shot super-resolution, the structures are transferred to higher resolutions, and corresponding electric fields are calculated with Code 1. Please refer to our Github repository for the script to generate the data.

**Table 2: Common hyperparameters**
| Name | Value |
| --- | --- |
| # of Epochs | 100 |
| Optimizer | AdamW |
| Learning rate | 1E-3 |
| LR scheduler | OneCycleLR |
| Base momentum | 0.85 |
| Max momentum | 0.95 |
| Activation | GELU |

##### Fourier neural operator (FNO)

We used 3,268,062 parameters for training FNO.
To serve as a baseline, we adhered closely to the architecture described in [], except for Tucker factorization.

**Table 3: FNO hyperparameters**
| Name | Value |
| --- | --- |
| # of modes | [24, 24] |
| Lifting channels | 32 |
| Hidden channels | 32 |
| Projection channels | 32 |
| # of layers | 10 |
| Domain padding | 0.015625 |
| Factorization | Tucker |
| Factorization rank | 0.5 |
| Normalization | BatchNorm |

##### UNet

We used vanilla UNet described in the original paper [], of which parameters counts up to 31,036,546.

##### Computational resource

Both FNO and UNet was trained on Beta server of Table 1.
FNO was trained for 3.80 hours, and UNet was trained for 1.86 hours.
Both algorithms used single GPU of Beta server and consumed most of the GPU memory.

##### Remark

We utilized the FNO code(^2^22[https://github.com/neuraloperator/neuraloperator](https://github.com/neuraloperator/neuraloperator)) of the original author, under MIT license.
Also, the widely used UNet implementation(^3^33[https://github.com/milesial/Pytorch-UNet](https://github.com/milesial/Pytorch-UNet)) available under the GPL-3.0 license.

### D.2 Test error

We train FNO on various losses shown in Table 4, and name it as {Model}-{Training loss}, e.g. FNO-L2 is a FNO trained with L2 loss. A model is trained specifically for single physical condition.

On metagrating, more intense and complex interactions occur around the grating area.
What makes this area more important is that, theoretically, the deflection efficiency can be calculated just with the field profile of grating area. Take a look at the supplementary of [].

Therefore, we derive a simple loss coined as region-wise (RW) L2 loss, which puts more weight on the grating area, s.t. $c_{1}+c_{2}=1$. We set $c_{1}=0.7$ and $c_{2}=0.3$.
Lastly, H1 loss is a norm in Sobolev space which integrates the norm of first derivative of the target. Training with H1 loss promotes smoother function [].

**Table 4: Loss functions. $\|\cdot\|$ is the Euclidean norm, RW is shorthand for region-wise. $grating,air,glass$ refers to the sets of indices for each region on matrix representation. All losses are relative error, i.e. normalized by the magnitude of the ground truth $y$.**
| Name | Notation | Definition |
| --- | --- | --- |
| L2 loss | $L_{2}$ | ${\|y-\hat{y}\|}\ /\ {\|y\|}$ |
| RW L2 loss | $L_{2,rw}$ | $c_{1}\cdot L_{2,grating}+c_{2}\cdot(L_{2,air}+L_{2,glass})$ |
| H1 loss | $H_{1}$ | $\sqrt{{\big{(}{\|y-\hat{y}\|}^{2}+{\|y^{\prime}-\hat{y}^{\prime}\|}^{2}\big{)}<br>}\ /\ {({\|y\|}^{2}+{\|y^{\prime}\|}^{2})}}$ |

Table 5 shows that FNO-RW L2 achieves slightly lower error than FNO-L2, but it is not significant enhancement compared to FNO-H1.
FNO-H1 shows best performance across all test metrics, L2, RW L2, and H1.
UNet is trained only with L2 loss, serving as a simple baseline. Comparing mean L2 values, 8.71 of FNO-L2 is 76% lower error than 34.80 of UNet-L2.

**Table 5: Test Error. Of the column names, top row is the name of the models and bottom row is the test metrics.**
|  | UNet-L2 | FNO-L2 | FNO-RW L2 | FNO-H1 |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Condition ($\lambda\$/ $\theta$) | L2$\downarrow$ | RW L2$\downarrow$ | H1$\downarrow$ | L2 | RW L2 | H1 | L2 | RW L2 | H1 | L2 | RW L2 | H1 |
| $1100nm\ /\ 70^{\circ}$ | 34.04 | 22.64 | 33.28 | 7.15 | 6.52 | 14.57 | 7.35 | 4.14 | 10.95 | 6.04 | 3.56 | 6.35 |
| $1100nm\ /\ 60^{\circ}$ | 41.61 | 41.86 | 47.82 | 14.57 | 17.37 | 26.65 | 16.03 | 14.7 | 24.11 | 11.09 | 10.98 | 14.62 |
| $1100nm\ /\ 50^{\circ}$ | 24.37 | 56.33 | 61.05 | 2.52 | 22.38 | 33.70 | 2.58 | 21.93 | 33.81 | 2.07 | 12.33 | 17.37 |
| $1000nm\ /\ 70^{\circ}$ | 43.44 | 22.17 | 29.55 | 15.15 | 5.7 | 12.16 | 15.19 | 4.93 | 11.91 | 9.02 | 3.35 | 5.42 |
| $1000nm\ /\ 60^{\circ}$ | 34.02 | 54.74 | 56.98 | 10.7 | 21.89 | 32.66 | 9.5 | 22.93 | 32.74 | 7.88 | 15.05 | 19.21 |
| $1000nm\ /\ 50^{\circ}$ | 28.46 | 39.62 | 44.28 | 2.88 | 12.34 | 22.51 | 2.25 | 11.66 | 21.50 | 2.19 | 8.26 | 12.15 |
| $900nm\ /\ 70^{\circ}$ | 40.78 | 27.21 | 34.25 | 15.14 | 8.37 | 15.05 | 13.63 | 6.51 | 12.67 | 10.8 | 5.03 | 7.31 |
| $900nm\ /\ 60^{\circ}$ | 31.36 | 30.53 | 34.07 | 6.07 | 11.10 | 17.27 | 5.47 | 9.08 | 14.61 | 4.85 | 7.26 | 9.24 |
| $900nm\ /\ 50^{\circ}$ | 35.11 | 51.64 | 51.59 | 4.23 | 22.87 | 30.79 | 3.77 | 19.89 | 27.33 | 3.29 | 14.91 | 17.75 |
| Mean | 34.80 | 38.53 | 43.65 | 8.71 | 14.28 | 22.81 | 8.42 | 12.86 | 21.07 | 6.36 | 8.97 | 12.16 |
| $\pm$Std | $\pm$5.95 | $\pm$12.81 | $\pm$10.79 | $\pm$5.95 | $\pm$6.58 | $\pm$7.93 | $\pm$5.12 | $\pm$6.92 | $\pm$8.47 | $\pm$3.32 | $\pm$4.31 | $\pm$5.00 |

## Appendix E Metasurface design

### E.1 Training RL agent

Budget refers to the total episode steps consumed for training an agent.

We limit the budget to 50,000 steps, which is 75% less than the budget used in [].

**Table 6: Fixed configurations for all algorithms**
| Parameter | Value |
| --- | --- |
| Budget | 50,000 steps |
| Initial structure | $g_{0}=[1,1,\ldots,1]$ |
| Episode length | $T=512$ |
| Replay buffer size | 25,000 (adjusted proportionally to budget) |
| Asynchronous environments | 8 |
| Number of cells | $k=256$ |
| Wavelength | $\lambda=1100\,\text{nm}$ |
| Desired deflection angle | $\theta=70^{\circ}$ |
| Fourier truncation order | 40 |

##### Dataset

Please refer to out Github repository for RL environment utilizing meent.

##### DQN

We mostly follow the previous work [].
The structure $g_{t}$ is encoded with shallow UNet, and the reward $r_{t}=\eta_{t}-\eta_{t-1}$ is received.
For fair comparison with DreamerV3 L [], following details were changed. Much larger number of parameters (70,711,873) was used, and physics-informed weight initialization was replaced by Pytorch’s default initialization []. Additionally, 1000 steps were used for warmup to fill empty replay buffer.

##### DreamerV3

DreamerV3 was trained with 99,789,440 parameters.
Most of the hyperparameters from original paper [] were reused, excluding: 1,024 steps were used for warmup, replay ratio was increased for sample efficiency, batch size and sequence length was adjusted due to our task’s relatively shorter episode length.

**Table 7: DreamerV3 hyperparameters**
| Name | Value |
| --- | --- |
| Model size | L |
| Replay ratio | 2 |
| Batch size | 8 |
| Sequence length | 32 |

As mentioned in the main text, DreamerV3 agent observes additional feature, the electric field $v_{t}$. The structure $g_{t}$ and electric field $v_{t}$ are encoded by MLP and CNN respectively, and concatenated to form a latent state.
With the input action and latent state, the dynamics model predicts next state, reward and done condition.
Simply put, the dynamics model functions as the environment.

Emprically, DreamerV3 failed the metasurface optimization with the reward of efficiency change $r_{t}=\eta_{t}-\eta_{t-1}$. From the experimental observation, we hypothesized that if the model truly understands the underlying physics, the reward predictor should directly predict efficiency $r_{t}=\eta_{t}$, not the change of efficiency. This hypothesis was inspired by the fact that efficiency can be derived from electric field as mentioned in D.2. With this hypothesized reward engineering, DreamerV3 agent successfully learnt to optimize the metasurface structure.

##### Computational resource

For servers in Table 1, DreamerV3 was trained for 10.88 hours on Softmax server with single GPU. DQN was trained for 1.6 hours on Alpha server.
Both algorithm used its server’s single GPU and consumed most of the GPU memory.

Despite the big difference in training time, when the device is expanded to high dimensionality, the simulation time can occupy the biggest portion of training time []. The main point of our example here is to show the dynamics model’s potential as surrogate solver in decision process, and we leave high dimensional problem as a future work.

##### Parallelization

Example benchmark on the axis of number of workers. The code for adapting RLlib wil be provided on our Github repository. Figure 9 shows that the calculation time sub-linearly decreases as the number of workers increases.

Figure: Figure 9: Parallelization of meent with Ray/RLlib.
Refer to caption: extracted/5658627/sections/applications/images/async-workers.png

##### Remark

The experimental code was adapted from SheepRL(^4^44[https://github.com/Eclectic-Sheep/sheeprl](https://github.com/Eclectic-Sheep/sheeprl)) [] for DreamerV3.
We utilized the previous version of DreamerV3, prior to the updated release on April 17, 2024.
Details of the training procedure and architecture are beyond the scope of this paper. For comprehensive information, please refer to [].

### E.2 High Impact Cell

We spotted the accurate prediction of dynamics model for scarcely happening transition. Introduced in [], a high impact cell refers to a cell that incurs abrupt change in some FoMs when flipped, which is very small change in the material distribution. To artificially create this case, a fully trained agent is run an episode and produces a trajectory. At the step of the trajectory when the efficiency reached about 75%, a high impact cell is manually found by flipping every cell of the structure at the step. With the found index of high impact cell, the action is fed into dynamics model, to predict next electric field.

As shown in Figure 10, dynamics model accurately captures the transition whereas FNO-H1 model, trained under same physical condition in Table 5, entirely fails to predict this phenomena. FNO-H1 might perform better if trained with similar distributions of data, but it is very difficult to draw similar patterns from extremely large design space size, $2^{256}/256$, where the denominator $256$ is for the periodicity.

Figure: Figure 10: High impact cell phenomena. Flipping a single $254^{th}$ cell from silicon to air (red arrow) results in completely different electric field and large decrease in deflection efficiency $\eta$ from about 75% to 17%. Our world model is able to capture the transition. Field intensity is clipped from 0.4 to 0.6 for clearer visualization. All of images are real part of the field.
Refer to caption: extracted/5658627/sections/applications/images/high-impact.png

## Appendix F OCD demonstration

To simulate a real-world scenario where we have the real devices and their spectra, we first determine values for ground truth of the design parameters denoted as P, and generate spectra S with simulation. These values are typically provided from domain experts. Our chosen values are in Table 8.

**Table 8: Design parameter information**
| Parameter | Variable name | Mean | STD | Ground Truth |
| --- | --- | --- | --- | --- |
| P1 | l1_o1_length_x | 100 | 3 | 101.5 |
| P2 | l1_o1_length_y | 80 | 3 | 81.5 |
| P3 | l1_o2_length_x | 100 | 3 | 98.5 |
| P4 | l1_o2_length_y | 80 | 3 | 81.5 |
| P5 | l2_o1_length_x | 30 | 2 | 31 |
| P6 | l2_o2_length_x | 50 | 1 | 49.5 |
| P7 | l1_thickness | 200 | 10 | 205 |
| P8 | l2_thickness | 300 | 10 | 305 |

Figure: Figure 11: Stack in experiment.
Refer to caption: extracted/5658627/sections/appendix/figures/ocd_model.png

Figure 11 depicts the stack utilized in the demonstration. Two layers are stacked on the silicon substrate, each containing objects within. Light is illuminated from the top, and the reflected light is acquired and processed into spectra.

To address this inverse problem of finding design parameters from spectra, initial values for optimization need to be determined. These conditions are also provided by domain experts. In this demonstration, these values are drawn from a normal distribution without correlation. The mean and standard deviation (STD) are listed in Table 8.

**Table 9: Design parameter information**
| Optimizer | Learning Rate | Other conditions |
| --- | --- | --- |
| Momentum | 1E2 | momentum: 0.9 |
| Adagrad | 1E0 |  |
| RMSProp | 1E-1 |  |
| Adam | 1E-1 |  |
| RAdam | 1E0 |  |

The hyperparameters utilized for the optimization demonstration are presented in Table 9. Default values from PyTorch are used for conditions not explicitly mentioned.
The learning rate was determined through a concise parameter-sweep test, which assessed three different values of the learning rate for each optimizer, as presented in Figure 12.

Figure: Figure 12: Hyperparameter sweep.
Refer to caption: extracted/5658627/sections/appendix/figures/ocd_hyperparameter_sweep.png

## Appendix G Background theory

RCWA is the sequence of the following processes: solving the Maxwell’s equations, finding the eigenmodes of a layer and connecting these layers including the superstrate and substrate to calculate the diffraction efficiencies.
Precisely, the electromagnetic field and permittivity geometry are transformed from the real space to the Fourier space (also called the reciprocal space or k-space) by Fourier analysis. Maxwell’s equations are then solved per layer through convolution operation, and a general solution of the field in each direction can be obtained. This general solution can be represented in terms of eigenmodes (eigenvectors) and eigenvalues with eigendecomposition, and used to calculate diffraction efficiencies by applying boundary conditions and connecting to adjacent layers.

This chapter provides a comprehensive explanation of the theories, formulations and implementations of meent in the following sections:

- 1.
Structure design: the device geometry is defined and modeled within meent framework.
- 2.
Fourier analysis of geometry: the device geometry is transformed into the Fourier space, allowing the decomposition of the structure into its corresponding spatial frequency components.
- 3.
Eigenmodes identification: RCWA identifies the eigenmodes that present within each layer of the periodic structure. These eigenmodes represent the possible electromagnetic field solutions that can exist within the system.
- 4.
Connecting layers: Rayleigh coefficients and diffraction efficiencies are determined using the transfer matrix method by connecting the layers. This step enables the determination of the overall electromagnetic response of the entire system.
- 5.
Enhanced transmittance matrix method: the implementation technique that avoids the inversion of some matrices which are possibly ill-conditioned.
- 6.
Topological derivative vs Shape derivative: two types of derivatives that meent supports are explained.

### G.1 Structure design

Figure: (a) raster modeling
Refer to caption: extracted/5658627/sections/manual/images/rcwa/rastermodelling.png

meent supports two distinct types of geometry modeling: the raster modeling and the vector modeling.
In the raster modeling, the device geometry is gridded and filled with the refractive index of the corresponding material as in Figure 13(a).
This approach is advantageous for solving optimization problems related to freeform metasurfaces.
The vector modeling (shown in Figure 13(b)), on the other hand, represents the geometry as an union of primitive shapes and each primitive shape is defined by edges and vertices like vector-type image.
Consequently, it is memory-efficient and has less parameters to optimize by not keeping the whole array as the raster-type does.
This feature is especially valuable in OCD metrology where semiconductor device comprises highly complex structures. raster-type methods may become impractical in such scenarios due to the limitations of grid-based representations. One of the key advantages provided by vector modeling is that the minimum feature size is not constrained by the grid size. This flexibility allows for more accurate and detailed representation of complex structures, making vector modeling essential for accurate simulation.

### G.2 Fourier analysis of geometry

In RCWA, the device geometry needs to be mapped to the Fourier space using Fourier analysis. To achieve this, the device is sliced into multi-layers so that each layer has Z-invariant (layer stacking direction) permittivity distribution. In other words, the permittivity can be considered as a piecewise-constant function that varies in X and Y but not Z direction in each layer. Then the geometry in real space can be expressed as a weighted sum of Fourier basis:

$$ $\displaystyle\varepsilon(x,y)=\sum_{m=-\infty}^{\infty}\sum_{n=-\infty}^{ \infty}c_{n,m}\cdot\exp{\left[j\cdot 2\pi\left(\frac{x}{\Lambda_{x}}m+\frac{y} {\Lambda_{y}}n\right)\right]},$ (6) $$

where $\Lambda_{x},\Lambda_{y}$ are the period of the unit cell and $c_{n,m}$ is the Fourier coefficients ($m^{th}$ in X and $n^{th}$ in Y).
However, due to the limitation of the digital computations, this has to be approximated with truncation:

$$ $\displaystyle\varepsilon(x,y)\simeq\sum_{m=-M}^{M}\sum_{n=-N}^{N}c_{n,m}\cdot \exp{\left[j\cdot 2\pi\left(\frac{x}{\Lambda_{x}}m+\frac{y}{\Lambda_{y}}n \right)\right]},$ (7) $$

where $M,N$ are the Fourier Truncation Order (FTO, the number of harmonics in use) in the X and Y direction, and these can be considered as hyperparamters that affects the simulation accuracy.

Here, $c_{n,m}$ is the permittivity distribution in the Fourier space which is our interest and can be found by one of these two methods:
Discrete Fourier Series (DFS) or Continuous Fourier Series (CFS). To be clear, CFS is Fourier series on piecewise-constant function (permittivity distribution in our case). This name was given to emphasize the characteristics of each type by using opposing words. The output array of DFS and CFS have the same shape and can be substituted for each other.

In DFS, the function $\varepsilon(x,y)$ to be transformed is sampled at a finite number of points, and this means it’s given in matrix form with rows and columns, $\varepsilon_{\mathtt{r},\mathtt{c}}$. The coefficients of DFS are then given by this equation:

$$ $c_{n,m}=\frac{1}{P_{x}P_{y}}\sum_{\mathtt{c}=0}^{P_{x}-1}\sum_{\mathtt{r}=0}^{ P_{y}-1}\varepsilon_{\mathtt{r},\mathtt{c}}\cdot\exp{\left[-j\cdot 2\pi\left( \frac{m}{P_{x}}\mathtt{c}+\frac{n}{P_{y}}\mathtt{r}\right)\right]},$ (8) $$

where $P_{x},P_{y}$ are the sampling frequency (the size of the array), $\varepsilon_{\mathtt{r},\mathtt{c}}$ is the ${(\mathtt{r},\mathtt{c})}^{th}$ element of the permittivity array.

There is an essential but easily overlooked fact: the sampling frequency ($P_{x},P_{y}$) is very important in DFS [, , ]. If this is not enough, an aliasing occurs: DFS cannot correctly capture the original signal (you can easily see the wheels of a running car in movies rotating in the opposite direction; this is also an aliasing and called the wagon-wheel effect).
In RCWA, this may occur during the process of sampling the permittivity distribution.
To resolve this, meent provides a scaling function by default - that is simply to increase the size of the permittivity array by repeatedly replicating the elements while keeping the original shape of the pattern.
This option improves the representation of the geometry in the Fourier space and results in more accurate RCWA simulations.

CFS utilizes the entire function to find the coefficients while DFS uses only some of them.
This means that CFS prevents potential information loss coming from the intrinsic nature of DFS, thereby enables more accurate simulation.
The Fourier coefficients can be expressed as follow:

$$ $\displaystyle c_{n,m}=\frac{1}{\Lambda_{x}\Lambda_{y}}\int_{x_{0}}^{x_{0}+ \Lambda_{x}}\int_{y_{0}}^{y_{0}+\Lambda_{y}}\varepsilon(x,y)\cdot\exp{\left[-j \cdot 2\pi\left(\frac{m}{\Lambda_{x}}x+\frac{n}{\Lambda_{y}}y\right)\right]}dydx.$ (9) $$

The information that CFS needs are the points of discontinuity and the permittivity value in each area sectioned by those points, whereas DFS needs the whole permittivity array as in Figure 13.

DFS and CFS have its own advantages and one can be chosen according to the purpose of the simulation. Basically, DFS is proper for Raster modeling since its operations are mainly on the pixels (array) and the input of the Raster modeling is the array. This is naturally connected to the pixel-wise operation (cell flipping) in metasurface freeform design.
CFS is suitable for Vector modeling because it deals with the graph (discontinuous points and length) of the objects and Vector modeling takes that graph as an input. Hence it enables direct and precise optimization of the design parameters (such as the width of a rectangle) without grid that severely limits the resolution. We will address this in section G.6

### G.3 Eigenmodes identification

Once the permittivity distribution is mapped to the Fourier space, the next step is to apply Maxwell’s equations to identify the eigenmodes of each layer.
In this section, we extend the mathematical formulation of the 1D conical incidence case described in [] to the 2D grating case as illustrated in Figure 14.
To ensure the consistency and clarity, we adopt the same notations and the sign convention of $(+jwt)$.

Figure: Figure 14: Geometry for the stack. Two-dimensional grating layers and incident ray.
Refer to caption: extracted/5658627/sections/manual/images/theory/geometry.png

We consider the normalized excitation wave at the superstrate to take the following form:

$$ $\displaystyle\mathbf{E}_{inc}=\mathbf{u}\cdot e^{-jk_{0}\mathtt{n}_{\text{I}}( \sin{\theta}\cdot\cos{\phi}\cdot x+\sin{\theta}\cdot\sin{\phi}\cdot y+\cos{ \theta}\cdot z)},$ (10) $$

where $\mathbf{u}$ is the normalized amplitudes of the wave in each direction:

$$ (11) $$

and $k_{0}=2\pi/\lambda_{0}$ with $\lambda_{0}$ the wavelength of the light in free space, $\mathtt{n}_{\text{I}}$ is the refractive index of the superstrate, $\theta$ is the angle of incidence, $\phi$ is the rotation (azimuth) angle and $\psi$ is the angle between the electric field vector and the plane of incidence.

The electric fields in the superstrate and substrate (we will designate these layers by I and II as in []) can be expressed as a sum of incident, reflected and transmitted waves as the Rayleigh expansion [, , ]:

$$ $\displaystyle\mathbf{E}_{\text{I}}$ $\displaystyle=\mathbf{E}_{inc}+\sum_{m=-M}^{M}\sum_{n=-N}^{N}\mathbf{R}_{n,m}e ^{-j(k_{x,m}x+k_{y,n}y-k_{\text{I},z,(n,m)}z)},$ (12) $\displaystyle\mathbf{E}_{\text{II}}$ $\displaystyle=\sum_{m=-M}^{M}\sum_{n=-N}^{N}\mathbf{T}_{n,m}e^{-j\{k_{x,m}x+k_ {y,n}y+k_{\text{II},z,(n,m)}(z-d)\}},$ (13) $$

where $M$ and $N$ are the Fourier Truncation Order (FTO) which is related to the number of harmonics in use, and the in-plane components of the wavevector ($k_{x,m}$ and $k_{y,n}$) are determined by the Bloch’s theorem (this has many names and one of them is Floquet condition) [, ],

$$ $\displaystyle k_{x,m}$ $\displaystyle=k_{0}\Big{(}\mathtt{n}_{\text{I}}\sin{\theta}\cos{\phi}-m\frac{ \lambda_{0}}{\Lambda_{x}}\Big{)},$ (14) $\displaystyle k_{y,n}$ $\displaystyle=k_{0}\Big{(}\mathtt{n}_{\text{I}}\sin{\theta}\sin{\phi}-n\frac{ \lambda_{0}}{\Lambda_{y}}\Big{)},$ (15) $$

where $\Lambda_{x}$ and $\Lambda_{y}$ are the period of the unit cell, and the out-of-plane wavevector is determined from the dispersion relation:

$$ $\displaystyle k_{\ell,z,(n,m)}=\begin{cases}+\ [(k_{0}\mathtt{n}_{\ell})^{2}-{ k_{x,m}}^{2}-{k_{y,n}}^{2}]^{1/2}&,\quad\text{if}\quad({k_{x,m}}^{2}+{k_{y,n}} ^{2})<(k_{0}\mathtt{n}_{\ell})^{2}\\ -j[{k_{x,m}}^{2}+{k_{y,n}}^{2}-(k_{0}\mathtt{n}_{\ell})^{2}]^{1/2}&,\quad\text {if}\quad({k_{x,m}}^{2}+{k_{y,n}}^{2})>(k_{0}\mathtt{n}_{\ell})^{2}\end{cases} ,\quad\ell=\text{I},\text{II}.$ (16) $$

Here, $k_{\ell,z,(n,m)}$ can be categorized into propagation mode and evanescent mode depending on whether it’s real or imaginary. $\mathbf{R}_{n,m}\text{ and }\mathbf{T}_{n,m}$ are the Rayleigh coefficients (also called the reflection and transmission coefficients): $\mathbf{R}_{n,m}$ is the normalized (3-dimensional) vector of electric field amplitude which is the ($m^{th}$ in X and $n^{th}$ in Y) mode of reflected waves in the superstrate and $\mathbf{T}_{n,m}$ is the normalized (3-dimensional) vector of electric field amplitude which is the ($m^{th}$ in X and $n^{th}$ in Y) mode of transmitted waves in the substrate.

Inside the grating layer, the electromagnetic field can be expressed as a superposition of plane waves by the Bloch’s theorem:

$$ $\displaystyle\mathbf{E}_{g}(x,y,z)$ $\displaystyle=\sum_{m=-M}^{M}\sum_{n=-N}^{N}\boldsymbol{\mathfrak{S}}_{g,(n,m) }\cdot e^{-j(k_{x,m}x+k_{y,n}y+k_{g,z}z)},$ (17) $\displaystyle\mathbf{H}_{g}(x,y,z)$ $\displaystyle=\sum_{m=-M}^{M}\sum_{n=-N}^{N}\boldsymbol{\mathfrak{U}}_{g,(n,m) }\cdot e^{-j(k_{x,m}x+k_{y,n}y+k_{g,z}z)},$ (18) $$

where $k_{g,z}$ is the wavevector in Z-direction (this is unique per layer hence the notation g was kept to distinguish) and $\boldsymbol{\mathfrak{S}}_{g,(n,m)}$ and $\boldsymbol{\mathfrak{U}}_{g,(n,m)}$ are the vectors of amplitudes in each direction at $(m,n)^{th}$ order:

$$ $\displaystyle\boldsymbol{\mathfrak{S}}_{g,(n,m)}$ $\displaystyle=\mathfrak{S}_{g,(n,m),x}\ \hat{x}+\mathfrak{S}_{g,(n,m),y}\ \hat {y}+\mathfrak{S}_{g,z}\ \hat{z},$ (19) $\displaystyle\boldsymbol{\mathfrak{U}}_{g,(n,m)}$ $\displaystyle=\mathfrak{U}_{g,(n,m),x}\ \hat{x}+\mathfrak{U}_{g,(n,m),y}\ \hat {y}+\mathfrak{U}_{g,z}\ \hat{z}.$ (20) $$

It is also possible to detach wavevector term on $z$ from exponent and combine with $\mathbf{\mathfrak{S}}_{g,(n,m)}$ and $\mathbf{\mathfrak{U}}_{g,(n,m)}$ in Equations 17 and 18 to make $\mathbf{S}_{g,(n,m)}(z)$ and $\mathbf{U}_{g,(n,m)}(z)$ which are dependent on $z$ as shown below:

$$ $\displaystyle\mathbf{S}_{g,(n,m)}(z)=\boldsymbol{\mathfrak{S}}_{g,(n,m)}\cdot e ^{-jk_{g,z}z},$ (21) $\displaystyle\mathbf{U}_{g,(n,m)}(z)=\boldsymbol{\mathfrak{U}}_{g,(n,m)}\cdot e ^{-jk_{g,z}z},$ (22) $$

then Equations 17 and 18 become

$$ $\displaystyle\mathbf{E}_{g}(x,y,z)$ $\displaystyle=\sum_{m=-M}^{M}\sum_{n=-N}^{N}\mathbf{S}_{g,(n,m)}(z)\cdot e^{-j (k_{x,m}x+k_{y,n}y)},$ (23) $\displaystyle\mathbf{H}_{g}(x,y,z)$ $\displaystyle=\sum_{m=-M}^{M}\sum_{n=-N}^{N}\mathbf{U}_{g,(n,m)}(z)\cdot e^{-j (k_{x,m}x+k_{y,n}y)}.$ (24) $$

Equations 17 and 18 are used in [, , ] and Equations 23 and 24 in [, ]. Whichever is used, the result is the same: we will show the development using ($\boldsymbol{\mathfrak{S}}_{g,(n,m)}$, $\boldsymbol{\mathfrak{U}}_{g,(n,m)}$) with the eigendecomposition and then come back to ($\mathbf{S}_{g,(n,m)}(z)$ and $\mathbf{U}_{g,(n,m)}(z)$) with the partial differential equations.

The behavior of the electromagnetic fields can be described by the formulae, called the Maxwell’s equations. Among them, we will use the third and fourth equations,

$$ $\displaystyle\nabla\times\mathbf{E}$ $\displaystyle=-j\omega\mu_{0}\mathbf{H},$ (25) $\displaystyle\nabla\times\mathbf{H}$ $\displaystyle=j\omega\varepsilon_{0}\varepsilon_{r}\mathbf{E},$ (26) $$

to find the electric and magnetic field inside the grating layer - $\mathbf{E}_{g}$ and $\mathbf{H}_{g}$.
Since RCWA is a technique that solves Maxwell’s equations in the Fourier space, curl operator in real space becomes multiplication and multiplication in real space becomes the convolution operator. For this convolution operation, the full set of the modes of the fields and the geometry are required so we introduce a vector notation in the subscript to denote it’s a vector with all the harmonics in use, i.e.,

$$ (27) $$

where $\boldsymbol{F}\in\{S,U,\mathfrak{S},\mathfrak{U}\}$ and $r\in\{x,y,z\}$. Some variables will be scaled by some factors:

$$ $\displaystyle\tilde{\mathbf{H}}_{g}=-j\sqrt{\varepsilon_{0}/{\mu_{0}}}\mathbf{ H}_{g},\quad\tilde{k}_{x}=k_{x}/k_{0},\quad\tilde{k}_{y}=k_{y}/k_{0},\quad \tilde{k}_{g,z}=k_{g,z}/k_{0},\quad\tilde{z}=k_{0}z.$ (28) $$

Substituting Equations 17 and 18 ($\mathbf{E}_{g}$ and $\tilde{\mathbf{H}}_{g}$ with $\boldsymbol{\mathfrak{S}}_{g}$ and $\boldsymbol{\mathfrak{U}}_{g}$) into Equations 25 and 26 (Maxwell’s equations) and eliminating Z-directional components ($\mathbf{E}_{g,z}$ and $\tilde{\mathbf{H}}_{g,z}$) derive the matrix form of the Maxwell’s equations composed of in-plane components $(\hat{x},\hat{y})$ in the Fourier space:

$$ $\displaystyle(-j\tilde{k}_{g,z})\begin{bmatrix}\boldsymbol{\mathfrak{S}}_{g, \vec{x}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\\ \boldsymbol{\mathfrak{S}}_{g,\vec{y}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\end{ bmatrix}=\boldsymbol{{\Omega}}_{g,L}\begin{bmatrix}\boldsymbol{\mathfrak{U}}_{ g,\vec{x}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\\ \boldsymbol{\mathfrak{U}}_{g,\vec{y}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\end{bmatrix}$ (29) $$

$$ $\displaystyle(-j\tilde{k}_{g,z})\begin{bmatrix}\boldsymbol{\mathfrak{U}}_{g, \vec{x}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\\ \boldsymbol{\mathfrak{U}}_{g,\vec{y}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\end{ bmatrix}=\boldsymbol{{\Omega}}_{g,R}\begin{bmatrix}\boldsymbol{\mathfrak{S}}_{ g,\vec{x}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\\ \boldsymbol{\mathfrak{S}}_{g,\vec{y}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\end{bmatrix}$ (30) $$

$$ $\displaystyle(-j\tilde{k}_{g,z})^{2}\begin{bmatrix}\boldsymbol{\mathfrak{S}}_{ g,\vec{x}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\\ \boldsymbol{\mathfrak{S}}_{g,\vec{y}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\end{ bmatrix}=\boldsymbol{{\Omega}}_{g,LR}^{2}\begin{bmatrix}\boldsymbol{\mathfrak{ S}}_{g,\vec{x}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\\ \boldsymbol{\mathfrak{S}}_{g,\vec{y}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\end{bmatrix}$ (31) $$

where

$$ $\displaystyle\boldsymbol{{\Omega}}_{g,L}=\begin{bmatrix}(-\tilde{\mathbf{K}}_{ x}\llbracket\varepsilon_{r,g}\rrbracket^{-1}\tilde{\mathbf{K}}_{y})&(\tilde{ \mathbf{K}}_{x}\llbracket\varepsilon_{r,g}\rrbracket^{-1}\tilde{\mathbf{K}}_{x }-\mathbf{I})\\ (\mathbf{I}-\tilde{\mathbf{K}}_{y}\llbracket\varepsilon_{r,g}\rrbracket^{-1} \tilde{\mathbf{K}}_{y})&(\tilde{\mathbf{K}}_{y}\llbracket\varepsilon_{r,g} \rrbracket^{-1}\tilde{\mathbf{K}}_{x})\end{bmatrix},$ (32) $$

$$ $\displaystyle\boldsymbol{{\Omega}}_{g,R}=\begin{bmatrix}(-\tilde{\mathbf{K}}_{ x}\tilde{\mathbf{K}}_{y})&(\tilde{\mathbf{K}}_{x}^{2}-\llbracket\varepsilon_{r ,g}\rrbracket)\\ (\llbracket\varepsilon_{r,g}^{-1}\rrbracket^{-1}-\tilde{\mathbf{K}}_{y}^{2})&( \tilde{\mathbf{K}}_{y}\tilde{\mathbf{K}}_{x})\end{bmatrix},$ (33) $$

$$ $\displaystyle\boldsymbol{{\Omega}}_{g,LR}^{2}=\begin{bmatrix}{\tilde{\mathbf{K }}_{y}}^{2}+(\tilde{\mathbf{K}}_{x}\llbracket\varepsilon_{r,g}\rrbracket^{-1} \tilde{\mathbf{K}}_{x}-\mathbf{I})\llbracket\varepsilon_{r,g}^{-1}\rrbracket^{ -1}&\tilde{\mathbf{K}}_{x}(\llbracket\varepsilon_{r,g}\rrbracket^{-1}\tilde{ \mathbf{K}}_{y}\llbracket\varepsilon_{r,g}\rrbracket-\tilde{\mathbf{K}}_{y})\\ \tilde{\mathbf{K}}_{y}(\llbracket\varepsilon_{r,g}\rrbracket^{-1}\tilde{ \mathbf{K}}_{x}\llbracket\varepsilon_{r,g}^{-1}\rrbracket^{-1}-\tilde{\mathbf{ K}}_{x})&{\tilde{\mathbf{K}}_{x}}^{2}+(\tilde{\mathbf{K}}_{y}\llbracket \varepsilon_{r,g}\rrbracket^{-1}\tilde{\mathbf{K}}_{y}-\mathbf{I})\llbracket \varepsilon_{r,g}\rrbracket\end{bmatrix},$ (34) $$

and

$$ $\displaystyle\tilde{\mathbf{K}}_{r}=\begin{bmatrix}\tilde{k}_{r,(-N,-M)}&0& \cdots&0\\ 0&\tilde{k}_{r,(-N,-M+1)}&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\tilde{k}_{r,(N,M)}\end{bmatrix},\quad r\in\{x,y\},$ (35) $$

and $\llbracket\ \rrbracket$ is the convolution (a.k.a Toeplitz) matrix: $\llbracket\varepsilon_{r,g}\rrbracket$ and $\llbracket\varepsilon_{r,g}^{-1}\rrbracket^{-1}$ are convolution matrices composed of Fourier coefficients of permittivity and one-over-permittivity (by the inverse rule presented in [] and []).

Equation 31 is a typical form of the eigendecomposition of a matrix. The vector [$\boldsymbol{\mathfrak{S}}_{g,\vec{x}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\quad
\boldsymbol{\mathfrak{S}}_{g,\vec{y}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}]^{T}$ is an eigenvector of $\boldsymbol{{\Omega}}_{g,LR}^{2}$ and $j\tilde{k}_{g,z}$ is the positive square root of the eigenvalues. This intuitively shows how the eigenvalues are connected to the Z-directional wavevectors.

It is also possible to use $\mathbf{S}_{g,\vec{x}}(\tilde{z})$ and $\mathbf{S}_{g,\vec{y}}(\tilde{z})$ instead of $\boldsymbol{\mathfrak{S}}_{g,\vec{x}}$ and $\boldsymbol{\mathfrak{U}}_{g,\vec{x}}$ because they satisfy the following relations:

$$ $\displaystyle\frac{\partial^{2}}{\partial(\tilde{z})^{2}}\begin{bmatrix} \mathbf{S}_{g,\vec{x}}(\tilde{z})\\ \mathbf{S}_{g,\vec{y}}(\tilde{z})\end{bmatrix}=\frac{\partial^{2}}{\partial( \tilde{z})^{2}}\begin{bmatrix}\boldsymbol{\mathfrak{S}}_{g,\vec{x}}\cdot e^{-j \tilde{k}_{g,z}\tilde{z}}\\ \boldsymbol{\mathfrak{S}}_{g,\vec{y}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\end{ bmatrix}=(-j\tilde{k}_{g,z})^{2}\begin{bmatrix}\boldsymbol{\mathfrak{S}}_{g, \vec{x}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\\ \boldsymbol{\mathfrak{S}}_{g,\vec{y}}\cdot e^{-j\tilde{k}_{g,z}\tilde{z}}\end{ bmatrix}.$ (36) $$

Hence it is just a matter of choice and we will use PDE form ($\mathbf{S}_{g}$ and $\mathbf{U}_{g}$) for the seamless connection to the 1D conical case in the previous work []. Then Equations 29, 30 and 31 become

$$ $\displaystyle\frac{\partial}{\partial(\tilde{z})}\begin{bmatrix}\mathbf{S}_{g, \vec{x}}(\tilde{z})\\ \mathbf{S}_{g,\vec{y}}(\tilde{z})\end{bmatrix}=\boldsymbol{{\Omega}}_{g,L} \begin{bmatrix}\mathbf{U}_{g,\vec{x}}(\tilde{z})\\ \mathbf{U}_{g,\vec{y}}(\tilde{z})\end{bmatrix},$ (37) $$

$$ $\displaystyle\frac{\partial}{\partial(\tilde{z})}\begin{bmatrix}\mathbf{U}_{g, \vec{x}}(\tilde{z})\\ \mathbf{U}_{g,\vec{y}}(\tilde{z})\end{bmatrix}=\boldsymbol{{\Omega}}_{g,R} \begin{bmatrix}\mathbf{S}_{g,\vec{x}}(\tilde{z})\\ \mathbf{S}_{g,\vec{y}}(\tilde{z})\end{bmatrix},$ (38) $$

$$ $\displaystyle\frac{\partial^{2}}{\partial(\tilde{z})^{2}}\begin{bmatrix} \mathbf{S}_{g,\vec{x}}(\tilde{z})\\ \mathbf{S}_{g,\vec{y}}(\tilde{z})\end{bmatrix}=\boldsymbol{{\Omega}}_{g,LR}^{2 }\begin{bmatrix}\mathbf{S}_{g,\vec{x}}(\tilde{z})\\ \mathbf{S}_{g,\vec{y}}(\tilde{z})\end{bmatrix},$ (39) $$

where Equation (39) is the second order matrix differential equation which has the general solution of the following form

$$ $\displaystyle\begin{bmatrix}\mathbf{S}_{g,\vec{x}}(\tilde{z})\\ \mathbf{S}_{g,\vec{y}}(\tilde{z})\end{bmatrix}=\boldsymbol{w}_{g,1}$ $\displaystyle(c_{g,1}^{+}e^{-q_{g,1}\tilde{z}}+c_{g,1}^{-}e^{+q_{g,1}\tilde{z} })+\cdots+\boldsymbol{w}_{g,\xi}(c_{g,\xi}^{+}e^{-q_{g,\xi}\tilde{z}}+c_{g,\xi }^{-}e^{+q_{g,\xi}\tilde{z}})$ (40) $\displaystyle=\sum_{i=1}^{\xi}\boldsymbol{w}_{g,i}(c_{g,i}^{+}e^{-q_{g,i} \tilde{z}}+c_{g,i}^{-}e^{+q_{g,i}\tilde{z}}),$ (41) $$

where $\xi=(2M+1)(2N+1)$, the total number of harmonics, and $\boldsymbol{w}_{g}$ is the eigenvector, $q_{g}$ is the positive square root of the corresponding eigenvalue ($j\tilde{k}_{g,z}$) and $c_{g}^{\pm}$ are the coefficients (amplitudes) of the mode in each propagating direction (+Z and -Z direction). This can be written in matrix form

$$ $\displaystyle\begin{bmatrix}\mathbf{S}_{g,\vec{x}}(\tilde{z})\\ \mathbf{S}_{g,\vec{y}}(\tilde{z})\end{bmatrix}$ $\displaystyle=\mathbf{W}_{g}\mathbf{Q}_{g}^{-}\mathbf{c}_{g}^{+}+\mathbf{W}_{g }\mathbf{Q}_{g}^{+}\mathbf{c}_{g}^{-}$ (42) $\displaystyle=\mathbf{W}_{g}\begin{bmatrix}\mathbf{Q}_{g}^{-}&\mathbf{Q}_{g}^{ +}\\ \end{bmatrix}\begin{bmatrix}{\mathbf{c}}_{g}^{+}\\ {\mathbf{c}}_{g}^{-}\\ \end{bmatrix},$ (43) $\displaystyle=\begin{bmatrix}\mathbf{W}_{g,11}&\mathbf{W}_{g,12}\\ \mathbf{W}_{g,21}&\mathbf{W}_{g,22}\end{bmatrix}\begin{bmatrix}{\mathbf{Q}_{g, 1}^{-}}&0&{\mathbf{Q}_{g,1}^{+}}&0\\ 0&{\mathbf{Q}_{g,2}^{-}}&0&{\mathbf{Q}_{g,2}^{+}}\end{bmatrix}\begin{bmatrix} \mathbf{c}_{g,1}^{+}\\ \mathbf{c}_{g,2}^{+}\\ \mathbf{c}_{g,1}^{-}\\ \mathbf{c}_{g,2}^{-}\end{bmatrix},$ (44) $$

where $\mathbf{Q}_{g}^{\pm}$ are the diagonal matrices with the exponential of eigenvalues

$$ $\displaystyle\mathbf{Q}_{g}^{\pm}=\begin{bmatrix}e^{\pm{q}_{g,1}}&&0\\ &\ddots&\\ 0&&e^{{\pm q_{g,\xi}}}\end{bmatrix},$ (45) $$

and $\mathbf{W}_{g}$ is the matrix that has the eigenvectors in columns and $\mathbf{c}_{g}^{\pm}$ are the vectors of the coefficients.

Now we can find the general solution of the magnetic field that shares same $\mathbf{Q}_{g}$ and $\mathbf{c}_{g}^{\pm}$ with the electric field in corresponding mode. It can be written in a similar form of Equation 42 as

$$ $\displaystyle\begin{bmatrix}\mathbf{U}_{g,\vec{x}}(\tilde{z})\\ \mathbf{U}_{g,\vec{y}}(\tilde{z})\end{bmatrix}$ $\displaystyle=-\mathbf{V}_{g}\mathbf{Q}_{g}^{-}\mathbf{c}_{g}^{+}+\mathbf{V}_{ g}\mathbf{Q}_{g}^{+}\mathbf{c}_{g}^{-}.$ (46) $$

The negative sign in the first term was given to adjust the direction of the curl operation, $E\times H$, to be in accordance with the wave propagation direction, $\tilde{k}_{g,z}$.
By substituting Equations 42 and 46 into Equation 38, we can get

$$ $\displaystyle\mathbf{V}_{g}=\boldsymbol{{\Omega}}_{g,R}\mathbf{W}_{g}\mathbf{q }_{g}^{-1},$ (47) $$

where $\mathbf{q}_{g}$ is the diagonal matrix with the eigenvalues. This can be written in matrix form

$$ $\begin{split}\mathbf{V}_{g}=\begin{bmatrix}\mathbf{V}_{g,11}&\mathbf{V}_{g,12} \\ \mathbf{V}_{g,21}&\mathbf{V}_{g,22}\end{bmatrix}=\begin{bmatrix}-\tilde{ \mathbf{K}}_{x}\tilde{\mathbf{K}}_{y}&\tilde{{\mathbf{K}}}_{x}^{2}-\llbracket \varepsilon_{r,g}\rrbracket\\ \llbracket\varepsilon_{r,g}^{-1}\rrbracket^{-1}-\tilde{{\mathbf{K}}}_{y}^{2}& \tilde{\mathbf{K}}_{y}\tilde{\mathbf{K}}_{x}\end{bmatrix}\begin{bmatrix} \mathbf{W}_{g,11}&\mathbf{W}_{g,12}\\ \mathbf{W}_{g,21}&\mathbf{W}_{g,22}\end{bmatrix}\begin{bmatrix}\mathbf{q}_{g,1 }&0\\ 0&\mathbf{q}_{g,2}\end{bmatrix}^{-1}\end{split}.$ (48) $$

### G.4 Connecting layers

Once the eigenmodes of each grating layer are identified, the transfer matrix method (TMM) can be utilized to determine the Rayleigh coefficients ($\mathbf{R}_{s},\mathbf{R}_{p},\mathbf{T}_{s},\mathbf{T}_{p}$) and the diffraction efficiencies.
TMM effectively represents this process as a matrix multiplication, where the transfer matrix is constructed by considering the interaction between the eigenmodes of neighboring layers.
This matrix accounts for the energy transfer and phase shift between the eigenmodes, and it is used to propagate the electromagnetic fields through the entire periodic structure.

From the boundary conditions, the systems of equations consisting of the in-plane (tangential) field components $(\mathbf{E}_{\mathtt{s}},\mathbf{E}_{\mathtt{p}},\mathbf{H}_{\mathtt{s}},
\mathbf{H}_{\mathtt{p}})$ can be described at each layer interface. We will first consider the case of a single grating layer cladded with the superstrate and substrate, then expand to multilayer structure. At the input boundary ($z=0$):

$$ (49) $$

and at the output boundary ($z=d$):

$$ $\displaystyle\begin{bmatrix}{\mathbf{W}_{g,ss}\mathbf{X}_{g,1}}&{\mathbf{W}_{g ,sp}\mathbf{X}_{g,2}}&{\mathbf{W}_{g,ss}}&{\mathbf{W}_{g,sp}}\\ {\mathbf{W}_{g,ps}\mathbf{X}_{g,1}}&{\mathbf{W}_{g,pp}\mathbf{X}_{g,2}}&{ \mathbf{W}_{g,ps}}&{\mathbf{W}_{g,pp}}\\ {\mathbf{V}_{g,ss}\mathbf{X}_{g,1}}&{\mathbf{V}_{g,sp}\mathbf{X}_{g,2}}&{- \mathbf{V}_{g,ss}}&{-\mathbf{V}_{g,sp}}\\ {\mathbf{V}_{g,ps}\mathbf{X}_{g,1}}&{\mathbf{V}_{g,pp}\mathbf{X}_{g,2}}&{- \mathbf{V}_{g,ps}}&{-\mathbf{V}_{g,pp}}\end{bmatrix}\begin{bmatrix}\mathbf{c}_ {g,1}^{+}\\ \mathbf{c}_{g,2}^{+}\\ \mathbf{c}_{g,1}^{-}\\ \mathbf{c}_{g,2}^{-}\\ \end{bmatrix}=\begin{bmatrix}\mathbf{I}&\mathbf{0}\\ \mathbf{0}&j\mathbf{Z}_{\text{II}}\\ j\mathbf{Y}_{\text{II}}&\mathbf{0}\\ \mathbf{0}&\mathbf{I}\end{bmatrix}\begin{bmatrix}\mathbf{T}_{s}\\ \mathbf{T}_{p}\end{bmatrix},$ (50) $$

where $\boldsymbol{\delta}_{00}$ is the Kronecker delta function that has 1 at the $(0,0)^{th}$ order and 0 elsewhere.

Here, the variables used above are defined: $\mathbf{X}_{g,1},\mathbf{X}_{g,2}$ are the diagonal matrices

$$ $\displaystyle\mathbf{X}_{g,1}=\begin{bmatrix}e^{-k_{0}q_{g,1,1}d_{g}}&&0\\ &\ddots&\\ 0&&e^{-k_{0}q_{g,1,\xi}d_{g}}\end{bmatrix},\quad\mathbf{X}_{g,2}=\begin{ bmatrix}e^{-k_{0}q_{g,2,1}d_{g}}&&0\\ &\ddots&\\ 0&&e^{-k_{0}q_{g,2,\xi}d_{g}}\end{bmatrix},$ (51) $$

where $d_{g}$ is the thickness of the grating layer, and
$\mathbf{Y}_{\text{I}}$ and $\mathbf{Z}_{\text{I}}$ are

$$ $\displaystyle\mathbf{Y}_{\text{I}}=\begin{bmatrix}\tilde{k}_{\text{I},z,(-N,-M )}&&0\\ &\ddots&\\ 0&&\tilde{k}_{\text{I},z,(N,M)}\end{bmatrix},\quad\mathbf{Z}_{\text{I}}=\frac{ 1}{(\mathtt{n}_{\text{I}})^{2}}\begin{bmatrix}\tilde{k}_{\text{I},z,(-N,-M)}&& 0\\ &\ddots&\\ 0&&\tilde{k}_{\text{I},z,(N,M)}\end{bmatrix},$ (52) $$

and $\mathbf{Y}_{\text{II}}$ and $\mathbf{Z}_{\text{II}}$ are

$$ $\displaystyle\mathbf{Y}_{\text{II}}=\begin{bmatrix}\tilde{k}_{\text{II},z,(-N, -M)}&&0\\ &\ddots&\\ 0&&\tilde{k}_{\text{II},z,(N,M)}\end{bmatrix},\quad\mathbf{Z}_{\text{II}}= \frac{1}{(\mathtt{n}_{\text{II}})^{2}}\begin{bmatrix}\tilde{k}_{\text{II},z,(- N,-M)}&&0\\ &\ddots&\\ 0&&\tilde{k}_{\text{II},z,(N,M)}\end{bmatrix}.$ (53) $$

Here, new set of $\mathbf{W}_{g}$ and $\mathbf{V}_{g}$ on SP basis $\{\hat{s},\hat{p}\}$ are introduced which are recombined from the set of $\mathbf{W}_{g}$ and $\mathbf{V}_{g}$ from XY basis $\{\hat{x},\hat{y}\}$:

$$ $\displaystyle\mathbf{W}_{g,ss}$ $\displaystyle=\mathbf{F}_{c}\mathbf{W}_{g,21}-\mathbf{F}_{s}\mathbf{W}_{g,11},$ $\displaystyle\mathbf{W}_{g,sp}$ $\displaystyle=\mathbf{F}_{c}\mathbf{W}_{g,22}-\mathbf{F}_{s}\mathbf{W}_{g,12},$ (54) $\displaystyle\mathbf{W}_{g,ps}$ $\displaystyle=\mathbf{F}_{c}\mathbf{W}_{g,11}+\mathbf{F}_{s}\mathbf{W}_{g,21},$ $\displaystyle\mathbf{W}_{g,pp}$ $\displaystyle=\mathbf{F}_{c}\mathbf{W}_{g,12}+\mathbf{F}_{s}\mathbf{W}_{g,22},$ (55) $\displaystyle\mathbf{V}_{g,ss}$ $\displaystyle=\mathbf{F}_{c}\mathbf{V}_{g,11}+\mathbf{F}_{s}\mathbf{V}_{g,21},$ $\displaystyle\mathbf{V}_{g,sp}$ $\displaystyle=\mathbf{F}_{c}\mathbf{V}_{g,12}+\mathbf{F}_{s}\mathbf{V}_{g,22},$ (56) $\displaystyle\mathbf{V}_{g,ps}$ $\displaystyle=\mathbf{F}_{c}\mathbf{V}_{g,21}-\mathbf{F}_{s}\mathbf{V}_{g,11},$ $\displaystyle\mathbf{V}_{g,pp}$ $\displaystyle=\mathbf{F}_{c}\mathbf{V}_{g,22}-\mathbf{F}_{s}\mathbf{V}_{g,12},$ (57) $$

with $\mathbf{F}_{c}$ and $\mathbf{F}_{s}$ being diagonal matrices with the diagonal elements $\cos\varphi_{(n,m)}$ and $\sin\varphi_{(n,m)}$, respectively, where

$$ $\displaystyle\varphi_{(n,m)}=\tan^{-1}(k_{y,n}/k_{x,m}).$ (58) $$

Equations 49 and 50 can be reduced to one set of equations by eliminating $\mathbf{c}^{\pm}_{1,2}$:

$$ $\displaystyle\begin{bmatrix}\sin\psi\ \boldsymbol{\delta}_{00}\\ \cos\psi\ \cos\theta\ \boldsymbol{\delta}_{00}\\ j\sin\psi\ \mathtt{n}_{\text{I}}\cos\theta\ \boldsymbol{\delta}_{00}\\ -j\cos\psi\ \mathtt{n}_{\text{I}}\ \boldsymbol{\delta}_{00}\\ \end{bmatrix}+\begin{bmatrix}\mathbf{I}&\mathbf{0}\\ \mathbf{0}&-j\mathbf{Z}_{I}\\ -j\mathbf{Y}_{I}&\mathbf{0}\\ \mathbf{0}&\mathbf{I}\end{bmatrix}\begin{bmatrix}\mathbf{R}_{s}\\ \mathbf{R}_{p}\end{bmatrix}=\begin{bmatrix}\mathbb{W}&\mathbb{WX}\\ \mathbb{V}&-\mathbb{VX}\end{bmatrix}\begin{bmatrix}\mathbb{WX}&\mathbb{W}\\ \mathbb{VX}&-\mathbb{V}\end{bmatrix}^{-1}\begin{bmatrix}\mathbb{F}\\ \mathbb{G}\\ \end{bmatrix}\begin{bmatrix}\mathbf{T}_{s}\\ \mathbf{T}_{p}\end{bmatrix},$ (59) $$

where

$$ (60) $$

This equation for a single layer grating can be simply extended to a multi-layer system as the following:

$$ (61) $$

where $L$ is the number of layers and

$$ $\displaystyle\mathbb{F}_{L+1}=\begin{bmatrix}\mathbf{I}&\mathbf{0}\\ \mathbf{0}&j\mathbf{Z}_{\text{II}}\end{bmatrix},\quad\mathbb{G}_{L+1}=\begin{ bmatrix}j\mathbf{Y}_{\text{II}}&\mathbf{0}\\ \mathbf{0}&\mathbf{I}\end{bmatrix}.$ (62) $$

Since we have four matrix equations for four unknown coefficients ($\mathbf{R}_{s},\mathbf{R}_{p},\mathbf{T}_{s},\mathbf{T}_{p}$), they can be derived and used for calculating diffraction efficiencies (also called the reflectance and transmittance).

The diffraction efficiency is the ratio of the power flux in propagating direction between incidence and diffracted wave of interest. It can be calculated by time-averaged Poynting vector [, , ]:

$$ $\displaystyle P$ $\displaystyle=\frac{1}{2}\operatorname{Re}{(E\times H^{*})},$ (63) $$

where ^∗ is the complex conjugate.
Now we can find the total power of the incident wave as a sum of the power of TE wave and TM wave:

$$ $\begin{split}P^{inc}&=P_{s}^{inc}+P_{p}^{inc}\\ &=\frac{1}{2}\operatorname{Re}\Bigg{[}{(E_{s}\times H_{s}^{*})+(E_{p}\times H_ {p}^{*})}\Bigg{]}\\ &=\frac{1}{2}\operatorname{Re}\Bigg{[}{(\sin\psi\ \cdot\sin\psi\ \mathtt{n}_{ \text{I}}\ \cos\theta)+(\cos\psi\ \cos\theta\ \cdot\cos\psi\ \mathtt{n}_{\text {I}})}\Bigg{]}\\ &=\frac{1}{2}\operatorname{Re}\Bigg{[}{(\sin^{2}\psi\ \mathtt{n}_{\text{I}} \cos\theta)+(\cos^{2}\psi\ \mathtt{n}_{\text{I}}\cos\theta)}\Bigg{]}\\ &=\frac{1}{2}\operatorname{Re}\Bigg{[}{(\mathtt{n}_{\text{I}}\cos\theta)}\Bigg {]}.\end{split}$ (64) $$

The power in each reflected diffraction mode is

$$ $\begin{split}P_{n,m}^{r}&=P_{nm,s}^{r}+P_{nm,p}^{r}\\ &=\frac{1}{2}\operatorname{Re}\Bigg{[}{(E_{nm,s}^{r}\times(H_{nm,s}^{r})^{*})+ (E_{nm,p}^{r}\times(H_{nm,p}^{r})^{*})}\Bigg{]}\\ &=\frac{1}{2}\operatorname{Re}\Bigg{[}{R_{nm,s}\cdot\frac{k_{\text{I},z,(n,m)} }{k_{0}}R_{nm,s}^{*}+\frac{k_{\text{I},z,(n,m)}}{k_{0}\mathtt{n}_{\text{I}}^{2 }}R_{nm,p}\cdot R_{nm,p}^{*}}\Bigg{]}\\ &=\frac{1}{2}\operatorname{Re}\Bigg{[}{R_{nm,s}R_{nm,s}^{*}\cdot\frac{k_{\text {I},z,(n,m)}}{k_{0}}+R_{nm,p}R_{nm,p}^{*}\cdot\frac{k_{\text{I},z,(n,m)}}{k_{0 }\mathtt{n}_{\text{I}}^{2}}}\Bigg{]},\end{split}$ (65) $$

and the power in each transmitted diffraction mode is

$$ $\begin{split}P_{n,m}^{t}&=P_{nm,s}^{t}+P_{nm,p}^{t}\\ &=\frac{1}{2}\operatorname{Re}\Bigg{[}{(E_{nm,s}^{t}\times(H_{nm,s}^{t})^{*})+ (E_{nm,p}^{t}\times(H_{nm,p}^{t})^{*})}\Bigg{]}\\ &=\frac{1}{2}\operatorname{Re}\Bigg{[}{T_{nm,s}\cdot\frac{k_{\text{II},z,(n,m) }}{k_{0}}T_{nm,s}^{*}+\frac{k_{\text{II},z,(n,m)}}{k_{0}\mathtt{n}_{\text{II}} ^{2}}T_{nm,p}\cdot T_{nm,p}^{*}}\Bigg{]}\\ &=\frac{1}{2}\operatorname{Re}\Bigg{[}{T_{nm,s}T_{nm,s}^{*}\cdot\frac{k_{\text {II},z,(n,m)}}{k_{0}}+T_{nm,p}T_{nm,p}^{*}\cdot\frac{k_{\text{II},z,(n,m)}}{k_ {0}\mathtt{n}_{\text{II}}^{2}}}\Bigg{]}.\end{split}$ (66) $$

Since the diffraction efficiency is the ratio between them $(P_{out}/P_{inc})$, we can get the efficiencies of reflected and transmitted waves:

$$ $\displaystyle DE_{r,(n,m)}$ $\displaystyle=|R_{s,(n,m)}|^{2}\operatorname{Re}{\bigg{(}\frac{k_{\text{I},z,( n,m)}}{k_{0}\mathtt{n}_{\text{I}}\cos{\theta}}\bigg{)}}+|R_{p,(n,m)}|^{2} \operatorname{Re}{\bigg{(}\frac{k_{\text{I},z,(n,m)}/{\mathtt{n}_{\text{I}}}^{ 2}}{k_{0}\mathtt{n}_{\text{I}}\cos{\theta}}\bigg{)}},$ (67) $\displaystyle DE_{t,(n,m)}$ $\displaystyle=|T_{s,(n,m)}|^{2}\operatorname{Re}{\bigg{(}\frac{k_{\text{II},z, (n,m)}}{k_{0}\mathtt{n}_{\text{I}}\cos{\theta}}\bigg{)}}+|T_{p,(n,m)}|^{2} \operatorname{Re}{\bigg{(}\frac{k_{\text{II},z,(n,m)}/{\mathtt{n}_{\text{II}}} ^{2}}{k_{0}\mathtt{n}_{\text{I}}\cos{\theta}}\bigg{)}}.$ (68) $$

### G.5 Enhanced Transmittance Matrix Method

As addressed in [, , ], solving Equation 61 may suffer from the numerical instability coming from the inversion of almost singular matrix when $\mathbb{X}_{\ell}$ has a very small and possibly numerically zero value. meent adopted Enhanced Transmittance Matrix Method (ETM) [] to overcome this by avoiding the inversion of $\mathbb{X}_{\ell}$.

The technique is sequentially applied from the last layer to the first layer. In Equation 61, the set of modes at the bottom interface of the last layer $(\ell=L)$ is

$$ $\begin{split}&\begin{bmatrix}\mathbb{W}_{L}&\mathbb{W}_{L}\mathbb{X}_{L}\\ \mathbb{V}_{L}&-\mathbb{V}_{L}\mathbb{X}_{L}\end{bmatrix}\begin{bmatrix} \mathbb{W}_{L}\mathbb{X}_{L}&\mathbb{W}_{L}\\ \mathbb{V}_{L}\mathbb{X}_{L}&-\mathbb{V}_{L}\end{bmatrix}^{-1}\begin{bmatrix} \mathbb{F}_{L+1}\\ \mathbb{G}_{L+1}\end{bmatrix}\begin{bmatrix}\mathbf{T}_{s}\\ \mathbf{T}_{p}\end{bmatrix}\\ &=\begin{bmatrix}\mathbb{W}_{L}&\mathbb{W}_{L}\mathbb{X}_{L}\\ \mathbb{V}_{L}&-\mathbb{V}_{L}\mathbb{X}_{L}\end{bmatrix}\begin{bmatrix}{ \mathbb{X}_{L}}^{-1}&\mathbb{0}\\ \mathbb{0}&{\mathbb{I}}\\ \end{bmatrix}{\begin{bmatrix}\mathbb{W}_{L}&\mathbb{W}_{L}\\ \mathbb{V}_{L}&-\mathbb{V}_{L}\end{bmatrix}}^{-1}\begin{bmatrix}\mathbb{F}_{L+ 1}\\ \mathbb{G}_{L+1}\end{bmatrix}\begin{bmatrix}\mathbf{T}_{s}\\ \mathbf{T}_{p}\end{bmatrix}.\\ \end{split}$ (69) $$

The matrix to be inverted can be decomposed into two matrices by isolating $\mathbb{X}_{L}$, which is the potential source of the numerical instability. The right-hand side can be shortened with new variables $\mathbb{A}_{L},\mathbb{B}_{L}$:

$$ $\begin{split}\begin{bmatrix}\mathbb{A}_{L}\\ \mathbb{B}_{L}\end{bmatrix}=\begin{bmatrix}{\mathbb{W}_{L}}&\mathbb{W_{L}}\\ \mathbb{V_{L}}&{-\mathbb{V}_{L}}\\ \end{bmatrix}^{-1}\begin{bmatrix}\mathbb{F}_{L+1}\\ \mathbb{G}_{L+1}\end{bmatrix},\end{split}$ (70) $$

then the right-hand side of Equation LABEL:eqn:etm-lastlayer becomes

$$ $\begin{split}\begin{bmatrix}\mathbb{W}_{L}&\mathbb{W}_{L}\mathbb{X}_{L}\\ \mathbb{V}_{L}&-\mathbb{V}_{L}\mathbb{X}_{L}\end{bmatrix}\begin{bmatrix}{ \mathbb{X}_{L}}^{-1}&\mathbb{0}\\ \mathbb{0}&{\mathbb{I}}\\ \end{bmatrix}\begin{bmatrix}\mathbb{A}_{L}\\ \mathbb{B}_{L}\end{bmatrix}\begin{bmatrix}\mathbf{T}_{s}\\ \mathbf{T}_{p}\end{bmatrix}.\\ \end{split}$ (71) $$

We can avoid the inversion of $\mathbb{X}_{L}$ by introducing the substitution $\mathbf{T}_{s}={\mathbb{A}_{L}}^{-1}\mathbb{X}_{L}\mathbf{T}_{s,L}$ and $\mathbf{T}_{p}={\mathbb{A}_{L}}^{-1}\mathbb{X}_{L}\mathbf{T}_{p,L}$. Equation 71 then becomes

$$ $\begin{split}&\begin{bmatrix}\mathbb{W}_{L}&\mathbb{W}_{L}\mathbb{X}_{L}\\ \mathbb{V}_{L}&-\mathbb{V}_{L}\mathbb{X}_{L}\end{bmatrix}\begin{bmatrix}{ \mathbb{X}_{L}}^{-1}&\mathbb{0}\\ \mathbb{0}&{\mathbb{I}}\\ \end{bmatrix}\begin{bmatrix}\mathbb{A}_{L}\\ \mathbb{B}_{L}\end{bmatrix}{\mathbb{A}_{L}}^{-1}{\mathbb{X}_{L}}\begin{bmatrix }\mathbf{T}_{s,L}\\ \mathbf{T}_{p,L}\end{bmatrix}\\ &=\begin{bmatrix}\mathbb{W}_{L}&\mathbb{W}_{L}\mathbb{X}_{L}\\ \mathbb{V}_{L}&-\mathbb{V}_{L}\mathbb{X}_{L}\end{bmatrix}\begin{bmatrix}{ \mathbb{X}_{L}}^{-1}&\mathbb{0}\\ \mathbb{0}&{\mathbb{I}}\\ \end{bmatrix}\begin{bmatrix}\mathbb{X}_{L}\\ \mathbb{B}_{L}\mathbb{A}_{L}^{-1}\mathbb{X}_{L}\end{bmatrix}\begin{bmatrix} \mathbf{T}_{s,L}\\ \mathbf{T}_{p,L}\end{bmatrix}\\ &=\begin{bmatrix}\mathbb{W}_{L}&\mathbb{W}_{L}\mathbb{X}_{L}\\ \mathbb{V}_{L}&-\mathbb{V}_{L}\mathbb{X}_{L}\end{bmatrix}\begin{bmatrix} \mathbb{I}\\ \mathbb{B}_{L}\mathbb{A}_{L}^{-1}\mathbb{X}_{L}\end{bmatrix}\begin{bmatrix} \mathbf{T}_{s,L}\\ \mathbf{T}_{p,L}\end{bmatrix}\\ &=\begin{bmatrix}\mathbb{W}_{L}(\mathbb{I}+\mathbb{X}_{L}\mathbb{B}_{L}\mathbb {A}_{L}^{-1}\mathbb{X})\\ \mathbb{V}_{L}(\mathbb{I}-\mathbb{X}_{L}\mathbb{B}_{L}\mathbb{A}_{L}^{-1} \mathbb{X})\end{bmatrix}\begin{bmatrix}\mathbf{T}_{s,L}\\ \mathbf{T}_{p,L}\end{bmatrix}\\ &=\begin{bmatrix}\mathbb{F}_{L}\\ \mathbb{G}_{L}\end{bmatrix}\begin{bmatrix}\mathbf{T}_{s,L}\\ \mathbf{T}_{p,L}\end{bmatrix}.\end{split}$ (72) $$

These steps can be repeated until the iteration gets to the first layer $(\ell=1)$, then the form becomes

$$ $\displaystyle\begin{bmatrix}\sin\psi\ \boldsymbol{\delta}_{00}\\ \cos\psi\ \cos\theta\ \boldsymbol{\delta}_{00}\\ j\sin\psi\ n_{\text{I}}\ \cos\theta\ \boldsymbol{\delta}_{00}\\ -j\cos\psi\ n_{\text{I}}\ \boldsymbol{\delta}_{00}\\ \end{bmatrix}+\begin{bmatrix}\mathbf{I}&\mathbf{0}\\ \mathbf{0}&-j\mathbf{Z}_{I}\\ -j\mathbf{Y}_{I}&\mathbf{0}\\ \mathbf{0}&\mathbf{I}\end{bmatrix}\begin{bmatrix}\mathbf{R}_{s}\\ \mathbf{R}_{p}\end{bmatrix}=\begin{bmatrix}\mathbb{F}_{1}\\ \mathbb{G}_{1}\end{bmatrix}\begin{bmatrix}\mathbf{T}_{s,1}\\ \mathbf{T}_{p,1}\end{bmatrix},$ (73) $$

where

$$ $\begin{bmatrix}\mathbf{T}_{s}\\ \mathbf{T}_{p}\end{bmatrix}=\mathbb{A}_{L}^{-1}\mathbb{X}_{L}\cdots\mathbb{A}_ {\ell}^{-1}\mathbb{X}_{\ell}\cdots\mathbb{A}_{1}^{-1}\mathbb{X}_{1}\begin{ bmatrix}\mathbf{T}_{s,1}\\ \mathbf{T}_{p,1}\end{bmatrix}.$ $$

### G.6 Topological derivative vs Shape derivative

AD enables the calculation of the gradient of the figure of merit (FoM) with respect to the design parameters of the device.
AD in meent can handle both modeling type - raster and vector - with two different forms: topological derivative and shape derivative.
If the raster modeling is utilized to obtain the geometry of the device, the gradients with respect to the refractive index of every pixel can be obtained through AD. This type of AD is known as the topological derivative, as the device design is updated pixel-wise and the topology is not conserved (Figure 15a).
On the contrary, a shape derivative is effective for vector modeling; the FoM derivative with respect to input dimensions is obtained as depicted in Figure 15b. The shape derivative is expected to be useful for cases where the device topology is known, but dimensions of specific structures, such as the radius of a cylinder in a layer or width and length of a cuboid, are to be found by optimization.

Figure: Figure 15: Topological and shape derivatives. A schematic diagram showing the difference between the (a) topological derivative and (b) shape derivative. The topological derivative results in the FoM derivative with respect to the permittivity changes of every cells in the grid and the shape derivative provides the derivative with respect to the deformations of a shape.
Refer to caption: x6.png

## Appendix H Program sequence

In this section, we will provide a detailed explanation of the functions in meent(^5^55for version 0.9.x) and discuss the simulation program sequence with examples.

### H.1 Initialization

A simple way to use meent is using ‘call_mee()’ function which returns an instance of Python class that includes all the functionalities of meent. Simulation conditions can be set by passing parameters as arguements (args) or keyword arguements (kwargs) in this function. It is also possible to change conditions after calling instance by directly assigning desired value to the property of the instance.

Figure: Code 3: Methods to set simulation conditions

Here are the descriptions of the input parameters in meent class:

: integer
meent supports three backends: NumPy, JAX, and PyTorch.

- •
0: NumPy (RCWA only; AD is not supported).
- •
1: JAX.
- •
2: PyTorch.

: integer
This parameter defines the simulation space.

- •
0: 1D grating without conical incidence $(\phi=0)$.
- •
1: 1D grating with conical incidence.
- •
2: 2D grating.

: integer or float
This parameter controls the linear polarization state of the incident wave by this definition: $\psi=\pi/2*(1-{\textit{pol}})$.
It can take values between 0 and 1. 0 represents fully transverse electric (TE) polarization, and 1 represents fully transverse magnetic (TM) polarization. Support for other polarization states such as the circular polarization state which involves the phase difference between TE and TM polarization will be added in the future updates.

: float
The refractive index of the superstrate.

: float
The refractive index of the substrate.

: float
The angle of the incidence in radians.

: float
The angle of rotation (or azimuth angle) in radians.

: float
The wavelength of the incident light in vacuum. Future versions may support complex type wavelength.

: integer or list of integers
Fourier truncation order (FTO). This represents the number of Fourier harmonics in use. If fourier_order = $N$, this is for 1D grating and meent utilizes $(2N+1)$ harmonics spanning from $-N$ to $N$:$-N,-(N-1),...,N$. For 2D gratings, it takes a sequence $[M,N]$ as an input, where $M$ and $N$ become FTO in $X$ and $Y$ directions, respectively. Note that 1D grating can also be simulated in 2D grating system by setting $N$ as 0 0.

: list of floats
The period of a unit cell. For 1D grating, it is a sequence with one element which is a period in X-direction. For 2D gratings, it takes a sequence [period in $X$, period in $Y$] as an input.

: integer
The datatype used in the simulation.

- •
0: complex128 (64 bit).
- •
1: complex64 (32 bit).

: integer
The selection of the device for the calculations: currently CPU and GPU are supported. At the time of writing this paper, the eigendecomposition, which is the most expensive step as $\mathcal{O}(M^{3}N^{3})$ where $M\text{ and }N$ are FTO, is available only on CPU. This means GPU may not as powerful as we expect as in deep learning regime.

- •
0: CPU.
- •
1: GPU.

: integer
This variable selects the type of Fourier series implementation. 0 and 1 are options for raster modeling and 2 is for vector modeling. 0 uses discrete Fourier series (DFS) while 1 and 2 use continuous Fourier series (CFS). Note that the name ‘fft_type’ may change since it is not correct expression.

- •
0: DFS for the raster modeling (pixel-based geometry). fft_type supports improve_dft option, which is True by default, that can prevent aliasing by increasing sampling frequency, and drives the result to approach to the result of CFS.
- •
1: CFS for the raster modeling (pixel-based geometry). This doesn’t support backpropagation. Use this option for debugging or in RCWA-only situation.
- •
2: CFS for the vector modeling (object-based geometry).

: list of floats
The sequence of the thickness of each layer from top to bottom.

: array of {floats, complex numbers}, shape is (i, j, k)
The input for the raster modeling. It takes a 3D array in ($Z$,$Y$,$X$) order, where $Z$ represents the direction of the layer stacking. In case of 1D grating, j is 1 (e.g., shape = (3,1,10) for a stack composed of 3 layers that are 1D grating).

### H.2 Structure design

meent provides two types of structure design methods: the vector modeling and the raster modeling.

#### H.2.1 vector modeling

Figure: (a) n_split: (1, 1)
Refer to caption: x7.png

Figure 16 shows rotated rectangles drawn on XY plane. meent decomposes the geometrical figures into the collection of sub-rectangles which of each side lies on the direction of either $\hat{x}$ or $\hat{y}$. Then CFS with the sinc function is used to find the Fourier coefficients. The degree of approximation can be determined by ‘n_split’ option in Code 4.

To add primitives to the simulation space, users can utilize ‘rectangle()’ or ‘rectangle_rotation()’ functions which allows the insertion of desired geometry. The ‘draw()’ function is then employed to create the final structure, taking into account any potential overlaps between the geometries. Code 4 is the example creating a layer that has rotated rectangle.

Figure: Code 4: vector modeling

Figure: (a) Red rectangle comes first, blue rectangle last

Code 5 and Figure 17 show how meent can handle the overlap of the shapes. Figure 17(a) and 17(b) have the same set of rectangles (red and blue) but they are placed in different order and this can be controlled by the function ‘layer_info_list’ in Code 5. It is the list that contains the base refractive index of the layer and the primitive shapes to be placed on the layer. In case of Figure 17(a), red rectangle comes first in the list and blue does for Figure 17(b).

Figure: Code 5: overlap

#### H.2.2 Raster modeling

Figure: (a) 1D grating
Refer to caption: extracted/5658627/sections/manual/sequence/figures/ucell_1d.png

We have 2 example structures of raster modeling as shown in Figure 18 and Code 6.
Figure 18(a) is a stack of 2 layers which has 1D grating. Note that 1D grating unit cell can be defined by setting the length of the second axis to 1 as (a) in Code 6. Figure 18(b) is a stack of single 2D grating layer.

Figure: Code 6: Raster modeling

### H.3 Electromagnetic simulation

Electromagnetic simulation (EM simulation) in meent can be divided into 3 main subcategories: convolution matrix generation, Maxwell’s equations computation and field calculation.
The method ‘conv_solve()’ does both convolution matrix generation and Maxwell’s equations computation sequentially. ‘conv_solve_field()’ method does the same and additionally calculates the field distribution of the structure. Code 7 is the example showing how to use those; ‘conv_solve()’ method returns the reflected and transmitted diffraction efficiencies and ‘conv_solve_field()’ does both and field distribution.

Figure: Code 7: Method call for EM simulation

#### H.3.1 Convolution matrix generation

Figure: (a) Coefficients matrix

The functions for convolution matrix generation are located in ‘convolution_matrix.py’ file for each backend. This part transforms the structure from the real space to the Fourier space and returns a convolution matrix (also called Toeplitz matrix) of the Fourier coefficients to apply convolution operation with the E and H fields. Figure 19 shows the Fourier coefficients matrix and convolution matrix made from the coefficient matrix.
Code 8 is the definition of ‘conv_solve()’ method and shows how the convolution matrix generation is integrated inside. As shown in the code, meent offers 3 different methods to get convolution matrix since each method has different input type and implementation. This can be chosen by the argument ‘fft_type’: 0 is for raster modeling with DFS, 1 for raster with CFS and 2 for vector with CFS.

Figure: Code 8: ‘conv_solve()’

#### H.3.2 Maxwell’s equations computation

After generating the convolution matrix, meent solves Maxwell’s equations and returns diffraction efficiencies with the method ‘solve()’. As in the Code 9, it is a wrapper of ‘_solve()’ method that actually does the calculations and returns the diffraction efficiencies with other information that is necessary for the field calculation.

Figure: Code 9: ‘solve()’

Input parameters:

: float
The wavelength of the incident light in vacuum.

: array of {float or complex}
A stack of convolution matrices of the permittivity array; this is $\llbracket\varepsilon_{r,g}\rrbracket$ in Chapter G. The order of the axes is the same as of ucell ($Z$ $Y$ $X$).

: array of {float or complex}
A stack of convolution matrices of the one-over-permittivity array; this is $\llbracket\varepsilon_{r,g}^{-1}\rrbracket$ in Chapter G. The order of the axes is the same as of ucell ($Z$ $Y$ $X$).

The diffraction efficiencies are 1D array for 1D and 1D-conical grating and 2D for 2D grating.

#### H.3.3 Field calculation

The ‘calculate_field()’ method in Code 10 calculates the field distribution inside the structure. Note that the ‘solve()’ method must be preceded.
This function returns 4 dimensional array that the length of the last axis varies depending on the grating type as shown in Code 11. 1D TE and TM has 3 elements (TE has Ey, Hx and Hz in order and TM has Hy, Ex and Ez) while the others have 6 elements (Ex, Ey, Ez, Hx, Hy and Hz) as in Figure 20.

Figure: Code 10: ‘calculate_field()’

Input parameters:

: integer
The field resolution in X direction (number of split which the period of x is divided by).

: integer
The field resolution in Y direction (number of split which the period of y is divided by).

: integer
The field resolution in Z direction (number of split in thickness of each layer).

: integer
The level of vectorization for the field calculation. Default is 2 which is fully vectorized for fast calculation while 1 is half-vectorized and 0 is none. Option 0 and 1 are remained for debugging or future development (such as parallelization).

- •
0: Non-vectorized
- •
1: Semi-vectorized: in X and Y direction
- •
2: Vectorized: in X, Y and Z direction

Figure: Code 11: the shape of returned array from ‘calculate_field()’

Figure: (a) Ex

## Appendix I Benchmark

In this section, we will address the 1D metasurface problem covered in our previous work [] with meent so that we can benchmark and analyze its capability and functionality.

### I.1 case application

Figure: Figure 21: The image of 1D diffraction metagrating on a silicon dioxide substrate.
Refer to caption: x19.png

This metagrating deflector is composed of silicon pillars placed on a silica substrate. The device period is divided into 64 cells, and each cell can be filled with either air or silicon. The Figure of Merit for this optimization is set to the deflection efficiency of the +1^st order transmitted wave when TM polarized wave is normally incident from the silica substrate as in Figure 21.

### I.2 Fourier series implementations

When the sampling frequency of permittivity distribution is not enough, Fourier coefficients from DFS is aliased. It can be resolved by increasing the sampling rate that is implemented in the way of duplicating the elements so the array is extended to have identical distribution but larger array size. We will call this Enhanced DFS, and it’s implemented in meent as a default option.

Figure: (a) Convergence test
Refer to caption: extracted/5658627/sections/manual/images/benchmarks/convergence_bad.png

Figure 22(a) illustrates the convergence tests of a particular structure with four different RCWA implementations.
Considering Reticolo as the reference, we can see CFS is well-matched but DFS shows different behavior. This is due to the insufficient sampling rate of permittivity distribution, which can be resolved by Enhanced DFS. Figure 22(b) is the histogram of the discrepancies from Reticolo result. About 600k structures were evaluated with 4 implementations and the errors of 3 meent implementations were calculated based on Reticolo. CFS shows the smallest errors and this is because Reticolo too uses CFS (CFS algorithms in meent are adopted from Reticolo). Enhanced DFS decreases the error about three orders of magnitudes (e.g., the median of DFS is 4.3E-4 and this becomes 1.4E-7).

### I.3 Computing performance

**Table 10: Performance test condition**
| backend | device | bit | alpha server | beta server | gamma server |
| --- | --- | --- | --- | --- | --- |
| NumPy | CPU | 64 | (A1) | (B1) | (C1) |
| NumPy | CPU | 32 | (A2) | (B2) | (C2) |
| JAX | CPU | 64 | (A3) | (B3) | (C3) |
| JAX | CPU | 32 | (A4) | (B4) | (C4) |
| JAX | GPU | 64 | - | (B5) | (C5) |
| JAX | GPU | 32 | - | (B6) | (C6) |
| PyTorch | CPU | 64 | (A7) | (B7) | (C7) |
| PyTorch | CPU | 32 | (A8) | (B8) | (C8) |
| PyTorch | GPU | 64 | - | (B9) | (C9) |
| PyTorch | GPU | 32 | - | (B10) | (C10) |

In this section, computing options to speed up the calculation - backend, device (CPU and GPU) and architecture (64bit and 32bit) - will be benchmarked. Table 1 is the hardware specification of the test server and Table 10 is the index of each test condition.

Figure: Figure 23: Performance test: calculation time with respect to FTO. Top row is the result from 64bit and bottom is from 32bit. The first column is the result from the test server alpha and the rest is beta and gamma in order.
Refer to caption: extracted/5658627/sections/manual/images/benchmarks/performance/result_all.png

The graphs in Figure 23 are calculation time vs FTO with all the data per machine and architecture. Before look into the details, we will briefly mention some notice in this figure. (1) JAX can’t afford large FTO regardless of device. We suspect that this is related to JIT compilation which takes much time and memory for the compilation at the first run. (2) GPU with JAX and PyTorch can’t accept large FTO even though GPU memory is more than needed for array upload. (3) if large amount of calculation is needed, Numpy or PyTorch on CPU is the option. (4) no golden option exists: it is recommended to find the best option for the test environment by doing benchmark tests.

We will visit these computing options one by one. The option C9 at FTO 1600 will be excluded in further analyses: this seems an optimization issue in PyTorch or CUDA.

#### I.3.1 Backend: NumPy, JAX and PyTorch

NumPy, JAX and PyTorch as a backend are benchmarked. NumPy is installed via PyPI which is compiled with OpenBLAS. There are many types of BLAS libraries and the most representative ones are OpenBLAS and MKL (Math Kernel Library).
As of now, PyPI provides NumPy with OpenBLAS while conda does one with MKL. This makes small discrepancy in terms of speed and precision hence pay attention when doing consistency test between machines.

Figure: Figure 24: Performance test: calculation time by FTO sweep. The result is normalized by NumPy case from the same options to compare the behavior of other backends. In these plots, black dashed line is $y=1$ and the results of NumPy cases lie on this line since they are normalized by themselves.
Refer to caption: extracted/5658627/sections/manual/images/benchmarks/performance/backend.png

Figure 24 is the relative simulation time per server and architecture normalized by the time of NumPy case in the same conditions to make comparison easy.
In small FTO regime, all the options were successfully operated and no champion exists. Hence it is strongly recommended to run benchmark test on your hardware and pick the most efficient one. In case of X7 (A7, B7 and C7), Alpha and Gamma show the same behavior - spike in 100 - while beta shows fluctuation around B1. One possible reason for this is the type of CPU. The CPUs of Alpha and Gamma belong to ‘Xeon Scalable Processors’ group but Beta is ‘Xeon E Processors’. Currently we don’t know if this actually makes difference or some other reason (such as the number of threads or BLAS implementation) does. This result may vary if MKL were used instead of OpenBLAS.
In large FTO, only two options are available: NumPy and PyTorch on CPU in 64 bit. In case of JAX, the tests were failed: we watched memory occupation surge during the simulation which seems unrelated to matrix calculation. This might be an issue of JIT (Just In Time) compilation in JAX. Between NumPy and PyTorch, PyTorch is about twice faster than NumPy in both architectures at Alpha and Gamma, but beta shows different behavior. This too, we don’t know the root cause but one notable difference is the family of CPU type.

#### I.3.2 Device: CPU and GPU

Figure: Figure 25: Performance test result. The calculation time of GPU cases are normalized by CPU cases from the same options to see the efficiency of GPU utilization. In these plots, black dashed line is $y=1$ where the capability of both are the same.
Refer to caption: extracted/5658627/sections/manual/images/benchmarks/performance/device.png

Figure 25 shows the relative simulation time of GPU cases normalized by CPU cases on the same backend and architecture. Note that it is relative time, so the smaller time does not mean it is a good option for the simulation experiments: the relative time can be small even if the absolute time of CPU and GPU are very large compared to other options.

JAX shows good GPU utilization throughout the whole range (except one point in beta) regardless of the architecture. Considering the architecture, the data trend in beta is not clear while the gamma clearly shows that GPU utilization can be more effective in 32bit operation. PyTorch data is a bit noisier than of JAX, but has the similar behavior per server. The data in beta is hard to conclude as the JAX cases and the gamma too shows ambiguous trend but we can consider GPU option is efficient with wide range of FTOs.

Up to date, eigendecomposition for non-hermitian matrix which is the most expensive step ($O(M^{3}N^{3})$) in RCWA, is not implemented on GPU in JAX and PyTorch hence the calculations are done on CPU and the results are sent back to GPU. As a result, we cannot expect great performance enhancement in using GPUs.

#### I.3.3 Architecture: 64 and 32 bit

Figure: Figure 26: Performance test result. The calculation time of 32bit cases are normalized by 64bit cases from the same options. In these plots, black dashed line is $y=1$ where the capability of both are the same.
Refer to caption: extracted/5658627/sections/manual/images/benchmarks/performance/archi.png

In Figure 26, calculation time of 32bit case is normalized by 64bit case with the same condition. With some exceptions, most points show that simulation in 32bit is faster than 64bit. Here are some important notes:
(1) From our understanding, the eigendecomposition (Eig) in NumPy operates in 64bit regardless of the input type - even though the input is 32bit data (float32 or complex64), the matrix operations inside Eig are done in 64bit but returns the results in 32bit data type. This is different from JAX and PyTorch - they provides Eig in 32bit as well as 64bit. Hence the 32bit NumPy cases in the figure approach to 1 as FTO increases because the calculation time for Eig is the same and it is the most time-consuming step.
(2) Keep in mind that 32bit data type can handle only 8 digits. This means that 1000 + 0.00001 becomes 1000 without any warnings or error raises. For such a reason, the accuracy of 32bit cases in the figures are not guaranteed - we only consider the calculation time.
(3) Eig in PyTorch shows interesting behavior: as FTO increases, calculation time in 32bit overtakes 64bit - see A8/A7, B8/B7 and C8/C7. This is counter-intuitive and we don’t have good explanation but cautiously guess that this might be related to the accuracy and precision in Eig or an optimization issue of PyTorch.

## Appendix J Applications

meent is expected to be useful for solving inverse design and optimization problems in OCD and metasurface design. In this section, we present some exemplary cases where meent proves its capabilities.
Leveraging the automatic differentiation function, we successfully carry out optimization for diffraction gratings and achieve inverse design of the geometric parameters.

### J.1 Inverse design of 1D diffraction grating

Figure: (a) Optimization curve
Refer to caption: extracted/5658627/sections/manual/images/benchmarks/optim_a.png

In this example, we will optimize 1D beam deflector that was used for benchmark in Chapter I using AD with these options - 256 cells, FTO $=100$, $\lambda_{0}=900$ nm and the deflection angle $=50^{\circ}$.
During the optimization, each cell can have non-binary refractive index values, leading to a gray-scale optimization. To obtain the final structure consisting of only silicon/air binary structures, an additional binary-push process is required.
The initial structure for optimization is randomly generated so the each cell can have the refractive index value between of air and silicon under uniform distribution.
The Figure of Merit for this optimization process is set to the $+1^{st}$ order diffraction efficiency, and the gradient is calculated by AD. The refractive indices are updated over multiple epochs using the ADAM optimizer [] with the learning rate of 0.5.

Figure 27(a) shows the deflection efficiency change by iteration. Two solid lines are averaged value of all the samples at the same iteration step. Shaded area is marked with $\pm$ standard deviation from the average. The blue line (Before binarization) is the result of device with any real number between two refractive indices (silicon and air), which is non-practical, and the orange line (After binarization) is the final device composed of silicon and air. The best result we found is 89.4%.

### J.2 Inverse design of 2D diffraction grating

Figure: Figure 28: Optimization result of 2D beam deflector. (a) The schematic of 2D beam deflector and the final structure after optimization. (b) Convergence test of the initial structure. (c) Learning curve of structure optimization for 110 epochs. Spatial blurring and binary push is applied on each epoch (d) The electric field distribution of the optimized structure in XZ plane.
Refer to caption: x20.png

Here, we demonstrate optimization of a 2D diffraction metagrating as shown in Figure 28a. Similar to the previous 1D diffraction metagrating, the 2D diffraction metagrating also consists of silicon pillars located on top of a silicon dioxide substrate. TM polarized wave with $\lambda=1000$ nm is normally incident from the bottom of the substrate and the device is designed to deflect the incident light with deflection angle $\theta=60^{\circ}$ in $X$-direction. The device has a rectangular unit cell of period $\lambda/\sin\theta\approx 1150$ nm and $\lambda/2=500nm$ for the x and y-axis, respectively. Moreover, the unit cell is gridded into $256\times 128$ cells which is either filled by air or silicon.
The convergence of RCWA simulation for different number of Fourier harmonics are plotted in Figure 28b. Considering the trade-off between simulation accuracy and time, we set $N_{x}=13$ and $N_{y}=10$.

After 110 epochs of optimization, the final structure achieves an efficiency of 92% and successfully deflects the incoming beam at a 60^∘ angle (Figure 28d). The optimized structure and the learning curve are presented in Figure 28a and Figure 28c, respectively.

### J.3 Inverse design of 1D grating color router

Figure: Figure 29: A schematic of a color router. The incoming light is guided to subpixels of corresponding wavelength.
Refer to caption: extracted/5658627/sections/manual/images/metacolorSchematic2.png

Until now, we have focused on the problems where the FoM was simply defined. However, in this example, we aim to demonstrate the optimization process of a meta color router, which involves a complex FoM.

A meta color router is an optical component designed for next-generation image sensors. It is designed to route the incoming light to the subpixel region of corresponding color, as depicted in Figure 29. In this exemplary case, we consider an RGB meta color router featuring a pixel pitch of $0.5\mu m$ that consists of vertically stacked 1D binary gratings. The constituent dielectrics are silicon dioxide and silicon nitride, with fixed refractive indices of 1.5 and 2.0, respectively.
The meta device region (width of $1.5\mu m$ and height of $2\mu m$) is sliced into 8 layers with 64 cells per layer.

The FoM for this meta color router is defined as the average of TE mode electric field intensity over the corresponding subpixel region, as given by Equation (74).

$$ $\displaystyle FoM=\frac{1}{N}\sum_{\lambda_{1}}^{\lambda_{N}}\frac{\int_{x_{1} }^{x_{2}}|\vec{\mathbf{E}}(\lambda)|^{2}dx}{\int_{0}^{P}|\vec{\mathbf{E}}( \lambda)|^{2}dx}\times T(\lambda)$ (74) $$

Here, $\mathbf{E}$ represents the electric field within the subpixel region, while $T$ represents transmittance. The parameter $x\in(x1,x2)$ determines the desired subpixel region corresponding to the incident beam wavelength. For simplicity, we define the wavelength ranges for R, G, and B as 600 nm - 700 nm, 500 nm - 600 nm, and 400 nm - 500 nm, respectively. Throughout the optimization process, optical efficiencies are averaged across 9 wavelength points to ensure a finely tuned broadband response.

The optimization procedure for the meta color router follows a similar approach to the previous examples, including random initialization, optimization via back-propagated gradients with or without binary push. The optimization curve and the final binarized device structure are shown in Figure 30.

Figure: Figure 30: Optimization result of 1D grating meta color router. (a) Optimization curve of greyscale device and binary-pushed device at each epoch. (b) Color sorting efficiency spectrum. (c) The electric field inside the final color router device.
Refer to caption: x21.png

## Appendix K Licenses

Development

- •
Numpy: BSD license
- •
JAX: Apache License 2.0
- •
PyTorch: BSD license (BSD-3)

Experiment

- •
Ray: Apache License 2.0
- •
SheepRL: Apache License 2.0
- •
neuraloperator: MIT license

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: - What is the technical novelty of this work?
- This work might be not fit to ICLR, because it proposes a simulation tool.  I would like to discuss this issue with the authors, reviewers, and area chairs.
soundness: 2
strengths: - It is generally well-written.
- It introduces machine learning researcher-friendly electromagnetic simulation software.
- It provides several interesting applications of meent.
summary: - This work emphasizes the significance of electromagnetic simulation in the analysis and design of photonic structures. 
- The integration of machine learning into electromagnetic simulation has emerged as a promising solution, with the optics research community increasingly leveraging machine learning algorithms.
- The paper introduces meent, a user-friendly electromagnetic simulation software developed in Python that utilizes rigorous coupled-wave analysis and automatic differentiation.
- It includes three key applications of meent.
weaknesses: - I don't know what the technical novelty of this work is.

---

confidence: 2
contribution: 4
initial_rating: 6
presentation: 3
questions: See weakness.
soundness: 3
strengths: - Combining ML with EM simulator is an important problem. 
- The paper is organized really well and well-written. Starting with the technical details of meent, the authors also present concrete applications. This helps significantly in illustrating the value of meent. 
- The contribution of the proposed framework has been clearly discussed and the comparison with existing packages is comprehensive.
summary: This work presents meent, a framework that aims to integrate EM simulators into the ML pipelines. In particular, meent contains a differentiable, Python-native EM simulator. The authors demonstrate the value of meent through three concrete applications: (1) generating datasets to train the neural operators; (2) enabling RL-based design of nanophotonic device; (3) constructing solutions for  inverse-problems
weaknesses: - Missing related work: *Benchmarking Data-driven Surrogate Simulators for Artificial Electromagnetic Materials.* Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2) 
- I think it would be helpful to present a set of experiments comparing the efficiency of meent. In particular, how fast does meent generate EM simulation comparing to existing Cpp-based methods? 
- In other words, does the whole ML+EM pipeline take longer time by incorporating differentiable EM simulation? If the extra time is considerable, then the benefits of meent should be evaluated more carefully. 
- Since meent has some simulation error as discussed in Section 3, what is the performance degrade comparing to ML + classical EM simulators (e.g., Reticolo)? I think **it is helpful conduct the same set of experiments with ML+ Reticolo** and compare the resultant performance with the ones under meent.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 3
questions: Could you either provide a comparison of your method to MaxwellNet (Lim & Psaltis, 2022)  for electric field prediction, or explain why MaxwellNet was not included as a baseline?"
soundness: 3
strengths: - Clear presentation
- Comprehensive set of applications, including investigating machine learning (ML) algorithms in optics problems, and on development of nanophotonic devices.
summary: A differentiable electromagnetic simulation framework (called "meent"), which is able to operate on a continuous space. Moreover, authors have presented six different applications for how to use "meent" as a tool to generate data for ML as well as a solver for inverse problems.
weaknesses: I have a concern about the contribution of this paper. While having access to a user-friendly and differentiable software for Physics applications (e.g., EM simulator) is important and definitely helps research communities to accelerate their ideas, I am not completely convinced that this conference is a right place and fit for this paper. The main contribution of this paper is to introduce a python-based software, making the use of other developed tools easier for solving Physics applications. There are certainly other great venues where readers can take advantage of reading this paper and might be more fit to the audience.

---

confidence: 5
contribution: 1
initial_rating: 3
presentation: 3
questions: See 'Weakness' section.
soundness: 2
strengths: The paper is easy to follow.
summary: This paper introduces Meent, a software designed for electromagnetic field simulation using RCWA, and highlights potential applications in machine learning (ML) through Meent. The paper is clear and easy to follow. See below for detailed comments.
weaknesses: A key contribution of this work is the development of Meent. However, it is important to note that Meent is purely software based on numerical algorithms, and its development does not involve any ML techniques. The only connection to ML is the authors' use of Meent as a simulation backbone for dataset generation or device optimization.

This raises a significant concern: **there is no substantial contribution to the field of ML in this work**. If I may respectfully provide some examples that would align with ICLR's focus on ML, especially in the context of ML for photonics or AI for science: (1) The authors could propose a new neural operator (or architecture) that outperforms existing approaches such as FNO in predicting electromagnetic fields, such as [1,2]. (2) The paper could introduce a novel reinforcement learning algorithm that demonstrates better performance than current RL methods or use RL addressing a problem no one has done before, such as AlphaFold, AlphaGo. (3) The authors could apply ML techniques to design a more compact photonic device with much better performance that no one can. 

Unluckily, this work doesn't fall into any of the categories above. While I recognize the potential value of this work, the primary contribution is not centered on ML. Therefore, I must recommend rejection for ICLR and suggest the authors consider alternative venues that might be more suitable, such as IEEE JLT, APL Photonics, Optica, Optica Express, Nature Photonics, or Nanophotonics.

[1] Jiaqi Gu et al., 'Neurolight: A physics-agnostic neural operator enabling parametric photonic device simulation,' Neurips 2022.

[2] Pingchuan Ma et al., 'PIC2O-Sim: A Physics-Inspired Causality-Aware Dynamic Convolutional Neural Operator for Ultra-Fast Photonic Device FDTD Simulation,' Arxiv.

Meta-review: The paper presents MEENT a differentiable electromagnetic simulation software, that can be used to generate synthetic datasets for ML surrogates training or for using the simulator as part of the training in an ML pipeline. 

Reviewers unanimously agreed that this is a great scientific contribution but they questioned the scope of this submission to ICLR, and found it would be more  appropriate submission to scientific journals. Given that the machine learning contribution of this paper is rather not highlighted in the current draft I concur with the authors.

Decision: Reject