Eigendecomposition-free inverse design of meta-optics devices

Yi Huang; Ziwei Zhu; Yunxi Dong; Hong Tang; Bowen Zheng; Viktor A. Podolskiy; Hualiang Zhang

doi:10.1364/OE.514347

1. Introduction

Metalenses and meta-atoms represent a paradigm shift in optical design, yielding devices with significantly reduced form factors and enhanced performance. Methods of inverse design driven by specific design goals have been demonstrated to be effective in optimizing high-performance meta-structures [1,2]. Despite these advancements, the inverse design of meta-atoms poses significant challenges primarily because of the requirements for high-resolution nanoscale simulation. Rigorous coupled-wave analysis (RCWA), initially introduced in [3], has gained prominence for simulating and optimizing meta-atoms as described in [4,5]. Although it is much faster compared to other fully discretized methods such as finite-difference time-domain (FDTD) [6], difficulties in parallelization curtail its efficiency because of the iterative eigendecomposition process [7], which hampers the ability to fully leverage modern computational resources such as graphics processing units (GPUs) to enhance performance.

The rigorous diffraction interface theory (R-DIT) [8,9] circumvents the eigendecomposition procedure in RCWA and, therefore, can be easily parallelized. It utilizes the Taylor series to approximate the matrix exponential and sidesteps the eigendecomposition process completely. Low-order R-DIT yields accurate results when the thickness is small compared to the wavelength and dramatically speeds up the computation. [8] By incorporating multiple terms in the Taylor expansion, R-DIT can adequately describe the photonic response of increasingly thick structures, as in the regime used in this work to design photonic devices with wavelength-scale thick meta-atoms.

Building on these insights, we introduce a fully-differentiable inverse design framework based on the R-DIT algorithm. There are two major contributions in this work: First, by reformulating the original R-DIT algorithm and implementing it with PyTorch [10], we achieve up to $16.2\times$ speedup over the traditional inverse design method based on RCWA on GPU. PyTorch also enables the computation of gradients of meta-atom performance featuring arbitrary design parameters and is free from extensive mathematical derivations. Secondly, we incorporate a novel neural network representation [11], which can implicitly encode the shapes of meta-atoms for topological optimization, providing a flexible and generalized approach to their design. We aim to extend the capabilities of our proposed framework to engineer innovative metasurfaces with multifunctional capabilities. These include but are not limited to actively reconfigurable or tunable beam deflectors [12,13], and spatiotemporal modulators [14], etc. The inverse design of such meta-atoms presents a formidable challenge, necessitating the creation of devices that offer high degrees of design freedom while fulfilling the desired performance on modulated materials that exhibit a high refractive index contrast. A significant challenge in this endeavor is devising a method to parameterize meta-atom shapes without restricting them to a limited set of shape candidates.

In this paper, we start by detailing our improved formalism of the R-DIT algorithm in Section 2. In Section 3, we exemplify an inverse design workflow based on our proposed TorchRDIT with demonstrations of two use cases, including parameter-constrained and free-form meta-atoms. Finally, Section 4. presents the benchmarks of the proposed inverse design workflow and evaluates its performance on both CPU and GPU platforms.

The source code of TorchRDIT has been made publicly available following GNU General Public License version 3 (GPLv3) [15].

2. Theoretical formulation of rigorous diffraction interface theory

In this section, we derive our improved formalism employed in TorchRDIT, as well as our previous work [9], based on the scattering matrices method [16] with buffered layers for the compatibility with existing semi-analytical methods.

2.1 Definition of scattering matrix

We start by defining the scattering matrix of the interested layer $i$. In this work, we denote $\mathbf {a}$ as incident wave coefficients of the light injecting from surrounding layers into layer $i$, such as $\mathbf {a}_{i-1}$ and $\mathbf {a}_{i+1}$ shown in Fig. 1(red arrows); similarly, $\mathbf {b}_{i-1}$ and $\mathbf {b}_{i+1}$ (green arrows) serve as the light transmitted out of layer $i$. The scattering matrix of layer $i$, describing its scattered properties, is bounded by these coefficients and can be formulated as:

(1)$$\begin{bmatrix} \mathbf{b}_{i-1} \\ \mathbf{b}_{i+1} \end{bmatrix} = \mathbf{S}^i \begin{bmatrix} \mathbf{a}_{i-1} \\ \mathbf{a}_{i+1} \end{bmatrix},\,\text{where}\, \mathbf{S}^i = \begin{bmatrix} \mathbf{S}^i_{11} & \mathbf{S}^i_{12} \\ \mathbf{S}^i_{21} & \mathbf{S}^i_{22} \end{bmatrix}.$$

Fig. 1. Definition of the scattering matrix of layered photonic structures. The scattering matrix of layer $i$ consistent along the $z$-direction is bounded by the incident and transmitted waves from its surrounding layers.

Download Full Size | PDF

$\mathbf {S}^i_{11}$ means the energy of modes reflected to layer $i-1$ and $\mathbf {S}^i_{21}$ means the transmitted energy from layer $i-1$ to layer $i+1$. Besides, for the convenience of the derivation of fields inside layer $i$, we define $\mathbf {a}_{i}$ as the light at the left interfaces $z=-h/2$ propagating to the center of layer $i$, and $\mathbf {b}_{i}$ as the scattered waves at the same interface. It is noteworthy that the choices of positions and directions of $\mathbf {a}_{i}$ and $\mathbf {b}_{i}$ do not affect the essence of the scattering matrix, and all field coefficients at any position inside layer $i$ can be effortlessly derived from the left interface.

