Nonlinear optics in plasmonic nanostructures

Although a relatively new area of nanoscience, nonlinear plasmonics has become a fertile ground for the development and testing of new ideas pertaining to light–matter interaction under intense field conditions, ideas that have found a multitude of applications in surface science, active photonic nanodevices, near-field optical microscopy, and nonlinear integrated photonics. In this review, we survey the latest developments in nonlinear plasmonics in three-dimensional (metallic) and two-dimensional (graphene) nanostructures and offer an outlook on future developments in this field of research. In particular, we discuss the main theoretical concepts, experimental methods, and computational tools that are used together in modern nonlinear plasmonics to explore in an integrated manner nonlinear optical properties of metallic and graphene based nanostructures.


Introductory remarks
Nonlinear optics at the nanoscale is becoming one of the major success stories of modern optics. Its beginnings can perhaps be traced to the advent of the laser [1], which has made possible the first demonstration of a nonlinear optical process, namely the second-harmonic generation (SHG) [2]. Ever since, we have witnessed a growing research interest in understanding the nonlinear optical interactions between light and matter. The scientific interest in nonlinear optical phenomena has been driven not only by the challenges they present at the fundamental science level, but also by the plethora of applications they have spawn in the commercial arena, including wavelength converters for telecommunications, optical signal processing, all-optical switching, ultrafast optical modulators, optical microscopy, surface science, optical tomography, and biosensing applications.
A nonlinear optical process can be viewed as the physical interaction between one or multiple optical beams and an optical medium, in which the response of the medium depends nonlinearly on the electric field of the beams. Much of the research focus during the early stages of nonlinear optics has been on nonlinear optical phenomena pertaining to the propagation of interacting near-monochromatic optical beams in bulk nonlinear optical media. To achieve strong nonlinear optical effects in this setting, a key prerequisite is to ensure that the interacting beams are phase-matched. In particular, the nonlinear medium plays a secondary rôle in determining the efficiency of the nonlinear optical processes, as it only facilitates the nonlinear coupling between the optical beams. As a result, the development of efficient phasematching techniques has been a central theme early on in nonlinear optics. Some of the most used techniques are birefringent phase matching [3], quasi-phase matching in periodic media [4,5], and modal phase matching in nonlinear optical waveguides [6,7].
In recent years, advances in nanofabrication techniques and a dramatic improvement of our ability to probe the optical field at deep-subwavelength scale have led to a gradual shift in the main research directions in nonlinear optics, from the study of coherent, phase-matched optical interactions that take place over many wavelengths to nonlinear optical processes in which the optical (near-) field interacts with matter over just a few wavelengths or even at subwavelength scale. In this emerging paradigm of nonlinear optics, the efficiency of the nonlinear optical processes is determined not by the quality of the phase matching between the interacting optical beams but by the degree of confinement and overlap between the optical near-field and nonlinear optical structures with subwavelength features. A representative example of such non-phase-matched nonlinear optical interaction is the nonlinear scattering of optical waves from dielectric or metallic particles, a process in which the excitation of Mie resonances [8,9] can lead to a significant enhancement of the amount of light generated at higher harmonics.
Extreme nonlinear optics, which occurs under the condition of very short laser pulse illumination and intense optical power, represent another research direction branching out from traditional nonlinear optics. In this regime, the influence of laser pulse carrier wave is dominant over that of laser pulse profile, and consequently many intriguing and unexpected nonlinear phenomena occur, including high-harmonic generation in gases, frequency doubling in inversion symmetric materials, and generation of attosecond electromagnetic pulses or pulse trains. Whereas this research field does not belong to the remit of this review, the interested reader is referred to several excellent references on this topic [10,11].
Since at the nanoscale the efficiency of nonlinear optical processes is primarily determined by the overlap between the optical structure and the optical near-field, significant theoretical and experimental research efforts in nonlinear nanooptics have been focused on designing optimized configurations in which not only that there is a good spatial matching between the optical near-field and the optical nanostructure but also the near-field is resonantly enhanced. For example, advances in optical near-field spectroscopy techniques make it possible not only to retrieve information about the intensity of the electric field but also to map full vectorial characteristics of both the electric and magnetic fields [12], both in the real-space [13] and the momentum-space (k-space) [14].
Importantly, these experimental developments have been facilitated by and spurred on research in computational electromagnetism. This has led to complex and highly accurate finite-difference [15,16] or finite-element [17] based numerical methods that can be used to describe the electromagnetic field both in the time and frequency domains, at deep-subwavelength and deep-suboptical-cycle scales.
One well-known and widely studied approach to significantly increase the optical near-field, and consequently the efficiency of nonlinear optical interactions at the nanoscale, is to employ metallic nanostructures, such as optical diffraction gratings, nano-apertures, and nanoparticles. This approach relies on the fact that metal-dielectric interfaces support surface plasmon polaritons (SPPs) [18][19][20][21][22], which are p-polarized strongly localized surface waves associated to oscillations of free electrons in metals. The generation of these localized or propagating modes is accompanied by strong spatial confinement and enhancement of the electric field at the interface between the two media. Because of this strong field enhancement and the large contrast between the permittivity of metals and that of the surrounding dielectric environment, most parameters defining the physical properties of SPP resonances, including the spatial distribution of the optical near-field and the resonance frequency, depend strongly on the geometry and optical constants of metallodielectric nanostructures. In particular, the tight confinement of the optical field of SPPs makes the surface effects at metaldielectric interfaces distinctly sensitive to variations in the shape of the plasmonic nanostructures and the dielectric characteristics of the embedding optical medium.
Despite the fact that SPPs are eminently surface phenomena, their generation in metals requires three-dimensional (3D) configurations as metallic nanostructures are 3D physical systems. The recent discovery of two-dimensional (2D) materials, however, has facilitated the study of SPPs and other nano-optics related phenomena in purely 2D physical systems. One such 2D material with salient physical properties is graphene, a monolayer of carbon atoms distributed in a hexagonal lattice [23,24]. This 2D nanomaterial has metallic characteristics and relatively low optical losses at terahertz and mid-IR frequencies, and thus it supports SPPs [25,26]. Equally important, due to the 2D nature of graphene, the optical field enhancement of graphene SPPs is much larger than that achievable in metallic systems. In conjunction with its unusually large nonlinear optical constants, this makes graphene an ideal material platform to study nonlinear optical interactions in 2D physical systems. More specifically, graphene, as a centrosymmetric material, exhibits large thirdharmonic generation (THG) [27,28], strong optical Kerr nonlinearity [29,30], and induced second-order nonlinearity [31][32][33] in a single atomic layer.

Scope and outline of this review
Although it is a relatively new nanoscience discipline, research in nonlinear plasmonics has produced a large body of results and thus has already been the topic of a few excellent review articles. Some notable examples are reviews of nonlinear metasurfaces [65,66], a review of second-and thirdorder nonlinear optical phenomena in plasmonic systems [67], and a review of nonlinear plasmonic antennas [68]. We considered, however, that a comprehensive review of nonlinear plasmonics, which would cover in-depth and in an integrated way the main experimental techniques, theoretical concepts, and computational methods that play the central rôle in modern nonlinear plasmonics is still missing. Equally important, we wanted to collect in the same place information about nonlinear plasmonic phenomena in the two most common physical systems where such effects are studied, namely metallic structures and graphene. Because of these characteristics, we hope that our review article will provide valuable knowledge both for experienced scientists who wants to keep abreast with the latest developments in nonlinear plasmonics as well as young researchers who ponder about entering this dynamic and exciting field of research.
Our review is organized as follows: first, we review in section 2 the theoretical concepts and tools used to investigate the nonlinear optical properties of plasmonic systems, with a particular focus on second-and third-order nonlinear optical processes. In section 3, we present experimental results in key areas of nonlinear plasmonics, including nonlinear SPPs at metal surfaces, plasmon-enhanced nonlinear response of single and arrays of metallic nanoparticles, phase-controlled nonlinear optics in plasmonic metasurfaces, mode-matching enhanced harmonics generation in plasmonic nanoantennas, and plasmon-enhanced nonlinear effects in nonlinear active materials. We then present in section 4 main results pertaining to nonlinear optical phenomena associated to the generation of SPPs in graphene. In section 5, we present the main numerical methods used to model nonlinear optical properties of plasmonic systems and illustrate how they can be used for specific nonlinear plasmonic nanostructures and active plasmonic devices. In particular, we discuss both time-and frequency-domain numerical methods. Future perspectives and the main conclusions of our review are summarized in section 6.

Nonlinear optical processes: theoretical background
In this section, we introduce the main theoretical concepts pertaining to nonlinear optics in plasmonic nanostructures. We start with a brief description of wave propagation in bulk nonlinear optical media, and continue with a more in-depth discussion of nonlinear light-matter interaction in the subwavelength regime, which is particularly relevant in the context of nonlinear plasmonics. We will focus on the two main types of nonlinear optical processes, namely secondand third-order nonlinear optical interactions.

