Abstract
Spatial differentiation is important in image-processing applications such as image sharpening and edge-based segmentation. In these applications, of particular importance is the Laplacian, the simplest isotropic derivative operator in two dimensions. Spatial differentiation can be implemented electronically. However, in applications requiring real-time and high-throughput image differentiation, conventional digital computations become challenging. Optical analog computing may overcome this challenge by offering high-throughput low-energy-consumption operations using compact devices. However, previous works on spatial differentiation with nanophotonic structures are restricted to either one-dimensional differentiation or reflection mode, whereas operating in the transmission mode is important because it is directly compatible with standard image processing/recognition systems. Here, we show that the Laplacian can be implemented in the transmission mode by a photonic crystal slab device. We theoretically derive the criteria for realizing the Laplacian using the guided resonances in a photonic crystal slab. Guided by these criteria, we show that the Laplacian can be implemented using a carefully designed photonic crystal slab with a non-trivial isotropic band structure near the point. Our work points to new opportunities in optical analog computing as provided by nanophotonic structures.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. INTRODUCTION
In image processing, spatial differentiation accomplishes image sharpening and edge-based segmentation, with broad applications ranging from microscopy and medical imaging to industrial inspection and object detection [1–5]. In these applications, of particular importance are isotropic derivative operators, whose responses are rotationally invariant. The simplest and most widely used isotropic derivative operator is the two-dimensional Laplacian [6].
Spatial differentiation can certainly be carried out with conventional digital electronic computation. However, there are many big-data applications that require real-time and high-throughput image processing, for which digital computations become challenging [7,8]. Optical analog computing may overcome this challenge by offering high-throughput derivative operation with almost no energy consumption other than the inherent optical loss associated with the differential operators [9–19]. Moreover, the recent development of nanophotonic structures such as metasurfaces and plasmonic structures have offered the possibility of doing optical differentiation using compact devices. However, most of the existing works on optical spatial differentiation using nanophotonic structures are restricted to one dimension, whereas most images are two-dimensional objects [20–26]. The two-dimensional Laplacian has been implemented using holograms [27,28] and phase-shifted Bragg gratings [29]. However, these implementations either rely upon bulky optical components or operate in the reflection mode.
In this paper, we show that a two-dimensional Laplacian can be implemented in the transmission mode, which is more compatible for image-processing applications, with the use of a photonic crystal slab, one of the most widely studied nanophotonic structures [30–32]. As shown in Fig. 1, the structure consists of a photonic crystal slab separated from a uniform dielectric slab by an air gap. The key in this design is to achieve an unusual band structure that is isotropic for the two polarizations.

