## Abstract

Intriguing and unusual physical properties of graphene offer remarkable potential for advanced, photonics-related technological applications, particularly in the area of nonlinear optics at the deep-subwavelength scale. In this study, we use a recently developed numerical method to illustrate an efficient mechanism that can lead to orders of magnitude enhancement of the third-harmonic generation in graphene diffraction gratings. In particular, we demonstrate that by taking advantage of the geometry dependence of the resonance wavelength of localized surface-plasmon polaritons of graphene ribbons and discs one can engineer the spectral response of graphene gratings so that strong plasmonic resonances exist at both the fundamental frequency and third-harmonic (TH). As a result of this double-resonant mechanism for optical near-field enhancement, the intensity of the TH can be increased by more than six orders of magnitude.

This article is part of the themed issue ‘New horizons for nanophotonics’.

## 1. Introduction

Graphene, a monolayer of carbon atoms arranged in a hexagonal lattice [1–3], has attracted over the last years a remarkable amount of interest, both at the fundamental science level and in the fields of engineering and technology, due chiefly to its outstanding and unique physical properties, as well as its potential for novel technological applications. In particular, this truly two-dimensional physical system possesses unprecedented electronic and thermal properties [4–7], which can be traced to its charge carriers having zero effective mass (the so-called Dirac fermions) and extremely large mean free path. This leads to very large thermal conductivity and high carrier mobility [8] and therefore provides an excellent materials platform for management of thermal processes in nanoelectronic circuits and high-frequency electronics [9–11]. Equally important, the physical properties of graphene can be readily and ultrafast tuned via chemical doping or electric gating, which makes it particularly appealing for technological applications in the area of ultrafast and active optoelectronic devices.

One field where the tunability of graphene could play a key role is photonics. This is primarily due to the fact that the dielectric constant of graphene is strongly dependent on the carriers density via the Fermi energy, which means that by chemical doping or electric gating one can change the dielectric constant of graphene, in relative terms, by much more than it can be done in the case of most other optical materials. Because generally the dielectric constant (index of refraction) is the main physical quantity that determines the optical properties of a material, the functionality of photonic devices that incorporate graphene in their structure can be actively and significantly modified. This feature of graphene has been recently exploited in a series of active photonic devices, including optical modulators, diffractive elements, optical switches, optical limiters, as well as photovoltaic and photoresistive devices [12–17].

Not only the linear optical properties of graphene can be exploited in photonic devices with new and advanced functionality, but also its nonlinear optical properties can play an equally important role. The main underlying reason why optical nonlinearities of graphene could play a key role both at the fundamental science level and in future graphene-based photonics technologies is the metallic character of graphene. To be more specific, owing to its metallic nature, graphene sheets or structured graphene can support propagating or localized surface-plasmon polaritons (SPPs), namely collective oscillations of the free carriers in graphene. The excitation of these highly localized SP modes is accompanied by a strong increase of the surface optical field and, owing to the nonlinear dependence of the intensity of the nonlinearly generated higher harmonics, an even stronger enhancement of nonlinear optical interactions. Importantly, the graphene tunability we alluded to and the strong dependence of the plasmon resonances on the geometrical configuration of graphene structures can be employed to spectrally shift over a large range the location of plasmon resonances, and thus one can greatly facilitate the design of new configurations of nonlinear interactions as well as active optical devices based on such nonlinear interactions. For example, enhanced second-order nonlinear optical effects have been predicted in hexagonal arrays of graphene patches [18], whereas increased third-order nonlinear optical interactions can be observed if graphene is incorporated in a photonic crystal [19].

In this paper, we demonstrate a mechanism for enhancing the third-harmonic generation (THG) in graphene-based optical diffraction gratings. Thus, we prove that by properly engineering the spectral location of plasmon resonances a pair of such resonances exist at the fundamental frequency (FF) and third harmonic (TH). This double-resonant excitation of localized graphene plasmons leads to strong enhancement of the local optical field, which in conjunction with the strong overlap between the optical fields of the plasmon modes, results in the increase of the intensity of the TH by several orders of magnitude in both one-dimensional and two-dimensional gratings. A brief discussion of the numerical method used in our study is also provided. Finally, we discuss the practical implications and potential technological applications of the ideas introduced in this paper, as well as other photonic structures based on graphene or other two-dimensional materials where our ideas could find applications.

## 2. Physical system and computational approach

