Second-order nonlinear optical response of graphene

# Second-order nonlinear optical response of graphene

Yongrui Wang Department of Physics and Astronomy, Texas A&M University, College Station, TX, 77843 USA    Mikhail Tokman Institute of Applied Physics, Russian Academy of Sciences    Alexey Belyanin Department of Physics and Astronomy, Texas A&M University, College Station, TX, 77843 USA
September 24, 2019
###### Abstract

Although massless Dirac fermions in graphene constitute a centrosymmetric medium for in-plane excitations, their second-order nonlinear optical response is nonzero if the effects of spatial dispersion are taken into account. Here we present a rigorous quantum-mechanical theory of the second-order nonlinear response of graphene beyond the electric dipole approximation, which includes both intraband and interband transitions. The resulting nonlinear susceptibility tensor satisfies all symmetry and permutation properties, and can be applied to all three-wave mixing processes. We obtain useful analytic expressions in the limit of a degenerate electron distribution, which reveal quite strong second-order nonlinearity at long wavelengths, Fermi-edge resonances, and unusual polarization properties.

## I Introduction

Nonlinear optical properties of graphene have attracted considerable interest in the community. The magnitude of the matrix element of the interaction Hamiltonian describing coupling of massless Dirac electrons to light scales in proportion to , i.e. it grows more rapidly with wavelength than in conventional materials with parabolic energy dispersion, where it scales roughly as . This promises a strong nonlinear response at long wavelengths. Unfortunately, graphene is also a centrosymmetric medium for low-energy in-plane excitations, which suppresses second-order nonlinear response in the electric dipole approximation. Therefore, most of the effort was concentrated on the third-order nonlinear processes that are electric dipole-allowed. Recent theoretical proposals and some experiments include third-harmonic generation Kumar et al. (2013); Hong et al. (2013), four-wave mixing Hendry et al. (2010); Gu et al. (2012); Sun et al. (2010) and current-induced second-harmonic generation Glazov and Ganichev (2014); Bykov et al. (2012); Cheng et al. (2014). In few-layer graphene, second-harmonic generation (SHG) arising from the interactions between layers, which breaks the inversion symmetry, has been observed Dean and van Driel (2009, 2010).

The aim of this paper is to show that monolayer graphene does demonstrate quite significant second-order nonlinearity at long wavelengths despite its inversion symmetry. Here and throughout the paper, we will discuss only the 2D (surface) nonlinearity due to in-plane motion of electrons. Like any surface, graphene exhibits anisotropy between in-plane and out-of-plane electron motion. However, the corresponding second-order nonlinearity is very small and we will not discuss it here.

We develop the full quantum-mechanical theory of the in-plane second-order nonlinear response beyond the electric dipole approximation. In this case one has to consider oblique or in-plane propagation of electromagnetic waves. A non-zero in-plane second-order susceptibility of monolayer graphene appears when one includes the dependence of on the in-plane photon wave vectors, i.e. the spatial dispersion. Physically, this means that the inversion symmetry of graphene is broken by the wave vector direction. The spatial dispersion in momentum space is of course equivalent to the nonlocal response in real space. Spatial dispersion effects turn out to be quite large because of a large magnitude of the electron velocity . A non-zero value of the nonlocal has been pointed out before for second-harmonic generation Mikhailov (2011); Glazov (2011); Smirnova et al. (2014) (which only included intraband transitions in a free-carrier model), difference-frequency generation Yao et al. (2014), and parametric frequency down-conversion Tokman et al. (2016). The latter two papers developed a quantum theory including both intraband and interband transitions and applied it to the nonlinear generation of surface plasmons. In the recent experiment Constant et al. (2016), evidence for the difference-frequency generation of surface plasmons in graphene was reported. Here we provide a systematic derivation of the second-order nonlinear conductivity tensor, valid for all second-order processes, all frequencies and doping densities, as long as the massless Dirac fermion approximation for a single-particle Hamiltonian is applicable. For graphene, this means the range of frequencies from zero (more precisely, from inverse scattering time) to the near-infrared. Our approach can be applied to any system of massless chiral Dirac fermions, for example surface states in topological insulators such as BiSe. The resulting nonlinear susceptibility tensor satisfies all symmetry and permutation properties, and predicts unusual polarization properties of the nonlinear signal. We also summarize main properties of the linear current as a necessary step in deriving the nonlinear response functions, and present a detailed discussion of its gauge properties and regularization.

