Theory of plasmonic effects in nonlinear optics: the case of graphene

# Theory of plasmonic effects in nonlinear optics: the case of graphene

Habib Rostami Istituto Italiano di Tecnologia, Graphene Labs, Via Morego 30, I-16163 Genova, Italy    Mikhail I. Katsnelson Radboud University, Institute for Molecules and Materials, NL-6525 AJ Nijmegen, The Netherlands    Marco Polini Istituto Italiano di Tecnologia, Graphene Labs, Via Morego 30, I-16163 Genova, Italy
###### Abstract

We develop a microscopic large- theory of electron-electron interaction corrections to multi-legged Feynman diagrams describing second- and third-order nonlinear response functions. Our theory, which reduces to the well-known random phase approximation in the linear-response limit, is completely general and is useful to understand all second- and third-order nonlinear effects, including harmonic generation, wave mixing, and photon drag. We apply our theoretical framework to the case of graphene, by carrying out microscopic calculations of the second- and third-order nonlinear response functions of an interacting two-dimensional (2D) gas of massless Dirac fermions. We compare our results with recent measurements, where all-optical launching of graphene plasmons has been achieved by virtue of the finiteness of the quasi-homogeneous second-order nonlinear response of this inversion-symmetric 2D material.

## I Introduction

Recent years have witnessed momentous interest koppens_nanolett_2011 (); grigorenko_naturephoton_2012 (); basov_rmp_2014 (); low_acsnano_2014 (); abajo_acsphoton_2014 () in the collective density oscillations of a doped graphene sheet, the so-called Dirac plasmons. The reason is partly related to the fact that the propagation of graphene plasmons has been directly imaged in real space by utilizing scattering-type scanning near-field optical microscopy fei_nature_2012 (); chen_nature_2012 (). In a series of pioneering experiments in the mid-infrared spectral range, Fei et al. fei_nature_2012 () and Chen et al. chen_nature_2012 () demonstrated that the plasmon wavelength can be - times smaller than the free-space excitation wavelength , allowing an extreme concentration of electromagnetic energy, and that Dirac plasmon properties are easily gate tunable.

Importantly, these figures of merit have been dramatically improved by utilizing van der Waals stacks geim_nature_2013 () comprising graphene encapsulated in boron nitride crystals mayorov_nanolett_2011 (); wang_science_2013 (); bandurin_science_2016 (). Mid-infrared plasmons in encapsulated graphene display woessner_naturemater_2015 () ultra-large field confinement (i.e. ), small group velocity, and a remarkably long lifetime, . In the Terahertz spectral range, acoustic plasmons with and a similar lifetime have been recently observed in hBN/graphene/hBN heterojunctions including a nearby metal gate alonsogonzalez_naturenano_2016 ().

Substantial theoretical efforts have also been devoted to understanding the nonlinear optical properties of graphene jafari_jpcm_2012 (); mikhailov_prb_2014 (); cheng_njp_2014 (); wehling_prb_2015 (); cheng_prb_2015 (); Habib_prb_2016 (); mikhailov_prb_2016 (); mikhailov_arxiv_2016 (); savostianova_arxiv_2016 (); cheng_arxiv_2016 (). Experimentally, Hendry et al. hendry_prl_2010 () demonstrated that the third-order optical susceptibility of graphene is remarkably large () and only weakly dependent on wavelength in the near-infrared frequency range. Third-harmonic generation (THG) from mechanically exfoliated graphene sheets has been measured by Kumar et al. kumar_prb_2013 () who extracted a value of the third-order susceptibility on the order of for an incident photon energy . Finally, Hong et al. hong_prx_2013 () reported strong THG in graphene grown by chemical vapor deposition, in the situation in which the incident photon energy is in three-photon resonance with the exciton-shifted van Hove singularity.

Since plasmons enable the concentration of electromagnetic energy into extremely small volumes, many groups have theoretically studied the interplay between plasmons and the nonlinear optical properties of graphene and its nanostructures Mikhailov_prb_2011 (); Mikhailov_prl_2014 (); Yao_prl_2014 (); Manzoni_njp_2015 (); Tokman_prb_2016 (); cox_acsnano_2016 (); Constant_arXiv_2016 (). An all-optical plasmon coupling scheme, which takes advantage of the intrinsic nonlinear optical response of graphene, has been implemented experimentally Constant_np_2016 (). Free-space, visible light pulses were used by the authors of Ref. Constant_np_2016, to launch Dirac plasmons in graphene. Difference-frequency wave mixing (see below) enabled the achievement of the energy- and wavevector-matching conditions. By carefully controlling the phase matching conditions, they also showed that one can excite Dirac plasmons with a definite wavevector and direction across a large frequency range, with an estimated efficiency approaching .

