1. Introduction

The design of optical elements made of quasi-periodic arrays of sub-wavelength scatterers, also known as metasurfaces, is a promising area of research. The miniatuarization of existing optical elements such as lenses [1–5], freeform optics [6], and retroreflectors [7] has already been shown using metasurfaces. Furthermore, multi-functional optical elements [8–14] and new point spread function engineering methods [15,16] have been demonstrated using metasurfaces. Until recently, however, these metasurfaces have generally been designed intuitively, termed here as forward design. Libraries of complex transfer coefficients of individual scatterers are pre-computed, and arranged in a periodic lattice to approximate a desired phase profile. The properties of these scatterers are computed with periodic boundary conditions and the metasurfaces are designed under the "local phase approximation": the scattering in any small region is taken to be the same as the scattering from a periodic surface [1,2,17]. This approximation neglects inter-scatterer coupling which is significant for metasurfaces composed of scatterers with rapidly varying geometries or with low refractive index [18]. Moreover, it is not always possible to know the phase-profile a priori, and in these cases forward design methods cannot be used.

Inverse design methods use a figure of merit (FOM) written in terms of adjustable geometric scatterer parameters to iteratively optimize the scatterers of a metasurface to implement a desired functionality. The process starts with an arbitrary initial scatterer configuration. Then the electric field scattered off the device, the FOM, and the gradient of the FOM with respect to the scatterer design parameters are computed. The scatterer geometries are then iteratively updated in the direction that optimizes the FOM. Thus, inverse design methods offer a clear path to create optical elements with unintuitive phase functions. Different optimization methods such as particle swarm optimization [19], genetic algorithm [20–22], and gradient based methods [23–31] have already been applied to design both integrated photonic elements and free space metasurface optics. One specific direction is to exploit Mie scattering of spherical scatterers to perform the inverse design [23,32]. This approach allows large-area metasurface design without relying on the local phase approximation, and thus accurately models the inter-scatterer coupling. Currently, this approach is restricted to spherical scatterers, for which the radii are the only free parameters available. We did not find a radius range over which these spherical scatterers smoothly span a $0-2\pi$ phase shift without suffering considerable optical losses, a common requirement when designing metaphotonic structures.

In this work, we present an inverse-design and optimization method for large-area ($\sim$ 40$\lambda$ in diameter) metasurfaces using transition matrix scattering theory, an extension of Generalized Multi-sphere Mie Theory (GMMT). Specifically, we extended a previously reported inverse design method [23] to include ellipsoidal scatterers. We first show the feasibility of using this method for inverse-designing single wavelength metasurfaces lenses. We then demonstrate the effectiveness of this method for designing non-intuitive devices without a known phase function by optimizing a polarization multiplexed lens that switches the location of the focal spots based on the incident light polarization. Finally we demonstrate the efficacy of inverse design techniques for optimization, by improving the efficiency of a high numerical aperture lens starting with a forward designed metalens as the initial condition.

2. T-matrix formalism

In our design, we employ the T-Matrix Method (TMM) for ellipsoidal scatterers. Since a rigorous treatment of TMM for solving electromagnetic scattering can be found elsewhere [33–35], we provide give a brief overview of this method. The TMM formalism allows for a faster and less memory intensive forward simulation compared to direct methods of solving Maxwell’s equations such as finite difference methods, at the cost of a restricted scatterer geometry.

In the case of a single scatterer $S_i$, the net electric field is a sum of the incident and scattered fields: $E(\vec {r}) = E_{in}^i(\vec {r}) + E_{scat}^i(\vec {r})$ where $E_{in}$ and $E_{scat}$ can be written as series expansions of the incident and scattered fields in the spherical vector wave function (SVWF) basis:

 

Fig. 1. A. Mie scattering schematic. Light is incident onto the set of ellipsoidal scatterers. Each scatterer has an associated T-matrix. The incident field onto each scatterer is described by the incident field $E_{in}$ and the scattered fields from all other scatterers. The inter-particle coupling is represented by the matrix $W_{n''n'}^{i''i'}$ which describes the coupling between spheres $i'$ and $i''$. B. Design parameters for ellipsoidal scatterers. The semi-major axes are taken to always be aligned with the particle frame: semi-major axis $a$ is aligned with the $x_{part}$ axis, $b$ with the $y_{part}$ axis, and $c$ with the $z_{part}$ axis. The rotation $\phi$ is about the z-axis, with the counterclockwise direction defined as a positive rotation.

Download Full Size | PDF

3. Adjoint optimization

Our design process begins with a set of scatterer locations and geometric properties of each individual ellipsoid as initial conditions. During the optimization process, each individual scatterer geometry is iteratively modified. To calculate the gradient, we use the adjoint method [23–25,27,36]. In our previous work [23], we calculated the general gradient of a FOM with respect to individual sphere radii R. Using a similar approach, we calculate the gradient of the FOM with respect to the free parameters of an elliptical scatterer (Fig. 1B). Given a set of design parameters $\{\textbf {P}\}$, we can write a FOM $f(\textbf {b(P)},\textbf {P})$, where b is the vector containing coefficients $b_n^i$. We want to calculate the FOM with respect to parameters $\textbf {P}$. The procedure of calculating the gradient with respect to the free parameters of the ellipsoidal scatterer is identical to that of a spherical scatterer, so we refer to Ref. [23] and write

Since, to the best of our knowledge, this is the first time these gradients have been calculated, we first numerically verify their validity. We denote the T-matrix of an ellipsoidal scatterer $S$ as a function of its geometry

 