## Ii Basic equations

Consider a 2D quantum system which in the absence of external fields can be described by the Dirac Hamiltonian

 ^H0(^\mathboldp)=vF^\mathboldσ⋅^\mathboldp, (1)

where , , , where are Pauli matrices. The spinor eigenfunctions of the Hamiltonian (1) are

 \mathboldΨ\mathboldk,s(\mathboldr)≡⟨\mathboldr|\mathboldk,s⟩=ei\mathboldk⋅\mathboldr√2A(seiθ(\mathboldk)), (2)

and the eigenenergies are where for conduction and valence bands, respectively; , is the angle between the electron momentum and -axis, and is the normalization area. This description is valid for carriers in monolayer graphene up to the energies of order 1 eV. For higher energies, quadratic and trigonal warping corrections become non-negligible.

Consider the most general light-matter interaction Hamiltonian utilizing both vector and scalar potentials: and . Following a standard procedure Landau and Lifshitz (2013); Gantmakher and Levinson (1987), we replace in the unperturbed Hamiltonian and add the potential energy operator assuming a particle with the charge . This gives

 ^H=^H0+^Hoptint,^Hoptint=evFc^\mathboldσ⋅\mathboldA−eφ⋅^1, (3)

where is a unit 22 matrix. The Hamiltonian in Eq. (3) leads to the von Neumann equation for the density matrix:

 iℏ∂∂tρmn=(Em−En)ρmn+∑l[(^Hoptint)mlρln−ρml(^Hoptint)ln], (4)

where .

We will consider a monochromatic electromagnetic field in plane of graphene,

 \mathboldE=12[\mathboldx0Ex(ω)+\mathboldy0Ey(ω)]e−iωt+iqx+C.C. (5)

or its bichromatic combinations. The field component can be ignored because neither this field component itself nor the magnetic field it generates can affect the 2D carrier motion. Furthermore, the component of the vector potential which generates the z-component of the electric field does not enter the Hamiltonian (3) because . The field described by Eq. (5) corresponds to the electromagnetic potentials

 φ =12ϕ(ω)e−iωt+iqx+C.C., \mathboldA =12[\mathboldx0Ax(ω)+\mathboldy0Ay(ω)]e−iωt+iqx+C.C.. (6)

Note that the P-polarized radiation can be defined through both the scalar potential,

 φ=12iEx(ω)qe−iωt+iqx+C.C., (7)

and the vector potential:

 \mathboldAP=12\mathboldx0cEx(ω)iωe−iωt+iqx+C.C. (8)

At the same time, the S-polarized radiation, can be defined only through the vector potential:

 \mathboldAS=12\mathboldy0cEy(ω)iωe−iωt+iqx+C.C. (9)

It is convenient to represent the surface current density generated in response to a harmonic field as a sum over spatial harmonics: , where ; the set of in-plane photon wave vectors is specified by appropriate conditions on the boundary of a large area . It is also convenient to choose the area to be a multiple of the normalization area , so that

 (2A)−1/2∫A\mathboldΨ∗n(\mathboldr)\mathboldΨm(\mathboldr)d2r=(2S)−1/2∫S\mathboldΨ∗n(\mathboldr)\mathboldΨm(\mathboldr)d2r. (10)

After calculating the matrix elements of the current density operator and solving independently the master equations (4), one can calculate the average amplitude of a given current density harmonic, which could be used as a source in Maxwell’s equations or to determine the conductivity tensor:

 \mathboldj(q)=∑mn\mathboldj(q)nmρmn. (11)

In order to evaluate we determine the velocity operator and define the current density operator as :

 ^\mathboldj=−evF^\mathboldσ. (12)

Next, we take into account a standard expression for the current density operator in a second-quantized formalism Landau and Lifshitz (2013): , where and are second-quantized operators, and and are fermion creation and annihilation operators. Treating and as Heisenberg operators and using , , we arrive at , which gives

 2−1\mathboldj(q)nm=⟨n|e−i\mathboldq⋅\mathboldr^\mathboldj|m⟩. (13)

