Comparison between the linear and nonlinear homogenization of graphene and silicon

In this paper, we use a versatile homogenization approach to model the linear and nonlinear optical response of two metasurfaces: a plasmonic metasurface consisting of a square array of graphene cruciform patches and a dielectric metasurface consisting of a rectangular array of photonic crystal (PhC) cavities in a silicon PhC slab waveguide. The former metasurface is resonant at wavelengths that are much larger than the graphene elements of the metasurface, whereas the resonance wavelengths of the latter one are comparable to the size of its resonant components. By computing and comparing the effective permittivities and nonlinear susceptibilities of the two metasurfaces, we infer some general principles regarding the conditions under which homogenization methods of metallic and dielectric metasurfaces are valid. In particular, we show that in the case of the graphene metasurface the homogenization method describes very well both its linear and nonlinear optical properties, whereas in the case of the silicon PhC metasurface the homogenization method is less accurate, especially near the optical resonances


I. INTRODUCTION
Metamaterials, whose emergence has opened up exciting new opportunities to create novel media with pre-designed physical properties, have been proving to have a significant impact on the development of new approaches and devices for controlling light interaction with matter and achieving key functionalities, including light focusing 1-3 , perfect lensing 4,5 , perfect absorption 6-10 , electromagnetic cloaking [11][12][13] , imaging with sub-diffraction resolution [14][15][16][17][18] , and optical sensing [19][20][21][22] . One of the most important functionalities provided by metamaterials is the enhancement of the local optical field [23][24][25][26][27][28][29][30] . This feature is particularly relevant to nonlinear optics since nonlinear optical interactions grow nonlinearly with the applied field. Promising applications of metamaterials can be found in a broad area of science and engineering, including optical filters, sensing and infrastructure monitoring, medical devices, remote aerospace applications, and smart solar power management.
As research in metamaterials advanced, it became clear that the two-dimensional (2D) counterpart of metamaterials, the so-called metasurfaces, would offer the fastest route to functional devices and applications. This is so chiefly because most nanofabrication techniques can conveniently be applied to the planar configuration of metasurfaces. These ultrathin and lightweight optical devices are generally made of sub-wavelength dielectric or metallic elements arranged in one-dimensional (1D) or 2D periodic arrays. Equally important from a practical perspective, the single-layer characteristics of photonic devices based on metasurfaces make them particularly amenable to system integration. Because of their small thickness, light-matter interaction occurs in a reduced volume and as such optical losses in metasurfaces are relatively small. Importantly in nonlinear optics applications, this reduced light propagation distance in meta-surfaces means that phase-matching requirements can be relaxed, which greatly reduces the design constraints of nonlinear optical devices based on metasurfaces [31][32][33][34][35][36] .
Metasurfaces can mainly be divided into two categories, namely metallic (plasmonic) 37,38 and dielectric metasurfaces 39,40 . Plasmonic metasurfaces, which exploit the resonant excitation of surface plasmons at specific frequencies [41][42][43] , can greatly enhance the local optical field, but this phenomenon is usually accompanied at optical frequencies by large dissipative losses. Dielectric metasurfaces, on the other hand, experience much smaller optical losses but only provide limited optical field enhancement. Moreover, another difference between the two classes of metasurfaces, which is directly related to the magnitude of the optical losses, is that whereas the resonances in the plasmonic metasurfaces are relatively broad, the (Mie) resonances associated to dielectric metasurfaces are particularly narrow. As a result, dielectric metasurfaces are usually much more dispersive than the plasmonic ones. One effective approach to study metasurfaces is using homogenization methods, which reduce the metasurface to a homogeneous material with specific linear and nonlinear retrieved optical coefficients. These effective physical quantities are determined in such a way that the metasurface and homogenous layer of material have the same optical response.
In this paper, we propose a homogenization method and investigate its accuracy when applied to plasmonic and dielectric metasurfaces. As plasmonic metasurface we consider a graphene metasurface [44][45][46] consisting of a square array of free-standing graphene cruciform patches, whereas the dielectric metasurface is made of a rectangular array of photonic crystal (PhC) cavities possessing high-Q optical modes, embedded in a silicon PhC slab waveguide. The rationale for our choice is that the two metasurfaces capture the general characteristics of the two main classes of metasurfaces, i.e. plasmonic and dielectric metasurfaces. In particular, the cruciform graphene patches possess strong plasmon resonances characterized by highly confined, enhanced optical near-field. As optical nonlinearity, we consider second-harmonic generation (SHG) by the nonlocal nonlinear polarization, as symmetry considerations imply that the SHG by the the local nonlinear polarization is zero 47,48 . Moreover, the silicon PhC cavities were designed so as to possess two high-Q optical cavity modes separated by the Raman frequency of silicon, which ensures a strong effective Raman nonlinearity of the metasurface 49 .
The remainder of the paper is organized as follows: In the next section we present the geometrical configurations and material parameters characterizing the two metasurfaces considered in this work. Then, in Sec. III, we introduce the linear and nonlinear homogenization method used to extract the constitutive parameters of the metasurfaces, whereas in Sec. IV we apply our homogenization approach to the two metasurfaces and derive general principles regarding the conditions in which the predictions of the homogenization method are valid. In particular, we extract the linear and nonlinear constitutive parameters of the metasurfaces and then compare the optical response of the metasurfaces with that of their homogenized counterparts. Finally, we conclude our paper by summarizing the main findings of our study and discussing some of their implications to future developments pertaining to metamaterials research.