Wave interactions in nonlinear optical media
The electromagnetic response of an optical medium is predominantly determined, at optical frequencies, by the polarization of the medium, t P r, ( ). In order to examine the nonlinear response of an optical medium, one expands in power series the polarization of the medium in terms of the electric field, is the ith-order bulk optical susceptibility, which is generally an (i+1)th-rank tensor, and ò 0 is the vacuum permittivity. Susceptibilities depend on the crystal structure of the medium, as in the case of graphene, whereas for homogeneous and isotopic materials, such as noble metals that play a central rôle in plasmonics, the susceptibilities are characterized by scalar quantities.
The polarization of the medium can be divided in a linear and nonlinear part, t t t P r P r P r , , , l n l = + ( ) ( ) ( ), according to its power dependence on the electric field, where r  P r  P r   ,  , :  ,  ,  , ,   c c In the equation above, the second-and third-order polarizations are, respectively: Let us now assume that the incident electromagnetic field is the superposition of two monochromatic plane waves. Under these circumstances, the incident electric field can be written as: where ω 1 and ω 2 are the frequencies of the interacting waves, k 1 and k 2 are the corresponding wave vectors, and 'cc' means complex conjugation.
Inserting ( ( ) The terms in the right hand side of this equation describe SHG at 2ω 1 , SHG at 2ω 2 , sum-frequency generation (SFG) at ω 1 +ω 2 , difference-frequency generation (DFG) at ω 1 −ω 2 , and optical rectification, respectively. Note that the last term in (6) describes a static polarization. The polarization t P r, 2 ( ) ( ) in (6) can be used to describe a multitude of important nonlinear optical processes in bulk optical media, including SHG, SFG, DFG, optical parametric amplification, and optical parametric oscillations. In all of these nonlinear interactions, the phase-matching between the interacting waves plays an important rôle. On the other hand, in nonlinear optics at the nanoscale, wave phase-matching between interacting optical fields does not affect the nonlinear optical interactions, and therefore in what follows the wave vector dependence of the optical fields is dropped.
If we express the second-order nonlinear polarization t P r, 2 ( ) ( ) as: the nonlinear polarizations in (6) describing different secondorder nonlinear optical processes can be written as: The most common second-order nonlinear optical processes of those described by the equations above, namely SHG, SFG, and DFG, are schematically illustrated in figure 1.
In practice, in bulk media it is not possible to observe simultaneously all nonlinear effects described by the nonlinear polarizations in the above equations because usually only one nonlinear process is phase-matched. In nonlinear nano-optics, however, it is no longer necessary to achieve phase-matching and therefore more frequencies can be generated as a result of nonlinear optical interactions.
The most common third-order nonlinear optical process investigated in nonlinear plasmonics and more generally nonlinear optics is the Kerr effect. It amounts to a change of the index of refraction of a nonlinear optical medium upon the interaction with an optical field. The third-order nonlinear polarization describing the Kerr effect can be expressed as: whereas the corresponding change of the index of refraction is: In these equations, 3 c ( ) is the third-order susceptibility and n 2 is the nonlinear-index coefficient. Sometimes one defines this coefficient in terms of the field intensity, I, via n n I I 2

D =
, the two nonlinear coefficients being related by n n n c I 2 2 0 0  = ( ), with n 0 being the index of refraction of the medium.

Quadratic nonlinear optical interactions in centrosymmetric media
Most metals used in plasmonics are centrosymmetric, namely their crystal lattice remains invariant upon inversion symmetry transformations. As far as quadratically nonlinear optical interactions are concerned, this property implies that the second-order susceptibility 2 c ( ) is identically equal to zero. This property can be easily demonstrated by considering that the field E is applied to a crystal with inversion symmetry. Under these conditions, the nonlinear polarization can be written as P E E : , which in conjunction with (4a) leads to 0 2 c º ( ) . Since the (local) dipole-allowed bulk nonlinear polarization cancels in a centrosymmetric optical medium, the quadratically nonlinear optical response in such nonlinear optical media is defined by higher-order effects. In particular, according to a model that is widely used in the study of the SHG in centrosymmetric media [69], the nonlinear polarization consists of two components. First, a (local) dipoleallowed surface nonlinear polarization, P r s 2 ( ) ( ) , defined within a surface layer several ngströms thin where the inversion symmetry property is broken. Second, a (nonlocal) bulk nonlinear polarization, P r b 2 ( ) ( ) , which is generated inside the nonlinear optical medium by electric quadrupoles and magnetic dipoles. In what follows, we briefly discuss the physical properties of these two nonlinear polarizations.
2.2.1. Surface contribution to the nonlinear optical response of centrosymmetric media. The surface contribution to the nonlinear optical response of centrosymmetric media is described by a nonlinear polarization, which can be written as: where γ, δ′, β, and ζ are material parameters. This (nonlocal) polarization originates from electric quadrupoles and magnetic dipoles located in the bulk of the medium [74]. The third term in (13) can usually be neglected, as in homogeneous media . Also, most theoretical models [75] predict that the second term is negligible, too. For example, in the case of plane wave propagation in a homogeneous medium this term exactly cancels due to the transverse character of the field at the (fundamental frequency) FF. Moreover, in the case of noble metals the ratio between δ′ and γ is of the order of n w [70], where ν is the damping frequency, and this ratio is negligible at optical frequencies. We note, however, that the degree to which the term proportional to δ′ influences the SHG is still a matter of debate [70]. Based on these considerations, in practice one sets δ′=0 and neglects the third term in (13). Furthermore, in the case of noble metals, the anisotropy parameter, ζ, has a negligible value, so that one usually sets ζ=0, too.

Experimental aspects of nonlinear plasmonics
By squeezing the incident electromagnetic radiation into a nanoscale near-field region, metallic nanostructures provide a prominent nanoplasmonic platform for boosting up various nonlinear optical effects. In addition, well-established nanofabrication techniques, including top-down approaches, such as electron-beam lithography and focused ion beam milling, and bottom-up methods, which include wet chemistry and colloidal self-assembly, allow for flexible design and preparation of plasmonic nanostructures with tailored nonlinear optical response. In the past two decades, a variety of nonlinear optical effects have been experimentally exploited with distinctive plasmonic nanostructures, including mainly the second-order nonlinear response such as second-harmonic generation, sum-/ difference-frequency generation, and third-order ones like third-harmonic generation, two-photon absorption induced luminescence (TPL), and four-wave mixing (FWM).
To clearly illustrate the role of plasmonic resonances in the nonlinear optical processes, as well as the opportunities they could offer for advanced nonlinear applications, in this section, we focus our attention mainly on the most representative and experimentally readily-performed ones-second-harmonic generation, third-harmonic generation, and FWM (see figure 2). The main text body will review many important and representative experimental studies on the nonlinear optical effects in various plasmonic nanostructures, and the content covers the following sub-topics: nonlinear SPPs at metal surfaces; localized surface-plasmon (LSP)enhanced nonlinear response of individual metallic nanostructures and arrays; phase-controlled nonlinear optics in plasmonic metasurfaces; mode-matching enhanced harmonics generation in plasmonic nanoantennas; and plasmonenhanced nonlinear effects in nonlinear active materials.

Nonlinear surface plasmon-polaritons at metal surfaces
Compared to conventional bulk nonlinear crystals, metals with considerably lower optical nonlinearities are seemingly not a suitable choice for nonlinear optics research and applications. In addition, the shallow penetration of electromagnetic waves in metals leads to the optical response of bulk metals mainly governed by the weak light-matter interaction in the region near the metal-dielectric interface. At the nanoscale, however, the light-matter interaction can be significantly enhanced in the presence of surface plasmons, which can leverage various nonlinear optical processes. Because of the centrosymmetric lattice structure of metals, the second-order nonlinear response of metals vanishes for the bulk and the dominant contribution originates from the lattice-constants-thick layer near the surface. In contrast, the third-order nonlinear response of metals is allowed in the bulk region and originates predominantly from the oscillating nonlinear dipoles in a surface layer with thickness of the order λ/2π, which is evidently much larger than the lattice constant. Importantly, excitation of SPPs can create surface waves traveling along the metal surface [76,77] and thus confine the free-space waves into the subwavelength-thick surface layer for enhanced nonlinear response.
Surface plasmon-polariton induced nonlinear polarization in metals was first reported for second harmonic generation in 1974, where by exciting the surface mode at a Agair interface enhanced the overall nonlinear polarization [78]. It has been found that the SPPs can modify the SHG response of the metal film in two different ways. First, the SPP excitation at the pump wavelength results in an enhanced electric field near the interface, which can significantly increase the nonlinear polarization in the interfacial region [79]. Second, the surface waves can be launched at the second-harmonic wavelength under the right phase-matching conditions [79]. This harmonic surface waves can then constructively interfere with other interfacial second-harmonic radiation, leading to an enhanced overall SH emission into the far-field [80].
Except for SHG, optical excitation of surface polaritons and SPPs at semiconductor and metal surfaces has also been demonstrated to affect the third-order nonlinear optical processes, such as FWM and THG. For instance, simultaneously exciting the surface polaritons in GaP at the pump frequencies ω 1 and ω 2 dramatically enhanced the nonlinear emission (c) A nano-patterned metal surface strongly scattering the SPP waves at the FWM frequency enables a nonlinear dark-field microscopy approach. Reprinted with permission from [85]. Copyright (2010) American Chemical Society. signal at the FMW frequency 2ω 1 − ω 2 [81]. In addition, launching surface polaritons at the FWM frequency can also be expected to generate strong nonlinear optical interaction if the wave-vector matching condition is fulfilled at the interface. Different from the SHG process, which only occurs at the metal surface, the surface polaritons mediated FWM emission can extend over the penetration depth of surface modes into the materials, which enables a sensitive spectroscopy probe for the third-order bulk nonlinear susceptibility, 3 c ( ) . In addition to the surface FWM waves, the intrinsic FWM originating from the metal bulk also contributes to the total emission, although the SPP-enhanced FWM can be dominant in proper experimental configurations [82]. Similar phenomena have also been observed for third-harmonic generation. More specifically, it has been demonstrated in recent experiments [80,83] that launching SPP modes at the FWM frequency was achievable by carefully selecting the incidence angles of the pump beams, as per figure 2(a). In this experimental configuration, the incidence waves were not directly coupled to the surface waves and thus enabled freespace excitation of SPP modes at the metal surface. In particular, the FWM signal can be further increased by coupling surface FWM waves with a nanostructured metal surface which is capable of more efficiently transferring the surface waves into outgoing, free-space propagating waves, as indicated in figure 2(b) [84]. Finally, the same group also demonstrated that a metal surface with locally structured features can dramatically enhance FWM, which has consequently inspired the development of nonlinear dark-field microscopes as depicted in figure 2(c) [85].

Localized surface plasmon-enhanced nonlinear response of individual metallic nanostructures and arrays
Distinct from the SPPs propagating at metal surfaces, localized surface plasmons supported by metallic nanostructures are non-propagating localized modes consisting primarily of evanescent waves. By efficiently confining the incident light into the nanoscale volume near the nanostructure surface, the LSP excitation creates substantially enhanced near-field intensities, which lays the ground for strongly leveraging various nonlinear effects both in plasmonic metals themselves and active nonlinear materials placed in their close vicinity.
As perhaps the simplest model to understand the nonlinear response from plasmonic nanostructures, a metallic nanosphere has been extensively investigated to reveal the origins of the nonlinearly generated signal, as illustrated in figure 3(a). According to the symmetry selection rule, the second-order nonlinear response originates from both the surface and bulk contributions [86]. Within the electrostatic approximation, the bulk nonlinear polarization, which dominates the total nonlinear emission intensity such as THG of individual plasmonic nanoparticles [87], has generally a dipolar character. The dipolar bulk second-order nonlinear response, however, is forbidden for most metals because of their structural inversion centrosymmetry. Thus, observable second-order nonlinear emission, like SHG in bulk metals, originates mainly from the nonlinear polarization at the metal surface where symmetry breaking occurs. Specifically, the total SH emission from a plasmonic nanoparticle can be viewed as the interference between all the nonlinear dipoles distributed on the particle surface.
For a perfectly spherical nanoparticle that is much smaller than the wavelength, the local nonlinear dipoles at the opposite sides of the sphere surface cancel out between each other, leading to negligible SH emission into the far-field. Indeed, the hyper-Rayleigh scattering results pertaining to spherical metal nanoparticles show a weak dipolar SH response for very small particles, as shown in figure 3(b), but a quadrupolar response for relatively large particles (see figures 3(c) and (d)) [88]. Close inspection reveals that the dipolar response for the small particles originates from geometry asymmetry associated with shape deviation from a perfect sphere [89,90], and that the retardation effect is responsible for the quadrupole SH emission observed in the larger particles [91,92]. For metal nanoparticles of other shapes and more complex geometries, the multipolar characters of their nonlinear response can be derived following the same approach as for the spherical nanoparticles and the signature of their mode characters can be more pronounced due to plasmonic amplification and modulation [93,94].
It is generally believed that the following two enhancement mechanisms are critical in typical nonlinear spectroscopy of plasmonic nanostructures. The first one is the plasmon induced excitation enhancement of the incident field at the fundamental frequency, which substantially strengthens the multi-waves coupling in the nonlinear conversion process. The second mechanism is related to the plasmon-amplified nonlinear emission efficiency, which is realized by coupling the outgoing nonlinear signals to radiative LSP modes or manipulating the local interference between the nonlinear sources at the subwavelength scale. In the first mechanism, complex metallic nanostructures having near-field plasmon coupling are favorable candidates for boosting the nonlinear effects (see figure 4). These coupled nanostructures having nanometer-sized gaps possess stronger field localization than individual nanoparticles and typically display a few orders of magnitude field enhancement within the gap region. Since the nonlinear polarization strength scales quadratically or cubically with the fundamental field intensity for second-and third-order nonlinear interactions, respectively, the plasmonic field intensity enhancement can result in at least tens of thousands times stronger nonlinear emission intensities [95,96].
When the near-field coupling across the nanoscale interparticle gap increases the nonlinear excitation rate, the diffractive long-range coupling occurring in metallic nanoparticle arrays can also modify the nonlinear response. As indicated in figure 5(a), two samples with subtle differences in ordering are found to have SH responses differing by a factor of up to 50 [97]. In these arrays, the linear optical properties are predominantly determined by the long-range diffractive coupling between the elemental particles. Since the nonlinear properties are strongly dependent on the linear optical properties, the orientational distribution and detailed ordering of the nanoparticles modify the plasmon resonance characteristics and thus affect the nonlinear excitation efficiency. Similar effects have also been demonstrated in periodically distributed nanodisk monomers and dimers [98]. The coexistence of localized plasmon modes and far-field photonic modes creates a Fano-type electromagnetic coupling in such arrays, which is responsible for the enhanced secondharmonic generation, as per figure 5 Another important representative manifestation of the geometry-ordering influence on the nonlinear response of plasmonic nanostructures is the presence of circular dichroism (CD) in second-harmonic generation (termed as SHG-CD) from structures with geometry chirality [99][100][101]. As depicted in figure 5(c), the chiral structure consisting of G-shaped gold elements exhibits a pronounced SHG-CD effect. The distinct SHG intensity patterns (see the four micrographic images rendered with false color) resolved by nonlinear microscopy with circularly polarized pump beams indicate that the CD effect originates from the enantiomerically sensitive plasmon modes. When all the G-shaped nano-elements in the two arrays are replaced with the same elements (the same geometry and orientation), the SHG-CD effect vanishes because of the identical plasmonic response of the reordered arrays.

Phase-controlled nonlinear optics in nonlinear plasmonic metasurfaces
Inspired by the concept of spin-orbit coupling of light in linear optics, we have recently witnessed growing research interest in nonlinear geometrical phase control, which opens up new ways for manipulating the harmonic generations at the nanoscale [102].
Upon excitation by a circularly polarized pump wave, the phase of the nonlinear polarization of a plasmonic nanostructure can be controlled by using the structure geometry and orientation-dependent spin-orbit coupling. Thus, for a circularly polarized pump beam propagating in the +z direction (see figure 6), the electric field can be described by  A single metal nanostructure itself or a nonlinear medium in the close vicinity of the metal nanoparticle locally forms an effective nonlinear dipole moment expressed by: where ω denotes the angular frequency of the pump wave, and a is the nth-order harmonic nonlinear polarizability tensor of the individual nanostructure with an orientation angle of θ.
Coordinate rotation can be employed to facilitate the analysis of orientation angle dependent nonlinear dipole moment of the nanostructure. As depicted in figure 6(a), the local coordinate frame is fixed to the nanostructure. When the local coordinate system (x′, y′) of the nanostructure is rotated by an angle of θ, the pump wave acquires a geometric phase through the spin-orbit coupling of the incident light, Here, the symbol 'L' denotes the local coordinate of the nanostructure.  In the local frame, the nth-order harmonic nonlinear polarizability is simply 0 0 a a = q q= | , and thus the corresponding nonlinear dipole moment can be written as: This nonlinear dipole moment can be decomposed into two in-plane rotating dipoles (according to the circular polarization states) as: After the coordinate transformation from the local frame to the laboratory frame, these two rotating dipole moments can be written in the following form: As a result, the nonlinear polarizabilities of the individual nanostructure can be expressed as: Obviously, a geometry phase equal to either n i is introduced in the nonlinear polarizability of the nth-order harmonic generation from the nanostructure with the same or opposite circular polarization to that of the pump wave. According to the selection rules for harmonic generation of circularly polarized pump waves, a single nanostructure with m-fold rotational symmetry allows only harmonic orders of n lm 1 =  , where l is an integer and the signs '+' and '−' correspond to harmonic generation of the same and opposite circular polarization, respectively [103,104]. For example, a single nanorod with two-fold rotational symmetry (symmetry group 2  ) allows for thirdharmonic generation with both the same and opposite circular polarization to that of the pump wave, as per figure 6(b), and the spin-dependent phases of the harmonic waves are 2σθand 4σθ, respectively. However, a nanostructure with four-fold rotational symmetry (symmetry group 4  ), like the nano-cross in figure 6(c), permits only third-harmonic signals with circular polarization state opposite to the incident polarization. According to equation (19b), the introduced geometry phase of this nonlinear wave is 4σθ.
Importantly, under circularly polarized illumination, the linear optical response of a single nanostructure with rotational symmetry is independent of its spatial orientation [102]. Thus, by properly assembling of rotation-symmetrical nanostructures with spatially varied orientation in a 3D or 2D lattice, one can form a nonlinear metamaterial or metasurface which shows homogeneous linear optical properties but well-defined local nonlinear polarizability distribution [102,103,105]. In particular, plasmonic nanostructure based nonlinear metasurfaces are the most suitable candidates for manipulating the local phase-controlled nonlinear response. This is because the resonant excitation of localized surface plasmons in the constituent nano-elements can substantially increase the pump efficiency that consequently leverages the generated harmonic signals from the nonlinear metasurface.
As indicated in figure 7(a), the measured diffractive TH intensity patterns for the orientation-varied nanorod ( 2  symmetry) array suggest that TH waves with the same circular polarization are mainly diffracted into the first-order grating angle while those with the opposite circular polarization are concentrated in the second-order diffraction angle. In comparison, the metasurface consisting of nano-crosses with 4  symmetry diffracts TH waves with opposite circular polarization into the first-order grating angle. Similarly, the efficient nonlinear responses for a 3  -type metasurface, such as the triangle-hole-arrayed 2D metamaterial shown in figure 7(b) [103], are limited to the second-harmonic generation with opposite circular polarization to that of the pump wave. Since this geometry-phase tunable nonlinear response is particularly sensitive to the pumping wavelength and polarization state, plasmonic metasurfaces with multiple-fold rotational symmetry can therefore be appealing candidates for applications in optical data storage, information encryption and background-free image reconstruction.
An illustrative such application is the image encoding technique realized with an ultrathin nonlinear metasurface, which has been recently demonstrated by controlling the phase difference of nonlinear polarization between two neighboring meta-atoms in a pixel, the details being presented in figure 7(c) [105]. Because of the plasmon-assisted nonlinear pumping efficiency and controlled local SH response of the 3  -nanostruture pair in a pixel, the image burned into the metasurface can only be read by a circularly polarized laser beam with wavelength at the plasmon resonance and then reconstructed through detection at the SH wavelength. For comprehensive reviews on nonlinear photonic metasurfaces, interested readers are referred to two excellent, recently published articles [65,66].

Mode-matching enhanced harmonics generation in single plasmonic nanoantennas
Nonlinear spectroscopy of single nanostructures is instrumental for fully understanding their nonlinear response [92,106]. In this context, it is more proper to treat a single plasmonic nanostructure as an optical antenna in the nonlinear regime, and the associated emission of harmonic waves can be viewed as a nonlinear scattering process. Note that many arrayed plasmonic nanostructures can be classified into this group, provided that the lattice distance between neighboring elements is properly designed such that they neither couple with each other in the near-field domain nor interfere in the far-field region [107,108].
The quest for efficient nonlinear light sources at the nanoscale in the past two decades has driven extensive efforts to enhance the nonlinear optical effects in plasmonic nanoantennas. It has been well known that plasmonic nanostructures possessing resonances are the preferred candidates for leveraging the nonlinear response by exploiting the huge near-field enhancement upon resonant optical excitation. In general, the plasmon resonance of an individual nanostructure can be tuned to spectrally overlap with either the excitation wavelength at the FF, to increase the nonlinear pump efficiency (see figure 8(a), [109]), or with the harmonic wavelength to enhance the emission cross-section (see figure 8(b), [110]).
In nonlinear frequency-conversion applications, nearinfrared pump lasers are frequently used to drive coherent frequency upconversion processes. Plasmonic nanostructures based on noble metals such as Au and Ag, possess plasmon resonances that are radiative dipolar modes, located within the near-infrared band. This indicates that significant radiation losses are present for the pump energy and hence only limited nonlinear conversion efficiencies can be achieved. To tackle this issue, Fano-type plasmon resonances can be created in metal nanostructures with controlled interference between a super-radiant dipolar mode and a sub-radiant higher-order plasmon mode [111]. As shown in figure 8(c), at the Fano resonance wavelength in a gold split-disk nanostructure, the radiative loss of the dipolar component is significantly suppressed due to near-field coupling mediated excitation of the dark plasmon mode, resulting in strong near-field localization and enhanced second-harmonic generation [108].
Plasmonic nanostructures with multiple Fano resonances are indeed required to achieve field enhancement at multiple frequencies, which particularly benefits the multiple-wave mixing process. In such systems, the plasmonic coupling between the constituent elements provides multiple subradiant high-order modes to destructively interfere with a broad dipolar mode [112,113]. Employing symmetry considerations, we have observed that a plasmonic crossbar nanostructure sustains multiple Fano resonances and polarization independent second-harmonic generation, as per figure 9(a) [107].
Multiple Fano resonances have also been employed to enhance third-order nonlinear processes in a metal nanodisk cluster [114]. As shown in figures 9(b) and (c), the two incident pumps with frequency ω 1 and ω 2 precisely match the two Fano  -type nonlinear metasurface. To visualize the image of characters META written on the metasurface, the read laser beam has to be circularly polarized and with wavelength at the plasmon resonance of the single pixel-a 3  nanostructure pair; in the meanwhile, it requires the detection wavelength at half that of the pump beam. Reprinted with permission from [105]. Copyright (2017) American Chemical Society. resonance energies of the nanodisk cluster, leading to readily detectable FWM signal at 2ω 2 −ω 1 . In general, multiple Fano resonances exist in a relatively large plasmonic nanostructure that typically possesses a broad dipolar bright mode. Therefore, a balanced trade-off between the increased dipolar radiation loss and the strong field localization of sub-radiant modes needs to be considered for optimized nonlinear enhancement.
Improving the nonlinear excitation efficiency with nearfield enhancement can also be employed to plasmonic nanocavities sized down to a few nanometers. A promising class of such plasmonic platforms consists of metal-film-coupled nanoparticles, which feature an extremely strong field confinement and enhancement within the particle-film gaps [115][116][117][118]. Benefiting from well-established fabrication technologies and their controllable plasmonic properties, important advances into the design of these nanostructure configurations have been made for nanoscale nonlinear optics [119][120][121].
More advanced plasmonic nanostructures have been designed in order to improve the nonlinear frequency-conversion efficiency by simultaneously increasing the pump efficiency at the fundamental wavelength and the scattering efficiency at the harmonic wavelength. For example, a wavelength-scaled optical antenna with multiple plasmon resonances across a broad spectral band has been fabricated for this purpose, the corresponding details being presented in figures 10(a) and (b) [122]. Due to the field enhancement within a wide spectral range covering both the fundamental and SH wavelengths, this plasmonic structure has demonstrated an effective second-order susceptibility value comparable to that of widelyused nonlinear inorganic crystals. However, such a broad plasmon band is typically not available for a subwavelength metallic nanostructure chiefly because of the reduced radiation decay.
To realize efficient nonlinear light sources at the nanoscale, recent studies have employed mode matching in multiresonant plasmonic nanostructures [123][124][125][126]. An early study with aluminum nanostructures has demonstrated a doublyresonant (DR) optical antenna in which two discrete plasmon resonances spectrally match the fundamental and SH wavelength, respectively [123]. As a result, both the second-harmonic generation and its re-emission into the far-field are significantly increased when compared to the case of a single dipolar antenna. In addition, the plasmonic mode character is crucially important for nonlinear DR nanostructures. For instance, with the same mode character at the emission wavelength, a Fano resonance at the fundamental wavelength can be more efficient to boost the nonlinear excitation efficiency than a dipolar one [124], a feature illustrated in figure 10(c). Moreover, the plasmon resonance at the harmonic wavelength of a DR nanostructure is typically a highorder mode, which is intrinsically sub-radiant and therefore, limits the scattering efficiency of the harmonic waves.
Thus far, the dipolar second-harmonic emission from DR nanoantennas has only been demonstrated in a few coupled plasmonic systems. The V-rod shaped nanostructures of this type shown in figures 10(d)-(f) can be a representative model to understand the enhanced dipolar SH response [125]. To be more specific, the plasmonic response of such nanostructure possesses a pronounced dipolar resonance V 1 , which originates from the dipolar plasmon bonding between the two constituents, and matches the fundamental wavelength for nonlinear excitation enhancement. The bonding and antibonding between the quadrupolar mode of the V element and the dipolar mode of the nanorod create the hybridized modes V B 2 and V A 2 , respectively. Through plasmon hybridization, both of these two modes possess a super-radiant dipolar character inherited from the nanorod plasmon resonance. By further matching the radiative V A 2 mode with the SH wavelength and providing an excellent spatial mode overlapping between the fundamental and SH mode, this coupled system demonstrates a SHG conversion efficiency as large as 5×10 −10 W −1 .
A few recent studies have shed light on the underlying mechanism responsible for harmonics generation enhancement in DR plasmonic systems and thus provided new design guidelines. It has been found that the maximal SHG enhancement in an aluminum DR nanostructure shown in figure 11(a) occurs when the quadrupolar and dipolar components of the SH mode are strongly coupled with each other through an intra-antenna gap [127]. More specifically, the SH quadrupolar mode, which is directly excited by the pump wave, can resonantly transfer its energy to the SH dipolar mode supported by the same antenna through near-field coupling, and then efficiently radiates into the far-field.
Another elegant configuration of a DR nanoantenna for efficient SHG, presented in figure 11(b) [128], consists of a gold nanorod as an ω-particle and a smaller aluminum nanorod as a 2ω-particle. The intra-antenna position dependent SH response of this coupled nanostructure renders a feasible approach to probe the near-field of second-harmonic light around the ω-particle with the dipolar 2ω-particle. Note that the probing of SHG directly in the near-filed of plasmonic nanoparticles has been demonstrated in earlier studies, too [129].
Similar conclusion has also been drawn for a gold DR nanoantenna consisting of a single bar-like ω-particle and two disk-shaped 2ω-particles, as per figure 11(c) [126], where the location of each 2ω-particle relative to the ω-particle determines the radiation phase at the second-harmonic wavelength. When the three particles are positioned in a non-centrosymmetric arrangement (CIII in figure 11(c)), the light emitted by the two 2ω-particles constructively interfere to generate a linearly polarized, super-radiant dipolar SH emission. By contrast, the centrosymmetric configuration CIV enforces destructive interference of the light emitted by the two 2ωparticles to leave only a subradiant remnant of quadrupole scattering at the SH wavelength.

Plasmon-enhanced nonlinear effects in nonlinear active materials
Besides the intrinsic nonlinear response of plasmonic metals, in many recent studies it has been explored the plasmonic enhancement of extrinsic nonlinear response in nonlinear optically active materials. Taking advantage of the plasmonmediated field confinement, the conventional nonlinear optical materials incorporated into the near-field amplification region of plasmonic nanostructures experience significantly enhanced nonlinear pump intensity and hence exhibit increased nonlinear conversion efficiency in harmonics generation and multiple-wave mixing. For example, it has been demonstrated that lithium niobate (LN) having a high intrinsic second-order nonlinearity generates 20 times stronger SH intensity when filled in a plasmonic nanoring resonator, as per figure 12(a) [130], and similar conclusions hold true if instead of the LN one uses GaAs [43,44].
The perovskite-type material barium titanate (BaTiO 3 ) has a non-centrosymmetric crystalline structure and thus a bulk second-order nonlinear response. Recent studies have reported a 15-fold enhancement of second-harmonic generation in a single BaTiO 3 nanoparticle coupled to a metal nanoparticle in a dimer configuration [131]. The metal nanoparticle in this hybrid dimer acts as a nonlinear nanoantenna for the BaTiO 3 nonlinear source, as illustrated in figure 12(c). In such a configuration, the SHG enhancement is only possible at the interface and a limited fraction of the bulk SH from the active source radiates through the plasmon channel.
In order to fully utilize the bulk nonlinearity for maximal nonlinear response, it has been proposed to embed BaTiO 3 into a core-shell nanocavity [132]. Using this approach, a SHG enhancement factor of 500 has been experimentally demonstrated (see figure 12(b)). Recently, this core-shell nanostructure configuration has also been adopted for investigating DFG at the nanoscale [133]. As depicted in figure 12(d), such elegantly engineered composite nanostructure supports plasmon resonances at pump, signal, and idler wavelengths, providing large field enhancement across the confined BaTiO 3 medium and efficient coupling of the frequency-converted idler emission to the far-field.
Although metals possess one of the largest intrinsic thirdorder susceptibilities in nature, their skin depth is generally small and the third-order nonlinear response is hence strongly reduced. Combining plasmonic metal nanostructures with non-metallic nonlinear optical materials represents a new paradigm for enhancing third-order nonlinear optical effects. For example, by feeding a semiconductor indium tin oxide (ITO) nanoparticle into the near-field hot spot of a gold dipolar antenna, as depicted in figures 12(e) and (f), the resultant hybrid dielectric-plasmonic nanostructure exhibits a third-harmonic conversion efficiency of up to 0.0007% [134]. Further analytical results reveal that the total third-harmonic response originates predominantly from the active ITO nanoparticle rather than the gold nanoparticles [135].
In recent years, we have witnessed an increasing research interest in nonlinear optics of nanostructured high-index dielectric materials [136][137][138][139][140], largely stimulated by their negligible dissipative losses, large nonlinearity, and strong magnetic multipolar resonances in the visible and near-infrared spectral band [141]. Compared to plasmonic metal nanostructures, resonant high-index dielectric ones can produce enhanced field over more extended spatial regions and thus, can mediate nonlinear optical processes in significantly larger volume. With specific arrangement of a plasmonic element and a high-index one in a hybrid subwavelength structure, as illustrated in figure 13 [142], the strongly confined near-field mediated by the plasmonic element can further excite the resonant high-index mode, resulting in significantly improved nonlinear conversion efficiency due to the enhanced nonlinear polarization over the whole high-index volume. In addition to enhanced nonlinearities, high-index dielectric nanostructures with rich spectra of Mie resonances can be conveniently manipulated for targeted nonlinear emission characteristics, including vector optical field generation and directional radiation [143,144]. In order to find out about other technological applications offered by dielectric nanostructures in nonlinear nano-optics, interested readers can refer to the publications [94,141] and other related works.
Similar to the plasmon-assisted coherent frequency conversion of active nonlinear optical materials, upconversion nanocrystals (UCNCs) have also been hybridized with plasmonic metal nanostructures for enhancing upconversion efficiency and modulating emission state. Similar to the aforementioned mode-matching configuration for enhanced harmonic generation in plasmonic nanostructures, simultaneous excitation and emission enhancement in upconversion luminescence can be expected for the UCNCs located in the close vicinity of a DR plasmonic nanostructure. To this end, it has been proposed that, as illustrated in figure 14(a), the longitudinal and transverse dipolar plasmon resonances of a gold nanorod can be geometrically tuned to respectively match the excitation (980 nm) and emission (540 nm) wavelength of lanthanide-doped UCNCs for double-resonance enhanced light upconversion [145]. Analytical results suggest that the distance between the UCNCs and the gold nanorod is crucially important to the upconversion efficiency, as the upconversion emission would suffer from energy quenching induced by non-radiative energy transfer to gold at very small metal-emitter separations [146].
Single gold nanorod@SiO @CaF : Yb , Er 2 2 3 3 + + hybrid core-shell-satellite nanostructures shown in figure 14(c) have been used to demonstrate the critical importance of controlling the distance between the plasmonic component and the UCNC emitters in maximizing plasmon-enhanced upconversion efficiency. It is found that the optimized distance is about 20 nm for the highest upconversion of the two emission bands (∼550 and 660 nm). Moreover, the two dipolar plasmon resonances of the gold nanorod can be geometrically tuned to simultaneously overlap with the two emission bands of the UCNCs, as per figure 14(b), which not only significantly enhances the emission intensity at both wavelengths, but also dramatically modifies their emission polarization character via optical antenna effects, as illustrated in figure 14(d). In addition, multiple lanthanides doping in UCNCs leads to rich physics for emission manipulation [147], where plasmon resonances can be used to selectively accelerate particular absorption and emission bands.
As illustrated in figures 14(e) and (f), the 808 nm excitation of the tri-doping Nd : Yb : Er 3 3 3 + + + UCNCs attached to Each nanobar (ω-particle) possesses a plasmon resonance at the fundamental wavelength and each nanodisk (2ω-particle) with a resonance at the SH wavelength. The non-centrosymmetric configuration CIII exhibits the strongest SH emission intensity (bottom). Reprinted with permission from [126]. Copyright (2016) American Chemical Society.
gold nanorods (with 808 nm longitudinal plasmon resonance) generates a much larger upconversion intensity than that of the 980 nm excitation seen in figure 14(f), although the latter was able to mediate an efficient Yb Er 3 3  + + energy transfer.
This is because the absorption rate of 808 nm photons by the Nd 3+ ions can be effectively accelerated by the longitudinal plasmon resonance of the gold nanorod. The absorbed photons are then non-radiatively transferred to the Yb 3+ ions after undergoing a fast relaxation to an intermediate state close to the excited state of Yb 3+ , leading to an enhanced sensitization efficiency. Because the transverse localized plasmon resonance of the gold nanorod matches the two green emission bands of Er 3+ , indicating an accelerated radiative decay rate, they can thus benefit from plasmon-induced simultaneous sensitization and emission enhancement while the blue and red emissions of Er 3+ undergo only an excitation enhancement. Such different plasmonic modifications to the three upconversion emission bands enable controlling their relative emission intensities and therefore realizing tunable upconversion colors [147].

Nonlinear graphene plasmonics
Graphene has attracted a tremendous research interest in science and engineering due primarily to its remarkable physical properties and potential for future technological  applications [23,24]. Regarding its optical properties, a key feature of graphene that makes it particularly appealing for photonics applications is that it supports SPPs [25,26]. This provides an unique material platform to study these excitations in truly 2D physical systems. In addition to the linear physical properties, nonlinear optical properties of graphene have been the topic of intense research, too. Graphene, as a centrosymmetric material, exhibits induced second-order nonlinearity, large third-harmonic generation, and strong optical Kerr nonlinearity in a single atomic layer. This makes it possible to employ graphene in active photonic devices with advanced functionality, including optical limiters, frequency converters, ultracompact modulators, and photovoltaic and photoresistive devices. The linear and nonlinear optical properties of graphene as well as its applications to active photonic devices are reviewed in this section.

Linear optical properties of graphene
The linear optical properties of graphene are described by the linear surface conductivity, whose frequency dependence is given by the Kubo's formula [148]. Within the random-phase approximation, this formula can be expressed as the sum of the inter-band and the intra-band contributions. The intraband part of the conductivity is given by where ω is the angular frequency, μ c is the chemical potential, τ is the relaxation time, T is the temperature, e indicates the electron charge, k B is the Boltzmann constant, and ÿ represents the reduced Planck's constant. Moreover, if μ c ?k B T, the inter-band contribution can be approximated as It can be seen that the intra-band contribution to the permittivity of graphene at THz frequencies is similar to that of noble metals, meaning that it is described by the Drude model. On the other hand, the inter-band part is similar to that of semiconductors, namely it can be represented by the  figure 15(a) and corresponds to ε F =0.6 eV, τ=0.25 ps/2π, and T=0.
From a computational point of view, it is of practical interest to introduce bulk quantities equivalent to the surface ones, as it is usually more convenient to work with bulk quantities rather than with surface ones. Consequently, one uses a bulk conductivity, h b s eff s s = , where h eff is the effective thickness of graphene, instead of the sheet conductance, σ s . This approach is ambiguous to some extent because the thickness of an atomic monolayer cannot be properly defined, especially in an electromagnetic context. Moreover, the optical properties of graphene can be described in an equivalent way using the electric permittivity, ò, which is related to the optical conductivity σ b via the following equation:

Nonlinear optical properties of graphene
The graphene lattice belongs to the h 6  space group (see the inset in figure 15(a)), and as such it is a centrosymmetric material. A direct implication of this property is that quadratically nonlinear interactions, such as SHG, are forbidden; however, THG and other cubic interactions are allowed and particularly strong in graphene [27,28]. The quadratic optical nonlinearity of graphene can be described using the nonlinear optical conductivity tensor, s 3 s ( ) , which relates the electric field, E, at the FF and the nonlinear surface current density, j nl , at the third-harmonic (TH). More specifically, assuming that the graphene sheet is located in the (x, y) plane at z=0, the nonlinear current density, J nl , can be expressed as [28]: where r t is the position vector lying in the graphene plane.
In these conditions, the nonlinear surface current density is given by: . A formula for s 3 s ( ) , derived under the assumptions that electron-phonon and electron-electron scattering processes together with thermal effects can be neglected, has been derived in the following form [28]: is the Fermi velocity, with γ 0 =2.7 eV being the nearestneighbor coupling constant and a 0 =1.42 Å the separation distance between nearest-neighbor carbon atoms in graphene. Figure 15 . This nonlinear conductivity is particularly useful in experimental studies of nonlinear optical phenomena in graphene because it is related to a physical quantity that is usually measured experimentally, namely the effective bulk third-order nonlinear susceptibility, b 3 c ( ) . Using the fact that in the case of harmonic fields i J r P r nl nl w = -( ) ( ), where P r nl ( ) is the nonlinear polarization, one can readily prove that where Ω t =3ω 0 is the frequency at the TH, with ω 0 being the FF.

Enhanced nonlinear optical response in optical gratings containing graphene
One approach to enhance the nonlinear optical response of graphene is to incorporate graphene structures in optical diffraction gratings that support optical modes at the frequencies of the graphene plasmons [149]. One such approach is illustrated by the structure schematically presented in figure 16. It consists of a slab waveguide that is periodically patterned and covered by a periodic array of graphene ribbons with width w=230 nm. The three graphene ribbons located onto the dielectric components of the slab waveguide with permittivity r  ¢ and ò r are called the inner and outer ribbons, respectively. This is inspired by the structure of the unit cell shown in figure 16. The center-to-center distance between both the inner and outer ribbons is w 3 0.69 m m = whereas the centerto-center distance between an inner ribbon an an adjacent outer ribbon is 1.37 μm. Moreover, the periodically structured waveguide consists of a slab of height, h=1.5 μm, placed in-between a substrate with relative permittivity 1.46 s 2  = (index of refraction, n s =1.46) and the graphene ribbons. The slab itself is periodically patterned, namely it consists of alternating domains with permittivity ò r =4 (n=2) and 3.24 r  ¢ = (n′=1.8) with a period Λ=5.5 μm. The interaction between the graphene plasmons and a waveguide mode was investigated using numerical simulations in which the angle of incidence, θ, was varied from 0°to 20°while keeping constant the azimuthal angle, j=0°. The fundamental wavelength values vary over the second-order plasmon peak, namely from 6 to 9 μm. The absorption map computed using this approach is shown in figure 17(a). The spectrum corresponding to θ=0°shows a narrow peak with maximum of A=0.33, due to the excitation of the TM 0 waveguide mode at 8.05 m TM 0 l m = , and two broad peaks with maximum values of A≈6×10 −3 , due to the excitation of surface plasmons. If the angle of incidence, θ, increases, the spectral location of the to plasmons remains unchanged, but the spectral location of the waveguide mode resonance varies. Moreover, for θ>0, the TM 0 mode exists at two wavelengths, TM . Their distance from the TM 0 l peak increases as θ increases, and of the two branches of the TM 0 mode the one on the left is more efficiently excited. For θ=4.96°and θ=10.71°the mode wavelengths, and θ=9°, that is where the inherited TM 0 band crosses the intrinsic TM 1 band. The corresponding conversion efficiency enhancement is significant and solely due to the excitation of the intrinsic nonlinear modes; however, it is orders of magnitude lower than in the case of inherited effects. For example, the same intrinsic TM 1 mode away from the simultaneous resonance, at λ TH =2. 16 μm and θ=20°, yields TH radiation of intensity I I 7.13 10 TH 9 0 =´-.

Double-resonant nonlinear graphene optical gratings
An alternative approach to enhance the nonlinear optical response of graphene gratings is to design them in such a way that they have plasmon resonances both at the FF and TH [150].
The geometrical structure of such a generic graphene grating is schematically depicted in figure 18. It consists of a 1D grating with period Λ=100 nm, grating vector K=2π/Λ, and width of the graphene ribbons, w=50 nm. The graphene grating lies in the z=z s plane and is placed onto a dielectric  substrate with relative permittivity, ò s , the relative permittivity of the cover being ò c . For specificity, it is assumed that the substrate is made of glass with ò s =2.25, whereas the cover is air, so that ò c =1. The linear optical properties of graphene are determined by its surface conductivity distribution, x z , s s s ( ), which means that one can describe the spatial dependence of the relative permittivity of the graphene grating by a piecewise continuous function, x z , r  ( ). The incident light is TM-polarized and is normally impinging onto this binary graphene grating.
The wavelength dependence of the reflectance, R, transmittance, T, and absorption, A, of this grating, as well as the corresponding spectrum of the TH intensity are presented in figure 19(a). The intensity of the incoming wave was set to I 10 W m 0 12 2 = -, the TH intensity being measured relative to this reference value. The spectra depicted in figure 19(a) show several important features. Thus, one can observe that both the absorption and TH intensity spectra possess a series of spectral peaks, the resonance wavelengths in both cases being the same. The magnitude of the corresponding spectral peaks, both in the case of absorption and TH intensity, decreases with the wavelength. Moreover, the spectral width of the resonances, and therefore the corresponding optical losses, decreases as the resonance wavelength decreases too. The origin of these spectral resonances can be traced to the enhancement of the local optical field because both absorption and TH intensity are directly related to the near-field at the FF. In addition, since graphene can be viewed as an extremely thin metallic film, one concludes that these resonances and the corresponding near-field enhancement at the FF are due to the excitation on the graphene ribbons of localized surface plasmons. The spectra depicted in figure 19(a) reveal a particularly effective way to design efficient graphene gratings for enhancement of the THG and other nonlinear optical processes in such diffraction gratings. To be more specific, let us assume that one can engineer the spectral resonances of the grating by changing the width of the ribbons in such a way that the resonance wavelength of the fundamental plasmon is equal to the wavelength of the incoming wave, whereas the wavelength at the TH is equal to the resonance wavelengths of one of the higher-order plasmons. This means that the resonance wavelength of one of the higher-order plasmons should be equal to a third of the resonance wavelength of the fundamental plasmon mode. In these conditions, the optical grating will be efficiently excited at the FF, which will induce a large field enhancement at this frequency, and will radiate efficiently at the TH because there is a plasmon resonance at the TH wavelength, too. In other words, such an optical grating can be viewed as a highly effective receiver at the FF and a strong emitter (an efficient antenna) at the TH. These ideas have recently been employed to design graphene optical gratings that can be used not only to enhance the THG but also to control the polarization state of the generated TH [151].
These theoretical arguments are validated by the results summarized in figure 19(b), where the dispersion map of the TH intensity of the diffraction grating is depicted, that is, the dependence of the nonlinear spectra on w. Also presented in the inset of this figure is the variation of the TH intensity with  the width of the graphene ribbons, calculated at the resonance wavelengths of the fundamental plasmon. An important conclusion illustrated by this figure is that, as we previously mentioned, the excitation of localized plasmons in graphene ribbons at the FF leads to a large increase of the TH intensity via local field enhancement. This can be clearly seen from the location of the spectral resonance bands in figure 19(b). However, more importantly, the plot presented in the inset of this figure demonstrates that a further enhancement of the TH intensity occurs when the double-resonance condition is verified. Specifically, the curve in the inset of figure 19(b) was determined by choosing the wavelength of the incident wave so as to coincide with the resonance wavelength of the fundamental plasmon (the plasmon band corresponding to the largest enhancement of the TH) and varying the width of the graphene ribbons. It can be seen that a maximum TH intensity is obtained when w=85 nm, namely exactly for the width at which plasmon modes exist at both the FF and TH.

Numerical methods for modeling nonlinear optical processes in plasmonic nanostructures
In the last decade, the bottom-up design of new optical nanostructures and the modeling of photonic nanodevices have been greatly accelerated by the availability of userfriendly, computationally powerful software tools. Thus, high fabrication costs of complex optical metamaterials make it imperative to have access to computational tools based on efficient numerical algorithms, which can greatly accelerate the device design process and reduce the design-fabricationtesting cycle. High-performance computing tools that are beginning to be used in nanophotonics are able to model a distributed multifunctional device 'continuum' in which the light is generated, processed, and collected in the same physical space. This device complexity implies that present and future computational algorithms for nonlinear plasmonics must fully integrate the description in the time and frequency domains of classical and quantum physical phenomena that take place over multiple spatial and temporal scales.

Time-domain methods and hydrodynamic model
Local-response approximation models, such as the Drude and Lorentz models, are commonly used to numerically model plasmonic nanostructures. However, these models do not consider the nonlinear and nonlocal optical characteristics of such nanostructures. At present, nanoscale fabrication techniques for plasmonic nanodevices are becoming increasingly more performant, which allows for investigating new phenomena occurring in unexplored physical regimes. When an optical field strongly couples to a plasmonic nanostructure, it interacts with free electrons in the system resulting in intrinsic nonlinear optical response, including higher-harmonic generation, sum-and DFG, and Kerr effects. Owing to the limitations of local-response models, physical models that incorporate both nonlinear and nonlocal effects are of particular importance.
The motion of free electrons in the conduction band can be described by the Boltzmann equation (BE) [152]: where f is a probability density function defined in the r k , ( )-phase space, E k ( ) is the energy of the electron with wave vector k, and the wave-particle duality is used, i.e.
The functions S and f t coll ¶ ¶ | in (27), are the source and collision terms, respectively.
When electrons are coupled to electromagnetic fields, the force F can be written as where E and H are the electric and magnetic fields, respectively, and v is the electron velocity.
Finding numerical solution to the BE is difficult and computationally expensive. Thus, in turn one solves its momentum balance equations. To this end, one defines the charge, current density, and kinetic energy density, respectively, as: The above equations can be written as a generalized source form, where f can be 1, −qv, etc.
Using the Galerkin scheme, we can multiply (27) by f and then integrate it over k, to finally get: å å can be rewritten as a series of balance equations: If one chooses f=1, the zeroth-order balance equation reads: where the generalized source term S f is assumed to be zero. In particular, the above is the current-continuity equation.
Moreover, for qv f = -, the first-order (momentum) balance equation reads: is the polarization current, τ m is the relaxation time, and m n vv is the kinetic energy density. In equation (36), the relaxation time approximation is used for the collision term [152]. This equation can be simplified by using the current-continuity equation (35), the result being the hydrodynamic equation: The nonlinear and nonlocal response of the conductionband electrons in metallic structures can be modeled by selfconsistently solving the hydrodynamic equation (38), currentcontinuity equation (35), and the Maxwell equations: where ò 0 and μ 0 are the vacuum permittivity and permeability, respectively. The second-harmonic polarization sources in isotropic and centrosymmetric metallic materials are classified into bulk and surface terms [70,74,75] (see also section 2.2). For ωτ m ? 1, the hydrodynamic model gives the following formulae for the bulk parameters: c ( ) is the linear bulk permittivity of the metal, n 0 is the equilibrium electron density, and ω p is the plasma frequency of free electrons. Moreover, in the same limit ωτ m ? 1, the following formulae for the components of the nonlinear surface susceptibility can be derived from the hydrodynamic model: is negligible. In a recent study [153], the second-harmonic response from metallic nanoparticles with different shapes has been analyzed perturbatively by solving the Maxwell-hydrodynamic model with the finite-difference time-domain (FDTD) method. The predictions of the theoretical model qualitatively agreed with experimental measurements. Figure 20 shows the SH signal strength emitted from C-shaped and T-shaped particles. The polarization state of the far-field SH signal is always y-polarized. Thus, a universal selection rule exists, i.e. mirror symmetry of the metallic particles in one direction completely forbids a polarization component of second-harmonic generation in that direction. Also, the presence of structural plasmonic resonances can also greatly enhance the SHG from metallic particles.
Plasmonic resonances of different order make different contributions to the SHG. Using the perturbative solution, it was found that the free-electron polarization P has the nonlinear sources (including the Coulomb term also referred to as quadrupole-like term), proportional to E P  ( · ), the magnetic Lorentz force contribution, P Ḃ , the convective terms P P  (˙· )˙and P P  ( ·˙)˙, the linear nonlocal pressure term described by P   ( · ), and the nonlinear pressure term P P    ( · ) ( · ) [154].
Using this formalism, SHG from a gold wire of diameter d=100 nm with the pump wavelength λ=1064 nm was modeled, the contributions from different nonlinear source terms being presented in figure 21. Each term defines a nonlinear scattering channel that could be amplified or suppressed by the other terms via interference effects. Consequently, it is not trivial to determine the contribution of an individual term to the generated field, which, in general, will depend on the particular geometrical configuration and other parameter choices. Additionally, in both [153,154], it was assumed that the fundamental fields are not affected by the generated harmonic, which is called the undepleted pump approximation.
The Maxwell-hydrodynamic model was rigorously solved by various numerical methods, such as FDTD method [155,156], finite-element time-domain method [157], and discontinuous Galerkin time-domain (DGTD) method [158]. The Maxwell equations or wave equation is coupled to the hydrodynamic model to produce a self-consistent solution going beyond the undepleted pump approximation. Thus, a time-split semi-implicit FDTD method to model SHG in metallic metamaterials has been developed [155]. Then, the finite-element time-domain method [157], has been used to investigate the nonlinear scattering spectra from an infinitely long gold cylinder of 100 nm radius. The power dependence of the nonlinear harmonic generation follows the well-defined powers of the pump, quadratic and cubic for SHG and THG, respectively, as expected for nonlinear optical interactions (see figure 22). In particular, it has been found that the straightforward perturbative hydrodynamic description (e.g. relying on phenomenological second-order nonlinear polarization) is inconsistent with the self-consistent Maxwell-hydrodynamic model.
The DGTD method has been employed to simulate a complex V-groove nanostructure, in particular to engineer a double-resonant scenario at the fundamental and the second-harmonic wavelengths. The two main resonances occur when exciting the structure with light polarized along the long axis of the groove. The first resonance ω 1 is located at 512 nm while the second one is found at half the wavelength, i.e. at 256 nm (marked with ω 2 ). This double-resonant scenario leads to an enhanced SHG signal that is slightly blue-shifted with respect to ω 2 , to 253 nm (marked with 2ω 1 , ω 2 in figure 23). In addition, it has been observed an additional SHG resonance at 223 nm, which is related to the resonance with light polarized along the short axis of the groove. This demonstrates that the SHG process is sensitive to all (linear) resonances and not only to those that occur for a specific excitation.
The selection rules of SHG from periodic metallic metamaterials have been recently investigated numerically by using an explicit FDTD method [156]. For a unit cell with N-fold rotational symmetry, the polarization selection rule for SHG is nN n m r -+ + = , where n is an arbitrary integer, ν is the spin angular momentum of the transmitted secondharmonic wave, and μ and ρ are the spin angular momenta of fundamental incident waves. The parameters ν, μ, and ρ are either 1 or −1. If the unit cell has three-fold rotational symmetry, namely N=3, only the combinations 1 n = , μ=−1, ρ=−1 or 1 n = -, μ=1, ρ=1 are allowed. It means that the metamaterial could convert left-circularly polarized fundamental waves to right-circularly polarized second-harmonic waves, and vice-versa. If N=1, both 1 n = -, μ=−1, ρ=−1 and ν=1, μ=1, ρ=1 are allowed. It means that the amplitude of the right-circularly polarized component is comparable to that of the left-circularly polarized component at the second-harmonic frequency. For N>3, no combination of ν, μ, ρ is allowed and thus the SHG is forbidden. Figure 24 shows the corresponding results for the triangle, L-shaped, and disk-shaped metamaterials.
Advances in nanofabrication techniques have made it possible to achieve ultrasmall nanoparticles and separation distance between two metallic elements of only a fraction of a nanometer. At such distances, nonlocal, electron spill-out and quantum tunneling effects become non-negligible. To consider these effects, the classical hydrodynamic equation should be modified using the Thomas-Fermi hydrodynamic equation [159] or the quantum hydrodynamic equation [160][161][162]. The quantum hydrodynamic equation can be written as: . In the quantum hydrodynamic model, on the other hand, n r 0 ( ) is the quantum-mechanical electron density of the system under consideration, obtained, for example, from a full ground-state Kohn-Sham density functional theory calculation.
Second harmonic generation in metallic Archimedean nanospirals has been numerically investigated [159], revealing the interplay between nonlocal and geometric effects. The quantum pressure term in the nonlinear hydrodynamic model is responsible for the emergence of fractional nonlinear   figure 25(c)), the nonlinear scattering intensity (with the linear scattering field subtracted) clearly shows the signature of generation of higher-harmonics up to the 3rd order, although no significant impact of the nonlocality is observed. The effect of nonlocality, however, is much more pronounced in the case of the spiral nanostructure (red lines in figure 25(c)). First, the generation of fractional harmonics, due to the quantum pressure term, is evident. A pronounced and remarkable difference between local and nonlocal scenarios manifests itself in the broadband supercontinuum generation. It can be traced to the collective electron-electron interactions, which gives rise to distinct fractional harmonics.
Recently, hydrodynamic equations governing the collective motion of massless fermions in graphene have been derived [163][164][165]. However, numerical solution of the quasirelativistic hydrodynamic equations for analyzing the nonlinear optical response of graphene has not been explored yet.

