A Photonic Crystal Slab Laplace Differentiator

We introduce an implementation of a Laplace differentiator based on a photonic crystal slab that operates at transmission mode. We show that the Laplace differentiator can be implemented provided that the guided resonances near the $\Gamma$ point exhibit an isotropic band structure. Such a device may facilitate nanophotonics-based optical analog computing for image processing.


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][2][3][4][5] In these applications, of particular importance are isotropic derivative operators, whose response is 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 highthroughput 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][10][11][12][13][14][15][16][17][18][19] Moreover, the recent developments of nanophotonic structures such as meta-surfaces 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][21][22][23][24][25][26] To date, there has been no optical realization of the two-dimensional Laplacian at transmission mode with a compact nanophotonic device.
In this article, we show that a two-dimensional Lapacian can be implemented with the use of a photonic crystal slab, one of the most widely studied nanophotonic structures. [27][28][29] As shown in Figure 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 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 = 12, the geometry parameters are: d = 0.55a, r = 0.111a, d s = 0.07a, d g = 0.21a. The plane above and below shows the input and output image respectively. The transmitted image (e.g. ring) through the device is the Laplacian of the incident image (e.g. disk) illuminated by a normally incident light with frequency ω 0 = 0.47656 × 2πc/a. (b) Coordinate system. (c) The Brillioun zone of the system. isotropic for the two polarizations.