II. DESCRIPTION OF THE GRAPHENE AND SILICON METASURFACES
In this section, we describe the geometrical configuration and material parameters of the two metasurfaces investigated in this work, namely the graphene cruciform metasurface illustrated in figure 1(a) and the silicon PhC nanocavity metasurface shown in figure 1(b). In addition, we explain the rationale for our choice of metasurfaces by presenting their main physical properties.
A. Geometrical configuration of the graphene and silicon metasurfaces The graphene metasurface lies in the x-y plane and consists of a square array of cruciform graphene patches. The symmetry axes of the array coincide with the x-and y-axes and are oriented along the arms of the crosses, as per figure 1(a). The length and width of the arms of the crosses along the two axes are L x (L y ) and w x (w y ), respectively, whereas the corresponding periods of the metasurface are P x and P y . In this, work, unless otherwise specified, L x = L y = 60 nm, w x = w y = 30 nm, and P x = P y = 100 nm.
The relative electric permittivity of graphene is given by the relation: where h g = 0.3 nm is the thickness of graphene, ω is the frequency, and the graphene surface conductivity, σ s , is described by the Kubo's formula 50 : Here, µ c , T , and τ are the chemical potential, temperature, and relaxation time, respectively. In this study, we use µ c = 0.6 eV, τ = 0.1 ps, and T = 300 K. The wavelength dependence of the real and imaginary parts of graphene permittivity, as described by Eq. (1) in conjunction with Eq. (2), are depicted in figure 2. It can be seen in this figure that Im(ε g ) > 0 and Re(ε g ) < 0, which are spectral characteristics shared by most noble metals.
The silicon metasurface, illustrated in figure 1(b), consists of a rectangular array of PhC cavities embedded in a PhC slab waveguide made of silicon (n Si = 3.4) 49 . The PhC slab waveguide comprises a 2D hexagonal lattice of air holes in a silicon slab with the hole radius r = 0.29a and slab thickness t = 0.6a, where a is the lattice constant. Moreover, the optical nanocavities are the so-called L5 PhC cavities, namely they are created by filling in 5 consecutive holes oriented along the ΓK symmetry axis of the hexagonal lattice. The periods of the PhC metasurface, defined as the center-to-center distance along the x-and y-axes between adjacent PhC cavities, are d l = 17a and d t = 6 √ 3a, respectively. In order to increase the Q-factor of the optical modes of this PhC nanocavity, we shifted outwardly the end-holes of the cavity by 0.15a 51 .
The PhC cavity is designed in such a way that it possesses two optical modes whose frequencies are separated by the Raman frequency of silicon, Ω = 2π × 15.6 THz 52 . This ensures a very strong nonlinear coupling between the two optical modes, both because of the large optical field enhancement inside the cavities and also due to favorable spatial overlap between the two optical modes. Consequently, one can achieve an efficient Raman interaction between the two optical modes. This means that the PhC nanocavities can be viewed as artificially engineered, strongly nonlinear "silicon meta-atoms", which when arranged in some spatial pattern give rise to photonic metasurfaces with large Raman nonlinearity. In particular, if one chooses the lattice constant a = 333 nm, the resonance frequency of the pump and Stokes modes are ω p = 1572.5 THz and ω S = 1474.6 THz, respectively, and therefore the condition ω p − ω S = Ω is fulfilled. Expressed in terms of normalized frequency of 2πc/a, the frequencies of the two optical modes are ω p = 0.2778 and ω S = 0.2605. Also, the Q-factors of the two modes are Q p = 1804 and Q S = 1.12 × 10 5 . Note that the photonic band structure of the PhC slab and the corresponding cavity modes were computed with RSoft's BandSOLVE 53 , with the cavity modes lying in the transverse-magnetic bandgap of the PhC slab waveguide (see figure 3).
The linear optical response of graphene and silicon metasurfaces is presented in figures 4(a) and 4(b), respectively, and correspond to an x-polarized plane wave normally incident onto the metasurface. This figure clearly shows that, in the case of the graphene metasurface, absorption, A, reflectance, R, and transmittance, T , have a series of spectral resonances occurring at common wavelengths. These resonances are due to the generation of localized surface plasmons in the graphene crosses, the resonance wavelength of the fundamental plasmon being λ = 5.88 µm. Moreover, it can be seen that the excitation of surface plasmons is accompanied by a large increase of the optical absorption, which suggests that at the corresponding resonance wavelengths the optical near-field is significantly enhanced. Moreover, the spectra of the silicon metasurface show two resonances at the frequencies of the two cavity modes, the width of these resonances being much smaller than that of plasmon resonances.