To solve the scattering matrix $\mathbf {S}^i$ in (1), we start from the $P-Q$ form of Maxwell equations in the Fourier domain. Using the same coordinates shown in Fig. 1, where $x-y$ plane is the cross-section of the given photonic structures and $\tilde {z}=k_0z$ is the normalized longitudinal coordinate ($k_0$ is the wave number in the free space), the Maxwell equation can be expressed as:

(2)$$\frac{\mathrm{d}}{\mathrm{d}\tilde{z}} \begin{bmatrix} \vec{\mathcal{E}}_{x,y} \\ \vec{\mathcal{H}}_{x,y} \end{bmatrix} = \begin{bmatrix} 0 & \hat{\mathbb{P}} \\ \hat{\mathbb{Q}} & 0 \end{bmatrix} \begin{bmatrix} \vec{\mathcal{E}}_{x,y} \\ \vec{\mathcal{H}}_{x,y} \end{bmatrix},\quad \vec{\mathcal{E}}_{x,y} = \begin{bmatrix} \vec{\mathcal{E}}_x \\ \vec{\mathcal{E}}_y \end{bmatrix},\quad \vec{\mathcal{H}}_{x,y} = \begin{bmatrix} \vec{\mathcal{H}}_x \\ \vec{\mathcal{H}}_y \end{bmatrix},$$

where $\vec {\mathcal {E}}_{x,y}$ is the Fourier components of the tangential electric field $\vec {E}_{x,y}$; similarly, $\vec {\mathcal {H}}_{x,y}$ is the Fourier components of the normalized tangential magnetic field $\vec {\tilde {H}}_{x,y}=-j\eta _0\vec {H}_{x,y}$ ($\eta _0$ is the free space impedance). Note that the negative sign convention $e^{-jkz}$ is used throughout this work. The matrices $\hat {\mathbb {P}}$ and $\hat {\mathbb {Q}}$ are block matrices formed by the distributions of permittivity and permeability at the cross-section of layer $i$, formulated as 3) and (4) below:

(3)$$\hat{\mathbb{P}} = \begin{bmatrix} \mathbf{K}_x\left[\kern-0.15em\left[{\varepsilon_r} \right]\kern-0.15em\right]^{{-}1}\mathbf{K}_y & \left[\kern-0.15em\left[{\mu_r} \right]\kern-0.15em\right] - \mathbf{K}_x\left[\kern-0.15em\left[{\varepsilon_r} \right]\kern-0.15em\right]^{{-}1}\mathbf{K}_x \\ \mathbf{K}_y\left[\kern-0.15em\left[{\varepsilon_r} \right]\kern-0.15em\right]^{{-}1}\mathbf{K}_y - \left[\kern-0.15em\left[{\mu_r} \right]\kern-0.15em\right] & - \mathbf{K}_y\left[\kern-0.15em\left[{\varepsilon_r} \right]\kern-0.15em\right]^{{-}1}\mathbf{K}_x \end{bmatrix},$$

(4)$$\hat{\mathbb{Q}} = \begin{bmatrix} \mathbf{K}_x\left[\kern-0.15em\left[{\mu_r} \right]\kern-0.15em\right]^{{-}1}\mathbf{K}_y & \left[\kern-0.15em\left[{\varepsilon_r} \right]\kern-0.15em\right] - \mathbf{K}_x\left[\kern-0.15em\left[{\mu_r} \right]\kern-0.15em\right]^{{-}1}\mathbf{K}_x \\ \mathbf{K}_y\left[\kern-0.15em\left[{\mu_r} \right]\kern-0.15em\right]^{{-}1}\mathbf{K}_y - \left[\kern-0.15em\left[{\varepsilon_r} \right]\kern-0.15em\right] & - \mathbf{K}_y\left[\kern-0.15em\left[{\mu_r} \right]\kern-0.15em\right]^{{-}1}\mathbf{K}_x \end{bmatrix},$$

where $\mathbf {K}_x$ and $\mathbf {K}_y$ are diagonal matrices describing the wave vectors in Fourier domain; assuming only isotropic materials are discussed in the text, $\left [\kern -0.15em\left [{\varepsilon _r} \right ]\kern -0.15em\right ]$ and $\left [\kern -0.15em\left [{\mu _r} \right ]\kern -0.15em\right ]$ are Toeplitz matrices computed from the permittivity and permeability distribution of the specified layer.

Following the same procedure in [16], the solution of the field inside the layer i can then be computed to a form as the following:

(5)$$\Psi_i(\tilde{z}) = \begin{bmatrix} \vec{\mathcal{E}}_{x,y}(\tilde{z}) \\ \vec{\mathcal{H}}_{x,y}(\tilde{z}) \end{bmatrix} = \begin{bmatrix} \mathbf{W}_i & \mathbf{W}_i \\ -\mathbf{V}_i & \mathbf{V}_i \end{bmatrix} \begin{bmatrix} \mathbf{X}_i & 0 \\ 0 & \mathbf{X}_i^{{-}1} \end{bmatrix} \begin{bmatrix} \mathbf{a}_{i} \\ \mathbf{b}_{i} \end{bmatrix},$$

where $\mathbf {W}_i$, $\mathbf {V}_i$ represent the eigen-vectors of the electric and magnetic coefficients, respectively, and $\mathbf {X}_i=\exp (-\mathbf {\Gamma } k_0 h)$ is the matrix exponential of the eigen-values $\mathbf {\Gamma }^2$. Instead of directly solving this eigen-problem as in RCWA, we propose R-DIT to calculate the scattering matrix shown in the next section.

2.2 Rigorous diffraction interface theory

The general idea of R-DIT is to analytically expand the fields at the center of the optically thin layer ($z=0$) at its two boundaries ($z=\pm h/2$) by using Taylor Expansion, which builds up the boundary conditions with expanded formulas. It has been illustrated that with enough order of Taylor series, it can compute structures whose thickness is close to half wavelength [8]. In this work, we define the order of the Taylor series used in the R-DIT algorithm as R-DIT order.