Frequency-domain methods
The classical hydrodynamic model only requires the relaxation time and electron density as physical parameters but only captures linear and nonlinear dynamics of electrons corresponding to intraband transitions. Thus, to capture the electromagnetic response due to interband transition processes, one can employ frequency-dependent surface and bulk susceptibilities obtained from experiments and model the nonlinear processes by using frequency-domain methods. In this approach, one first solves for the electric fields at the pump frequencies ω 1 and ω 2 using the linear bulk susceptibilities. Then, one calculates the surface and bulk nonlinear polarization sources by using the nonlinear surface and bulk susceptibilities. Finally, employing the nonlinear polarizations as the excitation sources, one can compute the sum-frequency or difference-frequency fields using the linear bulk susceptibilities at the frequencies of ω 1 +ω 2 or   nanostructures. The boundary element method can model efficiently the second-harmonic generation from arbitrarily shaped plasmonic structures since it only requires surface discretization. In addition, experimentally measured optical constants can be readily incorporated in the calculations.
Second-harmonic process can be described by the following two coupled linear wave equations: is a pump source current producing the incident field at the fundamental frequency. The nonlinear polarizations P E E : and P E E : are the nonlinear source terms for the fundamental and second-harmonic fields, respectively. Here, sd 2 c ( ) and su 2 c ( ) are the surface second-order nonlinear susceptibility tensors related to the down-conversion and up-conversion processes, respectively, E w ( ) and E W ( ) are the fundamental and secondharmonic electric fields, respectively, whereas k w ( ) and k W