In this Article, we present a formal theory of second- and third-order nonlinearities, which treats quantum effects, intra- and inter-band contributions, and electron-electron interactions on equal footing. Our theory starts from an equilibrium Matsubara approach and related Feynman diagrams for the non-interacting nonlinear susceptibilities Habib_prb_2016 () and includes electron-electron interactions via a large- approach Coleman (). Here, refers to the number of fermion flavors. In this approximation, only diagrams with the largest number of fermion loops are kept Coleman (), with the idea that each fermion loop (bubble) carries a factor . Our large- theory reduces to the ordinary Bohm-Pines random phase approximation (RPA) Pines_and_Nozieres (); Giuliani_and_Vignale () in the case of linear response theory. This Article therefore naturally generalizes RPA theory to the case of nonlinear response functions, capturing screening and plasmons.

While our approach is completely general, we carry out detailed microscopic calculations for the case of a system of two-dimensional (2D) massless Dirac fermions Katsnelson (); kotov_rmp_2012 () interacting via long-range Coulomb interactions. Large- theories are known to work very well for weakly correlated materials, like graphene, in which long-range Coulomb interactions (rather than on-site Hubbard-type interactions) play a major role, while they fail to describe e.g. excitonic effects in semiconductors. In this case, vertex corrections need to be taken into account. This is well beyond the scope of the present Article and is left for future works.

Our Article is organized as following. In Sect. II we present our large- theory of the second-order susceptibility, while the case of the third-order response is analyzed in Sect. III. In Sect. IV we detail the derivation of long-wavelength expressions for the second- and third-order density response functions, which are extremely useful to understand nonlinear optics experiments. Symmetry considerations that apply to the case of homogeneous and isotropic 2D systems are summarized in Sect. V. In Sect. VI we present explicit analytical and numerical calculations of the second-order conductivity of a 2D system of massless Dirac fermions. For the sake of clarity, the case of harmonic generation and sum/difference wave mixing are independently analyzed in Sects. VII and VIII, respectively. A comparison between our theory and available experimental results Constant_np_2016 () is reported in Sect. VIII.1. A summary of our main results and a brief set of conclusions are reported in Sect. IX. Three appendices report a wealth of technical details. In particular, in Appendix A, we show that the main formal results of our Article, Eqs. (II), (III), and  (III), can also be independently derived by using the time-dependent Hartree approximation Giuliani_and_Vignale (). This derivation highlights a pathway to transcend the large- approximation, by suggesting a straightforward approach to include exchange and correlation effects in the spirit of density functional theory Giuliani_and_Vignale ().

## Ii Second-order density response in the large-N limit

We start by considering the bare second-order density response function. This is diagrammatically represented in Fig. 1(a). In this diagram (usually termed “triangle” diagram), solid lines are non-interacting Green’s functions while filled circles represent density vertices Habib_prb_2016 ().

In the large- approximation Coleman (), electron-electron interactions are captured by the diagrams reported in Figs. 1(b)-(d). In Fig. 1(b), all the three density vertices of the bare diagram are dressed by the infinite RPA series of bubble diagrams shown in Fig. 1(e). Similarly, in Fig. 1(c) and (d), only two vertices (one vertex) are (is) dressed by the RPA series of bubble diagrams.

The logic of keeping only the diagrams in Fig. 1(b)-(d) is the following. In the large- approximation, only diagrams that dominate in the limit are retained, where stands for the number of fermion flavors (, for example, in graphene). As discussed in Ref. Coleman, , each fermion loop (i.e. bubble) brings a factor , while each electron-electron interaction line brings a factor . Therefore, a diagram with bubbles and electron-electron interaction lines scales like , where . In the limit , all diagrams with are negligible, while diagrams with dominate. Diagrams in Fig. 1(b)-(d) have , while the diagram in Fig. 1(f) has . We therefore conclude that the latter diagram must be discarded in the large- approximation.