Fig. 1. (a) Geometry of the photonic crystal slab differentiator, which consists of a photonic crystal slab separated from a uniform dielectric slab by an air gap. For , the geometry parameters are: , , , , where is the lattice constant. The plane above and below shows the input and output field profiles, respectively. The transmitted field profile (e.g., ring) through the device is the Laplacian of the incident field profile (e.g., disk) illuminated by a normally incident unpolarized light with frequency . (b) The coordinate system. (c) The Brillouin zone of the system.
2. THEORETICAL ANALYSIS
Our objective is to realize a Laplacian for a two-dimensional input field. Specifically, for a normally incident light beam along the axis with a transverse field profile , we would like to design an optical device for which the transmitted beam has a profile , where is the Laplacian. The task of realizing an ideal Laplacian in the real space is equivalent to designing an optical system with a response function in the wavevector space (-space) [33]:
where the subscripts denote the and polarizations [Fig. 1(b)]. Thus, to realize an ideal Laplacian, the device should have identical response to both polarizations, and have no polarization mixing in its response function.Instead of realizing an ideal Laplacian as discussed in Eq. (1), in many practical situations it is sufficient to achieve a response function of the form
where . In Section 3 of Supplement 1, we show that a response of the form of Eq. (2) is sufficient to achieve ideal Laplacian operation when the incident light is either unpolarized or circularly polarized. In the case where the incident light is linearly polarized, a device described by Eq. (2) does not perform an ideal Laplacian operation since the decomposition of linearly polarized light into the and polarization is wavevector dependent. Nevertheless, in practice a device with such a response function can still perform detection for edges with arbitrary orientations with linearly polarized incident light, as shown in Section 6 of Supplement 1.In this paper, we show how transmission coefficients in the form of Eq. (2) can be achieved using a photonic crystal slab. Equation (2) requires that at . In a lossless photonic crystal slab, the transmission for normally incident light can vanish near a guided resonance [34]. Thus we consider in more detail the transmission coefficient of a photonic crystal slab near normal incidence. Consider a single photonic band of guided resonances, as characterized by -dependent resonant frequencies and radiative linewidths , with the in-plane wavevector. Near the resonant frequencies, the transmitted amplitude is expressed as [34]
where is the incident light frequency, is the direct transmission coefficient, and is related to the complex decaying amplitude of the resonance to the transmission side of the slab.In general, is constrained by the direct process due to energy conservation and time-reversal symmetry [34,35]. In particular, if , that is, the direct pathway has a 100% transmission coefficient, then
even for structures without z-mirror symmetry, such as shown in Fig. 1, as has been derived in Ref. [35]. In this special case,We denote , . Using Eq. (5), zero transmission occurs at the point when the incident wave has the frequency
At , we perform an expansion of the transmission coefficient near :
where therefore,In this special case, is simply proportional to the band dispersion near . If , then , as well. Moreover, as we will prove, if , then and polarizations are decoupled, i.e., . In this case, Eq. (10) applies to both and , thus Eq. (2) will be realized. Therefore the transmission through a photonic crystal slab can be used to achieve the Laplacian.
The analysis above indicates that, to use a photonic crystal slab as a two-dimensional Laplace operator, it is sufficient that the slab satisfies the following three conditions:
To satisfy the first condition, we note that the direct transmission coefficient is related to the non-resonant transmission pathway [34]. Hence, it is possible to realize by e.g., changing the thickness of the slab. In the structure as shown in Fig. 1, we achieve by placing a uniform dielectric slab in the vicinity of the photon crystal slab, and by tuning the distance between the slabs. This has the advantage that we can tune without significantly affecting the band structure of the photonic crystal slab.
To satisfy the second and third conditions above, one will need to design the band structure for the photonic crystal slab. This design is in fact highly non-trivial due to the vectorial nature of electromagnetic waves. Since the Laplacian is isotropic in -space, it is natural to consider a photonic crystal slab structure that possesses rotational symmetry. As an illustration, here we consider a slab structure with a square lattice of air holes that has symmetry. For concreteness, we consider circular holes; the same analysis will apply to other hole shapes that preserve symmetry. For such a slab, it is known that, at the point, which corresponds to , the only modes that can couple to external plane waves must be twofold degenerate, belonging to a two-dimensional irreducible representation of the group [34,36]. As we consider only small wavevector range near the point, it is sufficient to include only the two modes at the point in our analysis. Near such modes, in the vicinity of the point, the band structure in general can be described by the following effective Hamiltonian (see Section 1 of Supplement 1):
where are three complex coefficients, and the ’s are the Pauli matrices. This Hamiltonian has two eigenvalues ofThe band structure as described by Eq. (12) in general does not satisfy the conditions for the Laplacian as outlined above. In this paper, for concreteness, we will fix the dielectric constant of the material for the slabs to be , which approximates that of Si or GaAs in the infrared wavelength range. In Figs. 2(a)–2(c), we plot the band structure for a slab with a thickness , and a radius of the holes, in the vicinity of the guided resonances at the point with a frequency . The band structure is numerically determined by using the guided-mode expansion method [37,38], and it agrees excellently with the analytic expression of Eq. (12) (see Section 2 of Supplement 1). At the point, the bands are twofold degenerate. Away from the point, the degeneracy is lifted, and the two bands are strongly anisotropic, as can be seen in Fig. 2(a), where the effective mass tensors along the and the directions are drastically different. The band anisotropy can also be visualized in the constant frequency contours plotted in Figs. 2(b) and 2(c) for the two bands. In general, it can be shown that, for a general choice of parameters and , the constant frequency contour as described by Eq. (12) is not a circle even at the limit.

