## Abstract

In inverse design, the design and background areas can be represented by different spatial resolutions; thus, adaptive meshes are more efficient than structured meshes. In this study, a second-order interpolation scheme is introduced to realize an inverse design process on an adaptive mesh. Experiment results show that the proposed scheme yields a 1.79-fold acceleration over that achieved using a structured mesh, aiding design time reduction or design area expansion. As the design area can be divided into multiple areas with different spatial resolutions, in future work, adaptive meshes can be combined with machine learning algorithms to further improve the inverse-design-process efficiency.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

The inverse design method, which is a combined simulation and optimization technique, has achieved remarkable results as regards the design of high-performance nanophotonic devices [1–6]. In this process, the processing speed and result accuracy are directly dependent on the simulation method. For nanophotonic device design using this approach, the finite-difference time-domain (FDTD) method or the finite difference frequency domain method (FDFD) can be employed to simulate the electromagnetic field [1,7]. Compared with the FDTD method, the FDFD method can be extremely fast and efficient, with lower memory consumption. In addition, the FDFD method exhibits superior numerical dispersion analysis and better ability to import oblique incidences and resolve sharp resonances [8–10]. As a result of these many advantages, the FDFD method is primarily employed for electromagnetic simulations in the context of inverse design.

In inverse design, optimization can be performed using the adjoint variable method, through which the gradient of an objective function for any arbitrary number of design parameters can be calculated via two electromagnetic simulations [11]. In this approach, a level-set function is usually used to design a more practical, discretely valued nanophotonic device [12–14]. When the level-set function is incorporated, adjoint optimization is usually referred to as topology optimization. Adjoint/topology optimization has been applied to linear, nonlinear, and active devices, and has exhibited significant ability as regards the design of nanophotonic devices at optical wavelengths [15–17].

In the inverse design process, the design time increases significantly with the design area; therefore, design of large-scale devices within a limited timeframe is difficult. Various studies have been conducted to reduce the design process time consumption for an enlarged design area, thereby accelerating the process. For example, researchers from Stanford University previously optimized the parameters of a perfectly matched layer (PML) [9] to improve the convergence of the iterative solver for an FDFD simulation. They then applied a 2.5-dimension FDFD method to the inverse design process to accelerate this process for a photonic crystal fiber [8]. They also employed the Schur complement method [18,19] to reduce the coefficient matrix size, and the conjugate orthogonal conjugate gradient method to reduce the iteration time of each step [20].

However, inverse design of large-scale devices remains challenging, as the design area is important and must be simulated at high spatial resolution. The background area containing waveguides and the PML can be simulated at low spatial resolution. As different spatial resolutions can be applied to different areas, an adaptive mesh is more efficient than a structured mesh.

Non-uniform gridding, octree grids, and subgridding/adaptive meshes, hereafter collectively referred to as “adaptive meshes,” are usually used to accelerate electromagnetic simulation. Several studies considering the FDFD method and adaptive meshes have been conducted to date [21–24]. For example, Kulas and Mrozowski proposed a novel three-dimensional subgridding scheme applicable to the finite-difference technique in the time and frequency domains [21]. When the refinement factor was set to 3, the field values on the coarse mesh overlapped with some of the field values on the fine mesh; hence, mapping of the field values from the coarse to fine mesh was no longer necessary. Therefore, the reflection between coarse and fine meshes could be reduced. Subsequently, Podwalski, Kowalczyk, and Mrozowski used a subgridding scheme to improve the efficiency of FDFD simulations for microwave and radiofrequency (RF) devices [22]. Further, Fotyga et al. presented a novel local and nested (multilevel) macromodel that could further improve the FDFD solver speed [23], while Luan et al. proposed a novel FDFD method for RF device simulation based on octree grids, combining a linear interpolation scheme and a second-order interpolation scheme to reduce reflection [24]. In those works, the simulation was accelerated through introduction of an adaptive mesh. However, the interpolation scheme on the boundary between the fine and coarse meshes introduced numerical dispersion and artificial reflections to the simulation results. In a small number of studies, the FDFD method has also been applied to an adaptive mesh for optical device simulation [23]. To the best of the authors’ knowledge, however, incorporation of an adaptive mesh with the FDFD method has not yet been applied to inverse design.