THEORETICAL ANALYSIS
Our objective is to realize a Lapacian for a two-dimensional input field. Specifically, for a normally incident light beam along the z-axis with a transverse field profile S in (x, y), we would like to design an optical device for which the transmitted beam has a profile S out (x, y) ∝ ∇ 2 S in (x, y), whre ∇ 2 = ∂ 2 x + ∂ 2 y is the Laplacian. The task of realizing the Laplacian ∇ 2 in the real space is equivalent to designing an optical system with response function in the wavevector space (k-space). [30] Equation (1) requires that t(k) = 0 at |k| = 0. In a lossless photonic crystal slab, the transmission for normally incident light can vanish near a guided resonance. [31] Thus we consider in more details the transmission coefficient of a photonic crystal slab near normal incidence. Consider a single photonic band of guided resonances, as characterized by k-dependent resonant frequencies ω(k) and radiative linewidths γ(k). (Here k = (k x , k y ) refers to the in-plane wavevector.) Near the resonant frequencies, the transmitted amplitude t is expressed as [31] t where ω is the incident light frequency, t d is the direct transmission coefficient, and f is related to the complex decaying amplitude of the resonance to the transmission side of the slab.
In general, f is constrained by the direct process due to energy conservation and time-reversal symmetry. [31,32] In particular, if t d = 1, that is, the direct pathway has a 100% transmission coefficient, then even for structure without z-mirror symmetry such as shown in Figure 1, as has been derived in Ref. [32]. In this special case, We denote ω 0 = ω(k = 0), γ 0 = γ(k = 0). Using Equation (4), zero transmission occurs at the Γ point when the incident wave has the frequency: At ω = ω 0 , we perform an expansion of the transmission coefficient t near k = 0 : where therefore, In this special case, t(k) is simply proportional to the band dispersion δω(k) near k = 0. If δω(k) ∝ |k| 2 , then t(k) ∝ |k| 2 as well. 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 spatial differentiator, it is sufficient that the slab satisfies the following three conditions: 2. Only one guided resonance band is coupled.
To satisfy the first condition above, we note that the direct transmission coefficient t d is related to the non-resonant transmission pathway [31]. Hence it is possible to realize t d = 1 by, e.g. changing the thickness of the slab. In the structure as shown in Figure 1, we achieve t d = 1 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 t d 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. The design here is in fact quite non-trivial due to the vectorial nature of electromagnetic waves. Since the Laplacian is isotropic in k-space, it is natural to consider a photonic crystal slab structure that has rotational symmetry. As an illustration, here we consider a slab structure with a square lattice of air holes that has C 4v symmetry. For such a slab, it is known that at the Γ point, which corresponds to |k| = 0, the only modes that can couple to external plane wave must be two-fold degenerate, belonging to a two-dimensional irreducible representation of the C 4v group. [31,33] Near such modes, in the vicinity of Γ point, the band structure in general can be described by the following 2 × 2 effective Hamiltonian: (See Supplementary Materials) where a, b, c are three complex coefficients, the σ's are the Pauli matrices. This Hamiltonian has two eigenvalues of The band structure as described by Equation (11) in general does not satisfy the conditions for ideal spatial differentiation as outlined above. In this paper for concreteness we will fix the dielectric constant of the material for the slabs to be = 12, which approximates that of Si or GaAs in the infrared wavelength range. In Figure 2(a-c), we plot the band structure for a slab with a thickness d = 0.55a, and a radius r = 0.111a of the holes, in the vicinity of a guided resonance at the Γ point with a frequency ω 0 = 0.38749 × 2πc/a. The band structure is numerically determined using guided-mode expansion method [34,35], and it agrees excellently with the analytic expression of Equation (11) (See Supplementary Materials). At the Γ point, the bands are two-fold degenerate. Away from the Γ point, the degeneracy is lifted, and the two bands are strongly anisotropic, as can be seen in Figure 2(a), where the effective masses along the Γ-X and the Γ-M direction are drastically different. The band anisotropy can also be visualized in the constant frequency contour plotted in Figure 2(b) and (c), for the two bands. In general, it can be shown that, for a general choice of parameters a, b and c, the constant frequency contour as described by Equation (11) is not a circle even at the |k| → 0 limit.
On the other hand, Equation (11) also indicates that when c = ±2b, both bands will be isotropic: To achieve such an isotropic band structure requires detailed tuning of the parameters. Numerically, we found that such an  Figure 2(d-f). In Figure 2(d), we see that the bands along the Γ-M and Γ-X direction have almost identical effective masses. And in Figure 2(e) and (f), we see that the constant frequency contours are almost completely circular. As we mentioned above, for a structure with C 4v symmetry, the modes that a normally-incident plane wave can couple to at the Γ point are always two-fold degenerate. Consequently, near the Γ point, there are always two bands of guided resonance present. Moreover, off the normal direction, if the direction of incident waves is away from the high symmetry planes, both S and P polarized lights may couple to both bands, leading to complex polarization conversion effects. Thus, in general, in addition to having an anisotropic band structure near Γ, a photonic crystal slab structure also does not satisfy the condition 2 above regarding the excitation of a single guided resonance band. Remarkably, however, below we show that once the condition for an isotropic band is satisfied, each of the two bands in fact only couple to one single polarization, for every direction of incidence. Below, we refer to this effect, where each polarization excites only a single band away from normal incidence, as the effect of single-band excitation.
To illustrate the effect of single-band excitation for every direction when the bands are isotropic, we first show that singleband excitation always occurs when the incident direction is in a mirror symmetry plane of the structure. Using the group theory notation in Ref. [36], as our system possesses C 4v symmetry, the doubly degenerate states at Γ are E modes. Along the Γ-X direction, this pair of E modes splits into two singly degenerate states of A and B modes which are even and odd, respectively, with respect to the reflection operator associated with the mirror plane y = 0. On the other hand, since the S and P polarized lights are also odd and even with respect to the same mirror plane, respectively, along the Γ-X direction, the S(P) polarized light can only couple to the B(A) modes. 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 y = 0 plane. [37] Next, we prove the following statement: if the two-band Hamiltonian is isotropic, that is, for every ϕ ∈ (0, 2π), then we have the effect of single-band excitation along all directions. Here, R(ϕ) andR(ϕ) are the rotation operators that describe the rotation around the z-axis by an angle ϕ, in the k space and the two-dimensional Hilbert space, respectively. To prove this, we denote the eigenstates of the two bands as |k, A and |k, B , which connect to the A and B modes along the Γ-X direction respectively. Equation (13) implies that: We denote the S and P polarized modes as |k, S and |k, P respectively. By definition, Then for k = R(ϕ)k x along any direction, using Equations (14) and (15), we also have As a numerical demonstration of the connection of the singleband excitation effect with the isotropic band structure, we consider the two-slab structure as illustrated in Figure 1. The photonic crystal slab structure has the same parameters as those of Figure 2. The uniform dielectric slab has a thickness of d s = 0.07a. The air gap between the two slabs has a width of d g = 2.1a. Figure 3 shows the transmission coefficients through the two slab system, as we vary the in-plane wavevector, while maintaining ϕ = 14 • , and hence the direction of the incident light is away from any symmetry plane. Near the frequency of ω = 0.38749 × 2πc/a, which is near the guided resonances as shown in Figure 2(a-c), the band structure is strongly anisotropic, and both the S and P polarizations excite both bands as shown in Figure 3(a-b). On the other hand, near the frequency of ω = 0.47656 × 2πc/a, which is near the guided resonances as shown in Figure 2(d-f), each polarization excites only one band as shown in Figure 3(c-d). In this structure, the thickness of the uniform dielectric slab and the gap size are chosen such that t d = 1 near ω = 0.47656 × 2πc/a. Therefore, we have designed a structure that satisfies all the sufficient conditions for achieving an ideal Laplacian spatial differentiator.
In Figure 4(a) and (b), we plot the transmission coefficients, as a function of in-plane wavevector k, for the S and P polarized incident light, at the frequency ω = 0.47656 × 2πc/a, for the same structure used in Figure 3. Both polarizations exhibit a transmission coefficient that is proportional to |k| 2 . In Figure 4(c), we plot the transmission coefficient for unpolarized light at the same frequency. The transmission coefficient for the unpolarized light can be derived as (See Supplementary Materials for a derivation) Here t s and t p are the transmission coefficients for S and P light respectively. The transmission coefficient for the unpolarized light also shows an isotropic response. In Figure 4(d), we plot |t u | as a function of |k|, which shows a quadratic dependency. To show this, we fit the curve with quadratic function |t u | = α|k| 2 , where α is the only fitting parameter. The fitting is almost perfect for |k| up to 6 × 10 −3 × 2π/a, which as we will show provides a sufficient range of wavevector for image differentiation. Thus, the transmission coefficients have all the required k-dependency in order to demonstrate a Laplacian.