III. LINEAR AND NONLINEAR HOMOGENIZATION METHOD
In this section, we present first the homogenization approach used to describe the effective optical response of the metasurfaces investigated in this paper and retrieve their effective permittivities and nonlinear susceptibilities and then compare the linear optical response of the original and homogenized metasurfaces so as to clarify the circumstances in which our homogenization approach is accurate.

A. Theory of the effective permittivity of metasurfaces
The general linear homogenization method presented here amounts to establishing a relationship between the averaged electric displacement field, D, and the electric field, E, and is known as the field averaging method. It will be expanded later on to the nonlinear case. As it is well known, the constitutive relation of a linear anisotropic material is expressed as D i = j ij E j , where ij is the permittivity tensor and the subscripts i, j = x, y, z. The field-average method relies on this relation, the effective permittivity of the medium being defined via a similar relation between the electric field and the electric displacement field: This method is particularly suitable for describing metasurfaces, because in this case the averaging domain can be naturally defined as the unit cell of the metasurface. Therefore, in the equations above, V is the volume of the unit cell of the metasurface. According to the field averaging method, for an isotropic medium whose permittivity tensor is diagonal, the effective permittivity is evaluated as i = D i /E i . This definition can be extended to anisotropic media as follows: First, one defines an auxiliary quantity, d ij = ij E j , and express the displacement field as D i = j d ij , and then calculate the averaged value of each component of the auxiliary quantity: If we assume that the averaged fields are related by a constitutive relation similar to that corresponding the local fields, namely D i = j ij E j , and by requiring that the average of the field D(r, ω) of the metasurface and the field D(ω) in the homogenized layer of material are termwise equal, the effective permittivity is determined by the following equation: This formula has been used to determine the effective permittivity of both metasurfaces. Before moving on to the calculation of the effective nonlinear susceptibilities of the metasurfaces, we note that in the case of the graphene metasurface the volume integrals can be reduced to surface integrals across the midsection of the graphene sheet because the fields across the ultrathin graphene layer vary only slightly.
B. Calculation of effective second-order susceptibility of graphene metasurfaces As we assume that the graphene crosses are free standing, symmetry considerations imply that the dipole (local) nonlinear polarization at the second harmonic (SH) exactly cancels. Consequently, the lowest order SHG is due to the nonlocal nonlinear polarization whose sources are magnetic dipoles and electric quadrupoles oscillating at the SH frequency. It should be noted that if the graphene patches lie onto a substrate the inversion symmetry is broken and consequently the generated SH is due to the dipole nonlinear polarization. The homogenization of such graphene metasurfaces has been recently studied 54 .
This nonlinear polarization and the associated nonlinear surface current in the case of graphene metasurfaces characterized by nonlocal nonlinear polarization can be expressed as 47 : J s (r, Ω) = σ (2) g (Ω; ω, ω) . . .E(r, ω)∇E(r, ω), where χ (2) g (Ω; ω, ω) and σ (2) g (Ω; ω, ω) are the bulk nonlinear second-order susceptibility and surface nonlinear second-order optical conductivity, respectively, and are related by the following formula: We stress that instead of a bulk nonlinear susceptibility one can use a surface one, defined as χ g , but we decided to use bulk quantities so that it is more convenient to compare these nonlinear optical coefficients of graphene with those of other centrosymmetric materials.
The surface nonlinear second-order optical susceptibility of graphene has been recently derived in 55 and is given by the following equation: where the scalar part of the surface second-order conductivity tensor is: Componentwise, the nonlinear polarization can be evaluated as: where we introduced a new nonlinear auxiliary quantity defined as q ijkl = 0 χ (2) g,ijkl E j ∇ k E l . Moreover, the spatial average of this quantity is: where χ (2) ijkl (r) = χ (2) g,ijkl if r corresponds to a point inside the graphene crosses and χ (2) ijkl (r) = 0 if r is in the air region.
If we express the nonlinear polarization in the homogenized metasurface as: where χ (2) ijkl is the effective nonlinear second-order susceptibility of the homogenized metasurface, and impose the condition that the spatial average of the nonlinear polarization described by Eq. (12) is termwise equal to the polarization in Eq. (14), we obtain the following formula for the effective nonlinear susceptibility: In this formula and in Eq. (14), the quantity ∇ i E j is defined as: C. Theory of effective Raman susceptibility of silicon metasurfaces The calculation of the effective Ramman susceptibility of the silicon based PhC metasurface described in Sec. 2 is similar to that of the effective second-order susceptibility presented in the preceding subsection, so that here we present only the main steps. A more detailed derivation can be found in 49 .
We start our analysis with the nonlinear Raman polarization at the Stokes frequency, P R (r, ω S ), which can be written as: R (r) is the Raman susceptibility and E(r, ω p ) and E(r, ω S ) are the optical fields at the pump and Stokes frequencies, respectively. For the sake of simplicity, we assume that the symmetry axes of the array of PhC cavities coincide both with the x-and y-axes and with the principal axes of silicon. Under these circumstances, the nonzero components of χ R,ijji , with i, j = x, y, z and i = j. The value at resonance of the only independent component is χ We then define the spatially averaged effective Raman polarization: where the volume integration is taken over the unit cell of the metasurface, together with the effective Raman polarization in a homogenized slab of nonlinear optical medium with the same thickness as that of the PhC slab: Here, χ R is the effective Raman susceptibility of the homogenized metasurface.
As in the case of the graphene metasurface, we cannot simply impose the condition that the components of the nonlinear polarizations in Eq. (18) and Eq. (19) are equal because in the general case the effective Raman susceptibility tensor, χ (3) R , has 81 independent components, so that the corresponding system of equations is overdetermined. Consequently, we impose the condition that the r.h.s. of equations Eq. (18) and Eq. (19) are termwise identical. Using this constraint, it can be seen that the components of χ (3) R are determined by the following relations: . (20) Note that the components of χ  R cancel for the same set of indices i, j, k, and l.