( )
are the corresponding wave numbers. In order to fully describe both the cross-coupling between the fundamental and second-harmonic fields and the pump depletion effect, equations (43) and (44) are solved self-consistently using a first-principle boundary element method with the initial The computational domain used for solving for the electromagnetic fields at the FF and SH is divided into the exterior background region and the interior of the metal, the two regions being separated by the interface, S. The object is illuminated by an incoming plane wave. By employing Love's equivalence principle, the corresponding currents lying on the outer side of the interface, S, generate the scattered field in the exterior region and no field inside the metallic object. Moreover, the equivalent currents on the inner side of S generate the total field inside the metallic object and no field in the exterior region. This is expressed mathematically by the following two equations: is the dyadic Green's function at the frequency ν (ν=ω, Ω=2ω represent the frequencies of the fundamental field and secondharmonic, respectively) and R r r = -¢ | | . Moreover, ℓ=e, i denotes the exterior and interior regions of the object, respectively, J ℓ The boundaries of the nanostructures are discretized using triangular meshes and the equivalent surface currents are expanded using Rao-Wilton-Glisson basis functions [166]. In the next step, a matrix system is constructed by exploiting the Galerkin testing procedure [167]. Then, a modified Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) formulation [167] is used to ensure accurate solutions even upon resonant conditions. The PMCHWT matrix equation can be expressed as: According to the nonlinear boundary conditions for the fields [69], the following relations hold: Here, ò′ is the permittivity in the selvedge region [69].
Equation (49) can be used to determine both the fundamental and second-harmonic fields by setting ν=ω, Ω, respectively. As indicated in (51), the source terms at the fundamental frequency include the incident excitation wave and the nonlinear polarization source back-coupled from the second-harmonic (down conversion process). The secondorder susceptibility mediates this interaction and acts as a coupling coefficient.
The interaction between the fundamental and secondharmonic fields is determined by solving iteratively the equations for the fundamental and second-harmonic fields, with the updated source terms computed in the previous iteration. In the numerical implementation of this procedure, one initially calculates the equivalent fundamental electric and magnetic currents induced by the incident excitation wave by setting the nonlinear polarization to P . In order to compute the initial value of the second-harmonic nonlinear polarization, one needs to evaluate the field at the fundamental frequency on the inner side of the surface, S. It is given by the equivalent surface currents via the following equations: Using this computed field at the fundamental frequency, the initial value of the nonlinear polarization at the secondharmonic is first calculated, so that subsequently the initial value of the second-harmonic field can be solved. After these initial values are calculated, one can iteratively solve for the field unknowns at the fundamental and second-harmonic frequencies until the iterative procedure converges asymptotically to the self-consistent solution.
A frequency-domain boundary element formulation for computing surface second-harmonic generation from nanoparticles of practically arbitrary shape and material properties has recently been developed [168,169]. By using this method, it has been shown that the particle geometry in conjunction with the polarization of the incident field play a central rôle in determining the characteristics of SHG from arrays of metallic nanoparticles, as per figure 26. Thus, for sample L, the x-polarized input field is nonresonant and the hot spots are weak, hence generating a weak SH signal. For sample T b , the y-polarized fundamental field is nonresonant and thus the allowed yyy-interaction is weak. Moreover, for sample T c , input y-polarized fields are near-resonant. The y-polarized hot spot is highly asymmetric in the y-direction, which leads to a strong signal produced by the yyy-interaction configuration. Consequently, the yyy-configuration SHG signal from sample T s is twice as large as that from sample L. Both the symmetry properties of the structure with respect to the polarization of the incident field and the existence of resonances determine the amplitude of the second-harmonic signal radiated into the far-field region.
The boundary element method has been extended to incorporate both local-surface and nonlocal-bulk sources [170]. Using this formalism, it has been demonstrated that the contribution of the nonlocal-bulk sources can be described using equivalent surface electric and magnetic currents. A similar method based on periodic Green functions has been used to study SHG from periodic metallic nanostructures [171]. Using the boundary element method, it has been shown that SHG in a heptamer plasmonic system can be increased when the structure possesses Fano resonances [172]. Thus, the geometry of the system of nanoparticles is engineered to simultaneously possess a Fano resonance at the fundamental frequency, leading to a strong localization of the fundamental field nearby the system, and a higher-order scattering peak at the second-harmonic frequency. Figure 27 depicts the second harmonic near-field intensity when the fundamental wavelength is λ=800 nm and λ=950 nm. A dramatic enhancement of the near-field at the second harmonic frequency in the Fano-resonant system can be clearly observed: in this resonant condition, the near-field intensity is about 25 times larger (λ=400 nm, figure 27(a)) than in the off-Fanoresonant case, namely at λ=475 nm; figure 27(b).
A self-consistent boundary element method that goes beyond the undepleted pump approximation has been introduced in a recent study [173]. Using this method, a plasmonic particle-in-cavity nanoantenna was designed to achieve strongly enhanced and directionally tunable secondharmonic radiation [174]. The efficiency of the particle-incavity nanoantenna has been compared to that of other designs reported in literature [172,175], as shown in figure 28. x-and y-polarization components of the near fields at the second-harmonic wavelength calculated at a plane 15 nm after the gold nanoparticles. The shown quantity corresponds to the real part of the field at a moment when the real part of the strongest spot is maximum. Sample L (top row), sample T b (middle row), and sample T s (bottom row) for x-polarized, y-polarized, and x+y polarized input beam at the fundamental wavelength, respectively. The numbers in the right upper corners show the maximum field amplitudes normalized to that for the yyy component of L. The input polarizations are indicated above the figures and the calculated SHG components are shown in the left upper corners. The labels in the right bottom corner show the tensor components that contribute to each panel. The near fields related to allowed tensor components are framed with thick blue, red, and green lines. Reprinted with permission from [169]. Copyright (2015) American Chemical Society.
In the same context, a theoretical study of the characteristics of the nonlinear spin-orbital angular momentum coupling induced by SHG in plasmonic (gold) and dielectric (silicon) clusters of nanoparticles made of centrosymmetric materials has been reported [176]. In particular, a general angular momentum conservation law has been formulated for the nonlinear spin-orbital angular momentum interaction, which includes the quasi-angular-momentum of chiral structures with different-order rotational symmetry. The conservation law can be expressed as j s l qN sca s = + + ( ) , where σ and l are the spin and orbital angular momentum numbers of the incident optical beam, respectively, s denotes the order of the generated harmonic, that is s=1 for linear processes, s=2 for SHG, etc, N is the quasi-angularmomentum number characterizing the nanostructure with N-fold rotational symmetry, q=0,±1,±2,±3,Kis an arbitrary integer, and j sca is the total angular momentum number of the scattered field.
The second-harmonic conversion efficiencies corresponding to plasmonic and dielectric chiral clusters illuminated by left-and right-circularly polarized Laguerre-Gaussian (LG) beams with orbital angular momentum number l=4 were calculated and compared. It should be noted that the conversion efficiency of the plasmonic and dielectric chiral clusters was normalized to the cross-sectional area of the plasmonic and dielectric nanospheres, respectively. This normalization ensures a per-particle quantification of the nonlinear conversion efficiency. The power of the incident LG optical beam for the dielectric and plasmonic nanostructures was chosen to be same, P in =2 W. As can be seen in figure 29, the dielectric chiral cluster made of (high-refractive index) silicon nanospheres supports both electric and magnetic Mie resonances, which leads to similar conversion efficiency when compared to that of the plasmonic chiral cluster. Considering that due to significantly lower optical losses in dielectrics, the silicon cluster can sustain a much larger optical power, the dielectric nanostructures can provide orders of magnitude larger second-harmonic conversion efficiency, as compared to their plasmonic counterparts.
It is also worth noting that, using the volume integral equation, it has numerically been studied the SHG from dipole gold nanoantennas by analyzing the different contributions of bulk and surface nonlinear terms [71]. In a different study [177], the relative contribution to the SHG in dielectric (silicon) and plasmonic (gold) of bulk and surface nonlinear effects has been compared, the main conclusion of this work being that in the case of plasmonic structures the optical power generated by nonlinear interactions at the surface surpasses by orders of magnitude that due to bulk effects, whereas in dielectric structures, in certain circumstances, these two contributions can become comparable.