We first formulate the Taylor series of the tangential electric and magnetic fields. The Fourier coefficients in (2) can be interpreted as follows:

(6)$$\vec{\mathcal{E}}_{x,y}(\tilde{z}+\tilde{\delta}) = \sum_{n=0}^\infty\frac{1}{n!}\frac{\mathrm{d}^{(n)}}{\mathrm{d}\tilde{z}^n}\vec{\mathcal{E}}_{x,y}\bigg|_{\tilde{z}}\cdot\tilde{\delta}^n = \vec{\mathcal{E}}_{x,y}(\tilde{z})+\tilde{\delta}\cdot\hat{\mathbb{P}}\cdot\vec{\mathcal{H}}_{x,y}(\tilde{z}) +\frac{1}{2}\tilde{\delta}^2\cdot\hat{\mathbb{P}}\hat{\mathbb{Q}}\cdot\vec{\mathcal{E}}_{x,y}(\tilde{z}) +\cdots,$$

(7)$$\vec{\mathcal{H}}_{x,y}(\tilde{z}+\tilde{\delta}) = \sum_{n=0}^\infty\frac{1}{n!}\frac{\mathrm{d}^{(n)}}{\mathrm{d}\tilde{z}^n}\vec{\mathcal{H}}_{x,y}\bigg|_{\tilde{z}}\cdot\tilde{\delta}^n = \vec{\mathcal{H}}_{x,y}(\tilde{z})+\tilde{\delta}\cdot\hat{\mathbb{Q}}\cdot\vec{\mathcal{E}}_{x,y}(\tilde{z}) +\frac{1}{2}\tilde{\delta}^2\cdot\hat{\mathbb{Q}}\hat{\mathbb{P}}\cdot\vec{\mathcal{H}}_{x,y}(\tilde{z}) +\cdots,$$

where $\tilde {\delta }=k_0\delta$ is the normalized length on $z$ direction. Substituting $\tilde {z}+\tilde {\delta }=0$ into the left-hand side and $\tilde {\delta }=\pm \tilde {h}/2$ to the right-hand side of (6), (7) separately, the Fourier series of the fields at the center expanded at the two boundaries of the layer can be formulated to the following two boundary conditions:

(8)$$\begin{bmatrix} \vec{\mathcal{E}}_{x,y} \\ \vec{\mathcal{H}}_{x,y} \end{bmatrix}\Bigg|_{\tilde{z}=0} = \hat{\mathbb{T}}^{(1)} \begin{bmatrix} \vec{\mathcal{E}}_{x,y} \\ \vec{\mathcal{H}}_{x,y} \end{bmatrix}\Bigg|_{\tilde{z}={-}\frac{\tilde{h}}{2}} , \hat{\mathbb{T}}^{(1)} = \begin{bmatrix} \hat{\mathbb{I}}+\frac{\tilde{h}^2}{8}\hat{\mathbb{P}}\hat{\mathbb{Q}}+\cdots & \frac{\tilde{h}}{2}\hat{\mathbb{P}}+\frac{\tilde{h}^3}{48}\hat{\mathbb{P}}\hat{\mathbb{Q}}\hat{\mathbb{P}}+\cdots \\ \frac{\tilde{h}}{2}\hat{\mathbb{Q}}+\frac{\tilde{h}^3}{48}\hat{\mathbb{Q}}\hat{\mathbb{P}}\hat{\mathbb{Q}}+\cdots & \hat{\mathbb{I}}+\frac{\tilde{h}^2}{8}\hat{\mathbb{Q}}\hat{\mathbb{P}}+\cdots \end{bmatrix},$$

(9)$$\begin{bmatrix} \vec{\mathcal{E}}_{x,y} \\ \vec{\mathcal{H}}_{x,y} \end{bmatrix}\Bigg|_{\tilde{z}=0} = \hat{\mathbb{T}}^{(2)} \begin{bmatrix} \vec{\mathcal{E}}_{x,y} \\ \vec{\mathcal{H}}_{x,y} \end{bmatrix}\Bigg|_{\tilde{z}={+}\frac{\tilde{h}}{2}} , \hat{\mathbb{T}}^{(2)} = \begin{bmatrix} \hat{\mathbb{I}}+\frac{\tilde{h}^2}{8}\hat{\mathbb{P}}\hat{\mathbb{Q}}+\cdots & -\frac{\tilde{h}}{2}\hat{\mathbb{P}}-\frac{\tilde{h}^3}{48}\hat{\mathbb{P}}\hat{\mathbb{Q}}\hat{\mathbb{P}}+\cdots \\ -\frac{\tilde{h}}{2}\hat{\mathbb{Q}}-\frac{\tilde{h}^3}{48}\hat{\mathbb{Q}}\hat{\mathbb{P}}\hat{\mathbb{Q}}+\cdots & \hat{\mathbb{I}}+\frac{\tilde{h}^2}{8}\hat{\mathbb{Q}}\hat{\mathbb{P}}+\cdots \end{bmatrix},$$

for simplicity, here we reduce $\hat {\mathbb {T}}^{(1)}$ and $\hat {\mathbb {T}}^{(2)}$ to:

(10)$$\hat{\mathbb{T}}^{(1)} = \begin{bmatrix} \underline{\mathcal{A}} & \underline{\mathcal{B}} \\ \underline{\mathcal{C}} & \underline{\mathcal{D}} \end{bmatrix},\quad \hat{\mathbb{T}}^{(2)} = \begin{bmatrix} \underline{\mathcal{A}} & -\underline{\mathcal{B}} \\ -\underline{\mathcal{C}} & \underline{\mathcal{D}} \end{bmatrix}\,.$$