To apply FDFD simulation to an adaptive mesh, an interpolation scheme between the coarse and fine meshes must be implemented on their boundary. Linear, high-order, and moving least-squares meshless interpolation have been used in previous studies [24–27]. However, the moving least-squares meshless interpolation, in which weights are applied to different simulation types, is not suitable for inverse design. Further, the simulation results obtained using linear interpolation are not sufficiently accurate. High-order interpolation yields greater accuracy; however, the simulation speed is also reduced.

In this study, to accelerate the inverse design process for large-scale devices, an adaptive mesh is introduced to an FDFD system and optimization is realized based on the adjoint variable method. As regards the interpolation scheme, experiments conducted in this study show that the second-order interpolation scheme is sufficiently accurate and, hence, a second-order interpolation scheme is employed. To implement the inverse design process, a linear coefficient matrix based on the FDFD method must be constructed. A less efficient approach to obtaining the coefficient matrix involves addition of several specific points to the coefficient matrix on the boundary between the coarse and fine meshes. However, the present study provides a general framework that allows effective construction of a coefficient matrix on an adaptive mesh. For a 1 × 2 beam splitter at 1550 nm as an inverse design target, a numerical experiment is performed and the FDFD simulation accuracy is verified through comparison with the FDTD module of the Lumerical commercial software package. Furthermore, to assess the speed of the design process, an FDFD simulation with a structured mesh is conducted for comparison.

The remainder of this paper is organized as follows. Section 2.1 presents the basic theory of the adjoint variable method and Section 2.2 illustrates the interpolation scheme. In Section 2.3, the general framework for constructing the FDFD coefficient matrix is presented. Section 3 analyzes the results of the inverse design numerical experiment and the limitations of the inverse design process based on the adaptive mesh. Finally, Section 4 concludes this study.

## 2. Basic theory and principles

#### 2.1 FDFD method and adjoint variable method

Before discussing acceleration of the inverse design process, we first briefly describe the FDFD and adjoint variable methods, which are used to realize the inverse design process.

The FDFD method based on the frequency-domain Maxwell equations is the simulation method typically applied for performance analysis of nanophotonic devices in inverse design. For this method, simulation is performed by solving a linear equation with form

where*Ce*and

*Ch*represent the curl-operator matrices on the electric and magnetic fields, respectively;

*ε*and

*μ*are the diagonal matrices of the permittivity and permeability, respectively;

*J*is the electric current source density;

*ω*is the frequency; and

*E*is the electric-field distribution.

Equation (1) can be simplified as

where*A*=

*Ce*×

*μ*

^{−1}×

*Ch*–

*ω*and

^{2}ε*b*= –

*iωJ*. By solving Eq. (2),

*E*can be calculated based on the permittivity and permeability.

In the inverse design process, the permittivity distribution with respect to the electric field can also be determined. This process can usually be realized using a convex optimization algorithm, with the adjoint variable method being used to calculate the optimization gradient.

For the latter, the figure of merit, described by a scalar objective function *L*, must be maximized for the design variables *φ*. The adjoint variable method then computes the sensitivity. Here,

Typically, ∂*L*/*∂E* can be solved analytically. In Eq. (3), *φ* represents the material distribution in the design region. To determine ∂*E*/*∂φ*, Eq. (2) can be derived as follows:

Hence, the following can be obtained:

Substituting Eq. (5) into Eq. (3) yields

*E’*is usually referred to as the adjoint field. In Eq. (6), ∂

*A*/

*∂φ*can be obtained analytically, and the vector

*E*can be obtained by solving the linear system

*AE = b*. A linear system must be constructed to calculate

*E’*, as follows:

In each optimization step, the design variables must be updated from *φ _{old}* to

*φ*as follows:

_{new}*α*is the learning rate. By using Eq. (8) to repeatedly update

*φ*, the final, required design variable can be obtained.

#### 2.2 Interpolation scheme

An adaptive mesh is usually a combination of several meshes having different spatial resolutions. In this study, only two meshes are used to construct the adaptive mesh, having low and high spatial resolutions and referred to as the coarse and fine meshes, respectively. Distributions of structured and adaptive meshes are illustrated in Fig. 1. Hereafter, the Yee cell sizes of the coarse and fine meshes are labeled *dc* and *df*, respectively.

