## Abstract

Optical second harmonic generation (SHG) from nanostructured graphene has been studied in the framework of classical electromagnetism using a surface integral equation method. Single disks and dimers are considered, demonstrating that the nonlinear conversion is enhanced when a localized surface plasmon resonance is excited at either the fundamental or second harmonic frequency. The proposed approach, beyond the electric dipole approximation used in the quantum description, reveals that SHG from graphene nanostructures with centrosymmetric shapes is possible when retardation effects and the excitation of high plasmonic modes at the second harmonic frequency are taken into account. Several SHG effects similar to those arising in metallic nanostructures, such as the silencing of the nonlinear emission and the design of double resonant nanostructures, are also reported. Finally, it is shown that the SHG from graphene disk dimers is very sensitive to a relative vertical displacement of the disks, opening new possibilities for the design of nonlinear plasmonic nanorulers.

© 2017 Optical Society of America

## 1. Introduction

Second-harmonic generation (SHG) is the nonlinear optical process corresponding to the conversion of two photons of an incident beam into a new photon at twice the fundamental frequency [1, 2]. The conversion rate for nonlinear optical phenomena is usually very low and a high intensity incident beam is consequently required to achieve a substantial nonlinear conversion. Although such a high intensity can be obtained with pulsed laser systems, researchers continue to seek new materials supporting strong nonlinear effects at low pump intensity [3–5], that would increase efficiencies of the devices and enable applications employing nonlinear optical phenomena, such as nanorulers [6], optoelectronics [7, 8], biological sensing [9] and signal processing [10].

Over the last decades, plasmonics has become an interesting field to design nonlinear nanostructures, giving birth to the field of nonlinear plasmonics [3]. Localized surface plasmons resonances (LSPR) correspond to the collective oscillations of the conduction electrons taking place in metallic nanoparticles [11]. LSPR dramatically increase the near-field intensity as they confine the electromagnetic field into sub-wavelength regions [12]. Furthermore, dimers of these nanoparticles enhance even more the field intensity between them, due to the coupling between the LSPR supported by each nanoparticle [13–15]. The field enhancement induced by LSPR has been widely used to enhance the nonlinear response of metallic nanostructures [16–21]. In addition, due to the small size of the nanoparticles compared to the wavelength, the paradigm of phase matching is transformed into mode matching [16,19]. Indeed, the LSPR frequency depends on the size, the shape and the composition of the nanoparticles but also on the surrounding medium and the potential coupling between plasmonic modes of non touching nanoparticles [12,22]. As a consequence, the resonance frequency of LSPR can be adjusted to match either the fundamental or the second harmonic frequency. When both the fundamental and the second harmonic frequencies correspond to LSPRs, the nonlinear response is significantly enhanced by this doubly resonant configuration [16,19,23]. However, the plasmon resonance of noble metal nanoparticles is hardly controllable after fabrication and it is then difficult to tune it in order to dynamically control their nonlinear response [8].

Recently, graphene has emerged as an interesting material for new developments in linear [8,24–38] and nonlinear plasmonics [39–50]. The electromagnetic response of this 2D material is tunable by applying a gate voltage, resulting in a modification of the plasmon frequency in the near infrared range [8, 24–26]. Moreover, graphene can sustain plasmons with relatively long lifetimes [25,26] and has a highly nonlinear response due to its linear dispersion relation [39,43]. In particular, nanostructured graphene, obtained by etching extended graphene with electron or focused ion beam [27–29], sustains LSPR with a level of confinement unattainable using metallic nanoparticles. The linear response of such nanostructured graphene has been widely studied both classically and quantum-mechanically [27–33, 36] while the nonlinear response of structures like nanotriangles and nanodisks of a few dozens nanometers has been mainly evaluated using quantum-mechanical calculation [40–44]. Studying larger nanostructures with this approach becomes challenging numerically due to the large number of involved electrons and a quantum description is not necessary as the quantum effects tends to vanish as the nanostructure size increases. Indeed, it was recently shown that the quantum effects in the nonlinear response of graphene nanostructures are negligible for nanostructures with typical dimensions larger than 40 nm [43]. In this study, J.D. Cox *et al* compared quantum simulations with classical ones in the electric dipole approximation. However, full-wave calculations of the nonlinear response are still lacking in the literature, despite their importance for large nanostructures and their demonstrated suitability in the linear regime.

SHG from graphene is often neglected due to the centrosymmetry of the crystal [41–43]. This argument is only valid in the electric dipolar approximation and for local dynamical conductivity. At local surface plasmon resonance, high local field and gradient lead to high multipolar response, as predicted for graphene [51]. The effect of the non-local conductivity of graphene on its second order optical response has been predicted [44]. Another source of SHG response of graphene is related to its very peculiar electronic structure and to the symmetry breaking associated with the valley polarization [52,53]. SHG response is also expected for nanostructured systems for which the symmetry is locally broken at the edge.