5.2.2.
Finite element method. The finite element method is another powerful tool for analyzing the optical response of Figure 27. Near-field distribution of the second harmonic intensity close to an ideal silver heptamer plotted on a logarithmic scale with SHG at (a) λ=400 nm (corresponding to the fundamental at λ=800 nm) and (b) 475 nm (corresponding to the fundamental at λ=900 nm), the angular intensity emission plots for the two corresponding cases of (c) SHG at λ=400 nm and (d) SHG at λ=475 nm. Note that in panel (d), the intensity has been scaled by a factor of 3 showing the much stronger SHG in the far-field for the Fano-resonant case. Reprinted with permission from [172]. Copyright (2013) American Chemical Society.  plasmonic structures, especially when the bulk effects cannot be neglected. In particular, it has been used to investigate the SHG of 150nm spherical gold nanoparticles [91]. Thus, it has been demonstrated that the interference effects between dipolar and octupolar plasmon modes can be used as a fingerprint to identify the local surface and nonlocal bulk contributions to the SHG. In a different study, it has been reported that experimentally measured SHG from gold nanoparticles with size of up to 100nm is in excellent agreement with predictions of finite element method simulations involving the normal surface term only in the nonlinear polarization source [178] The second-harmonic emitted by multiresonant plasmonic nanoantennas has also been analyzed [179] using simulations performed using a perturbative approach based on the undepleted-pump approximation, that is, assuming that the SHG field does not back-couple to the excitation field, and considering only the dominant contribution due to the freeelectron currents normal to the surface of the metal. Figure 30(a) shows a confocal map of the SHG from the nanostructure array. The plasmonic structure was excited with 100fs pulses centered at 1560 nm and the nonlinear signal was measured only within a narrow spectral window located at 780nm. A high emission occurs when the length of the V-shaped structure possesses optical resonances lying close to the wavelength of the excitation pulses. This is valid both in the case of isolated V-shaped particles (SEM image presented in figure 30(b)) and for coupled structures (right and left panels of figure 30(a), respectively). However, when the second-harmonic intensity collected from the coupled structures is compared to the one from the isolated V-shaped particles, a pronounced intensity modulation induced by the presence of the rod in the array of coupled antennas can be observed. A peak intensity for SHG, which is more than twice as large as the value obtained from a single resonant V-shaped antenna, is achieved for rods whose resonance wavelength matches the SH wavelength. The 2D SHG map computed numerically is presented in figure 30(c), and shows very good qualitative agreement with the corresponding map obteined experimentally and presented in figure 30(a). A full polarization investigation of the SHG was performed on the DR structure, and is illustrated in the angular plot depicted in figure 30(d). As expected from the SHG selection rules under strongly focused beam operation, the emission characteristics are those of an electric dipole. Further analysis of the spectrum in figure 30(e) reveals that the intensity of THG in these nonlinear plasmonic devices, which in nanostructured systems is often by far the dominant nonlinear optical interaction, is comparable to that of SHG.

