Quantum Effects in the Nonlinear Response of Graphene Plasmons

TL;DR

This paper demonstrates that finite-size quantum effects significantly enhance the nonlinear optical response of nanostructured graphene, surpassing classical predictions and highlighting its potential for advanced nanophotonic applications.

Contribution

It provides a quantum-mechanical analysis showing the importance of finite-size effects on nonlinear responses in graphene, which classical models underestimate.

Findings

Finite-size effects increase nonlinear response beyond classical predictions.

Quantum and classical models agree for linear response under certain conditions.

Third-order nonlinearities are generally underestimated by classical theories.

Abstract

The ability of graphene to support long-lived, electrically tunable plasmons that interact strongly with light, combined with its highly nonlinear optical response, has generated great expectations for application of the atomically-thin material to nanophotonic devices. These expectations are mainly reinforced by classical analyses performed using the response derived from extended graphene, neglecting finite-size and nonlocal effects that become important when the carbon layer is structured on the nanometer scale in actual device designs. Here we show that finite-size effects produce large contributions that increase the nonlinear response of nanostructured graphene to significantly higher levels than those predicted by classical theories. We base our analysis on a quantum-mechanical description of graphene using tight-binding electronic states combined with the random-phase approximation. While classical and quantum descriptions agree well for the linear response when either the plasmon energy is below the Fermi energy or the size of the structure exceeds a few tens of nanometers, this is not always the case for the nonlinear response, and in particular, third-order Kerr-type nonlinearities are generally underestimated by the classical theory. Our results reveal the complex quantum nature of the optical response in nanostructured graphene, while further supporting the exceptional potential of this material for nonlinear nanophotonic devices.

Tables

Table 1. Table 1: Numerical Values for the parameters used to calculate the classical polarizabilities.

	${\eta }_1$	${\xi }_1$	${\zeta }_{2\omega }^{\left(2\right)}$	${\zeta }_{3\omega }^{\left(3\right)}$	${\zeta }_{\omega }^{\left(3\right)}$
Triangle	-0.0933	0.541	-1.90	1.57	1.57
Ribbon	-0.0709	0.951	N/A	1.46	1.46

Equations

$${\alpha }_{\omega }^{(1)}(\omega )=\frac{{\eta }_{\omega }^{(1)}{\xi }_1^2{D}^3}{1-{\eta }_{\omega }^{(1)}(\omega )/{\eta }_1},$$

$${\xi }_1=-\left(1/D^2\right){\int }_S{d}^2\mathbf{R}{\epsilon }_{1,x}$$

$${\alpha }_{2\omega }^{(2)}(\omega )$$

$${\alpha }_{3\omega }^{(3)}(\omega )$$

$${\alpha }_{\omega }^{(3)}(\omega )$$

$${\zeta }_{2\omega }^{(2)}=\left(1/D\right){\int }_S{d}^2\mathbf{R}\left[{\epsilon }_{1,x}\left({\partial }_x{\epsilon }_{1,x}+5{\partial }_y{\epsilon }_{1,y}/3\right)+{\epsilon }_{1,y}\left({\partial }_y{\epsilon }_{1,x}/3-{\partial }_x{\epsilon }_{1,y}\right)\right],$$

$${\zeta }_{3\omega }^{(3)}=-\left(1/D^2\right){\int }_S{d}^2\mathbf{R}\left({\vec{\epsilon }}_1\cdot {\vec{\epsilon }}_1\right){\epsilon }_{1,x},$$

$${\zeta }_{\omega }^{(3)}=-\left(1/3D^2\right){\int }_S{d}^2\mathbf{R}\left(2|{\vec{\epsilon }}_1{|}^2{\epsilon }_{1,x}+{\vec{\epsilon }}_1\cdot {\vec{\epsilon }}_1{\epsilon }_{1,x}^{\ast }\right)$$

$$\frac{\partial f_{\mathbf{k}}(\mathbf{R},t)}{\partial t}-\frac{e}{\hslash }\text{E}\cdot {\nabla }_{\mathbf{k}}{f}_{\mathbf{k}}(\mathbf{R},t)\pm v_{\mathrm{F}}\frac{\mathbf{k}}{k}\cdot {\nabla }_{\mathbf{R}}{f}_{\mathbf{k}}(\mathbf{R},t)=-\frac{1}{\tau }\left[f_{\mathbf{k}}(\mathbf{R},t)-f_{\mathbf{k}}^0\right],$$

$${\sigma }_{\omega }^{(1)}(\omega )=\frac{\mathrm{i}{e}^2{E}_F}{\pi {\hslash }^2(\omega +\mathrm{i}{\tau }^{-1})},$$

$$j_{2\omega ,i}^{(2)}(\omega )={\sigma }_{2\omega }^{(2)}(\omega )\sum _{jkl}\left(\frac{5}{3}{\delta }_{ij}{\delta }_{kl}-{\delta }_{ik}{\delta }_{jl}+\frac{1}{3}{\delta }_{il}{\delta }_{jk}\right)E_{\omega ,j}{\partial }_k{E}_{\omega ,l},$$

$${\sigma }_{2\omega }^{(2)}(\omega )=\frac{3\mathrm{i}{e}^3{v}_F^2}{8\pi {\hslash }^2(\omega +\mathrm{i}{\tau }^{-1}{)}^3}$$

$${\sigma }_{3\omega }^{(3)}(\omega )=\frac{3\mathrm{i}{e}^4{v}_F^2}{4\pi {\hslash }^2{E}_F(\omega +\mathrm{i}{\tau }^{-1})(2\omega +\mathrm{i}{\tau }^{-1})(3\omega +\mathrm{i}{\tau }^{-1})},$$

$${\sigma }_{\omega }^{(3)}(\omega )=\frac{9\mathrm{i}{e}^4{v}_{\mathrm{F}}^2}{4\pi {\hslash }^2{E}_{\mathrm{F}}(\omega +\mathrm{i}{\tau }^{-1})(-\omega +\mathrm{i}{\tau }^{-1})(2\omega +\mathrm{i}{\tau }^{-1})}.$$

$$\frac{\partial\rho}{\partial t}=-\frac{{\rm i}}{\hbar}\left[H_{TB}-e\phi,\rho\right]-\frac{1}{2\tau}\big{(}\rho-\rho^{0}\big{)},$$

$$|j,k⟩=\frac{1}{\sqrt{N}}\sum _{l,m}{a}_{jl,k}{\mathrm{e}}^{\mathrm{i}kmb}|l,m⟩,$$