In this article, the linear and the second harmonic (SH) responses of a single graphene nanodisk with a 100 *nm* diameter, as well as nanodimers composed of the same two graphene nanodisks, are investigated using a surface integral equations formalism (SIE). For each nanostructure, the linear response is considered first, in order to confirm the suitability of the SIE for computing the optical response of nanostructured graphene and then the SHG is discussed in details. The main goal of this study is to investigate the SHG from nanostructured graphene beyond the electric dipole approximation, although a quantum description is not conducted in this case, since it was emphasized that full-wave computations are in general required to accurately describe the SHG from centrosymmetric nanostructures [54].

## 2. Numerical method

#### 2.1. Surface integral equations methods for the linear and second harmonic responses

The linear optical response of graphene nanostructures is computed using a SIE formulation considering a plane wave excitation. The typical triangle mesh used here for graphene is 30 nm^{2}, corresponding to approximately 1200 carbon atoms, and the thickness of the graphene particle is set to 0.5 nm. It was shown in the linear regime that this thickness is small enough to mimic a single graphene layer [25]. The dielectric constants for graphene at the fundamental and SH frequencies are evaluated as described in the next section. The surface integrations have been performed with an optimized approach in order to improve the numerical accuracy [55]. This improvement was mandatory due to the very short distance between triangles, a direct consequence of the ultra-thin thickness of the meshes with respect to the average triangle size. The convergence of the computations has been carefully assessed.

For SHG computations, the linear surface currents, which are expanded on Rao-Wilton-Glisson (RWG) basis functions [56], are used for the computation of the fundamental electric field inside the graphene layer in order to evaluate the surface SH polarization. This approach does not take into account the edge terminations. However, it was shown that classical computations agree better with the optical response of graphene nanostructures with arm-chair edges than with zig-zag edges [36]. For the SHG, a SH polarization vector is associated to each mesh triangle. The nonlinear polarization is computed considering the component ${\chi}_{nnn}^{(2)}$ of the surface tensor, where *n* denotes the normal to the surface mesh triangles. In this case, the SH sources standing on the top and bottom surfaces of the graphene nanodisks cancel out in the electric dipole approximation, as expected for a centrosymmetric material. Furthermore, the triangles on the disk edges reproduce the SHG induced by the centrosymmetry breaking at the edges. Although values for the tensor elements of the second order nonlinear susceptibility of graphene are only sparsely reported so far, recent studies seems to confirm the validity of this approximation in the case of graphene [50,57,58]. Furthermore, the implementation of all the tensor elements in the SIE formalism, which is straightforward, is not mandatory to assess the multipolar nature of the second harmonic wave. Moreover, the frequency dependence for ${\chi}_{nnn}^{(2)}$ which is set to 1 in the present work, is not considered here, allowing to obtain qualitative results [54].

The combination of a quantum evaluation of the nonlinear susceptibility with a full-wave evaluation of the radiation properties is probably the best method for a complete description of the SHG from nanostructured graphene, although this complete description is beyond the scope of this article [59,60]. For example, it will be interesting to investigate the influence of both chiral edge currents [61] and valley polarization [52,53] on the SHG from nanostructured graphene, although the comparison between the SHG induced by graphene edges and valley polarization is beyond the scope of the present study. The SH surface currents are obtained by solving the SIE formulation taking into account the nonlinear polarization and enforcing the boundary conditions at the nanostructure surfaces [62]. Similarly to the linear surface currents, the SH surface currents are expanded on RWG basis functions. The expanding coefficients are found by applying the method of moments with Galerkin’s testing [56]. A Poggio-Miller-Chang-Harrington-Wu-Tai formulation is used to ensure accurate solutions even at resonant conditions [56]. The SH electric field is then deduced from the SH surface currents using a two-term subtraction method for the evaluation of the Green’s functions [63].

The eigenmodes are found by seeking for the eigenvectors of the SIE matrix associated to vanishing eigenvalues without any external excitation [23]. This is practically done by scanning the complex frequency plane and computing the eigenvalues of the SIE matrix. When the complex frequency corresponds to an eigenfrequency, the eigenvalue vanishes and the corresponding eigenvector is then associated to an eigenmode.

#### 2.2. Dielectric function of graphene

The permittivity of graphene has been calculated within the random-phase approximation (RPA). In this model, the conductivity of graphene is given by the Kubo formula [25,64,65]:

*E*= 0.4 eV, a relaxation time

_{F}*τ*such that

*ħτ*

^{−1}= 1.6 meV and a temperature of

*T*= 300 K. The electron charge is

*e*,

*ħ*is the reduced planck constant,

*ω*is the frequency and

*k*is the Boltzman constant. The first term in Eq. (1) dominates for an energy much smaller than 2

_{B}*E*and can be reduced to a Drude-like conductivity proportional to the square root of the charge carrier density

_{F}*n*, when

*E*≫

_{F}*k*. As interband transitions become non-negligible for energy close to 2

_{b}T*E*, one has to consider this term since the calculation of the SHG requires the permittivity at twice the frequency, which can be close to the Fermi energy.

_{F}Furthermore, a third term has been added to the conductivity to take into account the contribution of the edge electron states, resulting from the finite nature of the nanostructures. For nanodisks, this term has been calculated as [66]:

*β*is the

_{ln}*n*zero of the Bessel function of order

^{th}*l*,

*l*= 100 and