IV. RESULTS AND DISCUSSIONS
In this section, we study the circumstances in which our method produces accurate results and use it to understand the main differences between the physical properties of graphene and silicon PhC metasurfaces. Based on our theoretical analysis, the effective permittivity of the two metasurfaces investigated in this work can be calculated using Eq. (6). The results of our calculations, corresponding to the graphene and silicon PhC metasurfaces, are presented in Fig. 5(a) and Fig. 5(b), respectively. These figures show some similarities between the two spectra but also significant differences. Thus, the effective permittivity of graphene metasurface displays a series of spectral resonances of Lorentzian nature, which occur at the plasmon resonance wavelengths of the graphene crosses. This means that the graphene crosses behave as meta-atoms that possess a series of resonant states, the overall optical response of the metasurface being primarily determined by these resonances. The main reason for this behavior can be traced to the size of the graphene crosses relative to the resonance wavelengths of the plasmons of the graphene crosses. Specifically, since the size of the crosses is much smaller than the plasmon wavelengths, the overall optical response of the graphene metasurface can be viewed as a superposition of the response of weakly interacting Lorentz-type oscillators.
The spectrum of the effective permittivity of the silicon PhC metasurface, on the other hand, presents a series of complex features, which are the result of several phenomena. Thus, the two main resonances of the effective permittivity are due the excitation of the two optical modes of the PhC cavity. The spectral separation between the frequencies of the two cavity modes is relatively small and this leads to interference features in the spectrum of the effective permittivity. The other, weaker spectral peaks are presumably due to leaky modes of the PhC slab.

