Many-Body Effects in Third Harmonic Generation of Graphene

The low-energy (intraband) range of the third harmonic generation of graphene in the terahertz regime is governed by the damping terms induced by the interactions. A controlled many-body description of the scattering processes is thus a compelling and desirable requirement. In this paper, using a Kadanoff-Baym approach, we systematically investigate the impact of many-body interaction on the third-harmonic generation (THG) of graphene, taking elastic impurity scattering as a benchmark example. We predict the onset in the mixed inter-intraband regime of novel incoherent features driven by the interaction at four- and five-photon transition frequencies in the third-harmonic optical conductivity with a spectral weight proportional to the scattering rate.We show also that, in spite of the complex many-body physics, the purely intraband term governing the limit $\omega \to 0$ resembles the constraints of the phenomenological model. We ascribe this agreement to the fulfilling of the conservation laws enforced by the conserving approach. The overlap with novel incoherent features and the impact of many-body driven multi-photon vertex couplings limit however severely the validity of phenomenological description.

The low-energy (intraband) range of the third harmonic generation of graphene in the terahertz regime is governed by the damping terms induced by the interactions. A controlled many-body description of the scattering processes is thus a compelling and desirable requirement. In this paper, using a Kadanoff-Baym approach, we systematically investigate the impact of many-body interaction on the third-harmonic generation (THG) of graphene, taking elastic impurity scattering as a benchmark example. We predict the onset in the mixed inter-intraband regime of novel incoherent features driven by the interaction at four-and five-photon transition frequencies in the third-harmonic optical conductivity with a spectral weight proportional to the scattering rate.We show also that, in spite of the complex many-body physics, the purely intraband term governing the limit ω → 0 resembles the constraints of the phenomenological model. We ascribe this agreement to the fulfilling of the conservation laws enforced by the conserving approach. The overlap with novel incoherent features and the impact of many-body driven multi-photon vertex couplings limit however severely the validity of phenomenological description.