NUMERICAL DEMONSTRATION
We now show that the structure used in Figure 3 indeed operates as a two-dimensional Laplacian differentiator. Figure 5(a) is 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 two-dimensional Lapacian, 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] We also test the spatial resolution of such spatial differentiator. Figure 5(c) is a series of slot patterns and Figure 5(d) is the calculated differentiated 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 30a. Considering a = 0.67µm corresponding to the resonant wavelength λ 0 = 1.4µm, the spatial resolution is 20µm, which is sufficient for most image processing application.

DISCUSSION AND CONCLUSION
In previous works, the two-dimensional Laplace operator has been implemented using holograms [38,39] and phase-shifted Bragg grating [40]. These implementations either rely upon bulky optical components, or operate in reflection modes. In contrast, here we show this operation can be achieved with a compact optical structure in transmission mode, which is more compatible for image processing applications.
In conclusion, we have shown that a Laplace differentiator can be implemented using a photonic crystal slab. Such a simple photonic device may have various applications involving image processing. Future research will be worthy to extend the Laplace differentiator to multi-frequency with multi-layer structures, enabling differentiation of color images.

FUNDING INFORMATION
This work is supported in part by Samsung Electronics, and by the U. S. Air Force Grant No. FA9550-17-1-0002.

DERIVATION OF THE EFFECTIVE HAMILTONIAN
In this section, we derive the 2 × 2 effective Hamiltonian (Equation (10) in the primary text) near the Γ point. We assume that the system has C 4v symmetry. And at the Γ point the system supports a pair of doubly degenerate states, denoted as |x and |y , respectively. With these states as bases, the 2 × 2 Hamiltonian in the vicinity of Γ has the following general form: whereÂ(k) andB(k) are both Hermitian and We are interested in the lowest-order non-vanishing terms in the functions f , g, h, r, s and t. By virtual of the two-fold degeneracy, we have Also, g, h, s and t are real sinceÂ(k) andB(k) are Hermitian. In general, the HamiltonianĤ(k) is constrained by the symmetric group G of the system: whereD(g) and D(g) are the representations of g in the Hilbert space and the k space, respectively. For our system, G = C 4v = {E, 2C 4 , C 2 , 2σ v , 2σ d }. As we choose |x and |y as bases for the Hilbert space and k x and k y for the k space, the representationsD(g) and D(g) have the same matrix forms. Now, we consider the constraint onĤ(k) from each symmetry element in G.