_{max}*n*= 20 are convergence parameters,

_{max}*ω*=

_{R}*v*with

_{F}/R*R*the disk radius,

*η*= 1/

*τ*and

*f*

_{ln,0}is the difference between the Fermi function at zero energy and at energy equal to

*ħω*. When the disk radius is larger than 50 nm, this term becomes quite small. For energies lower than

_{R}β_{ln}*E*/2, the relative difference between the conductivity of extended graphene and of 50 nm radius graphene disks is less than a percent and this contribution is negligible. However, the conductivity of extended graphene becomes smaller for energies closer to

_{F}*E*and the relative difference increases. Thus, this term has been taken into account in this work, although it does produce only a negligible modification of the computed spectra, with a maximum relative redshift of the order of one percent. The resulting conductivity is then

_{F}*σ*(

*ω*) =

*σ*(

_{k}*ω*) +

*σ*(

_{E}*ω*) and the relative permittivity of graphene is [25]: where

*t*is the thickness of graphene, which has been fixed at 0.5 nm in the present work. This value assures that computation converges properly.

*ε*is the relative permittivity of the surrounding medium, hence

_{b}*ε*= 1 in vacuum. In this article, the spectra are evaluated for energies ranging from 0.05 to 0.32 eV. Both real and imaginary parts of graphene permittivity are plotted in Fig. 1(a), top panel, for this energy range.

_{b}## 3. Results and discussion

#### 3.1. Single graphene disk

For all computations reported in this article the nanodisk diameter is fixed to *d* = 100 nm. The incident plane wave is propagating perpendicularly to the disk. In both cases (linear regime and SHG) the scattered intensity is integrated over a sphere of radius *R* = 50 *μ*m. The scattering cross-section from a single graphene nanodisk is considered first (Fig. 1(a), middle panel). The scattering spectrum reveals two peaks, one at 0.153 eV and one at 0.307 eV. It was demonstrated that, if only the Drude-like conductivity of graphene is considered (first term of Eq. (1) when *E _{F}* ≫

*k*), the resonance frequency

_{b}T*ω*of the dipolar mode for a single nanodisk of suspended graphene (without substrate), depends on its diameter and the Fermi energy and is given by [28,29]:

_{d}Contrary to noble metals, the properties of graphene give the possibility to tune the resonance frequency of the dipolar mode by applying a gate voltage, which modifies the Fermi energy. Though a full study of the influence of the Fermi energy is not conducted here, we note that this effect allows the tuning of the linear and nonlinear response of nanoparticles after fabrication. For a 100 nm diameter disk and a Fermi energy of 0.4 eV, Eq. (5) gives an resonant energy at 0.151 eV. The maximum of the first peak in the scattered field spectrum occurs at 0.153 eV, indicating that this resonance corresponds to the dipolar mode. This observation is confirmed by the eigenmode analysis (Fig. 1(b)) in excellent agreement with Eq. (5). The small difference is attributed to the modification of the conductivity caused by the interband transition and the electronic edge states, which are not considered in Eq. (5). A second resonance is observed at 0.307 eV (see eigenmode on Fig. 1(b)). A multipolar analysis reveals that this mode corresponds to an octupole (data not shown).

Now, we turn our attention to the SHG from the same graphene nanodisk (Fig. 1(a), bottom panel). In this case, the spectrum reveals three peaks. The two highest peaks in the SH spectrum correspond to a resonant excitation of the graphene nanostructures at the pump wavelength at respectively 0.153 and 0.308 eV. Indeed, these peaks are at the same energies as the peaks observed in the scattering spectrum and correspond to the resonant excitation of the dipolar and octupolar modes, respectively. However, a third peak is observed at 0.102 eV. The eigenmode analysis shows that a quadrupolar mode exists at 0.203 eV. This mode does not contribute to the linear scattering spectrum due to its symmetry and the normal incidence of the wave but explains the peak at low energy observed in the SH spectrum: an incident beam at 0.102 eV is converted into a SH wave at 0.204 eV, the resonant energy of the quadrupolar mode.

The observation of a quadrupolar emission is characteristic of the SHG from centrosymmetric nano-objects [67, 68]. Indeed, the SHG is forbidden in centrosymmetric media only in the electric dipole approximation, and the description of the SH wave must include high order modes. Furthermore, in the case of the graphene nanostructures studied here, the thickness is much smaller than the incident wavelength, meaning that the retardation effects are negligible at the pump excitation stage. In this case, the nonlinear sources distributed over the top and bottom mesh surfaces indeed cancel out each other and do not contribute to the far-field radiation. At the same time, this limitation does not stand for the SH sources located at the graphene edges, allowing for the generation of SH light. Due to its parity, the spatial distribution of the SH sources cannot induce a dipolar moment along the direction of polarization, but can indeed excite quadrupolar and other SH modes with even parity. This explains why there is no peak in the SHG spectrum at 0.153/2 *eV*.

#### 3.2. Dimers of graphene nanodisks