I. INTRODUCTION
A simple dimensional analysis shows that the zerotemperature third-harmonic response of clean graphene scales as σ (3) ∼ 1/(ω 3 |µ|) [1][2][3][4][5] where µ is the chemical potential and ω is the frequency of the incident field. Comparing to the scaling of linear conductivity at low frequency, σ (1) ∼ |µ|/ω, we can thus expect a huge enhancement of the nonlinear optical response at low frequencies. Accordingly, there is surge of experimental interest in exploring nonlinear optics of graphene and other two-dimensional materials in the terahertz frequency range. The effective three-dimensional THG susceptibility of graphene with an effective thickness ∼ 0.1 nm has been measured as 10 −19 − 10 −16 m 2 /V 2 [6][7][8][9] in the near-infrared/visible frequency range ω ∼ 100 − 500THz, whereas in the terahertz range ( ω ∼ 1THz) a susceptibility as large as 10 −9 m 2 /V 2 is obtained [10,11]. In these works however the nonlinear response was not related to multi-photon absorption/emission but rather to the underlying nonlinear dependence of the low-frequency intraband linear conductivity on the electronic temperature ruled by the magnitude of the incident laser power [9,12,13]. The intrinsic role of multi-photon processes is thus completely missing in the interpretation of this experimental observation.
In common graphene, the terahertz frequency range, 1THz ∼ 4meV, is much below the inter-band transition, ω 2|µ|, and therefore the optical response is mainly governed by intra-band processes. Even in lowdoping regime, the THz-range inter-band transitions will be essentially Pauli-blocked owing to temperature effects, charge-puddle formation and many-body induced bandbroadening. The intra-band response is very sensitive to the scattering processes which result in the momentum * habib.rostami@su.se relaxation of quasiparticles. Therefore, the many-body interaction is expected to have an unavoidable impact on the intra-band optical response of the Dirac fermions in graphene. A compelling study of nonlinear responses is usually a delicate and cumbersome task. Therefore, the number of studies exploring many-body effects on the nonlinear optical response in a systematic way is quite limited. [5,[14][15][16][17][18]. A practical short-cut for including the impact of the scattering in the analysis is through a phenomenological relaxation rate Γ independent from energy and field [2,3]. This rough approximation may work qualitatively well in the high-doping and high-frequency (|µ|, | ω| Γ) regime. However, even in high quality graphene samples a scattering rate Γ not less than Γ ≈ 2 − 5 THz was estimated in a wide range µ ∼ 0 − 200 meV [19,20]. Therefore, in the terahertz range ( ω ∼ 1 THz) we are rather in the regime ω < Γ, questioning the validity of a constant-Γ model . Furthermore, in lowfrequency regime vertex corrections to the linear ones are expected to be extremely relevant due to the proximity to the multi-photon self-generation [18].
Aim of the present work is to provide a compelling many-body approach for the nonlinear response of graphene and other two-dimensional Dirac systems in the terahertz frequency range, where the role of scattering is fundamentally important. We derive a conserving quantum theory based on the Kadanoff-Baym framework [21] which contains both the self-energy and vertex correction effects on the same footing. Our theory is therefore consistent with conservation laws and the gauge invariance. To further preserve the gauge invariance and not violate the Ward's identity, we employ a dimensional regularisation [22,23] in the evaluation of momentum integrals in Dirac model of graphene for which a high-energy cut-off is unavoidable. For the sake of simplicity, we only focus on the impact of disorder scattering within the selfconsistent Born approximation (SCBA). However, our formal and technical development is valid also for scattering with phonons and other electrons in the system. Most of our qualitative predictions are thus expected to be valid also for other types of scattering processes.
Our main results can be summarized as: (i) we predict the onset of four-and five-photon transitions in the third harmonic response of graphene. As a result of such four-and five-photon transitions, the THG of graphene is strongly enhanced in the intra-band regime where ω is smaller than three-photon transition edge. These novel transitions are intrinsically driven by the interaction with a spectral weight that scales with the magnitude of the one-particle scattering rate, revealing the intrinsic incoherent character; (ii) a strong impact of vertex corrections is revealed owing the presence of many-bodyinduced two-and three-photon current vertices which are absent in non-interacting Dirac fermions in graphene. A crucial role in this regards is played in particular by the occurring of the two-photon vertex self-generation in the intra-band terahertz regime, close to the dc limit; (iii) in the extreme low-frequency limit, we find a good agreement of the pure intraband term in the presence of manybody effects with the phenomenological models, as resulting of enforcing a conserving approach. Departures from this modelling are however observed also at relatively small scattering due to the onset of incoherent four and five-photon transitions and the two-photon vertex selfgeneration. Our theoretical modeling can be generalized in a straightforward way to explain intra-band THG in other two-dimensional materials such as transition-metal dichalcogenides, homo-and hetero-bilayer systems, etc.
The paper is structured in five Sections. In Section II, we introduce the conserving Baym-Kadanoff derivation employed to evaluate the nonlinear current within the Dirac model, and we introduce the elastic impurityinduced self-energy. In Section III, we formally derive the many-body-induced multi-photon vertices based on selfconsistent Bethe-Salpeter equations within the Kadanoff-Baym method. In Section IV we provide all of the analytical relations for the third-order response function using a diagrammatic quantum theory for THG in graphene and for generic two-dimensional Dirac systems. In Section V, we present our numerical results for the real and imaginary parts of the third-harmonic conductivity in graphene and we discuss the onset of novel incoherenttransition peaks, the impact of vertex renormalization and the spectral features of pure intraband processes.. Eventually in Section VI we provide a summary and conclusion.

II. MODEL AND METHOD
We use the Dirac Hamiltonian of low-energy carriers in graphene asĤ where v ∼ 10 6 m/s is the Fermi velocity and where the Hamiltonian includes the dependence both on the pseudospin sublattice degree of freedom and on the valley.
More explicitly we writeσ = (τσ x ,σ y ), whereσ x ,σ y stand for the Pauli matrices in the sublattice basis and τ = ± accounts for the two inequivalent valleys in the Brillouin zone of graphene.
In dipole approximation we can model light-matter interaction by applying minimal coupling transformation k → k + eA(t) where A(t) stands for an external vector potential and the corresponding electric field is given by E(t) = −∂ t A(t). For sake of shortness that the speed of light is set c = 1. The inverse of the bare Green's function in the presence of the external vector potential reads thus: where we use the shorthand notation 1 ≡ (r 1 , t 1 ) for the space-time coordinate. In the presence of many-body scattering, it is useful to introduce an interacting Green's function:Ĝ where . . . stands for the thermodynamical average, T for the time-ordering operation, andψ H (r, t) denotes the field operator in the Heisenberg picture in the basis of full Hamiltonian H which contains kinetics, light-matter and many-body interaction terms. Using a standard quantum-field formalism, the effects of the many-body interacting can be conveniently cast in terms of the many-body self-energyΣ(1, 2). Using Dyson recursive relation, the full field-dependent and interacting Green's function is given in terms of a field-dependent self-energy,Σ, and of a bare Green's function,Ĝ 0 , as followŝ where we denoted 1 + ≡ (r 1 , t 1 + 0 + ) and "tr" stands for the "trace" operation over all spinor indexes i.e. tr[ÂB] = ss [A ss B s s ]. The bare one-photon current vertex is hence obtained in terms of variational derivatives of non-interacting Green's function versus the gauge field:Λ Because of the linear dispersion of the Dirac model, the two-and three-photon current vertices are of course null at the non-interacting level. However, as we are going to see, the presence of a field-dependent self-energŷ Σ(1, 1 ; A) in Eq. (5) breaks the linear dependence of the inverse Green's functionĜ −1 (1, 1 ; A) on the external field, and it is expected to give rise to higher order n-photon current vertices. More precisely, such nonlinear processes can be computed in terms of multi-point correlation functions for n-photon vertex operators : The possibility of obtaining a computationally affordable expression forΛ n (1 , 1 ; 1, . . . , n) depends of course on the specific characteristics of the scattering source.
In the following we focus on the role of elastic disorder/impurity scattering which, along with a particularly simple structure allowing for a direct computation, preserves in Dirac materials the fundamental frequencydependence of the self-energy, which is a crucial ingredient in determining the nonlinear electromagnetic response. A. One-particle impurity self-energy 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, so that V imp (r) = 0 and the effective scattering potential reads where the average . . . imp is meant over all the impurity configurations, and where V imp parametrizes the strength of impurity scattering.
In the absence of external fields, the lowest-order selfconsistent Born self-energy reads: where γ imp = n imp V 2 imp , S cell is the two-dimensional unitcell area, and the variable z lies in the complex frequency space. Due to the isotropic impurity scattering the selfenergy spinor structure is trivial asΣ(z) = Σ(z)Î and therefore the Green's function can be explicitly written as followsĜ where where S(z) = z + µ 0 − Σ(z). The introduction of a high-energy (ultra-violet) cut-off is an unavoidable requirement of Dirac models. There is however a relative large degree of freedom in the way how to introduce it, and particular care is needed in order to avoid spurious results and to preserve physical consistencies, like Ward's identities, and gauge invariance. Dimensional regularization has proven to be a formidable tool to ensure that physical correctness is preserved [22,23]. We consider the evaluation of the disorder self-energy which displays a primary diverging integral. In arbitrary D dimensions, we have thus Note that the above integral in D dimensions can be solved in terms of Euler's Gamma-function, Γ E (z), by utilizing the following identity [23] where W is a proper ultra-violet energy cut-off, and where we use the prescription lim →0 1/ ≡ ln[W ]. Eventually, we obtain the following self-consistent formula for the self-energy where U is a dimensionless parameter representing the electron-impurity scattering strength: In a direct comparison with graphene, we have S cell = √ 3a 2 /2 and v = 3at 0 /2, where a ≈ 0.246 nm is the lattice constant and t 0 ∼ 3eV is the nearest neighbor hopping energy. In order to preserve the number of states, S cell defines also an effective 2D Brillouin zone, V BZ = 4π 2 /S cell which, employing the isotropic symmetry, defines a momentum cut-off k c , πk 2 c = V BZ , and a natural ultra-violet energy cut-off W = vk c for the Dirac linear dispersion. Using the above parameters for graphene we get W = 7.2 eV.
The above procedure of dimensional regularization is employed in similar way in the evaluation of all the momentum integrals in the present work.

III. MANY-BODY-DRIVEN MULTI-PHOTON VERTEX GENERATION
The analytical expression of the functional dependence of the one-particle self-energy on the external field allows, within the spirit of a Baym-Kadanoff approach, the derivation of a closed set of self-consistent equations governing the transport properties at the chosen (linear or nonlinear) order. For Dirac materials, like graphene, where the second order response vanishes by symmetry, a particular interest is paid on the third-order response, and, within this framework, on the third-harmonic generation. At the non-interacting level, vertex corrections are null and the third-harmonic generation is governed by the well-known "square" diagram with bare one-photon current vertices at the corners [4].
Things are much more complex in the presence of many-body interactions where the intrinsic dependence of the self-energy on the frequency and on external fields triggers in novel nonlinear effects which are not predictable at the non-interacting level or within a phenomenological model using a constant (frequencyindependent and external-field-independent) one-particle scattering rate.
A careful investigation of the many-body effects driven by disorder scattering, at the lowest-order self-consistent Born level, is remarkably enlightening since it preserves all the relevant nonlinearity but with a particular simple expression for the self-energy which results to depend linearly on the fully interacting Green's function in the presence of external field: Here, V (1, 2) stands for the many-body interaction potential where for the impurity-driven the interaction potential is given by Eq. (9). On this basis, after performing a lengthy but straightforward algebra, we can construct a diagrammatic theory for the third-order response function of graphene as a two-dimensional Dirac material. Feynman diagrams for the third harmonic response function χ (3) (1; 2, 3, 4) are depicted in Fig. 1. Here solid lines represent Green's functions, wavy lines the incoming/outcoming photons, and the empty/filled symbols the bare/renormalized nphoton vertices, where n can be identified by the number of attached wavy-lines (photons).
The corresponding third-harmonic optical conductivity σ THG (ω) can be hence obtained after analytical continuation iω m → ω + i0 + as: After a straightforward algebra one can obtain thus [18] For sake of compactness, we use here a short notation where P ν0ν1ν2ν3 = P ( 0,ν0 , 1,ν1 , 2,ν2 , 3,ν3 ), where j,νj = + jω + jη νj (j = 0, 1, 2, 3) and where η νj is a vanishingly small quantity with η νj > 0 if ν j = R and η νj < 0 if ν j = A. Notice also that P AAAA = (P RRRR ) * . As depicted in Fig. 1, the P -function contains four different contributions, P = P 1 + P 2 + P 3 + P 4 , associated   respectively with square (P 1 ), triangles (P 2 , P 3 ) and bubble (P 4 ) diagrams. More explicitly we can write: where α = e 4 v 2 N f /2π 2 . The sum over spin and valley indices just leads to an overall degeneracy factor N f = N s N v where N s = 2 and N v = 2. The function Ω 1 (z 0 , z 1 , z 2 , z 3 ) represents the square diagram neglecting vertex renormalization, and represents thus the onephoton Bethe-Salpeter renormalization factor which is discussed in details in Appendix A. In similar way we can write the contributions of the two triangles diagrams as  where λ 2 (z 0 , z 1 , z 2 ) is the lowest order two-photon current vertex (Fig. 3a), and Q 2 (z i , z j ) = Λ 2 (z i , z k , z j )/λ 2 (z i , z k , z j ) is twophoton Bethe-Salpeter renormalization factor (see Appendix B). Furthermore, we can also express the triangle diagram as: where Finally we can express the bubble term as: where Q 3 (z i , z k , z l , z j ) = Λ 3 (z i , z j )/λ 3 (z i , z k , z l , z j ) is three-photon Bethe-Salpeter renormalization factor (see Appendix C), with λ 3 (z 0 , z 1 , z 2 , z 3 ) being the lowest order three-photon vertex function (Fig. 3b), and A close inspection of the topological structure of diagrams for the three-photon vertex (see Fig. 3b and Appendix C) permits further simplifications as: Using Q 3 = Q 1 = 1/[1 − U X 1 ] (see Appendix C), we find the following result for the total P -function We can see thus that the net impact of three-photon vertex diagram in the third-order response function simply leads to the appearance of the three-photon renormalization factor Q 3 (z 0 , z 3 ) on the contribution of the other diagrams. The explicit expressions of λ n , and Q n , and Ω n are provided in great details in Appendix sections for the one-, two-, and three-photon vertex Bethe-Salpeter renormalizations.
Equipped with the all analytical expressions needed for the computation of optical properties of third-harmonic generation response, in the following Section we present numerical results for the low-energy intraband thirdharmonic conductivity of two-dimensional Dirac modelling of graphene.

V. RESULTS AND DISCUSSION
Based on a mere dimensional analysis, we can conveniently express the zero-temperature third harmonic conductivity of graphene in terms of a dimensionless function f 3 (x, y, z): where µ = µ 0 − ReΣ(0) is the renormalized chemical potential, Γ = Γ(0) = −ImΣ(0) is the Fermi surface scattering rate, σ 0 = e 2 /4 is the universal conductivity including the spin and valley degeneracy, and E 0 = πt 0 / √ 3ea ≈ 22.0 V/nm is a characteristic electricfield scale determined by the inter-atomic hopping energy t 0 and by the lattice constant a. In the dc limit ω → 0, lim x→0 f 3 (x, y, z)/x 4 , we recover the transport regime discussed in Ref. [18]. The appearance of the bare chemical potential µ 0 in the definition of x, and of the renormalized one µ in y is dictated by the different role of µ 0 and µ in governing the dc and optical properties, as discussed in more details below. Besides the obvious role of the real part of the third harmonic generation, the imaginary part of the conductivity bares also a strong relevance. As a matter of fact, the THG efficiency η THG scales indeed as At high-frequency regime, the third-harmonic generation response function χ (3) THG (ω) is described by purely interband transitions, χ (3) THG,inter (ω), giving rise to steplike functions at ω ≈ 2µ 0 /n corresponding to n = 1, 2, 3 photon resonances where the three-photon resonance at ω ≈ 2µ 0 /3 defines the interband optical edge. The low frequency ω µ corresponds to the purely intraband regime of third-harmonic optical conductivity. In addition to this structure, in the intermediate frequency range, here, we obtain novel incoherent four and fivephoton transitions with mixed intra-interband characters.

A. Four and five-photon incoherent transitions
The full quantum treatment of many-body interaction with Green's functions and Kubo formalism predicts intriguing novel spectral features. In Fig. 4a,b we plot the real and imaginary parts of the THG function THG versus the laser-frequency ω in the whole intra-to inter-band frequency ranges for different values of the scattering strength U . In the almost clean limit (U = 0.001) the dominant features are the the multi-photon interband resonances at ω ≈ 2µ 0 /n with n = 1, 2, 3, which appear as smeared structures by the interaction in the real and imaginary parts. Note that on this scale the purely intraband features are not visible since their magnitude scales as ω 4 . Upon increasing the scattering strength, besides the obvious smearing of the interband features, we notice the appearance of novel resonances below the interband edge ω ≈ 2µ 0 /3, i.e. in the intraband range. A closer look reveals that such spectral features occurs at frequencies ω = 2µ 0 /n with n = 4, 5 (purple vertical dashed lines in Fig. 4), namely at energies corresponding to four-and five-photon resonances, respectively.
The very possibility of observing four-and five-photon transitions in the third order conductivity is quite surprising and it calls for a further deeper investigation. Also worth being noticed is the fact that, although a full many-body treatment is here enforced, the n-photon resonances (n = 1, . . . , 5) occurs at the energies ω ≈ 2µ 0 /n dictated by the bare chemical potential µ 0 , rather than by the effective renormalized one µ. This puzzling result has not been detected previously in literature since in non-interacting models as well as in phenomenological models where only a constant scattering rate Γ (imaginary part of self-energy) is included, no renormalization of the chemical potential is operative, and µ = µ 0 .
The appearance of four-and five-photon transitions in the third order conductivity can be rationalized by considering at the simplest level the square diagrams depicted in Fig.  1 furthermore neglecting the Bethe-Salpeter vertex renormalization.
In this case, according Eq. (22), the kernel response function P in Eq, (21) will read simple P (z 0 , z 1 , z 2 , z 3 ) = P 1 (z 0 , z 1 , z 2 , z 3 ) ∝ Ω 1 (z 0 , z 1 , z 2 , z 3 ). The explicit expression of Ω 1 (z 0 , z 1 , z 2 , z 3 ) is long and cumbersome and it is provided in Appendix C. The main feature in regards with the present issue is that it depends as: .
Accordingly Eq. (34), multi-photon transitions can occur when Re[S(z i )+S(z j )] = 0 with z j = +j ω. Neglecting for the moment the contribution of the self-energy, this implies that = −µ 0 − (i + j) ω/2. We assume for simplicity µ > 0, T = 0 and we focus on the RRRR channel. Enforcing the boundary conditions −3 ω ≤ ≤ 0 which originates from the factor "n F ( )−n F ( +3ω)" in Eq. (21), we get that possible n-photon inter-band transitions can occur at ω = 2µ 0 /n with n = 6 − (i + j), and, considering i, j = 0, 1, 2, 3 with i = j, we obtain the possible values n = 1, 2, 3, 4, 5. Similar analysis can be applied for other channels where the possibility of detecting nphoton transition is dictated by the integration window over (see Fermi function prefactors in Eq. (21) for different channels) and by the retarded/advanced character of the complex frequencies z i involved is the transition. More explicitly, it can be shown than n-photon transitions in each channel have a finite spectral weight only when z i , z j have the same (retarded or advanced) character, in similar way as it occurs in the linear optical response. With such roadmap, we can analyze theoretically the possible appearance of multi-photon transitions in each separate contribution P RRRR , P ARRR , P AARR , and P AAAR . Our theoretical predictions are summarized in Table. I, and the numerical results are shown in Fig. 5a,b, in excellent agreement with each other.
TABLE I. Predicted n-photon transitions ω = 2µ0/n in the third-order optical conductivity. The transitions marked with asterix are expected to have null spectral weight the for elastic scattering here considered but they might gain a finite spectral weight in the presence of inelastic scattering (see discussion in the text).
Channel -range at T = 0 n-photon trans. RRRR −3 ω ≤ ≤ 0 n = 1, 2, 3, 4, 5 ARRR − ω ≤ ≤ 0 1 * AARR −2 ω ≤ ≤ − ω n = 1, 2 * , 3 AAAR −3 ω ≤ ≤ −2 ω n = 1, 2, 3, 4, 5 The total nonlinear response is determined by the sum of all the channels. As we can see in RRRR and AAAR channels show step-like transitions at the four-and five-photon resonances, i.e. at ω ≈ µ 0 /2 and ω ≈ 2µ 0 /5, but with opposite sign. The sum of these two contributions would exactly cancel out in the clean limit U = 0. Such cancellation is however just partial for finite U (or finite Γ), leaving finite spectral structures at ω ≈ µ 0 /2 and ω ≈ 2µ 0 /5, as seen in Fig. 5c. The behavior as a function of the scattering strength is shown in Fig. 5d. As mentioned above, in the clean limit U → 0 the multi-photon resonances at n = 4, 5 in the individual channels P RRRR , P ARRR , P AARR and P AAAR cancel out exactly and they are thus absent in the total response. However, such cancelation is not perfect in the presence of a finite electron-impurity scattering, leaving residual spectral structures at frequencies corresponding to the four-and five-photon resonances. The spectral weight of these multi-photon structures scales with the impurity scattering itself. The absence of four-and five-photon transitions in the clean (non-interacting) limit reveal thus that these transitions are in fact incoherenttransitions with mixed inter-and intra-band characters. The partial intra-band character is highlighted by the transition weight being proportional to the relaxation rate Γ, whereas the partial inter-band character of such features is pointed out by their lying at finite frequency, i.e. at exactly the four-and five-photon resonance energy with the peak positions not affected by of U . It is worth emphasizing that these incoherent-transitions at finite frequency pinned by µ 0 are a peculiar property of nonlinear optical conductivity and they do not emerge in the linear optical conductivity. Now, we focus on the role of the many-body self-energy renormalization in determining the spectral features of the optical third-harmonic generation response. As discussed above, n-photon transitions can be theoretically identified by enforcing the condition Re[S(z i ) + S(z j )] = 0 together with = min , max , where min , max are the lower and upper energy integration limits for each channel. With such prescription, neglecting the manybody self-energy, the n-photon transitions occur at ω ≈ 2µ 0 /n, where µ 0 is the bare chemical potential. It is however straightforward to check that the same holds true also in the presence of (frequency-dependent) elastic scattering driven by disorder/impurity preserving the mirror symmetry with respect to the Dirac point. The analysis of linear optical conductivity is enlightening on this point. In similar way as in the third-order response function, at T = 0 the edge of the interband optical transitions is determined by the conditions Re[S( ) + S( + ω)] = 0 and Re[S( ) + S * ( + ω)] = 0, respectively, together with the constraints = 0, = − ω determined by the window of energy integration over . We obtain thus that the edge of optical interband transitions is determined by with a spectral weight I sw that scales as: where ν, ν = A,R. The elastic impurity scattering selfenergy respect the following symmetry relation owing to symmetry of Dirac dispersion: These relations imply in a direct way that: (i) the interband optical edge in the linear optical conductivity is determined only by µ 0 and not by the renormalized chemical potential µ; (ii) only the retarded/retarded channel is responsible for the n = 1 photon transition observed thus in linear optics. Similar argumentations holds true in a straightforward way in nonlinear optics where we conclude that: (iii) the n-photon transitions are determined only by µ 0 and not by the renormalized chemical potential µ; (iv) only retarded/retarded or advanced/advanced transitions show a sizable spectral weight and can be observed thus in the optical features. The strict validity of these symmetry relations is affected in the presence of inelastic scattering where the even/odd symmetries with respect to the Dirac points are lifted. Considering however the realistic case of lowenergy inelastic scattering (e.g, phonons), the breaking of the symmetry relations in the imaginary part in Eqs. (38)-(39) is limited to a narrow shell around the Fermi level, affecting thus only the ω ≈ 2µ 0 resonance (n = 1). Inelastic scattering might thus reduce a bit the spectral weight of n = 1 photon transitions in the RRRR, AARR, AAAR channels and it might induce a finite spectral weight in the ARRR channel which is predicted to be null under the above symmetry conditions valid however only for elastic scattering. Inelastic scattering would affect as well the symmetry relation for the real part of the self-energy, as it can be obtained by Kramers-Kronig transform of the imaginary part of the self-energy. Since the imaginary part of the self-energy is typically affected only in a narrow region around the Fermi level, we expect deviations from the symmetry relation in the real part in Eq. (37) to be relatively small and scaling with the strength of the inelastic coupling. Multi-photon transitions might be just slightly shifted from the edge determined by µ 0 .

B. Intraband THG conductivity
The analysis of the THG optical conductivity allows to investigate in details also the low-energy spectral features associated with the pure intraband processes, which are not easily visible in the dimensionless function f 3 . A phenomenological model previously obtained [2,3] based on density-matrix formalism using a constant scattering rate can qualitatively describe low-frequency profile of THG conductivity in the Boltzmann regime where µ Γ, ω: where C > 0 is a constant and where Γ stands for a constant (frequency-independent) scattering rate. Eq. (40) captures a similar physics as the Drude term in the linear optics. In similar way as the Drude term, and unlike the interband counterpart ruled by µ 0 , the spectral properties of Eq. (40) are governed only by Γ and by the ratio ω/Γ. Intraband transitions occur at the Fermi surface and accordingly the intraband THG conductivity depends on the renormalized chemical potential µ. Remarkably, the above model predicts for any strength of Γ   Fig. 6 for the graphical clarification of parameters defined above.
In Fig. 7a we show the real part of the THG optical conductivity as obtained by our quantum conserving approach for different impurity scattering strengths. In the weak scattering limit (red curves), the four-and fivephoton edges are visible in the bump and in the dip at ω = 0.5µ 0 and ω = 0.4µ 0 , respectively. Besides these spectral structures, with mixed intra-interband character, we are able now to reveal a purely intraband low frequency structure with a strong resemblance with the phenomenological model of σ THG,Γ given in Eq. (40). Giving the complex nature of the nonlinear response, even in the intraband dc limit where the multiphoton vertex renormalization has shown to be highly relevant, [18] the robustness of the phenomenological constant-Γ model appears by no way trivial. This issue is pointed out in Fig. 7b,c where we plot the THG optical conductivity at two levels of (non-conserving) approximations: neglecting the two-photon vertex renormalization (X 2 = 0), that has been shown to be dominant in the dc limit, see panel (b); and neglecting all the vertex renormalization processes, see panel (c). For both cases we find a completely different scenario at low-frequency with respect to the phenomenological constant-Γ model and to the compelling numerical results. We remind that the constant-Γ model is itself a conserving approx-imation, where to a field-independent and frequencyindependent self-energy correspond unrenormalized vertices. We rationalize the agreement at weak scattering between the constant-Γ model and the fully many-body theory on the basis of the requirement of a conserving analysis. Quantum many-body theories, when based on non-conserving approaches, might severely give spurious results.
For stronger scattering strength U (see blue curves in Fig. 7) we enter to the quantum regime where the scattering rate Γ is comparable (or larger) than the chemical potential µ. For this case, we observe a full merging of all spectral features (pure intraband peak merges with four and five-photon incoherent structures) and morphing into a single peak as well as the nonlinear dc conductivity changes sign. This sign change is a result of twophoton vertex self-generation which was discussed previously [18]. Obviously, the phenomenological model given in Eq. (40) is no longer valid in this strong interacting regime.
In order to investigate in a deeper way the comparison of our many-body theory with the phenomenological model in the weak scattering limit, we analyze the numerical output of the fully quantum conserving theory in the low-frequency range within the ansatz of the phenomenological model of Eq. (40). More precisely, we determine in the full quantum theory, as a function of the scattering strength U , the dc limit σ  THG (ω = 0); the frequency ω 0 where the THG has the first low-frequency zero; the frequency ω max where the THG conductivity has its first maximum; the value σ   conductivity at ω = ω max . In Fig. 8a-d we compare the mutual dependence of the such characteristic spectral parameters as obtained from the numerical results and as estimated from the model in Eq. (40). We find an excellent agreement in the very weak scattering regime, but a quick deviation by increasing the scattering strength U . Keeping in mind that in the phenomenological model ω max ≈ Γ, we estimate a breakdown of the constant-Γ model for Γ ≈ 0.01µ, e.g. see Fig. 8d. On the basis of Fig. 7a, we rationalize this failure as the overlap of the purely intraband term with the five-and four-photon edge transitions, starting from ω ≈ 0.4µ 0 − 0.5µ 0 .

VI. SUMMARY AND CONCLUSION
In this paper we have explored the effects of elastic impurity scattering effect on the third-harmonic generation of graphene by using a conserving diagrammatic method that includes self-consistently both self-energy and vertex renormalization contributions. As a result of the field dependence of the self-energy, we showed that interactioninduced multi-photon vertex diagrams are relevant. In particular we have predicted the onset of novel incoherent resonances at four and five-photon transition energy with mixed intra-and inter-band character. Furthermore, we have shown that the main features of the purely intraband contribution in the full quantum theory are qualitatively comparable with phenomenological models that assume a frequency/field independent self-energy, whereas non-conserving approaches give rise to spurious results. In the terahertz regime the proximity of the fivephoton transition edge at ω = 2µ 0 /5 might affect the interpretation of experimental data in terms of a purely intraband term.
Our study, shed light on the importance of many-body interaction on the qualitative explanations of nonlinear optics in the terahertz and infrared (intra-band) regimes. Although the numerical study is performed for the thirdharmonic generation, the formalism is quite general and applicable for two-photon absorption, nonlinear Kerr effect, etc. Moreover, our study can be simply generalized to explain terahertz nonlinear response in other twodimensional materials such as transition-metal dichalcogenides and twisted bilayer graphene.

VII. ACKNOWLEDGEMENTS
H.R. acknowledges the support from the Swedish Research Council (VR 2018-04252).