To calculate the matrix elements and we will need the following useful relationships:

 (eiqx)mn =12(smsn+ei(θn−θm))δ\mathboldkm,\mathboldkn+\mathboldq, (14) (^\mathboldσeiqx)mn =12[(\mathboldx0−i\mathboldy0)smeiθn+(\mathboldx0+i\mathboldy0)sne−iθm]δ\mathboldkm,\mathboldkn+\mathboldq. (15)

The above general equations should allow one to calculate the conductivity in any order with respect to the external optical field. There is however a complication related to the fact that the model described by the effective Hamiltonian Eq. (1) contains a ”bottomless” valence band with electrons occupying all states to . Therefore, only the converging integrals make sense:

 ∑mn\mathboldj(q)nmρmn⇒g∑ss′\lx@underaccentset∞′∫d2k′4π2\lx@underaccentset∞∫d2k4π2\mathboldj(q)\mathboldk′\mathboldks′sρ\mathboldk\mathboldk′ss′, (16)

where is the degeneracy factor. Otherwise the optical response could be determined by the electron dispersion far from the Dirac point where the effective Hamiltonian Eq. (1) is no longer valid. It turns out that the convergence of the linear current depends on the choice of the gauge, whereas for the second-order nonlinear current the integral in Eq. (16) converges for any gauge. The divergence of the linear response can be regularized as discussed in the next section. In addition, the gauge dependence of the linear response violates gauge invariance, which is a consequence of the fact that the density matrix corresponding to the bottomless Hamiltonian in Eq. (1) has an infinite trace. In the next section we discuss this issue in more detail.

## Iii The linear response of massless Dirac fermions

The perturbation expansion of the nonlinear response functions implies that the second-order nonlinear terms depend on the first-order linear response. Therefore, in this section we outline the derivation of the linear current. The nontrivial aspect of this derivation is an apparent violation of gauge invariance and divergence of the linear current. We address these issues in this section and related Appendix sections.

The solution of the density matrix equation (4) in the linear approximation with respect to the field is

 (17)

where we defined .
Here and is a complex-valued amplitude of the linear perturbation of the density matrix. For a monochromatic current we have

 \mathboldj(q)(ω)=∑mn\mathboldj(q)mnρ(1)nm(ω). (18)

The expression (18) is evaluated in Appendix A. The most straightforward derivation is for a P-polarized field defined through a scalar potential, Eq. (7), since in this case the integral (16) converges. If we keep only the terms of the lowest order in (i.e. the linear terms since ), the resulting 2D (surface) conductivity tensor is independent of . In the limit of strong degeneracy or low temperatures, the relevant terms are
(i) intraband conductivity, which has a Drude-like form:

 σ(intra)xx(ω)=ie2vFkFπℏ(ω+iγ), (19)

(ii) and the interband term:

 σ(inter)xx(ω)=ie24πℏln[2vFkF−(ω+iγ)2vFkF+(ω+iγ)]. (20)

Here is Fermi momentum, and we also added the relaxation terms by replacing in Eq. (17); in the limit one can obtain from Eq. (20) the well known result for the interband conductivity Nair et al. (2008): , where is the Heaviside step function.

If we define the optical field with a vector potential, the same calculation will lead to divergent integrals. In this case the finite, and at the same time gauge-invariant, expression for the linear current at frequency can be obtained by subtracting the same current evaluated at zero frequency Falkovsky and Varlamov (2007):

 \mathboldj(q)(ω)=∑mn\mathboldj(q)mn[ρ1,Anm(ω)−ρ1,Anm(ω→0)]. (21)

Here is Eq. (17) with in the interaction Hamiltonian. This prescription cancels the divergent term and leads to the Kubo formula for the linear response. In our case Eq. (21) is equal to the sum of Eqs. (19) and (20) for the diagonal conductivities , and gives . The procedure in Eq. (21) can be justified by considering the graphene Hamiltonian with a small quadratic term in the energy dispersion:

 E=sℏvFk+ϵℏ2k22, (22)

where is a small parameter. Adding this term provides a bottom to the valence band. As shown in Appendix B, the linear current for such a system approaches Eq. (21) when .