$$\rho =\sum _{j{j}^{\prime },k{k}^{\prime }}{\rho }_{j{j}^{\prime },k{k}^{\prime }}|j,k⟩⟨j^{\prime },k^{\prime }|=\sum _{j{j}^{\prime },k{k}^{\prime },l{l}^{\prime }}{\rho }_{j{j}^{\prime },k{k}^{\prime }}{a}_{jl,k}{a}_{j^{\prime }{l}^{\prime },k^{\prime }}^{\ast }|l,k⟩⟨l^{\prime },k^{\prime }|,$$

$$|l,k⟩=\frac{1}{\sqrt{N}}\sum _m{\mathrm{e}}^{\mathrm{i}kmb}|l,m⟩.$$

$$\frac{\partial {\rho }_{j{j}^{\prime },k{k}^{\prime }}}{\partial t}$$

$$\displaystyle\quad-\frac{1}{2\tau}\big{(}\rho_{jj^{\prime},kk^{\prime}}-\rho^{0}_{jj^{\prime},kk^{\prime}}\big{)},$$

$$\rho =\sum _{n,s}{\rho }^{ns}{\mathrm{e}}^{\mathrm{i}s(qmb-\omega t)}.$$

$${\varphi }_{j{j}^{\prime },k{k}^{\prime }}=\sum _{l,m}{a}_{jl,k}^{\ast }{a}_{j^{\prime }l,k^{\prime }}{\varphi }_l{\mathrm{e}}^{\mathrm{i}\left(q-k+k^{\prime }\right)mb}/N={\delta }_{k-q,k^{\prime }}\sum _l{a}_{jl,k}^{\ast }{a}_{j^{\prime }l,k}{\varphi }_l,$$

$${\varphi }_l={\varphi }_l^{\text{ext}}-2e\sum _{l^{\prime }}{\bar{v}}_{l{l}^{\prime }}{\rho }_{l^{\prime }{l}^{\prime }},$$

$${\rho }_{j{j}^{\prime },k{k}^{\prime }}^{ns}{\mathrm{e}}^{\mathrm{i}sqmb}=-\frac{e}{\hslash }\sum _{l{l}^{\prime }}\frac{{\varphi }_l^{ns}{\rho }_{l{l}^{\prime },k-sq,k^{\prime }}^0-{\varphi }_{l^{\prime }}^{ns}{\rho }_{l{l}^{\prime },k,k^{\prime }+sq}^0}{s\omega +\mathrm{i}{\tau }^{-1}/2-\left({\epsilon }_{j,k}-{\epsilon }_{j^{\prime },k^{\prime }}\right)}{a}_{jl,k}^{\ast }{a}_{j^{\prime }{l}^{\prime },k^{\prime }}+{\eta }_{j{j}^{\prime },k{k}^{\prime }}^{ns},$$

$${\eta }_{j{j}^{\prime },k{k}^{\prime }}^{ns}=-\frac{e}{\hslash }\sum _{n^{\prime }=1}^{n-1}\sum _{s^{\prime }=-n^{\prime }}^{n^{\prime }}\sum _{l{l}^{\prime }}\frac{{\varphi }_l^{n^{\prime }{s}^{\prime }}{\rho }_{l{l}^{\prime },k-s^{\prime }q,k^{\prime }}^{n-n^{\prime },s-s^{\prime }}-{\varphi }_{l^{\prime }}^{n^{\prime }{s}^{\prime }}{\rho }_{l{l}^{\prime },k,k^{\prime }+s^{\prime }q}^{n-n^{\prime },s-s^{\prime }}}{s\omega +\mathrm{i}{\tau }^{-1}/2-\left({\epsilon }_{j,k}-{\epsilon }_{j^{\prime },k^{\prime }}\right)}{a}_{jl,k}^{\ast }{a}_{j^{\prime }{l}^{\prime },k^{\prime }}$$

$${\varphi }_l^{ns}={\varphi }_l^{\text{ext}}{\delta }_{n,1}\left({\delta }_{s,-1}+{\delta }_{s,1}\right)-2e\sum _{l^{\prime },m}{v}_{l{l}^{\prime },0m}{\mathrm{e}}^{\mathrm{i}sqmb}{\rho }_{l^{\prime }{l}^{\prime }}^{ns}.$$

$${\rho }_{ll,q}^{ns}=\frac{-1}{2e}\sum _{l^{\prime }}{\chi }_{l{l}^{\prime },q}^0(s\omega ){\varphi }_{l^{\prime }}^{ns}+{\beta }_{l,q}^{ns}$$

$${\chi }_{l{l}^{\prime },q}^0(\omega )=\frac{2e^2}{\hslash }\frac{b}{2\pi }{\int }_{-\pi /b}^{\pi /b}dk\sum _{j{j}^{\prime }}\left(f_{j^{\prime },k-sq}-f_{j,k}\right)\frac{a_{jl,k}{a}_{j^{\prime }l,k-sq}^{\ast }{a}_{j{l}^{\prime },k}^{\ast }{a}_{j^{\prime }{l}^{\prime },k-sq}}{\omega +\mathrm{i}{\tau }^{-1}/2-\left({\epsilon }_{j,k}-{\epsilon }_{j^{\prime },k-sq}\right)}{\mathrm{e}}^{-\mathrm{i}sqmb}$$

$${\beta }_{l,q}^{ns}=\sum _{j{j}^{\prime },k{k}^{\prime }}{a}_{jl,k}{a}_{j^{\prime }l,k^{\prime }}{\eta }_{j{j}^{\prime },k{k}^{\prime }}^{ns}{\mathrm{e}}^{-\mathrm{i}sqmb}.$$

$${\alpha }_{s\omega }^{(n)}=-\frac{2e}{(E_0{)}^s}\sum _l{\rho }_{ll,0}^{ns}{\text{R}}_l\cdot \hat{\text{e}}.$$

$${\alpha }_{s\omega }^{(n)}(\omega )=\frac{1}{(E_0{)}^n}\int d^2\mathbf{R}\left(\hat{\text{e}}\cdot \mathbf{R}\right){\rho }_{s\omega }^{\text{ind},(n)}(\mathbf{R},\omega ),$$

$${\rho }_{s\omega }^{\text{ind},(n)}(\mathbf{R},\omega )=-\frac{\mathrm{i}}{s\omega }{\nabla }_{\mathbf{R}}\cdot {\mathbf{j}}_{s\omega }^{(n)}(\mathbf{R},\omega ),$$

Full text

I Introduction