Fig. 2. Verification of the analytical T-matrix derivatives. A shows the error between the analytical T-matrix derivative and the numerical derivative with respect to semi-major axis $a$, B with respect to $b$, C with respect to $c$. Figure 2D shows the T-matrix derivative with respect to the azimuthal rotation of the ellipsoid $\phi$. As the step size of the numerical approximation to the derivative gets smaller, the mean error between numerical and analytical derivative gets closer to 0, which implies that the analytical derivatives are valid.

Download Full Size | PDF

4. Inverse design of optical elements

Using the aforementioned TMM formalism and adjoint optimization method, we present the design and optimization of two optical elements: a lens with numerical aperture (NA) of $\sim 0.83$, and a lens that switches focal lengths based on the polarization of light incident onto the lens. To design these devices, we must specify a FOM that encompasses its performance. In both cases, the lenses were designed for $915 nm$ incident wavelength, with scatterers having refractive index of 3.56 surrounded by vacuum. Each device was designed by starting off with identical ellipsoids, and minimizing a specified FOM.

4.1 High numerical aperture lens

Figure 3A shows the distribution of ellipsoidal scatterers for a high NA lens with a diameter of 30 $\mu m$ and a focal length of 10 $\mu m$. The FOM for the lens can be written as:

 

Fig. 3. A Final distribution of scatterers with periodicity $450$ nm for the inverse designed lens. semi-major axes $a$ and $b$ are allowed to range between $40$ and $150$ nm. Semi-major axis $c$ is allowed to range between $40$ and $300$ nm. B the field cross-section in the x-z plane at $y=0 \mu m$, C the cross-section in the x-y plane at $z=10\mu m$. D shows the Gaussian fit to the field at the focal spot $z=10 \mu m$ along $x=0$. In order to calculate the lens efficiency, the full-width at half-maximum (FWHM) was calculated for the fitted Gaussians. The integral of the field intensity around the disk $d=3 \times FWHM$ about the center of the focal spot was calculated, and then divided by the total incident field intensity. The units of all plots are arbitrary light intensity units. The efficiency of the inverse designed lens was calculated to be $3.38\%$.

Download Full Size | PDF

Here, $I^T$ is some arbitrary intensity value at the focal spot of the lens, and $I^A$ is the actual intensity at that spot calculated via TMM. For this problem, we want to minimize the FOM over parameters $a,b,c$. We choose to optimize only over these parameters as we found very little dependence of the lens performance on the scatterer rotation. We initialize a grid of identical ellipsoidal scatterers with nominal semi-major axis radii $a=b=100nm$ and $c=300nm$, and a lattice periodicity of $450nm$. The radii are allowed to vary between $40 nm$ and $150 nm$ for axes $a$ and $b$, and between $40nm$ and $300 nm$ for semi-major axis $c$. These parameters were picked in such a way that the scatterer phase responses span a $0-2\pi$ phase shift over this parameter range. The maximum radii of the semi-major axes are chosen in such a way that a circumscribing sphere with the radius of the largest semi-major axis radius of one particle cannot overlap with the surface of any neighboring particle. This is due to limitations in our method, as we compute the inter-scatterer coupling by assuming the incident electric field onto each particle is composed of the incident field and the field scattered from the surfaces of spheres inscribing the ellipsoidal particles [35,37]. We also choose to cutoff our field expansions at the multiple order of $l$=3 [23].

Figure 3 shows the performance of the designed lens. Figure 3B shows that there is a clear focal spot at $10 \mu m$. All electric field in this work were calculated by using our extension of the CELES code [35]. The efficiency of the lens was calculated by fitting a Gaussian shape to the field profile at the focal spot $z=10 \mu m$, for $x=0$ as shown in Fig. 3D. Then we found the full-width at half maximum of the Gaussian, and integrated the intensity of the field at that focal spot, and divided it over the total intensity of the incident light. This quantity is defined as the efficiency of the lens $\eta$ and for a lens with focal length $F$ is defined as:

 

Fig. 4. A Scatterer distribution of the polarization multiplexed lens. Lattice periodicity is 650 nm, radii were limited to range from 40 nm to 292.5 nm for the $a$ and $b$ axes, and $0$ to $357.5 nm$ for the c axis. For the initial condition, all of the semi-major axis radii were set to $250 nm$, and the rotations were set to $0$ radians. In the final parameter distribution, the scatterers look very similar, and indeed, the minimum semi-major axis radius in the design is $\sim 205 nm$ and the maximum is $\sim 289 nm$. B-D Are field distributions correspond to x-polarized light, and E-G correspond to y-polarized light. B,E are scattered field slices in the x-z plane at $y=0\mu m$. C,F are x-y profiles at each focal spot. C is a slice at $z=20 \mu m$, and E is a slice at $z=30 \mu m$. D,G are Gaussian fits at each focal spot.

Download Full Size | PDF

4.2 Inverse design of polarization switched focal length lens

Then we designed a lens with a diameter of 40 $\mu m$, and focal lengths of 20 $\mu m$ (NA $\sim 0.71$) and 30 $\mu m$ (NA $\sim 0.55$) for the x and y polarizations respectively. The lattice constant for this lens was taken to be $650 nm$. Semi-major axis radii $a$ and $b$ were allowed to range between $40 nm$ and $292.5 nm$. Semi-major axis radius $c$ was allowed to range between $40 nm$ and $357.5 nm$. The azimuthal rotation around the $z$ axis of the scatterers, was allowed to range from $-\pi /2$ to $\pi /2$.

The optimization problem was framed as a min-max optimization problem [29]. For this optimization, we write the total FOM as a sum of FOM’s for each polarization, given by

