Optical spatial differentiation, which overcomes the speed and energy limitations of digital image analysis, plays a pivotal role in all-optical information processing, such as image reconstruction and analog optical computing [1–3]. In recent years, optical spatial differentiators with different structures have been proposed both theoretically [4,5] and experimentally [6–12], such as multilayer slabs [13], subwavelength gratings [14,15], photonic crystals [16,17], and metasurfaces [18–23]. Embeddable optical spatial differentiation metasurfaces [24–26] have been utilized to detect image edges induced by both amplitude and phase variation [27,28], offering a compact footprint and high flexibility due to their abilities to manipulate electromagnetic waves at subwavelength scale without Fourier transform lenses [29], such as quasi-bound states metasurface [30] and laplace differential metasurface [31]. By manipulating differentiation-gains [32], a high-contrast differentiation metasurface effectively extracts edges associated with weak intensity variations, such as those found in medical images. In contrast, its high-resolution [33,34] counterpart is utilized to achieve high-pass filtering of the object image and detect fine structural details, for instance, the precise segmentation of tissues and cells. A trade-off exists between the high-contrast and resolution functionalities of embeddable optical spatial differentiation metasurfaces, making it challenging to be optimized by forward design strategies. The existing inverse design methods, such as a heuristic algorithm [35,36], topology optimization [37–40], and neural networks [41,42], mainly focus on the spatial domain to obtain optimized performance of nanophotonic devices, but the structural design of spatial differentiation metasurfaces with spatial-frequency domain functionalities is rarely studied.

In this paper, two kinds of embeddable optical spatial differentiation metasurfaces with functionalities of high root-mean-square (RMS) contrast and high resolution were inversely designed by the spatial-frequency Trust-Region algorithm, which experimentally extracts boundaries induced by weak intensity variations and fine structures of unstained cell images without a Fourier transform lens. The spatial-frequency Trust-Region algorithm framework was built by combining the Fourier modal method (FMM) and Trust-Region algorithm to iterate the nanostructure of the metasurface until its spatial-frequency transfer function (SFTF) converges to customized differentiation-gain. The designed differentiation metasurfaces were fabricated on a quartz substrate by plasma enhanced chemical vapor deposition, electron-beam lithography (EBL), and reactive ion etching (RIE). Without a Fourier transform lens, the experimental resolutions of high-contrast and high-resolution differentiation metasurfaces were respectively measured as 3.9 μm (NA 0.26) and 1.2 μm (NA 0.87) by a standard resolution test chart. The boundary images of unstained cells and quasi-transparent cells were successfully extracted by the fabricated metasurfaces, which have the ability of detecting both amplitude and phase variation. Meanwhile, dynamic image edge detections of paramecia were also implemented without motion smear. By integrating with optical microscopes, the fabricated high-contrast differentiation metasurface was first utilized to rapidly locate the position of live cells or tissues owing to high transmission efficiency; then the high-resolution counterpart was employed to detect the fine-structure information due to high-spatial-frequency components of SFTF. Although the spatial-frequency Trust-Region algorithm cannot guarantee to give a globally optimal solution, the structure-optimized metasurfaces still implemented the customized functionality optical spatial differentiation, which will potentially apply in microsurgery and magnetic resonance imaging.

Embeddable metasurface geometric structures with desired differentiation-gains were optimized by the spatial-frequency Trust-Region algorithm, which is realized by integrating the Fourier modal method with the Trust-Region algorithm. The Fourier modal method [43,44] is employed for the purpose of attaining a spatial-frequency response of periodic optical metasurfaces due to its simplicity and versatility. By solving the algebraic Maxwell eigenvalue equation in the spatial-frequency domain, FMM is possible to mathematically construct a complete set of eigenmodes for metasurfaces spanning all possible optical fields within the finite dimensional numerical framework.