The sum of diagrams in Fig. 1(b)-(d) renormalizes the bare second-order response shown in Fig. 1(a). We find that the second-order nonlinear response function in the large- limit, which will be denoted by the symbol , is given by the following expression:

 χ(2)=χ(2)0R2 , (1)

where

 1R2 ≡ 1+∑i=1,2,Σviχ(1)(i) (2) + ∑i=1,2,Σv1χ(1)(1)v2χ(1)(2)vΣχ(1)(Σ)viχ(1)(i) + v1χ(1)(1)v2χ(1)(2)vΣχ(1)(Σ) .

In Eq. (2) we have used the following shorthand: and . Here, is the 2D Fourier transform of the Coulomb interaction potential and

 χ(1)(−q,q,−ω,ω)≡χ(1)0(q,ω)1−vqχ(1)0(q,ω) , (3)

is the usual RPA series of bubble diagrams Pines_and_Nozieres (); Giuliani_and_Vignale (), where is the frequency- and wavevector-dependent first-order non-interacting density response function Pines_and_Nozieres (); Giuliani_and_Vignale (). Finally, for , we have and .

Carrying out the sums in Eq. (2), we find

 R2=ϵ(Σ)ϵ(2)ϵ(1) , (4)

where is a shorthand for

 ϵ(qi,ωi)=1−vqiχ(1)0(qi,ωi) , (5)

which is the dynamical RPA screening function Pines_and_Nozieres (); Giuliani_and_Vignale ().

Therefore, the second-order density-density response function in the large- limit is given by

 χ(2)(−qΣ,q1,q2,−ωΣ,ω1,ω2)= χ(2)0(−qΣ,q1,q2,−ωΣ,ω1,ω2)ϵ(qΣ,ωΣ)ϵ(q2,ω2)ϵ(q1,ω1) . (6)

In the harmonic case, , Eq. (II) reduces to a result that has been obtained earlier by using a self-consistent density-matrix approach Mikhailov_prb_2011 (); Mikhailov_prl_2014 ().

## Iii Third-order density response in the large-N limit

In this Section we lay down a large- theory for the third-order response function. In this case, the situation is more subtle. The point is that the bare third-order response function (square diagram) contains four density vertices, see Fig. 2(a). One can therefore create two families of large- diagrams that contribute to the third-order response. The first family, which is based on square-type diagrams, is shown in Fig. 2. These large- series just renormalizes the bare third-order response, as in the case of the second-order response in Fig. 1. The second family is topologically distinct and based on triangle-type diagrams, as illustrated in Fig. 3. The idea of the second family is that you can create Feynman diagrams with four density vertices by “glueing” together two second-order triangular diagrams via an electron-electron interaction line. The sum of the diagrams in the first family will be denoted by the symbol , while the sum of the diagrams in the second family will be denoted by . The full third-order response function in the large- approximation is given by: .

In analogy with the second-order response function (II), the sum of the diagrams in Fig. 2 can be written as

 χ(3)a=χ(3)0R3 , (7)

where

 1R3≡1+∑iviχ(1)(i)+ 12∑i,jv1χ(1)(1)v2χ(1)(2)v3χ(1)(3)vΣχ(1)(Σ)viχ(1)(i)vjχ(1)(j)(1−δij)+ ∑iv1χ(1)(1)v2χ(1)(2)v3χ(1)(3)vΣχ(1)(Σ)viχ(1)(i)+ v1χ(1)(1)v2χ(1)(2)v3χ(1)(3)vΣχ(1)(Σ) . (8)

Here, and is the Kronecker delta. Carrying out the sums in Eq. (III), we find

 R3=ϵ(Σ)ϵ(3)ϵ(2)ϵ(1) (9)

and

 χ(3)a(−qΣ,q1,q2,q3,−ωΣ,ω1,ω2,ω3)= χ(3)0(−qΣ,q1,q2,q3,−ωΣ,ω1,ω2,ω3)ϵ(qΣ,ωΣ)Π3i=1ϵ(qi,ωi) . (10)

The situation is quite different for the second family of Feynman diagrams shown in Fig. 3. The sum of these diagrams can be written as

 χ(3)b=∑i=1,2,3χ(2)(i)viχ(2)0(i)Ki , (11)