In inverse design, the designed device is usually a thin slice with fewer Yee cells in the *z*-direction. Thus, it is not necessary to set different spatial resolutions for the coarse and fine meshes in the *z*-direction, and only the interpolation scheme for the other two directions is discussed in this section. To expand the application potential of the proposed approach, the interpolation scheme in three directions is presented in Appendix 1. To perform FDFD simulation on an adaptive mesh, field values are obtained from the fine mesh with the field values from the coarse mesh, and field values are obtained from the coarse mesh with the field values from the fine mesh.

First, we discuss acquisition of field values from the fine mesh to the coarse mesh. In the FDFD simulation, the field component is mapped onto a staggered mesh, which is usually referred to as the Yee mesh. To apply an adaptive mesh to an FDFD system, the coarse and fine meshes must obey the Yee theory, which states that different field components in a single Yee cell have different positions. Therefore, if one entire Yee cell is split into multiple sub-Yee cells, an interpolation scheme is required to directly obtain the field values from the sub-Yee cells for the single entire Yee cell. In the context of the present work, if the entire Yee cell is regarded as the Yee cell on the coarse mesh and the multiple sub-Yee cells are taken as the Yee cells on the fine mesh, the field values must be interpolated from the fine mesh to the coarse mesh. However, previous research has shown that each field value on the coarse mesh can be spatially aligned with one of the field values on the fine mesh, as shown in Fig. 2, if the refinement factor (*dc/df*) is set to 3 [21]. Then, the field values on the coarse mesh can be obtained by determining the spatially aligned field values on the fine mesh. Based on this theory, the interpolation scheme can be simplified and the reflection on the boundary between meshes can be reduced.

Next, interpolation to the fine mesh with field values from the coarse mesh is considered. To simplify the explanation, only the interpolation scheme for the *Hz* field component is discussed, where *H* represents the magnetic field. Using a similar approach, interpolation functions can be constructed for the *Ex, Ey, Ez, Hx*, and *Hy* components. The interpolation function for *Hz* in the Yee cell can be defined as

*a, b, c, d*, and

*e*are unknown values. The spatial arrangement of

*Hz*(

*x*,

_{0}*y*– 1),

_{0}*Hz*(

*x*– 1,

_{0}*y*),

_{0}*Hz*(

*x*),

_{0},y_{0}*Hz*(

*x*+ 1,

_{0}*y*), and

_{0}*Hz*(

*x*,

_{0}*y*+ 1) is shown in Fig. 3, where

_{0}*Hz*(

*x*) is the central point and

_{0},y_{0}*x*and

_{0}*y*are its coordinates.

_{0}To calculate *a, b, c, d,* and *e*, it is assumed that *Hz* can be expressed as a second-order function of the spatial coordinates such that

By substituting the field values and spatial positions of the five *Hz* components into Eq. (10), the following can be obtained:

Then, solving for *a’, b’, c’, d’*, and *e*’ and substituting the solutions into Eq. (10), the interpolation function can be written as

Compared with Eq. (9), *a, b, c, d,* and *e* can be expressed as

Based on the interpolation function for each component, several sub-interpolation matrices can be constructed for each component, here labeled *D _{Ex}, D_{Ey}, D_{Ez}, D_{Hx}, D_{Hy},* and

*D*. Details of the sub-interpolation matrix construction are provided in Appendix 2. By combining

_{Hz}*D*, and

_{Ex}, D_{Ey}, D_{Ez}, D_{Hx}, D_{Hy}*D*, interpolation matrices labeled

_{Hz}*D*,

_{inp}*D*, and

_{E}*D*can be obtained, such that

_{H}*E*and

_{fine}*H*represent the electric and magnetic fields on the fine mesh,

_{fine}*E*,

_{coarse}*H*represent the electric and magnetic fields on the coarse mesh and

_{coarse}*D*represents the interpolation matrix from the electric and magnetic fields on the coarse mesh to the electric and magnetic fields on the fine mesh.

