## Abstract

The development of metasurfaces has enabled unprecedented portability and functionality in flat optical devices. Spaceplates have recently been introduced as a complementary element to reduce the space between individual metalenses, which will further miniaturize entire imaging devices. However, spaceplates necessitate an optical response which depends on the transverse spatial frequency component of a light field — therefore making it challenging both to design them and to assess their ultimate performance and potential. Here, we employ inverse-design techniques to explore the behaviour of general thin-film-based spaceplates. We observe a tradeoff between the compression factor $\mathcal{R}$ and the numerical aperture NA of such devices; we obtained a compression factor of $\mathcal{R} = 5.5$ for devices with an NA = 0.42, and up to a record $\mathcal{R} = 340$ with NA of 0.017. Our work illustrates that even simple designs consisting of realistic materials (*i.e.*, silicon and glass) permit capable spaceplates for monochromatic applications.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

## 1. Introduction

Nanostructured surfaces known as metasurfaces have been gaining much attention for enabling the creation of flat optics [1–3]. Among the most notable devices are metalenses, optical elements that promise to yield ultra-thin imaging systems by replacing the relatively thick glass of lenses with a much thinner optic [4–6]. Unfortunately, these recent advances in miniaturizing lenses do not address the space *between* the numerous optics, sources, and sensors, which make up the largest part of most imaging systems.

Recently, a new optical element has been proposed that can simulate free-space propagation to further reduce imaging systems, called a ‘spaceplate’ [7–9]. This device consists of a plate of thickness $d$; however a transmitted lightfield will experience an effective propagation of length $d_\mathrm {eff}>d$ (Fig. 1). This expanded propagation is accomplished through a non-local operation [10–16] that imparts an angle-dependent phase onto an incident beam. By replacing free-space, spaceplates could one day lead to ultra-thin telescopes, microscopes, and cameras, for example eliminating the camera bump on the back of smartphones.

Past works have experimentally demonstrated spaceplates in the form of a low-index medium or a birefringent crystal [7], and have theoretically demonstrated spaceplates in the form of thin-film multilayer stacks [7], a two-dimensional photonic crystal slab [8], and coupled cavities [9]. These devices featured compression factors $\mathcal{R}\equiv (d_\mathrm {eff}/d)$ that ranged from $\mathcal{R}=1.1$ to $\mathcal{R}=144$, with other metrics of varying performance, such as numerical aperture (NA), insertion loss, spectral operating bandwidth, and total simulated space $d_\mathrm {eff}$. Aside from the special case of resonator-based spaceplates [9], a general design methodology and identification of the fundamental limitations of high-performance spaceplates has remained a challenge. This stands in contrast to other more mature devices, such as metalenses [17] or invisibility cloaks [18], for which there are established design prescriptions. As a result, it is as of yet unknown what the fundamental tradeoffs or limitations on these performance metrics may be for this new type of device.

In this paper, we expand on previous inverse-design work on a class of spaceplates based on multilayer thin-film stacks. The main goal is to explore the design parameter space for this class of spaceplates. In doing so, we find high-performance structures and observe trade-offs between their various performance metrics and design constraints. Thin-film stacks are a mature technology that can consist of structures of over a thousand layers. It is the basis of standard commercial dielectric mirrors that exhibit a NA of 0.35 and reflectivity >99% across the visible spectrum. Commercial thin-film stack bandpass-filters are available that have bandwidths of 350 nm and a transmission of >90%. Multilayer stacks provide a complex and large parameter space conducive to optimization.

Whereas previous work featuring this multilayer platform used a genetic algorithm [7] and produced a maximum $\mathcal{R}=4.9$, here we develop a gradient descent optimization that produces much higher performance spaceplates. We switched to gradient descent due to its ability to find local maxima quickly [19]. We kept the global search aspect of the genetic algorithms by concurrently optimizing a set of random starting structures. In this paper, we do two types of optimizations. The first is aimed at exploring the space of realistic devices made solely from silica and silicon. It optimizes the layer thickness. The second type of optimization aims to reveal the role of refractive index in spaceplate performance. We optimize a continuously variable index while keeping the layer thicknesses fixed. However, the focus of this paper is not the optimization algorithm but rather the exploration of the parameter space for multilayer spaceplates.