where

 χ(2)(i)≡χ(2)(−qΣ,qi,˜qi,−ωΣ,ωi,˜ωi) , (12)
 χ(2)0(1)≡χ(2)0(−˜q1,q2,q3,−˜ω1,ω2,ω3) , (13)
 χ(2)0(2)≡χ(2)0(−˜q2,q3,q1,−˜ω2,ω3,ω1) , (14)
 χ(2)0(3)≡χ(2)0(−˜q3,q1,q2,−˜ω3,ω1,ω2) , (15)

and

 1Ki = 1+∑j=1,2,3vjχ(1)(j)(1−δij) (16) + v1χ(1)(1)v2χ(1)(2)v3χ(1)(3)viχ(1)(i) .

In Eqs. (12)-(15), and .

From Eq. (16), one can show that , with similar expressions holding for and , provided that suitable cyclic permutations of the ,, and indices are carried out.

After lengthy but straightforward algebra we conclude that

 χ(3)b(−qΣ,q1,q2,q3,−ωΣ,ω1,ω2,ω3)= 3∑i=1vsc(˜qi,˜ωi)χ(2)0(−qΣ,qi,˜qi,−ωΣ,ωi,˜ωi)ϵ(qΣ,ωΣ)Π3l=1ϵ(ql,ωl)× χ(2)0(−˜qi,qj,qk,−˜ωi,ωj,ωk) , (17)

where for and so on and so forth, in a cyclic manner, and the dynamically screened interaction is defined by Giuliani_and_Vignale ()

 vsc(q,ω)≡vqϵ(q,ω) . (18)

An alternative derivation of Eqs. (II), (III), and (III), which is based on the time-dependent Hartree approximation, is offered in Appendix A.

## Iv Long-wavelength expansion of nonlinear density response functions

In this Section we present a long-wavelength expansion of the nonlinear response functions introduced in the previous Sections. To this end, we take advantage of the gauge invariance principle and introduce nonlinear conductivity tensors.

Using gauge invariance habib_gaugeinv_2016 (), we obtain the following relation between the -th order nonlinear density response function and the corresponding nonlinear conductivity:

 χ(n)=(−i)nωΣ∑ℓ,{αi}qΣ,ℓΠni=1qi,αiσ(n)ℓα1…αn , (19)

where and are the Cartesian components of the vectors and , respectively. In writing Eq. (19) we have dropped for simplicity the argument of the nonlinear functions and : and .

Using Eq. (19), we can first express the dynamical screening function in terms of the linear-response conductivity tensor:

 ϵ(q,ω)=1+ivq∑ℓαqℓqαωσ(1)ℓα(−q,q,−ω,ω) . (20)

Using Eqs. (1), (7), (11), and (19), we obtain the following formal relations for the second- and third-order conductivities:

 σ(2),eeℓα1α2=σ(2)ℓα1α2R2 , (21)

and

 σ(3),eeℓα1α2α3=σ(3)ℓα1α2α3+˜σ(3)ℓα1α2α3R3 . (22)

In Eqs. (21) and (22), and denote the second- and third-order conductivities of the interacting electron system, while and denote their non-interacting counterparts.

Also, in Eq. (22) we have introduced

 ˜σ(3)ℓα1α2α3(−qΣ,q1,q2,q3,−ωΣ,ω1,ω2,ω3) = −i3∑i=1{vsc(˜qi,˜ωi)˜ωi∑β,β′˜qi,β˜qi,β′σ(2)ℓα1β(−qΣ,qi,˜qi,−ωΣ,ωi,˜ωi) (23) × σ(2)β′α2α3(−˜qi,qj,qk,−˜ωi,ωj,ωk)} .

The contribution denoted by the symbol stems from the family of diagrams shown in Fig. 3. As we have discussed earlier, a similar contribution does not exist in the case of the second-order response—cf. Eq. (21). Physically, represents an interaction-induced nonlocal contribution to the third-order conductivity. A caveat is now in order. Some care must be exercised when adding and in Eq. (22). In principle, indeed, one should expand in powers of wavevectors in the long-wavelength limit, up to the same order that appears in Eq. (23). This calculation is very cumbersome and will be not carried out in this work. The numerical results in Fig. 5(c) below have been calculated by neglecting nonlocal corrections to .