Graphene, the two-dimensional monolayer of carbon atoms arranged in a honeycomb lattice, has proven to be an ideal material platform for numerous photonic and opto-electronic applications Bonaccorso et al. (2010); Bao and Loh (2012). Owing to its gapless, linear charge-carrier dispersion relation Wallace (1947); Castro Neto et al. (2009), graphene is capable of relatively strong, broadband coupling with light Nair et al. (2008); Mak et al. (2008), and offers a large electrically-tunable optical response Liu et al. (2011). Valence and conduction electrons in this material present a uniform velocity that clearly emphasizes their anharmonic response to external fields, which has stimulated considerable interest in the graphene nonlinear optical properties Mikhailov (2007); Hendry et al. (2010); Dean and van Driel (2010); Wu et al. (2011); Zhang et al. (2012); Gu et al. (2012); Kumar et al. (2013); Hong et al. (2013); Glazov and Ganichev (2014); Peres et al. (2014); Yao et al. (2014); Politano and Chiarello (2014); Nasari and Abrishamian (2014); Cheng et al. (2014, 2015); M (15). Further motivation to study optical nonlinearities in monolayer graphene is provided by the availability of interband optical transitions within a continuous range of low photon energies, accompanied by a high electrical mobility Novoselov et al. (2004, 2005). The large third-order susceptibilities corresponding to four-wave mixing Hendry et al. (2010); Gu et al. (2012), the optical Kerr effect Wu et al. (2011); Zhang et al. (2012), and third-harmonic generation Kumar et al. (2013); Hong et al. (2013) (THG) that have been measured in graphene are indeed attributed to these electronic properties. However, the superiority of graphene with respect to conventional nonlinear optical semiconductors for ultrafast all-optical switching is not yet fully established Khurgin (2014), although this atomically thin material has been shown to be an excellent broadband saturable absorber for laser mode locking Bao and Loh (2012).

The intrinsically-high nonlinear optical response of graphene can be further enhanced by plasmonic excitations supported by the carbon layer when it is electrically doped. Unlike traditional plasmons in noble metals, whose resonance frequencies are determined by geometric and intrinsic properties, graphene plasmons can be actively tuned by modifying the doping density Ju et al. (2011); Fei et al. (2011); Chen et al. (2012); Fei et al. (2012); Yan et al. (2012a, b); Fang et al. (2013); Brar et al. (2013); Fang et al. (2014); Yan et al. (2013); Freitag et al. (2014); Grigorenko et al. (2012); García de Abajo (2014). Because of its superior electrical conductivity, plasmons in graphene are also longer-lived and couple more efficiently to light than their noble-metal counterparts. Additionally, they exhibit extreme subwavelength confinement Jablan et al. (2009); Koppens et al. (2011), so that the excitation of propagating plasmons in extended graphene often relies on near-field coupling to satisfy energy-momentum conservation Fei et al. (2011); Chen et al. (2012); Fei et al. (2012). This constraint is relaxed in nanostructured graphene, enabling optical excitation of localized plasmons from free space Yan et al. (2012c); Fang et al. (2013); Brar et al. (2013); Yan et al. (2013). The strong near electric fields generated by plasmons in noble metals have been widely exploited to enhance nonlinear optical processes Danckwerts and Novotny (2007); Kauranen and Zayats (2012); Butet et al. (2015), and thus, even larger enhancements are expected from graphene plasmons, as recently predicted by several theoretical studies Mikhailov (2011); Gullans et al. (2013); Gorbach (2013); Scholz et al. (2013); Nesterov et al. (2013); Yao et al. (2014); Cox and García de Abajo (2014, 2015); Smirnova et al. (2015); Jablan and Chang (2015); Christensen et al. (2015); Manzoni et al. (2015).

For extended graphene, the nonlinear response is commonly described using a classical nonlinear conductivity derived from the Boltzmann transport equation (BTE) Mikhailov (2007); Peres et al. (2014) assuming that intraband transitions dominate the optical response (i.e., for photon energies $\hbar\omega$ roughly below the Fermi energy ${E_{\rm F}}$). At higher energies, a quantum-mechanical (QM) treatment is required to account for the linear and nonlinear response arising from interband electron transitions Cheng et al. (2014, 2015); M (15). This approach has been used to predict a large nonlinear optical susceptibility in extended graphene at low doping (i.e., without plasmonic enhancement), for which reasonable agreement with experimental observations has been reported Hendry et al. (2010).

To describe plasmon-enhanced nonlinear optical processes in nanostructured graphene, a finite graphene structure is typically assigned a local nonlinear conductivity, and its response simulated within the framework of classical electrodynamics Gullans et al. (2013); Smirnova et al. (2015); Christensen et al. (2015); Manzoni et al. (2015). However, this classical approach neglects nonlocal, finite-size, and atomistic (e.g., edge termination) effects, all of which have been found to play crucial roles in describing both the linear Thongrattanasiri et al. (2012); Christensen et al. (2014); Wang et al. (2015) and nonlinear Cox and García de Abajo (2014, 2015) optical response of graphene structures with small ($\sim 10$ nm) features. Here we reveal strong finite-size and atomistic effects in the nonlinear response of graphene nanoribbons and nanoislands, predicted to take place from a realistic, QM description of these structures, beyond classical theory, which can produce dramatically different results for structure sizes up to a few tens of nanometers. The discrepancy between classical and QM descriptions is particularly large for the Kerr effect/nonlinear absorption (i.e., the complex third-order susceptibility oscillating at the fundamental frequency, which from here onwards will be referred to simply as the “Kerr nonlinearity”) where classical theory underestimates the strength of the third-order response by several orders of magnitude even for $>20\,$nm structures. This is at odds with the conclusions previously drawn by examining the linear response regime, where the plasmon frequencies and strengths were found to be similar within classical and QM approaches either when the structures were a few tens of nanometers in size or when they did not possess zigzag edges. The QM effects here reported for such large sizes represent a rather unusual scenario in plasmonics, while they support the use of doped nanographene structures as plasmon-driven nonlinear enhancers, where they perform much better than previously estimated from the study of extended graphene.

II Results and discussion

II.1 Graphene nanoribbons