In Section 2, we introduce the theory behind the spaceplate and give an overview of the structure optimization method. In Section 3, we present the results of parameter space exploration. In Section 4, we discuss the implications of these results including design tradeoffs and then conclude.

## 2. Theory

A spaceplate will replace a slab of a medium of length $d_\mathrm {eff}$ if the spaceplate acts on every possible incident planewave in the same manner as the slab. In particular, each planewave incident at $z=0$ should be transmitted through the spaceplate at $z=d$ without loss and receive an angle-dependent phase of the form

where $\theta$ is the incident angle of the light relative to the interface normal, and $\mathbf {k}=(k_x,k_y,k_z)$ is the wavevector in the background medium [7]. Throughout this publication, the background medium is air ($n=1$).The spaceplates that we consider are made from layers of materials with parallel interfaces along the $x$-$y$ plane. As such, they will conserve momentum along this plane. Moreover, the imparted phase will be independent of position. These equivalent three aspects describe that a ray exiting a spaceplate will not change its propagation angle, as illustrated in Fig. 1. This angle-conservation is contrary to the behaviour of a ray at a lens, where the ray angle changes due to refraction. In this way, lenses and spaceplates are complementary devices.

A spaceplate is a non-local optic in the sense used in the theory of metamaterials [10,11]. A general linear optic is associated with the function $G(x,x')$ and transforms an input field as $E_\mathrm {out}(x)=\int E_\mathrm {in}(x')G(x,x')dx'$. If the output field $E_\mathrm {out}$ at $x$ only depends on the input field local to $x$, *i.e.*, $E_\mathrm {in}(x)$, then $G(x,x')=G(x)$ and the optic is classified as ‘local’. An example is a simple phase-mask or thin lens, $G(x)=\mathrm {exp}[-i\beta (x)]$. All optics for which $G(x,x')\neq G(x)$ are ‘non-local’. In particular, the effect of the spaceplate is to transform the input field as $E_{\mathrm {out}}(x)=\frac {1}{\sqrt {2\pi }}\int E_\mathrm {in}(k_{x})H(k_{x})e^{-ik_{x}x}dk_{x}$, where $H(k_x)=\mathrm {exp}(-i \phi )$ is the Fourier transfer function of an ideal spaceplate. By the convolution theorem, $E_{\mathrm {out}}(x)=\int E_\mathrm {in}(x')h(x-x')dx'$, where $h(x)$ is the inverse Fourier Transform of $H(k_x)$. Through $h(x-x')$ the output field at $x$ will depend on the input field at all positions, *i.e.*, non-locally. The spaceplate is an extreme example of a non-local optic since its action solely depends on the input angle $\theta$.

In order to quantify the performance of the spaceplates that we generate through inverse design, we use the Strehl ratio $S$ [20]. This quantity is commonly used to measure the quality of lenses. It is the ratio of the peak intensity of a focus of a Gaussian beam that travels through a real (*i.e.*, aberrated) lens to the peak intensity created by an ideal lens. $S$ quantifies aberration performance from poor to ideal by ranging 0 to 1, respectively. Most importantly for this work, a Strehl ratio of $S \geq 0.8$ is considered ‘diffraction limited’ [21]. That is, with an image system containing a lens with performance $S \geq 0.8$, diffraction will limit image resolution more than aberrations will. Similarly, a real *(i.e.*, designed) spaceplate will impart phase errors to a beam that is already focusing through it and so the Strehl ratio is also a relevant metric. When calculating a high Strehl ratio, the following approximation [22] holds:

Here, $\sigma$ is the root-mean-square-error (RMSE) between the phase of an ideal spaceplate $\phi$ and a designed spaceplate $\phi '$ for the same nominal $d_\mathrm {eff}$,