We are now in the position to expand the nonlinear density response functions in the long-wavelength limit. To this end, we just need to expand the conductivity tensors , keeping the leading contributions. We also need to specify the functional dependence of on . For long-range Coulomb interactions in a free-standing graphene sheet, the 2D Fourier transform of the Coulomb potential is given by , with and the vacuum permittivity. In the long-wavelength limit the RPA dynamical screening function can be expanded as

 ϵ(q,ω)=1+iq2ϵ0ωσ(1)L(ω)+… , (24)

where we have introduced the longitudinal linear conductivity

 σ(1)L(ω)≡∑ℓ,αqℓqαq2σ(1)ℓα(0,0,−ω,ω) . (25)

In Eq. (24) and below, “” denote higher-order corrections, which vanish faster that the leading term in the long-wavelength limit.

Similarly, the long-wavelength expansion of the second-order density response function requires an expansion of the second-order conductivity up to linear order in :

 σ(2)ℓα1α2(−qΣ,q1,q2,−ωΣ,ω1,ω2)=σ(2)ℓα1α2(−ωΣ;ω1,ω2) +∑i=1,2∑βqi,βd(2)ℓα1α2β,i(−ωΣ;ω1,ω2)+… , (26)

where the zeroth-order term, , is the second-order optical conductivity, while its dipole is defined by

 d(2)ℓα1α2β,i(−ωΣ;ω1,ω2)≡∂σ(2)ℓα1α2∂qi,β∣∣{q1,q2}→0 . (27)

Using Eqs. (II), (19) and (IV), we can rewrite the large- second-order density response in the long-wavelength limit as

 χ(2)(−qΣ,q1,q2,−ωΣ,ω1,ω2)=(−i)2q1q2qΣωΣ1R2 ×[σ(2)L(−ωΣ;ω1,ω2)+∑iqid(2)L,i(−ωΣ;ω1,ω2)]+…

where

 σ(2)L(−ωΣ;ω1,ω2)≡ ∑ℓ,α1,α2q1,α1q2,α2qΣ,ℓq1q2qΣ× (29) σ(2)ℓ,α1,α2(−ωΣ;ω1,ω2) ,

and

 d(2)L,i(−ωΣ;ω1,ω2)≡ ∑ℓ,α1,α2,βq1,α1q2,α2qΣ,ℓqi,βq1q2qΣqi× (30) d(2)ℓα1α2β,i(−ωΣ;ω1,ω2) .

In Eq. (IV), we have introduced the following long-wavelength expansion of the factors:

 Rn = [1+iqΣ2ϵ0ωΣσ(1)L(ωΣ)]Πnj=1[1+iqj2ϵ0ωjσ(1)L(ωj)] (31) + …

with .

Similarly, we can expand the third-order nonlinear density response functions. In the long-wavelength limit Eq. (7) reduces to

 χ(3)a(−qΣ,q1,q2,q3,−ωΣ,ω1,ω2,ω3)=(−i)3q1q2q3qΣωΣ ×1R3σ(3)L(−ωΣ;ω1,ω2,ω3)+… (32)

where the longitudinal third-order optical conductivity is given by

 σ(3)L(−ωΣ;ω1,ω2,ω3) = ∑ℓ,{αi}q1,α1q2,α2q3,α3qΣ,ℓq1q2q3qΣ× (33) σ(3)ℓα1α2α3(−ωΣ;ω1,ω2,ω3) .

In the same limit, Eq. (11) reduces to

 χ(3)b(−qΣ,q1,q2,q3,−ωΣ,ω1,ω2,ω3)=q1q2q3qΣωΣ(−i)42ϵ0R3 3∑i=1˜qi/˜ωi[1−i˜qi2ϵ0˜ωiσ(1)L(˜ωi)][σ(2)L(−ωΣ;ωi,˜ωi)+ qid(2)L,i(−ωΣ;ωi,˜ωi)+˜qid(2)L,2(−ωΣ;ωi,˜ωi)][σ(2)L(−˜ωi;ωj,ωk) +qjd(2)L,1(−˜ωi;ωj,ωk)+qkd(2)L,2(−˜ωi;ωj,ωk)]+…. (34)

Notice again that for , and so on and so forth, in a cyclic way.

## V Symmetry considerations for homogeneous and isotropic 2D systems