For a P-polarized field which can be represented through both scalar and vector potentials the renormalization procedure in Eq. (21) is equivalent to the gauge transformation of the density matrix from the A-gauge (8) to the -gauge (7). Indeed, let the function correspond to the solution of Eq. (17) for the field defined in the gauge given by Eq. (8), whereas the function correspond to the gauge of Eq. (7). Since we just found that the sum is finite, it makes sense to try the transformation . The gauge transformation from and to and corresponds to the unitary transformation of the density matrix (see Appendix C)

 ~ρnm=∑qp(e−iefℏc)nqρqp(e+iefℏc)pm, (23)

where the scalar function determines the gauge transformation of the potentials

 (24)

In particular, the transformation from the vector potential (8) to scalar potential (7) is

 ∇f=−\mathboldAP. (25)

Within the linear approximation with respect to we obtain from Eq. (23):

 ρ1,APnm⇒ρ1,APnm−ieℏcfnm(nm−nn). (26)

Next, we will use the general relationship (see e.g. Tokman (2009))

 fnm=−iℏEn−Em(∇f⋅^\mathboldv+^\mathboldv⋅∇f2)nm, (27)

from which we obtain from that

 fnm=−iℏvF(^\mathboldσ⋅∇f)nmEn−Em. (28)

As a result, from Eqs. (26), (28) and (25) one gets

 ρ1,APnm(ω)⇒ρ1,APnm(ω)+evFc[^σxAx(ω)]nm(nm−nn)En−Em. (29)

Taking into account Eq. (17), Eq. (29) can be represented as , which is identical to Eq. (21).

The structure of transformation (23) makes it clear why the density matrix with an infinite trace can give rise to the divergent current. Consider the density matrix in the form , where is a small perturbation. The sum can converge in a certain gauge even if the trace diverges. However, the transformation (23) to a different gauge projects the diagonal of the matrix with an infinite trace onto off-diagonal elements, which can lead to the divergence in Eq. (16). The inverse is also true: the divergence can be eliminated by the transformation (23) as we have just shown above.

It is also clear that the separation of the response into intraband and interband components depends generally on the choice of the gauge since the transformation (23) mixes different contributions. At the same time, a correctly defined current has to be gauge-invariant.

## Iv Second-order nonlinear response

Now we consider the second-order nonlinear response to the bichromatic field which we will represent through the vector potential in order to describe both P- and S-polarized fields with the same formalism. We will write the in-plane field components at frequencies directed along unit vectors as

 \mathboldA=12\mathboldη1A(ω1)ei(\mathboldq1⋅\mathboldr∥−ω1t)+12\mathboldη2A(ω2)ei(\mathboldq2⋅\mathboldr∥−ω2t)+c.c. (30)

We need to calculate the perturbation of the density matrix at the sum frequency . The term quadratic with respect to the field can be written as

 =ρ(2)mn(ω1+ω2)=(e2c)1ℏ(ω1+ω2)−(ϵm−ϵn) ×∑l≠m,n[((^\mathboldv⋅\mathboldη1)ei\mathboldq1⋅\mathboldr)mlA(ω1)ρ(1)ln(ω2)−ρ(1)ml(ω1)((^\mathboldv⋅\mathboldη2)ei\mathboldq2⋅\mathboldr)lnA(ω2)]+{1↔2} =12(ec)2A(ω1)A(ω2)ℏ(ω1+ω2)−(ϵm−ϵn)×∑l≠m,n((^\mathboldv⋅\mathboldη1)ei\mathboldq1⋅\mathboldr)ml((^\mathboldv⋅\mathboldη2)ei\mathboldq2⋅\mathboldr)ln ×[(ρnn−ρll)ℏω2−(ϵl−ϵn)−(ρll−ρmm)ℏω1−(ϵm−ϵl)]+{1↔2}. (31)

The trace of the corresponding Fourier harmonic of the induced current can be then calculated as

 \mathboldj(q1+q2)(ω1+ω2)=−e∑mn(^\mathboldve−i(\mathboldq1+\mathboldq2)⋅\mathboldr)nmρ(2)mn(ω1+ω2). (32)