where $\left \langle f\right \rangle \equiv \left ({\textstyle \sum _{i=1}^{N}}f(\theta _{i})\right )/N$ denotes an average over angle across the NA of the device. Here, $N$ is the number of points between $\theta =0$ and the acceptance half-angle $\theta _{\mathrm {NA}}=\sin ^{-1}{(\mathrm {NA})}$. Thus, the RMSE and, in turn, the Strehl ratio quantify the aberrations of a designed spaceplate within its numerical aperture.In our exploration, we focus on multilayer stacks. These can be easily and accurately modeled using the transfer-matrix method (TMM) [23]. This method is used to obtain the phase $\phi$ and the transmittance $T$ of a stack over NA. We consider $p$-polarized light at an operating wavelength of $\lambda =1550$ nm, thus restricting our study to monochromatic operation. Monochromatic spaceplates could be used in the many applications that utilize narrow bands of frequencies. The requirements of commercial applications are diverse and often proprietary. While it is beyond the scope of this paper to target a specific application, monochromatic spaceplates could be used to reduce the size of sensors (*e.g.*, LIDAR or retina scanners), monochromatic displays, or spatial mode mux/demux devices in telecommunications.

We use gradient descent as our optimization technique, where we define the figure of merit (FOM) that we maximize to be $1/\sigma (\phi ', \phi _\mathrm {target})$. The RMSE $\sigma$ is taken between the phase of a candidate spaceplate and an ideal spaceplate with $d_\mathrm {eff}=\mathcal{R}_{\mathrm {target}}d$. Each device optimization thus targets a $\mathcal{R}_{\mathrm {target}}$ and NA. The result of the optimization is the designed device. Its aberration performance is quantified by the Strehl Ratio (Eq. (2)) calculated using RMSE $1/\sigma (\phi ', \phi _\mathrm {fit})$, where $\phi _\mathrm {fit}$ is the phase of an ideal spaceplate where $\mathcal{R}_{\mathrm {fit}}$ is found from a fit of Eq. (1) to the device phase $\phi '$ over the target NA.

It is unknown whether this FOM establishes a convex optimization for multilayer structures. If it is not convex, the globally optimal solution may not be reached. In order to temper this potential limitation, we use a randomization procedure to globally search the parameter space. We describe this procedure in more detail in the following section.

## 3. Optimizations and results

We explore the optimization of two different parameters. The first explores the space of realistic devices by using solely silica and silicon and optimizes the layer thicknesses. The goal of this optimization is to investigate the role of layer thickness on spaceplate performance and also to identify realistic structures with high $\mathcal{R}$ and/or NA. The second optimization type uses fixed layer thicknesses and optimizes each layer’s refractive index. The goal is to observe the influence of index contrast on the spaceplate performance.

Both optimization types use gradient descent. For both types, we pre-define the number of layers, the targeted compression factor $\mathcal{R}$, and acceptance half-angle $\theta _\mathrm {NA}$ (*i.e.*, the $\mathrm {NA}=\sin {(\theta _\mathrm {NA})}$). It is unknown whether either optimization type is convex. To mitigate this possibility we introduce randomness, with each run starting from a structure comprising random initial parameters. In this sense, we are performing a global search to increase the likelihood of finding global or near-global optima. For the first optimization type, we randomize the layer thicknesses by simultaneously optimizing a ‘swarm’ of 200 gradient descent optimization runs. In the second type, we randomize the refractive index profile; however, due to the increased number of degrees of freedom, the optimization is much more numerically cumbersome, and so we optimize a swarm of only 10 devices.

#### 3.1 Optimizing layer thicknesses