_{inp}It must be emphasized that, in the present study, all the field values are mapped from the coarse mesh to the fine mesh to simplify the implementation. However, only parts of *D _{inp}* are used on the boundary to construct the FDFD system, such that the interpolation scheme is only applied on the boundary.

#### 2.3 Implementation

In this section, we discuss construction of an FDFD system based on the interpolation scheme. To simplify the construction, two FDFD systems are built: one on the coarse mesh and the other on the fine mesh. We define *Ce _{coarse}* and

*Ch*, and

_{coarse}*Ce*and

_{fine}*Ch*as the difference matrices of the two FDFD systems, respectively. To simplify the discussion, the region of the design area in which the Yee cells on the fine mesh spatially align with the Yee cells on the coarse mesh is defined as the overlapping area. The regions of the design area containing the other Yee cells that are not spatially aligned are defined as the non-overlapping areas. The spatial distribution of the overlapping areas, non-overlapping areas, and background area of the fine mesh is shown in Fig. 4.

_{fine}Based on this spatial distribution, *Ce _{fine}, Ce_{coarse}*, and

*D*can be expressed as

_{E}*Cec*and

_{background}*Cec*represent the difference matrices in the background and design areas on the coarse mesh, respectively; and

_{design}*Cef*,

_{overlap}*Cef*, and

_{non-overlap}*Cef*represent the difference matrices in the overlapping, non-overlapping, and background areas on the fine mesh, respectively. The remaining matrices in Eq. (16) are coupling matrices.

_{background}To build the FDFD system on an adaptive mesh, the following matrices should be constructed: two difference matrices representing the curl operator on the electric and magnetic fields, *Ce _{combine}* and

*Ch*, respectively; two diagonal matrices representing the permittivity and permeability,

_{combine}*ε*and

_{combine}*μ*, respectively; and one vector representing the electric current source density,

_{combine}*J*. Of all the required matrices and vectors, construction of

_{combine}*Ce*and

_{combine}*Ch*is particularly challenging.

_{combine}First, we consider construction of the difference matrix. Only construction of *Ce _{combine}* is discussed in detail, as

*Ch*can be constructed using a similar method. To simplify the construction,

_{combine}*Ce*is obtained by combining

_{combine}*Ce*and

_{coarse}, Ce_{fine},*D*. This combination can be performed because the design area has higher spatial resolution than the background area in the inverse design process. Therefore, only the portions of the difference matrix related to the background area on the coarse mesh and the design area on the fine mesh are considered.

_{E}For detailed discussion of the difference matrix construction, we briefly review the interpolation scheme. The *Hz* distribution shown in Fig. 5 is used as an example. On the coarse mesh, the three field components marked with hollow circles are required in order to realize the forward difference process. As shown in Fig. 5(a), the hollow-circle field components overlap some of the field components on the fine mesh. In this case, the hollow-circle field components can be obtained by selecting the overlapping components on the fine mesh. On the fine mesh, the nine field components marked with hollow circles are required in order to realize the backward difference process. As shown in Fig. 5(b), the hollow-circle field components can be obtained by implementing an interpolation scheme.

First, a scheme for interpolation from the fine to coarse mesh is considered. Because the Yee cells in the overlapping area on the fine mesh are spatially aligned with those in the design area on the coarse mesh, the coupling matrix *Cec _{d}*

_{→b}can be used directly. Combining

*Cec*,

_{background}*Cec*

_{d}_{→b}, and

*Ce*yields the following:

_{fine}Next, interpolation from a coarse to fine mesh is considered. To introduce the interpolation scheme, *D _{E}* is incorporated into Ce

_{add}_{1}to yield

*I*is the unit matrix. The rank of Ce

_{add}_{2}in Eq. (18) exceeds the size of Ce

_{add}_{2}, i.e.,

*n*. Therefore, Eq. (18) can be eliminated, and the required matrix

*Ce*can be obtained, where

_{combine}Further, *ε _{combine}*,

*μ*, and

_{combine}*J*can be constructed by extracting and combining components of the permittivity, permeability, and electric current source from the FDFD systems on the coarse and fine meshes having the spatial distribution shown in Fig. 4.

_{combine}Using the matrices *Ce _{combine}, Ch_{combine}, ε_{combine}*, and

