Topology optimization of electromagnetic metamaterials

Filtering

• convolution: \(f(\B r) = \frac{1}{A}{\rm exp}(-|\B r|^2 /R_f^2)\), with \(\int_{\Omega_{\rm des}} f(\B r) =1\) \[\begin{equation*} \densf(\B r) = p * f = \int_{\Omega_{\rm des}} p(\B r') f(\B r -\B r') {\rm d} \B r' \label{eq:gaussian_filt} \end{equation*}\]
• PDE: \(-R_f^2 \B\nabla ^2 \densf + \densf = p {\quad\rm on \,}\Omega_{\rm des}, \grad\densf\cdotp\B n = 0 {\quad\rm on \,}\partial\Omega_{\rm des}\)⁴

Filtering

Projection

\[\densp(\densf) = \frac{\tanh\left[\beta\nu\right] + \tanh\left[\beta(\densf-\nu)\right] }{\tanh\left[\beta\nu\right] + \tanh\left[\beta(1-\nu)\right]}\] with \(\nu=1/2\) and \(\beta>0\) increased during the course of the optimization.⁵

Projection

Interpolation

\(\varepsilon(\densp)=(\varepsilon_{\rm max}-\varepsilon_{\rm min})\,\densp^m + \varepsilon_{\rm min}\)⁶

Finite differences

\[\begin{equation} \frac{\mathrm{d}\B G}{\mathrm{d}p_i} \approx \frac{\B G(\B p + h \B e_i) - \B G(\B p)}{h} \end{equation}\] where \(\B e_i\) is the vector with \(0\) in all entries except for \(1\) in the \(i^{th}\) entry.

• numerical inaccuracy
• expensive for large \(M\) and/or \(N\)

Tangent linear equation

Explicitly, the gradient can be computed applying the chain rule: \[\begin{equation} \frac{\mathrm{d}\B G}{\mathrm{d}\B p} = \frac{\partial \B G}{\partial \B p} + \frac{\partial \B G}{\partial \B u} \frac{\mathrm{d}\B u}{\mathrm{d}\B p}. \label{eq:gradG} \end{equation}\] Taking the total derivative of Eq. (\(\ref{eq:genimpeq}\)) we obtain the tangent linear equation: \[\begin{equation} {\frac{\partial \B F(\B u, \B p)}{\partial \B u}} {\frac{\mathrm{d}\B u}{\mathrm{d}\B p}} = {-\frac{\partial \B F(\B u, \B p)}{\partial \B p}}. \end{equation}\]

Adjoint equation

Assuming the tangent linear system is invertible, we can rewrite the Jacobian as: \[\begin{equation} \frac{\mathrm{d}\B u}{\mathrm{d}\B p} = - \left(\frac{\partial \B F(\B u, \B p)}{\partial \B u}\right)^{-1} \frac{\partial \B F(\B u, \B p)}{\partial \B p}. \end{equation}\] After substituting this value in (\(\ref{eq:gradG}\)) and taking the adjoint (Hermitian transpose, denoted by \(\dagger\)) we get: \[\begin{equation} \frac{\mathrm{d}\B G}{\mathrm{d}\B p}^{\dagger} = \frac{\partial \B G}{\partial \B p}^{\dagger} - \frac{\partial \B F(\B u, \B p)}{\partial \B p}^{\dagger} \left(\frac{\partial \B F(\B u, \B p)}{\partial \B u}\right)^{-\dagger} \frac{\partial \B G}{\partial \B u}^{\dagger} . \end{equation}\]

Adjoint equation

Defining the adjoint variable \(\B \lambda\) as: \[\begin{equation} \B \lambda = \left(\frac{\partial \B F(\B u, \B p)}{\partial \B u}\right)^{-\dagger} \frac{\partial \B G}{\partial \B u}^{\dagger} \end{equation}\] we obtain the adjoint equation \[\begin{equation} \left(\frac{\partial \B F(\B u, \B p)}{\partial \B u}\right)^{\dagger} \B \lambda = \frac{\partial \B G}{\partial \B u}^{\dagger}. \end{equation}\]

Automatic differentiation (AD)

\(f: \mathbb{R}^n \rightarrow \mathbb{R}^m\)

Example: \(f(x_1,x_2) = x_1x_2 + \sin(x_1)\)

Finite Element Method (FEM)

• The two polarizations decouple and Maxwell’s equations can be combined into the following scalar propagation equation: \[\begin{equation} \mathcal {M}_{\tens{\xi},\chi}(u) = \ddiv \left[\tens{\xi}\grad u\right] + k_0^2 \chi u= 0,\label{eq:propa} \end{equation}\]
  • TE polarization: \(u=H_z\), \(\tens{\xi} = \tens{\e_\parallel}^T/\det\tens{\e_\parallel}\), \(\chi = \muzz\)
  • TM polarization: \(u=E_z\), \(\tens{\xi} = \tens{\mu_\parallel}^T/\det\tens{\mu_\parallel}\), \(\chi = \ezz\)

Finite Element Method

• Excitation \(u_0\), the diffracted field \(u_s = u - u_0\) must satisfy an Outgoing Waves Condition.
  • Plane wave: \(u_0=\exp(i\B k\cdot\B r)\)
  • Line source: Green’s function \(u_0=G(\B r,\B r',k)=-\frac{i}{4}H_0^{(1)}(k|\B r- \B r'|)\)
• Scattering problem (\(\ref{eq:propa}\)) is recast into a radiation problem with a source \(\mathcal {S}\) localized inside the object: \[\begin{equation} \mathcal {M}_{\tens{\xi},\chi}(u_s) = -\mathcal {M}_{\tens{\xi}-\tens{\xi_b},\chi-\chi_b}(u_0) = \mathcal {S}.\label{eq:rad} \end{equation}\]
• Weak form: \[\begin{equation} \mathcal {R}_{\tens{\xi},\chi}(u,v) = -\int_\Omega \,\tens{\xi}\grad u \grad v^{\star} {\dd\B r}+ \int_{\Omega} \, \chi u v^\star{\dd\B r} + \int_{\D\Omega}\, (\tens{\xi}\grad u) v^\star {\dd s}. \label{eq:residual} \end{equation}\] Surface term = 0 since \((\tens{\xi}\grad u)\cdotp \B n = 0\) on the outer boundaries of the PMLs.
• Find \(u\in\LSO\) such that: \[\begin{equation} \mathcal {R}_{\tens{\xi},\chi}(u,v) = 0, \quad \forall v \in \LSO \label{eq:weak} \end{equation}\]

Application: cloaking

¹⁰

Cloaking

Objective: minimize scattered field in the background \[\begin{equation*} \min_{p(\B r)} \quad \phi=\gamma\phi_1 + (1-\gamma) \phi_2 \end{equation*}\] \[\begin{equation*} \phi_1=\frac{1}{W}\int_{\Omega_{\rm B}}\left|E_{\rm d}^{\rm cloak}\right|^2 {\rm d} \B r\end{equation*}\] \[\begin{equation*} \phi_2=\int_{\Omega_{\rm D}}\frac{h_{\rm m}^2}{S}\left| \B \nabla\,\varepsilon \right|^2 {\rm d} \B r\end{equation*}\]

• \(W=\int_{\Omega_{\rm B}}\left|E_{\rm d}^{\rm object}\right|^2{\rm d} \B r\) (uncloaked case)
• \(\phi_2\): alternative approach for filtering, \(S=(R_{\rm c}-R_{\rm s})^2/2\) is the area of the design domain and \(h_{\rm m}\) is the maximum size of a mesh element
• weight parameter \(\gamma\) between 0 and 1

Cloaking

Sensitivity to source

Bandwidth enhancement with multi-frequency optimization

Modal analysis and the origin of cloaking effect

Illusion

Objective: make a metallic cylinder look like an arbitrarily shaped dielectric \[\begin{equation*} \min_{p(\B r)} \quad \phi=\gamma\phi_1 + (1-\gamma) \phi_2 \end{equation*}\] \[\begin{equation*} \phi_1=\frac{1}{W}\int_{\Omega_{\rm B}}\left|E^{\rm ill} - E^{\rm ref}\right|^2 {\rm d} \B r\end{equation*}\] \[\begin{equation*} \phi_2=\int_{\Omega_{\rm D}}\frac{h_{\rm m}^2}{S}\left| \B \nabla\,\varepsilon \right|^2 {\rm d} \B r\end{equation*}\]

• \(W=\int_{\Omega_{\rm B}}\left|E^{\rm ref}\right|^2{\rm d} \B r\) (reference object only)

Illusion

¹¹

Experimental setup

Experimental results

Experimental results

FEM: implementation

Open source libraries with bindings for the python programming language using a custom code gyptis.

• Geometry and mesh generation: gmsh
• FEM library: fenics using second order Lagrange basis functions
• Gradient calculations: dolfin-adjoint library with automatic differentiation

^{1213 1415}

Bi-focal lens

Objective: focal point at two different locations depending on the excitation frequency \[\begin{equation*} \max_{p(\B r)} \quad \Phi = \left|E_1(\omega_1,\B r_1)\right| + \left|E_2(\omega_2,\B r_2)\right| \end{equation*}\]

Optimization history

Bi-focal lens

¹⁶

Measurements

Measurements

Ferrorelectric metamaterials

¹⁷

Nonlinear permittivity

¹⁸

Nonlinear permittivity coupled with electrostatics

Solve for the potential \(V\) satisfying law for a given permittivity distribution \(\e\): \[\begin{equation} \ddiv (\e(\B E) \grad V) = 0 \label{eq:elstat}\end{equation}\] biased by a constant uniform electric field \(\B E_{\rm B} =\Eb \B x\). The electric field \(\B E=-\grad V\) derived from the solution depends on the permittivity distribution, which itself depends on the electric field.

Homogenized properties

Two scale expansion \(\Hm(\B r)=\Hm_0(\B r,\B r/\nu) + \nu \Hm_1(\B r,\B r/\nu) + \nu^2 \Hm_2(\B r,\B r/\nu) + ...\) for magnetic field \(\Hm\)

Two annex problems \(\mathcal P_j\), \(j={x, y}\) on the unit cell : \[\begin{equation} \ddiv \left[ \e(\B E) \grad(\psi_j + r_j) \right] = 0, \label{eq:hom_annex}\end{equation}\] The homogenized permittivity \(\ehom\) is then obtained as: \[\begin{equation} \ehom(\B E) = \langle \e(\B E)\rangle \B I + \B \phi(\B E), \label{eq:hom_epsi}\end{equation}\] where \(\mn{.}\) denotes the mean value over the unit cell, \(\B I\) is the \(2\times 2\) identity matrix and \(\B \phi_{ij}(\B E) = \langle \e(\B E) \grad \psi_i(\B E) \rangle_j\) are correction terms.

¹⁹

Enhancement of tunability

Optimized ferroelectric metamaterials

Objective: maximize the tunability of the composites along the direction of the applied electric field defined as \(\tunh(E)=\exxh(0)/\exxh(E)\) \[\begin{equation} \begin{aligned} \max_{p(\B r)} \quad & \sigma(E) = \frac{\tunh(E)}{\eta(E)}-1\\ \textrm{s.t.} \quad & f_{\rm min} < f < f_{\rm max}\\ \end{aligned} \label{eq:max_tuna} \end{equation}\] where \(\sigma\) is the tunability enhancement compared to the bulk BST tunability \(\eta(E)=\ef(0)/\ef(E)\).

²⁰

Optimized ferroelectric metamaterials

Optimized ferroelectric metamaterials

Bi-objective optimization

Loss reduction factor: \[\begin{equation} \theta(E) = 1- \frac{\tdhxx(E)}{\tdf(E)} \label{eq:loss_red_factor} \end{equation}\] where the homogenized loss tangent \(xx\) component is \(\tdhxx(E) = -\im \exxh(E) / \re\exxh(E)\).

Objective: high tunability and low loss \[\begin{equation} \begin{aligned} \min_{p(\B r)} \quad & \gamma\,\sigma(E) + (1-\gamma)\,\theta(E)\\ \textrm{s.t.} \quad & f_{\rm min} < f < f_{\rm max}\ \end{aligned} \label{eq:biobj} \end{equation}\]

Bi-objective optimization

Bi-objective optimization

Inverse design of superscatterers

• Design domain: circular rod of diameter \(D=2R=\lambda_0=600\)nm, materials: \(\rm SiO_2\) (\(\varepsilon_{\rm min}=2.17\)) and silicon (\(\varepsilon_{\rm max}=15.59-0.22i\)).
• Scattering width computed on a circle \(\Gamma\) enclosing the object: \[\begin{equation} \sigma_s=\frac{1}{|\B S_i|}\int_\Gamma \B n \cdotp \B S_s {\rm d} \Gamma \label{eq:scs} \end{equation}\] with Poynting vectors \(\B S_i = \frac{1}{2} {\rm Re}[\B E_i \times \B H_i^\star]\) and \(\B S_s = \frac{1}{2} {\rm Re}[\B E_s \times \B H_s^\star]\) associated with the incident and scattered field respectively, \(\B n\) is the unit vector normal to \(\Gamma\).

Objective: maximize the normalized scattering width: \[\begin{equation} \max_{p(\B r)} \quad \Phi = \sigma_s/2R \end{equation}\]

History

Superscatterers: results

Superscatterers: modal analysis

FMM

• Rescaled magnetic field \(\tilde{\B H} = \sqrt{\mu_0/\varepsilon_0}\B H\) expanded as: \[\begin{equation} \tilde{\B H}(\rpara, z)=\sum_{\B {G}} \B{H}_{\B {G}}(z) {\rm e}^{i(\B {k}+\B {G}) \cdot \rpara}, \label{eq:Hexpansion} \end{equation}\] with \(\B {G} = n L_{k,1} + m L_{k,2}\) a reciprocal lattice vector and \(\rpara\) the in plane component of the position vector
• Fourier transform: \[ \varepsilon_{\B{G}}=\frac{1}{\left|L_{r}\right|} \int_{\Omega} \varepsilon(\B{r}) e^{-i \B{G} \cdot \B{r}} d \B{r} \] where \(\Omega\) denotes one unit cell of the lattice
• Form block Toepliz matrices \(\hat{\varepsilon}_{m n}=\varepsilon_{\left(\B{G}_{m}-\B{G}_{n}\right)}\) for each components of \(\B \varepsilon\).

FMM

Fourier coefficients \(d_{x}\) and \(d_{y}\) of the displacement field \(\B{D}\) related to the electric field by: \[\begin{equation} \left[\begin{array}{c}-d_{y}(z) \\ d_{x}(z)\end{array}\right]=\mathcal{E}\left[\begin{array}{c}-e_{y}(z) \\ e_{x}(z)\end{array}\right] \label{eq:DE} \end{equation}\]

• Early formulation:²³ Laurent’s rule \[\mathcal{E}=\left[\begin{array}{ll}\hat{\epsilon}_{x x} & \hat{\epsilon}_{x y} \\ \hat{\epsilon}_{y x} & \hat{\epsilon}_{y y}\end{array}\right].\]
• Improved convergence: proper factorization Rules (Li, 1997):²⁴ decompose the electric field in its normal and tangential components at material interfaces (normal vector field generation)²⁵

FMM

• Anzatz for eigenmode: \[\begin{equation*}\B{H}(z)=\sum_{\B{G}}\left[\phi_{\B{G}, x} \B{x}+\phi_{\B{G}, y} \B{y}-\frac{\left(k_{x}+G_{x}\right) \phi_{\B{G}, x}+\left(k_{y}+G_{y}\right) \phi_{\B{G}, y}}{q} \B{z}\right] e^{i(\B{k}+\B{G}) \cdot \B{r}+i q z}\end{equation*}\]
• Fourier transform + eliminating \(z\) components, we get the EVP: \[\begin{equation*} M \Phi = \left(\mathcal{E}\left(k_0^{2}-\mathcal{K}\right)-K\right) \Phi= Q^{2}\Phi, \quad \Phi=\left[\begin{array}{l} \phi_{x} \\ \phi_{y} \end{array}\right] \end{equation*}\] \(Q^{2}={\rm diag\,} q_n^2\), eigenvalues \(q_{n}^2\), \(\left[\phi_{x, n}, \phi_{y, n}\right]^{\rm T}\) Fourier coefficients of the eigenmodes \[\begin{equation*} \mathcal{K}=\left[\begin{array}{cc}\hat{k}_{y} \hat{\varepsilon}_{z}^{-1} \hat{k}_{y} & -\hat{k}_{y} \hat{\varepsilon}_{z}^{-1} \hat{k}_{x} \\ -\hat{k}_{x} \hat{\varepsilon}_{z}^{-1} \hat{k}_{y} & \hat{k}_{x} \hat{\varepsilon}_{z}^{-1} \hat{k}_{x}\end{array}\right], \quad K=\left[\begin{array}{ll}\hat{k}_{x} \hat{k}_{x} & \hat{k}_{x} \hat{k}_{y} \\ \hat{k}_{y} \hat{k}_{x} & \hat{k}_{y} \hat{k}_{y}\end{array}\right], \end{equation*}\] where \(\hat{k}_{\nu}\), \(\nu\in\{x,y\}\), are diagonal matrices with entries \(\left(k_{\nu}+G_{1 \nu}, k_{\nu}+G_{2 \nu}, \ldots\right)\).

FMM

• Compute fields in each layer
• \(S\)-matrix algorithm is then used starting from the incident medium to recursively find the \(S\)-matrix of each layer and form the total \(S\)-matrix.
• Computation of relevant electromagnetic quantities of interest, e.g. diffraction efficiencies in transmission and reflection for each order by calculating the power flux through the outermost layer (more efficient in Fourier space)²⁶

Implementation

FMM and PWEM in python with various numerical backends for core linear algebra operations and array manipulation

• numpy²⁷
• scipy²⁸
• autograd (AD)²⁹
• pytorch (AD + GPU)³⁰
• jax (AD + GPU)³¹

FMM benchmark

Deflective metasurface

• \(\lambda_t=732\)nm, periods \(L_x=800\)nm and \(L_y=400\)nm, plane wave normally incident from the superstrate (silica, \(\varepsilon=2.16\)).
• Design domain: metasurface layer of thickness \(350\)nm, restricted to one unit cell, \(\varepsilon_{\rm min}=1\) (air) and \(\varepsilon_{\rm max}=14.06-0.074i\) (silicon)

Objective: maximize the average of the transmission coefficient in the \((1,0)\) diffracted order for both polarizations: \[\begin{equation} \max_{p(\B r)} \quad \Phi = \frac{1}{2} \left( T^{\rm TE}_{(1,0)} + T^{\rm TM}_{(1,0)}\right) \end{equation}\]

Photonic crystals: maximizing bandgaps

• TE modes, square array with enforced \(C_4\) symmetry on the unit cell, \(\varepsilon_{\rm min}=1\) (air) and \(\varepsilon_{\rm max}=9\)

Objective: open and maximize a bandgap between the \(5^{th}\) and \(6^{th}\) eigenvalues:\[\begin{equation*}\max_{p(\B r)} \quad \Phi = \min_{\B k} \omega_{6}(\B k) - \max_{\B k} \omega_{5}(\B k)\end{equation*}\]

• Final distribution in agreement with simple geometrical rules: the walls of an optimal centroidal Voronoi tessellation with \(n=5\) points³³

Photonic crystals: dispersion engineering

• TM modes, symmetry with respect to \(y\)

Objective: obtain a prescribed dispersion curve for the \(6^{th}\) band \[\begin{equation*}\min_{p(\B r)} \quad \Phi = \left\langle\left|\omega_{6}(k_x) - \langle \omega_{6}\rangle - \omega_{\rm tar}(k_x) \right|^2\right\rangle \end{equation*}\] with \[\begin{align*}\omega_{\rm tar}(k_x) =& -0.02 \cos(k_x a) + 0.01 \cos(2 k_x a) \\ &+ 0.007 \cos(3 k_x a)\end{align*}\] \(\langle f\rangle =\frac{1}{M}\sum_{m=0}^M f_m\)

Get the code

• Development on gitlab: continuous integration for testing and documentation deployment
• Install / fork it / run it online / report bugs!

FEM

pip install gyptis

conda install -c conda-forge gyptis

FMM

pip install nannos

conda install -c conda-forge nannos

PWEM

pip install protis