The fields at the two boundaries given in (8) and (9) can be written as the fields of the surrounding layers $\Psi _{i-1}$ and $\Psi _{i+1}$ (shown in Fig. 1). In this work, we construct each layer to be solved as a layer surrounded by two gap layers (layers $i-1$ and $i+1$) with zero thickness, and then relate all solved layers by the Redheffer star product [16–18]. Under these circumstances, (8) and (9) can be further rewritten as:

(11)$$\begin{bmatrix} \vec{\mathcal{E}}_{x,y} \\ \vec{\mathcal{H}}_{x,y} \end{bmatrix}\Bigg|_{\tilde{z}=0} = \begin{bmatrix} \underline{\mathcal{A}} & \underline{\mathcal{B}} \\ \underline{\mathcal{C}} & \underline{\mathcal{D}} \end{bmatrix} \begin{bmatrix} \mathbf{W_0} & \mathbf{W_0} \\ -\mathbf{V_0} & \mathbf{V_0} \end{bmatrix} \begin{bmatrix} \mathbf{a_{i-1}}\\ \mathbf{b_{i-1}} \end{bmatrix},$$

(12)$$\begin{bmatrix} \vec{\mathcal{E}}_{x,y} \\ \vec{\mathcal{H}}_{x,y} \end{bmatrix}\Bigg|_{\tilde{z}=0} = \begin{bmatrix} \underline{\mathcal{A}} & -\underline{\mathcal{B}} \\ -\underline{\mathcal{C}} & \underline{\mathcal{D}} \end{bmatrix} \begin{bmatrix} \mathbf{W_0} & \mathbf{W_0} \\ -\mathbf{V_0} & \mathbf{V_0} \end{bmatrix} \begin{bmatrix} \mathbf{b_{i+1}}\\ \mathbf{a_{i+1}} \end{bmatrix},$$

where $\mathbf {W_0}$ and $\mathbf {V_0}$ are eigen-vectors of the gap layers which can be easily derived analytically [16].

Combining (11) and (12), we can obtain the scattering matrix of layer $i$ after a few derivation steps as follows:

(13)$$\begin{aligned} \mathbf{S}^i_{11} & =\mathbf{S}^i_{22}=\frac{1}{2}[\mathbf{M}_1-\mathbf{M}_2],\\ \mathbf{S}^i_{12} & =\mathbf{S}^i_{21}=\frac{1}{2}[\mathbf{M}_1+\mathbf{M}_2],\\ \mathbf{M}_1 & =(-\underline{\mathcal{C}}\mathbf{W_0}-\underline{\mathcal{D}}\mathbf{V_0})^{{-}1}(\underline{\mathcal{C}}\mathbf{W_0}-\underline{\mathcal{D}}\mathbf{V_0}),\\ \mathbf{M}_2 & =(\underline{\mathcal{A}}\mathbf{W_0}+\underline{\mathcal{B}}\mathbf{V_0})^{{-}1}(\underline{\mathcal{A}}\mathbf{W_0}-\underline{\mathcal{B}}\mathbf{V_0})\,.\\ \end{aligned}$$

To conclude, our proposed derivation of the R-DIT algorithm greatly simplifies the computations of the scattering matrices, requiring only basic matrix operations such as multiplication and inversion. Since the Taylor approximation is used, the approximation error can be estimated from the R-DIT order and the thickness of the interested layer (See Supplement 1 Section 1 for more analysis on the convergence and error estimations of our proposed R-DIT algorithm).

3. Inverse design workflows and numerical results

In this section, we begin by elucidating the inverse design workflow incorporating our proposed TorchRDIT. To exemplify this, we present case studies including the inverse design of meta-atoms with parameter constraints and the creation of meta-atoms with a free-form structure.

Fig. 2 shows the schematic of the workflow, where TorchRDIT functions as a differentiable electromagnetic solver here, taking patterns of meta-atoms as pixelated matrices. Here, we utilize topology optimization techniques [19–21] to construct the permittivity distributions of our devices from a density matrix $\rho _i\in [0,1]$ as follows:

(14)$$\varepsilon_{r,i} = \varepsilon_{r,\mathrm{min}}^i + \rho_i\cdot(\varepsilon_{r,\mathrm{max}}^i - \varepsilon_{r,\mathrm{min}}^i),$$

where $\varepsilon _{r,\mathrm {min}}^i$ and $\varepsilon _{r,\mathrm {max}}^i$ are user defined values according to the materials of layer $i$. When implementing this workflow with PyTorch, $\rho _i$ is initialized as a tensor variable in PyTorch, and all the operations on it are tracked by torch.autograd, which is the automatic differentiation engine of PyTorch [10].

Fig. 2. Inverse design workflows. The proposed workflow constructs computational graphs that transition from initial input vectors to a user-defined loss function utilizing the PyTorch framework, enabling a seamless end-to-end optimization process. In the process of generating optical geometries for resolution, the approach is divided into two kinds: (a) design parameters are initially synthesized by a neural network and then form the geometries by topology operations; (b) for the representation of free-form geometries, a neural network with sinusoidal activation, known as SIREN[11], is employed. These generated geometries, alongside their associated computational graphs, are subsequently fed into TorchRDIT and optimized through iterative minimization of the loss calculated from the electromagnetic responses obtained. This iterative process is designed to converge on the optimally designed geometries that align with the predefined objectives.

Download Full Size | PDF

3.1 Inverse design of structures constrained by parameters