We study graphene nanoribbons illuminated with light polarized perpendicular to their direction of translational symmetry. As second-order processes are forbidden in this geometry Boyd (2008), we concentrate on first- and third-order phenomena. In Fig. 1 we examine the incident-frequency dependence of the linear and nonlinear polarizabilities for ribbons roughly 10 nm in width. In particular, we present atomistic QM simulations for ribbons of armchair (carbon-to-carbon width 9.97 nm) and zigzag (width 9.81 nm) edges, compared with classical theory (width 10 nm). Details of the calculation procedures are given in the Appendix. The linear response (Fig. 1b) is characterized by a prominent plasmon feature, for which the classical calculations are in reasonable agreement with the quantum model for both types of edges, provided that the Fermi energy ${E_{\rm F}}$ is larger than the plasmon energy. The only discrepancy in this regime is a small systematic plasmon redshift and broadening in the QM model relative to the classical simulations. In contrast, when the plasmon energy is above ${E_{\rm F}}$, the quantum model predicts a strong plasmon attenuation and the shift becomes dramatic.

Remarkably, similar agreement between classical and quantum models is observed in THG at high doping (Fig. 1c), which deteriorates as the Fermi energy is lowered, again producing an overestimate of the classical response compared with the QM simulations. This is unlike the Kerr nonlinearity (Fig. 1d), which we find to be substantially smaller in the classical description, a result that could originate in the inability of this model to account for the simultaneous plasmonic enhancement of the frequency-degenerate input and output waves. The values of $\chi^{(3)}(\omega)$ predicted in the QM description for graphene nanoribbons are on the order of $10^{-5}$ esu ($10^{-12}$ m2/V2), to be compared with those measured in gold, ranging from $10^{-7}$–$10^{-13}$ esu ($10^{-14}$–$10^{-20}$ m2/V2) Boyd et al. (2014) (values for the nonlinear refractive index/nonlinear absorption coefficient are provided in Supporting Information, SI). We also note that these values depend strongly on the inelastic scattering decay time, $\tau$, for which we have used a highly conservative value $\hbar\tau^{-1}=50$ meV, or $\tau\simeq 13$ fs, corresponding to moderate values of the DC mobility less than 660 cm${}^{2}/($V s$)$, as estimated from the Drude model Ashcroft and Mermin (1976) for the Fermi energies under consideration.

An overview of the size and doping dependence of the linear and nonlinear ribbon polarizabilities is presented in Fig. 2. Again, the linear absorption of armchair and zigzag ribbons in the quantum description displays very similar plasmon frequencies and strengths when the plasmon energy is below ${E_{\rm F}}$ (i.e., for widths $\sim 5-15\,$nm, depending on the actual value of ${E_{\rm F}}$ under consideration), and they in turn agree with the classical simulations. For larger plasmon energies (smaller sizes) quantum effects become important, reflecting in particular a strong coupling to electronic edge states in zigzag ribbons Silveiro et al. (2015). We extract similar conclusions regarding the ability of classical theory to describe THG (Fig. 2), in good qualitative agreement with QM simulations, except at low doping levels and small ribbon widths. However, classical theory underestimates the polarizability associated with the Kerr nonlinearity by roughly an order of magnitude, even for ribbon widths as large as 40 nm.

II.2 Graphene nanoislands

In Fig. 3, we examine the spectral response of the linear and nonlinear polarizability for equilateral triangular nanographenes roughly 10 nm in side length and illuminated by light polarized in a direction perpendicular to one of the triangle sides (see Fig. 3 top). Despite the centrosymmetry of the graphene crystal structure, even-ordered, dipolar nonlinear processes are enabled by the symmetry-breaking of finite structures lacking inversion symmetry along the direction of the induced nonlinear dipole moment. The chosen triangular shape allows us to study islands with approximately equal dimensions, but containing exclusively either armchair- or zigzag-terminated edges (see Fig. 3a, bottom). The linear absorption cross-sections, normalized to the nanoisland areas, are presented in Fig. 3b,c, where, for lower doping, we find dramatic differences between atomistic and classical results, and also strong discrepancies among the two types of edge terminations considered in the atomistic simulations. With increasing Fermi level, the numerous plasmon resonances that appear due to the finite-size and nonlocal effects captured by the QM description tend to coalesce towards the prominent dipolar plasmon mode predicted by the classical simulations. We also note that the spectral features of the zigzag island are typically smaller in magnitude than those of the armchair island, which is attributed to the presence of zero-energy electronic states confined to the zigzag edges Manjavacas et al. (2013); Christensen et al. (2014).

Our results employing the classical and quantum descriptions of the second-order nonlinear polarizability corresponding to second-harmonic generation (SHG) (Fig. 3d,e) reveal that the classical theory provides good qualitative agreement with QM simulations, in particular for armchair-edged islands, when plasmon resonances are at energies below $E_{F}$. The classical description relies on the second-order conductivity for SHG in extended graphene, for which the leading contribution is nonlocal due to photon momentum transfer Manzoni et al. (2015). Conversely, the atomistic simulations are not limited by any such assumed symmetry, and also account for the responses at both the fundamental and second harmonic frequencies, which can significantly enhance SHG at low doping levels Cox and García de Abajo (2014). This is similar to surface-dipole-induced SHG, which is known to be generally stronger than bulk quadrupole Dean and van Driel (2010); Butet et al. (2015). Nevertheless, our findings suggest that the nonlocal electric field enhancement associated with a plasmon resonance, which is well described by our classical model when $E_{F}$ is larger than the resonance frequency, provides the leading contribution to SHG in nanoislands.

The third-order nonlinear polarizabilities associated with THG and the Kerr nonlinearity are presented in Fig. 3f,g and Fig. 3h,i, respectively, for the triangular nanoisland. We find fairly good agreement between the THG polarizabilities calculated using the classical and quantum approaches at high doping, while for the Kerr nonlinearity the classical treatment underestimates the polarizability by one order of magnitude compared with the quantum approach, except for very low doping levels. Conversely, the Kerr nonlinearity involves the mixing of three waves at the fundamental frequency to produce a fourth wave, also oscillating at the fundamental frequency. At the dominant plasmon mode, all four waves are then simultaneously enhanced by the near electric field of the plasmon. This enhancement at the output frequency is not taken into account in the classical description, resulting in a significantly weaker Kerr response.

Figure 4 presents an overview of our results for nanotriangles. As expected, the linear response shows that the dipolar plasmon resonances predicted in the classical treatment generally agree more closely with the atomistic simulations for armchair-edged nanoislands. For SHG and THG, the nonlinear polarizabilities obtained from the classical treatment are in good agreement with the atomistic simulations for the armchair islands, presumably for similar reasons as in the linear response, while the quantum simulations predict a much larger nonlinear response for the Kerr nonlinearity, which we attribute to the simultaneous enhancement at the fundamental and mixed frequencies –an effect that is not captured in the classical description.

III Conclusions