*μ*, the coefficient matrix for the FDFD system on the adaptive mesh can be constructed as

_{combine}The field value of the adaptive mesh can be obtained by solving a linear system in the form

Following introduction of the interpolation matrix to the FDFD system, the coefficient matrix is no longer the Hermitian matrix. However, for the experiment conducted in this study, the biconjugate gradient (BICG), which is usually applied to Hermitian system solution, could provide an appropriate solution for the FDFD system on the adaptive mesh. Therefore, the BICG solver was chosen for application to the FDFD system.

Finally, as shown in Eq. (20), *ε _{combine}*, which incorporates the design variable

*φ*, is on the diagonal line of

*A*. Therefore,

_{combine}*A*can be directly introduced to the adjoint variable method.

_{combine}## 3. Numerical experiment

This section describes the numerical experiment and discusses the results. All computations were performed on a computer with two graphics processing units (NVIDIA Tesla P100 16G) and one central processing unit (Intel Xeon E5-2630). The software was based on the g++ platform with the CUDA package.

Because inverse design acceleration is more advantageous for large-scale problems, only an FDFD system containing more than 180 × 180 × 50 Yee cells was considered. Thus, the Yee cell dimensions were set to 75 and 25 nm for the coarse and fine meshes, respectively, where the dimensions were identical in both the *x*- and *y*-directions; the *z*-direction dimension for both meshes was set to 25 nm. On the structured mesh, the Yee cell dimensions in both the *x*- and *y*-directions were set to 25 nm. Therefore, all simulation areas for the inverse design process exceeded 4.5 μm × 4.5 μm × 1.25 μm.

To illustrate the advantage of the adaptive mesh, a 1 × 2 beam splitter was used as the experimental target. The initial scheme of the beam splitter is illustrated in Fig. 6. The waveguide width and Si layer thickness were 500 and 220 nm, respectively. In the inverse design process, the Si and SiO_{2} refractive indices were set to 3.5 and 1.5, respectively. To normalize the comparison, the distance between the design area and the boundary of the simulation area was set to 1.125 μm.

#### 3.1 Acceleration

In this study, the simulation area dimensions in both the *x*- and *y*-directions were used as the length fraction. Before analyzing the acceleration, the size of the coefficient matrix for the adaptive-mesh FDFD system and was compared with that for the structured-mesh FDFD system under different length fractions, as shown in Fig. 7. The coefficient matrix size for the adaptive-mesh FDFD system was far smaller than that of the structured-mesh FDFD system. As the length fraction increased, the difference in coefficient matrix size between the structured- and adaptive-mesh FDFD systems increased. For the BICG solver, the speed of each step depends on the coefficient matrix size; thus, the above result implies that, for each step of the BICG solver, the calculation for the adaptive-mesh FDFD system can be faster than that for the structured mesh.

For the BICG solver, the convergence also directly affects the computational time. Here, the number of iterative steps required for convergence was taken to represent convergence. A comparison of the convergence for different length fractions is shown in Fig. 8, and reveals that the iterative solver convergence on the adaptive mesh was poorer than that on the structured mesh for a simulation area smaller than 6.1 μm × 6.1 μm × 1.25 μm. As the simulation area increased, the iterative-solver convergence degraded for both cases. However, for a simulation area exceeding 6.1 μm × 6.1 μm × 1.25 μm, the iterative-solver convergence on the adaptive mesh exceeded that on the structured mesh.

Each iterative-solver step was conducted more quickly on the adaptive mesh than on the structured mesh; thus, the BICG solver and, hence, the inverse design process, can be implemented more quickly on an adaptive mesh for large-scale problems.

To analyze the performance of the adaptive-mesh FDFD simulation for a large-scale problem, the total computational times and matrix-vector-multiplication times of both the adaptive- and structured-mesh FDFD simulations were compared for different length fractions, as shown in Fig. 9. Even though the simulation area was quite large, the adaptive-mesh FDFD simulation was faster. Furthermore, the comparative results for the total computational time and vector multiplication time had a similar form, which confirms that the iterative-solver time consumption depends on the vector multiplication time. As the vector multiplication time depends on the coefficient matrix size, the time-consumption comparison shown in Fig. 9 matches the size comparison shown in Fig. 7.