To illustrate the efficacy of TorchRDIT in the inverse design of geometries bound by parameters, we construct a model of Huygens meta-atoms operating at a wavelength, $\lambda _0=5.2\,{\mu}\textrm{m}$, with a periodicity of $p=2.5\,{\mu}\textrm{m}$, on the PbTe-CaF$_2$ structures as characterized in [22]. Within this framework, we specify four distinct shape categories as depicted in Fig. 3(a), all maintaining a consistent PbTe layer height of $650\,\mathrm {nm}$. The primary objective is to inverse design meta-devices that achieve our predetermined phase targets and boast transmission efficiencies surpassing the 0.6 threshold. Following the process flow illustrated in Fig. 2(a), we consequently establish our loss function as follows:

(15)$$\mathcal{L}={-}\frac{\boldsymbol{t}_{xn}}{|\boldsymbol{t}_{xn}|}\exp{({-}j\phi_{\mathrm{target}})}^*\cdot w_1+\mathcal{L}_{\mathrm{MSE}}(\mathbf{Trn}_x, 1)\cdot w_2,$$

where the first term is tasked with calculating the discrepancy between the phase of the transmitted x-polarized electric fields $\boldsymbol {t}_{xn}$, and the user-stipulated target phase $\phi _{\mathrm {target}}$, expressed in radians; the second term assesses the deviation (where $\mathcal {L}_{\mathrm {MSE}}$ symbolizes the mean squared error) between the transmission efficiency for x-polarization $\mathbf {Trn}_x$ and the desired efficiency. These two loss functions are moderated by two individual weights, $w_1, w_2$, offering the capacity to adjust them as part of the hyperparameters within the operational workflow. For this demonstration, we maintain the truncated Fourier modes on a $21\times 21$ scale and set the R-DIT order to 15.

Fig. 3. Schematics and optimized results of parameter constrained devices. To demonstrate the case of parameter-constrained devices, we define four types of geometries illustrated in (a). (b) shows a 3D structure of an H-shape PbTe meta-atom array on the PbTe-CaF$_2$ substrate. The PbTe layer are model with a two-layer model ($n_{\mathrm {top}}=4.8; n_{\mathrm {bottom}}=5.4, k_{\mathrm {bottom}}=0.01$) as in [22], and $n_{\mathrm {CaF_2}=1.4}$ for the substrate. The simulated transmission efficiencies and phases using TorchRDIT (marked as solid lines with filled plus markers) and CST Microwave Studio (marked as hollow squares) are shown in (d). (c) illustrates the number of iterations cost when each device is optimized. Insets of (d): Grid view of the inverse designed shapes of PbTe meta-atoms (See Supplement 1 Section 2 for the dimension of each meta-atom).

Download Full Size | PDF

Fig. 3(d) illustrates all inversely optimized meta-atoms spanning the phase from $0$ to $2\pi$ and achieving transmission efficiencies in excess of 0.6. A significant proportion of the devices attains efficiencies that exceed 0.75. By amalgamating diverse geometric forms, our proposed framework also successfully addresses the issue of low-efficiency gaps, as previously highlighted in the study by [22]. Crucially, the workflow we propose can autonomously generate these outcomes within a span of one thousand iterations for each meta-atom, as evidenced in Fig. 3(c). This convergence rate is notably expeditious when juxtaposed with the conventional manual/empirical design techniques. Furthermore, the neural network at the core of this workflow operates effectively without requiring a pre-training phase using any dataset.

3.2 Inverse design of reconfigurable free-form mata-atoms

This section focuses on the application of free-form meta-atoms to realize sophisticated functionalities, specifically utilizing optical phase change materials (O-PCMs). The choice of O-PCMs is strategic, as they amplify the capability of dynamically modulating light while preserving their slim profile. We aim to conceptualize a reconfigurable transmitting beam-deflector operating at a central wavelength of $\lambda _0=4\,{\mu}\textrm{m}$ based on Ge$_2$Sb$_2$Se$_4$Te$_1$ (GSST) for its ability to facilitate the reversible transition between amorphous and crystalline states through electrical control. This transition yields a significant refractive index contrast, thus enabling versatile optical devices with dynamic responses [23,24].

3.2.1 Definition of the loss function

The challenge of designing reonfigurable/tunable optical devices with high efficacy has traditionally been formidable using conventional optimization methods. However, this challenge is adeptly overcome by our proposed inverse design workflows, as illustrated in Fig. 2. We postulate the functionality of this particular meta-atom, with a focus on the transverse magnetic (TM) mode (x-polarized) diffraction along the x-axis, to encompass two states: (1) in the amorphous state of GSST, it should direct the normal incident light towards the +1 order diffraction, while concurrently minimizing the energy directed into the -1 and 0 diffraction orders; (2) in the crystalline state, it aims to minimize the transmitted energies at the operating wavelength across the -1, 0, and +1 diffraction orders.

Given that our consideration is limited to x-directional diffraction, we can simplify our geometries by enforcing a symmetry constraint about the x-axis and reducing the $L_y$ dimension to eliminate y-directional diffraction. Using the operational parameters defined above, we subsequently articulate our loss function accordingly:

(16)$$\begin{aligned} \mathcal{L}_{\mathrm{sum}}= & \mathcal{L}_{\mathrm{MSE}}(\boldsymbol{T}_{A,({+}1,0)}, 1)\cdot w_1+ \mathcal{L}_{\mathrm{MSE}}(\boldsymbol{T}_{A,(0,0)}, 0)\cdot w_2+\\ & (\boldsymbol{T}_{A,(0,0)} - \boldsymbol{T}_{A,({+}1,0)})\cdot w_3+ (\boldsymbol{T}_{A,({-}1,0)} - \boldsymbol{T}_{A,({+}1,0)})\cdot w_4+\\ & \mathcal{L}_{\mathrm{MSE}}(\boldsymbol{T}_{C,({+}1,0)}, 0)\cdot w_5+ \mathcal{L}_{\mathrm{MSE}}(\boldsymbol{T}_{C,(0,0)}, 0)\cdot w_6,\\ \end{aligned}$$