Based on FMM, both the permittivity and electromagnetic field distribution of a single unit in a periodic metasurface are uniformly discretized (ϵp′q′, Ep,q), and then converted into the spatial-frequency domain by discrete Fourier transform: {ϵp′,q′=∑m′=04M∑n′=04Nϵm′,n′exp(iπp′2Mm′+iπq′2Nn′),Ep,q=∑m=02M∑n=02NEm,nexp(iπpMm+iπqNn),(1)where p′∈[0,4M], q′∈[0,4N], p∈[0,2M], q∈[0,2N], and M, N∈Z. ϵm′n′ and Em,n are Fourier coefficients of ϵp′q′ and Ep,q, respectively. By substituting Eq. (1) into Maxwell’s equations and solving it in the spatial-frequency domain, all the Bloch eigenmodes Eη(kx,ky) in a metasurface structure are derived. The electric field in a metasurface is represented as a linear superposition of Bloch eigenmodes: Emeta(kx,ky)=∑ηCηEη(kx,ky),(2)where Cη is linear coupling coefficients; kx,ky are transverse incident wavevectors. Based on the plane-wave condition and transverse field continuation conditions at the interface of air and a metasurface, the transmitted electric field Etransmit(kx,ky) is calculated; then the spatial frequency transfer function (SFTF) of the metasurface is derived: $$\begin{gathered} h(kx,ky)=Etransmit(kx,ky)Eincident(kx,ky) \\ =(Extransmit(kx,ky),Eytransmit(kx,ky),Eztransmit(kx,ky))(Exincident(kx,ky),Eyincident(kx,ky),Ezincident(kx,ky)), \end{gathered}$$(3)where Eincident(kx,ky) and Etransmit(kx,ky) are incident and transmitted electric fields, respectively. The permittivity and electromagnetic field distribution are represented by discrete Fourier transform instead of the pseudo-Fourier series used in traditional FMM, which avoids the Gibbs phenomenon to improve the convergence of FMM [45,46].

The spatial-frequency Trust-Region algorithm is then built by merging FMM with the Trust-Region algorithm to design optical spatial differentiation metasurfaces with the desired Fourier spatial-frequency response. The framework of the spatial-frequency Trust-Region algorithm for a spatial differentiation metasurface is shown as the following five steps and diagrammed in Fig. 1.

Step 1: Set target SFTF of the metasurface with desired differentiation-gains [47]: h(kx,ky)=ak(kx2+ky2),(4)where ak is the differentiation-gain of the metasurface; kx,ky are transverse incident wavevectors. A large differentiation-gain ak in Eq. (4) indicates high transmission efficiency of metasurfaces leading to differentiation operation characterized by a high-contrast functionality, while small differentiation-gain ak allows more high-spatial-frequency components in SFTF, achieving differentiation featured by a high-resolution functionality to finer image structures.

Step 2: Starting with an initial metasurface structure, derive its SFTF by Eqs. (1)–(3).

Step 3: Calculate the Pearson correlation coefficient O(s) between the iterating SFTF and target counterpart: O(s)=coff[hiterating(kx,ky),hdesired(kx,ky)],(5)where s is independent variables representing the structural parameters of the metasurface to be optimized. Then construct an optimization subproblem: {min1−O(si+Δsi)s.t.si+Δsi∈Ω,(6)where metasurface structural parameter increments Δsi and O(si+Δsi) are the optimization variables and objectives, respectively; subscript i stands for the i-th iteration.

Step 4: The optimization subproblem of Eq. (6) is a multivariate, nonlinear optimization problem of metasurface inverse design; therefore the Trust-Region algorithm is applied due to its fast global convergence: the optimal shift vector of independent variables Δsimin is solved by the following subproblem to ensure the largest increment of objective function O(si+Δsimin) in the range of the trust region: {min1−Ui(Δsi)=1−(fi+giTΔsi+12ΔsiTGiΔsi)s.t.‖Δsi‖≤hi,(7)where Ui(Δsi) is a second-order Taylor expansion of objective function O(si+Δsi) at point si; hi is the trust region radius; fi, giT, and Gi are zeroth-, first- and second-order Taylor expansion coefficients of O(si+Δsi), respectively. Then the ratio ri is calculated by the solved optimal shift vector Δsimin in Eq. (7): ri=fi(si)−fi(si+Δsimin)Ui(0)−Ui(Δsimin),(8)based on which the independent variables and trust region radius of the next iteration are adjusted.

Step 5: Update the metasurface structural parameters s, and then repeat Steps 2–4 until the Pearson correlation coefficient O(s) converges to one.

The Trust-Region algorithm is fundamentally a local optimization algorithm. In order to approach the global optimum, the independent optimization processes were conducted starting from different initial parameter sets. When the optimized structures converge to a similar configuration, the solution is considered to be approaching the global optimum.

The material of an embeddable differentiation metasurface is amorphous silicon (α-Si) with a refractive index of nα-Si=3.75 at 808 nm. The α-Si metasurface is placed on a quartz substrate (refractive index is nquartz=1.453). The laser is incident normally with 45° polarization. The initial metasurface structure is set as periodic α-Si polygon voxels of which the side lengths are optimization variables. Limited by actual fabrication capability, the minimum polygon side length is set to 20 nm. Through multiple iterations, the schematics of the optimized high-contrast and high-resolution metasurface are shown in Figs. 2(a) and 2(b), respectively, and structural parameters are marked.

The customized functionality of the designed embeddable differentiation metasurfaces by the spatial-frequency Trust-Region algorithm was simulated by the finite-difference time-domain (FDTD) method. By simulating the transmission of plane-waves at different transverse incident wavevectors kx,ky, the SFTF curves of two kinds of metasurfaces were obtained and are shown by the solid line in Figs. 2(e) and 2(f), respectively, of which the quadratic trend indicates the functionality of second-order optical spatial differentiation. By performing quadratic curve fitting for the SFTFs, the simulated differentiation-gains of high-gain and high-resolution differentiation metasurfaces were fitted with ak=12.8k0−2 and 1.1k0−2 in Eq. (4), and the corresponding numerical apertures (NAs) are 0.28 and 0.94, respectively, realizing customized functionalities of the metasurfaces.

The designed two-dimensional second-order optical spatial differentiation metasurfaces were fabricated experimentally. The α-Si film was deposited on a quartz substrate by plasma enhanced chemical vapor deposition (PECVD). Then metasurface structures were fabricated by a single electron-beam lithography (EBL) and reactive ion etching (RIE) step. The scanning electron microscope (SEM) images of fabricated high-contrast and high-resolution differentiation metasurfaces are shown in Figs. 2(c) and 2(d), respectively. The fabricated structural parameters of two kinds of metasurfaces are also shown in Figs. 2(c) and 2(d).

The SFTF curves of the fabricated embeddable differentiation metasurfaces were measured by the experimental set-up shown in Fig. 2(g). The metasurfaces were positioned on a high-precision rotary table to finely adjust the incident angle of a linearly polarized laser (808 nm). Transmitted laser powers of metasurfaces at different incident angles, namely, different transverse incident wavevectors kx,ky, were recorded by an optical power meter. Based on these measurements, the experimental data points of SFTF for the fabricated differentiation metasurfaces were obtained and shown in Figs. 2(e) and 2(f), while the dashed line in each figure shows a second-order polynomial fit to these experimental data points. Both experimental data points act as quadratic dependence given by Eq. (4), illustrating the functionality of second-order optical spatial differentiation of metasurfaces. By conducting quadratic curve fitting for the SFTFs, the experimental differentiation-gains of the high-contrast and high-resolution differentiation metasurfaces were determined as ak=14.8k0−2 and 1.3k0−2, and the corresponding numerical apertures (NAs) are 0.26 and 0.87, respectively, allowing for direct discrimination of the image edges with customized differentiation-gains in different situations. The discrepancies between experiment and simulation of differentiation-gains and numerical apertures are mainly caused by the fabrication deviations of metasurfaces.

The edge detection of the fabricated embeddable differentiation metasurfaces was experimentally studied by the experimental set-up shown in Fig. 3(a). A standard resolution test chart was illuminated by a linearly polarized laser (808 nm), forming the images to serve as the input of fabricated metasurfaces [Figs. 3(b)–3(d)]. The output images of the fabricated differentiation metasurfaces were captured by a charge-coupled device (CCD) camera. The captured transmitted images of the high-contrast and high-resolution differentiation metasurfaces are shown in Figs. 3(e)–3(g) and Figs. 3(h)–3(j), respectively. The normalized field profiles of incident and transmitted images at the position of the white dashed line in Figs. 3(b), 3(e), and 3(h) are extracted and shown in Figs. 3(k)–3(m), respectively. The two-peak positions of the experimental transmitted intensity profiles exactly match the edges of the incident intensity profile, namely, two kinds of differentiation metasurfaces both clearly exhibiting all edges of the input images with higher transmission or finer structure. The minimum test chart sizes that image edges can be clearly distinguished, namely, experimental resolutions are 3.9 and 1.2 μm for high-contrast and high-resolution differentiation metasurfaces, respectively, which are sufficient to detect boundaries of not only tissues but also cells.

Edge detections of cells were studied by the experimental set-up shown in Fig. 3(a). First, two types of unstained cells, onion epidermis [Fig. 4(a)] and tomato cells [Fig. 4(b)], were used as the imaging specimens. The specimen images were implemented by high-contrast and high-resolution differentiation metasurfaces, and the transmitted images are displayed in Figs. 4(e) and 4(f) and Figs. 4(i) and 4(j), respectively. The high-contrast differentiation metasurfaces extracted cell edges with high brightness, while the high-resolution counterparts detected more high-spatial-frequency information, such as outlines of the nucleus. The root-mean-square (RMS) contrast of boundary images in Figs. 4(e), 4(f), 4(i), and 4(j) is 45.5, 27.9, 30.3, and 20.8, respectively, revealing an obvious improvement by the high-contrast metasurface compared to the high-resolution one.

The transparent cells show almost uniform intensity distribution over the entire field-of-view of the microscope, which is insufficient to be directly distinguished. The optical spatial differentiation metasurfaces were utilized to detect not only intensity changes but also phase-variation usually existing at the cell boundary, which enables to visualize the transparent cells. Optical spatial differentiation of quasi-transparent mouse hepatocytes [Fig. 4(c)] and mouse spleen cells [Fig. 4(d)] was implemented by the fabricated metasurfaces as shown in Figs. 4(g), 4(k) and Figs. 4(h), 4(l), respectively. High-contrast differentiation metasurfaces are applied for weak phase-variation detection relying on its high differentiation-gain, while the high-resolution counterparts prefer accurately detecting small-size transparent cells due to the extraction of high-spatial-frequency information.

Dynamic medical imaging, such as in microsurgery, requires clear and timely edge detection to remove irrelevant background and preserve important information. Generally, image edge detection realized by computers usually has a delay of about ∼10−2s [48], leading to the motion smear of images. The computation time of all-optical spatial differentiable metasurfaces depends on the ratio of the optical path length to the speed of light, enabling “real-time” edge detection with a delay of about ∼10−11s [48], which overcomes the motion smear in dynamic biomedical imaging. Dynamic biomedical edge detection of live paramecia images was experimentally operated by the fabricated embeddable optical spatial differentiation metasurfaces. The paramecia edge-extracted results by the fabricated metasurfaces are shown in Visualization 1 in which frames were extracted as in Figs. 5(a)–5(c). Both fabricated metasurfaces detected live paramecia edges in real-time without motion smear. In comparison, high-contrast edge detection of the fabricated metasurface realized quick dynamic screening of target cells, while the high-resolution feature of another metasurface is more suitable for real-time biomedical image segmentation such as in magnetic resonance imaging (MRI) [49,50].

A two-step edge detection framework of dynamic biomedical images was proposed to maximize the advantages of both fabricated embeddable metasurfaces. As shown in Fig. 5(d), the high-contrast and high-resolution differentiation metasurfaces were embedded into 20× and 40× objective lenses of a commercial optical microscope (Leica DM750M, 10× eyepiece), respectively. Step 1 of the framework: the embedded high-contrast differentiation metasurface was utilized to rapidly locate the position and state of live cells or tissues, which reduces the impact of background noise on the target by high transmission efficiency of derivative images. Step 2 of the framework: the embedded high-resolution differentiation metasurface was utilized to segment the fine boundaries of the located cells or tissues, which provides more high-spatial-frequency information of images. As an example, a live Brachionus was observed by the two-step edge detection framework, and its observation video is attached in Visualization 2, in which the frames at Steps 1 and 2 are shown in Figs. 5(e) and 5(f), respectively.

The high-performance optical spatial differentiation metasurfaces optimized by the spatial frequency inverse design method have been successfully used in transparent and dynamic biomedical image edge detection in real-time. Compared to optical differentiators designed by forward methodology, the spatial-frequency Trust-Region algorithm counterpart realizes customizable differentiation-gains, which have potential applications for image edge detection in different situations, such as high-contrast dynamic images in microsurgery and high-resolution biomedical image segmentation in MRI. In addition to differential metasurface optimization, the spatial-frequency Trust-Region algorithm can also be employed to inversely design metasurface structures with optical functionalities in the spatial-frequency domain, such as metasurfaces for spatial filtering, image reconstruction, noise suppression, and analog optical computing. Although the wavelength used in our experiments is 808 nm, which offers superior transmittance for biological tissues, the inverse design method can also be employed to design metasurfaces with alternative materials at other wavelengths.

In conclusion, two-dimensional second-order optical spatial differentiation metasurfaces with customized differentiation-gains were inversely designed by the spatial-frequency Trust-Region algorithm and fabricated through a single lithography step, which implemented biomedical image edge detection with high-contrast or high-resolution functionalities. The fabricated embeddable metasurface with a high-resolution functionality achieved a resolution of 1.2 μm and numerical apertures (NAs) of 0.87, while the high-contrast one obtained a root-mean-square (RMS) contrast higher than that of the first with a numerical aperture of 0.26. With a large differentiation-gain, a high-contrast differentiation metasurface effectively extracts edge images associated with weak intensity variations, which can be applied in microsurgery to detect the boundaries of dynamic biomedical images without motion smear. With the small differentiation-gain allowing more high-spatial-frequency components in SFTF, the high-resolution counterpart is utilized to detect fine structural details of unstained, even transparent cell images, offering a potential avenue for precise segmentation. Despite a globally optimal solution not being guaranteed, the spatial-frequency Trust-Region algorithm has optimized metasurface structures realizing optical spatial differentiation, which can potentially be employed to inversely design metasurfaces with optical functionalities in the spatial frequency domain, such as spatial filtering, image reconstruction, noise suppression, and analog optical computing.