We have demonstrated that the nonlinear optical response of nanostructured graphene is strongly influenced by nonlocal and finite-size effects, which cannot easily be ignored in structures with dimensions of $<20$ nm, and quite possibly in even larger systems. The classical electrodynamic description fails to capture these effects, which we have found to be particularly important when dealing with the optical Kerr nonlinearity. Interestingly, THG in nanoribbons is very well described by the classical approach, which also provides reasonable agreement with quantum-mechanical results for graphene nanoislands, though better agreement is obtained for predominantly armchair-edged islands. Qualitatively similar agreement between quantum-mechanical and classical descriptions is also found for SHG in nanoislands lacking inversion symmetry, indicating that the second-order nonlocal bulk nonlinear conductivity, combined with the dipolar electric field associated with a plasmon mode, is sufficient to describe SHG in such structures. Barring edge-effects in nanoislands, it is then apparent that SHG and THG are reasonably well described by the classical electrodynamic simulations when plasmon energies are much less than the chemical potential gap $2{E_{\rm F}}$, yet the optical Kerr nonlinearity is lacking significant contributions from the nonlocal response at the fundamental frequency, for which four waves are involved; in this case, a more fundamental quantum-mechanical description predicts a much higher nonlinear response than classical theory, thus providing further support for the use of highly doped graphene as an excellent nonlinear element in nanophotonic devices.

Appendix A Classical electrodynamical simulations

Our classical calculations are based upon the linear and nonlinear conductivities for extended graphene $\sigma^{(1)}_{\omega}$, $\sigma^{(2)}_{2\omega}$, $\sigma^{(3)}_{3\omega}$, and $\sigma^{(3)}_{\omega}$, which describe linear, second-harmonic, third-harmonic, and Kerr-type processes, respectively. As explained in more detail in the SI, we adopt an eigenmode expansion of the graphene optical response García de Abajo (2013) and concentrate on a spectral region dominated by a single plasmon, whose in-plane electric field is $\vec{\varepsilon}_{1}({\bf R})$. Assuming a local dielectric function $\epsilon=1+4\pi{\rm i}\sigma^{(1)}_{\omega}(\omega)/\omega t$, where $t$ the thickness of the graphene layer (here we take $t=0.5$ nm, which is found to be well converged with respect to the $t\rightarrow 0$ limit), we calculate the field $\vec{\varepsilon}_{1}({\bf R})$ for equilateral triangles using a finite-element method, and normalize it in such a way that $\int_{S}d^{2}{\bf R}|\vec{\varepsilon}_{1}({\bf R})|^{2}=D^{2}$, where $S$ denotes the 2D region occupied by the graphene structure and $D$ is a characteristic length (e.g., the triangle side length). We now define $\eta^{(n)}_{s\omega}={\rm i}\sigma^{(n)}_{s\omega}/s\omega D$, and consider, without loss of generality, a plasmon polarized along a symmetry direction $x$. Then, the linear polarizability is given by

[TABLE]

where

[TABLE]

and $\eta_{1}$ is a plasmon eigenvalue García de Abajo (2013), which is determined from the calculated linear extinction spectrum. Similarly, the nonlinear optical polarizabilities associated with SHG, THG, and the Kerr nonlinearity are given as

[TABLE]

respectively, where

[TABLE]

[TABLE]

and

[TABLE]

are dimensionless coefficients, for which numerical values are provided in Table 1. Nonlinear susceptibilities are estimated by normalizing the polarizabilities to the nanostructure volume, assuming an effective thickness of $0.33$ nm for the graphene layer (i.e., the interplane distance in graphite) Hendry et al. (2010); Cox and García de Abajo (2014).

For one-dimensional graphene nanoribbons, we use the boundary-element method García de Abajo and Howie (2002) to calculate the in-plane electric field $\vec{\varepsilon}_{1}({\bf R})$, taking $D$ as the ribbon width. Following the same procedure as for finite structures, we recover Eqs. (1–4) with the left-hand side of each equation divided by the (infinite) ribbon length, so that $\alpha^{(n)}_{s\omega}$ instead denotes the polarizability per unit length along the ribbon for an $n^{\text{th}}$-order process at harmonic $s$ of the illumination frequency.

Appendix B Nonlinear conductivities in extended graphene

The linear and nonlinear surface currents in graphene are obtained via iterative solution of the Boltzmann transport equation Peres et al. (2014),

[TABLE]

where $f_{\bf k}^{0}$ is the Fermi-Dirac distribution, to which the system relaxes at the phenomenological scattering rate $\tau^{-1}$, and we consider the interaction with an ac field ${\bf E}(t)={\bf E}_{\omega}({\rm e}^{-{\rm i}\omega t}+{\rm e}^{{\rm i}\omega t})$ that is polarized parallel to the graphene plane with polarization across the ribbon (i.e., the $x$-$y$ plane). At first-order, in the long-wavelength limit, the linear surface current is ${\bf j}^{(1)}_{\omega}(\omega)=\sigma^{(1)}_{\omega}(\omega){\bf E}_{\omega}$, with the conductivity

[TABLE]

which coincides with the Drude model, also resulting from the local limit of the random-phase approximation (RPA) after neglecting interband transitions.

The dominant contribution to the second-harmonic current ${\bf j}^{(2)}_{2\omega}(\omega)$ arises from the lowest-order expansion in the photon momentum (i.e., a nonlocal contribution originating in the $\nabla_{\bf R}$ term of Eq. (52)), which leads to Manzoni et al. (2015)

[TABLE]

where the scalar part of the second-harmonic conductivity tensor is given by

[TABLE]

and the $i$ subindex in Eq. (7) refers to the direction of the second-harmonic current.

The leading contribution to the third-harmonic and Kerr-type nonlinear currents is local (i.e., we can neglect the spatial variations in $f_{\bf k}({\bf R},t)$). For THG, the nonlinear surface current is ${\bf j}^{(3)}_{3\omega}(\omega)=\sigma^{(3)}_{3\omega}(\omega)\left({\bf E}_{\omega}\cdot{\bf E}_{\omega}\right){\bf E}_{\omega}$, where the conductivity reads

[TABLE]

while for the Kerr nonlinearity we have ${\bf j}^{(3)}_{\omega}(\omega)=\sigma^{(3)}_{3\omega}(\omega)\left[2|{\bf E}_{\omega}|^{2}{\bf E}_{\omega}+({\bf E}_{\omega}\cdot{\bf E}_{\omega}){\bf E}_{\omega}^{*}\right]/3$, with

[TABLE]

Note that the above expressions for the linear and nonlinear conductivities are obtained in the zero-temperature limit, which should be valid for $T=300$ K provided $\hbar\omega\gg k_{\rm B}T$ (see SI for a comparison of results obtained at $T=0$ and 300 K).