In this section, dimers composed of graphene nanodisks with gap distances *g* and vertical shifts *h* are considered, Fig. 2(a). The *z*-axis is defined perpendicular to the graphene plane and the *x*-axis and *y*-axis are defined respectively parallel and perpendicular to the dimer axis when there is no vertical shift. Then *g* is measured along *x* and *h* along *z*. The incident wave is propagating along the *z*-axis and is polarized along the *x*-axis. The case of the dimer with different gap distances but without vertical shift (*h* = 0 nm) is considered first (Fig. 3 and 4).

Three peaks are observed in the scattering spectra (Fig. 3(a), top panel). The two peaks at low energy blueshift as the gap distance increases. On the contrary, the spectral position of the highest energy peak is only slightly influenced by the gap distance. Here, the eigenmode analysis reveals two bonding and two anti-bonding modes. Bonding modes and anti-bonding modes correspond to modes with the charges of the two discs in the gap having respectively the opposite and the same signs. The dipolar and quadrupolar bonding modes correspond to the two modes at low energy observed in the spectra, while the two anti-bonding modes are not observed (Fig. 3(b)). Interestingly, the scattering intensity does not decrease significantly with the gap distance for the dipolar mode (inset Fig. 3(a)).

Nonlinear optical processes require high field intensity to be substantial. It is well known that due to the coupling between the LSPR, the near field intensity is enhanced in the nanogap. The spectra for the field intensity enhancement (the ratio between the field intensity at the considered point and the incident field intensity) evaluated at the center of the nanogap for different gap distances (Fig. 3(a), middle panel), clearly indicate that there are other modes between 0.2 eV and 0.3 eV, not visible in the far-field linear response.

These modes are identified by plotting the imaginary part of the charge distributions at the peak frequency (Fig. 4(a)). Indeed, when a resonant system is excited at its eigenmode frequency, its response is 90° out of phase relatively to the excitation [69]. In the frequency domain, this phase relation translates to complex valued quantities (field and charges oscillations for example) where the argument of the complex number gives the phase difference between the excitation and the considered physical quantity. The imaginary part thus gives the resonant part of the system response [69].The first two modes at low energy have already been identified as the bonding dipolar and bonding quadrupolar modes. The bonding dipolar mode permits to enhance the intensity in the nanogap by up to seven orders of magnitude, in agreement with computations done with a finite element method [36]. The other modes also have a bonding nature, but between higher order modes. Their bonding nature results in the adequate symmetry for these modes being excited by an incident planewave at normal incidence. It is worth mentioning that an anti-bonding nature results in the adequate symmetry for SH emission. Interestingly, for the last peak at 0.305 eV, the real part of the charge distribution appears to be different than the imaginary part, (Fig. 4(a)) indicating that two eigenmodes spectrally overlap at this energy.

For all the observed modes, the near-field enhancement increases as the distance between the graphene nanodisks decreases. However, the SHG spectra indicate that SHG induced by the fundamental dipolar mode increases slowly with the gap distance (Fig. 3(a), bottom panel). This behavior is explained by the silencing of the SHG emitted from the gap due to the symmetry of the nonlinear sources. The same effect has been recently reported in the case of dimers composed of gold nanorods [70]. Indeed, the SH sources standing at each side of the nanogap oscillate out of phase, as evidenced by the longitudinal components (along the *x*-axis) of the SH near-field distribution (Fig. 4(b)). As a consequence, the SH waves coming from each side of the nanogap tend to annihilate each other in the far-field, suppressing a high contribution to far-field SH intensity. This silencing effect increases as the distance between the nanodisks, i.e. between the out-of-phase SH sources, decreases. This tends to compensate the fundamental intensity enhancement occurring in the gap, resulting in an attenuation of the SHG. To conclude this analysis of the SHG from graphene dimer with different gap distances, it is worth noting the steep increase in the SH intensity from the graphene dimer with a gap of 80 nm. It occurs when the bonding dipolar mode is resonantly excited at the fundamental frequency and is due to the resonant excitation of a second mode at the SH frequency, around 0.3 eV resulting in a doubly resonant configuration [16,19,23].

#### 3.3. Influence of a vertical shift on the SHG from graphene dimers

To overcome the limitation in the SHG yield due to the silencing effect, it is necessary to reduce the destructive interference between the SH field coming from the nonlinear sources at each side of the nanogap. One simple way to achieve this is to vertically shift, i.e. along the propagation direction *z* of the incident wave, one of the nanodisks relatively to the other one (see Fig. 2).

The spectra for the scattering, near-field enhancement and SHG are shown for vertical shifts *h* ranging from 0 nm to 40 nm (Fig. 5(a)) with a constant horizontal gap distance of 2 nm. As the vertical shift increases, the effect on the scattering and the enhancement of the fundamental electric field between both nanodisks is similar to the effect of an increase of the distance between the two nanodisks (Fig. 3(a)): for the bonding dipolar mode observed at low energy, a blueshift is observed as the vertical shift increases and the intensity increases slightly. However, the field enhancement in the nanogap decreases quickly when the vertical shift increases.

On the other hand, the SH intensity in the far-field first increases with the vertical shift, as expected from the symmetry breaking, but then decreases because the coupling between the two graphene nanodisks vanishes. The maximum occurs for a vertical shift close to 20 nm. In Fig. 5(b), the near-field intensity at the fundamental and the second harmonic frequencies are plotted for a graphene dimer without vertical shift (left panels) and for a dimer with a vertical shift *h* = 20 nm (right panels).

