Nonlinear graphene plasmonics

The rapid development of graphene has opened up exciting new fields in graphene plasmonics and nonlinear optics. Graphene's unique two-dimensional band structure provides extraordinary linear and nonlinear optical properties, which have led to extreme optical confinement in graphene plasmonics and ultrahigh nonlinear optical coefficients, respectively. The synergy between graphene's linear and nonlinear optical properties gave rise to nonlinear graphene plasmonics, which greatly augments graphene-based nonlinear device performance beyond a billion-fold. This nascent field of research will eventually find far-reaching revolutionary technological applications that require device miniaturization, low power consumption and a broad range of operating wavelengths approaching the far-infrared, such as optical computing, medical instrumentation and security applications.


Introduction
The last decade has witnessed the discovery of a wonder material called graphene. Graphene is an atomically thin layer of carbon found in lamellar-structured graphite. Although the band theory of graphene was postulated over 70 years ago [1], the successful isolation of a monolayer or a few layers of graphene was relatively recent [2]. Since then, there has been an explosion of graphene research in broad scientific and technological fields, ranging from electronic, photonic, mechanical to biological applications [3]. The breakthroughs achieved in graphene research are largely built on graphene's exceptional properties, such as its extreme mechanical strength and flexibility, high electron mobility and thermal conductivity, chemical purity and easy functionalization, to name but a few.
The existence of graphene's supreme properties, which were hidden in its graphite-crystal formation, is a manifestation of a naturally occurring, unique, two-dimensional (2D) crystal that physicists used to fantasize about. The special arrangement of the carbon atoms to form a monolayer of hexagonal lattices gave rise to its unique 2D band structure. This honeycomb network of carbon atoms is a result of the sp 2 hybridization between the s-and p-orbitals, which led to the formation of strong σ-bands that bestowed its robust mechanical strength. The remaining valence electron from each carbon atom then forms the half-filled π-bands, which enable the long-range π-conjugation that leads to graphene's extraordinary electronic, thermal and optical properties [1,4,5]. This review focuses on nonlinear graphene plasmonics, which has exceptional performance, as embodied by its linear and nonlinear optical properties. Graphene nonlinear optics and graphene plasmonics have been two separate popular research areas during the past decade [6,7]. The first experimental verification of graphene's ultrahigh nonlinear optical response was in 2010 [8] (figure 1), while graphene plasmons were first observed in 2012 through excitation by scattering at atomic force microscopy (AFM) tips [9,10] (figure 2). The rapid development in both fields in recent years has culminated in the study of nonlinear graphene plasmonics. We will learn here that the synergy between the graphene plasmon and nonlinear susceptibility unlocks graphene's enhanced nonlinear optical processes. We also review some of the proposed applications of nonlinear graphene plasmonics, which have leveraged nonlinear enhancement for lower power consumption and device miniaturization.  and zero doping, graphene has a minimum finite optical conductivity of σ min = e 2 /4h, which confirms prior theoretical predictions [12,13]. More generally, taking into account finite temperature and doping conditions, graphene's 2D optical conductivity could be represented by the Kubo formula; in the low-photon momentum and electron-scattering limit, Falkovsky derived a neat analytical formula [14,15] σ (1) g (ω) = ie 2 (2k B T) where e is the electronic charge, k B is the Boltzmann constant, T is the temperature,h is the reduced Planck's constant, ω is the angular frequency, E F is the Fermi level of graphene, and ν 1 and ν 2 are the scattering frequencies of the intraband and interband electrons, respectively. The first term on the right-hand side of equation (2.1) is the intraband conductivity, which is cast in a semi-classical Drude form, while the second term is the interband conductivity. When E F and T are zero, we can easily recover the minimum optical conductivity, σ min = e 2 /4h. While T is always assumed to be either at absolute zero for theoretical simplicity or at room temperature for practicality, E F , on the other hand, is leveraged for its reconfigurability. Since graphene is only an atomically thin layer, carrier doping of graphene is easily tuned by electrostatic gating [2]. An important implication for this dynamic reconstruction of graphene's optical conductivity is the ability to tailor the same graphene devices for situational uses; for example, operation at a wide range of optical frequencies from terahertz to the visible [16], or even electro-optic modulation by the electrical gate signals [17]. Recent innovative methods like surface carrier transfer have managed to tune graphene's E F to over 0.8 eV, which enables the coupling of near-infrared graphene plasmons [18].  Figure 3. Surface plasmons are the collective motion of free electrons that induces charge density oscillations. When a photon is coupled to the surface plasmon, they form surface plasmon polaritons that propagate along the metal surface. (Reproduced with permission from Zia et al. [19] (Copyright © 2006 Elsevier Ltd).) (Online version in colour.) electrons induces charge density oscillations that move the coupled photon energy along the metal surface, forming plasmonic waves, as shown in figure 3. As these are surface waves, they are inherently 2D waves, having oscillating electric fields only in the direction along the wave propagation path and perpendicular to the surface. An important consequence of this phenomenon is the absence of the concept of the diffraction limit: optical waves can now be squeezed to subwavelength scales, opening up a slew of new device concepts in nano-optics that will eventually revolutionize the field of nanotechnology [21].
The first observed plasmonic phenomenon on a metal surface was in the 1957 electron bombardment experiment by Ritchie [22]. From the electron energy loss data, Ritchie deduced the formation of surface plasmons with a surface plasmon resonance frequency equivalent to where ω p is the bulk plasma frequency of metal, given by ω p = Ne 2 /ε 0 m * e . From this, we see that the plasmonic response is entirely dependent on the material properties of the metal, i.e. the free-electron density, N, and the electron effective mass, m * e . One can also fit ω p to the Drude model to find the permittivity of metals, ε m = 1 − ω 2 p /ω 2 . Then, the plasmon dispersion is written as follows [19,20]: where ε d is the permittivity of the dielectric background. By nature, graphene is a semimetal. It is apparent from the solution of equation (2.1) that graphene behaves as a metal under a wide range of optical frequencies and Fermi levels. Hence, plasmonic behaviour on graphene surfaces is intuitively anticipated. Theoretical studies of graphene plasmons were carried out by Ryzhii [28]. (We will learn that the mid-infrared regime is particularly important due to graphene's Fermi level by substrate doping, usually oxide silicon, that sits in the range of 0.2-0.3 eV [29].) But, it was not until 2012 that graphene plasmons were successfully excited and observed [9,10].
In a similar vein, the dispersion of the 2D graphene plasmon is derived as follows [23][24][25][26][27][28]: In the long-wavelength limit, graphene's electronic susceptibility can be written as follows:   with N in accordance with E F =hv F √ π N (v F ≈ 10 6 ms −1 is the Fermi velocity and N is the 2D free-electron density of graphene) [30]; thus, we see a scaling dependence of ω p ∝ N 1 / 4 ; thirdly, it contains Planck's constant instead of the electron effective mass, which connotes the Dirac nature of the graphene plasmons.

(c) Extreme confinement of graphene plasmons
In figure 4, we plot the SPP dispersion for graphene (E F = 0.2 eV, assuming a silica substrate) and several representative metals: aluminium [31], copper [32], gold and silver [33], according to equations (2.2) and (2.3). Noticeably, graphene's plasmon dispersion levels off very early in the frequency range around 30 THz, or the wavelength equivalent of 10 µm, i.e. in the mid-infrared spectrum as stated earlier. Important information to be gleaned from the dispersion curves is the effective waveguide index, n wg , which is a good indicator of the confinement factor. The effective waveguide index is related to q and ω by the expression n wg = qc/ω. It is apparent from the dispersion curves in figure 4 that, for the same q, graphene has a much lower ω than metals by two orders of magnitude. For example, if we fix q = 30 µm −1 , a rough estimation of n wg has metals in the range of 5-20, while graphene can reach 300. In terms of wavelength compression, we see that, when a 10 µm free-space photon is coupled to the graphene plasmon, the plasmon wavelength is scaled down by a factor of 300 to only 33 nm.
The high confinement of graphene is also associated with the trade-off of high linear optical losses. Linear optical losses of the graphene plasmons are governed primarily by the intraband scattering frequency, given by ν 1 = ev 2 F /μ e E F [27], which is a function of the electron mobility of graphene and the Fermi level. A recent ultrafast optical experiment by Ni et al. [34] confirms that high-quality, high-mobility graphene samples do reduce the plasmonic losses of graphene. Even with the assumption of a high mobility of 10 4 cm 2 Vs −1 , the theoretical graphene plasmon losses are considerably larger than those of the metal plasmons [28,[35][36][37][38]. Nevertheless, this issue is less severe in nonlinear graphene plasmonics due to the unusually high saturable absorption property of graphene, as will be discussed in the subsequent sections.  photon. Therefore, it is of interest to the experimental optics community to model graphene as a 3D material in order to measure its optical response. An ingenious working model was first devised by Vakil & Engheta [28] using the effective medium theory. Graphene physically exists as a 3D structure-its atomic thickness is measured to be 0.33 nm [39,40]. Therefore, if we assume that the free electrons are bound to graphene's physical structure, it is instructive to model the effective '3D' optical conductivity as σ g−3D = σ g /d eff , where d eff = 0.33 nm is graphene's thickness. Hence, the dielectric permittivity could be easily written as ε g = 1 + iσ g /ε 0 ωd eff , and the refractive index is simply the square root of the dielectric permittivity. The original intent of Vakil and Engheta's model was to simulate graphene in commercial electromagnetic software; it was emulated with tremendous success by researchers working on device simulations thereafter [37,38,[41][42][43]. Later, experimental optics researchers also successfully adopted this model to fit their experimental measurements of graphene's refractive index [44][45][46][47][48][49]. It should be noted that the rapid development of graphene research has prompted improvement in modern Maxwell solvers to incorporate surface current boundary conditions for simulation of graphene's 2D optical conductivity in 3D devices [50][51][52].
For illustration, in figure 5, we plot the complex refractive index of graphene using the effective medium theory, for various Fermi levels at 1550 nm wavelength. For this wavelength, note that, at E F = 0.5 eV, there is a crossing point between the real and imaginary refractive indices. This is the so-called epsilon-near-zero (ENZ) point, where the insulator nature of graphene transits to a metallic one. However, care has to be exercised when designing ENZ-based devices, because graphene's refractive index is anisotropic. The anisotropic refractive index of graphene was recently discussed by Kwon [53] and de Oliveira & de Matos [54], simulated through ab initio density functional theory calculations by Cheon et al. [48] and experimentally confirmed by Chang & Chiang [55]. If we lay a graphene sheet on the x-y plane, we can write the refractive index tensor as follows: where the refractive index perpendicular to the graphene plane is usually taken as a constant with zero extinction coefficient. The anisotropic refractive index of graphene has a prominent effect on the design considerations for graphene optical devices that could lead to the device's performance being severely overestimated if the graphene's material model is applied incorrectly [54]. Finally, another quantity that could be characterized by the 3D model of graphene is the effective mode area of the waveguide, A eff . A eff can be used as another measure of the confinement strength of plasmonic waveguides. The expression for the A eff of graphene plasmonic waveguides   is simply written as A eff = ( |E| 2 dx dy) 2 / slab 1 2 |E z | 4 dx dy, where integration of the electric fields can be done through mode simulation. Figure 6 shows an illustrated comparison of gold versus graphene plasmonic waveguides [56], where matching parameters are the thickness of the waveguide (0.33 nm) and excitation wavelength (1.55 µm). Graphene is doped to 0.7 eV to ensure that it is metallic at the studied wavelength. With this configuration, it is found that graphene's A eff is smaller by two orders of magnitude than gold, a result which corroborates the findings for n wg . The high confinement of the optical plasmonic modes to the graphene ensures high light-matter interaction that enhances nonlinear optical processes, as we will cover later in this review.

Nonlinear optical properties of graphene (a) Theoretical studies of nonlinear optics of graphene
There is a wealth of literature authored by a number of research groups on the theoretical derivations of the nonlinear optical coefficients of graphene [6]. The series of studies on graphene's numerous nonlinear optical processes, as shown in figure 7, can be traced back to contributions by Mikhailov   where v g is the group velocity derived from the semi-classical dispersion relation v g = ∂ε k (p)/∂k, N(p) is the carrier population found through the Boltzmann transport equation and the product of both terms is integrated over momentum space p. While the Boltzmann transport model captures the physics of the intraband response, it does not sufficiently describe the interband response. Perturbative models are developed to supplement the interband physics, for example the quantum mechanical (QM) Dirac equation model [61,62,64,66,67,89], and also the semiconductor Bloch equations [73,75,76,85]. The QM model replaces v g in equation (3.1) with v g = ψ|∂ H/∂k|ψ , where ψ is the wave function of graphene and H is the Dirac Hamiltonian, which replaces N(p) with N(p) = f (−ε k ) − f (ε k ), the Fermi distribution difference between the conduction and valence bands. These models could solve for the interband response of graphene for high E F , yielding reasonably simple analytical results after approximations. However, they suffer from divergences as E F → 0. Later, a non-perturbative quantum mechanical model is developed by casting the time-dependent Dirac equation into the Bloch equations [70,71,88]. This model describes the interplay between the intraband and interband response, which is a more accurate description of the electron dynamics in graphene and does not diverge at very low E F .
The theoretical models mentioned above usually cast the optical Kerr nonlinearity of graphene in the form of the 2D third-order nonlinear optical conductivity, σ g . However, to compare with experimental measurements, the recasting of the 2D nonlinear conductivity into the 3D nonlinear refractive index via graphene's thickness is applied, as in the case of the modelling of graphene's linear refractive index. This is prominently practised by Cheng et al. [73,75,76] in several of their landmark papers on the theoretical derivations of graphene's nonlinear optical response. Quantitatively, the theoretical models predict a Kerr-type nonlinear coefficient of the order of n 2 ∼ 10 −11 -10 −13 m 2 W −1 around the visible and telecommunication wavelengths. This nonlinear coefficient is vastly larger than the reported average values of conventional silicon photonics at n 2 ∼ 10 −18 m 2 W −1 [90], and on a par with gold (at visible wavelengths) at n 2 ∼ 10 −12 -10 −16 m 2 W −1 [91]. Across wavelengths, using Cheng's model (with scattering frequencies ν 1 = ev 2 F /μ e E F [27], where μ e is graphene's electron mobility and ν 2 ≈ 0.8 THz [92]), we find that n 2 scales inversely proportional with E F at the mid-and far-infrared wavelengths, and converges at the near-infrared spectrum, as illustrated in figure 8. For k 2 , there are low variations of the magnitude with E F , but the sign swings from positive at the far-infrared to negative at the near-infrared, and the swing occurs at a distinct wavelength for each E F . The earliest device that exploits the giant nonlinear phenomenon of graphene would be the graphene saturable absorber, as shown in figure 9, which enabled ultra-efficient mode-locking of femtosecond laser pulses [92,101,102]. (A wealth of literature on graphene saturable absorbers can be found in [101,102] and thus will not be cited here.) More recently, when graphene was integrated into waveguides, a number of proposals to use graphene in waveguide-based devices include regenerative oscillation and FWM as shown in figure 10a [98,103,104], soliton propagation as shown in figure 10b [105], enhancement of nonlinearity and third-harmonic generation in a graphene-clad fibre [50,106,107] as well as improvements to third-harmonic generation and nonlinear absorption through field localization in a photonic crystal waveguide [108,109].

(c) Negative nonlinear refractive index of graphene
The early theoretical derivations of graphene's nonlinearity did not specifically mention the sign of the nonlinear coefficients. In 2014, a simple model of the relationship of graphene's nonlinearity with E F was derived by Khurgin [110] to be n 2,g = χ (1) η 0 (ev F ) 2 /2(ωE F ) 2 . When the non-retarded linear Drude susceptibility of graphene is applied, the equation can be reduced to n 2,g = 4πα 2 cv 2 F /E F ω 4 d eff , which indicates a positive nonlinear index of refraction. Ooi et al. [111] substituted the full expression of the linear susceptibility, which includes the electronic scattering frequencies, and found that graphene may exhibit negative nonlinear refraction for certain E F values and wavelengths. Subsequent Bloch derivations by Cheng et al. [75] also yield negative nonlinear conductivities for both the real and imaginary parts, when appropriate scattering frequencies and room temperature are applied, as shown in figure 8.  Similarly, while early experimental investigations did not explicitly distinguish between the signs of graphene's nonlinear coefficient, very recent works have unanimously found negative nonlinear refractive indices in the visible and telecommunication wavelengths. Three independent works in 2016, published during June, July and August, respectively, verified the negative index of refraction for graphene [112][113][114]. Demetriou et al. [112] measured n 2 to be of the order of −10 −13 m 2 W −1 at wavelengths between 1150 and 2400 nm with a 100 fs source using the Z-scan technique. Dremetsika et al. [113] also found n 2 to be −1.1 × 10 −13 m 2 W −1 using the optically heterodyne-detected optical Kerr effect (OHD-OKE), as shown in figure 11b, using a 1600 nm Ti : sapphire laser with 180 fs pulses, concurring with the previous Z-scan measurements. Finally, Vermeulen et al. [114] conducted SPM measurements, as shown in figure 11c, and also reported n 2 ∼ −10 −13 m 2 W −1 at 1550 nm for 3 ps pulses.

Nonlinear graphene plasmonics (a) Nonlinear graphene plasmonic waveguides
There are generally two classes of nonlinear graphene plasmonic waveguides. The first class uses graphene only for its robust linear optical properties, while the nonlinear optics is handled by an external active dielectric Kerr medium [115][116][117][118][119][120][121]. This class of waveguide uses the extreme optical confinement property of graphene plasmons to concentrate the fields into a small effective area, thereby increasing the local optical intensity in the Kerr substrate and hence increasing the overall nonlinear parameter of the waveguide. Relatively large shifts of the dispersion and resonant modes with optical intensity can be observed in the graphene plasmonic waveguides   bounded by Kerr substrates [115,119], as shown in figure 12, which increase the switching contrast of all-optical modulation. In addition, 'designer waveguides' can be engineered by using different waveguide configurations around the Kerr medium. For example, in parallel-plate graphene plasmonic waveguides, one could select either linear or nonlinear materials for each of the bounding media, and the nonlinear propagation characteristics vary with the configuration of the waveguide [116,118]. However, the analyses of these waveguides are not complete due to graphene's high nonlinearity, which would certainly overshadow the effects from the dielectric Kerr nonlinearities. The second class of waveguides exploits both the tight confinement of graphene plasmons and the ultrahigh Kerr nonlinearity of graphene [51,56,111,[122][123][124][125][126][127][128][129][130]. Here, squeezing of the plasmonic modes into a small effective area enables a surface-induced nonlinear enhancement (SINE) effect [122,131], as will be elaborated in the following section.

(b) Surface-induced nonlinear enhancement in graphene plasmonic waveguides
There are a number of theoretical papers discussing the effect of plamonics on graphene's nonlinearity [122,126,[132][133][134]. Pertaining to waveguide structures, as early as 2010, Skryabin et al. [131] derived a set of expressions from the nonlinear Schrödinger equations (NLSEs) to describe the surface effects of various subwavelength plasmonic waveguide structures and coined the term SINE. Later on, this derivation was applied to graphene [122,126], and, owing to graphene's ultraconfined plasmonic modes, the SINE effect could potentially contribute to the enhancement by a factor in excess of billions.
While the full formalism for the SINE effect is available in [122,131], it is possible to perform a snap evaluation for a simple graphene plasmonic waveguide bounded by a homogeneous background. From Gorbach [122], the total effective nonlinear parameter is given as Γ = gγ , where γ is the standard nonlinear parameter given by γ = ωn 2 /cA eff and g is the SINE factor given by g = (1 + η) −2 , where η = −2/Pq d · q d is the wavevector of the background medium, given where β is the graphene plasmon wavevector, while P is the Poynting vector, given as For very large β, which is true for graphene plasmons, we find β ε d (ω/c) 2 and hence q d = β. Substituting this into P and η, we find that the SINE factor is reduced to the following expression: From equation (4.2), it is very clear that the SINE factors for graphene plasmons, which scale as n 4 wg , are very much higher than those for metals. Previously, we have estimated the n wg for metals to be in the range of 5-20, thus having SINE factors from 2.5 × 10 3 to 6.4 × 10 5 . Meanwhile, n wg for graphene goes from 100 to 300, which translates to SINE factors from 4 × 10 8 to 3.2 × 10 10 , which are, on average, greater than those for metals by five orders of magnitude. This simple estimation of the SINE factors can be validated by comparison of values between metal and graphene plasmonics, calculated using full vectorial analysis of the NLSE, which are taken from Gorbach [122] and Skryabin et al. [131] and reproduced in figure 13. It can be observed from figure 13a that the various configurations of metallic plasmonic waveguides can only yield SINE factors up to 100, while in figure 13b SINE factors of graphene range from 10 7 to 10 13 . In figure 14, the total effective nonlinear parameters Γ of various graphene-based waveguides are compared [126]. As expected, graphene-on-silicon waveguides shown in figure 14a have the lowest Γ due to the principally dielectric waveguide modes. In figure 14b, graphene-on-metal plasmonic waveguides have much higher Γ , arising from the surface effects of the metal plasmons. Finally, graphene plasmonic waveguides have the highest Γ as a result of its extreme optical confinement and high n wg .

Applications of nonlinear graphene plasmonics (a) Nonlinear generation of graphene plasmons
Graphene plasmons are difficult to excite due to the large mismatch between the momenta of freespace photons and plasmons. Moreover, graphene plasmons exist at the mid-and far-infrared wavelengths, the spectra where only a few excitation sources had been developed. The currently viable excitation method is through optical scattering at AFM tips, as shown in figure 2 [9,10], albeit with low efficiency. Recently, nonlinear optical processes were also studied as a way to generate graphene plasmons. For example, the difference frequency generation (DFG) method makes use of ubiquitous visible light sources and converts them to far-infrared plasmons via graphene's ultrahigh second-order nonlinearity [135,136]. This method uses two pump photons that create a plasmon at their difference frequency, as shown in figure 15. Generation efficiencies of 10 −5 have been recorded [136], and theoretical efficiencies are predicted to reach 1%.
Another proposed method is through the FWM process that exploits the third-order Kerr nonlinearities [137][138][139]. With a similar concept to DFG, the FWM nonlinear process takes two pump photons and converts them into two signal photons through ω GP = 2ω pump − ω signal , in   which one of the signals, ω signal , is at a higher frequency than the pump. The other signal, ω GP , is at a lower frequency, which is easily coupled to the graphene plasmon if the frequency is tuned to the mid-and far-infrared spectrum.

(b) Enhancement of nonlinearities through graphene nanostructures and metamaterials
While graphene plasmons on graphene plasmonic waveguides can already reach optical confinement levels above two orders, patterned graphene nanostructures could push the confinement to even more extreme levels to realize single-photon, quantum-level nonlinear optics. An early theoretical proposal by Gullans et al. [140] designed the graphene plasmon cavity as small as five hexagonal carbon-atom rings. As the size of the graphene is so small, it induces a large dispersive nonlinearity where only a single photon can resonantly excite the graphene plasmon on the cavity. Jadidi et al. [141] studied nonlinear absorption of graphene plasmons and found a nonlinear enhancement of two orders of magnitude when graphene is patterned into nanoribbons. Meanwhile, a quantum dot and graphene nanoflake system was studied by Cox et al. [142] for its configurable two-photon absorption (TPA) efficiency. It was found that, due to its highly resonant structure, the TPA process could be switched on or off depending on the Fermi level of the graphene dot. In later works, Cox and co-workers also studied triangularshaped graphene nanoislands [133,143] and found that the nonlinear response varies with the size, Fermi level as well as the edge structures (armchair or zigzag) of graphene, as shown in figure 16. You et al. [144] studied the nonlinear enhancement of graphene ribbons and discs on the third-harmonic generation efficiency.
On a different note, patterned metal plasmonic nanostructures and metamaterials can also be combined with graphene to enhance nonlinearities. Khorasaninejad et al. [145] experimentally found that Raman scattering intensity is enhanced by a factor of approximately 900 when graphene is put on top of silver plasmonic nanostructures. Smirnova et al. [146] designed gold metamaterial structures placed on top of graphene to excite the asymmetric Fano resonance modes that enhance second-harmonic generation in graphene.

(d) Enhanced second-, third-and high-harmonic generation
In the previous discussion, harmonic generation of photons in graphene can be used to excite graphene plasmons. These harmonic generation methods are, however, weak due to a lack of strong light-matter interaction in thin-film graphene. Proposals to exploit the strong confinement of graphene plasmons would improve the light-matter interaction to increase the efficiencies of second-, third-and even higher-harmonic generation. Mikhailov [132] has studied the effect of graphene plasmons on the second-harmonic generation and found the efficiency to increase by several orders of magnitude. Using plasmonic resonance, Nasari & Abrishamian [149,150] achieved five orders of improvement to the efficiency of third-harmonic generation in graphene.
In hybrid graphene plasmonic waveguides, Smirnova & Solntsev [125] simultaneously phasematched the cascaded third, second and fundamental graphene plasmon harmonics to improve the third-harmonic generation by a factor of five compared with the non-cascaded regime. Finally, Cox et al. [151] studied high-harmonic generation in graphene plasmons and found that harmonics as high as the 13th order could be theoretically generated with optical intensity as low as 100 MW cm −2 . Several other interesting nonlinear phenomena can be seen or enhanced in graphene plasmons. One of them is the plasmonic bistability phenomenon, which takes advantage of the angle sensitivity in the Kretschmann excitation method of graphene plasmons [152][153][154], as shown in figure 18a. The coupling of photons to the graphene plasmons is dependent on resonant phase-matching and therefore dependent on both the incident angle of the photon and the optical conductivity of graphene. When the photons are successfully coupled to plasmons, the plasmon field intensity slowly rises, activating the nonlinear response that alters graphene's optical conductivity, to the point where the photon-plasmon coupling is no longer resonant, as shown in figure 18b. The plasmon fields would therefore drop in intensity to the point where resonant coupling is available again, forming a hysteresis loop.  Nonlinear wave-mixing is enhanced by graphene plasmons in graphene nanoislands [155], which includes sum frequency generation (SFG), DFG and degenerate FWM. The graphene plasmonic nanoislands exhibit nonlinearities two orders of magnitude higher than their metallic counterparts of the same size. These nanoislands are also proposed to be used as nonlinear plasmonic biosensors, as the presence of an individual molecule is sufficient to trigger a large change in the nonlinear response of the graphene plasmons [156]. Meanwhile, plasmonsolitons are supported by graphene waveguides and metamaterial structures [128,[157][158][159]. The characteristics of the plasmon-solitons can be tuned by the input optical intensity as well as graphene's Fermi level. Graphene plasmonic waveguides have also been proposed for supercontinuum generation that could cover a broad range of wavelengths, from 6 to 13 µm, using peak optical powers as low as 0.5 mW [51]. On top of that, plasmonic gain in graphene has also been investigated; however, there have been lacklustre results due to the ultrafast spontaneous emission of the inverted carriers in graphene [160]. Finally, swift electrons can interact strongly with graphene to create nonlinear graphene plasmons [161,162]. It has been theoretically calculated that a single electron moving at 1% of the speed of light is sufficient to generate high-intensity graphene plasmons that alter the plasmon dispersion curve.

Summary and outlook
We have reviewed the recent explosive developments in graphene nonlinear optical properties and graphene plasmonics. Graphene has ultrahigh nonlinear coefficients originating from its peculiar 2D electronic bandstructure. On the other hand, graphene plasmons exhibit extreme optical confinement at long wavelengths due to its unique 2D linear optical conductivity. The synergy between the linear and nonlinear optical properties of graphene plasmons is driving the development of a new field called nonlinear graphene plasmonics, which amplifies further the nonlinear efficiencies of graphene-based devices by over a billion-fold.
It has been shown that the SINE effect, which scales as n 4 wg , augments graphene's nonlinear parameter by up to 10 10 times and beyond, built entirely upon graphene's high plasmonic waveguide indices. This allows all-optical switching devices to be constructed small, and the requirements of optical pump power are reduced to the pW range. Furthermore, when graphene is patterned into nanostructures that enhance the optical confinement ability, further enhancement of various nonlinear optical processes is observed, which includes harmonic generation, nonlinear absorption, SFG and DFG, FWM and plasmonic soliton propagation. On top of that, interesting nonlinear optical phenomena such as optical bistability can only be observed in graphene plasmonics owing to the highly nonlinear resonant coupling of the photon and plasmon.
The birth of nonlinear graphene plasmonics could revolutionize the field of all-optical control, enabling further miniaturization of nonlinear optical components with reduced power requirements and improved switching contrast. The ultra-efficient high-harmonic generation, predicted up to 13 orders, could technically generate photons at any desired frequency. Besides that, the study of supercontinuum generation puts the wavelength span as broad as 7 µm, allowing for chip-scale broadband sources to be built. Graphene as a wonder material in nonlinear plasmonics is expected to play an important role in various optical applications for years to come, including optical computing, biosensing, medical apparatus and security, to name but a few.