Appendix C Quantum-mechanical simulations

We employ a tight-binding Hamiltonian $H_{TB}$ to describe the $\pi$-band electronic structure of graphene nanoislands and nanoribbons, where each carbon atom is represented by a single $p$ orbital oriented perpendicular to the graphene plane, and a hopping energy of 2.8 eV connects nearest-neighbor carbon sites. To simulate the optical response of these systems we solve the equation of motion for the single-electron density matrix,

[TABLE]

where $\phi$ is the self-consistent electric potential. The last term in Eq. (11) describes inelastic electron scattering with a phenomenological decay time $\tau$, taking $\rho^{0}$ as the $t\rightarrow-\infty$ equilibrium density matrix. The polarizability spectra of Figs. 1 and 3 are obtained with $\hbar\tau^{-1}=50$ meV, while we take $\hbar\tau^{-1}=20$ meV for the polarizability maps of Figs. 2 and 4. These damping rates correspond to moderate values of the DC mobility up to 1646 cm${}^{2}/($V s$)$, as estimated from the Drude model Ashcroft and Mermin (1976) for the Fermi energies under consideration.

For a graphene nanoribbon consisting of $N$ unit cells with period $b$ along its direction of translational symmetry, electronic Bloch states are constructed as

[TABLE]

where $j$ denotes the band index, $k$ is the in-plane Bloch wave vector along the ribbon, and $l$ labels carbon sites $\textbf{R}_{l}$ within each unit cell $m$. We then write the density matrix as

[TABLE]

where we have defined

[TABLE]

Eqs. (12) and (14) yield the relations $\rho_{ll^{\prime},kk^{\prime}}=\sum_{jj^{\prime}}\rho_{jj^{\prime},kk^{\prime}}a_{jl,k}a^{*}_{j^{\prime}l^{\prime},k^{\prime}}$ and $\rho_{jj^{\prime},kk^{\prime}}=\sum_{ll^{\prime}}\rho_{ll^{\prime},kk^{\prime}}a^{*}_{jl,k}a_{j^{\prime}l^{\prime},k^{\prime}}$, which connect site and state representations, allowing us to express the equation of motion in the latter as

[TABLE]

where we have used $H_{TB}|j,k\rangle=\hbar\varepsilon_{j,k}|j,k\rangle$. The equilibrium density matrix $\rho^{0}_{jj^{\prime},kk^{\prime}}=f_{j,k}\delta_{jj^{\prime}}\delta_{kk^{\prime}}$ is constructed by filling electron states according to the Fermi-Dirac distribution occupation numbers $f_{j,k}$ Hedin and Lundqvist (1970).

We consider excitation of a graphene nanoribbon by continuous-wave (cw) illumination of frequency $\omega$ and in-plane wave vector $q$. An iterative solution of the above equation of motion is now facilitated by expanding the density matrix as

[TABLE]

In general, the external potential is a function of the position along the nanoribbon direction, according to

[TABLE]

where

[TABLE]

$\phi^{\text{ext}}_{l}=-\textbf{R}_{l}\cdot\textbf{E}(t)$ is the external potential describing the incident field $\textbf{E}(t)$, $\bar{v}_{ll^{\prime}}=\sum_{m}v_{l0,l^{\prime}m}$, and $v_{l0,l^{\prime}m}{\rm e}^{iqmb}$ is the Coulomb interaction between atoms $l$ and $l^{\prime}$ separated by $m$ unit cells.

Inserting Eqs. (16) and (17) into Eq. (15), and identifying terms with the same ${\rm e}^{-{\rm i}s\omega t}$ dependence on both sides of the equation, we find

[TABLE]

where

[TABLE]

and

[TABLE]

Moving from $|j,k\rangle$ (state) to $|l,m\rangle$ (site) representation, we obtain the diagonal density matrix elements (within a unit cell of the nanoribbon, i.e. for a single $m$) as

[TABLE]

where

[TABLE]

is the noninteracting RPA susceptibility, and

[TABLE]

Here we focus on normally incident light relative to the graphene plane, and thus the external potential does not depend on the position along the nanoribbon direction, i.e., we take $q=0$ in the above expressions. Following a procedure described elsewhere Cox and García de Abajo (2014), we compute the density matrices $\rho^{ns}_{ll^{\prime},0}$, from which we obtain the $n$th-order induced charge at site $l$ oscillating with harmonic $s$ as $\rho^{\text{ind}}_{l}=-2e\rho^{ns}_{ll,0}$. The polarizability per unit length at order $n$ and harmonic $s$ is then given by

[TABLE]

Using this formalism, we compute the linear polarizability $\alpha^{(1)}_{\omega}$, as well as the nonlinear polarizabilities corresponding to SHG $\alpha^{(2)}_{2\omega}$, THG $\alpha^{(3)}_{3\omega}$, and the optical Kerr nonlinearity $\alpha^{(3)}_{\omega}$.

For finite graphene nanoislands, the linear and nonlinear optical response is simulated following the procedure previously described for the analysis of SHG and THG Cox and García de Abajo (2014), employing a fast Fourier transform (FFT) method to expedite the calculation of the noninteracting RPA susceptibilities Thongrattanasiri et al. (2012) (see SI for details on the convergence of the FFT method with direct calculations).

Appendix D Classical electrostatic theory for nanostructured graphene

We consider a graphene nanostructure with a characteristic size $D$, corresponding to the width of a nanoribbon or the side length of an equilateral triangle that is much less than the wavelength of the incident illumination, for which the electric field is given by ${\bf E}^{\text{ext}}=E_{0}\hat{\textbf{e}}\left({\rm e}^{-{\rm i}\omega t}+\text{c.c.}\right)$, where $\hat{\textbf{e}}$ is the polarization unit vector. We quantify the linear and nonlinear optical response of a finite nanostructure by the $n^{\text{th}}$-order dipole induced along $\hat{\textbf{e}}$ oscillating at harmonic $s$ of the excitation frequency $\omega$, for which the polarizability is given by

[TABLE]

where the induced charge density is obtained from the surface currents ${\bf j}^{(n)}_{s\omega}({\bf R},\omega)$ using the continuity equation,

[TABLE]

and ${\bf R}=(x,y)$ are 2-D coordinate vectors in the $x$-$y$ plane.