The computational approach employed in this paper [20] can be used for a variety of physical settings, namely periodic photonic structures containing both bulk materials and graphene. Specifically, our numerical method incorporates the main linear and nonlinear optical effects in three-dimensional media and graphene, thus accurately describing the physics of such photonic structures. In particular, because the method describes the optical nonlinearities via nonlinear currents associated with various nonlinear polarizations in bulk materials and graphene, our computational approach can be used to investigate a multitude of nonlinear optical interactions, including second-harmonic generation, sum- and difference-frequency generation, THG, Kerr effects and four-wave mixing [21–24].

To be more specific, our numerical method for the calculation of the linear and nonlinear optical fields in the grating and far-field regions is an extension of the well-known rigorous coupled-wave analysis method [25–27], which incorporates the nonlinear optical processes via the corresponding nonlinear currents and the optical effects of the infinitely thin graphene sheets via modified (nonlinear) boundary condition at those interfaces where graphene layers are located. In addition, we used a recently introduced -matrix algorithm [20] for propagation of nonlinear fields generated by graphene sheets contained in multilayered periodic structures. In a very brief description, the numerical method consists of the following three steps: in the first step (*linear calculations*), one computes the field at the FF, *ω*_{0}. Then, in the second step (*polarization evaluation*), one evaluates the nonlinear polarization, **P**^{(nl,Ω)}(**E**^{(Ω)};**r**), at the TH, *Ω*. Finally, in the third step (*nonlinear calculation*), one calculates the generated nonlinear electric field, **E**^{(Ω)}.

The geometrical structure of the graphene gratings studied in this paper is schematically presented in figure 1. Thus, we considered one-dimensional gratings with period *Λ*, grating vector *K*=2*π*/*Λ*, and width of the graphene ribbons, *W*, and two-dimensional gratings with periods *Λ*_{x} and *Λ*_{y}, corresponding to grating vectors **K**_{x} and **K**_{y} lying in the (*x*,*y*) plane, the radius of the graphene discs being, *r*. In both cases, the graphene structures are assumed to be located at *z*=*z*_{s} and placed on a substrate with relative permittivity, *ϵ*_{s} (for specificity, assumed to be glass, *ϵ*_{s}=2.25), the relative permittivity of the cover being *ϵ*_{c} (for simplicity, assumed to be air, *ϵ*_{c}=1). The linear optical properties of graphene are characterized by its surface conductivity distribution, *σ*_{s}(*x*,*y*,*z*_{s}), so that we can describe the spatial distribution of the relative permittivity of the graphene grating by a single, piecewise continuous function, *ϵ*_{r}(*x*,*y*,*z*).

The graphene grating is excited by an incident harmonic plane-wave with the electric field expressed as
2.1where **r**_{t} is the in-plane component of the position vector, and are the in-plane and normal components of the wavevector of the incident wave in the cover, respectively, with being the wavenumber in the cover region, *ω* is the optical frequency and *c* is the speed of light. The polarization state of the incoming plane-wave is described by the relative field components in the transverse electric (TE) and transverse magnetic (TM) configurations, namely **E**_{TE} and **E**_{TM}, respectively, such that the field amplitude, , is defined by the polarization angle, *ψ*. In this study, *θ*=0 and *ψ*=*π*/2, so that the incident electric field is TM-polarized.

Because graphene is a single atomic layer physical system, one can conveniently describe its optical properties via surface quantities. In particular, the interaction of light with graphene can be readily described by the sheet conductance (sometimes also called conductivity), which in the random-phase approximation and zero temperature conditions is expressed as [28,29]
2.2where is the universal dynamic conductivity of graphene, *e* is the elementary charge, *ε*_{F} is the Fermi level, is the reduced Planck constant, *τ* is the relaxation time and *θ* is the Heaviside step function. The values of the Fermi level and relaxation time used in our investigations are *ε*_{F}=0.6 eV and *τ*=0.25 ps/2*π*, respectively. The dispersion of *σ*_{s}(λ) for graphene is presented in figure 2*a*. The sheet conductance contains the intraband (Drude) part, described by the first term in the r.h.s. of equation (2.2), and the interband component given by the next two terms. A more complete model valid at finite temperature can be used as well (see references [28–30]); however, the model described by equation (2.2) already describes the main features of graphene conductivity, namely the possibility of tuning *σ*_{s}(λ) by changing the Fermi level via chemical doping or by applying a gate voltage.