Fig. 2. Band structure of the photonic crystal slab with the dielectric constant , thickness , and radius of the holes. (a)–(c) Band structure near as a typical example of conventional anisotropic bands that are doubly degenerate at . (a) Band dispersion diagram along and . (b) Constant frequency contours of the lower band with respect to . (c) Constant frequency contours of the upper band. (d)–(f) Band structure at , which is nearly isotropic. (d) Band dispersions along and . (e) Constant frequency contours of the lower band. (f) Constant frequency contours of the upper band.
On the other hand, Eq. (12) also indicates that when , both bands will be isotropic:
To achieve such an isotropic band structure requires detailed tuning of the parameters. Numerically, we found that such an isotropic band structure can be approximately achieved with the same slab as indicated above, but near the guided resonances at with a different frequency , where the complex coefficients , , , thus (see Section 2 of Supplement 1 for details). The resulting isotropic band structures are shown in Figs. 2(d)–2(f). In Fig. 2(d), we see that the bands along the and directions have almost identical effective masses, and in Figs. 2(e) and 2(f), we see that the constant frequency contours are almost completely circular.
As we mentioned above, for a structure with symmetry, the modes that a normally incident plane wave can couple to at the point are always twofold degenerate. Consequently, near the point, there are always two bands of guided resonances present. Moreover, for incident wavevectors that are not along a high-symmetry plane, both - and -polarized light may couple to both bands, leading to complicated polarization conversion effects. Thus, in general, in addition to having an anisotropic band structure near , a photonic crystal slab structure also does not satisfy Condition 2 regarding the excitation of a single guided resonance band. Remarkably, however, below we show that once the condition for isotropic bands is satisfied, each of the two bands in fact couple to only one single polarization, for every direction of incidence. Below, we refer to this effect as single-band excitation, and, because of it, and polarizations are decoupled, i.e., in this special case, as required by Eq. (2) (see Section 4 of Supplement 1).
To illustrate the effect of single-band excitation for every direction when the bands are isotropic, we first show that single-band excitation always occurs when the incident direction is in a mirror symmetry plane of the structure. Using the group theory notation in Ref. [39], as our system possesses symmetry, the doubly degenerate states at are modes. Along the direction, this pair of modes splits into two singly degenerate states of and modes, which are even and odd, respectively, with respect to the reflection operator associated with the mirror plane . On the other hand, since the - and -polarized lights are odd and even with respect to the same mirror plane, respectively, along the direction, the -polarized light can couple to only the mode. Thus, in general, we have the effect of single-band excitation when the direction of incidence is in a high-symmetry plane such as the plane [40].
Next, we prove the following statement: if the two-band Hamiltonian is isotropic, that is,
for every , then we have the effect of single-band excitation along all directions. Here, and are the rotation operators that describe the rotation around the axis by an angle , in the -space and the two-dimensional Hilbert space, respectively. To prove this, we denote the eigenstates of the two bands as and , which connect to the and modes along the direction, respectively. Equation (14) implies thatWe denote the - and -polarized modes as and , respectively. By definition,
For along the direction, the mode does not couple to the polarization, i.e.,
Then, for along any direction, using Eqs. (15) and (16), we also have
As a numerical demonstration of the connection of the single-band excitation effect with the isotropic band structure, we consider the two-slab structure as illustrated in Fig. 1. The photonic crystal slab has the same structural parameters as those in Fig. 2. The uniform dielectric slab has a thickness of . The air gap between the two slabs has a width of . Figure 3 shows the transmission coefficients through the two slab system, as we vary the in-plane wavevector, while maintaining ; hence, the direction of the incident light is away from any symmetry plane. Near the frequency of , which is near the guided resonances as shown in Figs. 2(a)–2(c), the band structure is strongly anisotropic, and both the and polarizations excite both bands, as shown in Figs. 3(a) and 3(b). On the other hand, near the frequency of , which is near the guided resonances as shown in Figs. 2(d)–2(f), the bands are isotropic, and each polarization excites only one band, as shown in Figs. 3(c) and 3(d). Moreover, due to the isotropic band structure, the results in Figs. 3(c) and 3(d) are the same for any angle of .