For the linear optical response (taking $n=s=1$ in the above expressions), the surface current is ${\bf j}^{(1)}_{\omega}({\bf R},\omega)=\sigma^{(1)}_{\omega}({\bf R},\omega){\bf E}({\bf R},\omega)$, where $\sigma^{(1)}_{\omega}({\bf R},\omega)$ is the linear conductivity of graphene (see Appendix) and ${\bf E}({\bf R},\omega)=-\nabla_{{\bf R}}\phi({\bf R},\omega)$ is the total electric field acting on the graphene nanostructure due to the self-consistent potential $\phi$, which is given in the electrostatic approximation by García de Abajo (2013)

[TABLE]

Following the method of Ref. García de Abajo (2013), we assume that the linear conductivity can be separated as $\sigma^{(1)}_{\omega}({\bf R},\omega)=f({\bf R})\sigma^{(1)}_{\omega}(\omega)$, where the occupation factor $f({\bf R})=1$ within the graphene structure and is zero everywhere else, and we express Eq. (27) in terms of a reduced 2-D coordinate vector $\vec{\theta}={\bf R}/D$ as

[TABLE]

In obtaining the above expression we have defined $\eta^{(n)}_{s\omega}(\omega)={\rm i}\sigma^{(n)}_{s\omega}(\omega)/s\omega D$ and introduced the normalized electric field $\vec{\varepsilon}(\vec{\theta},\omega)=-\sqrt{f(\vec{\theta})}\nabla_{\vec{\theta}}\phi(\vec{\theta},\omega)$, which is expanded in a complete set of eigenmodes $\vec{\varepsilon}_{j}$ with real eigenvalues $1/\eta_{j}$ as García de Abajo (2013)

[TABLE]

where the expansion coefficients are

[TABLE]

and the eigenmodes are orthogonal, i.e.,

[TABLE]

Now, assuming that the optical response is dominated by the lowest-order dipolar mode (the $j=1$ term in Eq. (30)), the induced charge density given in Eq. (29) becomes

[TABLE]

We then express Eq. (26) in terms of normalized coordinates $\vec{\theta}$ and use Eq. (33) to write the linear polarizability as

[TABLE]

where we have defined

[TABLE]

Moving back to ${\bf R}$ space, we write $\xi_{1}=-\hat{\textbf{e}}\cdot\int_{S}d^{2}{\bf R}\,\vec{\varepsilon}_{1}({\bf R})/D^{2}$, as given above in the Appendix, by using the function $f({\bf R})$ to restrict the integration to the surface of the graphene nanostructure, $S$, while Eq. (32) guarantees that $\int_{S}d^{2}{\bf R}|\vec{\varepsilon}_{1}({\bf R})|^{2}=D^{2}$.

To describe the nonlinear response associated with second-harmonic generation, we express the second-harmonic current as Manzoni et al. (2015)

[TABLE]

where we assume that the nonlinear conductivity can be written as $\sigma^{(2)}_{2\omega}({\bf R},\omega)=f({\bf R})\sigma^{(2)}_{2\omega}(\omega)$, and we have isolated its tensorial part, ${\bf{\Delta}}^{(2)}_{ijkl}=5\delta_{ij}\delta_{kl}/3-\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}/3$. We then express Eq. (27), for $n=s=2$, in reduced coordinates as

[TABLE]

Keeping the total electric field to linear order, we again use only the $j=1$ term in Eq. (30) and write

[TABLE]

Using the above expression in Eq. (26), expressed in terms of $\vec{\theta}$, then yields the nonlinear polarizability for second-harmonic generation,

[TABLE]

where we introduce the unitless parameter

[TABLE]

Note that $\varepsilon_{1,i}(\vec{\theta})/\sqrt{f(\vec{\theta})}$ is proportional to the physical electric field associated with the dipolar plasmon mode, and therefore is continuous across the graphene edge, even in the limit where $f(\vec{\theta})$ jumps from 1 to 0 at the edge itself; the leading factor of $\sqrt{f(\vec{\theta})}$ then limits the integration in the above expression to the region occupied by the graphene, $S$. By moving to ${\bf R}$ space, and considering (without loss of generality) the response along the $x$-direction, i.e., for $\hat{\textbf{e}}_{x}$, we may express $\zeta^{(2)}_{2\omega}$ in the form

[TABLE]

which corresponds to the expression provided in the Appendix.

The extension to third-order nonlinearities follows straightforwardly: For third-harmonic generation, we start with the third-order surface current

[TABLE]

where the tensor part is ${\bf{\Delta}}^{(3)}_{ijkl}=\left(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)/3$. We isolate the spatial dependence of the third-order conductivity according to $\sigma^{(3)}_{3\omega}({\bf R},\omega)=f({\bf R})^{2}\sigma^{(3)}_{3\omega}(\omega)$, and so Eq. (27) (for $n=s=3$) becomes

[TABLE]

Using Eq. (30), we find that

[TABLE]

from which Eq. (26) provides the nonlinear polarizability corresponding to third-harmonic generation,

[TABLE]

where

[TABLE]

The above expression for the unitless parameter $\zeta^{(3)}_{3\omega}$ can be simplified by moving to ${\bf R}$ coordinates and using the fact that $f({\bf R})=1$ for ${\bf R}\in S$, and is zero otherwise, to write

[TABLE]

corresponding to the expression provided in the Appendix.

In an analogous manner, the nonlinear polarizability corresponding to the Kerr nonlinearity (Eq. (26) for $n=3$, $s=1$) is obtained from the third-order surface current

[TABLE]

as

[TABLE]

where

[TABLE]

After performing the summation in the above expression, we obtain

[TABLE]

which coincides with the expression provided in the Appendix.

Appendix E Temperature-dependence of the nonlinear response

Through iterative solution of the Boltzmann transport equation (in the local limit),

[TABLE]

we take ${\bf E}(t)={\bf E}_{\omega}({\rm e}^{-{\rm i}\omega t}+{\rm e}^{{\rm i}\omega t})$ and find that, to first-order,

[TABLE]

with

[TABLE]

so that Eq. (53) leads to Eq. (5) in the $T\rightarrow 0$ limit. Actually, the linear result is virtually unchanged by considering $T=300$ K, for which we find $F^{(1)}_{\bf k}\simeq E_{F}$ for the Fermi energies considered in this work, i.e. for $0.2\leq E_{F}/\text{eV}\leq 2.0$. At third-order, we obtain expressions for the local intraband conductivity of graphene as

[TABLE]

where

[TABLE]