#### 3.2 Accuracy verification

To verify the accuracy of the simulation results, the FDTD module of the Lumerical commercial optical simulation software was used to simulate a nanophotonic device, which was designed with an inverse design process on an adaptive mesh. As shown in Fig. 10, the simulation result obtained from the adaptive-mesh FDFD method matched that of the Lumerical FDTD module. However, for several steps, the transmission obtained with the adaptive mesh differed slightly from the Lumerical result. This phenomenon arose because the interpolation scheme provided an approximate result only, and because reflection was introduced on the boundary between the coarse and fine meshes. The design and simulation results are presented in Fig. 11.

Usually, improved FDFD simulation accuracy corresponds to reduced speed, and vice versa. Here, the accuracy and acceleration were compared to appropriately analyze the FDFD simulation efficiency. Table 1 compares the transmission in the top and bottom ports for the designed 1 × 2 beam splitters and the total time consumption of the inverse design process with different mesh resolutions. Note that the transmission was calculated using the Lumerical FDTD module. For condition normalization, the target coupling efficiency at both output ports was set to 0.49–0.51. The maximum optimization step was defined as 30. The inverse design process with the adaptive mesh achieved a 1.79-fold acceleration over the structured-mesh case (25 nm). In addition, for a mesh resolution of approximately 40 nm, the time consumption (i.e., speed) of the inverse design process was similar to that of the adaptive mesh. However, the performance of the nanophotonic device obtained through the inverse design process on the structured mesh was lower than that obtained through the inverse design process on the adaptive mesh.

#### 3.3 Discussion

Inverse design on an adaptive mesh can accelerate the design process and increase the design area. However, the adaptive mesh introduces a small reflection at the boundary between the coarse and fine meshes, which generates a difference in transmission. To verify the simulation results, transmissions in the top and bottom ports given by the adaptive-mesh FDFD simulation and Lumerical FDTD simulation were compared for different frequencies. As shown in Fig. 12, both the resonant frequency and transmission differed slightly between the two simulations. This phenomenon arose because the interpolation scheme introduced a reflection at the boundary between the coarse and fine meshes. Further, the simulations with different spatial resolutions may have had different resonant frequencies. This problem can be solved by employing a high-order interpolation scheme.

For inverse design, the designed device performance degrades as the spatial resolution decreases. To analyze the limitations of the adaptive mesh with regard to the spatial resolution, the designed-device transmissions for different spatial resolutions were compared. As detailed in Table 2, when the coarse mesh resolution was 105 nm, the top- and bottom-port transmissions were 0.345 and 0.347, respectively; these values are low for the design of a nanophotonic device. However, for a 90-nm coarse mesh resolution, the top- and bottom-port transmissions were acceptable, at 0.413 and 0.414, respectively. Therefore, we suggest that the coarse mesh resolution must exceed 90 nm for effective implementation of the inverse design process.

## 4. Conclusion

In this study, the inverse design process computation was accelerated using an FDFD simulation on an adaptive mesh. Experiment results showed that 1.79-fold acceleration was achieved by the adaptive-mesh FDFD simulation over the structured-mesh case for an FDFD system simulation area of 4.5 μm × 4.5 μm × 1.25 μm. For verification, the simulation results were compared with those given by the Lumerical FDTD module, with good agreement being obtained. Further, a 2.25 μm × 2.25 μm beam splitter at 1550 nm was obtained through inverse design using an FDFD simulation on an adaptive mesh, with 0.458-dB insertion loss and low time consumption.

However, FDFD simulation on an adaptive mesh usually yields a system error at the boundary between the fine and coarse meshes, which in turn causes artificial reflection and numerical dispersion in the inverse design process. Therefore, in the context of inverse design, we believe that an adaptive mesh is more suitable for initial high-precision designs of nanophotonic devices rather than final designs.

