Dominant Role of Two-Photon Vertex in Nonlinear Response of Dirac Materials
DDominant Role of Two-Photon Vertex in Nonlinear Response of Dirac Materials
Habib Rostami ∗ and Emmanuele Cappelluti Nordita, KTH Royal Institute of Technology and Stockholm University,Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden Istituto di Struttura della Materia (ISM), CNR, 34149 Trieste, Italy
Using a conserving Baym-Kadanoff approach, we present a fully compelling theory of nonlineardc response of a Dirac system to electric fields in the presence of disorder scattering. We showthat the nonlinear terms are strikingly ruled by the appearance of a dominant two-photon vertexwhich is absent at the bare level and finite even in the weak-coupling limit. Such two-photon vertexself-generation highlights the crucial role of the frequency and field dependence of the scatteringrates in the nonlinear regime. Our study reveals a novel many-body mechanism in the nonlinearresponse of Dirac materials whose effects are predicted to be observable.
Due to their linear dispersion, (cid:15) k ∼ ∣ k ∣ , and to the un-derlying chiral structure, Dirac materials show a varietyof exotic features that makes them a versatile platformfor theoretical investigations of new physics and for ap-plication purposes. Despite the complex physics, manyproperties of these materials are often rationalized us-ing concepts of non-interacting particle or semi-classicalmodel [1–4]. For instance, a standard transport model isconventionally applied for the dc conductivity in highly-doped graphene ( Boltzmann regime), where the mobilityis evaluated at the non-interacting level, and the interac-tions enter only through the effective parameter known astransport scattering rate Γ tr [5]. At odds with the abovescenario, there is a wide consensus that the quantum regime (low-energy transitions in undoped Dirac model)is much more complex and it might be significantly af-fected by many-body effects [6].The nonlinear electromagnetic response of Dirac ma-terials has attracted recently a considerable interest intwo [7–14] and three dimensions [15–20]. Widely investi-gated are the nonlinear optical properties and in partic-ular the appearing of four-wave mixing, nonlinear Kerreffect, second and third-harmonic-generation in single-layer graphene, with remarkable technological interest[21–30]. Peculiar of Dirac material is, due to the lineardispersion, the absence of the bare two-photon-electroncoupling, which should give rise to the so-called diamag-netic term. The lack of such term prompts several widelydebated issues, as the validity of optical sum rules [31–33]. Most of the theoretical descriptions of nonlinear ef-fects rely at the moment upon non-interacting analyses,or semiclassical approaches [23, 24, 34–36] where, in asimilar way as in the Boltzmann theory of linear response,the dominant transition processes resemble the ones ofthe non-interacting case and the scattering sources areaccounted through effective parameters as the scatteringrate Γ (or equivalently through the mean-free path l , thelifetime τ , etc.).In this Letter, we show that the a compelling analy- ∗ Electronic address: [email protected] sis of the many-body physics, beyond the semi-classicalapproaches, can drastically change the above scenario,pointing out that different physical processes can be re-sponsible for the relevant properties of the nonlinear dctransport. Analyzing the case of disorder scattering as abasilar benchmark example, we show how non-conservingphenomenological models of scattering intrinsically failand high-order vertex processes must properly taken intoaccount. More in particular, we show that, despite thebare diamagnetic two-photon vertex (TPV) being nullin Dirac materials, the many-body renormalized TPVis finite and relevant and it can play a dominant role.Our results, besides providing a consistent frameworkfor a proper analysis of nonlinear transport and opti-cal response in realistically interacting Dirac materials,open novel perspectives for understanding and predict-ing new functional properties of these complex promisingsystems.We consider the two-dimensional (2D) Dirac Hamilto-nian ˆ H k = ̵ hv ˆ σ ⋅ k − µ ˆ I , where µ is the bare chemicalpotential ruling the charge doping. For realistic purposeswe consider the paradigmatic case of graphene [37] whereˆ σ = ( τ ˆ σ x , ˆ σ y ) ,where ˆ σ i stand for the Pauli matrices inthe spinor space, and τ = ± stands for valley index in theBrillouin zone of graphene. In the dipole approximationthe light-matter interaction can be modelled by applyingthe minimal coupling transformation ̵ h k → ̵ h k + e A ( t ) where A ( t ) stands for an external vector potential. Thecorresponding electric field is given by E ( t ) = − ∂ t A ( t ) .Due to the linear dispersion, the electron-photon cou-pling does not present a diamagnetic (two-photon) bareterm but only the linear coupling: H light − matter = ̵ hev ∫ d r ˆ ψ † ( r ) ˆ σ ⋅ A ( t ) ˆ ψ ( r ) . (1)Without the loss of generality, we assume an electric fieldalong the y axis. As scattering source we consider long-range impurity centers with standard Born impurity cor-relations [38–40]. Within this framework we can writethe Born impurity self-energy in the complex frequencyspace: ˆΣ ( z ) = γ imp ∑ k ˆ G ( k , z ) where the Green’s func-tion follows ˆ G ( k , z ) = [ z − ˆ H k − ˆΣ ( z )] − . For isotropicscattering we get a diagonal self-energy in the spinor basis a r X i v : . [ c ond - m a t . m e s - h a ll ] J u l RPA = + RPARPA = + RPA RPA=+RPA RPA=+RPA = RPA=+RPA RPA=+RPARPA = + RPARPA = + RPA
RPA=+RPA R P A =+ R P A R P A = + R P A RPA = + RPA
RPA = + RPA RPA=+RPA R P A =+ R P A R P A = + R P A RPA = + RPA
RPA = + RPA R P A = + R P A R P A = + R P A RPA=+RPA +++ (a)
RPA = + RPARPA = + RPARPA = + RPA = + (b) = + R P A = + R P A R P A = + R P A R P A = + R P A R P A = + R P A R P A = + R P A R P A = + R P A (c) = + RPA = + RPA R P A = + R P A R P A = + R P A RPA = + RPA R P A = + R P A R P A = + R P A RPA = + RPA R P A = + R P A R P A = + R P A (d) RPA = + RPARPA = + RPA R P A = + R P A R P A = + R P A = (e) RPA = + RPA R P A = + R P A R P A = + R P A RPA = + RPARPA = + RPARPA = + RPA
RPA = + RPA R P A = + R P A R P A = + R P A RPA = + RPA R P A = + R P A R P A = + R P A = (f) ++ FIG. 1:
Diagrammatic representation of nonlinear e.m. response in Dirac materials. (a) nonlinear responseexpressed in terms of renormalized one-, two- and three-photon vertices; (b)-(d) self-consistent Bethe-Salpeter eqautions forone-, two- and three-photon vertices; (e)-(f) many-body definition of unrenormalized two- and three-photon vertices in terms oflower order vertices. Note that solid and wavy lines stand for fermion propagator and external photons, respectively. Dashedline indicate impurity interaction line. Void and filled circle stand for bare and renormalized one-photon vertex, respectively.Void and filled square stand for unrenormalized and renormalized two and three-photon vertices. as ˆΣ ( z ) = Σ ( z ) ˆ I . It is well known that under these con-ditions the impurity self-energy, as well the Coulomb andother scattering ones, depends intrinsically on the ultra-violet energy cut-off W representing the range of validityof the Dirac model. In order to provide a conserving ap-proach, this issue needs to be cured by means of a properregularization [41, 42]. As detailed in the SupplementaryMaterial (SM), we employ standard dimensional regular-ization leading to:Σ ( z ) = − U S ( z ) ln [ − W / S ( z )] , (2)where S ( z ) = z + µ − Σ ( z ) , and U is a dimensionlessparameter characterizing the strength of impurity scat-tering [38–40].Conserving approaches, based for instance on a Baym-Kadanoff derivation [43, 44], are fundamental in theo-retical physics to ensure that compelling results are ob-tained. This aim is particularly important in nonlinearresponse since an arbitrary selection of diagrams can eas-ily lead to spurious conclusions. The choice of the vector-potential gauge, within the paradigmatic Born impurityscattering we consider here, permits us an exact deriva-tion of self-consistent equations (see SM for details [45])for all the high order processes relevant in the third-order response function which is the leading nonlinearterm in centrosymmetric Dirac materials. The diagram-matic expression of the third-order response function is provided in Fig. 1 where, roughly speaking, empty sym-bols represent n -photon vertices ( n = ,
3) expressed interms of the renormalized lower-order vertices (Fig. 1e,f),whereas filled symbols represent the solution of a Bethe-Salpeter-like (BS) self-consistent resummation for a given n -photon vertex (Fig. 1b-d). Leaving aside the com-plexity of the self-consistent set of equations, few rele-vant things are worth to be underlined here. First ofall, we notice that an effective multi-photon coupling isinduced by the disorder scattering source even if it isabsent in the Hamiltonian at the bare (non-interacting)level (Fig. 1e,f). Second, that the relevance of each n -photon vertex is largely governed by the BS many-bodyresummation as depicted in Fig. 1b-d. This might lead toa reduction (screening) or to an enhancement of differentmulti-photon scattering depending on the Pauli structureof the corresponding photon vertex, as we discuss moreextensively later.The diagrammatic expressions in Fig. 1 represent infull generality the optical frequency-dependent third-order response function in Dirac materials, includingthird-harmonic generation, four-wave mixing, etc. Fora generic interaction, the effective solution of such cou-pled equations on the real-frequency axis is a formidabletask that does not allow for a practical solution. Thefocus on the isotropic disorder scattering is on the otherhand particularly suitable to investigate many-body ef-fects in nonlinear electromagnetic response since it allowsfor a set of equations in the Matsubara space which canbe generalized in a rigorous way on the real frequencyaxis, using the well-known procedure of multiple branchcuts in the complex frequency space. The derivation islengthly and cumbersome but compelling and it is sum-marized in the SM [45]. We consider first the dc transportlimit. Without loss of generality, it is possible to expressthe linear and the third-order dc conductivity in termsof two dimensionless quantities: σ ( ) dc = σ f ( µ Γ ( µ ) ; U ) , (3) σ ( ) dc = σ E [ t Γ ( µ ) ] f ( µ Γ ( µ ) ; U ) , (4)where µ is the effective chemical potential µ = µ − ReΣ ( ω = ) and Γ ( µ ) the scattering rate Γ ( µ ) = − ImΣ ( ω = ) , σ ∝ e / ̵ h the universal conductivity unitand E ∝ t / ea is a characteristic electric field scale de-termined by inter-atomic hopping energy t and by thelattice constant a [37].It is worth to stress again that Eqs. (3)-(4) are tied to-gether since they must descend in a compelling way froma common approximation for the self-energy. Hereto-fore, although many approaches for the self-energy havebeen discussed for the linear response, the third-order re-sponse has been analyzed only in the simplistic case ofa phenomenological constant scattering rate Γ ( µ ) = Γ.Since such phenomenological self-energy does not de-pend on the applied external field, the third-order re-sponse function reduces to the first “square” diagram ofFig. 1a dropping all the vertex renormalization processes,i.e. replacing the filled circles with empty ones (= bareelectron-photon coupling). A similar scheme can as wellbe employed for the linear response. Under these ultra-simplified conditions, one can see that the linear andthird-order dc transport depend uniquely on the semi-classical parameter x = µ / Γ ( µ ) , i.e. f ( x ; y ) = f ( x ) , f ( x ; y ) = f ( x ) . An analytical expression for the func-tions f ( x ) , f ( x ) is obtained in the SM [45]. In partic-ular, in the Boltzmann regime one gets results f ( ∞ ) ≈ µ / π Γ, f ( ∞ ) ≈ − π Γ / µ , implying that nonlinear ef-fects lead to a reduction of the dc conductivity in theBoltzmann regime. A similar analysis is performed in thequantum regime, giving [45] f ( ) = / π , f ( ) = / enhance-ment of the dc conductivity in the quantum regime.The above predictions, based on the phenomenologi-cal model of a constant scattering rate Γ, are challengedwhen many-body effects are computed in a compellingconserving scheme. In Fig. 2 we show the characteristicdc transport [46] function f as a function of the semi-classical parameter x = µ / Γ ( µ ) , from the extreme quan-tum limit ( x ≪
1) to the Boltzmann regime ( x ≫ W , the quantumlimit can be conveniently investigated by fixing U andvarying µ , whereas the Boltzmann regime is more eas- - - f ( x ; U ) - - x f ( x ; U ) - - - - - - x f ( x ; U ) (a) x x x U=0.06U=0.12 x ⇥ µ = 2[THz] µ = 26[THz] (b) x U=0.06U=0.12
FIG. 2: Characteristic nonlinear dc transport function f as a function of the dimensionless parameter x = µ / Γ ( µ ) forthe many-body conserving scheme. Also shown is f for thephenomenological model (dashed line). Curves of f vs. x are plotted for different U ’s in the quantum regime (panel a),and for different µ ’s in the Boltzmann regime in THz unit(panel b). ily spanned by fixing µ and varying U . Let’s discussfirst the Boltzmann regime (Fig. 2b). We notice thata compelling many-body analysis recovers qualitatively(but with an increase in the magnitude of factor 10) thepredictions of the phenomenological constant-Γ modelwith f ( x, y ) ≈ f ( x ) ∝ − / x in the Boltzmann regime( x ≫ x ∗ ≈ x ≪ x ∗ (Fig. 2a), f shows a significant depen-dence on µ , signalizing that the nonlinear dc transportproperties are no more governed uniquely by the semi-classical parameter x but that the detailed value of µ (or conversely, of U ) starts playing a relevant role. Somestriking things are worth being pointed out: ( i ) counter-intuitively, the third-order contribution to the dc trans-port appears to be magnified approaching the clean limit U →
0; ( ii ) there are two isosbestic points (i.e. x -pointswhere f does not depend on U ) coinciding in a verygood approximation with the zeroes of the f function;( iii ) whereas the phenomenological constant-Γ modellingpredicts a a well-determined positive sign of the nonlin-ear dc correction in the quantum regime (implying anincrease of the total conductivity), the sign of the non-linear terms of the full conserving many-body theory inthe quantum regime is not univocally determined, pre-senting a positive region in the crossover range and anegative sign in the extreme quantum limit.We can rationalize points ( i )-( ii ) by assuming that inthe quantum regime the nonlinear characteristic dc trans-port function f ( x, U ) can be factorized as [45]: f ( x ; U ) ≈ C ( x ) U γ , (5)(with γ > U of the interactionrules governs the intensity of the third-order dc trans-port, while the semiclassical parameter x seems to dic-tate the sign of the third-order correction. To assess themeaningfulness of such description we plot in Fig. 3a ina log-log scale the absolute value of the function f ( x, U ) versus U for a representative case x = .
03 in the quan-tum regime. We find a perfect agreement with a scal-ing behavior f ( x, U ) ∝ / U γ with γ slightly smallerthan 2 ( γ = .
82) signalizing that in the clean limit thethird-order dc transport is expected to be dominant withrespect to the linear one. As detailed in [45], a similaranalysis is valid in the whole quantum regime.In order to gain a full understanding of these novelfeatures, we analyze separately in Fig. 3a,b the rele-vance of each family of diagrams contributing to the to-tal third-order conductivity as depicted in Fig. 1a. Wecan thus realize that the contribution of the conventional“square” diagram (which is the only one present in thenon-interacting case and for the phenomenological damp-ing model), is essentially marginal, as well as the contri-bution of the last “bubble” associated with the renor-malized three-photon vertex. The dominant role is in-stead played by the “triangle” diagrams containing therenormalized two -photon vertex. A quantitative analysisshows that each family of diagrams obeys Eq. (5) withan approximately integer exponent (i.e. γ ■ ≈ γ ▲ ≈ γ • ≈ γ ≈ γ ▲ ). The self-consistent BS renormalization of the TPV (Fig. 1c) isa crucial ingredient in such novel scenario. This can beassessed in Fig. 3a,b where one can see that, once ne-glected the BS renormalization, the contribution of thetriangle diagrams results to be of the same order (evensmaller) of that of the conventional square diagram. Thedominant role of the TPV renormalization appears evenmore evident by investigating the scaling of the charac-teristic third-order dc transport function f ( x, U ) versus U . As depicted in Fig. 3a, once replaced the BS renor-malized TPV (Fig. 1c) with the “bare” one (Fig. 1e), thethird-order dc conductivity scales as 1 / U ( γ △ ≈ ∼ / U that diverges in the clean limit. The n -photonvertex matrix structure, which reads ˆΛ n = ( − ev ˆ σ y ) n Λ n ,plays a crucial role in the relevance of the BS renormal-ization effect. The impressive effect is peculiar of theTPV renormalization Λ = Λ ( ) /[ − U X ] and doesnot appear in the BS renormalization of the one- three-photon vertex (Λ n = Λ ( ) n /[ − U X n ] , with n = ,
3) [45].This different impact can be traced down to the differentstructure in the Pauli space. As detailed in [45], we getindeed X = X ∝ Tr [ ˆ σ y ˆ G ˆ σ y ˆ G ] , X ∝ Tr [ ˆ G ˆ G ] . In thequantum regime, one can thus show that in the dc limit U X ≈ U X ∝ U , whereas U X ≈ + O ( U ) , so thatresulting in an effective divergence of Λ / Λ ( ) in the dclimit at zero temperature and in the clean limit ( U → SquareTrianglesTriangles, X = - PhotonFull - - - - - μ [ THz ] σ d c ( ) E / σ d c ( ) (b) OOVVXX O x = V x = X x = U | f t r i ang l e s ( x ; U ) | - U | f ( x ; U ) | (a) . NegativePositiveNegative (c) FIG. 3:
Nonlinear conductivity of graphene at zerotemperature. (a) Log-log scale plot for the absolute valueof the characteristic third-order transport function f ( µ, U ) versus U at x = .
03. Different lines correspond to the indi-vidual contribution of diagrams in Fig. 1a as mentioned theshared plot-legend in panel (b). (b) f ( µ, U ) as a functionof the chemical potential for U = .
11. Similar to the panel(a), different curves correspond to different diagram’s contri-bution as mentioned the plot-legend. (c) Colormap plot for g = ∣ σ ( ) dc E / σ ( ) dc ∣ factor with E = / nm versus chemi-cal potential µ and relaxation rate Γ ( µ ) in the full conservingmodels. The sign of σ ( ) dc is written on the plot where two sign-switch borders are highlighted by dashed red lines. Green andblue dashed lines stand for the contour lines with g = g = .
1, respectively. Similar colormap plot for the constant-Γmodel is given in the SM [45]. photon renormalization factor, [45] Q ( z , z ) = − U X ( z , z ) = S ( z ) − S ( z ) z − z , (6)where X ( z , z ) ∝ ∑ k Tr [ ˆ G ( k , z ) ˆ G ( k , z )] and where z and z are the electronic frequencies in the complexplane. Particularly enlightening is the analysis of theretarded-retarded (RR) channel. In the dc limit ( ω → Q RR2 ( (cid:15), (cid:15) + ω ) is simply given bylim ω → Q RR2 ( (cid:15), (cid:15) + ω ) = dS ( (cid:15) ) d(cid:15) = S ( (cid:15) ) U S ( (cid:15) ) + µ + (cid:15) , (7)where ω is the photon energy and (cid:15) is the electronic en-ergy from the Fermi surface. For low-energy excitations (cid:15) = Q RR2 = S ( )/[ U S ( ) + µ ] . TheBoltzmann regime is achieved as µ ≫ U S ( ) . In theclean limit U → S ( ) = µ and Q RR2 = µ ≪ U S ( ) , and we get Q RR2 = / U , leadingthus to a huge enhancement (divergence) in the cleanlimit. A similar behavior Q ∝ / U appears also in theretarded-advance (RA) channel, although is a more deli-cate way. In the dc limit, we find indeed the leading termin Q RA2 ( (cid:15), (cid:15) + ω ) ≈ − i ( (cid:15) )/ ω which shows a divergenceas a function of the photon energy ω →
0. Once pluggedthis behavior in the response function associated withthe “triangle” diagrams containing the renormalized BStwo-phonon vertex, χ ren . triangles ( ω ) ∼ Γ ( (cid:15) ) χ unren . triangles ( ω )/ ω the divergence is ω implies that χ unren . triangles ( ω ) must beexpanded to a higher order in ω , involving the second derivative S ′′ ( ) = − /[ U Γ ( )] and higher orders. Asa net result, the divergence 1 / ω in Q RA2 ( (cid:15), (cid:15) + ω ) is re-flected in a consequent divergence − / U in the responsefunction, similarly as for the RR channel, but with a negative sign. The balance between RR and RA termsdetermines the change of sign of the third-order dc con-ductivity as a function of x n the quantum regime. Wemust stress that the huge enhancement of the third-orderdc transport is governed by the dominant role of the two-phonon vertex renormalization. Since Λ ( ) scales as U and 1 /[ − U X ] scales as 1 / U , such enhancement canbe regarded as TPV self-generation (Λ ≠
0) which sur-vives in the weak-coupling (clean) limit U → E = g = ∣ σ ( ) dc E / σ ( ) dc ∣ , in the physical space of the effective chem-ical potential µ and scattering rate Γ ( µ ) , as they can beobtained directly in an experimental way. The Boltz-mann regime corresponds thus to the right-lower cornerwhereas the extreme quantum regime ( x →
0) is recov-ered in the left-upper corner. As noticed before, at oddswith the predictions of the phenomenological model, wefind that the third-order conductivity σ ( ) dc is negative notonly in the Boltzmann regime, but also in the quantum regime. Note that the zeroes of the third-order dc con-ductivity, see x and x in Fig. 2, appear in this plot asstraight dashed lines. This is a consequence of the fac-torizable expression for the characteristic third-order dctransport function as shown in Eq. (5). We mark withtiny dotted in this plot the regions where the third-orderterms start to be relevant g ≈ . g ≈ µ → ∞ , one should correspondingly need infinitesimally smallΓ → µ → E max above which third-order corrections to the dc transport become of the sameorder of the linear term, ∣ σ ( ) dc ∣ E ∼ σ ( ) dc . In the phe-nomenological constant-Γ model we obtain E max = αE with α ≈ . × µ Γ / t in the Boltzmann regime and α ≈ . × Γ / t in the quantum limit. These values canbe compared with the estimates for the full conservingtheory that gives α ≈ . × µ Γ ( µ )/ t in the Boltzmannregime and α ≈ . U × Γ ( ) / t in the quantum limitwith U = / ln [ W / Γ ( )] . In Ref. [48] a roughly con-stant value Γ ≈
15 meV was estimated in the wide range µ ∼ −
200 meV. With these values the phenomenologicalmodel would estimate a critical field E max ≈ . µ ≈
200 meV and E max ≈ . µ V/nm for µ = µ ∗ ≈
150 meV,whereas the full conserving theory predicts E max ≈ . µ ≈
200 meV and E max ≈ .
051 mV/nm for µ =
0, in a more observable range.In conclusion, in this Letter, we have presented a fullyconserving theory of nonlinear transport response in 2DDirac materials. Our results show that the previous anal-yses in literature, based on phenomenological scatteringmodels, can be qualitatively (but not quantitatively) reli-able in the Boltzmann regime but they completely fail inthe quantum regime. We have shown that, in a wide re-gion of the phase diagram, close to the neutral point, thenonlinear dc transport response is dominated by novelphysical processes where the two-photon vertex, absentin the bare Dirac Hamiltonian, plays a relevant role. Itshould be furthermore stressed that our results, focusedon the dc limit, implies that the current knowledge aboutthe nonlinear optical response in the terahertz regimeshould be deeply revised. Our work opens new scenariosfor a deep understanding of the electromagnetic responseof Dirac systems, whose relevance ranges from condensedmatter to high-energy physics.
Acknowledgements–.
H.R. acknowledges the supportfrom the Swedish Research Council (VR 2018-04252). [1] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. ,109 (2009), URL https://link.aps.org/doi/10.1103/RevModPhys.81.109 .[2] T. Wehling, A. Black-Schaffer, and A. Bal-atsky, Advances in Physics , 1 (2014),https://doi.org/10.1080/00018732.2014.927109, URL https://doi.org/10.1080/00018732.2014.927109 .[3] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev.Mod. Phys. , 015001 (2018), URL https://link.aps.org/doi/10.1103/RevModPhys.90.015001 .[4] M. Katsnelson, Graphene: Carbon in Two Dimen-sions (Cambridge University Press, 2012), ISBN9780521195409. [5] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi,Rev. Mod. Phys. , 407 (2011), URL https://link.aps.org/doi/10.1103/RevModPhys.83.407 .[6] N. M. R. Peres, Rev. Mod. Phys. , 2673 (2010), URL https://link.aps.org/doi/10.1103/RevModPhys.82.2673 .[7] E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko,and S. A. Mikhailov, Phys. Rev. Lett. , 097401(2010), URL https://link.aps.org/doi/10.1103/PhysRevLett.105.097401 .[8] H. Zhang, S. Virally, Q. Bao, L. Kian Ping, S. Mas-sar, N. Godbout, and P. Kockaert, Optics Letters ,1856 (2012), URL http://ol.osa.org/abstract.cfm?URI=ol-37-11-1856 .[9] N. Kumar, J. Kumar, C. Gerstenkorn, R. Wang, H.-Y. Chiu, A. L. Smirl, and H. Zhao, Phys. Rev. B , 121406 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.87.121406 .[10] R. I. Woodward, R. T. Murray, C. F. Phelan, R. E. P.de Oliveira, T. H. Runcorn, E. J. R. Kelleher, S. Li, E. C.de Oliveira, G. J. M. Fechine, G. Eda, et al., 2D Mate-rials , 011006 (2016), URL https://doi.org/10.1088%2F2053-1583%2F4%2F1%2F011006 .[11] G. Soavi, G. Wang, H. Rostami, D. G. Purdie,D. De Fazio, T. Ma, B. Luo, J. Wang, A. K. Ott,D. Yoon, et al., Nature Nanotechnology , 583 (2018),URL https://doi.org/10.1038/s41565-018-0145-8 .[12] H. A. Hafez, S. Kovalev, J.-C. Deinert, Z. Mics, B. Green,N. Awari, M. Chen, S. Germanskiy, U. Lehnert, J. Te-ichert, et al., Nature , 507 (2018), URL https://doi.org/10.1038/s41586-018-0508-1 .[13] G. Soavi, G. Wang, H. Rostami, A. Tomadin, O. Balci,I. Paradisanos, E. A. A. Pogna, G. Cerullo, E. Lidorikis,M. Polini, et al., ACS Photonics , 2841 (2019), URL https://doi.org/10.1021/acsphotonics.9b00928 .[14] Q. Ma, S.-Y. Xu, H. Shen, D. MacNeill, V. Fatemi, T.-R.Chang, A. M. Mier Valdivia, S. Wu, Z. Du, C.-H. Hsu,et al., Nature , 337 (2019), URL https://doi.org/10.1038/s41586-018-0807-6 .[15] Q. Ma, S.-Y. Xu, C.-K. Chan, C.-L. Zhang, G. Chang,Y. Lin, W. Xie, T. Palacios, H. Lin, S. Jia, et al., Na-ture Physics , 842 (2017), URL https://doi.org/10.1038/nphys4146 .[16] F. de Juan, A. G. Grushin, T. Morimoto, and J. E.Moore, Nature Communications , 15995 (2017), URL https://doi.org/10.1038/ncomms15995 .[17] H. Rostami and M. Polini, Phys. Rev. B ,195151 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.195151 .[18] J. Ma, Q. Gu, Y. Liu, J. Lai, P. Yu, X. Zhuo, Z. Liu,J.-H. Chen, J. Feng, and D. Sun, Nature Materi-als , 476 (2019), URL https://doi.org/10.1038/s41563-019-0296-5 .[19] B. Cheng, N. Kanda, T. N. Ikeda, T. Matsuda,P. Xia, T. Schumann, S. Stemmer, J. Itatani, N. P.Armitage, and R. Matsunaga, Phys. Rev. Lett. ,117402 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.117402 .[20] J. Cheng, J. Sipe, and S. Wu, arXiv preprintarXiv:2005.13693 (2020), URL https://arxiv.org/pdf/2005.13693 .[21] S. A. Mikhailov, Europhysics Letters (EPL) , 27002(2007), URL https://doi.org/10.1209%2F0295-5075%2F79%2F27002 . [22] J. L. Cheng, N. Vermeulen, and J. E. Sipe, New Journalof Physics , 053014 (2014), URL https://doi.org/10.1088%2F1367-2630%2F16%2F5%2F053014 .[23] J. L. Cheng, N. Vermeulen, and J. E. Sipe, Phys. Rev.B , 235320 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.91.235320 .[24] S. A. Mikhailov, Phys. Rev. B , 085403 (2016),URL https://link.aps.org/doi/10.1103/PhysRevB.93.085403 .[25] H. Rostami and M. Polini, Phys. Rev. B ,161411 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.161411 .[26] H. Rostami, M. I. Katsnelson, and M. Polini, Phys. Rev.B , 035416 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.95.035416 .[27] S. A. Mikhailov, Phys. Rev. B , 115416 (2019),URL https://link.aps.org/doi/10.1103/PhysRevB.100.115416 .[28] A. Principi, D. Bandurin, H. Rostami, and M. Polini,Phys. Rev. B , 075410 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.99.075410 .[29] H. A. Hafez, S. Kovalev, K.-J. Tielrooij, M. Bonn,M. Gensch, and D. Turchinovich, Advanced Optical Ma-terials , 1900771 (2020), URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adom.201900771 .[30] H. Rostami and V. Juriˇci´c, Phys. Rev. Research ,013069 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.013069 .[31] S. P. Goldman and G. W. F. Drake, Phys. Rev. A , 2877 (1982), URL https://link.aps.org/doi/10.1103/PhysRevA.25.2877 .[32] R. Cenni, Nuclear Physics A , 605 (2001),ISSN 0375-9474, URL .[33] J. Sabio, J. Nilsson, and A. H. Castro Neto, Phys. Rev.B , 075410 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.78.075410 .[34] D. E. Parker, T. Morimoto, J. Orenstein, and J. E.Moore, Phys. Rev. B , 045121 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.99.045121 .[35] Y. Shen, The Principles of Nonlinear Optics (Wiley, NewYork, 1984).[36] R. Boyd,
Nonlinear Optics (Academic Press, Cambridge,2008).[37] For numerical calculations we consider a two-dimensionalDirac material with realistic parameters for graphene. Inparticular we consider a honeycomb lattice with inter-atom distance b = .
42 ˚A, with lattice parameter a = √ b = .
46 ˚A and nearest-neighbor hopping t ≈ v = √ t a / ̵ h ≈ m/s. The size of the unit-cell reads thus S cell = √ a / N s = N v =
2. In order to preserve the size of the Brilloiunzone V BZ = π / √ a , we introduce a momentum cut-off k c = .
09 ˚A − , so that πk c = V BZ , defining thus anultra-violet energy cut-off W = ̵ hvk c = . σ = e / ̵ h and E = πt /√ ea ≈ . / nm is an ultra-strong charac-teristics electric field (see Supplementary Materials [45]).[38] N. H. Shon and T. Ando, Journal of thePhysical Society of Japan , 2421 (1998),https://doi.org/10.1143/JPSJ.67.2421, URL https://doi.org/10.1143/JPSJ.67.2421 .[39] A. Tsuneya, Y. Zheng, and S. Hidekatsu, Journal of the Physical Society of Japan , 1318 (2002),https://doi.org/10.1143/JPSJ.71.1318, URL https://doi.org/10.1143/JPSJ.71.1318 .[40] H. Rostami and E. Cappelluti, Phys. Rev. B ,054205 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.054205 .[41] G. Leibbrandt, Rev. Mod. Phys. , 849 (1975), URL https://link.aps.org/doi/10.1103/RevModPhys.47.849 .[42] M. E. Peskin, An Introduction To Quantum Field Theory (CRC Press; 1 edition, 1995).[43] G. Baym and L. P. Kadanoff, Phys. Rev. ,287 (1961), URL https://link.aps.org/doi/10.1103/PhysRev.124.287 .[44] L. P. Kadanoff and G. Baym,
Quantum Statistical Me-chanics (W.A. Benjamin; 1st edition, 1962).[45] Supplementary Materials.[46] We estimate the third-order dc conductivity as σ ( ) dc = lim ω → σ ( ) THG ( ω ) . In Supplementary Materials [45] weplot the frequency dependence of the third-harmonicgeneration (THG) response function σ ( ) THG ( ω ) for µ =
17 THz and few representative values of U .[47] The condition g =
1, in the regions where the third-orderdc conductivity σ ( ) dc is negative should not be regardedas an onset of negative total dc conductivity, rather asa sign that the expansion at the third-order in E startsto be a poor approximation and higher order terms inpowers of E must be included in the analysis.[48] J. Horng, C.-F. Chen, B. Geng, C. Girit, Y. Zhang,Z. Hao, H. A. Bechtel, M. Martin, A. Zettl, M. F. Crom-mie, et al., Phys. Rev. B , 165113 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.83.165113 .[49] G. D. Mahan, Many-Particle Physics (Springer US; NewYork; 3 edition, 2000).
Supplemental Materials:“Dominant Role of Two-Photon Vertex in Nonlinear Response of Dirac Materials”
Habib Rostami, Emmanuele Cappelluti Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden Istituto di Struttura della Materia, CNR, 34149 Trieste, Italy
Contents
References S1. Ultra-violet cut-off and dimensional regularization S2. Baym-Kadanoff derivations
S3. Analytical expressions for nonlinear conductivity
S4. Linear conductivity S5. Nonlinear dc conductivity f ( x ) function in the constant-Γ model 20B. Numerical evaluation of f ( x, U ) function in the full quantum theory 21 S6. Analytical continuation for the third order response function S7. Numerical evaluation of the nonlinear optical and dc conductivities S1. ULTRA-VIOLET CUT-OFF AND DIMENSIONAL REGULARIZATION
The introduction of a high-energy (ultra-violet) cut-off is an unavoidable requirement of Dirac models. There ishowever a relative large degree of freedom in the way how to introduce it, and particular care is needed in order toavoid spurious results and to preserve physical consistencies, like Ward’s identities, and gauge invariance. Dimensionalregularization has proven to be a formidable tool to ensure that physical correctness is preserved [41, 42]. Here weshow how such approach provides a consistent framework for evaluating self-energy and susceptibilities. In our workwe employ dimensional regularization that endures gauge invariance. As a benchmark example, and for the sake ofsimplicity, we consider the evaluation of the disorder self-energy which displays a primary diverging integral.We consider scattering on local impurity centers with density n imp and potential V imp ( r ) = ∑ i V i δ ( r − R i ) where R i are the coordinates of the lattice sites. We assume standard Born impurity correlations as ⟨ V imp ( r )⟩ = V ( , ) = ⟨ V imp ( r ) V imp ( r )⟩ imp = γ imp δ ( r − r ) . (S1)Note that the average ⟨ . . . ⟩ imp is meant over all the impurity configurations and γ imp = n imp V in which n imp = N imp / N cell (number of impurity centers per number of unit-cells) stands for the density of scattering centers and V imp is the average strength of the scattering potential energy. The lowest-order self-consistent Born self-energy readsˆΣ ( ik n ) = γ imp ∑ k ˆ G ( k , ik n ) = γ imp S ∫ d k ( π ) ˆ G ( k , ik n ) . (S2)Note that S = N cell S cell is the system area with S cell being the unit-cell area. Due to the isotropic impurity scatteringthe self-energy spinor structure is trivial as ˆΣ ( ik n ) = Σ ( ik n ) ˆ I and therefore the Green’s function can be explicitlywritten as follows ˆ G ( k , ik n ) = [ z + µ − Σ ( z )] ˆ I + ̵ hv ˆ σ ⋅ k [ z + µ − Σ ( z )] − ( ̵ hvk ) (S3)where z = ik n . Accordingly, we find Σ ( z ) = γ imp N cell S cell ( ̵ hv ) ∫ d (cid:96) ( π ) S ( z ) S ( z ) − (cid:96) (S4)where S ( z ) = z + µ − Σ ( z ) . In arbitrary D dimensions, we haveΣ ( z ) = γ imp N cell S cell ( ̵ hv ) D ∫ d D (cid:96) ( π ) D S ( z ) S ( z ) − (cid:96) . (S5)Note that the above integral in D dimensions can be solved in terms of Euler’s Gamma-function, Γ E ( z ) , by utilizingthe following identity [42] ∫ d D (cid:96) ( π ) D ( (cid:96) + ∆ ) n = ( π ) D Γ E ( n − D ) Γ E ( n ) ( ) n − D . (S6)Therefore, we find Σ ( z ) = − U S ( z ) ( π ) D − ( ̵ hv ) D − Γ E ( − D ) ( − S ( z ) ) D − (S7)in which [37] U = γ imp N cell S cell π ( ̵ hv ) = N imp π √ ( V imp t ) . (S8)Now we set D = d − (cid:15) where d = (cid:15) →
0. Note that Γ E ( (cid:15) / ) ≈ / (cid:15) for (cid:15) → (cid:15) → ( X ) − (cid:15) / (cid:15) / = lim (cid:15) → (cid:15) − ln [ X ] = ln [ W ] − ln [ X ] = ln [ W X ] . (S9) U (cid:1) U (cid:1) U (cid:1) U (cid:1) U (cid:1) - - - - - μ [ THz ] Γ ( μ ) [ T H z ] U (cid:1) U (cid:1) U (cid:1) U (cid:1) U (cid:1) - - - - - - μ [ THz ] Δ ( μ ) [ T H z ] U (cid:1) U (cid:1) U (cid:1) U (cid:1) U (cid:1) - - - - - μ [ THz ] Γ ( μ ) [ T H z ] (a) (b) (c) FIG. S1:
Numerical result for the imaginary and real parts of the self-energy. (a) The imaginary part of theself-energy at the Fermi surface Γ = − Im [ Σ ( ω = )] is shown versus bare chemical potential µ . (b) The imaginary part ofthe self-energy at the Fermi surface ∆ = Re [ Σ ( ω = )] is shown versus bare chemical potential µ . (c) Γ is shown versusrenormalized chemical potential µ = µ − ∆ ( µ ) . Different curves correspond to different values of U . (cid:15) → (cid:15) ≡ ln [ W ] where W is the ultra-violet energy cut-off. Eventually, we obtainthe following self-consistent formula for the self-energyΣ ( z ) = − U S ( z ) ln [ − W S ( z ) ] . (S10)After solving the above relation self-consistently, the real, ∆ = Re [ Σ ( ω = )] , and imaginary, Γ = − Im [ Σ ( ω = )] ,parts of the self-energy at the Fermi surface are depicted in Fig. S1. For the undoped regime and at the Dirac point,we have Σ ( ) = − i Γ ( ) where Γ ( ) = W exp { − U } . (S11)The above procedure of dimensional regularization is employed in similar way in the evaluation of other momentumintegrals in this study. S2. BAYM-KADANOFF DERIVATIONS
Within the lowest-order self-consistent Born approximation, the self-energy correction induced by elastic impurityscattering reads ˆΣ ( , ) = V ( , ) ˆ G ( , ) . (S12)We use a shorthand notation 1 ≡ ( r , t ) for the space-time coordinate. Using Dyson recursive relation, the fullfield-dependent and interacting Green’s function is given in terms a field-dependent self-energy, ˆΣ, and a bare Green’sfunction, ˆ G , as follows ˆ G ( , ′ ; A ) = ˆ G ( , ′ ; A ) + ∫ ¯2 , ¯3 ˆ G ( , ¯2; A ) ˆΣ ( ¯2 , ¯3; A ) ˆ G ( ¯3 , ′ ; A ) (S13)or equivalently we have ˆ G − ( , ′ ; A ) = ˆ G − ( , ′ ; A ) − ˆΣ ( , ′ ; A ) . (S14)The inverse of bare Green’s function readsˆ G − ( , ′ ; A ) = [ i∂ t − v ˆ σ ⋅ ( − i ̵ h ∇ + e A ( ))] δ ( − ′ ) . (S15)We assume an external gauge field along y axis, A ( ) = A ( ) ˆ y . The one-photon current vertex is given in terms ofvariational derivative of bare Green’s function versus the gauge field:ˆΛ ( ) ( , ′ ; 1 ′′ ) = δ ˆ G − ( , ′ ) δA ( ′′ ) (cid:187)(cid:187)(cid:187)(cid:187)(cid:187)(cid:187) A → = − ev ˆ σ y δ ( − ′′ ) δ ( − ′ ) . (S16)The thermodynamic physical current, J ( A ) = J ( A ) ˆ y , in Dirac systems reads J ( A ) = − i ∫ ′ , ′′ tr [ ˆΛ ( ) ( , ′ ; 1 ′′ ) ˆ G ( , ′+ ; A )] . (S17)Note that 1 ′+ ≡ ( r ′ , t ′ = t + + ) and “tr” stands for the “trace” operation over all spinor indexes i.e. tr [ ˆ A ˆ B ] = ∑ ss ′ [ A ss ′ B s ′ s ] . Note the field operator ˆ ψ H ( r , t ) in the Heisenberg picture of ˆ ψ ( r ) in the basis of full Hamiltonian H which contains kinetics, light-matter and many-body interaction terms. The Baym-Kadanoff (or contour) Green’sfunction [43, 44] follows ˆ G ( , ′ ; A ) = − i ⟨ T [ ˆ ψ H ( ) ˆ ψ † H ( ′ )]⟩ . (S18)where ⟨ . . . ⟩ stands for the thermodynamical average and T is for the time-ordering operation.1 A. Conserving linear response theory in Dirac systems
We define the linear susceptibility as follows χ ( ) ( , ) = δJ ( A ) δA ( ) (cid:187)(cid:187)(cid:187)(cid:187)(cid:187)(cid:187) A → , (S19)Using Eq. (S17), the linear susceptibility reads χ ( ) ( , ) = − i ∫ ′ , ′′ tr [ ˆΛ ( ) ( , ′ ; 1 ′′ ) δ ˆ G ( , ′+ ; A ) δA ( ) ] A → (S20)Using ˆ G ˆ G − =
1, we have δ ˆ G ( , ′ ; A ) δA ( ) = − ∫ ¯2 , ¯3 ˆ G ( , ¯2; A ) δ ˆ G − ( ¯2 , ¯3; A ) δA ( ) ˆ G ( ¯3 , ′ ; A ) (S21)Using Eq. (S14), we find δ ˆ G − ( , ′ ; A ) δA ( ) = δ ˆ G − ( , ′ ; A ) δA ( ) − δ Σ ( , ′ ; A ) δA ( ) (S22)Since the self-energy depends on the external potential only through its dependance on the Green’s function, we canwrite down δ ˆΣ ( , ′ ; A ) δA ( ) = ∫ ¯3 , ¯4 ˆΞ ( , ′ ; ¯3 , ¯4; A ) δ ˆ G ( ¯3 , ¯4; A ) δA ( ) (S23)We define Bethe-Salpeter kernel as bellow ˆΞ ( ,
2; 3 , A ) = δ ˆΣ ( , A ) δ ˆ G ( , A ) (S24)Note that ˆ G ( , ′ ; A )∣ A → = ˆ G ( , ′ ) and ˆΞ ( ,
2; 3 , A )∣ A → = ˆΞ ( ,
2; 3 , ) . For our self-energy model within theself-consistent Born approximation given in Eq. (S12), we haveˆΞ ( ,
2; 3 , ) = V ( , ) δ ( − ) δ ( − ) (S25)Therefore, using Eq. (S21) the self-consistent Bethe-Salpater relation for the one-photon vertex function followsˆΛ ( , ′ ; 2 ) = ˆΛ ( ) ( , ′ ; 2 ) + V ( , ′ ) ∫ ¯5 , ¯6 ˆ G ( , ¯5 ) ˆΛ ( ¯5 , ¯6; 2 ) ˆ G ( ¯6 , ′ ) (S26)Note that the dressed one-photon vertex function are defined asˆΛ ( , ′ ; 2 ) = δ ˆ G − ( , ′ ; A ) δA ( ) (cid:187)(cid:187)(cid:187)(cid:187)(cid:187)(cid:187) A → (S27)For the sake of simplicity, we extend the definition of space-time parameter to include also spinor indexes as ¯1 = ( r ¯1 , t ¯1 , s ¯1 ) and, for instance, we can drop “ˆ” symbol in ˆΛ ˆ G instead use Λ G since the spinor multiplication is takeninto account when we replace ∫ d ¯1 ∑ s ¯1 → ∫ d ¯1. Moreover we use shorthand notation for C = ∫ d ¯1 A ( . . . , ¯1 ) B ( ¯1 , . . . ) as C = AB . We use this compact notation from now on χ ( ) αβ ( , ) = i Tr [ Λ ( ) ( ) G Λ ( ) G ] (S28)Λ ( ) = Λ ( ) ( ) + V G Λ ( ) G (S29)Note that Tr [ . . . ] in the above formula stands for the sum over un-contracted spinor index. In the compact notation1 and 2 are space-time symbols.2 B. Conserving third order response theory in Dirac systems
Third-order response function is given by χ ( ) ( , , , ) ≡ δ J ( A ) δA ( ) δA ( ) δA ( ) (cid:187)(cid:187)(cid:187)(cid:187)(cid:187)(cid:187) A → = − i
3! Tr [ Λ ( ) ( ) δ G ( , + ; A ) δA ( ) δA ( ) δA ( ) ] A → (S30)Using GG − = δ GδA ( ) δA ( ) δA ( ) = − δ GδA ( ) δA ( ) δG − δA ( ) G − δ GδA ( ) δA ( ) δG − δA ( ) G + G δG − δA ( ) G δ G − δA ( ) δA ( ) G − δ GδA ( ) δA ( ) δG − δA ( ) G + G δG − δA ( ) G δ G − δA ( ) δA ( ) G + G δG − δA ( ) G δ G − δA ( ) δA ( ) G − G δ G − δA ( ) δA ( ) δA ( ) G . (S31)Using Eq. (S27), we have δ GδA ( ) δA ( ) δA ( ) = − δ GδA ( ) δA ( ) Λ ( ) G − δ GδA ( ) δA ( ) Λ ( ) G + G Λ ( ) G δ G − δA ( ) δA ( ) G − δ GδA ( ) δA ( ) Λ ( ) G + G Λ ( ) G δ G − δA ( ) δA ( ) G + G Λ ( ) G δ G − δA ( ) δA ( ) G − G δ G − δA ( ) δA ( ) δA ( ) G . (S32)Once more we use GG − = δ GδA ( ) δA ( ) = − δGδA ( ) δG − δA ( ) G − δGδA ( ) δG − δA ( ) G − G δ G − δA ( ) δA ( ) G . (S33)The two-photon vertex function is defined as followsΛ ( , ) = δ G − δA ( ) δA ( ) . (S34)Using Eq. (S21), Eq. (S27) and Eq. (S34), we find δ GδA ( ) δA ( ) = − G Λ ( , ) G + ∑ P ( ) G Λ ( ) G Λ ( ) G (S35)where P (
2; 3 ) stands for the permutation between 2 ↔
3. The three-photon vertex function readsΛ ( , , ) = δ G − δA ( ) δA ( ) δA ( ) . (S36)By plugging Eq. (S34) , Eq. (S35), and Eq. (S36) in Eq. (S32), we obtain δ GδA ( ) δA ( ) δA ( ) = − G Λ ( , , ) G − ∑ P ( ) G Λ ( ) G Λ ( ) G Λ ( ) G + ∑ P ( ) G Λ ( , ) G Λ ( ) G + ∑ P ( ) G Λ ( ) G Λ ( , ) G . (S37)where P (
2; 3; 4 ) stands for all six permutations among 2,3, and 4 space-time coordinates. Therefore, the third-orderresponse function reads χ ( ) ( , , , ) = i
3! Tr [ Λ ( ) ( ) G Λ ( , , ) G ] + i ∑ P ( ) Tr [ Λ ( ) ( ) G Λ ( ) G Λ ( ) G Λ ( ) G ] − i ∑ P ( ) Tr [ Λ ( ) ( ) G Λ ( , ) G Λ ( ) G ] − i ∑ P ( ) Tr [ Λ ( ) ( ) G Λ ( ) G Λ ( , ) G ] . (S38)3
1. Interaction induced two-photon vertex
The two-photon vertex function is defined as followsΛ ( , ) = δ G − δA ( ) δA ( ) = δ G − δA ( ) δA ( ) − δ Σ δA ( ) δA ( ) . (S39)The second derivative of the bare Green’s function vanishes in Dirac systems which implies the absence of bare two-photon vertex. However, an interaction induced two-photon vertex function is obtained owing to the field-dependentself-energy. Since the self-energy depends on the external field only through the dependence on the Green’s function,we have Λ ( , ) = − δ Σ δA ( ) δA ( ) = − V δ GδA ( ) δA ( ) . (S40)Using Eq. (S35), we obtain Λ ( , ) = Λ ( ) ( , ) + V G Λ ( , ) G (S41)where Λ ( ) ( , ) = − ∑ P ( ) V G Λ ( ) G Λ ( ) G . (S42)
2. Interaction induced three-photon vertex
The three-photon vertex function is defined as followsΛ ( , , ) = δ G − δA ( ) δA ( ) δA ( ) = δ G − δA ( ) δA ( ) δA ( ) − δ Σ δA ( ) δA ( ) δA ( ) . (S43)The third derivative of the bare Green’s function vanishes in Dirac systems which implies the absence of bare three-photon vertex. However, an interaction induced three-photon vertex function is obtained owing to the field-dependentself-energy. Since the self-energy depends on the external field only through the dependence on the Green’s function,we have Λ ( , , ) = − δ Σ δA ( ) δA ( ) δA ( ) = − V δ GδA ( ) δA ( ) δA ( ) . (S44)Using Eq. (S37), we obtain Λ ( , , ) = Λ ( ) ( , , ) + V G Λ ( , , ) G (S45)whereΛ ( ) ( , , ) = ∑ P ( ) V [ G Λ ( ) G Λ ( ) G Λ ( ) G − G Λ ( , ) G Λ ( ) G − G Λ ( ) G Λ ( , ) G ] . (S46)To summarise, the formal derivation given in the current Section guides us in constructing a conserving diagram-matic theory for the third-order response function in Dirac systems. Lengthy mathematical relations for differ-ent contributions to the nonlinear response function, i.e. Eq. (S38), and the multi-photon vertex functions, i.e.Eqs. (S29),(S41),(S42),(S45), and (S46), are graphically illustrated in Feynman diagrams depicted in Fig. 1 of themain text. Quantitative evaluation of these diagrams are explicitly discussed with great details in the next Section.4 S3. ANALYTICAL EXPRESSIONS FOR NONLINEAR CONDUCTIVITY
We present here the analytical expression of the third-order optical response function which in the Matsubara spacecan be formally written as: χ ( ) THG ( m , m , m ) = β ∑ n P ( n, n + m , n + m + m , n + m + m + m ) , (S47)where n stands for a short notation n = iω n ( ω n being a fermionic frequency), and m = iω m ( ω m being a bosonicfrequency). For the third-harmonic generation we have m = m = m = m and χ ( ) THG ( m ) = β ∑ n P ( n, n + m, n + m, n + m ) . (S48)After a straightforward algebra we perform the Matsubara summation and analytic continuation as iω m → ̵ hω + i + ,we obtain (see Section S6) χ ( ) THG ( ω ) = ∫ +∞−∞ d(cid:15) πi { n F ( (cid:15) ) P RRRR − n F ( (cid:15) + ̵ hω ) P AAAA + ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) )) P ARRR + ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) + ̵ hω )) P AARR + ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) + ̵ hω )) P AAAR } . (S49)Note that P RRRR ( (cid:15) , (cid:15) , (cid:15) , (cid:15) ) with (cid:15) j = (cid:15) + jω + iη j means that all frequency arguments are in the retarded channel(i.e. ( η j → + ) while P ARRR implies that the first argument is in the advanced channel, η → − + , but the other areretarded, η j ≠ → + . Therefore, we have P AAAA = ( P RRRR ) ∗ . By knowing the response function, the third-harmonicoptical conductivity reads σ ( ) THG ( ω ) = i χ ( ) THG ( ω ) ω . (S50) A. Different diagram contributions
Using Baym-Kadanoff analysis summarised in the previous section, we construct diagrams for the third-orderresponse function in Dirac materials as illustrated in Fig. 1 of the main text. As depicted in Fig. 1a of the main text,the P -function contains four main different contributions, P = P + P + P + P , associated respectively with square( P ), triangles ( P , P ) and bubble ( P ) diagrams. More explicitly we can write: P ( z , z , z , z ) = ( e v N f π ̵ h ) Q ( z , z ) Q ( z , z ) Q ( z , z ) Ω ( z , z , z , z ) (S51)The sum over spin and valley indices just leads to an overall degeneracy factor N f = N s N v where N s = N v = ( z , z , z , z ) is the bare square diagram in the absence of vertex renormalization,Ω ( z , z , z , z ) = γ imp U ∑ k Tr [ ˆ σ y ˆ G ( k , z ) ˆ σ y ˆ G ( k , z ) ˆ σ y ˆ G ( k , z ) ˆ σ y ˆ G ( k , z )] , (S52)and the one-photon vertex ˆΛ ( z i , z j ) = ( − ev ˆ σ y ) Λ ( ) Q ( z i , z j ) where Λ ( ) = Q ( z i , z j ) is the one-photon Bethe-Salpeter renormalization factor which is given in the following subsection. Insimilar way we can write the contributions of the two triangles diagrams, namely P ( z , z , z , z ) = ( e v N f π ̵ h ) Q ( z , z ) Q ( z , z ) Λ ( ) ( z , z , z ) Ω ( z , z , z ) , (S53)where Ω ( z , z , z ) = − γ imp U ∑ k Tr [ ˆ σ y ˆ G ( k , z ) ˆ G ( k , z ) ˆ σ y ˆ G ( k , z )] , (S54)5Note that ˆΛ ( z , z , z ) = ( − ev ˆ σ y ) Λ ( ) ( z , z , z ) Q ( z , z ) is the two-photon vertex function where Λ ( ) ( z , z , z ) is the unrenormalized two-photon vertex function and Q ( z , z ) is the two-photon Bethe-Salpeter renormalizationfactor (see subsections below). We have also the further triangle diagram: P ( z , z , z , z ) = ( e v N f π ̵ h ) Q ( z , z ) Q ( z , z ) Λ ( ) ( z , z , z ) Ω ( z , z , z ) . (S55)where Ω ( z , z , z ) = − γ imp U ∑ k Tr [ ˆ σ y ˆ G ( k , z ) ˆ σ y ˆ G ( k , z ) ˆ G ( k , z )] , (S56)Finally we make the bubble term explicit: P ( z , z , z , z ) = ( e v N f π ̵ h ) Λ ( ) ( z , z , z , z ) Q ( z , z ) X ( z , z ) , (S57)where ˆΛ ( z , z , z , z ) = ( − ev ˆ σ y ) Λ ( ) ( z , z , z , z ) Q ( z , z ) is the three-photon vertex function withΛ ( ) ( z , z , z , z ) being the unrenormalized three-photon vertex function and Q ( z , z ) is the three-photon Bethe-Salpeter renormalizationfactor (see subsections below). The explicit expressions of Λ ( ) ( z , z , z ) , Λ ( ) ( z , z , z , z ) ,and Q n = , , ( z i , z j ) , Ω ( z , z , z , z ) , Ω ( z , z , z ) , and Ω ( z , z , z ) are provided in the next subsections. B. Renormalization of the one-photon vertex
The one-photon vertex renormalization is depicted diagrammatically in Fig. 1b of the main text and it readsˆΛ ( p , p + q ; n, n + m ) − ˆΛ ( ) ( p , p + q ; n, n + m ) = γ imp ∑ k ˆ G ( k , n ) ˆΛ ( k , k + q ; n, n + m ) ˆ G ( k + q , n + m ) , (S58)where m ≡ iq m and n ≡ ik n = ip n stand for the bosonic and fermionic Matsubara frequencies, respectively. Notethat in the integrand we have shifted the dummy momentum k as k + p → k and therefore we can see that vertexcorrection does not depends on the fermion momentum p . For the optical (or dipole) approximation we have q = .Therefore, the Bethe-Salpeter relation for the one-photon vertex function readsˆΛ ( n, n + m ) = ˆΛ ( ) + γ imp ∑ k ˆ G ( k, n ) ˆΛ ( n, n + m ) ˆ G ( k, n + m ) . (S59)Note that we have ˆΛ ( ) = δ ˆ G − / δA ∣ A → = − ev ˆ σ y where ˆ G stands for the non-interacting Green’s function. Weassume the following ansatz for the vertex functionˆΛ ( n, n + m ) = a ˆ I + bσ x + ( c + v ) σ y + dσ z . (S60)Using the fact that the integral of odd-function of k is zero, we obtain a = b = d = = ( − ev ˆ σ y ) Λ whereΛ ( n, n + m ) = Λ ( ) − U X ( n, n + m ) . (S61)Note that Λ ( ) = X ( n, n + m ) = γ imp U ∑ k Tr [ ˆ σ y ˆ G ( k , n ) ˆ σ y ˆ G ( k , n + m )] . (S62)Using dimensional regularization, we find the following formula for the X function: X ( n, n + m ) = S ( n ) S ( n + m ) S ( n ) − S ( n + m ) ln [ S ( n + m ) S ( n ) ] . (S63)The above equation can be straightforwardly generalized in the generic complex space by replacing n → iω n → z , n + m → iω n + iω m → z . Eq. (S61) defines the one-photon vertex renormalization factor: Q ( z , z ) = − U X ( z , z ) . (S64)Obviously X ( z , z ) = X ( z , z ) .6 C. Renormalization of the two-photon vertex
Similar to the one-photon vertex case, it can be shown that the two-photon vertex function is independent of thefermionic momentum, p . Moreover, for the optical limit we can neglect the photon momentum q . The self-consistentBethe-Salpeter relation for the two-photon vertex function is depicted in Fig. 1c of the main text and it readsˆΛ ( n, n + m, n + m ) = ˆΛ ( ) ( n, n + m, n + m ) + γ imp ∑ k ˆ G ( k , n ) ˆΛ ( n, n + m, n + m ) ˆ G ( k , n + m ) . (S65)In the non-interacting Dirac system the “bare” two-photon vertex function is zero, ˆΛ ( ) ∝ δ ˆ G − / δA ∣ A → =
0, dueto the linear momentum dependence of the Hamiltonian. However, due to interaction the unrenormalized two-photonvertex ˆΛ ( ) is finite given by the following relation (see Fig. 1e of the main text)ˆΛ ( ) ( n, n + m, n + m ) = − γ imp ∑ k ˆ G ( k , n ) ˆΛ y ( n, n + m ) ˆ G ( k , n + m ) ˆΛ y ( n + m, n + m ) ˆ G ( k , n + m ) . (S66)From now on we adopt the short-hand notation z j = n + jm with j = , , ,
3. We find ˆΛ ( ) = ( − ev ˆ σ y ) Λ ( ) withˆ σ y = ˆ I and Λ ( ) ( z , z , z ) = Q ( z , z ) Q ( z , z ) U Z ( z , z , z ) , (S67)in which Q ( z i , z j ) is the one-photon renormalization factor defined in the previous subsection, and where Z ( z , z , z ) = − γ imp U ∑ k Tr [ ˆ G ( k , z ) ˆ σ y ˆ G ( k , z ) ˆ σ y ˆ G ( k , z )] . (S68)By performing the momentum integration using the dimensional regularization, we obtain Z ( z , z , z ) = S ( z ) S ( z ) − S ( z ) { S ( z ) S ( z ) − S ( z ) ln [ S ( z ) S ( z ) ] + S ( z ) S ( z ) − S ( z ) ln [ S ( z ) S ( z ) ] } . (S69)In a compact form we can write Z ( z , z , z ) = X ( z , z ) − X ( z , z ) S ( z ) − S ( z ) . (S70)By solving the the self-consistent Bethe-Salpeter relation for the two-photon vertex given in Eq. (S65), we obtainˆΛ = ( − ev ˆ σ y ) Λ with Λ ( z , z , z ) = Q ( z , z ) Λ ( ) ( z , z , z ) , (S71)in which Q ( z , z ) is the two-photon Bethe-Salpeter renormalization factor Q ( z , z ) = − U X ( z , z ) , (S72)and where X ( z, z ′ ) = γ imp U ∑ k Tr [ ˆ G ( k , z ) ˆ G ( k , z ′ )] . (S73)We explicitly obtain X ( z, z ′ ) = S ( z ) − S ( z ′ ) { S ( z ) ln [ − W S ( z ) ] − S ( z ′ ) ln [ − W S ( z ′ ) ] } . (S74)Note that X ( z, z ′ ) = X ( z ′ , z ) . Using the self-energy relation Eq. (S10) we have S ( z ) ln [ − W S ( z ) ] = − Σ ( z ) U . (S75)7Therefore, the two-photon renormalizationfactor reads (Eq. (6) of the main text) Q ( z, z ′ ) = − U X ( z, z ′ ) = S ( z ) − S ( z ′ ) z − z ′ . (S76)It is useful to evaluate of Q RR2 ( (cid:15), (cid:15) ) which is given in Eq. (7) of the main text. Using the above relation in theretarded-retarded (RR) channel and employing the self-energy relation Eq. (S10), we obtain Q RR2 ( (cid:15), (cid:15) ) = S ′ ( (cid:15) ) = dS ( (cid:15) ) d(cid:15) = − d Σ ( (cid:15) ) d(cid:15) = + U S ′ ( (cid:15) ) ln [ − W S ( (cid:15) ) ] − U S ( (cid:15) ) d ln [ S ( (cid:15) ) ] d(cid:15) = + U S ′ ( (cid:15) ) ln [ − W S ( (cid:15) ) ] − U S ′ ( (cid:15) ) . (S77)Therefore, we have [ + U − U ln [ − W S ( (cid:15) ) ]] S ′ ( (cid:15) ) = . (S78)Consequently, we arrive at S ′ ( (cid:15) ) = + U − U ln [ − W S ( (cid:15) ) ] . (S79)Using Eq. (S10), we obtain 1 − U ln [ − W S ( (cid:15) ) ] = + Σ ( (cid:15) ) S ( (cid:15) ) = µ + (cid:15)S ( (cid:15) ) . (S80)Therefore, we find the following relation which is given in Eq. (7) of the main text. Q RR2 ( (cid:15), (cid:15) ) = S ′ ( (cid:15) ) = S ( (cid:15) ) U S ( (cid:15) ) + µ + (cid:15) . (S81) D. Renormalization of the three-photon vertex
Similar to the case of two-photon case, the impurity scattering induces a finite three-photon vertex as defined inFig. 1f of the main text. Accordingly we find ˆΛ ( ) = ( − ev ˆ σ y ) Λ ( ) withΛ ( ) ( z , z , z , z ) = M ( z , z , z , z ) + M ( z , z , z , z ) + M ( z , z , z , z ) (S82)where M ( z , z , z , z ) = U Ω ( z , z , z , z ) Q ( z , z ) Q ( z , z ) Q ( z , z ) , (S83) M ( z , z , z , z ) = U Ω ( z , z , z ) Λ ( ) ( z , z , z ) Q ( z , z ) Q ( z , z ) , (S84) M ( z , z , z , z ) = U Ω ( z , z , z ) Λ ( ) ( z , z , z ) Q ( z , z ) Q ( z , z ) . (S85)Here Q ( z i , z j ) , Q ( z i , z j ) are the Bethe-Salpeter one- and two-photon renormalization functions, respectively. Theexplicit expression for Ω function is given byΩ ( z , z , z , z ) = u ( z , z , z , z ) ln [ S ( z ) S ( z ) ] + u ( z , z , z , z ) ln [ S ( z ) S ( z ) ] + u ( z , z , z , z ) ln [ S ( z ) S ( z ) ] (S86)8where u ( z , z , z , z ) = S ( z ) S ( z ) S ( z ) S ( z ) + [ S ( z ) S ( z ) + S ( z ) S ( z )] S ( z ) [ S ( z ) − S ( z ) ][ S ( z ) − S ( z ) ][ S ( z ) − S ( z ) ] , (S87) u ( z , z , z , z ) = S ( z ) S ( z ) S ( z ) S ( z ) + [ S ( z ) S ( z ) + S ( z ) S ( z )] S ( z ) [ S ( z ) − S ( z ) ][ S ( z ) − S ( z ) ][ S ( z ) − S ( z ) ] , (S88) u ( z , z , z , z ) = S ( z ) S ( z ) S ( z ) S ( z ) + [ S ( z ) S ( z ) + S ( z ) S ( z )] S ( z ) [ S ( z ) − S ( z ) ][ S ( z ) − S ( z ) ][ S ( z ) − S ( z ) ] . (S89)Similarly, one can obtain Ω ( z , z , z ) = Z ( z , z , z ) , (S90)Ω ( z , z , z ) = Z ( z , z , z ) . (S91)Finally, the Bethe-Salpeter renormalization of the three-photon vertex function gives:Λ ( z , z , z , z ) = Q ( z , z ) Λ ( ) ( z , z , z , z ) . (S92)where Q ( z , z ) = − U X ( z , z ) (S93)Note that X ( z , z ) = X ( z , z ) . S4. LINEAR CONDUCTIVITY
Linear response function is obtained after performing a Matsubara summation as follows χ ( ) ( m ) = β ∑ n P ( n, n + m ) . (S94)Note that β = / k B T where k B is the Boltzmann constant and T stands for the electronic temperature. The P -function is analytical in the complex plain except two branch cuts at (cid:15) and (cid:15) − m where (cid:15) span over whole real axes.After performing the summation and an analytical continuation as m → ̵ hω + i + , we find [49] χ ( ) ( ω ) = ∫ ∞−∞ d(cid:15) πi {[ n F ( (cid:15) ) − n F ( (cid:15) + ̵ hω )] P AR ( (cid:15), (cid:15) + ̵ hω ) − n F ( (cid:15) ) P RR ( (cid:15), (cid:15) + ̵ hω ) + n F ( (cid:15) + ̵ hω ) P AA ( (cid:15), (cid:15) + ̵ hω )} (S95)where n F ( x ) = /( e βx + ) is the Fermi-Dirac distribution function. Note that “R” and “A” superscripts stand for theretarded and advanced, respectively. Accordingly, we have P RR ( (cid:15), (cid:15) + ̵ hω ) = [ P AA ( (cid:15), (cid:15) + ̵ hω )] ∗ = P ( (cid:15) + i + , (cid:15) + ̵ hω + i + ) and P AR ( (cid:15), (cid:15) + ̵ hω ) = P ( (cid:15) − i + , (cid:15) + ̵ hω + i + ) . The linear optical conductivity reads σ ( ) ( ω ) = i χ ( ) ( ω ) ω . (S96)For the 2D Dirac model, the P -function reads P ( n, n + m ) = N f e π ̵ h X ( n, n + m ) − U X ( n, n + m ) . (S97)Note that the sum over spin and valley index just leads to an overall degeneracy factor N f = N s N v where N s = N v =
2. The linear dc -conductivity follows σ ( ) dc = σ π { X AR1 ( , ) − U X
AR1 ( , ) − X RR1 ( , ) − U X
RR1 ( , ) } (S98)9where σ = e / ̵ h . The retarded self-energy can be decomposed into its real and imaginary parts Σ ( (cid:15) ) = ∆ ( (cid:15) ) − i Γ ( (cid:15) ) with ∆ ( (cid:15) ) and Γ ( (cid:15) ) > X RR1 ( , ) = − X AR1 ( , ) = w ( x ) with x = µ / Γ ( µ ) in which the renormalized chemical potential is given by µ = µ − ∆ ( µ ) , and w ( x ) reads w ( x ) = + x x arctan [ − x + x , x + x ] . (S99)Therefore, we find σ ( ) dc = σ f ( µ Γ ( µ ) ; U ) (S100)where f ( x ; U ) = π { + U + w ( x ) − U w ( x ) } . (S101)The functional dependence of f ( x ; U ) on x and U is illustrated in Fig. S2. In the constant-Γ model, Σ = − i Γ, we U U U U U f ( x ; U ) FIG. S2: Universal f ( x ; U ) function versus x for several values of U . have f ( x ) = [ + w ( x )] π . (S102)The asymptotic form of f ( x ) for small and large x follows f ( x ) ≈ π ( + x ) , x ≪ , (S103) f ( x ) ≈ π ( x + x ) , x ≫ . (S104)0 S5. NONLINEAR DC CONDUCTIVITYA. Analytical derivation f ( x ) function in the constant- Γ model Nonlinear dc conductivity is given by σ ( ) dc = − lim ω → Im [ χ ( ) THG ( ω )] ω (S105)where the imaginary part of nonlinear (third-harmonic) response function followsIm [ χ ( ) THG ( ω )] = ∫ +∞−∞ d(cid:15) π { ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) )) Re [ P RRRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )] − ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) )) Re [ P ARRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )] − ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) + ̵ hω )) Re [ P AARR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )] − ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) + ̵ hω )) Re [ P AAAR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )]} . (S106)For small frequency we can expand the Fermi function as follows n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) ) ≈ ̵ hω∂ (cid:15) n F ( (cid:15) ) + ( ̵ hω ) ∂ (cid:15) n F ( (cid:15) ) + ( ̵ hω ) ∂ (cid:15) n F ( (cid:15) ) + . . . = − ( ̵ hω ) δ ( (cid:15) ) − ( ̵ hω ) ∂ (cid:15) δ ( (cid:15) ) − ( ̵ hω ) ∂ (cid:15) δ ( (cid:15) ) + . . . . (S107)Notice that at zero temperature we have n F ( (cid:15) ) = Θ ( − (cid:15) ) that is the Heaviside step function and its first derivative ∂ (cid:15) n F ( (cid:15) ) = − δ ( (cid:15) ) which stands for the Dirac delta function. After integration over (cid:15) , we findIm [ χ ( ) THG ( ω )] ≈ − ( ̵ hω ) h ( ω ) + ( ̵ hω ) h ( ω ) − ( ̵ hω ) h ( ω ) (S108)where we define h ( ω ) = [ P RRRR ( , ̵ hω, ̵ hω, ̵ hω )] − Re [ P ARRR ( , ̵ hω, ̵ hω, ̵ hω )] − Re [ P AARR ( , ̵ hω, ̵ hω, ̵ hω )] − Re [ P AAAR ( , ̵ hω, ̵ hω, ̵ hω )] , (S109) h ( ω ) = lim (cid:15) → ∂∂(cid:15) { [ P RRRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )] − Re [ P ARRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )] − [ P AARR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )] − [ P AAAR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )]} , (S110) h ( ω ) = lim (cid:15) → ∂ ∂(cid:15) { [ P RRRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )] − Re [ P ARRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )] − [ P AARR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )] − [ P AAAR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )]} . (S111)In the absence of vertex correction, we have P ( z , z , z , z ) = ( e v N f π ̵ h ) Ω ( z , z , z , z ) . (S112)We assume Σ ( (cid:15) ) = ∓ i Γ where Γ > − / + stands for the retarded (R) andadvanced (A) channel, respectively. For the case of ̵ hω ≪ µ, Γ, it is legitimate to expand the integrand for small ω .In the constant-Γ model, the contribution from h ( ω ) exactly cancels that of h ( ω ) . Eventually, the nonlinear dcconductivity in the constant-Γ model reads σ ( ) dc = ( e v N f ̵ h π ̵ h t ) t Γ f ( µ Γ ) . (S113)1Considering ̵ hv = √ t a / N f =
4, we have ( e v N f ̵ h ( π ) ̵ h t ) = σ E (S114)with σ = e / ̵ h and E = π √ t ea (S115)where with t ≈ a ≈ . E ≈ / nm. The universal f ( x ) function reads f ( x ) = [ x − x ] + [ + x ] + [ x − + x − x ] w ( x ) . (S116)There is a sign change in f ( x ) at x ≈ . f ( x ) ≈ ( − x ) , x ≪ , (S117) f ( x ) ≈ − π ( x − x ) , x ≫ . (S118) B. Numerical evaluation of f ( x, U ) function in the full quantum theory The third-harmonic generation (THG) conductivity is given by σ ( ) THG ( ω ) = i χ ( ) THG ( ω ) ω = ( e v N f π ̵ h ) ω ∫ ∞−∞ d(cid:15) π Re [ K ( (cid:15), ω )] = ( e v N f ̵ h ( π ) ̵ h t ) [ t ( ̵ hω ) ∫ ∞−∞ d(cid:15) Re [ K ( (cid:15), ω )]] (S119)where K ( (cid:15), ω ) = ( e v N f π ̵ h ) − { n F ( (cid:15) ) P RRRR − n F ( (cid:15) + ̵ hω ) P AAAA + ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) )) P ARRR + ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) + ̵ hω )) P AARR + ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) + ̵ hω )) P AAAR } . (S120)Therefore, we obtain σ ( ) THG ( ω ) = σ E t ( ̵ hω ) ∫ ∞−∞ d(cid:15) Re [ K ( (cid:15), ω )] . (S121)We evaluate the third-order dc conductivity as dc limit of the third-harmonic conductivity, σ ( ) dc = lim ω → σ ( ) THG ( ω ) = σ E [ t Γ ( µ ) ] f ( µ Γ ( µ ) ; U ) (S122)where the universal f function can be evaluated numerically by using the following relation f ( µ Γ ( µ ) ; U ) = lim ω → Γ ( µ ) ( ̵ hω ) ∫ ∞−∞ d(cid:15) Re [ K ( (cid:15), ω )] . (S123)2 z = ✏ z = ✏ m z = ✏ mz = ✏ m FIG. S3: The Matsubara summation is performed by utilizing an integration on a contour enclosing whole complex plane exceptfour branch cuts which are shown by dashed lines on the contour with radius R → ∞ . S6. ANALYTICAL CONTINUATION FOR THE THIRD ORDER RESPONSE FUNCTION
The summation of the fermionic Matsubara frequency n is performed by the contour integration technique. B ( m ) = β ∑ n P ( n, n + m, n + m, n + m ) . (S124)The P -function contains four brach cut in complex plane. The Matsubara summation is performed on a contour withfour cuts at (cid:15) , (cid:15) − m , (cid:15) − m , and (cid:15) − m , see Fig. S3. Note that (cid:15) runs over the entire real axes. Therefore, we write B ( m ) = ∮ dz πi n F ( z ) P ( z, z + m, z + m, z + m ) = ∫ +∞−∞ d(cid:15) πi { n F ( (cid:15) + iη ) P ( (cid:15) + iη, (cid:15) + iη + m, (cid:15) + iη + m, (cid:15) + iη + m ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi { n F ( (cid:15) − m + iη ) P ( (cid:15) − m + iη, (cid:15) + iη, (cid:15) + iη + m, (cid:15) + iη + m ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi { n F ( (cid:15) − m + iη ) P ( (cid:15) − m + iη, (cid:15) − m + iη, (cid:15) + iη, (cid:15) + iη + m ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi { n F ( (cid:15) − m + iη ) P ( (cid:15) − m + iη, (cid:15) − m + iη, (cid:15) − m + iη, (cid:15) + iη ) − ( η → − η )} (S125)where n F ( z ) = /( + e βz ) is the Fermi-Dirac distribution function. Note that m stands for external bosonic Matsobara(imaginary) frequency. Since η ≪ ∣ m ∣ , we have m + iη → m which implies B ( m ) = ∫ +∞−∞ d(cid:15) πi { n F ( (cid:15) + iη ) P ( (cid:15) + iη, (cid:15) + m, (cid:15) + m, (cid:15) + m ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi { n F ( (cid:15) − m ) P ( (cid:15) − m, (cid:15) + iη, (cid:15) + m, (cid:15) + m ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi { n F ( (cid:15) − m ) P ( (cid:15) − m, (cid:15) − m, (cid:15) + iη, (cid:15) + m ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi { n F ( (cid:15) − m ) P ( (cid:15) − m, (cid:15) − m, (cid:15) − m, (cid:15) + iη ) − ( η → − η )} . (S126)3Note that for bosonic frequency m we have n F ( (cid:15) − m ) = n F ( (cid:15) ) and therefore we find B ( m ) = ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ) { P ( (cid:15) + iη, (cid:15) + m, (cid:15) + m, (cid:15) + m ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ) { P ( (cid:15) − m, (cid:15) + iη, (cid:15) + m, (cid:15) + m ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ) { P ( (cid:15) − m, (cid:15) − m, (cid:15) + iη, (cid:15) + m ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ) { P ( (cid:15) − m, (cid:15) − m, (cid:15) − m, (cid:15) + iη ) − ( η → − η )} . (S127)We do analytical continuation as m → ̵ hω + iη : B ( ω ) = ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ) { P ( (cid:15) + iη, (cid:15) + ̵ hω + iη, (cid:15) + ̵ hω + i η, (cid:15) + ̵ hω + i η ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ) { P ( (cid:15) − ̵ hω − iη, (cid:15) + iη, (cid:15) + ̵ hω + iη, (cid:15) + ̵ hω + i η ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ) { P ( (cid:15) − ̵ hω − ̵ hω − i η, (cid:15) − ̵ hω − iη, (cid:15) + iη, (cid:15) + ̵ hω + iη ) − ( η → − η )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ) { P ( (cid:15) − ̵ hω − i η, (cid:15) − ̵ hω − i η, (cid:15) − ̵ hω − iη, (cid:15) + iη ) − ( η → − η )} . (S128)Since η → + , we have B ( ω ) = ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ){ P RRRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω ) − P ARRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ){ P ARRR ( (cid:15) − ̵ hω, (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω ) − P AARR ( (cid:15) − ̵ hω, (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ){ P AARR ( (cid:15) − ̵ hω, (cid:15) − ̵ hω, (cid:15), (cid:15) + ̵ hω ) − P AAAR ( (cid:15) − ̵ hω, (cid:15) − ̵ hω, (cid:15), (cid:15) + ̵ hω )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ){ P AAAR ( (cid:15) − ̵ hω, (cid:15) − ̵ hω, (cid:15) − ̵ hω, (cid:15) ) − P AAAA ( (cid:15) − ̵ hω, (cid:15) − ̵ hω, (cid:15) − ̵ hω, (cid:15) )} . (S129)Note that “R” and “A” superscript stand for the retarded and advanced, respectively. We shift (cid:15) in such a way thatall P -function arguments is ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω ) : B ( ω ) = ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) ){ P RRRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω ) − P ARRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) + ̵ hω ){ P ARRR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω ) − P AARR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) + ̵ hω + ̵ hω ){ P AARR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω ) − P AAAR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )} + ∫ +∞−∞ d(cid:15) πi n F ( (cid:15) + ̵ hω ){ P AAAR ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω ) − P AAAA ( (cid:15), (cid:15) + ̵ hω, (cid:15) + ̵ hω, (cid:15) + ̵ hω )} . (S130)4Eventually, we obtain B ( ω ) = ∫ +∞−∞ d(cid:15) πi { [ n F ( (cid:15) ) P RRRR − n F ( (cid:15) + ̵ hω ) P AAAA ] + ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) )) P ARRR + ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) + ̵ hω )) P AARR + ( n F ( (cid:15) + ̵ hω ) − n F ( (cid:15) + ̵ hω )) P AAAR } . (S131) S7. NUMERICAL EVALUATION OF THE NONLINEAR OPTICAL AND DC CONDUCTIVITIES
In Fig. S4, we plot the frequency dependence of the third-harmonic generation (THG) response function σ ( ) THG ( ω ) for µ =
17 THz and few representative values of U . For comparison we show also the THG optical response σ ( ) THG ( ω ) for non-interacting electrons which is null for ̵ hω ≤ µ /
3. Note the change of sign of σ ( ) dc varying the scatteringstrength. ⇥ ⇥ U = = = = - - - - ω / μ R e [ σ T H G ( ) ( ω ) ]/ σ ( ) FIG. S4: Real part of third-harmonic optical conductivity versus frequency in comparison with the non-interacting result. Notethat the chemical potential is set µ =
17 THz, and σ ( ) = σ / E . In Fig. S5, we illustrate the universal scaling of f versus U in the quantum regime for different values of x = µ / Γ ( µ ) <
1. As seen the slop of the curves in the log-log scale plot does not strongly depends on the value of x whichsupport the validity of Eq. (5) given in the main text.In Fig. S6, we show the phase diagram for the constant-Γ model. As it is seen this phase diagram is completelydifferent from that of the full quantum theory which is given in Fig. 3c of the main. text. We can see only onesign-change in the constant-Γ model in contrast to that of full quantum theory which gives two sing-changes. Unlikethe full quantum theory, the constant-Γ model predicts a positive nonlinear correction in the quantum regime.5 SquareTrianglesTriangles, X = - PhotonFull - - - - - μ [ THz ] σ d c ( ) E / σ d c ( ) OOVVXX O x = V x = X x = U | f t r i ang l e s ( x ; U ) | OOVVXX O x = V x = X x = U | f t r i ang l e s ( x ; U ) | (c) OOVVXX O x = V x = X x = U | f t r i ang l e s ( x ; U ) | (a) x=0.1 x=0.7 (b) x=0.4 - - U | f ( x ; U ) | - - U | f ( x ; U ) | - - U | f ( x ; U ) | FIG. S5: Log-log scale plot for the absolute value of the universal f ( µ, U ) function versus U at x = . , . . NegativePositive
FIG. S6: Colormap plot for g = ∣ σ ( ) dc E / σ ( ) dc ∣ factor with E = / nm versus chemical potential µ and scattering rate Γ inthe constant-Γ model, Σ = − i Γ. The sign of σ ( ) dc is written on the plot where the sign-switch border is highlighted by a dashedred line. Green and blue dotted lines stand for the contour lines with g = g = ..