Fig. 3. Transmittances of the two-slab structure as illustrated in Fig. 1. The photonic crystal slab has the same structural parameters as those in Fig. 2. The thickness of the uniform dielectric slab is , and the air gap width is . (a), (b) Transmittance near as a function of and frequency along a general direction for (a) light and (b) light . Both bands couple to and light. (c), (d) Transmittance near along for (c) light and (d) light . light couples with only the upper band while with the lower band, showing the single-band excitation effect.
In this structure, the thickness of the uniform dielectric slab and the gap size are chosen such that near . Therefore, we have designed a structure that satisfies all the sufficient conditions for achieving the Laplace operator.
In Figs. 4(a) and 4(b), we plot the transmission coefficients and , as a function of the in-plane wavevector at the frequency , for the same structure used in Fig. 3. Both and are proportional to . In Fig. 4(c), we plot the transmission coefficient for unpolarized light at the same frequency, which can be derived as (see Section 3 of Supplement 1)
also shows an isotropic response. In Fig. 4(d), we plot as a function of along the direction; due to the isotropic transmittance, the result is the same for any . To show the quadratic dependency, we fit the curve with the quadratic function , where is the only fitting parameter. The fitting is almost perfect for up to , which corresponds to an acceptance angle , and a numerical aperture . As we will show below, this provides a sufficient range of wavevectors for image differentiation. Thus, the transmission coefficients have all the required dependency in order to demonstrate a Laplacian [Eq. (2)].
Fig. 4. Contour plots of transmittance as a function of and at frequency for (a) light , (b) light , and (c) unpolarized light , which are all isotropic near . (d) as a function of along the direction, and the quadratic fitting where is the only fitting parameter.
3. NUMERICAL DEMONSTRATION
We now show that the structure used in Fig. 3 indeed operates as a two-dimensional Laplacian. In the numerical demonstration below, the incident beam is unpolarized light. The transmitted image is calculated in the following way, which is standard in image processing [1]. (1) Compute the Fourier transform of the incident field profile . Note the incident image is . (2) Compute the Fourier transform of the output field profile, . Here is the transmittance for the unpolarized light, which is the transfer function here. (3) Obtain the output field profile by inverse Fourier transform. Calculate the output image .
Figure 5(a) shows the Stanford emblem as the incident image. Figure 5(b) is the calculated transmitted image, which clearly shows all the edges with the same intensity despite different edge orientations. The ability to detect edges with different orientations is a significant advantage of an isotropic differential operator such as the Laplacian, as compared with the one-dimensional differential operators. The Laplacian also highlights the sharp corners, exhibiting its higher sensibility to fine details as a second-order derivative operator [1]. Here we note that a perfect Laplacian, where the response function for all including , would generate double contours for edges. Our device instead maintains only in a wavevector range near , as is upper-bounded by unity due to energy conservation. Thus, the double contours are not observed for the sharp edges in the input image. We also note that there is no diffraction effect, because the periodicity of the photonic crystal slab is smaller than the free space wavelength, and incident waves are near the normal incidence.