B. Validation of the homogenization approach
In order to validate the conclusions drawn in the preceding subsection and to investigate the situations in which our homogenization method is accurate, we calculated the absorption, A, reflectance, R, and transmittance, T of both metasurfaces and their homogenized counterparts. The main results of these calculations are summarized in Fig. 6 and they reveal several important ideas. Thus, it can be seen in this figure that in the case of the graphene metasurface the resonances of the transmittance occur at the same wavelengths as the resonances of the effective permittivity, whereas for the silicon PhC metasurface the the two sets of resonances differ to some extent. In order to explain these results, one should note that generally the resonances of the transmittance of a planar optical system correspond to excitation of bound modes of the optical system. In the cases investigated here, the bound modes are localized surface plasmons of graphene crosses for the graphene metasurface and optical modes of the PhC cavities and leaky modes of the PhC slab. Moreover, the excitation of these resonances induces a resonant response of the medium, via the polarization of the medium, which in turn leads to resonances in the spectrum of the effective permittivity of the homogenized metasurfaces. Fig. 6(b) also shows that whereas the linear optical responses of the graphene metasurface and its homogenized counterpart are almost identical, in the case of the silicon PhC metasurface they markedly differ from each other. The main reason for this dichotomy is that the graphene crosses are much smaller than the operating wavelength, which makes them respond to the incident optical field as if they were point-like resonators. By contrast, the size of the PhC cavities is comparable to the resonance wavelength of the cavity modes, which renders our homogenization method to be less accurate in this case. It should also be noted, however, that although the homogenization approach for silicon PhC metasurface cannot provide extremely accurate quantitative values for the effective permittivity, it can still provide us valuable qualitative insights into the governing physics.
These ideas are further illustrated by the dependence of the linear optical response of the graphene metasurface on the angle of incidence, θ, of the incoming plane wave, which is presented in Fig. 7. Thus, it can be observed that the spectral resonances of the graphene metasurface, calculated for θ = 0 • , 30 • and 60 • , only slightly varies with θ, whereas the values of A, R, and T at the resonance wavelengths depend more pronouncedly on θ. These findings are explained by the fact that the plasmon resonances depend chiefly on the shape of the graphene nano-patches, and thus are independent on θ, whereas the particular values of A, R, and T depend on the coupling between the incident field and graphene crosses, more specifically on the spatial overlap between the incident wave and the optical field of the graphene plasmons, which is obviously θ-dependent. mons at the fundamental frequency (FF) induces a strong optical near-field and consequently enhanced nonlinear polarization, which is the source of the generated SH. This phenomenon can be clearly seen in Fig. 8, where we present the absorption spectra at the FF and SH. In particular, it can be observed in this figure that the occurrence of a resonance at the FF is accompanied by a resonance at half of its wavelength in the SH spectrum. For example, the plasmon resonance at the FF of λ F F = 5.878 µm has a correspondent in the SH spectrum at λ SH = λ F F /2 = 2.939 µm.
As can be easily inferred from Eq. (9) and Eq. (10), the dominant components (in absolute value) of the second-order susceptibility of graphene are χ (2) g,xxyy = χ (2) g,yyxx . The value of this component, determined for a FF wavelength λ = 1 µm, is χ (2) g,xxyy = (−8.37 + 0.133i) × 10 −19 m 2 V −1 . Moreover, as a consequence of our approach to the calculation of the effective second-order susceptibility, it is equal to zero for the same set of indices for which the graphene secondorder susceptibility is equal to zero. Therefore, in order to quantify the enhancement of the nonlinear optical response of the graphene metasurface, we computed the enhancement factor, η SH = |χ (2) xxyy /χ (2) g,xxyy |, for several values of the angle of incidence, θ = 0 • , 30 • and 60 • . We summarize the results of these calculations in Fig. 9.
The most important conclusion of this analysis is that, at the wavelength of the fundamental plasmon, the effective second-order susceptibility of the homogenized graphene metasurface is enhanced by almost 200×. More precisely, the maximum achievable enhancement factor is η SH = 175, which means that the maximum value of the effective second-order susceptibility of the homogenized graphene metasurface is |χ (gold at λ = 810 nm) 57 and γ = 1.3 × 10 −19 m 2 V −1 (silicon at λ = 800 nm) 58 . Moreover, it can be inferred from Fig. 9 that the enhancement factor is smaller for higherorder plasmons, as in this case the plasmon-induced field enhancement decreases. Another feature revealed by Fig. 9 is that the enhancement factor decreases as the angle of incidence increases, a finding explained by the fact that when θ increases the spatial overlap between the incident wave and the plasmon mode becomes less favorable and thus the enhancement of the local optical field decreases.