The first approach considers a structure of alternating layers of silicon ($n=3.48$) and silica ($n=1.44$). The structure is surrounded on both sides by air ($n=1$). Here, we set a fixed number of thin-film layers and the thicknesses of the individual layers are used as the optimization parameters. Our optimization starts by initializing a group of 200 spaceplate candidate designs, each with randomly generated layer thicknesses. For each optimization case, the code then applies the gradient descent algorithm to each layer of the spaceplate, discretely differentiating by the layer thickness to find the gradient. For a given candidate design, the algorithm starts by simultaneously optimizing the first three layers, and once they are fully optimized, it includes the next layer into the simultaneous optimization process. The $\mathrm {FOM}=1/\sigma (\phi ', \phi _\mathrm {target})$. Once all of the layers have been sequentially included and optimized, the spaceplate’s actual compression $\mathcal{R}_{\mathrm {fit}}$ and Strehl ratio are assessed.

Figure 2 illustrates a resulting device consisting of 27 layers optimized for operation for an acceptance half-angle of $\theta _\mathrm {NA}=1^{\circ }$ ($\mathrm {NA}=0.017$). The device is ${7.6}\;\mathrm{\mu} \textrm {m}$ thick and replaces propagation for 2.6 mm, corresponding to a compression factor of $\mathcal{R} = 340$. This value is the largest observed in any study. This device has a maximum insertion loss of only 3.3 dB over its designated angular operating bandwidth. Figure 2(b) gives the imparted phase of the device (blue line), which fits well to the ideal spaceplate curve (dashed line, given by Eq. (1)). This device has a Strehl ratio of $S=0.9992$, far above diffraction-limited performance.

To further probe the potential performance of the thin-film multilayer stack platform as a spaceplate, we explored the parameter space in a thorough series of simulations. Though a direct search of this type is not in any sense exhaustive (*i.e.*, completely comprehensive), its results can be instructive in drawing certain conclusions, as we shall see. In our study, for a given number of layers ranging from 9 – 29, we swept over targeted compression factors from $\mathcal{R}_\mathrm {target}=2$ – 60 and explored devices with five NAs: 0.087, 0.17, 0.26, 0.34, and 0.42, corresponding to acceptance half-angles of $\theta _\mathrm {NA}=5^{\circ }$, $10^{\circ }$, $15^{\circ }$, $20^{\circ }$ and $25^{\circ }$, respectively.

The results from this study are summarized in Fig. 3. Figure 3(a) shows the maximum $\mathcal{R}$ obtained using the optimization procedure outlined above as a function of the device’s NA independent of the number of layers. As mentioned above, the largest compression factor, $\mathcal{R}=340$, is obtained for devices with the smallest NA. Conversely, for the largest NA ($\mathrm {NA}=0.43$, corresponding to an acceptance half-angle of $\theta _\mathrm {NA}=25^{\circ }$), we obtain $\mathcal{R}=5.5$ in a ${1.61}\;\mathrm{\mu} \textrm {m}$-thick device. Thus, we observe a tradeoff between NA and $\mathcal{R}$. Moreover, we observe that for large NA, the highest $\mathcal{R}$ occurs for structures with fewer layers.

Figure 3(b) contains a summary of the performance of multilayer stacks comprising only 13 layers as a function of $\mathcal{R}_\mathrm {target}$ for a range of NAs. We observe a clear tradeoff between $S$ and $\mathcal{R}$ of a spaceplate for a given NA. Conversely, for a given $\mathcal{R}$, a higher NA results in lower device Strehl ratio. Though we do not reproduce it here, this trend holds for all optimized devices for every number of layers we explored (9, 11,…, 27, 29 layers).

We observe threshold values for $\mathcal{R}_\mathrm {target}$ above which $S$ starts to decrease from 1. At these thresholds, the optimization procedure is obtaining devices with increasing aberration (*i.e.*, decreasing $S$) as it is more challenging to reach larger values of $\mathcal{R}_\mathrm {target}$

In Fig. 3(d), we explore the impact of the total number of layers on the performance of the device. Here, we fix the NA to 0.17, and optimize devices for a range of $\mathcal{R}_\mathrm {target}$ and for a number of layers ranging from 9 to 29. Increasing the number of layers enables devices with better performance and higher compression factors, which makes intuitive sense as doing so expands the space of possible devices.

In Fig. 3(c), we show that for a spaceplate designed for a given target compression factor $\mathcal{R}$, when increasing the thickness of any one layer, the compression factor is oscillatory. The multiple maxima suggest that the optimization is non-convex and that the results of optimizations may not be global maxima.

Our search yielded a few other qualitative observations of interest. First, our optimization procedure favoured devices with an odd number of layers, and where the outermost layers had a higher refractive index. Additionally, devices optimized for operation with p-polarized light tended to yield higher compression factors.

#### 3.2 Optimizing layer indices

In our second approach to optimizing the spaceplate structures, we fix the total number of layers and their respective thicknesses and use the refractive indices of the layers as the optimization parameters. This approach is similar to the one explored in previous work featuring a different non-local operation in a thin-film multilayer stack [11]. Here, we set the thickness of 101 individual layers to 5 nm for a total device thickness of 505 nm, and we constrain the refractive index to a chosen range. These ranges define two optimization sub-types: the first is a realistic range for dielectric materials at optical frequencies ($n=1$ – $4$), and the second is an expanded range of any positive value ($n>0$). This latter range explicitly includes materials for which $n<1$, so-called epsilon-near-zero media [24,25]. Unlike the dielectric materials in the first sub-type, these tend to have high absorption, which makes them less feasible for realistic devices.

In our first attempts at this second type of optimization, we observed that the resulting structures exhibited impractically small transmission. To counter this effect, we also maximized the transmission by including it in the FOM$=1/(\sigma (1 - \langle 10\log _{10}(T) \rangle ))$. Here, $T$ is the normalized transmittance, and therefore $\langle 10\log _{10}(T) \rangle$ is the magnitude of the insertion loss in dB, where $\langle \rangle$ denotes an average over angle across the NA of the device, as in Eq. (3). Since this quantity is located in the denominator of the FOM, the optimization works to reduce loss, increasing overall transmittance. In this approach, the total number of layers is substantially larger than in the devices explored previously. This optimization approach is much more numerically cumbersome and so we are forced to iterate over fewer parameters. Consequently we only optimize for a swarm of 10 candidate devices at a time.

For $\mathrm {NA}=0.26$, the first sub-type ($n=1$ – $4$) yielded a maximum compression factor of $\mathcal{R}=6.4$ (Fig. 4(a)) when targetting $\mathcal{R}_\mathrm {target}=7$. The second sub-type ($n>0$) yielded a device with refractive indices ranging from $n=0.2$ to $n=14$, and a maximum compression factor of $\mathcal{R}=80$ for a target value $\mathcal{R}_\mathrm {target}=80$ (Fig. 4).

Though this second type of optimization yields less feasible devices due to the continuous refractive index values, it suggests that the ultimate compression of multilayer spaceplates is likely tied to the refractive index of the constituent materials. To test this aspect more thoroughly, we conducted a new optimization study of the first type (*i.e.*, optimizing the layer thicknesses only) where we fix the number of layers to 17, and explored the indices of the two materials used. We found that as the index contrast (the ratio of the refractive indices, $n_H/n_L$) increases, the achieved compression factor increases as well (Fig. 4(e)). Additionally, spaceplates with similar index contrast perform similarly even though the materials used may widely differ (*e.g.*, the orange and light green curves in Fig. 4(e), corresponding to $n_H/n_L = 4/2$ and 2/1, respectively).

## 4. Discussion and conclusion

Our study with the first optimization type suggests that the maximum compression factor is limited by the NA of the device. For example, for devices with an identical $S$ and number of layers, reducing the NA from 0.26 to 0.17 had the effect of nearly doubling $\mathcal{R}$ from 9.6 to 18.0. This observation echoes a recent finding that the total effective propagation distance $d_\mathrm {eff}$ of resonator-based spaceplates is limited by the NA [9]. Despite this limit, we find that spaceplates based on this platform are shown to enable diffraction-limited operation with record high compression factors for large numerical apertures.

We had expected the second optimization type to yield devices with larger compression factors since it considered a larger parameter space with access to many more degrees of freedom. However, the best compression factor obtained by this approach ($\mathcal{R}=6.4$) fell short of the solution found in the previous search for $\mathrm {NA}=0.26$, which yielded a maximum compression factor of $\mathcal{R}=11.5$. We suspect that the dramatic increase in optimization parameters hampered locating the global maximum for this class of devices. Nevertheless, increasing the range of allowable indices further increased the compression factor from $\mathcal{R}=6.4$ to 80. Paired with the study where we explored making spaceplates out of two different indices using the first optimization type, we can comfortably claim that $\mathcal{R}$ is also limited by index contrast.

In applications, it may be desirable to constrain $d_\mathrm {eff}$ to a target length in the optimization. Our designs result in a $d_\mathrm {eff}$ ranging from about ${10}\;\mathrm{\mu} \textrm {m}$ to 1 mm, a factor of 100, which demonstrates an impressive range of performance. In the second type of optimization, we have constrained the physical size of the spaceplate $d=505$ nm while targeting a series of values of $\mathcal{R}$, which is equivalent to targeting specific values of $d_\mathrm {eff}$. We demonstrated a range of performance with this $d$. This shows that it is possible to target both $d_\mathrm {eff}$ and $\mathcal{R}$. Moreover, recent work has proposed a periodic design that allows for any chosen $d_\mathrm {eff}$ [9]. Future work will focus on the possibility of cycling our optimization procedures in order to achieve a desired $d_\mathrm {eff}$.

In summary, we have successfully inversely designed high-performing spaceplates based on multilayer thin-film stacks. We observed tradeoffs between the number of layers, NA, compression factor $\mathcal{R}$, and the degree of aberration (*i.e.*, the Strehl ratio). While this work was focused on monochromatic spaceplates, these tradeoffs suggest there may be similar tradeoffs with the spaceplate’s spectral bandwidth. Broadband operation and achromatic design are topics for future research. We observe signs that our gradient descent is non-convex which suggests that other inverse-design schemes might be even more fruitful than gradient descent and, thus, may fully capitalize on the enormous design parameter space of the thin-film platform. Finally, while previous spaceplate designs have achieved very high $\mathcal{R}$ [8], this work demonstrates that even higher $\mathcal{R}$ can be achieved in thin-film stacks, a mature commercial fabrication platform.

## Funding

Canada First Research Excellence Fund (Transformative Quantum Technologies); Canada Research Chairs (Canada Research Chairs Program).

## Acknowledgments

The authors thank Marc-André Renaud for programming advice, and Víctor J. López-Pastor and Ali H. Alhulaymi for helpful discussions.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request. The gradient descent code used to design spaceplates can be found online [26].

## 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. **N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nat. Mater. **13**(2), 139–150 (2014). [CrossRef]

**3. **W. T. Chen, A. Y. Zhu, and F. Capasso, “Flat optics with dispersion-engineered metasurfaces,” Nat. Rev. Mater. **5**(8), 604–620 (2020). [CrossRef]

**4. **M. Khorasaninejad, W. T. Chen, R. C. Devlin, J. Oh, A. Y. Zhu, and F. Capasso, “Metalenses at visible wavelengths: Diffraction-limited focusing and subwavelength resolution imaging,” Science **352**(6290), 1190–1194 (2016). [CrossRef]

**5. **M. Khorasaninejad and F. Capasso, “Metalenses: Versatile multifunctional photonic components,” Science **358**(6367), eaam8100 (2017). [CrossRef]

**6. **H. Liang, Q. Lin, X. Xie, Q. Sun, Y. Wang, L. Zhou, L. Liu, X. Yu, J. Zhou, T. F. Krauss, and J. Li, “Ultrahigh Numerical Aperture Metalens at Visible Wavelengths,” Nano Lett. **18**(7), 4460–4466 (2018). [CrossRef]

**7. **O. Reshef, M. P. DelMastro, K. K. M. Bearne, A. H. Alhulaymi, L. Giner, R. W. Boyd, and J. S. Lundeen, “An optic to replace space and its application towards ultra-thin imaging systems,” Nat. Commun. **12**(1), 3512 (2021). [CrossRef]

**8. **C. Guo, H. Wang, and S. Fan, “Squeeze free space with nonlocal flat optics,” Optica **7**(9), 1133–1138 (2020). [CrossRef]

**9. **A. Chen and F. Monticone, “Dielectric Nonlocal Metasurfaces for Fully Solid-State Ultrathin Optical Systems,” ACS Photonics **8**(5), 1439–1447 (2021). [CrossRef]

**10. **G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Nonlocal transformation optics,” Phys. Rev. Lett. **108**(6), 063902 (2012). [CrossRef]

**11. **A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alu, and N. Engheta, “Performing Mathematical Operations with Metamaterials,” Science **343**(6167), 160–163 (2014). [CrossRef]

**12. **H. Kwon, D. Sounas, A. Cordaro, A. Polman, and A. Alù, “Nonlocal Metasurfaces for Optical Signal Processing,” Phys. Rev. Lett. **121**(17), 173004 (2018). [CrossRef]

**13. **A. Cordaro, H. Kwon, D. Sounas, A. F. Koenderink, A. Alù, and A. Polman, “High-index dielectric metasurfaces performing mathematical operations,” Nano Lett. **19**(12), 8418–8423 (2019). [CrossRef]

**14. **M. Gerken and D. A. B. Miller, “Multilayer thin-film structures with high spatial dispersion,” Appl. Opt. **42**(7), 1330–1345 (2003). [CrossRef]

**15. **F. Monticone, C. A. Valagiannopoulos, and A. Alù, “Parity-time Symmetric nonlocal metasurfaces: All-angle negative refraction and volumetric imaging,” Phys. Rev. X **6**(4), 041018 (2016). [CrossRef]

**16. **A. C. Overvig, S. C. Malek, and N. Yu, “Multifunctional Nonlocal Metasurfaces,” Phys. Rev. Lett. **125**(1), 017402 (2020). [CrossRef]

**17. **F. Presutti and F. Monticone, “Focusing on bandwidth: achromatic metalens limits,” Optica **7**(6), 624–631 (2020). [CrossRef]

**18. **R. Fleury, F. Monticone, and A. Alù, “Invisibility and Cloaking: Origins, Present, and Future Perspectives,” Phys. Rev. Appl. **4**(3), 037001 (2015). [CrossRef]

**19. **S. Boyd and L. Vandenberghe, * Convex Optimization* (Cambridge university press, 2009).

**20. **A. V. D. Bos, “Aberration and the Strehl ratio,” J. Opt. Soc. Am. A **17**(2), 356–358 (2000). [CrossRef]

**21. **A. Maréchal, “Étude des effets combinés de la diffraction et des aberrations géométriques sur l’image d’un point lumineux,” Éditions de la Revue d’optique théorique et instrumentale (1948).

**22. **V. N. Mahajan, “Strehl ratio for primary aberrations in terms of their aberration variance,” J. Opt. Soc. Am. **73**(6), 860–861 (1983). [CrossRef]

**23. **A. Yariv and P. Yeh, * Photonics: Optical electronics in modern communications* (Oxford University Press, New York, NY, 2007), 6th ed.

**24. **A. Alù, M. G. Silveirinha, A. Salandrino, and N. Engheta, “Epsilon-near-zero metamaterials and electromagnetic sources: Tailoring the radiation phase pattern,” Phys. Rev. B **75**(15), 155410 (2007). [CrossRef]

**25. **O. Reshef, I. D. Leon, M. Z. Alam, and R. W. Boyd, “Nonlinear optical effects in epsilon-near-zero media,” Nat. Rev. Mater. **4**(8), 535–551 (2019). [CrossRef]

**26. **MagicIncubator, “Inverse design spaceplate,” Github (2021), https://github.com/MagicIncubator/Inverse_Design_Spaceplate.