5.2.3.
Mode matching schemes. The mode-matching method is a commonly used technique for the analysis of the electromagnetic properties of optical systems, especially for structures consisting of two or more piecewise optically homogeneous and separated regions. It is based on expanding the fields in a certain basis function and matching them at the boundaries between different regions. Among the most representative of the mode-matching methods, we mention here the rigorous coupled-wave analysis, scattering matrix method, plane-wave expansion method, and T-matrix method. Requiring relatively modest computer resources, these numerical methods are particularly useful in characterizing the linear and nonlinear optical response of periodic, multilayered, and multi-component structures. However, mode-matching methods are not particularly well suited for describing metallic scatterers with complex geometries because a very large number of modes are required in the field expansion in order to accurately describe plasmonic phenomena.
In the first step of the numerical implementation of such a method, one solves the fundamental field. To this end, the incident and the scattered electromagnetic fields are expanded in Fourier series using a properly chosen basis of functions, the choice being usually guided by the symmetry of the problem. Then, the boundary conditions (see equations (47), (51), and (52)) at the surface of the scatterers are used to (e) Plot of emission lines of the doubly-resonant nanoantenna in the visible-near infrared wavelength region, featuring a broad twophoton photoluminescence band between the SHG and THG peaks. Spectral analysis reveals that while the THG is centered around 519 nm, the SHG is emitted around 776 nm. Inset: overlap between the SHG peak (dark blue line) and the excitation laser band (red line). The theoretical SHG band obtained by self-convoluting the laser spectrum is also sketched (light blue line). Horizontal scales are expressed in hertz and the experimental and theoretical SHG peaks (FWHM ≈7.8 THz and 10.3 THz, respectively) are plotted on a frequency scale that is double the scale of the laser peak (FWHM ≈10.5 THz). In all measurements the excitation power is set to 50 μW, which is low enough to exclude any photodamage. Reprinted by permission from Springer Nature: Nature Nanotechnology [179], (2015).
construct a system of linear equations whose solution provides the Fourier coefficients of the series expansion of the scattered field. Once these coefficients are computed, by solving the corresponding system of linear equations, the electromagnetic field can be found at any point in space.
The second step consists in the calculation of the electromagnetic field at the second harmonic under the nondepleted pump approximation. One first determines the source of the field at the second harmonic, namely, the nonlinear polarization (surface and bulk polarizations) induced by the field at the fundamental frequency. Then, following similar procedures used in solving the field at the fundamental frequency, the field at the second harmonic can be readily obtained.
The multiple-scattering matrix method has been recently employed to rigorously describe the SHG from a set of cylinders with arbitrary electric permittivity, which are made of centrosymmetric (metal or dielectric) materials [180]. Regarding the optical properties of such structures, figure 31 shows that, depending on the frequency of the excitation wave, a chain of optically coupled plasmonic cylinders supports either modes that only propagate at the FF (panels A and B) or modes that propagate at both the fundamental frequency and second harmonic (panels C and D). Importantly, the latter type of plasmonic modes can find important applications to subwavelength active nanodevices for generation and transport of optical power at subwavelength scale. Another notable effect illustrated in panel A is the formation at the end of the chain of plasmonic cylinders of an optical beam with width of about λ/3, a so-called optical nanojet. This phenomenon can be employed to achieve subwavelength light focusing or optical nanoprobes. It has also been demonstrated that the surface second-harmonic generation can lead to the formation of nonlinear plasmonic whisperinggallery modes in microcavities made of metallic nanowires [181]. An unusual feature of these modes is that they possess fractional azimuthal modal numbers.
The multiple-scattering matrix method has also been employed to model SHG from 3D structures consisting of arbitrary distributions of metallic spheres made of centrosymmetric materials [182]. In particular, the interaction between circularly polarized LG optical beams and clusters of metallic nanoparticles has been investigated [183], the results being summarized in figure 32. Thus, it can be seen in this figure that the Stokes parameters S 1 and S 2 for the fields at the fundamental frequency and second harmonic have a symmetric spatial profile at the frequency of the resonant collective mode, which is one of the typical characteristics of the spin-orbit interaction of light. Moreover, the S 3 Stokes parameter characterizes the local density of the spin angular momentum. Therefore, the spin-dependent characteristics of the interaction of light with the cluster of spheres can be captured by the spatial distribution of the S 3 parameter. In particular, this is one of the signatures of the spin-Hall effect of light. Indeed, the spatial splitting of spin states for both the fundamental and second-harmonic fields has been observed for this type of clusters. The calculations show that the spin-orbit interaction can also cause the giant spin-Hall effect for the second-harmonic field.
Recently, a highly effective rigorous coupled-wave analysis for efficient and accurate description of both linear and nonlinear optical phenomena in nanostructured 2D materials embedded in nanopatterned photonic structures containing regular 3D optical materials, such as diffraction gratings and periodic metasurfaces, has been developed [149]. As discussed previously, graphene is a lossy conductor at THz and optical frequencies and therefore supports surface waves called plasmons. Examples of such plasmon modes of graphene disks are presented in figures 33(a)-(c). More specifically, in these figures, the dominant field component E x | | of the linear electric field corresponding to the first three plasmon modes is depicted. These field profiles exhibit distinct mode shapes with one, three, and five maxima, which demonstrates that they represent the first three localized plasmon modes of graphene disks. Strongly enhanced and Figure 31. The top two panels show the spectra of the scattering cross section corresponding to a chain of N=12 metallic cylinders. The radius is R=100 nm, the angle of incidence is 0°, and the separation distance is d=20nm. The spatial profile of the amplitude of the electric field, calculated at λ=313 nm (A and B) and λ=229 nm (C and D), is presented in the bottom panels. The panels A and C correspond to the fundamental, whereas the panels B and D correspond to the second harmonic. Reprinted figure with permission from [180], Copyright (2010) by the American Physical Society.
highly confined optical field resulting from the excitation of such localized plasmon modes gives rise to enhanced THG. This fact is supported by the field profiles plotted in figures 33(d)-(f), where the dominant electric field component at the third-harmonic, E x | |, is presented.