where $\boldsymbol {T}_{A, (x, 0)}$ and $\boldsymbol {T}_{C, (x, 0)}$ correspond to the transmitted power levels in the amorphous and crystalline states, respectively, for the $x$th order of diffraction; The term $w_i$ ($i=1, 2,\ldots, 6$) denotes the weights assigned to each loss term. Notably, the third and fourth terms in (16) have the potential to yield negative values. Consequently, the scalar loss values derived from the function are not bounded as positive and vary according to the assigned weights.

Although possibilities exist to define an extensive array of loss terms to closely align the inversely designed behaviors with the predefined conditions, imposing such strict constraints can instead complicate the optimization process. Specifically, stringent bounds may hinder the convergence of the optimization algorithm and necessitate an increase in hyperparameters, which would then require more extensive manual tuning workloads.

3.2.2 Implicit neural representations of free-form geometries

In this work, we use neural networks to implicitly represent the geometries of the free-form meta-atoms, as shown in Fig. 2(b). Unlike classical image processing models that treat geometries as discrete pixels with certain resolutions [25], implicit neural networks map arbitrary spatial coordinates to a level set function, which is then binarized into fabricable permittivity distribution. Here, we use sinusoidal representation networks (SIREN) [11] to perform this mapping. We establish an implicit representation that intakes $64\times 64$ uniformly sampled spatial coordinates ($\mathbf {x}=(x, y)$) and outputs the corresponding density distributions ($\rho _i=\Phi _i(\mathbf {x}_i)$) as described in (14). In this work, we use a SIREN with 16 hidden layers with 6 hidden features of each. The frequency factors are set to 12 for the first and hidden layers. Specifically, TorchRDIT manages the mapping between the implicit representation $\Phi _i$ and the loss function $\mathcal {L}_\mathrm {sum}(\Phi _i)$, thereby facilitating the optimization of the topological forms to meet the defined requirements.

Beyond the matrices defining the geometries, TorchRDIT can also treat additional simulation parameters as differentiable entities, which extends the flexibility of the design process to encompass a broader set of physical attributes, such as the layer thicknesses illustrated in Fig. 2(b). Within this context, the parameters $h_1$, $h_2$, and $h_3$ are identified as variables subject to optimization. Following a similar topology optimization method, we define the thickness as follows:

(17)$$h_{i} = h_{\mathrm{min}}^i + \theta_i\cdot(h_{\mathrm{max}}^i - h_{\mathrm{min}}^i),$$

where $h_{\mathrm {max}}^i$ and $h_{\mathrm {min}}^i$ are defined by users, and $\theta _i\in [0,1]$ are vectors tracked and updated iteratively by PyTorch.

3.2.3 Results

Adhering to the inverse design workflow delineated in Fig. 2(b), which utilizes identical spatial modes and orders discussed in the preceding section, we successfully optimized a meta-atom aligning with our defined goals. A graphical representation of this optimized meta-atom is sketched in Fig. 4(c) and (d). Moreover, the simulated spectral response of this meta-atom covering from $3\,{\mu}\textrm{m}$ to $5\,{\mu}\textrm{m}$ is shown in Fig. 5. It provides insight into the optical performance of the meta-atom across the specified range.

Fig. 4. Structures and materials of the GSST meta-atoms. (a) illustrates the device configuration we use for the beam deflector. The protective SiN$_x$ layer covers the GSST meta-atom resting on the doped-silicon heater. The permittivity profiles of GSST and the doped-silicon are shown in (b). The permittivity profiles of Si, SiO2, and SiN$_x$ are offered by [26], [27,28] and [29] respectively. (c) and (d) present the top view and 3D perspective of the optimized meta-atom. (Thickness: $h_0=1.70\,\mu m$, $h_1=87.7\,nm$, $h_2=3.20\,\mu m$ and $h_3=0.33\,\mu m$; Dimensions: $L_x=4.25\,\mu m$ and $L_y=1.8\,\mu m$).

Download Full Size | PDF

Fig. 5. Simulated results of the optimized meta-atom in Fig. 4. Utilizing the TorchRDIT, we conduct simulations on the optimized meta-atom and validate the simulations with the assistance of S$^4$, an open-source software package for electromagnetic analysis [30], with a maximum spatial mode of 400. (a), (b) show the simulated TM mode (x-polarized) spectrum under amorphous and crystalline states, and (c), (d) present the TE mode (y-polarized) results.

Download Full Size | PDF

Fig. 5(a) illustrates the optimized response in TM mode: (1) the engineered meta-atom directs approximately 60${\% }$ of energy into the +1 diffraction order while effectively suppressing other mode transmissions; (2) the bandwidth corresponding to the +1 order diffraction extends to approximately $1\,{\mu}\textrm{m}$. On the flip side, although we have not imposed specific constraints on the transverse electric (TE) mode behaviors, its TE spectrum also indicates a higher transmission at the +1 order diffraction in the amorphous state (see Fig. 5(c), (d)).

The discrepancies observed in the simulation outcomes between TorchRDIT and S$^4$, as illustrated in Fig. 5, are attributable to boundary extraction and edge smoothing procedures, which are necessary when translating the pixelated patterns TorchRDIT generated into polygonal shapes compatible with S$^4$. The current geometrical matrices utilized in TorchRDIT have a resolution of $64\times 64$. An enhancement in the resolution within TorchRDIT would be conducive to minimizing the mismatch between the simulated results of the two systems, thereby yielding a more accurate representation of the meta-atom’s performance as predicted by the simulation.

4. Benchmark