The near-field maps do not show a large influence of the vertical shift, even though there is a slight symmetry breaking in the SH near-field intensity, not present in the fundamental intensity distribution. The effect of the vertical shift is obvious when looking at the emission pattern of the dimer (Fig. 5(b), bottom panel). Without vertical shift, the dimer does not emit SH waves along the *x*-axis but only along the *y*-direction. The symmetry breaking induced by the vertical shift clearly enables the SHG along the dimer axis with an intensity approximately 6.2 times higher than the intensity along the y-axis, which is not enhanced by the vertical shift.

The maximum of SH intensity is plotted as a function of the vertical shift in Fig. 5(c). Although a small vertical shift of 1 or 2 nm already dramatically enhances the SH intensity in the near-field (on Fig. 5(c), left panel, the intensity is taken 1 nm away from the edge of a disk), the symmetry breaking is not high enough to stop the destructive interference occurring in the far-field (Fig. 5(c) right panel). A vertical shift as high as 20 nm is required to maximize the SHG in the far-field. Indeed, the symmetry breaking along with a strong enhancement of the near-field intensity allows an improvement in the yield of the SHG. For a gap distance of 2 nm, the vertical shift permits to enhance the SHG up to a factor 10. This phenomenon can be used to determine the vertical shift between two adjacent graphene nanodisks, with a resolution close to 1 nm. This high resolution is directly related to the limited thickness of the graphene layer.

## 4. Conclusion

In conclusion, full wave calculations of the SH response of graphene disks have been performed using a surface integral equation method. These results confirm that the nonlinear conversion is enhanced when a LSPR is excited at either the fundamental or second harmonic frequency. The role played by the different plasmonic modes has been determined, in connection with their symmetry properties, using an eigenmode analysis. In the case of dimers, silencing of the SHG has been observed and the influence of the symmetry breaking induced by a vertical shift of one of the disks has been investigated in details. Especially, it was pointed out that the far-field SH intensity reaches a maximal value for a vertical shift of 20 nm. This optimal geometry results from the competition between fundamental intensity enhancement and symmetry breaking. It can be noted that a non-normal incidence would have some effects. Firstly, the coupling between the pump wave and the graphene nanostructures would change as well as the field enhancement would decreases. Secondly, the symmetry would be broken and the selection rules for the SHG would be different. Other methods for breaking the symmetry in nanostructured graphene have been reported in the literature. For example, it is possible to control the charge doping in each disk [32], and it is anticipated – on the basis of the results reported in this article – that this method is promising for the dynamic control of the nonlinear optical conversion down to the nanoscale [71]. Furthermore, it will be interesting to investigate the influence of both chiral edge currents [61] and valley polarization [52,53] on the SHG from nanostructured graphene.

## Funding

European Research Council (ERC-2015-AdG-695206 Nanofactory); Swiss National Science Foundation (project 200020_153662).

## Acknowledgments

This work was completed while Dr Michaël Lobet was a recipient of a Fellowship of the Belgian American Educational Foundation.

## References and links

**1. **R. W. Boyd, *Nonlinear Optics* (Academic, 2008).

**2. **Y.R. Shen, *The Principles of Nonlinear Optics* (Wiley-Blackwell, 2002).

**3. **M. Kauranen and A. V. Zayats, “Nonlinear plasmonics,” Nat. Photon. **3**, 737–748 (2012). [CrossRef]

**4. **J. B. Khurgin and G. Sun, “Plasmonic enhancement of the third order nonlinear optical phenomena: Figures of merit,” Opt. Express **21**(22), 27460–27480 (2013). [CrossRef] [PubMed]

**5. **J. Butet, P.-F. Brevet, and O. J. F. Martin, “Optical Second Harmonic Generation in Plasmonic Nanostructures: From Fundamental Principles to Advanced Applications,” ACS Nano , **9**(11), 10545–10562 (2015) [CrossRef] [PubMed]

**6. **J. Butet and O. J. F. Martin, “Nonlinear plasmonic nanorulers,” ACS Nano **8**(5), 4931–4939 (2014). [CrossRef] [PubMed]

**7. **F. Bonaccorso, Z. Sun, T. Hasan, and A. C. Ferrari, “Graphene Photonics and Optoelectronics,” Nat. Photon. **3**, 611–622 (2010). [CrossRef]

**8. **Q. Bao and K. P. Loh, “Graphene Photonics, Plasmonics, and Broadband Optoelectronic Devices,” ACS Nano , **6**(5), 3677–3694 (2012). [CrossRef] [PubMed]

**9. **P. Pantazis, J. Maloney, D. Wu, S. E. Fraser, P. Pantazis, J. Maloney, D. Wu, and S. E. Fraser, “Second harmonic generating (SHG) nanoprobes for in vivo imaging,” PNAS **107**(33), 14535–14540 (2010). [CrossRef] [PubMed]

**10. **E. Garmire, “Nonlinear optics in daily life,” Opt. Express **21**(25), 30532–44 (2013). [CrossRef]