Conclusions and future perspectives
In this review, we have discussed previous and present research in nonlinear plasmonics in 3D (metallic) and 2D (graphene) nanostructures, a field that is now growing rapidly due to increasing interest in many exciting optical phenomena and a multitude of important emerging applications in surface science, active photonic nanodevices, nonlinear integrated photonics, and bio-medicine. In particular, we discussed the main theoretical tools, experimental techniques, and computational methods that are used in modern nonlinear plasmonics to study in an integrated manner nonlinear optical properties of metallic and graphene based nanostructures. The focus of this work is on nonlinear optical processes at the subwavelength scale and thus it goes significantly beyond the scope of traditional nonlinear optics in bulky nonlinear optical structures, in which case phasematching considerations play a central rôle. By contrast, in this new paradigm of nonlinear optics, the enhancement of the interacting optical near-fields and the strength of the overlap between optical modes with size comparable to or smaller than the optical wavelength are the defining factors that determine the efficiency of nonlinear optical interactions. The material presented in this review article is roughly equally divided among a background section covering the basic theoretical tools and concepts commonly used in nonlinear plasmonics, an experimental section that presents the latest developments in several key areas of nonlinear plasmonics, including nonlinear SPPs at metal surfaces, plasmon-enhanced nonlinear response of individual metallic nanostructures and arrays, phase-controlled nonlinear optics in plasmonic metasurfaces, mode-matching enhanced harmonics generation in plasmonic nanoantennas, and plasmon-enhanced nonlinear effects in nonlinear active materials, and a chapter in which we introduce the main time-and frequency-domain numerical methods used to model nonlinear optical properties of plasmonic systems and illustrate how they can be used to investigate computationally specific nonlinear plasmonic nanostructures and active plasmonic devices. Gathering in this way in the same place the main tools employed in nonlinear plasmonics research, namely theoretical modeling, experimental techniques, and computer simulations, makes it easy for the interested reader to form an informed opinion and readily gain valuable information about this dynamic field of science.
Nonlinear plasmonics is still a maturing field of research, so that its potential, both at a fundamental science level as well as for technological applications, is yet to be fully fulfilled. Therefore, at this point, it might be relevant to speculate briefly about future directions of research in nonlinear plasmonics. Thus, present and future technologies will allow one to fabricate metallic and graphene nanostructures with size significantly smaller than the optical wavelength and experimental techniques will make it feasible to individually probe the optical response of these plasmonic structures. As a result of these developments, nonlinear plasmonics will open up exciting new avenues to unique applications, such us subwavelength optical sources, ultra-small and highly efficient sensors, and optical imaging at deep subwavelength scale. Moreover, as the size of the building blocks of complex plasmonic systems that can be fabricated with available technologies decreases, quantum effects begin to play an increasingly important rôle in defining the physical properties of such plasmonic systems. Specifically, due to quantum coherence effects, the interaction between plasmonic structures in the quantum regime is strongly dependent on the frequency, optical power, and structural configuration of plasmonic quantum systems, which makes these physical systems to be highly nonlinear. An immediate consequence of this fact is that  quantum plasmonic systems possess remarkably large optical nonlinearity. Finally, it is expected that future numerical methods for nonlinear plasmonics will be increasingly more effective and versatile by incorporating not only the classical electromagnetic response described by optical constants such as linear and nonlinear optical susceptibilities but also by being able to describe in a unified manner the dynamics of the optical field coupled with the classical or quantum dynamics of the plasma carriers. All these mean that exciting new physics pertaining to nonlinear plasmonic phenomena are yet to be discovered, but of course, the most exciting developments will be those that we presently do not anticipate. ORCID iDs N C Panoiu https:/ /orcid.org/0000-0001-5666-2116 W E I Sha https://orcid.org/0000-0002-7431-8121 D Y Lei https:/ /orcid.org/0000-0002-8963-0193 G-C Li https:/ /orcid.org/0000-0002-9903-8900