In this section, we evaluate the single-precision performance of TorchRDIT on the CPU (AMD Ryzen Threadripper 1950X) and GPU (NVIDIA A100) platforms. For comparative analysis, we have integrated the differentiable RCWA method [5] into TorchRDIT to ensure that the respective algorithms compute only layers with inhomogeneous geometries. Furthermore, the tested photonic structure consists of only one inhomogeneous layer, minimizing the impact of other simulation stages on the performance evaluation.

Our analysis of TorchRDIT begins with its performance on CPU platforms. The RCWA algorithm, requiring eigen-decomposition for forward computation, contrasts with R-DIT’s reliance on matrix multiplications of $\hat {\mathbb {P}}$ and $\hat {\mathbb {Q}}$. Experimentally, R-DIT’s simulation time remains relatively flat as its order increases (see Fig. 6(b)) and achieves a 1.20 times speed-up over RCWA with high spatial modes and high R-DIT orders (Fig. 6(d)). When it comes to the backward computation, the R-DIT algorithm performs twice as fast with high spatial modes, detailed in (Fig. 6(e)) because gradient calculations required by RCWA are impeded by the matrix solution computations, which become particularly burdensome for large dimensional matrices [5]. Considering the complete optimization process consisting of both forward and backward computations, the R-DIT algorithm demonstrates a speed-up ranging from 1.2 to 2.2 times when compared to the differentiable RCWA algorithm, as evidenced in Fig. 6(f).

Fig. 6. Benchmark of R-DIT solver on the CPU platform. (a) Simulation time (per iteration) of the forward and backward computations of the differentiable RCWA. (b)-(c): The forward computation time (b) and the backward computation time (c) of R-DIT per iteration. (d)-(f) Speed-up of the forward computations (d), backward computations (e), and the total computation time (forward+backward) (f) of R-DIT over RCWA.

Download Full Size | PDF

Our evaluation of TorchRDIT’s performance on GPU platforms identifies a critical limitation with RCWA: eigen-decomposition on PyTorch is CPU-bound, leading to increased computation time and data transfer delays between the CPU and GPU. As Fig. 6(a) and 7(a) demonstrate, for high spatial mode simulations, the GPU-executed RCWA can be even slower than its CPU counterpart. On the other hand, GPUs excel in parallel processing tasks. As the spatial modes increase, the GPU performance of RCWA is highly restricted due to the large dimensions of the matrices, resulting in a dramatic peaking of the speed-up of R-DIT as shown in Fig. 7(e). When considering complete iterations, R-DIT achieves up to a 16.2 speed-up over the RCWA (see Fig. 7(f)), underscoring R-DIT’s superior efficiency and scalability on GPU platforms for high spatial mode scenarios.

Fig. 7. Benchmark of R-DIT solver on the GPU platform. (a) Simulation time (per iteration) of the forward and backward computations of the differentiable RCWA. (b)-(c): The forward computation time (b) and the backward computation time (c) of R-DIT per iteration. (d)-(f) Speed-up of the forward computations (d), backward computations (e), and the total computation time (forward+backward) (f) of R-DIT over RCWA.

Download Full Size | PDF

5. Conclusion

In conclusion, our research presented an enhanced R-DIT algorithm that serves as the core of the proposed differentiable solver TorchRDIT, significantly boosting the inverse design speed. Through TorchRDIT, we illustrated the inverse design of parameter-determined meta-atoms to yield high-efficiency Huygens meta-atoms with precise phase control. Furthermore, we have utilized SIREN to facilitate the design of reconfigurable free-form GSST-based meta-atoms. Validations were achieved using CST Microwave Studio and S$^4$.

The benchmark results confirmed that the differentiable R-DIT algorithm surpasses the differentiable RCWA speed on both CPU and GPU platforms. Our method effectively bypasses the complexities of eigendecomposition, enhancing the numerical stability and computational efficiency, particularly for GPU architectures. Combined with emerging techniques of the machine learning community and rapidly evolving computing hardware, this work further lays the foundation for realizing fast inverse design of large-scale, three-dimensional, high-performance photonic structures with complex functionalities. Meanwhile, both TorchRDIT and the proposed inverse design workflow contribute to the machine learning community as one form of a physics-informed neural network (PINN) [31,32].

Funding

National Science Foundation (2118787, 2132929).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper may be obtained from the authors upon reasonable request.

Supplemental document

See Supplement 1 for supporting content.

References

1. M. M. R. Elsawy, S. Lanteri, R. Duvigneau, et al., “Numerical Optimization Methods for Metasurfaces,” Laser Photonics Rev. 14(10), 1900445 (2020). [CrossRef]

2. Z. Li, R. Pestourie, J.-S. Park, et al., “Inverse design enables large-scale high-performance meta-optics reshaping virtual reality,” Nat. Commun. 13(1), 2409 (2022). [CrossRef]

3. M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc. Am. 71(7), 811–818 (1981). [CrossRef]

4. Z. Zhu and C. Zheng, “Differentiable scattering matrix for optimization of photonic structures,” Opt. Express 28(25), 37773–37787 (2020). [CrossRef]

5. S. Colburn and A. Majumdar, “Inverse design and flexible parameterization of meta-optics using algorithmic differentiation,” Commun. Phys. 4(1), 65 (2021). [CrossRef]

6. Yee Kane, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14(3), 302–307 (1966). [CrossRef]

7. Z. Bai, J. Demmel, J. Dongarra, et al., Templates for the Solution of Algebraic Eigenvalue Problems (Society for Industrial and Applied Mathematics, 2000), chap. 7.

8. C. M. Roberts and V. A. Podolskiy, “Rigorous diffraction interface theory,” Appl. Phys. Lett. 110(17), 171108 (2017). [CrossRef]