In a homogeneous and isotropic 2D system, the properties of the nonlinear conductivity tensors are highly constrained by rotational, translational, and inversion symmetries.

We start by recalling that, due to mirror () symmetry, all elements of the third-order conductivity with an odd number of () Cartesian indices are identically zero. Full rotational symmetry implies

 σ(3)xxxx=σ(3)xxyy+σ(3)xyxy+σ(3)xyyx . (35)

Moreover, mirror symmetry with respect to diagonal in the - plane provides an exchange symmetry between and Cartesian indices: we therefore have

 σ(3)yyyy=σ(3)xxxx,   σ(3)yyxx=σ(3)xxyy σ(3)yxyx=σ(3)xyxy,   σ(3)yxxy=σ(3)xyyx . (36)

By using Eqs. (33), (35), and (V) and for , we get

 σ(3)L = σ(3)xxyy+σ(3)xyxy2cos(θ1+θ2−θ3−θΣ) (37) + σ(3)xxyy+σ(3)xyyx2cos(θ1+θ3−θ2−θΣ) + σ(3)xyxy+σ(3)xyyx2cos(θ2+θ3−θ1−θΣ) .

Here, is the azimuthal angle of , i.e.

 cos(θΣ−θ1)=q1+q2cos(θ21)+q3cos(θ31)√∑3i=1q2i+2∑i>jqiqjcos(θij) (38)

with . All conductivity tensor elements in Eq. (37) have argument . Eq. (37) reduces to for the particular case of THG.

In an inversion symmetric system, we have . However, a non-vanishing dipole of the second-order conductivity is expected. Since is a rank- tensor, it obeys the same symmetry properties of the third-order nonlinear conductivity.

Because of the intrinsic permutation symmetry of the second-order conductivity tensor, we have . Moreover, as demonstrated in Appendix B, .

We can therefore write the following result for :

 d(2)L,1(−ωΣ;ω1,ω2)=d(2)xyyx,1cos(2θ1−θ2−θΣ) +[d(2)xyxy,1+d(2)xyyx,1]cos(θΣ−θ2) . (39)

All tensor elements on the right-hand side of Eq. (V) have argument . For the case of , we have (see Appendix B)

 d(2)L,2(−ωΣ;ω1,ω2)=d(2)xyyx,1cos(2θ2−θ1−θΣ) +(d(2)xyxy,1+d(2)xyyx,1)cos(θΣ−θ1) . (40)

From now on, we use the symbol as a shorthand for .

## Vi Second-order optical conductivity and its dipole in 2D Dirac materials

We are now ready to specialize our general results to the case of a specific material.

We will consider a 2D Dirac material with two valleys, like graphene. The low-energy Hamiltonian reads as following Katsnelson (); kotov_rmp_2012 (): , where is the Fermi velocity, stands for the valley (, ) index, and represent ordinary Pauli matrixes acting in sublattice space. The eigenstates (i.e. bands) of the this Hamiltonian are given , where indicates conduction and valence bands. The corresponding eigenvectors are , where is the polar angle of the vector . The wavefunctions in real space are where is the 2D electron system area. The matrix elements of the charge density () and charge current () operators are given by and where stands for the anti-commutation and the paramagnetic current operator in the first quantization picture reads Therefore, .

In the scalar potential gauge, the second-order charge current reads as follows Habib_prb_2016 ():

 J(2)ℓ(qΣ,ωΣ) = Π(2)ℓ(−qΣ,q1,q2,−ωΣ,ω1,ω2) (41) × V(q1,ω1)V(q2,ω2) .

Here, denotes the Fourier transform of the external scalar potential and is the second-order response function that establishes a link between the current response and the product of two external scalar potentials, and . The quantity is diagrammatically represented by a triangular diagram similar to the one in Fig. 1a), with two density vertices and one current vertex, i.e.

 Π(2)ℓ(−qΣ,q1,q2,−ωΣ,ω1,ω2)=−e3vFS∑k∑{λi}′∑P Fℓ,λ1λ2λ3(k,q1,q2)ℏωΣ+Eλ1,k−Eλ3,k+qΣ[nF(Eλ1,k)−nF(Eλ2,k+q1)ℏω1+Eλ1,k−Eλ2,k+q1− nF(Eλ2,k+q1)−nF(Eλ3,k+q</