For computational purposes, it is usually more convenient to work with bulk rather than to surface quantities and therefore it is of practical interest to introduce bulk equivalents of the surface quantities. To this end, instead of the sheet conductance, *σ*_{s}, one uses a bulk conductivity, *σ*_{b}=*σ*_{s}/*h*_{eff}, where *h*_{eff} is the effective thickness of graphene. This approach can be confusing because of the ambiguity contained in the definition of the thickness of an atomic monolayer. Moreover, the optical properties of graphene can be equally well described by the electric permittivity, *ϵ*, which is related to the optical conductivity *σ*_{b} via the following equation:
2.3

The graphene lattice belongs to the space group, meaning that it is a centrosymmetric material. Therefore, quadratically nonlinear interactions such as second-harmonic generation are forbidden; however, THG is allowed and particularly strong in graphene [31,32]. The nonlinear optical response of graphene can be described by using the nonlinear optical conductivity tensor, *σ*^{(3)}_{s}, which relates the nonlinear surface current density, **j**^{nl}, at the TH and the electric field, **E**, at the FF. In this description, assuming that the graphene sheet lies in the (*x*,*y*) plane at *z*=0, the nonlinear current density, **J**^{nl}, can be written as [32]
2.4where **r**_{t} is the position vector lying in the graphene plane. Then, the nonlinear surface current density can be written as
2.5where *α*=*x*,*y*,*z* and **E**_{t} is the electric field component lying in the plane of graphene. Equation (2.5) shows that the nonlinear current density lies in the plane of graphene and only depends on the tangential field components. As a result, and . A formula for *σ*^{(3)}_{s}, derived under the assumptions that electron–electron and electron–phonon scattering as well as thermal effects can be neglected, has been recently derived in reference [32] in the following form:
2.6In this equation, *T*(*x*)=17*G*(*x*)−64*G*(2*x*)+45*G*(3*x*), with , and is the Fermi velocity, with *a*_{0}=1.42 Å being the separation distance between nearest-neighbour carbon atoms in graphene and *γ*_{0}=2.7 eV is the nearest-neighbour coupling constant. Figure 2*b* depicts *σ*^{(3)}_{s}(*ω*), whereas the inset of the figure highlights the lowest spectral peak at λ=1.033 μm.

As in the linear case, we introduce a ‘bulk’ nonlinear conductivity, . This nonlinear conductivity is particularly useful in experimental investigations of nonlinear optics of graphene because it is related to an effective bulk third-order nonlinear susceptibility, , the physical quantity that is usually measured in experiments. Using the fact that for harmonic fields **J**^{nl}(**r**,*t*)=−i*ω***P**^{nl}(**r**,*t*), where **P**^{nl}(**r**,*t*) is the nonlinear polarization, one can prove that
2.7where *Ω*_{t}=3*ω*_{0} is the frequency at the TH, with *ω*_{0} being the FF.

## 3. Results and discussion

To begin our analysis of THG in graphene gratings, let us consider the one-dimensional periodic array of graphene ribbons depicted in figure 1*a*. The period of the grating is *Λ*=100 nm and the width *W*=50 nm, that is a grating with duty cycle of 50%. The incident light is normally impinging onto this binary graphene grating and is TM polarized. For this grating, we have calculated using the method described in the previous section its linear spectral response, quantified by the reflectance, *R*, transmittance, *T* and absorption, *A*, as well as the spectrum of the TH intensity. The results of these calculations are presented in figure 3. In this and all our subsequent calculations, the intensity of the incoming beam was assumed to be , the TH intensity being measured relative to this reference value.

The spectra presented in figure 3 show a series of important features. Thus, it can be seen that both the optical absorption and TH intensity spectra have a series of spectral resonances, which in both cases are located at the same resonance wavelengths. The intensity of the corresponding spectral peaks decreases as the wavelength decreases, both in the case of absorption and TH intensity. In addition, the spectral width of the resonances, and, consequently, the associated optical losses, also decreases as the resonance wavelength decreases. Because absorption and TH intensity are directly related to the local optical field at the FF, the origin of these spectral resonances can be traced to the enhancement of this local optical field. Moreover, because graphene can be viewed as an infinitely thin metallic sheet, we can conclude that these resonances and the associated enhancement of the local field at the FF are owing to the excitation of localized SPs on the graphene ribbons.

In order to further validate these conclusions, we have determined the spatial field distribution of the local optical field at both the FF and TH, the results of these calculations being summarized in figure 4. These plots unambiguously demonstrate the excitation of localized graphene SPs at the resonance wavelengths, as the main features of such surface optical modes are clearly illustrated by these field profiles. Thus, it can be seen in figure 4 that the linear and nonlinear optical fields are strongly confined at the surface of the graphene ribbons, the degree of confinement increasing as the resonance wavelength decreases. Moreover, it can be seen that as the order of plasmon mode (resonance) increases, the number of local maxima of the field profile increases, which is one of the characteristics of a family of optical modes supported by a confined optical structure. Additionally, if we compare the field profiles at the FF (top panels in figure 4) and the TH (bottom panels in figure 4), we can see that they have very similar features. The reason for this is that the nonlinear optical response is primarily determined by the distribution of the linear optical near-field, as described by equation (2.5).

These conclusions are not specific to one-dimensional graphene diffraction gratings but they also remain valid, with certain changes, in the case of two-dimensional diffraction gratings. In order to illustrate this idea, we considered a two-dimensional periodic array of graphene discs, as depicted in figure 1*b*. The periods of the diffraction grating in this case are *Λ*_{x}=*Λ*_{y}=*Λ*=200 nm, the radius of the discs being *r*=50 nm. As in the one-dimensional case, the incident light is normally impinging onto the graphene grating. In order to characterize the spectral properties of this graphene structure, we have calculated the linear spectral response of the diffraction grating, namely the reflectance, transmittance and absorption, as well as the spectrum of the TH intensity, the main conclusions of these computations being outlined in figure 5. The results of these calculations show that, as in the one-dimensional case, both the optical absorption spectrum and the spectrum of the TH intensity exhibit a series of resonances, which are associated with the excitation of localized SPs on the graphene discs. In addition, the corresponding optical field is strongly confined and enhanced at the surface of the graphene discs, which explains the large increase in both optical absorption and nonlinear optical interaction at the surface of the graphene discs when graphene SPs are excited.

The optical spectra presented in figures 3 and 5 suggest there is a convenient way to design diffraction gratings based on graphene for enhancement of the THG and other nonlinear optical interactions in such diffraction gratings. Thus, let us assume that by varying the width of the graphene ribbons (or the diameter of the graphene discs) we can engineer the spectral resonances of the grating in such a way that the resonance wavelength of the fundamental plasmon coincides with the wavelength of the incoming beam, whereas the wavelength at the TH coincides with one of the resonance wavelengths of the higher-order plasmons. In other words, the resonance wavelength of one of these higher-order plasmon modes should be exactly one-third of the resonance wavelength of the fundamental plasmon mode. Under these circumstances, the diffraction grating will be efficiently excited at the FF, which will lead to a strong field enhancement at this frequency, and will radiate effectively at the TH, because there is also a plasmon resonance at this wavelength. In effect, such a diffraction grating would act as a highly effective receiver at the FF and a strong emitter (efficient antenna) at the TH.

In order to see whether it is possible to design such a configuration of plasmon resonances, we first determined the dispersion map of the absorption of the one-dimensional diffraction grating, that is the dependence of the absorption spectra on the width of the graphene ribbons, *W*. The results of these calculations, presented in figure 6, suggest that it is indeed possible to engineer a diffraction grating with the desired property. Thus, in addition to the plasmon absorption bands that in the linear spectra appear as spectral peaks, it can be seen in figure 6 that there are certain values of the graphene ribbon width, *W*, for which a graphene plasmon mode exists at both the FF and TH wavelengths. To be more specific, for *W*=85 nm, the graphene diffraction grating supports a (fundamental) plasmon mode at the FF corresponding to λ_{FF}=9.03 μ*m* and a third-order plasmon mode at λ_{TH}=λ_{FF}/3=3.01 μm. Another double-resonance occurs for *W*=97 nm in that case the plasmons excited at the FF and TH are a fundamental plasmon and a second-order one, respectively. We note that, because of the field periodicity, in figures 4 and 5 we have only shown the near-field distribution corresponding to the unit cell.

The next step of our analysis was devoted to the investigation of the nonlinear optical response of the one-dimensional graphene diffraction grating, so as to verify that enhanced THG is achieved when the double-resonance condition is satisfied. To this end, we have computed the dispersion map of the TH intensity of the diffraction grating, namely the dependence of the TH intensity spectra on *W*. This dispersion map of the nonlinear optical response of the diffraction grating is presented in figure 7, together with the dependence of the TH intensity on the width of the graphene ribbons, determined for the resonance wavelengths of the fundamental plasmon, which is shown in the inset of this figure. An important result illustrated in figure 7 is that, as previously discussed, the excitation of graphene-localized plasmons at the FF induces a strong increase of the intensity at the TH via local field enhancement, as it is clear from the location of the spectral resonance bands in figure 7. However, more importantly, the plot shown in the inset of figure 7 proves that a further enhancement of the TH intensity occurs when the double-resonance condition holds. To be more specific, the plot in the inset of figure 7 was determined by choosing the wavelength of the incident beam to be equal to the resonance wavelength of the fundamental plasmon (the plasmon band corresponding to the strongest enhancement of the TH) and varying the width of the graphene ribbons. It can be seen that a maximum intensity of the TH is achieved when *W*=85 nm, that is exactly at the width at which there is a plasmon mode at both the FF and TH. Importantly, from the results presented in figure 6, we find that the double resonance condition is fulfilled when the graphene ribbons are in relative close proximity to their nearest neighbours. Therefore, despite the fact that the excitation of localized surface plasmons on the graphene ribbons plays the major role in the observed enhancement of the TH intensity, the optical coupling between neighbouring ribbons and other diffractive effects could also affect the optical response of the graphene structure.

The conclusions valid for a periodic distribution of graphene ribbons also hold in the two-dimensional case. To illustrate this, we calculated the linear and nonlinear optical response of a two-dimensional array of square graphene patches with side length, *W*. This particular choice of the shape of graphene patches was guided by the fact that numerical simulations converge much faster in the case of square patches when compared to discs. From the results presented in figure 8, it can also be seen that in the two-dimensional case there are certain values of *W* for which a graphene plasmon mode exists at the FF and TH wavelengths. To be more specific, for *W*=84 nm, the graphene diffraction grating supports a (fundamental) plasmon mode at the FF corresponding to λ_{FF}=9.69 μm and a third-order plasmon mode at λ_{TH}=λ_{FF}/3=3.23 μm.

To verify that enhanced THG is achieved when the double-resonance condition is satisfied in the two-dimensional case, we have computed the dependence of the TH intensity spectra on the side length, *W*, of the square patches. This dispersion map of the nonlinear optical response is presented in figure 8. In the inset of figure 8, we also show the dependence of the TH intensity on *W*, determined for the resonance wavelengths of the fundamental plasmon together. Similar to the nonlinear spectra computed in the one-dimensional case, the spectral resonance bands in figure 8 show that the intensity at the TH is enhanced upon the excitation of graphene-localized plasmons at the FF. In addition, the plot in the inset of figure 8 shows that a maximum intensity of the TH is achieved when *W*=84 nm, that is exactly at the patch size at which there is a plasmon mode at the FF and TH. Thus, our results prove that the TH intensity can also be resonantly enhanced in two-dimensional graphene diffraction grating.

## 4. Conclusion

Here, we have investigated TH generation in one- and two-dimensional graphene diffraction gratings, paying special attention to the relations between the excitation of localized SPs on the graphene structures and the nonlinear optical response of the diffraction gratings. Our analysis has revealed that, similar to the case of plasmonic structures, the generation of localized SPs in graphene ribbons or discs leads to a strong enhancement of the local optical field. This, in turn, is accompanied by a large increase of the intensity of both the nonlinear near-field and the nonlinear radiation emitted into the far-field. We have also demonstrated that the intensity of the TH can be further enhanced by engineering the graphene diffraction gratings so as plasmon modes exist both at the frequency of the excitation field as well as at the frequency of the TH.

It should be noted that this dual-resonance approach to enhance the linear and nonlinear optical response of graphene structures is not specific to the gratings studied here but can also be extended to other configurations. For example, one expects that diffraction gratings containing graphene ribbons or discs with different sizes can generate an even stronger nonlinear optical response, as in this case, one can design the graphene structures in such a way that a fundamental plasmon exists at the fundamental frequency and TH (these plasmons correspond to the two graphene structures of different sizes). Finally, the ideas presented in this paper can be extended to other two-dimensional materials, such as transition-metal dichalcogenide monolayers, but the corresponding analysis is left for a future study.

## Authors' contributions

J.W.Y. and J.Y. performed the numerical simulations. M.W. developed the computer code used to perform the numerical simulations. N.C.P. conceived and designed the study, and drafted the manuscript. All authors analysed the data, read the manuscript and approved it.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by the European Research Council, grant agreement no. ERC-2014-CoG-648328 and the Royal Society’s International Exchanges Scheme. The work done by M.W. was partly supported through a UCL Impact Award graduate studentship supported by UCL and Photon Design, Ltd, whereas the work done by J.Y. was supported by a China Scholarship Council Studentship.

## Acknowledgements

M. Weismann and N. C. Panoiu wish to acknowledge the hospitality of Y. S. Kivshar and the Nonlinear Physics Centre of the Australian National University. The authors acknowledge the use of the UCL Legion High Performance Computing Facility (Legion@UCL) and associated support services in the completion of this work.

## Footnotes

One contribution of 15 to a theme issue ‘New horizons for nanophotonics’.

- Accepted November 23, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.