**11. **S. A. Maier, *Plasmonics: Fundamentals and Applications* (Springer2007).

**12. **M. I. Stockman, “Nanoplasmonics: The physics behind the applications,” Phys. Today **64**(2), 39–44. (2011). [CrossRef]

**13. **P. Muhlschlegel, H.-J. Eisler, O. J. F. Martin, B. Hecht, and D. W. Pohl, “Resonant Optical Antennas,” Science **308**(5728), 1607–1609 (2005). [CrossRef] [PubMed]

**14. **I. Romero, J. Aizpurua, G. W. Bryant, and F. J. G. De Abajo, “Plasmons in nearly touching metallic nanoparticles: singular response in the limit of touching dimers,” Opt. Express **14**(21), 9988–9999 (2006). [CrossRef] [PubMed]

**15. **P. Nordlander, C. Oubre, E. Prodan, K. Li, and M. I. Stockman, “Plasmon Hybridization in Nanoparticle Dimers,” Nano Lett. **4**(5), 899–903 (2004). [CrossRef]

**16. **K. Thyagarajan, S. Rivier, A. Lovera, and O. J. F. Martin, “Enhanced second-harmonic generation from double resonant plasmonic antennae,” Opt. Express **20**(12), 12860–12865 (2012). [CrossRef] [PubMed]

**17. **Z. Fang and X. Zhu, “Plasmonics in nanostructures,” Adv. Mater. **25**(28), 3840–3856 (2013). [CrossRef] [PubMed]

**18. **J. Butet and O. J. F. Martin, “Fano resonances in the nonlinear optical response of coupled plasmonic nanostructures,” Opt. Express **22**(24), 29693 (2014). [CrossRef]

**19. **M. Celebrano, X. Wu, M. Baselli, S. Gromann, P. Biagioni, A. Locatelli, C. De Angelis, G. Cerullo, R. Osellame, B. Hecht, L. Duo, F. Ciccacci, and M. Finazzi, “Mode matching in multiresonant plasmonic nanoantennas for enhanced second harmonic generation,” Nature Nanotech. **10**(5), 412–417 (2015). [CrossRef]

**20. **B. Metzger, M. Hentschel, T. Schumacher, M. Lippitz, X. Ye, C. B. Murray, B. Knabe, K. Buse, and H. Giessen, “Doubling the Efficiency of Third Harmonic Generation by Positioning ITO Nanocrystals into the Hot-Spot of Plasmonic Gap-Antennas,” Nano Lett. **14**(5), 2867–2872 (2014). [CrossRef] [PubMed]

**21. **H. Aouani, M. Rahmani, M. Navarro-Cia, and S. A. Maier, “Third-harmonic-upconversion enhancement from a single semiconductor nanoparticle coupled to a plasmonic antenna,” Nature Nanotech. **9**(4), 290–294 (2014). [CrossRef]

**22. **N. J. Halas, S. Lal, W. S. Chang, S. Link, and P. Nordlander, “Plasmons in strongly coupled metallic nanostructures,” Chem. Rev. **111**(6), 3913–3961 (2011). [CrossRef] [PubMed]

**23. **G. D. Bernasconi, J. Butet, and O. J. F. Martin, “Mode analysis of second-harmonic generation in plasmonic nanostructures,” J. Opt. Soc. Am. B **33**(4), 768 (2016). [CrossRef]

**24. **E. H. Hwang and S. D. Sarma, “Dielectric function, screening, and plasmons in two-dimensional graphene,” Phys. Rev. B **75**(20), 205418 (2007). [CrossRef]

**25. **F. H. L. Koppens, D. E. Chang, and F. J. G. de Abajo, “Graphene Plasmonics: A Platform for Strong Light-Matter Interactions,” Nano Lett. **11**(8), 3370–3377 (2011). [CrossRef] [PubMed]

**26. **F. J. G. de Abajo, “Graphene Plasmonics: Challenges and Opportunities,” ACS Photonics , **1**(3), 135–152 (2014). [CrossRef]

**27. **H. Yan, F. Xia, Z. Li, and P. Avouris, “Plasmonics of coupled graphene micro-structures,” New J. Phys. **14**(12), 125001 (2012). [CrossRef]

**28. **Z. Fang, S. Thongrattanasiri, A. Schlather, Z. Liu, L. Ma, Y. Wang, P. M. Ajayan, P. Nordlander, N. J. Halas, and F. J. G. de Abajo, “Gated Tunability and Hybridization of Localized Plasmons in Nanostructured Graphene,” ACS Nano , **14**(3), 2388–2395 (2013). [CrossRef]

**29. **Z. Fang, Y. Wang, A. E. Schlather, Z. Liu, P. M. Ajayan, F. J. G. De Abajo, P. Nordlander, X. Zhu, and N. J. Halas, “Active tunable absorption enhancement with graphene nanodisk arrays,” Nano Lett. **14**(1), 299–304 (2014). [CrossRef]

**30. **H. Yan, X. Li, B. Chandra, G. Tulevski, Y. Wu, M. Freitag, W. Zhu, P. Avouris, and F. Xia, “Tunable infrared plasmonic devices using graphene/insulator stacks,” Nature Nanotech. , **7**(5), 330–334 (2012). [CrossRef]