Fig. 5. (a) Incident Stanford emblem with a size of . (b) Calculated transmitted image with unpolarized light, which clearly shows the edges with different orientations. (c) Incident slot patterns with length and widths 100, 50, 30, and . (d) Calculated transmitted images with unpolarized light, which show that the spatial resolution of our design is around .
We also test the spatial resolution of such a spatial differentiator. Figure 5(c) shows a series of slot patterns and Fig. 5(d) is the calculated transmitted image. The performance of edge detection degrades as the slot width decreases. The spatial resolution of our differentiator, which is the minimum separation between the two edges that can be resolved, is around . Considering , corresponding to a resonant wavelength , the spatial resolution is 20 μm, which is sufficient for most image-processing applications. (Note ).
4. DISCUSSION AND CONCLUSION
Here we would like to stress that our device operates in the transmission mode, whereas previous works using Fabry–Perot etalon, such as [29], operate in the reflection mode. Operating in the transmission mode is important for image processing, since the device can be directly used as the first layer of an image recognition system. It is more complicated to use a device operating in the reflection mode for this purpose. The Fabry–Perot etalon, intrinsically, cannot operate in the transmission mode: the Laplacian requires the transmission to vanish at a certain , while for an etalon, the transmission does not vanish at any . Our use of a photonic crystal slab thus generally points to new opportunities in optical analog computing as provided by nanophotonic structures.
In conclusion, we have shown that the Laplacian can be implemented using a photonic crystal slab device. Such a simple photonic device may have various applications involving image processing. Future research will be worthy to extend the Laplacian to multi-frequency with multi-layer structures, enabling differentiation of color images.
Funding
Samsung Electronics; U.S. Air Force (USAF) (FA9550-17-1-0002); Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (SNF) (P2ELP2_165174).
Acknowledgment
The authors thank Yu Guo and Dr. Alexander Cerjan for helpful discussions.
See Supplement 1 for supporting content.
REFERENCES
1. R. Gonzales and R. Woods, Digital Image Processing, 3rd ed. (Prentice Hall, 2008).
2. R. Markham, S. Frey, and G. Hills, “Methods for the enhancement of image detail and accentuation of structure in electron microscopy,” Virology 20, 88–102 (1963). [CrossRef]
3. M. D. Abràmoff, P. J. Magalhães, and S. J. Ram, “Image processing with ImageJ,” Biophoton. Int. 11, 36–42 (2004).
4. T. Brosnan and D. W. Sun, “Improving quality inspection of food products by computer vision—a review,” J. Food Eng. 61, 3–16 (2004). [CrossRef]
5. N. Dalal and B. Triggs, “Histograms of oriented gradients for human detection,” in IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR) (IEEE, 2005), pp. 886–893.
6. A. Rosenfeld and A. C. Kak, Digital Picture Processing (Academic, 1982).
7. D. L. Pham, C. Xu, and J. L. Prince, “Current methods in medical image segmentation,” Annu. Rev. Biomed. Eng. 2, 315–337 (2000). [CrossRef]
8. R. J. Holyer and S. H. Peckinpaugh, “Edge detection applied to satellite imagery of the oceans,” IEEE Trans. Geosci. Remote Sens. 27, 46–56 (1989). [CrossRef]
9. D. R. Solli and B. Jalali, “Analog optical computing,” Nat. Photonics 9, 704–706 (2015). [CrossRef]
10. D. Görlitz and F. Lanzl, “A holographic spatial filter for direction independent differentiation,” Jpn. J. Appl. Phys. 14, 223–228 (1975). [CrossRef]
11. J. Eu, C. Liu, and A. Lohmann, “Spatial filters for differentiation,” Opt. Commun. 9, 168–171 (1973). [CrossRef]
12. E. R. Reinhardt and W. H. Bloss, “Optical differential operation processor,” Opt. Eng. 17, 69–72 (1978). [CrossRef]
13. R. Sirohi and V. R. Mohan, “Differentiation by spatial filtering,” Opt. Acta Int. J. Opt. 24, 1105–1113 (1977). [CrossRef]
14. H. Kasprzak, “Differentiation of a noninteger order and its optical implementation,” Appl. Opt. 21, 3287–3291 (1982). [CrossRef]
15. S. K. Yao and S. H. Lee, “Spatial differentiation and integration by coherent optical-correlation method,” J. Opt. Soc. Am. 61, 474–477 (1971). [CrossRef]
16. K. S. Nesteruk, I. P. Nikolaev, and A. V. Larichev, “Image differentiation with the aid of a phase knife,” Opt. Spectrosc. 91, 295–299 (2001). [CrossRef]
17. C. Warde and J. Thackara, “Operating modes of the microchannel spatial light modulator,” Opt. Eng. 22, 695–703 (1983). [CrossRef]
18. J. A. Davis, W. V. Brandt, D. M. Cottrell, and R. M. Bunch, “Spatial image differentiation using programmable binary optical elements,” Appl. Opt. 30, 4610–4614 (1991). [CrossRef]
19. J. Lancis, T. Szoplik, E. Tajahuerce, V. Climent, and M. Fernández-Alonso, “Fractional derivative Fourier plane filter for phase-change visualization,” Appl. Opt. 36, 7461–7464 (1997). [CrossRef]
20. T. Zhu, Y. Zhou, Y. Lou, H. Ye, M. Qiu, Z. Ruan, and S. Fan, “Plasmonic computing of spatial differentiation,” Nat. Commun. 8, 15391 (2017). [CrossRef]
21. N. V. Golovastikov, D. A. Bykov, L. L. Doskolovich, and V. A. Soifer, “Spatiotemporal optical pulse transformation by a resonant diffraction grating,” J. Exp. Theor. Phys. 121, 785–792 (2015). [CrossRef]
22. N. V. Golovastikov, D. A. Bykov, and L. L. Doskolovich, “Resonant diffraction gratings for spatial differentiation of optical beams,” Quantum Electron. 44, 984–988 (2014). [CrossRef]
23. A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alu, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343, 160–163 (2014). [CrossRef]
24. A. Youssefi, F. Zangeneh-Nejad, S. Abdollahramezani, and A. Khavasi, “Analog computing by Brewster effect,” Opt. Lett. 41, 3467–3470 (2016). [CrossRef]
25. S. AbdollahRamezani, K. Arik, A. Khavasi, and Z. Kavehvash, “Analog computing using graphene-based metalines,” Opt. Lett. 40, 5239–5242 (2015). [CrossRef]
26. A. Pors, M. G. Nielsen, and S. I. Bozhevolnyi, “Analog computing using reflective plasmonic metasurfaces,” Nano Lett. 15, 791–797 (2015). [CrossRef]
27. C.-S. Guo, Q.-Y. Yue, G.-X. Wei, L.-L. Lu, and S.-J. Yue, “Laplacian differential reconstruction of in-line holograms recorded at two different distances,” Opt. Lett. 33, 1945–1947 (2008). [CrossRef]
28. J. P. Ryle, D. Li, and J. T. Sheridan, “Dual wavelength digital holographic Laplacian reconstruction,” Opt. Lett. 35, 3018–3020 (2010). [CrossRef]
29. D. A. Bykov, L. L. Doskolovich, E. A. Bezus, and V. A. Soifer, “Optical computation of the Laplace operator using phase-shifted Bragg grating,” Opt. Express 22, 25084–25092 (2014). [CrossRef]
30. E. Chow, S. Y. Lin, S. G. Johnson, P. R. Villeneuve, J. D. Joannopoulos, J. R. Wendt, G. A. Vawter, W. Zubrzycki, H. Hou, and A. Alleman, “Three-dimensional control of light in a two-dimensional photonic crystal slab,” Nature 407, 983–986 (2000). [CrossRef]
31. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light (Princeton University, 2011).
32. W. Zhou, D. Zhao, Y.-C. Shuai, H. Yang, S. Chuwongin, A. Chadha, J.-H. Seo, K. X. Wang, V. Liu, Z. Ma, and S. Fan, “Progress in 2D photonic crystal Fano resonance photonics,” Prog. Quantum Electron. 38, 1–74 (2014). [CrossRef]
33. R. N. Bracewell, The Fourier Transform and Its Applications, 3rd ed. (McGraw-Hill, 1999).
34. S. Fan and J. D. Joannopoulos, “Analysis of guided resonances in photonic crystal slabs,” Phys. Rev. B 65, 235112 (2002). [CrossRef]
35. K. X. Wang, Z. Yu, S. Sandhu, and S. Fan, “Fundamental bounds on decay rates in asymmetric single-mode optical resonators,” Opt. Lett. 38, 100–102 (2013). [CrossRef]
36. T. Ochiai and K. Sakoda, “Dispersion relation and optical transmittance of a hexagonal photonic crystal slab,” Phys. Rev. B 63, 125107 (2001). [CrossRef]
37. L. C. Andreani and D. Gerace, “Photonic-crystal slabs with a triangular lattice of triangular holes investigated using a guided-mode expansion method,” Phys. Rev. B 73, 235114 (2006). [CrossRef]
38. M. Minkov and V. Savona, “Automated optimization of photonic crystal slab cavities,” Sci. Rep. 4, 5124 (2015). [CrossRef]
39. K. Sakoda, Optical Properties of Photonic Crystals, Vol. 80 of Springer Series in Optical Sciences (Springer, 2005).
40. S. G. Tikhodeev, A. L. Yablonskii, E. A. Muljarov, N. A. Gippius, and T. Ishihara, “Quasiguided modes and optical properties of photonic crystal slabs,” Phys. Rev. B 66, 045102 (2002). [CrossRef]