The second-order response at the difference frequency, can be obtained by replacing

 ω2⇒−ω2,q2⇒−q2,A(ω2)⇒A∗(ω2). (33)

Next, we transform from summation to integration over -states, introduce the corresponding occupation numbers of the momentum states in each band, apply the momentum conservation in a three-wave mixing process, and take into account spin and valley degeneracy. Note that the integral over the electron momenta converges, as opposed to the linear response calculations where one needs to regularize the integral by either subtracting the contribution at zero frequency or adding a term to the Hamiltonian, as discussed above. The result is

 =\mathboldJ(2)(ω1+ω2) =−e3v3F16π2c2ℏ2A(ω1)A(ω2)∑sm,sn,sl∫d2\mathboldk1(ω1+ω2)−vF(sm|\mathboldk+\mathboldq1|−sn|\mathboldk−\mathboldq2|) ×[f(sn,|\mathboldk−\mathboldq2|)−f(sl,|\mathboldk|)ω2−vF(sl|\mathboldk|−sn|\mathboldk−\mathboldq2|)−f(sl,|\mathboldk|)−f(sm,|\mathboldk+\mathboldq1|)ω1−vF(sm|\mathboldk+\mathboldq1|−sl|\mathboldk|)] ×[(η1x−iη1y)smeiθ(\mathboldk)+(η1x+iη1y)sle−iθ(\mathboldk+\mathboldq1)] ×[(η2x−iη2y)sleiθ(\mathboldk−\mathboldq2)+(η2x+iη2y)sne−iθ(\mathboldk)] ×[(\mathboldx0+i\mathboldy0)sme−iθ(\mathboldk−\mathboldq2)+(\mathboldx0−i\mathboldy0)sneiθ(\mathboldk+\mathboldq1)] +{1↔2}. (34)

This equation can be integrated numerically for any given geometry of incident fields and electron distribution. We consider the limit of the Fermi distribution with a strong degeneracy, direct all in-plane photon wave vectors along x-axis, and expand the integrand in Eq. (34) in powers of . The integral over the term of zeroth-order in vanishes, as expected from symmetry. We will keep the terms linear in . Also we have to evaluate separately the intraband contribution and all types of mixed interband-intraband contributions: , , and . After performing this procedure, we find the following nonzero components of the second-order nonlinear conductivity tensor, while all other components are zero:

 =σ(2)xxx(ω1+ω2;ω1,ω2) =−s(ϵF)e3v2F2πℏ21ω21ω22(ω1+ω2)1(ω21−4v2Fk2F)(ω22−4v2Fk2F)((ω1+ω2)2−4v2Fk2F) ×[−4v4Fk4F(q1ω32(2ω1+ω2)+q2ω31(ω1+2ω2))+16v6Fk6F(q1ω2(2ω1+ω2)+q2ω1(ω1+2ω2))], (35) =σ(2)xyy(ω1+ω2;ω1,ω2) =−s(ϵF)e3v2F2πℏ21ω21ω22(ω1+ω2)1(ω21−4v2Fk2F)(ω22−4v2Fk2F)((ω1+ω2)2−4v2Fk2F) ×[4(vFkF)2ω1ω2(ω1+ω2)2(q1ω22+q2ω21) =+4(vFkF)4(q1ω42−(6q1+4q2)ω1ω32−8(q1+q2)ω21ω22−(4q1+6q2)ω31ω2+q2ω41) =+16(vFkF)6(q1ω2(2ω1−ω2)+q2ω1(2ω2−ω1))], (36) =σ(2)yxy(ω1+ω2;ω1,ω2) =−s(ϵF)e3v2F2πℏ21ω21ω22(ω1+ω2)1(ω21−4v2Fk2F)(ω22−4v2Fk2F)((ω1+ω2)2−4v2Fk2F) ×[4(vFkF)2ω21ω2(ω1+ω2)(q1ω22−q2ω1(ω1+2ω2)) =+4(vFkF)4(q2ω1(ω1+2ω2)3−q1ω2(4ω31+4ω21ω2+2ω1ω22+3ω32)) =+16(vFkF)6(q1ω2(2ω1+3ω2)−q2ω1(ω1+2ω2))], (37) =σ(2)yyx(ω1+ω2;ω1,ω2) =−s(ϵF)e3v2F2πℏ21ω21ω22(ω1+ω2)1(ω21−4v2Fk2F)(ω22−4v2Fk2F)((ω1+ω2)2−4v2Fk2F) ×[4(vFkF)2ω1ω22(ω1+ω2)(q2ω21−q1ω2(2ω1+ω2)) =+4(vFkF)4(q1ω2(2ω1+ω2)3−q2ω1(3ω31+2ω21ω2+4ω1ω22+4ω32)) =+16(vFkF)6(q2ω1(3ω1+2ω2)−q1ω2(2ω1+ω2))]. (38)