Inversion symmetry
From Equation (S4) we have: Thus, the lowest order terms in f , g, h, r, s and t must all be quadratic.
2. Four-fold rotational symmetry C 4 .D From Equation (S4) we have: Here we consider the mirror plane perpendicular to the y-axis. From Equation (S4) we have: The other mirror symmetry σ d doesn't provide new constraints as it can be realized with a combination of σ v and C 4v . Combining all the symmetry requirements, we determine the forms of the lowest order terms in all the functions. As an example, we consider the function f (k), which can be expanded as where all the C coefficients are in general complex. Since f (k) = − f * (Rk) = − f (M y k), we have C x = C y = 0, C xy = C * xy ≡ C. Thus f (k) = Ck x k y , where C is real. Repeating the procedures for all the functions, we summarize the results below: where all the coefficients are real. Therefore, we obtain the Hamiltonian in Equation (10) of the main text.: where a = A − iA , b = B − iB , c = C − iC are the complex coefficients.

VALIDATION OF THE EFFECTIVE HAMILTONIAN
In this section we provide quantitative validation of the effective Hamiltonian. Figure (S1) plots the band structure of the photonic crystal slab near frequency ω 0 = 0.38749 × 2πc/a. Figure (S1) (a-d) exhibit contour plots of (ω − ω 0 ) and (γ − γ 0 ) for the two bands numerically calculated using guided-mode expansion method, while (e-h) show the corresponding analytical results from the effective Hamiltonian (Equation (S13)). The complex coefficients a = 0. 33

TRANSMITTANCE OF ARBITRARY POLARIZED STATES
In this section we compute the transmission of an arbitrarily polarized state, given the transmission response of the system for S and P polarized light. The transmission response of the system is described by the S matrix: where t ss , t sp , t ps and t pp are complex. Here we adopt |s and |p as the basis states. For a general input state |e , the output state is S |e . Consider an input beam with its polarization described by a polarization density matrix ρ, normalized such that I = tr ρ is the incident light intensity. The transmitted state is then described by ρ = SρS † , with the transmission defined as:   S2. Band structure of the photonic crystal slab near frequency ω 0 = 0.47656 × 2πc/a. The structure has the dielectric constant = 12, the thickness d = 0.55a, and the radius r = 0.111a of the holes. (a-d) Contour plots of the band structure numerically calculated using guided-mode expansion method. (a) (ω − ω 0 )-(k x , k y ) for the lower band. (b) (ω − ω 0 )-(k x , k y ) for the upper band. (c) (γ − γ 0 )-(k x , k y ) for the lower band. (d) (γ − γ 0 )-(k x , k y ) for the upper band. γ 0 = 1.3 × 10 −4 × 2πc/a. The (ω − ω 0 ) contours for both bands are almost completely circular, while the (γ − γ 0 ) contours are anisotropic. Nonetheless, (γ − γ 0 ) are much smaller than (ω − ω 0 ) and thus don't affect the isotropy of the band structure much. (e-h) Corresponding contour plots of the band structure analytically calculated from the effective Hamiltonian (Equation (S13)). The complex coefficients a = 0.68 − 0.14i, b = 1.11 − 0.12i, c = 2.24 + 0.15i are fitted from the band dispersions along Γ-X and Γ-M direction. The numerical results agree with the analytical expressions excellently. In all the plots k x and k y are in the units of 10 −3 × 2π/a, (ω − ω 0 ) and (γ − γ 0 ) are in the units of 10 −4 × 2πc/a.
Note that for our system with an isotropic band structure, S is diagonal: Then the interference term in Equation (S20) disappears and we have |t l | = |t u |. So the differentiation performance will be the same under illumination of unpolarized or circularly polarized light.

TRANSMITTANCE AT LARGER WAVEVECTORS
Figure (S3) plots the transmittance |t| as a function of k x and k y for S, P and unpolarized light in a larger wavevector region as compared to Figure 4 of the primary text. All |t| are isotropic near the Γ point. They become anisotropic at larger wavevector k. The background transmittances are all 1. The occurrence of transmittance dips at larger k are due to other guided resonance modes. We calculate the performance of our differentiator based on the transmittance shown here, which covers all the relevant wavevector regions for the incident images that we assumed.  S3. Transmittance |t| as a function of k x and k y at frequency ω 0 = 0.47656 × 2πc/a for (a) S, (b) P and (c) unpolarized light, which are all isotropic near Γ, which is the major relevant wavevector region for image processing, but anisotropic at larger k. k x and k y are in the units of 10 −2 × 2π/a.