In addition to acceleration and simulation area expansion, FDFD simulation on an adaptive mesh can be used to develop machine-learning-based optimization methods to identify the more important parts of a nanophotonic device. Hence, components requiring individual simulations with high and low resolutions could be determined. The acceleration achieved in this work could then be further improved by optimizing the device on multilevel meshes. In addition, inverse design on an adaptive mesh could be used to determine the most important characteristic in the design area, for a high-resolution, regularly shaped area. Indeed, a regular shape could be adopted as a characteristic for nanophotonic devices. The authors intend to develop these techniques in future work.

## Appendix 1. Interpolation scheme in three directions

The interpolation scheme in three directions differs from that in two directions. First, instead of 3 × 3 field values on the fine mesh, there are 3 × 3 × 3 field values. Second, the interpolation function has three inputs. The interpolation function can be written as

The spatial distribution of *Hz*(*x, y, z*) is shown in Fig. 13. To calculate *a, b, c, d, e, f*, and *g*, spatial coordinates are used to represent *Hz*, where

The interpolation function must satisfy

From Eq. (25), *a, b, c, d, e, f*, and *g* can be calculated such that

## Appendix 2. Construction of sub-interpolation matrix

Based on Eqs. (9) and (14), the interpolation matrix for a Yee cell on a coarse mesh can be constructed as

*x*

_{1}, …

*x*

_{9}and

*y*

_{1}, …

*y*

_{9}are the spatial coordinates of each Yee cell on the fine mesh within the Yee cell on the coarse mesh.

The interpolation matrix *DH _{zone_cell}* is the same for each Yee cell on a coarse mesh. Therefore, the sub-interpolation matrix

*DHz*can be constructed through point-by-point setting of

*DH*.

_{zone_cell}To simplify the discussion, the simulation area shown in Fig. 14 is taken as an example. Then, *DHz* can be expressed as

To accelerate the construction, the indexes of each Yee cell on the fine and coarse meshes can be constructed. Faster construction of *DHz* can be achieved based on these indexes and *DHz _{one_cell}*. Further, the sub-interpolation matrices

*DEx, DEy, DEz, DHx,*and

*DHy*can be constructed based on a similar method.

## Funding

National Natural Science Foundation of China (61935003).

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

**1. **A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nat. Photonics **9**(6), 374–377 (2015). [CrossRef]

**2. **Z. Ye, J. Qiu, C. Meng, L. Zheng, Z. Dong, and J. Wu, “Inverse design of a SOI T-junction polarization Beamsplitter,” J. Phys.: Conf. Ser. **844**, 012009 (2017). [CrossRef]

**3. **N. V. Sapra, D. Vercruysse, L. Su, K. Y. Yang, J. Skarda, A. Y. Piggott, and J. Vuckovic, “Inverse design and demonstration of broadband grating couplers,” IEEE J. Select. Topics Quantum Electron. **25**(3), 1–7 (2019). [CrossRef]

**4. **N. Mohammadi Estakhri, B. Edwards, and N. Engheta, “Inverse-designed metastructures that solve equations,” Science **363**(6433), 1333–1338 (2019). [CrossRef]

**5. **B. Kailkhura, B. Gallagher, S. Kim, A. Hiszpanski, and T. Y. Han, “Reliable and explainable machine-learning methods for accelerated material discovery,” npj Comput. Mater. **5**(1), 108–109 (2019). [CrossRef]

**6. **L. Su, R. Trivedi, N. V. Sapra, A. Y. Piggott, D. Vercruysse, and J. Vučković, “Fully-automated optimization of grating couplers,” Opt. Express **26**(4), 4023–4034 (2018). [CrossRef]

**7. **B. Shen, P. Wang, R. Polson, and R. Menon, “An integrated-nanophotonics polarization beamsplitter with 2.4 × 2.4 μm2 footprint,” Nat. Photonics **9**(6), 378–382 (2015). [CrossRef]

**8. **J. Lu, S. Boyd, and J. Vučković, “Inverse design of a three-dimensional nanophotonic resonator,” Opt. Express **19**(11), 10563–10570 (2011). [CrossRef]

**9. **W. Shin and S. Fan, “Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers,” J. Comput. Phys. **231**(8), 3406–3431 (2012). [CrossRef]

**10. **W. Shin and S. Fan, “Accelerated solution of the frequency-domain Maxwell’s equations by engineering the eigenvalue distribution of the operator,” Opt. Express **21**(19), 22578–22595 (2013). [CrossRef]