Here depending on whether the Fermi level is in the conduction or valence band. A sketch of the second-order nonlinear process for an obliquely incident light of mixed polarization is shown in Fig. 1. Note that when both pump fields have either S- or P-polarization, the generated nonlinear current has only the -component (along the in-plane direction of propagation of the pumps). When the polarizations are mixed, the -component of the nonlinear current appears due to and components of the nonlinear conductivity (they are different only by permutation of indices 1 and 2 referring to the two pump fields).

Apparent “non-reciprocity” of the expressions for (P-in, S-out channel) and (S-in, P-out channel) has a simple physical explanation: a P-polarized incident field cannot create a current orthogonal to the electric field, whereas an incident S-polarized field creates such a current via the magnetic field component normal to the layer.

The expressions for the nonlinear conductivity tensor that we obtained pass all symmetry and gauge invariance tests. Indeed, one can verify that the value of agrees with the one derived using scalar potential in the interaction Hamiltonian. Furthermore, after converting the nonlinear conductivity to the nonlinear susceptibility according to

 χ(2)ijk(ω1+ω2;ω1,ω2)=iσ(2)ijk(ω1+ω2;ω1,ω2)ω1+ω2,

one can verify that all components of the nonlinear susceptibility tensor satisfy proper permutation relations; see e.g. Ch. 2.9 in Il’inskii and Keldysh (1994):

 χ(2)ijk(ω3=ω1+ω2)=χ(2)jik(−ω1=−ω3+ω2)=χ(2)kji(−ω2=−ω3+ω1), (39)

where in-plane wave vectors have to be permuted together with frequencies.

The second-order response goes to zero when the Fermi energy goes to zero, and has maxima at resonances when one of the three frequencies involved in three-wave mixing is close to . Far from these resonances and for high frequencies or low doping, , expressions for the nonlinear conductivity are greatly simplified (we will give only the expressions for and for brevity):

 σ(2)xxx=s(ϵF)2e3v2Fπℏ2v4Fk4F[q1ω32(2ω1+ω2)+q2ω31(2ω2+ω1)]ω41ω42(ω1+ω2)3, (40) σ(2)xyy=−s(ϵF)2e3v2Fπℏ2v2Fk2F(q1ω22+q2ω21)ω31ω32(ω1+ω2). (41)

An interesting and surprising result contained in these expressions is that the nonlinear frequency conversion of S-polarized radiation into P-polarized radiation is much more efficient at high frequencies as compared to the P-in, P-out channel: . The dominance of the S-in, P-out channel is due to the chiral nature of Dirac fermion states. In particular, for the second-harmonic generation process and , and the dominant component of the nonlinear conductivity tensor is simply

 σ(2)xyy=−s2e3πℏ2v4Fk2Fqω5. (42)

In the opposite limit of low frequencies or high doping, , we also obtain simplified expressions:

 σ(2)xxx=se3v2F4πℏ2ω1ω2(q1+q2ω1+ω2+q1ω1+q2ω2), (43) σ(2)xyy=−se3v2F4πℏ2ω1ω2[q1+q2ω1+ω2+ω1−ω2ω1+ω2(q1ω1+q2ω2)]. (44)

We verified that Eqs. (43) and (44) can be derived independently from the single-band kinetic equation, i.e. in the quasiclassical approximation described in Appendix D. This provides another test of our general expressions, since one should indeed expect that the single-band physics emerges in the limit of a strong doping and low frequencies, when all interband transitions become suppressed by Pauli blocking. Note that although Eqs. (43) and (44) do not depend on