9. Y. Huang, H. Tang, B. Zheng, et al., “Inverse Design of Photonic Structures Using Automatic Differentiable Rigorous Diffraction Interface Theory,” in CLEO 2023, (Optica Publishing Group, 2023), p. JTu2A.119.

10. A. Paszke, S. Gross, S. Chintala, et al., “Automatic differentiation in PyTorch,” in NIPS 2017 Workshop Autodiff, (2017).

11. V. Sitzmann, J. N. Martel, A. W. Bergman, et al., “Implicit neural representations with periodic activation functions,” in Proc. NeurIPS, (2020).

12. M. Y. Shalaginov, S. D. Campbell, S. An, et al., “Design for quality: Reconfigurable flat optics based on active metasurfaces,” Nanophotonics 9(11), 3505–3534 (2020). [CrossRef]

13. K. Du, H. Barkaoui, X. Zhang, et al., “Optical metasurfaces towards multifunctionality and tunability,” Nanophotonics 11(9), 1761–1781 (2022). [CrossRef]

14. A. M. Shaltout, V. M. Shalaev, and M. L. Brongersma, “Spatiotemporal light control with active metasurfaces,” Science 364(6441), eaat3100 (2019). [CrossRef]

15. Y. Huang and Z. Zhu, “TorchRDIT,” https://github.com/yi-huang-1/torchrdit (2023).

16. R. C. Rumpf, “Improved Formulation of Scattering Matrices for Semi-Analytical Methods That Is Consistent with Convention,” Prog. Electromagn. Res. B 35, 241–261 (2011). [CrossRef]

17. R. Redheffer, “Difference equations and functional equations in transmission-line theory,” Modern mathematics for the engineer 12, 282–337 (1961).

18. M. G. Moharam and A. B. Greenwell, “Efficient rigorous calculations of power flow in grating coupled surface-emitting devices,” in Photon Management, vol. 5456 (SPIE, 2004), pp. 57–67.

19. Y. Tsuji, K. Hirayama, T. Nomura, et al., “Design of optical circuit devices based on topology optimization,” IEEE Photonics Technol. Lett. 18(7), 850–852 (2006). [CrossRef]

20. J. Jensen and O. Sigmund, “Topology optimization for nano-photonics,” Laser Photonics Rev. 5(2), 308–321 (2011). [CrossRef]

21. D. Sell, J. Yang, S. Doshay, et al., “Large-Angle, Multifunctional Metagratings Based on Freeform Multimode Geometries,” Nano Lett. 17(6), 3752–3757 (2017). [CrossRef]

22. L. Zhang, J. Ding, H. Zheng, et al., “Ultra-thin high-efficiency mid-infrared transmissive Huygens meta-optics,” Nat. Commun. 9(1), 1481 (2018). [CrossRef]

23. Y. Zhang, C. Fowler, J. Liang, et al., “Electrically reconfigurable non-volatile metasurface using low-loss optical phase-change material,” Nat. Nanotechnol. 16(6), 661–666 (2021). [CrossRef]

24. M. Y. Shalaginov, S. An, Y. Zhang, et al., “Reconfigurable all-dielectric metalens with diffraction-limited performance,” Nat. Commun. 12(1), 1225 (2021). [CrossRef]

25. Y. Dong, S. An, B. Zheng, et al., “Data Collection and Network Design for Deep Learning Based Metasurface Design,” in 2023 International Applied Computational Electromagnetics Society Symposium (ACES), (2023), pp. 1–2.

26. E. Shkondin, O. Takayama, M. E. A. Panah, et al., “Large-scale high aspect ratio Al-doped ZnO nanopillars arrays as anisotropic metamaterials,” Opt. Mater. Express 7(5), 1606–1627 (2017). [CrossRef]

27. I. H. Malitson, “Interspecimen Comparison of the Refractive Index of Fused Silica*,†,” J. Opt. Soc. Am. 55(10), 1205–1209 (1965). [CrossRef]

28. C. Z. Tan, “Determination of refractive index of silica glass for infrared wavelengths by IR spectroscopy,” J. Non-Cryst. Solids 223(1-2), 158–163 (1998). [CrossRef]

29. K. Luke, Y. Okawachi, M. R. E. Lamont, et al., “Broadband mid-infrared frequency comb generation in a Si₃N₄ microresonator,” Opt. Lett. 40(21), 4823–4826 (2015). [CrossRef]

30. V. Liu and S. Fan, “S⁴ : A free electromagnetic solver for layered periodic structures,” Comput. Phys. Commun. 183(10), 2233–2244 (2012). [CrossRef]

31. S. Cuomo, V. S. Di Cola, F. Giampaolo, et al., “Scientific Machine Learning Through Physics–Informed Neural Networks: Where we are and What’s Next,” J. Sci. Comput. 92(3), 88 (2022). [CrossRef]

32. Z. K. Lawal, H. Yassin, D. T. C. Lai, et al., “Physics-Informed Neural Network (PINN) Evolution and Beyond: A Systematic Literature Review and Bibliometric Analysis,” Big Data Cogn. Comput. 6(4), 140 (2022). [CrossRef]

Eigendecomposition-free inverse design of meta-optics devices

Abstract

1. Introduction

2. Theoretical formulation of rigorous diffraction interface theory

2.1 Definition of scattering matrix

2.2 Rigorous diffraction interface theory

3. Inverse design workflows and numerical results

3.1 Inverse design of structures constrained by parameters

3.2 Inverse design of reconfigurable free-form mata-atoms

3.2.1 Definition of the loss function

3.2.2 Implicit neural representations of free-form geometries

3.2.3 Results

4. Benchmark

5. Conclusion

Funding

Disclosures

Data availability

Supplemental document

References

Supplementary Material (1)

Data availability

Cited By

Figures (7)

Equations (17)

Optics Express