D. Effective Raman susceptibility of the silicon photonic crystal metasurface
Due to the symmetry properties of the silicon PhC metasurface and the orientation of the cavity array with respect to the principal axes of silicon, the only non-zero component of the effective Raman susceptibility of the metasurface is χ R . Therefore, similarly to the case of the graphene metasurface, we define the enhancement factor η R = |χ is the dominant component of the Raman susceptibility of silicon. The parameter η R quantifies the enhancement of the Raman nonlinearity of the silicon PhC metasurface. Moreover, in order to investigate the dependence of the enhancement factor on the angle of incidence, these calculations were performed for θ = 0 • , 30 • and 60 • .
Following the procedure we just described, we found out that for the values of the incidence angle of 0 • , 30 • and 60 • , the enhancement factor was 2.29 × 10 4 , 3.19 × 10 3 and 1.99 × 10 3 , respectively. Thus, it can be seen that a giant enhancement of the effective Raman susceptibility of the silicon PhC metasurface of more than 4 orders of magnitude can be achieved at normal incidence. In order to understand the main reason for this remarkable nonlinearity enhancement, one should note that due to the large Q-factor of the pump and Stokes cavity modes the field is significantly enhanced as compared to the amplitude of the incident wave, which in conjunction with the fact that the Raman intensity is proportional to the local field to the power of 6, leads to the extremely large resonant enhancement of the Raman response of the silicon PhC metasurface. Moreover, the cavity field enhancement decreases when the angle of incidence increases, due to a weaker coupling between the incoming wave and the cavity modes, which results in reduced nonlinearity enhancement at larger θ.
As a final remark, it should be noted that the specific design of our metasurface ensures a particularly efficient Raman amplification. To be more specific, let us compare the spectral width of Raman interaction in silicon, ∆ν R = 105 GHz 59 , with the spectral width of the cavity mode at the Stokes frequency, ∆ν S = ω S /(2πQ S ) = 2.1 GHz. Thus, since ∆ν S ∆ν R , an efficient Raman interaction can be achieved.

V. CONCLUSION
In summary, two generic metasurfaces, a graphene metasurface based on graphene cruciform patches and a silicon metasurface with photonic crystal cavities as building blocks, are studied using a versatile and powerful homogenization method. In particular, in order to quantify the linear and nonlinear optical response of the two metasurfaces, we computed their effective permittivities and nonlinear susceptibilities. Our calculations revealed that, in both cases, the nonlinear optical response of the metasurfaces was enhanced by several orders of magnitude at the resonances of the metasurface building blocks. Moreover, by comparing the optical response of the metasurfaces and their homogenized counterparts, we showed that the homogenization approach is more suitable for graphene-based metasurfaces, because in this case the size of the resonant graphene nanostructures is much smaller than the operating wavelength. Even though the homogenization approach for silicon-based metasurfaces appeared to be less accurate, it could still provide valuable qualitative insights into their nonlinear optical response.
It should be noted that our homogenization approach is rather general, in that it can be readily extended to metasurfaces of different configurations or made of optical media with various dispersive and nonlinear optical properties. Moreover, the ideas presented in this paper have broad applicability, as they can be easily extended to other nonlinear optical interactions of practical interest, including third-harmonic generation, four-wave mixing, and sum-and difference frequency generation.