**31. **F. Ramirez, B. Liu, and S. Shen, “Extreme blueshift of surface plasmon resonance frequency in graphene nanodisk stacks,” J. Quant. Spectrosc. Radiat. Transf. **158**, 27–35 (2015). [CrossRef]

**32. **G. Rosolen and B. Maes, “Asymmetric and connected graphene dimers for a tunable plasmonic response,” Phys. Rev. B , **92**(20), 205405 (2015). [CrossRef]

**33. **S. Thongrattanasiri, A. Manjavacas, and F. J. G. De Abajo, “Quantum finite-size effects in graphene plasmons,” ACS Nano , **6**(2), 1766–1775 (2012) [CrossRef] [PubMed]

**34. **W. Wang, P. Apell, and J. Kinaret, “Edge plasmons in graphene nanostructures,” Phys. Rev. B **84**(8), 085423 (2011). [CrossRef]

**35. **A. Manjavacas, P. Nordlander, and F. J. G. de Abajo, “Plasmon Blockade in Nanostructured Graphene,” ACS Nano **6**(2), 1724–1731 (2012). [CrossRef] [PubMed]

**36. **S. Thongrattanasiri and F. J. G. de Abajo, “Optical field enhancement by strong plasmon interaction in graphene nanostructures,” Phys. Rev. Lett. **110**(18), 187401 (2013). [CrossRef] [PubMed]

**37. **J. Christensen, A. Manjavacas, S. Thongrattanasiri, F. H. L. Koppens, and F. J. G. De Abajo, “Graphene plasmon waveguiding and hybridization in individual and paired nanoribbons,” ACS Nano **6**(1), 431–440 (2012) [CrossRef]

**38. **Z. Wang, T. Li, K. Almdal, N. Asger Mortensen, S. Xiao, and S. Ndoni, “Experimental demonstration of graphene plasmons working close to the near-infrared window,” Opt. Lett. **41**(22), 5345 (2016). [CrossRef] [PubMed]

**39. **S. A. Mikhailov, “Non-linear electromagnetic response of graphene,” EPL **79**(2), 27002 (2007). [CrossRef]

**40. **J. D. Cox, M. R. Singh, M. a. Anton, and F. Carreño, “Plasmonic control of nonlinear two-photon absorption in graphene nanocomposites,” J. Phys. Condens. Matter **25**(38), 385302 (2013). [CrossRef] [PubMed]

**41. **J. D. Cox and F. J. G. de Abajo, “Electrically tunable nonlinear plasmonics in graphene nanoislands,” Nature Comm. **5**, 5725 (2014). [CrossRef]

**42. **J. D. Cox and F. J. G. de Abajo, “Plasmon-Enhanced Nonlinear Wave Mixing in Nanostructured Graphene,” ACS Photonics , **2**(2), 306–312 (2015). [CrossRef]

**43. **J. D. Cox, I. Silveiro, and F. J. G. de Abajo, “Quantum Effects in the Nonlinear Response of Graphene Plasmons,” ACS Nano , **10**(2), 1995–2003 (2016). [CrossRef] [PubMed]

**44. **M. T. Manzoni, I. Silveiro, F. J. G. de Abajo, and D. E. Chang, “Second-order quantum nonlinear optical processes in single graphene nanostructures and arrays,” New J. Phys. **17**(8), 83031 (2015). [CrossRef]

**45. **E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, and S. A. Mikhailov, “Coherent nonlinear optical response of graphene,” Phys. Rev. Lett. **105**(9), 1–4. (2010). [CrossRef]

**46. **S. A. Mikhailov, “Theory of the giant plasmon-enhanced second-harmonic generation in graphene and semiconductor two-dimensional electron systems,” Phys. Rev. B **84**(4), 045432 (2011). [CrossRef]

**47. **M. Gullans, D. E. Chang, F. H. L. Koppens, F. J. G. de Abajo, and M. D. Lukin, “Single-photon nonlinear optics with graphene plasmons,” Phys. Rev. Lett. **111**(24), 247401 (2013). [CrossRef]

**48. **M. Jablan and D. E. Chang, “Multiplasmon Absorption in Graphene,” Phys. Rev. Lett. **114**(23), 236801 (2015). [CrossRef] [PubMed]

**49. **D. A. Smirnova, I. V. Shadrivov, A. E. Miroshnichenko, A. I. Smirnov, and Y. S. Kivshar, “Second-harmonic generation by a graphene nanoparticle,” Phys. Rev. B **90**(3), 035412 (2014). [CrossRef]

**50. **M. Lobet, M. Sarrazin, F. Cecchet, N. Reckinger, A. Vlad, J. F. Colomer, and D. Lis, “Probing Graphene Chi(2) Using a Gold Photon Sieve,” Nano Lett. **16**(1), 48–54 (2016). [CrossRef]

**51. **J. D. Cox, A. Marini, and F. J. G. de Abajo, “Plasmon-assisted high-harmonic generation in graphene,” Nature Comm. **8**, 14380 (2017). [CrossRef]