**11. **M. G. Damigos and E. M. Papoutsis, “kgIachagias, Giannakoglou K C, adjoint variable-based shape optimization with bounding surface constraints,” Int. J. Numer. Meth. Fluids (2020).

**12. **N. Lebbe, ““Contribution in topological optimization and application to nanophotonics,” Analysis of PDEs,” NNT 2019GREAM047, tel-02518286 (2019).

**13. **R. E. Christiansen, F. Wang, and O. Sigmund, “Topological insulators by topology optimization,” Phys. Rev. Lett. **122**(23), 234502 (2019). [CrossRef]

**14. **S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, “Inverse design in nanophotonics,” Nat. Photonics **12**(11), 659–670 (2018). [CrossRef]

**15. **T. W. Hughes, M. Minkov, I. A. D. Williamson, and S. Fan, “Adjoint method and inverse design for nonlinear nanophotonic devices,” ACS Photonics **5**(12), 4781–4787 (2018). [CrossRef]

**16. **J. Wang, Y. Shi, T. Hughes, Z. Zhao, and S. Fan, “Adjoint-based optimization of active nanophotonic devices,” Opt. Express **26**(3), 3236–3248 (2018). [CrossRef]

**17. **L. Su, A. Y. Piggott, N. V. Sapra, J. Petykiewicz, and J. Vučković, “Inverse design and demonstration of a compact on-chip narrowband three-channel wavelength demultiplexer,” ACS Photonics **5**(2), 301–305 (2018). [CrossRef]

**18. **N. Z. Zhao, S. Boutami, and S. Fan, “Accelerating adjoint variable method based photonic optimization with Schur complement domain decomposition,” Opt. Express **27**(15), 20711–20719 (2019). [CrossRef]

**19. **N. Zhao, S. Verweij, W. Shin, and S. Fan, “Accelerating convergence of an iterative solution of finite difference frequency domain problems via Schur complement domain decomposition,” Opt. Express **26**(13), 16925–16939 (2018). [CrossRef]

**20. **L. Su, D. Vercruysse, J. Skarda, N. V. Sapra, J. A. Petykiewicz, and J. Vučković, “Nanophotonic inverse design with SPINS: Software architecture and practical considerations,” Appl. Phys. Rev. **7**(1), 011407 (2020). [CrossRef]

**21. **L. Kulas and M. Mrozowski, “Low-reflection subgridding,” IEEE Trans. Microwave Theory Tech. **53**(5), 1587–1592 (2005). [CrossRef]

**22. **J. Podwalski, P. Kowalczyk, and M. Mrozowski, “Efficient multiscale finite difference frequency domain analysis using multiple macromodels with compressed boundaries,” Prog. Electromagn. Res. **126**(463), 479 (2012). [CrossRef]

**23. **G. Fotyga, P. Kowalczyk, L. Kulas, K. Nyka, J. Podwalski, and M. Mrozowski, “Reduced order models in computational electromagnetics (in memory of Ruediger Vahldieck),” ASIA PAC. INTERNATIONAL SYMPOSIUM ON ELECTROMAGNETIC COMPATIBILITY (APEMC), 2012.

**24. **H. Luan, Y. Geng, Y. Yu, and S. Guan, “Three-dimensional transient electromagnetic numerical simulation using FDFD based on octree grids,” IEEE Access **7**, 161052–161063 (2019). [CrossRef]

**25. **X. Zhu, Y. Shu, C. Luo, F. Huo, and W. Zhu, “Transient electromagnetic voltage imaging of dense UXO-like targets based on improved mathematical morphology,” IEEE Access **8**, 150341–150349 (2020). [CrossRef]

**26. **R. Goldade, Y. Wang, M. Aanjaneya, and C. Batty, “An adaptive variational finite difference framework for efficient symmetric octree viscosity,” ACM Trans. Graph. **38**(4), 1–14 (2019). [CrossRef]

**27. **A. Fedeli, C. Montecucco, and G. L. Gragnani, “Open-source software for electromagnetic scattering simulation: The case of antenna design,” Electronics **8**(12), 1506 (2019). [CrossRef]