We calculated the efficiency of this lens for each polarization using the method described in the previous section, and found values of $\eta = 2.31\%$ for the x-polarization and $\eta = 3.38 \%$ for y polarization. Other relevant quantities to characterize the performance of this device are the contrast ratios of the focal spots. We define and report two different contrast ratios. The first one is the ratio between the intensity at the focal spot, where light should be maximized, to the intensity at the focal spot of the orthogonal polarization. We found these ratio values to be $\frac {I(0,0,30\mu m)}{I(0,0,20 \mu m} = 8.75$ for y polarized light and $\frac {I(0,0,20\mu m)}{I(0,0,30 \mu m} = 5.11$ for x polarized light. The second one is the ratio between the intensity at the focal spot for both polarizations. We found these values to be $\frac {I_x(0,0,20\mu m}{I_y(0,0,20)} = 5.58$ and $\frac {I_y(0,0,30\mu m}{I_x(0,0,30)} = 5.92$.

5. Metasurface lens optimization

Finally, we discuss the optimization of a forward designed metalens using our inverse design method. By using rigorous coupled wave analysis (RCWA), we computed the phase and amplitude response of a library of ellipsoidal scatterers with periodic boundary conditions [38]. The ellipsoidal scatterers we chose for this design have identical range of geometric and material properties to those described in section 4.1. Then we discretized the design space in the x-y plane, using a scatterer periodicity of 450 nm, and by using the phase equation for a lens given by

 

Fig. 5. Figs. A-D correspond to the forward designed lens, and Figs. E-H to the optimized lens. A,E are the scatterer distributions. E,F are the x-z slices of the resulting field profile at $y=0\mu m$. C,G correspond to the x-y field slice at $z=10\mu m$. D,H are the Gaussians fitted to the field profiles at their focal spot with $y=0\mu m$. The forward design lens efficiency was determined to be $25.59\%$, and the optimized efficiency was calculated to be $32.00\%$

Download Full Size | PDF

To optimize this device, we started off with the scatterer distribution given by the forward design as the initial condition and maximize the light intensity at a the focal spot. The performance of the optimized device is summarized in Fig. 5. The efficiency of this device was calculated to be 32.00%, which is a 6.41% improvement over the forward design lens. On average, each individual scatterer was changed by approximately $3.03 nm$ along the $a$ axis, $4.8 nm$ along the $b$ axis, and $0.17 nm$ along the c axis. The standard deviations for each axis are $3.53 nm$, $4.89 nm$, $0.42 nm$ respectively. The maximum changes for each axis were $33.42 nm$, $33.43 nm$ and $5.57 nm$ respectively. It is worth noting that this improvement implies that lenses designed by the conventional forward design methods are not necessarily globally optimal, even for high index materials. We can also see that the initial conditions are very important for the final design, as starting with identical ellipsoids, the final design provides very low efficiency. In fact, based on our analysis, we believe that our inverse design method will be more suitable for optimization type of problem, where the initial conditions are developed based on intuition and prior knowledge.

6. Discussion

We demonstrated a new optimization method for designing large area dielectric metasurfaces made of ellipsoidal scatterers based on the adjoint method and a generalization of GMMT. Starting from an array of identical ellipsoidal scatterers, we designed a high NA ($\sim 0.83$) lens and a polarization multiplexed lens that focuses light at $30\mu m$ and $20 \mu m$ based on the polarization of incident light. Although polarization was the parameter that we chose to optimize over, the same approach can be used to design angle or wavelength multiplexed devices. Furthermore, by modifying the FOM to finely sample the wavelength or angle of the incident light, broad-band and broad-angle devices can be designed at the expense of simulation time. We have also shown that starting with a forward-designed lens as an initial condition, a higher efficiency design can be obtained via optimization.

We note that all the reported devices were designed at refractive index $n = 3.56$. As our method requires the bounding spheres of ellipsoidal scatterers not to overlap other ellipsoidal scatterers, we are limited by the aspect ratio and density of the ellipsoids, and only with high index ellipsoids we can maintain low density of scatterers while spanning the whole $0-2\pi$ phase. Unfortunately, there is currently no straightforward way to fabricate these structures with such a high index. One solution could be to use a high index resin in additive manufacturing [39]. It is also possible to fabricate cylindrical scatterers at high refractive indices by using traditional lithography. This would require a further generalization of the T-Matrix method to expand the incident and scattered fields in terms of spheroidal wave-functions instead of SVWF or by using the plane wave coupling method [40].

Appendix A. T-matrix derivatives

The T-Matrix of the ellipsoid requires the computation of the $Q$ and $RgQ$ matrices that represent the coupling between the scattered field and the incident field to the surface fields respectively. It depends only on the geometric and material properties of particle itself, and the wavelength of excitation. The T-matrix is then given by [41–43]:

(18)$$T = RgQ (Q)^{-1}$$

The $Q$ matrix is square, and is composed of four square submatrices $\bar {P}$, $\bar {R}$, $\bar {S}$, and $\bar {U}$, given by:

(19)$$Q = \begin{bmatrix} \bar{P} & \bar{R} \\ \bar{S} & \bar{U} \end{bmatrix},$$

These individual square matrices are given by [37]:

(20)$$\bar{P}_{lml'm'} = -i k k_s J^{(21)}_{lml'm'} - i k^2 J^{(12)}_{lml'm'},$$

(21)$$\bar{R}_{lml'm'} = -i k k_s J^{(11)}_{lml'm'} - i k^2 J^{(22)}_{lml'm'},$$

(22)$$\bar{S}_{lml'm'} = -i k k_s J^{(22)}_{lml'm'} - i k^2 J^{(11)}_{lml'm'},$$

(23)$$\bar{U}_{lml'm'} = -i k k_s J^{(12)}_{lml'm'} - i k^2 J^{(21)}_{lml'm'},$$

where the $J$ terms represent integrals over the surface of the particle, and are given by:

(24)$$J^{(pq)}_{lml'm'} = (-1)^m \int_S dS \hat{n}(\textbf{r}) \cdot \boldsymbol{\Psi}_{p,l',m'}^{(1)} (k_s r, \theta,\phi) \times \boldsymbol{\Psi}_{q,l,-m}^{(3)}(kr,\theta,\phi),$$

where $S$ is the surface bounding the particle, $dS$ is infinitesimal surface area, and $\hat {n}$ is a outward pointing unit normal at $dS$. The SVWFs $\boldsymbol{\Psi }^{(1)}$ and $\boldsymbol{\Psi }^{(3)}$ are given by [42]:

(25)$$\boldsymbol{\Psi}^{(\nu)}_{1lm}(\textbf{r}) = \frac{e^{i m \phi}}{\sqrt{2l(l+1)}}b_l(kr)\big[im\pi_{lm}(\theta)\hat{\theta}-\tau_{lm}(\theta)\hat{\phi}\big],$$

(26)$$\begin{aligned} \boldsymbol{\Psi}^{(\nu)}_{2lm}(\textbf{r}) = & \frac{e^{i m \phi}}{\sqrt{2l(l+1)}}\bigg\{l(l+1)\frac{b_l(kr)}{kr}P^{|m|}_l(cos\theta)\hat{r} \\ & +\frac{1}{kr}\frac{\partial (kr b_l(kr))}{\partial (kr)}\big[\tau_{lm}(\theta)\hat{\theta}+im\pi_{lm}(\theta)\hat{\phi}\big]\bigg\}, \end{aligned}$$

Here we have defined:

(27)$$\pi_{lm}(\theta) = \frac{P^{|m|}_l(cos\theta)}{sin\theta}, \tau_{lm}(\theta) = \frac{\partial P^{|m|}_l(cos\theta)}{\partial \theta}.$$

$P^m_l(x)$ is the associated Legendre polynomial. $j_l$ is the spherical Bessel function of order $l$, and $b_l$ is either a spherical Bessel function ($j_l$) for $\nu = 1$ or spherical Hankel function of the first kind ($h^{(1)}_l$) of order $l$ for $\nu = 3$, depending on whether $RgQ$ or $Q$ is being computed. In spherical coordinates, the product of the unit normal and the infinitesimal area is:

(28)$$dS \hat{n}(\textbf{r}) = r^2 sin(\theta)\boldsymbol{\sigma}(\textbf{r}) d\theta d\phi,$$

and $\boldsymbol{\sigma }$ is given by:

(29)$$\boldsymbol{\sigma}(\textbf{r}) = \hat{r} - \hat{\theta}\frac{1}{r}\frac{\partial r}{\partial \theta} - \hat{\phi}\frac{1}{r sin\theta}\frac{\partial r}{\partial \theta}.$$

In this case, $r$ is parameterizing the surface of a particle, and for an ellipsoid in spherical coordinates, $r$ is given by:

(30)$$r(\theta,\phi) = \left[ sin^2\theta \left(\frac{cos^2\phi}{a^2}+\frac{sin^2\phi}{b^2}\right) + \frac{cos^2\theta}{c^2}\right]^{-1/2}$$

To compute $RgQ$ rather than $Q$, we simply need to replace $\boldsymbol{\Psi }^{(3)}$ in the $J$ integrals with $\boldsymbol{\Psi }^{(1)}$.

The derivative of the T-matrix of a particle with respect to some parameter $p$ is given by:

(31)$$\frac{\partial T}{\partial p} = \left(\frac{\partial RgQ}{\partial p} - T \frac{\partial Q}{\partial p}\right)Q^{-1},$$

Hence we need to find the derivatives of the sub-matrices $\bar {P}$, $\bar {R}$, $\bar {S}$, and $\bar {U}$ with respect to $p$. This requires us to take the derivatives of the surface integrals $J$ with respect to the parameter $p$. In general, our parameter of interest $p$ will be some geometric quantity that determines the shape of the surface of integration $S$. In the specific case of ellipsoidal scatterers, they will be the three independent axes $a$, $b$, and $c$ along the $x$, $y$, and $z$ axes respectively.

The expressions for the derivatives with respect to a spatial variable ($a$, $b$, $c$) are as follows where $p$ represents any of the ellipsoid axes:

(32)$$\begin{aligned} \frac{\partial J^{(11)}_{lml'm'}}{\partial p} = & -i \iint \alpha_{lml'm'} \left( m'\pi_{l'm'} \tau_{lm}+m \pi_{lm} \tau_{l'm'}\right) \\ & \left[r\left(k\frac{\partial b_l}{\partial p} j_{l'} + k_s b_l \frac{\partial j_{l'}}{\partial p}\right)+2b_l j_{l'}\right]r sin\theta d\theta d\phi, \end{aligned}$$

(33)$$\begin{aligned} \frac{\partial J^{(12)}_{lml'm'}}{\partial p} = & \iint \alpha_{lml'm'} \bigg\{\bigg[ \frac{\partial R^{(12)}_{lml'm'}}{\partial r}+(\Theta^{(12)}_{lml'm'} E_\theta+\Phi^{(12)}_{lml'm'} E_\phi)\frac{\partial \rho_{l,l'}}{\partial r}\bigg] \frac{\partial r}{\partial p} \\ & +\bigg( \Theta^{(12)}_{lml'm'} \frac{\partial E_\theta}{\partial p}+\Phi^{(12)}_{lml'm'}\frac{\partial E_\phi}{\partial p}\bigg) \rho_{l,l'}\bigg\}d\theta d\phi, \end{aligned}$$

(34)$$\begin{aligned} \frac{\partial J^{(21)}_{lml'm'}}{\partial p} = & \iint \alpha_{lml'm'} \bigg\{\bigg[ \frac{\partial R^{(21)}_{lml'm'}}{\partial r}+(\Theta^{(21)}_{lml'm'} E_\theta+\Phi^{(21)}_{lml'm'} E_\phi)\frac{\partial \rho_{l,l'}}{\partial r}\bigg] \frac{\partial r}{\partial p} \\ & +\bigg( \Theta^{(21)}_{lml'm'} \frac{\partial E_\theta}{\partial p}+\Phi^{(21)}_{lml'm'}\frac{\partial E_\phi}{\partial p}\bigg) \rho_{l,l'}\bigg\}d\theta d\phi, \end{aligned}$$

(35)$$\begin{aligned} \frac{\partial J^{(22)}_{lml'm'}}{\partial p} = & \iint \alpha_{lml'm'} \bigg\{\bigg[\frac{\partial R^{(22)}_{lml'm'}}{\partial r}+\frac{\partial \Theta^{(22)}_{lml'm'}}{\partial r}E_\theta+\frac{\partial \Phi^{(22)}_{lml'm'}}{\partial r}E_\phi\bigg]\frac{\partial r}{\partial p} \\ & +\Theta^{(22)}_{lml'm'} \frac{\partial E_\theta}{\partial p}+\Phi^{(22)}_{lml'm'} \frac{\partial E_\phi}{\partial p} \bigg\} d\theta d\phi, \end{aligned}$$

where we have defined:

(36)$$\alpha_{lml'm'} = \frac{(-1)^m(1+(-1)^{m'-m})(1+(-1)^{l'+l+1}}{2\sqrt{l(l+1)l'(l'+1)}}e^{i(m'-m)\phi}$$

$k$ and $k_s$ are the k vectors of light in the medium surrounding the particle, and in the particle itself. Then we define:

(37)$$E_\theta = \frac{cos^2\phi}{a^2} + \frac{sin^2\phi}{b^2} - \frac{1}{c^2}$$

(38)$$E_\phi = \frac{1}{b^2}-\frac{1}{a^2},$$

(39)$$\rho_{l,l'} = r^3 j_{l'}b_l,$$

Now, we can define the specific terms used to construct each $J$ surface integral. For $J^{(12)}$, we define:

(40)$$ \frac{\partial R^{(12)}_{lml'm'}}{\partial r} = \frac{sin\theta}{k} (mm' \pi_{l'm'}\pi_{lm}+\tau_{l'm'}\tau_{lm})\bigg(j_{l'}\frac{\partial (kr b_l)}{\partial (kr)}\\+r\bigg(k_s \frac{\partial j_{l'}}{\partial r}\frac{\partial (kr b_l)}{\partial (kr)}+k j_{l'} \frac{\partial}{\partial r}\left(\frac{\partial (kr b_l)}{\partial (kr)}\right)\bigg)\bigg), $$

(41)$$\Theta^{(12)}_{lml'm'} = -\frac{sin\theta}{k}l(l+1)P^{|m|}_l\tau_{l'm'},$$

(42)$$\Phi^{(12)}_{lml'm'} = -i\frac{sin\theta}{k} l(l+1)m'P^{|m|}_l\pi_{l'm'}.$$

For $J^{(21)}$, we define:

(43)$$ \frac{\partial R^{(21)}_{lml'm'}}{\partial r} = - \frac{sin\theta}{k_s} (mm' \pi_{l'm'}\pi_{lm}+\tau_{l'm'}\tau_{lm})\bigg(\frac{\partial (k_s r j_{l'})}{\partial(k_s r)}b_l \\+r\bigg(k_s \frac{\partial}{\partial r}\left(\frac{\partial (k_s r j_{l'})}{\partial (k_s r)}\right)b_l + k \frac{\partial b_l}{\partial r}\frac{\partial (k_s r j_{l'})}{\partial(k_s r)}\bigg)\bigg), $$

(44)$$\Theta^{(21)}_{lml'm'} = \frac{sin\theta}{k_s}l'(l'+1)P^{|m'|}_{l'}\tau_{lm},$$

(45)$$\Phi^{(21)}_{lml'm'} = -i\frac{sin\theta}{k_s} l'(l'+1)mP^{|m'|}_{l'}\pi_{lm}.$$

Finally, for $J^{(22)}$ we define:

(46)$$\begin{aligned} \Theta^{(22)}_{lml'm'}= & i \frac{r^2 sin\theta}{k k_s} \bigg(m'l(l+1) \frac{\partial (k_s r j_{l'})}{\partial (k_s r)} b_l P^{|m|}_l \pi_{l'm'} \\ & +ml'(l'+1)j_{l'}\frac{\partial (kr b_l)}{\partial (kr)}P^{|m'|}_{l'}\pi_{lm}\bigg) \end{aligned}$$

(47)$$\begin{aligned} \Phi^{(22)}_{lml'm'} = & \frac{r^2 sin\theta}{k k_s}\bigg(l'(l'+1)j_{l'}P^{|m'|}_{l'}\frac{\partial (kr b_l)}{\partial (kr)}\tau_{lm} \\ & -l(l+1)\frac{\partial (k_s r j_{l'})}{\partial (k_s r)}\tau_{l'm'}b_lP_{lm}\bigg) \end{aligned}$$

and the three derivative terms:

(48)$$\begin{aligned} \frac{\partial R^{(22)}_{lml'm'}}{\partial r} = & -i\frac{sin\theta}{k k_s}( m'\pi_{l'm'} \tau_{lm}+m \pi_{lm} \tau_{l'm'}) \\ & \bigg(k \frac{\partial}{\partial r}\left(\frac{\partial (kr b_l)}{\partial (kr)}\right)\frac{\partial (k_s r j_{l'})}{\partial (k_s r)} \\ & + k_s\frac{\partial}{\partial r} \left(\frac{\partial (k_s r j_{l'})}{\partial (k_s r)}\right)\frac{\partial (kr b_l)}{\partial (kr)}\bigg) \end{aligned}$$

(49)$$\begin{aligned} \frac{\partial \Theta^{(22)}_{lml'm'}}{\partial r} = & i\frac{sin\theta}{k k_s}\bigg(ml'(l'+1)P^{|m'|}_{l'}\pi_{lm}\bigg(2r \frac{\partial (kr b_l)}{\partial (kr)} j_{l'} \\ & +r^2\bigg(k \frac{\partial}{\partial r}\left(\frac{\partial (kr b_l)}{\partial (kr)}\right) j_{l'}+k_s \frac{\partial j_{l'}}{\partial r} \frac{\partial (kr b_l)}{\partial (kr)}\bigg)\bigg) \\ & +m'l(l+1)P^{|m|}_l\tau{l'm'}\bigg(2rb_l\frac{\partial (k_s r j_{l'})}{\partial (k_s r)} \\ & +r^2\bigg(k \frac{\partial b_l}{\partial r}\frac{\partial (k_s r j_{l'})}{\partial (k_s r)}+k_sb_l\frac{\partial}{\partial r} \left(\frac{\partial (k_s r j_{l'})}{\partial (k_s r)}\right)\bigg)\bigg)\bigg) \end{aligned}$$

(50)$$\begin{aligned} \frac{\partial \Phi^{(22)}_{lml'm'}}{\partial r} = & \frac{sin\theta}{k k_s} \bigg(l'(l'+1)P^{|m'|}_{l'}\tau_{lm}\bigg(2r\frac{\partial (kr b_l)}{\partial (kr)}j_{l'} \\ & +r^2\bigg( \frac{\partial}{\partial r}\left(\frac{\partial (kr b_l)}{\partial (kr)}\right)j_{l'}+k_s \frac{\partial j_{l'}}{\partial r}\frac{\partial (kr b_l)}{\partial (kr)}\bigg)\bigg) \\ & -l(l+1)P^{|m|}_l\tau_{l'm'}\bigg(2r b_l \frac{\partial (k_s r j_{l'})}{\partial (k_s r)} \\ & +r^2\bigg(k\frac{\partial b_l}{\partial r}\frac{\partial (k_s r j_{l'})}{\partial (k_s r)}\frac{\partial}{\partial r} \left(\frac{\partial (k_s r j_{l'})}{\partial (k_s r)}\right)\bigg)\bigg)\bigg). \end{aligned}$$

Now with these $J$ integrals, we can compute the quantity $\frac {\partial T}{\partial p}$ for a given axis of an ellipsoid in its own particle frame where $a$, $b$, and $c$ are aligned along the $x_{part}$, $y_{part}$, and $z_{part}$ axes.

In addition to computing the response of the T-matrix of the ellipsoid to the contraction or extension of one of its axes, we are also interested in its response to rotations about the $z_{part}$ axis. To do this we will first define the transformation of the T-matrix or a derivative matrix from the particle frame to some rotated lab frame that has new axes $x_{lab}$ and $y_{lab}$, but shares $z_{lab} = z_{part}$. Given some rotation angle $\phi _{rot}$, we can then define our new axes:

(51a)$$ x_{lab} = x_{part}cos(\phi_{rot})+y_{part}sin(\phi_{rot}) $$

(51b)$$ y_{lab} = -x_{part}sin(\phi_{rot})+y_{part}cos(\phi_{rot}) $$

(51c)$$ z_{lab} = z_{part} $$

The general form of this orthogonal transformation in three dimensions can be represented by the Euler angles $\alpha$, $\beta$, and $\gamma$. The general transformation of an element of an operator $O$ from the particle frame to the lab frame can be written as [42]:

(52)$$O^{lab}_{plmp'l'm'}(\alpha,\beta,\gamma) = \sum^l_{m_1 = -l}\sum^{l'}_{m_2 = -l'}D^{l}_{mm_1}(\alpha,\beta,\gamma)O^{particle}_{plm_1p'l'm_2}D^{l'}_{m_2m'}(-\gamma,-\beta,-\alpha),$$

where the $D$ operator is a Wigner D-function. It can be represented as:

(53)$$D^l_{m'm}(\alpha,\beta,\gamma) = e^{-im'\alpha}d^l_{m'm}(\beta)e^{-im\gamma},$$

where $d^l_{m'm}(\beta )$ is Wigner’s (small) d-matrix given by:

(54)$$d^l_{m'm}(\beta) = \langle{l,m'}|e^{-i\beta J_y}|{l,m}\rangle.$$

However, as we are only concerned with rotations about the $z$ axis, we can simplify our expressions knowing that $\alpha$ is our only nonzero angle, and equation (56) becomes:

(55)$$O^{lab}_{plmp'l'm'}(\alpha,0,0) = \sum^l_{m_1 = -l}\sum^{l'}_{m_2 = -l'}D^{l}_{mm_1}(\alpha,0,0)O^{particle}_{plm_1,p'l'm_2}D^{l'}_{m_2m'}(0,0,-\alpha).$$

In this case, our $D$ operator has a much simplified form:

(56a)$$ D^{l}_{m'm}(\alpha,0,0) = e^{-im'\alpha}\delta_{m'm} $$

(56b)$$ D^{l}_{m'm}(0,0,\gamma) = e^{-im\gamma}\delta_{m'm}. $$

Combining equations (46), (47a), and (47b), we obtain a simple expression transforming $O$ from the particle frame to the lab frame:

(57)$$O^{lab}_{plmp'l'm'}(\alpha) = e^{i(m'-m)\alpha}O^{particle}_{plmp'l'm'}.$$

 

Fig. 6. Transmission of individual scatterers with periodic boundary conditions as a function of the radius of the ellipsoids (semi-major axes $a = b$). A is the plot of the complex transmission of the ellipsoids used from section 4.1 and 5. Ellipsoid height is fixed to be $600 nm$. The lattice constant is $450 nm$. B transmission response for ellipsoids outlined in section 4.2. Ellipsoidal height is fixed at $715nm$ with a lattice constant of $650nm$.

Download Full Size | PDF

Equation (61) is applicable to for transforming both T-matrices and the derivative matrices computed in the particle frame into the lab frame. It also gives us a prescription for computing the derivative matrix with respect to the particle’s angular orientation. We already have derivatives characterizing the response of the particle to contractions and extensions of its principal axes, and can now rotate these to a lab frame where the particle has an arbitrary angular orientation relative to the $z$ axis. We can now compute the derivative with respect to the particle’s angular orientation $\alpha$ as:

(58)$$\frac{\partial T^{lab}_{plmp'l'm'}(\alpha)}{\partial \alpha} = i(m'-m)e^{i(m'-m)\alpha}T^{particle}_{plmp'l'm'}.$$

With equations (31) and (58), we have characterized the derivatives of the T-matrix representing an ellipsoid with respect to its axes and orientation.

These integrals are implemented in MATLAB, and performed using Gaussian quadrature.

Appendix B. Scatterer electromagnetic response

Figure 6 shows the transmission of of individual scatterers with periodic boundary conditions as a function of the radius of the ellipsoids.

Appendix C. Machine specifications

Ubuntu 16.04

MATLAB v9.5.0 R2018b with parallel computing toolbox v2.4 2x Intel E5-2620 at 2.1 GHz

NVIDIA Tesla K40 12 GB Memory running CUDA 9.1

64 GB DDR3 Memory

Our inverse and forward design methods are solved using a modified version of CELES. More details about CELES are available from Ref. [35].

CELES is available free of charge, and our implementation of the optimization algorithm is available upon request.

Funding

NSF (1825308).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References

1. N. Yu, P. Genevet, M. A. Kats, F. Aieta, J.-P. Tetienne, F. Capasso, and Z. Gaburro, “Light propagation with phase discontinuities: Generalized laws of reflection and refraction,” Science 334(6054), 333–337 (2011). [CrossRef]  

2. S. Sun, Q. He, S. Xiao, Q. Xu, X. Li, and L. Zhou, “Gradient-index meta-surfaces as a bridge linking propagating waves and surface waves,” Nat. Mater. 11(5), 426–431 (2012). [CrossRef]  

3. W. T. Chen, A. Y. Zhu, V. Sanjeev, M. Khorasaninejad, Z. Shi, E. Lee, and F. Capasso, “A broadband achromatic metalens for focusing and imaging in the visible,” Nat. Nanotechnol. 13(3), 220–226 (2018). [CrossRef]  

4. A. Zhan, S. Colburn, R. Trivedi, T. K. Fryett, C. M. Dodson, and A. Majumdar, “Low-contrast dielectric metasurface optics,” ACS Photonics 3(2), 209–214 (2016). [CrossRef]  

5. A. Arbabi, Y. Horie, A. J. Ball, M. Bagheri, and A. Faraon, “Subwavelength-thick lenses with high numerical apertures and large efficiency based on high-contrast transmitarrays,” Nat. Commun. 6(1), 7069 (2015). [CrossRef]  

6. A. Zhan, S. Colburn, C. M. Dodson, and A. Majumdar, “Metasurface freeform nanophotonics,” Sci. Rep. 7(1), 1673 (2017). [CrossRef]  

7. A. Arbabi, E. Arbabi, Y. Horie, S. M. Kamali, and A. Faraon, “Planar metasurface retroreflector,” Nat. Photonics 11(7), 415–420 (2017). [CrossRef]  

8. A. Arbabi, Y. Horie, M. Bagheri, and A. Faraon, “Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission,” Nat. Nanotechnol. 10(11), 937–943 (2015). [CrossRef]  

9. S. M. Kamali, E. Arbabi, A. Arbabi, Y. Horie, M. Faraji-Dana, and A. Faraon, “Angle-multiplexed metasurfaces: Encoding independent wavefronts in a single metasurface under different illumination angles,” Phys. Rev. X 7(4), 041056 (2017). [CrossRef]  

10. S. Xiao, F. Zhong, H. Liu, S. Zhu, and J. Li, “Flexible coherent control of plasmonic spin-hall effect,” Nat. Commun. 6(1), 8360 (2015). [CrossRef]  

11. B. Yang, T. Liu, H. Guo, S. Xiao, and L. Zhou, “High-performance meta-devices based on multilayer meta-atoms: interplay between the number of layers and phase coverage,” Sci. Bull. 64(12), 823–835 (2019). [CrossRef]  

12. H. Guo, J. Lin, M. Qiu, J. Tian, Q. Wang, Y. Li, S. Sun, Q. He, S. Xiao, and L. Zhou, “Flat optical transparent window: mechanism and realization based on metasurfaces,” J. Phys. D: Appl. Phys. 51(7), 074001 (2018). [CrossRef]  

13. X. Shiyi, W. Jiarong, L. Fu, Z. Shuang, Y. Xiaobo, and L. Jensen, “Spin-dependent optics with metasurfaces,” Nanophotonics 6, 215–234 (2016). [CrossRef]  

14. S. Xiao, H. Mühlenbernd, G. Li, M. Kenney, F. Liu, T. Zentgraf, S. Zhang, and J. Li, “Helicity-preserving omnidirectional plasmonic mirror,” Adv. Opt. Mater. 4(5), 654–658 (2016). [CrossRef]  

15. S. Colburn, A. Zhan, and A. Majumdar, “Metasurface optics for full-color computational imaging,” Sci. Adv. 4(2), eaar2114 (2018). [CrossRef]  

16. S. Colburn and A. Majumdar, “Simultaneous varifocal and broadband achromatic computational imaging using quartic metasurfaces,” (2019).

17. L. Hsu, M. Dupré, A. Ndao, J. Yellowhair, and B. Kanté, “Local phase method for designing and optimizing metasurface devices,” Opt. Express 25(21), 24974–24982 (2017). [CrossRef]  

18. E. Bayati, A. Zhan, S. Colburn, M. V. Zhelyeznyakov, and A. Majumdar, “Role of refractive index in metalens performance,” Appl. Opt. 58(6), 1460–1466 (2019). [CrossRef]  

19. K. Wang, J. Zhao, Q. Cheng, D. S. Dong, and T. J. Cui, “Broadband and broad-angle low-scattering metasurface based on hybrid optimization algorithm,” Sci. Rep. 4(1), 5935 (2015). [CrossRef]  

20. F. Aieta, M. A. Kats, P. Genevet, and F. Capasso, “Multiwavelength achromatic metasurfaces by dispersive phase compensation,” Science 347(6228), 1342–1345 (2015). [CrossRef]  

21. V. Egorov, M. Eitan, and J. Scheuer, “Genetically optimized all-dielectric metasurfaces,” Opt. Express 25(3), 2583–2593 (2017). [CrossRef]  

22. M. Donelli, “Design of broadband metal nanosphere antenna arrays with a hybrid evolutionary algorithm,” Opt. Lett. 38(4), 401–403 (2013). [CrossRef]  

23. A. Zhan, T. K. Fryett, S. Colburn, and A. Majumdar, “Inverse design of optical elements based on arrays of dielectric spheres,” Appl. Opt. 57(6), 1437–1446 (2018). [CrossRef]  

24. P. Hansen and L. Hesselink, “Accurate adjoint design sensitivities for nano metal optics,” Opt. Express 23(18), 23899–23923 (2015). [CrossRef]  

25. A. Y. Piggott, J. Petykiewicz, L. Su, and J. Vuckovic, “Fabrication-constrained nanophotonic inverse design,” Sci. Rep. 7(1), 1786 (2017). [CrossRef]  

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

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

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

29. R. Pestourie, C. Pérez-Arancibia, Z. Lin, W. Shin, F. Capasso, and S. G. Johnson, “Inverse design of large-area metasurfaces,” Opt. Express 26(26), 33732–33747 (2018). [CrossRef]  

30. Z. Lin, V. Liu, R. Pestourie, and S. G. Johnson, “Topology optimization of freeform large-area metasurfaces,” Opt. Express 27(11), 15765–15775 (2019). [CrossRef]  

31. J. A. Fan, “High performance metasurfaces based on inverse design (conference presentation),” (2017).

32. A. Zhan, R. Gibson, J. Whitehead, E. Smith, J. R. Hendrickson, and A. Majumdar, “Controlling three-dimensional optical fields via inverse mie scattering,” Sci. Adv. 5(10), eaax4769 (2019). [CrossRef]  

33. D. W. Mackowski and M. I. Mishchenko, “Calculation of the t matrix and the scattering matrix for ensembles of spheres,” J. Opt. Soc. Am. A 13(11), 2266–2278 (1996). [CrossRef]  

34. Y. lin Xu, “Electromagnetic scattering by an aggregate of spheres,” Appl. Opt. 34(21), 4573–4588 (1995). [CrossRef]  

35. A. Egel, L. Pattelli, G. Mazzamuto, D. S. Wiersma, and U. Lemmer, “Celes: Cuda-accelerated simulation of electromagnetic scattering by large ensembles of spheres,” J. Quant. Spectrosc. Radiat. Transfer 199, 103–110 (2017). [CrossRef]  

36. C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch, “Adjoint shape optimization applied to electromagnetic design,” Opt. Express 21(18), 21693–21701 (2013). [CrossRef]  

37. J. B. Schneider and I. C. Peden, “Differential cross section of a dielectric ellipsoid by the T-matrix extended boundary condition method,” IEEE Trans. Antennas Propag. 36(9), 1317–1321 (1988). [CrossRef]  

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

39. S. C. Ligon, R. Liska, J. Stampfl, M. Gurr, and R. Mülhaupt, “Polymers for 3d printing and customized additive manufacturing,” Chem. Rev. 117(15), 10212–10290 (2017). [CrossRef]  

40. D. Theobald, A. Egel, G. Gomard, and U. Lemmer, “Plane-wave coupling formalism for T-matrix simulations of light scattering by nonspherical particles,” Phys. Rev. A 96(3), 033822 (2017). [CrossRef]  

41. P. C. Waterman, “Symmetry, unitarity, and geometry in electromagnetic scattering,” Phys. Rev. D 3(4), 825–839 (1971). [CrossRef]  

42. M. I. Mishchenko, L. D. Davis, and A. A. Lacis, Scattering, Absoprtion, and Emission of Light by Small Particles (NASA Goddard Institute for Space Studies, 2002).

43. A. Doicu, T. Wriedt, and J. A. Eremin, Light Scattering by Systems of Particles: Null-Field Method with Discrete Sources (Springer, 2006).