**52. **T. O. Wehling, A. Huber, A. I. Lichtenstein, and M. I. Katsnelson, “Probing of valley polarization in graphene via optical second-harmonic generation,” Phys. Rev. B **91**(4), 041404 (2015). [CrossRef]

**53. **L. E. Golub and S. A. Tarasenko, “Valley polarization induced second harmonic generation in graphene,” Phys. Rev. B **90**(20), 201402 (2014). [CrossRef]

**54. **C. Forestiere, A. Capretti, and G. Miano, “Surface integral method for second harmonic generation in metal nanoparticles including both local-surface and nonlocal-bulk sources,” J. Opt. Soc. Am. B **30**(9), 2355 (2013). [CrossRef]

**55. **T. V. Raziman, W. R. C. Somerville, O. J. F. Martin, and E. C. Le Ru, “Accuracy of surface integral equation matrix elements in plasmonic calculations,” J. Opt. Soc. Am. B **32**(3), 485(2015). [CrossRef]

**56. **A. M. Kern and O. J. F. Martin, “Surface integral formulation for 3D simulations of plasmonic and high permittivity nanostructures,” J. Opt. Soc. Am. A **26**(4), 732 (2009). [CrossRef]

**57. **V. a. Margulis, E. E. Muryumin, and E. a. Gaiduk, “Optical second-harmonic generation from two-dimensional hexagonal crystals with broken space inversion symmetry,” J. Phys. Condens. Matter **25**(19), 195302 (2013). [CrossRef] [PubMed]

**58. **Y. Q. An, J. E. Rowe, D. B. Dougherty, J. U. Lee, and A. C. Diebold, “Optical second-harmonic generation induced by electric current in graphene on Si and SiC substrates,” Phys. Rev. B **89**(11), 115310 (2014). [CrossRef]

**59. **J. D. Cox, R. Yu, and F. J. G. de Abajo, “Analytical description of the nonlinear plasmonic response in nanographene,” Phys. Rev. B **96**(4), 045442 (2017). [CrossRef]

**60. **J. L. Cheng, N. Vermeulen, and J. E. Sipe, “Second order optical nonlinearity of graphene due to electric quadrupole and magnetic dipole effects,” Sci. Rep. **7**, 43843 (2017). [CrossRef] [PubMed]

**61. **J. Karch, C. Drexler, P. Olbrich, M. Fehrenbacher, M. Hirmer, M. M. Glazov, S. A. Tarasenko, E. L. Ivchenko, B. Birkner, J. Eroms, D. Weiss, R. Yakimova, S. Lara-Avila, S. Kubatkin, M. Ostler, T. Seyller, and S. D. Ganichev, “Terahertz Radiation Driven Chiral Edge Currents in Graphene,” Phys. Rev. Lett. **107**(27), 276601 (2011). [CrossRef]

**62. **T. Heinz, *Nonlinear Surface Electromagnetic Phenomena* (Ponath and Stegeman, 1991)

**63. **J. Butet, B. Gallinet, K. Thyagarajan, and O. J. F. Martin, “Second-harmonic generation from periodic arrays of arbitrary shape plasmonic nanostructures: a surface integral approach,” J. Opt. Soc. Am. B **30**(11), 2970 (2013). [CrossRef]

**64. **L. A. Falkovsky and A. A. Varlamov, “Space-time dispersion of graphene conductivity,” EPJB , **56**(4), 281–284 (2007). [CrossRef]

**65. **L. A. Falkovsky, “Optical properties of graphene,” J. Phys.: Conf. Ser. **129**(11), 12004 (2008).

**66. **T. Christensen, W. Wang, A.-P. Jauho, M. Wubs, and N. A. Mortensen, “Classical and quantum plasmonics in graphene nanodisks: Role of edge states,” Phys. Rev. B **90**(24), 241414 (2014). [CrossRef]

**67. **J. I. Dadap, J. Shan, K. B. Eisenthal, and T. F. Heinz, “Second-Harmonic Rayleigh Scattering from a Sphere of Centrosymmetric Material,” Phys. Rev. Lett. **83**(20), 4045–4048 (1999). [CrossRef]

**68. **J. I. Dadap, J. Shan, and T. F. Heinz, “Theory of optical second-harmonic generation from a sphere of centrosymmetric material: small-particle limit,” J. Opt. Soc. Am. B **21**(7), 1328–1347 (2004). [CrossRef]

**69. **T. V. Raziman and O. J. F. Martin, “Does the real part contain all the physical information,” J. Opt. **18**(9), 95002 (2016). [CrossRef]

**70. **J. Berthelot, G. Bachelier, M. Song, P. Rai, G. Colas des Francs, A. Dereux, and A. Bouhelier, “Silencing and enhancement of second-harmonic generation in optical gap antennas,” Opt. Express **20**(10), 10498 (2012). [CrossRef] [PubMed]

**71. **Y. Q. An, F. Nelson, J. U. Lee, and A. C. Diebold, “Enhanced Optical Second-Harmonic Generation from the Current-Biased Graphene/SiO 2 /Si(001) Structure,” Nano Lett. **13**(5), 2104–2109 (2013). [CrossRef] [PubMed]