and $D_{3\omega}(\omega)=(\omega+{\rm i}\tau^{-1})(2\omega+{\rm i}\tau^{-1})(3\omega+{\rm i}\tau^{-1})$ for third-harmonic generation ($s=3$) or $D_{\omega}(\omega)=(\omega+{\rm i}\tau^{-1})(2\omega+{\rm i}\tau^{-1})(-\omega+{\rm i}\tau^{-1})$ for the Kerr nonlinearity ($s=1$). Now, at zero temperature, we have $\partial f^{0}_{\bf k}/\partial\varepsilon_{k}=-\delta\left(E_{F}-\varepsilon_{k}\right)$, and Eq. (55) reduces to $F_{\bf k}=-3\pi/4$, after having used the general relation

[TABLE]

Thus, we recover

[TABLE]

and

[TABLE]

corresponding to Eqs. (9) and (10) in the Appendix, and in agreement with the result of Ref. Peres et al. (2014). These third-order expressions differ in multiplicative factors from those reported in a previous study M (15), where the phenomenological decay is introduced in the Fourier integrals; if in Eq. (52) we had instead defined ${\bf E}(t)=\int^{\infty}_{-\infty}d\omega{\bf E}_{\omega}{\rm e}^{-i\omega t+t/\tau}$, we would recover the purely intraband contributions to the third-order currents obtained from Ref. M (15):

[TABLE]

and

[TABLE]

For the intraband third-order nonlinear conductivities at non-zero temperatures, we find that Eq. (55) reduces to

[TABLE]

Note that the integrand in the first term above has a singularity at $\varepsilon_{k}=0$, and so $F^{(3)}_{\bf k}$ diverges. This suggests that the perturbation theory (up to third order) is inadequate for dealing with finite temperatures, in the same way it fails as $E_{F}\rightarrow 0$ (consider the $1/E_{F}$ dependence of $\sigma^{(3)}_{s\omega}$). However, we expect that the $T=0$ description of the third-order, purely intraband conductivities should also describe the $T=300$ K case reasonably well, considering the negligible change in the linear response (see Eq. (53)), as well as those in the nonlinear response for graphene nanoribbons predicted by atomistic simulations (see Fig. 5).

Appendix F Regarding interband contributions to the nonlinear optical response

Computing the full third-order nonlinear optical conductivity of extended graphene, i.e., including contributions from both intraband and interband optical transitions, is presently an area of active study, with various approaches in the literature that apparently yield different results (see for example the work of Cheng et al. in Refs. Cheng et al. (2014, 2015) and by Mikhailov in Ref. M (15)). In this work, we are primarily interested in the nonlinear optical response of graphene nanostructures enhanced by plasmonic exctations, and thus we restrict our investigation to the regime where only intraband contributions will contribute significantly to the nonlinear conductivities. Incidentally, it has been argued that interband transitions only become important at energies $\hbar\omega\geq 2E_{F}$, and even then their effect is orders of magnitude smaller than the contribution of intraband transitions M (15).

A relatively simple, straightforward improvement to the classical description of the plasmon-enhanced nonlinear response in nanostructured graphene is to include the effect of interband optical transitions in the linear response, which can be accomplished by using the linear conductivity for graphene obtained from the random-phase approximation in the local limit García de Abajo (2014),

[TABLE]

The linear conductivity enters the factors $\eta^{(1)}_{\omega}$ in Eqs. (1-4) of the Appendix, effectively describing the local field generated by the graphene nanostructure with greater accuracy. An accurate description of the local field enhancement is arguably more important when describing the nonlinear response at plasmon resonances than to include the effect of interband optical transitions in the nonlinear conductivities. Indeed, as we show in Fig. 6 for graphene nanoribbons, the peak intensities for the linear and nonlinear responses are in better agreement with those predicted in quantum-mechanical simulations when we use Eq. (60) to describe the linear conductivity in graphene, where the most noticeable improvements appear for resonances above or near the Fermi level.

Appendix G Nonlinear refractive index and nonlinear absorption

The complex third-order susceptibility oscillating at the fundamental frequency of illumination, $\chi^{(3)}_{\omega}$, which we have referred to here as the Kerr nonlinearity, is related to the nonlinear refractive index $n_{2}$ and nonlinear absorption coefficient $\beta$ of a medium. These quantities are typically what are measured directly in nonlinear optical experiments, from which values for the nonlinear susceptibility are inferred. The relations among a medium’s $\chi^{(3)}_{\omega}$, $n_{2}$, and $\beta$ are provided in Ref. del Coso and Solis (2004) as

[TABLE]

and

[TABLE]

where $n_{0}=\text{Re}\{(1+\chi^{(1)}_{\omega})^{1/2}\}$ and $k_{0}=\text{Im}\{(1+\chi^{(1)}_{\omega})^{1/2}\}$, with all of the above quantities in SI units. The above expressions are intended to properly account for the interplay between the complex first- and third-order susceptibilities in highly-absorbing media del Coso and Solis (2004), and we use them here to show the frequency-dependence of $n_{2}$ and $\beta$ for the $\sim 10$ nm graphene nanoribbons considered in Fig. 1.

Appendix H Convergence of the fast-Fourier transform method for nanoislands

In this work, the nonlinear response for finite graphene nanoislands is described using the perturbative expansion procedure for the single-electron density matrix outlined in Ref. Cox and García de Abajo (2014). In this approach, we compute several noninteracting RPA susceptibilities, $\chi^{0}_{ll^{\prime}}(s\omega)$, which have the same form as Eq. (23) of the Appendix (but lack any momentum dependence), to obtain the nonlinear polarizabilities for a single nanoisland. To expedite these computations, which require summing $\sim N^{4}$ terms, $N$ being the number of carbon atoms in a nanoisland, we employ a fast-Fourier transform (FFT) method, as described in Ref. Thongrattanasiri et al. (2012), for which only $\sim N^{3}$ operations are required. To perform the FFT, a finite grid of $N_{\omega}$ equally-spaced frequencies must be defined, and it is this discretization that determines the convergence of the $\chi^{0}_{ll^{\prime}}(s\omega)$ computed using the FFT method with that obtained from a direct evaluation. In Fig. 8, we show that while the FFT method can provide excellent convergence for $N_{\omega}=10^{12}$ frequencies, particularly in the linear response, more satisfactory results for the nonlinear polarizabilities are obtained with $10^{13}$ frequencies. Note that while this convergence is rather independent of a graphene nanoisland’s size, more frequencies are required as the phenomenological relaxation rate is reduced.

Acknowledgements.

This work has been supported in part by the Spanish MINECO (MAT2014-59096-P and SEV-2015-0522) and the European Commission (Graphene Flagship CNECT-ICT-604391 and FP7-ICT-2013-613024-GRASP). The authors thank Andrea Marini and Renwen Yu for useful discussions.