Scaling laws for nonlinear electromagnetic responses of Dirac fermions
SScaling laws for nonlinear electromagnetic responses of Dirac fermions
Takahiro Morimoto
1, 2 and Naoto Nagaosa
2, 3 Department of Physics, University of California, Berkeley, CA 94720 RIKEN Center for Emergent Matter Science (CEMS), Wako, Saitama, 351-0198, Japan Department of Applied Physics, The University of Tokyo, Tokyo, 113-8656, Japan (Dated: October 9, 2018)We theoretically propose that the Dirac fermion in two-dimensions shows the giant nonlinearresponses to electromagnetic fields in terahertz region. A scaling form is obtained for the currentand magnetization as functions of the normalized electromagnetic fields
E/E ω and B/B ω , where thecharacteristic electric (magnetic) field E ω ( B ω ) depends on the frequency ω as (cid:126) ω /ev F ( (cid:126) ω /ev F ),and is typically of the order of 80 V/cm ( 8 mT) in the terahertz region. Applications of the presenttheory to graphene and surface state of a topological insulator are discussed. I. INTRODUCTION
Nonlinear electromagnetic responses (NLEMRs) [1–9]are essential for creating efficient electric/magnetic de-vices because their functions are closely related to thenonlinearity. In particular, the NLEMRs in the low en-ergy region, e.g., terahertz frequency region, is the focusof keen attention due to the development of the inten-sive terahertz excitations, and the strong coupling be-tween the free carriers and electromagnetic field (EMF)in semiconductors [10–13]. To achieve the strong non-linearity in low energy, two conditions must be satisfied.One is that there are enough excited states in that region,and the other is that the strong EMF can exist in thematerial. Nonlinear responses are enhanced when thereexist resonant states corresponding to the frequency ofthe external fields and the electric field within the mate-rial can be large by avoiding the screening effect. If wefocus on the low frequency region which is our interestin this paper, metals could be a good candidate to showthe NLEMRs for the first point, because there are largedensity of states in the low energy region. However, theelectromagnetic field is screened by the conducting elec-trons and hence cannot be strong enough to show thenonlinearity. In insulators, on the other hand, the energygap E g leads to the stability of the system, while the en-ergy denominator (cid:126) ω − E g appearing in the perturbationtheory reduces the higher order effects, and hence strongnonlinear effects. The resonance at ω ∼ E g is one wayto enhance the NLEMRs, but it does not work for thelow energy responses in the THz regime. Another sourceof the nonlinearity is the symmetry breaking such as theferroelectricity, which produces the giant responses bythe small external field. For example, photovoltaic cur-rents have recently been studied for ferroelectric materi-als [8, 9]. However, the switching speed of ferroelectricityis usually very slow, which deteriorates its applicabilityfor realistic devices.From this viewpoint, the Dirac fermion having lin-ear dispersions (Fig. 1) is located between the met-als and insulators, and has an advantage of havinglow energy excitations and poor screening. Therefore,it deserves to explore theoretically the possible giant FIG. 1. Schematic picture of nonlinear responses of a Diracfermion.
NLEMRs of two-dimensional Dirac fermions which arerealized in condensed matters at the surface of three-dimensional topological insulators (TIs) [14, 15], and alsoin graphene [16, 17]. The screening of external fields areweak for Dirac fermions since the density of states van-ish at zero energy. At the same time, the linear disper-sion of Dirac fermions enables that a light of any fre-quency matches the resonance condition, especially inthe terahertz region. Actually, it has been proposed thatgraphene is a promising candidate for the NLEMRs [18–24]. In addition, Dirac fermions in TIs show a non-trivial topological nature such as magneto-electric re-sponses [25, 26], which makes TIs the attractive systemfor the NLEMRs. Since Dirac fermions in TIs are sub-ject to the spin-orbit coupling, NLEMRs in TIs indicatenonlinear responses of spins. In this regard, nonlinearresponses of TIs will lead to novel spintronics functionsof TIs.In this paper, we study nonlinear responses of a two-dimensional Dirac fermion and derive a scaling form forthem. This provides a unified description of various non-linear phenomena supported by two-dimensional Diracfermions which are realized in graphene and at the sur-face of TIs. a r X i v : . [ c ond - m a t . m e s - h a ll ] M a r FIG. 2. Two regimes of nonlinear responses. Strong nonlin-earity appears when
E > E ω and E > E c . Weak nonlinearityappears otherwise. Similar phase diagram applies for ω − B plane. II. TWO REGIMES OF NONLINEARRESPONSES
We consider the Dirac fermion at the surface of a topo-logical insulator or in a graphene (neglecting the valleyand spin degrees of freedom) which is described by theDirac Hamiltonian, H = (cid:126) v F ( k x σ y − k y σ x ) , (1)where h = 2 π (cid:126) is the Planck constant, σ are Pauli matri-ces acting on spin degrees of freedom, and v F is the Fermivelocity of the Dirac fermion. The Green’s function forthe Dirac fermion is given by G ( iω ) = ( i (cid:126) ω − H ) − = (cid:126) − iω + v F ( k x σ y − k y σ x )( iω ) − v F k , (2)and the current operator is given by j x = ev F σ y . (3)We study NLEMRs of the Dirac fermion focusing onresponses of uniform current J and magnetization M ofthe Dirac fermion to external ac EMFs E and B . TheNLEMRs of J and M obey the following scaling forms J = eω v F f (cid:18) EE ω , BB ω (cid:19) , (4) M = eωg (cid:18) EE ω , BB ω (cid:19) , (5)because the Dirac fermion is gapless and the only energyscale is introduced by the frequency ω of the externalEMFs. Here, E ω and B ω are the typical scales of EMFsdetermined by the frequency of EMFs, and f and g are scaling functions. We note that the spatially uniformlimit q → E ( B ) is largerthan E ω and E c ( B ω and B c ), and (ii) weak nonlinearityotherwise, as summarized in Fig. 2. The typical order ofelectric fields E ω is determined by eE ω v F t = ht , (6)which equates the energy that a Dirac fermion driven by E gains by the time t and the energy scale correspondingto t . The typical order of magnetic fields B ω is deter-mined by ( v F t ) B ω = h e , (7)which equates the magnetic flux in an area ( v F t ) toa flux quantum. These two equations indicate thatthe nonlinear effects are qualitatively different when theEMFs are larger or smaller than E ω = hev F t = hω ev F , B ω = h ev F t = hω ev F , (8)where we set the typical time scale to be t = 1 /ω : (i)the NLEMRs are non-perturbative with respect to theapplied EMFs when E and B are stronger than thesetypical scales, (ii) the NLEMRs can be obtained by aperturbation theory with respect to E and B in the op-posite limit. If we consider the static limit ( ω → t = 1 /ω is replaced by the relaxationtime of electrons τ which is in the order of THz regime,i.e., τ (cid:39) E c and B c corresponding to the relaxation time τ are estimated to be E c (cid:39)
80 V/cm , B c (cid:39) . (9)Figure 2 shows the two typical regions of nonlinear re-sponses in the parameter space spanned by ω and E .NLEMRs show strong nonlinearity in a region boundedby E > E ω and E > E c . NLEMRs show weak non-linearity in the perturbative region E < E ω or in therelaxation-dominant region E < E c . III. STRONG NONLINEARITY
An insight into the strong nonlinearity regime is ob-tained by considering the static limit ω = 0 and E > E c .In this case, the effective action for static EMFs is knownfrom (2 + 1)D QED. The action is given by [28] S = (cid:90) d xL, L = A (cid:0) B − E (cid:1) , (10)where A = ζ (cid:0) (cid:1) e / (2 π ) is a numerical constant with ζ ( x ) being the zeta function and we adopted the conven-tion (cid:126) = v F = 1 [Note that E in this unit corresponds to( c/v F ) E in the natural unit ( c = 1)]. This effective actionwas obtained for the gauge fields with constant E and B fields coupled to a Dirac fermion. This (2 + 1)D QEDaction in Eq. (10) is understood from the energy den-sity of Landau levels for Dirac fermions as explained inAppendix A. By using this effective action, the nonlinearresponses and current and magnetization are summarizedas J = iω ∂L∂E = − i A ωE (cid:0) B − E (cid:1) − (11) M = − ∂L∂B = 3 A B (cid:0) B − E (cid:1) − , (12)which can be translated to the asymptotic behavior ofthe scaling functions f and g as f ( x, y ) ∼ x (cid:0) y − x (cid:1) − , g ( x, y ) ∼ y (cid:0) y − x (cid:1) − in the limit of large | x | and | y | .These can be interpreted by the more directly measurablequantities, i.e., the dielectric constant (cid:15) and magneticsusceptibility χ as (cid:15) = − d LdE = 3 A (cid:0) B − E (cid:1) − (cid:0) B − E (cid:1) , (13) χ = − d LdB = − A (cid:0) B − E (cid:1) − (cid:0) B − E (cid:1) . (14)If we focus on the magnetic susceptibility χ , we can read-ily see that χ shows a singularity for B → E = 0as χ = − A B − , (15)which indicates that the Taylor expansion with respect to B fields around B = 0 is not applicable and this is a non-perturbative result. The above formula for χ = dM/dB satisfies the scaling form for M in Eq. (5) with M (cid:39) ω ( B/B ω ) / . The magnetic susceptibility χ is subject toa nonlinear suppression with respect to small E as χ = − A B − (cid:18) − E B (cid:19) . (16)When the magnitude of E exceeds that of B ( | E | > | B | ),the action in Eq. (10) is no longer real and indicates thatthe system becomes unstable. This corresponds to theSchwinger process where applied electric fields E induceZener tunneling of an electron between adjacent Landaulevels in a successive manner and the energy distributionof the Dirac electrons is no longer stable. However, thescaling forms in Eqs. (5), (12) together with Eq. (10) stillwork. A theoretical study [29] on the Schwinger processof the Dirac fermion concluded the crossover from thelinear response regime J ∝ E to the nonlinear regime J ∝ τ E / at E ∼ = E c = E ω =1 /τ with τ being the relaxationtime. Both are consistent with Eqs. (5), (12) as f ( x, ∼ x ( | x | (cid:28) f ( x, ∼ x / ( | x | (cid:29)
1) when one replace ω by 1 /τ . In the case of the surface Dirac fermion in TIs,the current J is tied to the spin polarization (cid:104) s (cid:105) due tothe spin-momentum locking in TIs. Thus the Schwingerprocess of the Dirac fermions in TIs leads to a strongnonlinear spin generation with applied electric fields; weexpect the similar crossover for spin generation from thelinear response regime (cid:104) s (cid:105) (cid:39) E to the nonlinear responseregime (cid:104) s (cid:105) (cid:39) τ E / with strong E ( (cid:29) E c ). IV. WEAK NONLINEARITY
When the magnitudes of the applied electromagneticfields do not exceed the typical scales (
E < E ω , B < B ω )or they are below the the scales determined by the relax-ation time ( E < E c , B < B c ), the nonlinear responsesshow weak nonlinearity where perturbative treatmentswith respect to E and B are valid.Weak nonlinear responses are described by the expan-sion in terms of the applied E and B fields as J = β n,m ( ω ) B n ( ω ) E m ( ω ) , (17)where we consider Fourier components of J , B , and E oscillating in frequencies of the order of ω . The behaviorof β n,m ( ω ) can be determined by the scaling forms Eq.(5)as β n,m ∝ e n + m +1 (cid:126) n + m v n + m − F ω n +2 m − . (18)In a similar manner, the nonlinear response of the orbitalmagnetization to the applied E and B fields is given by M = χ n,m ( ω ) B n E m , (19)and χ n,m ∝ e n + m +1 (cid:126) n + m v n + mF ω n +2 m − . (20) V. NONLINEAR RESPONSE FUNCTIONS
Now we study nonlinear response functions in the weaknonlinearity regime. We obtain the response functions bythe diagram approach and confirm that they are consis-tent with the scaling form in Eqs. (18) and (20). (Thedetails of the calculations are given in Appendices.) Westart with a nonlinear suppression of the orbital magneticsusceptibility with electric fields. The orbital magneticsusceptibility χ ( ω ) given by M = χ , ( ω ) B, (21)is obtained from the current-current correlation function χ , ( ω ) = − ∂ ∂q y (cid:104) j x j x (cid:105) ( ω, q ) (cid:12)(cid:12)(cid:12)(cid:12) q = , (22)where q y is the wavenumber in the y -direction. Its realpart is written asRe[ χ , ( ω )] = − πe v F h δ ( ω ) . (23)(For details, see Appendix B.) The magnetic susceptibil-ity of Dirac fermions is diverging in the dc limit wherecutoffs are introduced by the relaxation time or the fi-nite temperature [30, 31]. Now we consider modulationof the orbital magnetic susceptibility by applying electricfields E ( ω ) which is captured by the nonlinear responsefunction M ( ω (cid:48) ) = χ , ( ω (cid:48) , ω ) B ( ω (cid:48) ) E ( ω ) E ( − ω ) . (24)Here, the magnetic susceptibility is measured in the fre-quency ω (cid:48) in the presence of laser light of the frequency ω , and χ , ( ω (cid:48) , ω ) quantifies the change of the magneticsusceptibility proportional to the laser intensity. In thelimit of static magnetic fields ( ω (cid:48) →
0) for a fixed fre-quency ω of the incident light, the nonlinear responsefunction is given byRe[ χ , ( ω (cid:48) , ω )] = πe v F h ω δ ( ω (cid:48) ) , (25)where the power of ω − is the one expected from thescaling analysis (for derivation, see Appendix B). Thusthe orbital magnetic susceptibility undergoes a nonlinearsuppression by applying the electric fields as M ( ω (cid:48) ) = − πe v F h δ ( ω (cid:48) ) (cid:18) − e v F h ω | E ( ω ) | (cid:19) B ( ω (cid:48) ) . (26)This is a nonlinear cross correlation effect between themagnetic properties and the electric fields. This nonlin-ear cross correlation effect is enhanced in the low fre-quency limit, typically in the terahertz regime.Next we discuss high harmonic generations which areschematically illustrated in Fig. 1. The N -th harmonicgeneration is a nonlinear process where the electron ab-sorbs N photons of frequency ω and emits one photonof frequency N ω . Here we focus on the third harmonicgeneration which is described by J (3 ω ) = β , ( ω ) E ( ω ) , (27)where J (3 ω ) is the induced current with frequency 3 ω ,and E ( ω ) is the electric field of incident light. [32]The response function β , ( ω ) is related to a correlationfunction of three current operators (as detailed in Ap-pendix C) and is given by β , ( ω ) = e (cid:126) v F ω . (28)This response function is usually denoted as χ (3) ( ω ) andtakes a value in the order of 10 − m / V for molecular liquids like CS in the terahertz regime [11]. In the caseof Dirac fermions in TIs, β , ( ω ) in Eq. (28) amounts tothe order of 10 − m / V for ω = 1 THz, if thin filmsof TIs are stacked in a separation of 10 nm. This isa significant enhancement compared to the conventionalmolecular liquids which is attributed to the gaplessness ofDirac fermions at the surface of TIs, as has been discussedin the case of graphene [19, 20]. We note that a similarorder of χ (3) was theoretically predicted to stacked layersof two-dimensional electron gas [33].Finally we discuss the nonlinear modulation of refrac-tive index by an incident light. This effect is describedby the response function J ( ω (cid:48) ) = ˜ β , ( ω (cid:48) , ω ) E ( ω (cid:48) ) E ( ω ) E ( − ω ) . (29)Here, ω (cid:48) is the frequency of transmitting light for whichwe are interested in the refractive index, and ω is thefrequency of incident light that modulates the refractiveindex. The response function ˜ β , ( ω (cid:48) , ω ) is a monotoni-cally decreasing function of ω (cid:48) for a fixed ω , and it takesthe asymptotic form ˜ β , ( ω (cid:48) → , ω ) = ( e / (cid:126) )( v F / ω )for the nonlinear modulation of the dc refractive indexwith a laser light of the frequency ω . (For details, see Ap-pendix D.) Again this power law behavior is consistentwith the scaling form Eq. (18), and hence, its significantenhancement is expected in the terahertz region. Fur-thermore, due to the spin-momentum locking in Diracfermions at the surface of TIs, the induced current in-dicates spin generation. Therefore, in the case of TIs,the nonlinear modulation of the refractive index leadsto a nonlinear Edelstein effect [34], i.e., nonlinear spingeneration induced by strong electric fields, as follows.When the weak electric field is applied to the surface ofa TI, there appears spin generation (cid:104) s (cid:105) proportional tothe applied field E . As E is increased, the above non-linear refractive index in Eq. (29) indicates the nonlinearspin generation (cid:104) s (cid:105) (cid:39) E . If E is further increased, thecrossover to the strong nonlinearity regime takes placeand the spin generation is governed by the Schwingerprocess with (cid:104) s (cid:105) (cid:39) tE / . VI. SUMMARY
To summarize, we have studied the nonlinear electro-magnetic responses of Dirac fermion from the viewpointof the scaling laws, which provide a unified view of var-ious nonlinear phenomena in this system. The nonlin-earity increases as the frequency ω decreases, and thestrong limit of the nonlinearity can be understood by aLagrangian for constant electric ( E ) and magnetic ( B )fields. The crossover strength E ω ( B ω ) of E ( B ) aredetermined by E ω = hω ev F ( B ω = hω ev F ), and is of theorder of E ω (cid:39)
80 V/cm ( B ω (cid:39) E = 120 V/cm [24] where non-linearity of thermalized excited electrons has been ob-served. Unfortunately, strong high harmonic generationhas not been observed in these experiments, which maybe attributed to some extrinsic effects such as disordereffects. Thus observation of scaling behavior of nonlinearresponses of Dirac fermions is expected in future experi-ments. Furthermore, given the recent advances of exper-iments on TIs [35, 36], Dirac fermions in TIs will offer aplatform for novel nonlinear electromagnetic responses. ACKNOWLEDGMENTS
We thank N. Ogawa and Y. Tokura for fruitful discus-sions. This work was supported by the EPiQS initiativeof the Gordon and Betty Moore Foundation and by JSPSGrant-in-Aid for Scientific Research (No. 24224009, andNo. 26103006) from MEXT, Japan.
Appendix A: (2+1)D QED action
The (2+1)D QED action for constant gauge strength can be obtained by considering Landau levels for Weylfermions. The Landau level energy and Landau level degeneracy are given by E N = sgn( N ) v F √ N (cid:126) eB, (A1)12 π(cid:96) = eB π (cid:126) . (A2)Then the energy density is given by H = (cid:88) N = −∞ E N π(cid:96) = ( eB ) v F π (cid:114) (cid:126) ζ (cid:18) − (cid:19) = v F ζ ( )2 π √ (cid:126) (cid:18) eB (cid:19) , (A3)where we used ζ ( − ) = ζ ( ) / π in the last line. If we replace B → ( B − E ) for the Lorenz invariance and set e = (cid:126) = v F = 1, we obtain the effective Lagrangian for (2 + 1)D QED. Appendix B: Orbital magnetic susceptibility
We study the orbital magnetic susceptibility, M = χ ( ω ) B. (B1)The magnetic susceptibility is obtained from the current-current correlation function as χ = − ∂ ∂q y (cid:104) j x j x (cid:105) , (B2) (cid:104) j x j x (cid:105) ( i Ω , q ) = 1(2 π ) (cid:90) d k (cid:88) iω n Tr[ j x G ( iω n + i Ω , k + q ) j x G ( iω n , k )] , (B3)where q y is the wavenumber in the y -direction. This is understood from M = ∇ × J and B = ∇ × A . The Green’sfunction and the current operator for the Weyl fermion are given by G ( iω, k ) = ( iω − H ) − , (B4) j i = ∂H∂k i , (B5)with H = k x σ y − k y σ x , (B6)where we used the convention (cid:126) = v F = 1. By using the relation ∂ q y G = − G∂ q y G − G = Gj y G , this is rewritten as χ ( i Ω) = − π ) (cid:90) d k (cid:88) iω n Tr[ j x G ( iω n + i Ω , k ) j y G ( iω n + i Ω , k ) j y G ( iω n + i Ω , k ) j x G ( iω n , k )] . (B7)The integration over the angle of k and the summation of Matsubara frequencies give χ ( ω ) = − π (cid:90) kdk k − k ω − ω k (4 k − ω ) , (B8)after the analytic continuation i Ω → ω . The imaginary part is obtained from the residue at k = ω asIm[ χ ( ω )] = − ω . (B9)By using the Kramers-Kronig relation, the real part is obtained asRe[ χ ( ω )] = 1 π P (cid:90) ∞−∞ dω (cid:48) Im[ χ ( ω (cid:48) )] ω (cid:48) − ω = − π δ ( ω ) . (B10)Now we proceed to the nonlinear effect in χ in the presence of electric fields by focusing on M ( ω (cid:48) ) = χ , ( ω (cid:48) , ω ) B ( ω (cid:48) ) E ( ω ) E ( − ω ) , (B11)which is a nonlinear modulation of χ ( ω (cid:48) ) proportional to an intensity of incident light | E ( ω ) | . From the dimensionanalysis, this modification is expected to take the form χ , ( ω (cid:48) , ω ) ∝ χ ( ω (cid:48) ) ω . (B12)We verify this by an explicit calculation in the following. The coefficient χ , ( ω (cid:48) , ω ) is obtained by computing fourpoint correlation functions of the current operator j x , which are illustrated in Fig. 3, and looking at the contributionsproportional to − q y . If we again use the relation, ∂ q y G = − G∂ q y G − G = Gj y G , this is written as χ , ( ω (cid:48) , ω ) = Re (cid:20) Q ( ω )( − iω )( iω ) (cid:21) (B13) Q ( i Ω) = − π ) (cid:90) d k (cid:88) iω n Tr[ A ( iω n , i Ω , i Ω (cid:48) ) + A ( iω n , i Ω , i Ω (cid:48) ) + A ( iω n , i Ω , i Ω (cid:48) )] , (B14)where A ( iω n , i Ω , i Ω (cid:48) ) = j x G ( iω ) j x G ( iω + i Ω) j x G ( iω ) j x G ( iω + i Ω (cid:48) ) j y G ( iω + i Ω (cid:48) ) j y G ( iω + i Ω (cid:48) ) , (B15) A ( iω n , i Ω , i Ω (cid:48) ) = j x G ( iω ) j x G ( iω − i Ω) j x G ( iω ) j x G ( iω + i Ω (cid:48) ) j y G ( iω + i Ω (cid:48) ) j y G ( iω + i Ω (cid:48) ) , (B16)and A ( iω n , i Ω , i Ω (cid:48) )= j x G ( iω ) j x G ( iω + i Ω) j x G ( iω + i Ω + i Ω (cid:48) ) j y G ( iω + i Ω + i Ω (cid:48) ) j y G ( iω + i Ω + i Ω (cid:48) ) j x G ( iω + i Ω (cid:48) )+ j x G ( iω ) j x G ( iω + i Ω) j x G ( iω + i Ω + i Ω (cid:48) ) j x G ( iω + i Ω (cid:48) ) j y G ( iω + i Ω (cid:48) ) j y G ( iω + i Ω (cid:48) )+ j x G ( iω ) j x G ( iω + i Ω) j x G ( iω + i Ω + i Ω (cid:48) ) j y G ( iω + i Ω (cid:48) + i Ω (cid:48) ) j x G ( iω + i Ω (cid:48) ) j y G ( iω + i Ω (cid:48) ) , (B17)The result for the imaginary part is given byIm[ χ ( ω )] = (cid:40) ω − ω ω (cid:48) +2 ωω (cid:48) + ω (cid:48) ω ω (cid:48) ( ω + ω (cid:48) ) , ( ω (cid:48) < ω ) − ω +4 ω ω (cid:48) − ωω (cid:48) − ω (cid:48) ωω (cid:48) ( ω + ω (cid:48) ) . ( ω (cid:48) ≥ ω ) (B18)In the limit of ω (cid:48) → ω , this reduces toIm[ χ , ( ω (cid:48) , ω )] = 116 ω ω (cid:48) , (B19)which is the form expected from the scaling analysis. This indicates that the real part of the orbital magneticsusceptibility behaves as Re[ χ ( ω (cid:48) )] = − π δ ( ω (cid:48) ) (cid:18) − ω | E ( ω ) | (cid:19) , (B20)which shows a nonlinear suppression of the orbital magnetic susceptibility with an applied electric field. jx(i Ω ,0) [E(i Ω ,0)]jx(-i Ω ,0) [E(-i Ω ,0)] jx(i Ω ’ ,q) [B(i Ω ’ ,q)]jx(-i Ω ’ ,-q) [M(i Ω ’ ,q)]i ω ,k i ω+ i Ω ,k i ω ,k i ω+ i Ω ’ ,k+q (a) jx(-i Ω ,0) [E(-i Ω ,0)]jx(i Ω ,0) [E(i Ω ,0)] jx(i Ω ’ ,q) [B(i Ω ’ ,q)]jx(-i Ω ’ ,-q) [M(i Ω ’ ,q)]i ω ,k i ω− i Ω ,k i ω ,k i ω+ i Ω ’ ,k+q (b) jx(i Ω ,0) [E(i Ω ,0)] jx(-i Ω ,0) [E(-i Ω ,0)]jx(i Ω ’ ,q) [B(i Ω ’ ,q)] jx(-i Ω ’ ,-q) [M(i Ω ’ ,q)]i ω ,k i ω+ i Ω ,k i ω+ i Ω+ i Ω ’ ,k+q i ω+ i Ω ’ ,k+q (c) FIG. 3. Diagrams for the nonlinear response χ , ( ω (cid:48) , ω ). Panels (a), (b), and (c) corresponds to A , A , and A in Eq. (B14),respectively. Appendix C: High harmonics generation
We study the high harmonic generations described by the nonlinear response, J ( N ω ) = χ N ( ω ) E ( ω ) N , (C1)where J ( N ω ) is the induced current with frequency
N ω , and E ( ω ) is the electric field of incident light. The coefficient χ N ( ω ) is given by an N -particle correlation function as χ N ( ω ) = Re (cid:20) Q ( ω )( − iω ) N (cid:21) , (C2) Q ( i Ω) = 1(2 π ) (cid:90) d k (cid:88) iω Tr (cid:20) j x G (cid:18) iω − i N −
12 Ω (cid:19) j x G (cid:18) iω − i N −
32 Ω (cid:19) . . . j x G (cid:18) iω + i N −
12 Ω (cid:19)(cid:21) . (C3)We notice that the coefficients χ N ( ω ) for even N vanish after integration over the angle of k because of the reflectionsymmetry. Here we study the third harmonic generation χ ( ω ) deduced from Q ( i Ω) = (cid:90) d k (cid:88) iω Tr (cid:20) j x G (cid:18) iω − i
32 Ω (cid:19) j x G (cid:18) iω − i
12 Ω (cid:19) j x G (cid:18) iω + i
12 Ω (cid:19) j x G (cid:18) iω + i
32 Ω (cid:19)(cid:21) . (C4)After integrating over the direction of k and summing over the Matsubara frequency, one obtains Q ( ω ) = 12 π (cid:90) kdk − ( k + 4 kω )16( k − ω )( k − ω )( k − ω ) . (C5)Since the imaginary part of Q ( ω ) is obtained from residues at k = ω , ω, ω (which are evaluated with ω shifted intothe upper half plane), the coefficient χ ( ω ) is given by χ ( ω ) = 1384 ω , (C6)which is consistent with the scaling law. Appendix D: Intensity dependent refractive index
The intensity dependent refractive index is related to the intensity dependent optical conductivity which is describedby J ( ω ) = χ (3) ( ω , ω ) E ( ω ) E ( ω ) E ( − ω ) . (D1)The coefficient χ (3) is obtained from a correlation function of four current operators as χ (3) ( ω , ω ) = Re (cid:20) Q ( ω , ω )( − iω ) ω (cid:21) , (D2)and Q ( i Ω , i Ω ) = 1(2 π ) (cid:90) d k (cid:88) iω Tr (cid:104) j x G ( iω + i Ω ) j x G ( iω ) j x G ( iω + i Ω ) j x G ( iω )+ j x G ( iω + i Ω ) j x G ( iω ) j x G ( iω − i Ω ) j x G ( iω )+ j x G ( iω + i Ω ) j x G ( iω + i Ω + i Ω ) j x G ( iω + i Ω ) j x G ( iω ) (cid:105) . (D3)After integrating over the direction of k and summing over Matsubara frequencies, one obtains Q ( ω , ω ) = 12 π (cid:90) kdk − k ( ω − k ) ( ω − k ) (4 k − ( ω − ω ) ) (4 k − ( ω + ω ) ) × (cid:104) k + 192 k (cid:0) ω + ω (cid:1) − k (cid:0) ω − ω ω + 7 ω (cid:1) + 4 k (cid:0) ω − ω ω − ω ω + ω (cid:1) + 5 ω ω (cid:0) ω − ω (cid:1) (cid:105) . (D4)The final result after the k integration reads χ (3) ( ω , ω ) = (cid:40) ω − ω ω ω ( ω + ω ) , ( ω > ω ) ω − ω ω ( ω + ω ) , ( ω ≤ ω ) (D5)The coefficient χ (3) ( ω , ω ) is a monotonically decreasing function of ω for fixed ω , where χ (3) (cid:39) / ω ω for ω (cid:29) ω . In the limit of ω → ω , one obtains χ (3) ( ω → , ω ) = 18 ω . (D6) [1] R. W. Boyd, Nonlinear optics (Academic press, London,2003).[2] H. Kishida, H. Matsuzaki, H. Okamoto, T. Manabe,M. Yamashita, Y. Taguchi, and Y. Tokura, Nature ,929 (2000).[3] H. Kishida, M. Ono, K. Miura, H. Okamoto, M. Izumi,T. Manako, M. Kawasaki, Y. Taguchi, Y. Tokura, T. To-hyama, K. Tsutsui, and S. Maekawa, Phys. Rev. Lett. , 177401 (2001).[4] J. McIver, D. Hsieh, H. Steinberg, P. Jarillo-Herrero, andN. Gedik, Nat. Nanotech. , 96 (2012).[5] A. Junck, G. Refael, and F. von Oppen, Phys. Rev. B , 075144 (2013).[6] N. Ogawa, M. S. Bahramy, H. Murakawa, Y. Kaneko,and Y. Tokura, Phys. Rev. B , 035130 (2013).[7] N. Ogawa, M. S. Bahramy, Y. Kaneko, and Y. Tokura,Phys. Rev. B , 125122 (2014).[8] S. M. Young and A. M. Rappe, Phys. Rev. Lett. ,116601 (2012).[9] I. Grinberg, D. V. West, M. Torres, G. Gou, D. M. Stein,L. Wu, G. Chen, E. M. Gallo, A. R. Akbashev, P. K. Davies, et al. , Nature , 509 (2013).[10] S. D. Ganichev and W. Prettl, Intense terahertz exci-tation of semiconductors (Oxford Univ. Press, Oxford,2006).[11] M. C. Hoffmann, N. C. Brandt, H. Y. Hwang, K.-L. Yeh,and K. A. Nelson, Appl. Phys. Lett. , 231105 (2009).[12] D. Turchinovich, J. M. Hvam, and M. C. Hoffmann,Phys. Rev. B , 201304 (2012).[13] M. Cornet, J. Degert, E. Abraham, and E. Freysz, J.Opt. Soc. Am. B , 1648 (2014).[14] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045(2010).[15] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. , 1057(2011).[16] A. K. Geim and K. S. Novoselov, Nat. Mat. , 183 (2007).[17] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109(2009).[18] M. Glazov and S. Ganichev, Physics Reports , 101(2014), high frequency electric field induced nonlinear ef-fects in graphene. [19] A. R. Wright, X. G. Xu, J. C. Cao, and C. Zhang, Ap-plied Physics Letters , 072101 (2009).[20] E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, andS. A. Mikhailov, Phys. Rev. Lett. , 097401 (2010).[21] S. Mikhailov, EPL (Europhysics Letters) , 27002(2007).[22] S. A. Mikhailov, Phys. Rev. B , 241301 (2014).[23] M. J. Paul, Y. C. Chang, Z. J. Thompson, A. Stickel,J. Wardini, H. Choi, E. D. Minot, B. Hou, J. A. Nees,T. B. Norris, and Y.-S. Lee, New Journal of Physics ,085019 (2013).[24] Z. Mics, K.-J. Tielrooij, K. Parvez, S. A. Jensen,I. Ivanov, X. Feng, K. M¨ullen, M. Bonn, and D. Turchi-novich, Nat. commun. , 7655 (2015).[25] X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B , 195424 (2008).[26] A. M. Essin, J. E. Moore, and D. Vanderbilt, Phys. Rev.Lett. , 146805 (2009).[27] D.-X. Qu, Y. S. Hor, J. Xiong, R. J. Cava, and N. P.Ong, Science , 821 (2010).[28] A. N. Redlich, Phys. Rev. D , 2366 (1984).[29] B. D´ora and R. Moessner, Phys. Rev. B , 165431(2010). [30] J. W. McClure, Phys. Rev. , 666 (1956).[31] M. Koshino and T. Ando, Phys. Rev. B , 085425(2007).[32] The high harmonic generation in graphene has been pre-viously studied by using the semiclassical equation [21]and by using the quantum perturbation theory [22]. Here,we apply our scaling arguments to the third harmonicgeneration and examine it by the quantum perturbationtheory.[33] D. Citrin, Optics letters , 554 (2001).[34] V. M. Edelstein, Solid State Commun. , 233 (1990).[35] C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang,M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, Z.-Q. Ji,Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S.-C.Zhang, K. He, Y. Wang, L. Lu, X.-C. Ma, and Q.-K.Xue, Science , 167 (2013).[36] J. G. Checkelsky, R. Yoshimi, A. Tsukazaki, K. S.Takahashi, Y. Kozuka, J. Falson, M. Kawasaki, andY. Tokura, Nat. Phys.10