Dominant Role of Two-Photon Vertex in Nonlinear Response of Dirac Materials

Using a conserving Baym-Kadanoff approach, we present a fully compelling theory of nonlinear dc response of a Dirac system to electric fields in the presence of disorder scattering. We show that the nonlinear terms are strikingly ruled by the appearance of a dominant two-photon vertex which is absent at the bare level and finite even in the weak-coupling limit. Such two-photon vertex self-generation highlights the crucial role of the frequency and field dependence of the scattering rates in the nonlinear regime. Our study reveals a novel many-body mechanism in the nonlinear response of Dirac materials whose effects are predicted to be observable.

Using a conserving Baym-Kadanoff approach, we present a fully compelling theory of nonlinear dc response of a Dirac system to electric fields in the presence of disorder scattering. We show that the nonlinear terms are strikingly ruled by the appearance of a dominant two-photon vertex which is absent at the bare level and finite even in the weak-coupling limit. Such two-photon vertex self-generation highlights the crucial role of the frequency and field dependence of the scattering rates in the nonlinear regime. Our study reveals a novel many-body mechanism in the nonlinear response of Dirac materials whose effects are predicted to be observable.
Due to their linear dispersion, k ∼ |k|, and to the underlying chiral structure, Dirac materials show a variety of exotic features that makes them a versatile platform for theoretical investigations of new physics and for application purposes. Despite the complex physics, many properties of these materials are often rationalized using concepts of non-interacting particle or semi-classical model [1][2][3][4]. For instance, a standard transport model is conventionally applied for the dc conductivity in highlydoped graphene (Boltzmann regime), where the mobility is evaluated at the non-interacting level, and the interactions enter only through the effective parameter known as transport scattering rate Γ tr [5]. At odds with the above scenario, 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 affected by many-body effects [6].
The nonlinear electromagnetic response of Dirac materials has attracted recently a considerable interest in two [7][8][9][10][11][12][13][14] and three dimensions [15][16][17][18][19][20]. Widely investigated are the nonlinear optical properties and in particular the appearing of four-wave mixing, nonlinear Kerr effect, second and third-harmonic-generation in singlelayer graphene, with remarkable technological interest [21][22][23][24][25][26][27][28][29][30]. Peculiar of Dirac material is, due to the linear dispersion, the absence of the bare two-photon-electron coupling, which should give rise to the so-called diamagnetic term. The lack of such term prompts several widely debated issues, as the validity of optical sum rules [31][32][33]. Most of the theoretical descriptions of nonlinear effects rely at the moment upon non-interacting analyses, or semiclassical approaches [23,24,[34][35][36] where, in a similar way as in the Boltzmann theory of linear response, the dominant transition processes resemble the ones of the non-interacting case and the scattering sources are accounted through effective parameters as the scattering rate Γ (or equivalently through the mean-free path l, the lifetime τ , etc.).
In this Letter, we show that the a compelling analy- * Electronic address: habib.rostami@su.se sis of the many-body physics, beyond the semi-classical approaches, can drastically change the above scenario, pointing out that different physical processes can be responsible for the relevant properties of the nonlinear dc transport. Analyzing the case of disorder scattering as a basilar benchmark example, we show how non-conserving phenomenological models of scattering intrinsically fail and high-order vertex processes must properly taken into account. More in particular, we show that, despite the bare diamagnetic two-photon vertex (TPV) being null in Dirac materials, the many-body renormalized TPV is finite and relevant and it can play a dominant role. Our results, besides providing a consistent framework for a proper analysis of nonlinear transport and optical response in realistically interacting Dirac materials, open novel perspectives for understanding and predicting new functional properties of these complex promising systems.
We consider the two-dimensional (2D) Dirac Hamilto-nianĤ k = ̵ hvσ ⋅ k − µ 0Î , where µ 0 is the bare chemical potential ruling the charge doping. For realistic purposes we consider the paradigmatic case of graphene [37] wherê σ = (τσ x ,σ y ),whereσ i stand for the Pauli matrices in the spinor space, and τ = ± stands for valley index in the Brillouin zone of graphene. In the dipole approximation the light-matter interaction can be modelled by applying the minimal coupling transformation ̵ hk → ̵ hk + eA(t) where A(t) stands for an external vector potential. The corresponding electric field is given by E(t) = −∂ t A(t). Due to the linear dispersion, the electron-photon coupling does not present a diamagnetic (two-photon) bare term but only the linear coupling: Without the loss of generality, we assume an electric field along the y axis. As scattering source we consider longrange impurity centers with standard Born impurity correlations [38][39][40]. Within this framework we can write the Born impurity self-energy in the complex frequency space:Σ(z) = γ imp ∑ kĜ (k, z) where the Green's function followsĜ(k, z) = [z −Ĥ k −Σ(z)] −1 . For isotropic scattering we get a diagonal self-energy in the spinor basis arXiv:2007.08282v1 [cond-mat.mes-hall] 16 Jul 2020  asΣ(z) = Σ(z)Î. It is well known that under these conditions the impurity self-energy, as well the Coulomb and other scattering ones, depends intrinsically on the ultraviolet energy cut-off W representing the range of validity of the Dirac model. In order to provide a conserving approach, this issue needs to be cured by means of a proper regularization [41,42]. As detailed in the Supplementary Material (SM), we employ standard dimensional regularization leading to: where S(z) = z + µ 0 − Σ(z), and U is a dimensionless parameter characterizing the strength of impurity scattering [38-40]. Conserving approaches, based for instance on a Baym-Kadanoff derivation [43,44], are fundamental in theoretical physics to ensure that compelling results are obtained. This aim is particularly important in nonlinear response since an arbitrary selection of diagrams can easily lead to spurious conclusions. The choice of the vectorpotential gauge, within the paradigmatic Born impurity scattering we consider here, permits us an exact derivation of self-consistent equations (see SM for details [45]) for all the high order processes relevant in the thirdorder response function which is the leading nonlinear term in centrosymmetric Dirac materials. The diagrammatic expression of the third-order response function is provided in Fig. 1 where, roughly speaking, empty symbols represent n-photon vertices (n = 2, 3) expressed in terms 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 complexity of the self-consistent set of equations, few relevant things are worth to be underlined here. First of all, we notice that an effective multi-photon coupling is induced by the disorder scattering source even if it is absent in the Hamiltonian at the bare (non-interacting) level (Fig. 1e,f). Second, that the relevance of each nphoton vertex is largely governed by the BS many-body resummation as depicted in Fig. 1b-d. This might lead to a reduction (screening) or to an enhancement of different multi-photon scattering depending on the Pauli structure of the corresponding photon vertex, as we discuss more extensively later.
The diagrammatic expressions in Fig. 1 represent in full generality the optical frequency-dependent thirdorder response function in Dirac materials, including third-harmonic generation, four-wave mixing, etc. For a generic interaction, the effective solution of such coupled equations on the real-frequency axis is a formidable task that does not allow for a practical solution. The focus on the isotropic disorder scattering is on the other hand particularly suitable to investigate many-body effects in nonlinear electromagnetic response since it allows for a set of equations in the Matsubara space which can be generalized in a rigorous way on the real frequency axis, using the well-known procedure of multiple branch cuts in the complex frequency space. The derivation is lengthly and cumbersome but compelling and it is summarized in the SM [45]. We consider first the dc transport limit. Without loss of generality, it is possible to express the linear and the third-order dc conductivity in terms of two dimensionless quantities: where µ is the effective chemical potential µ = µ 0 − ReΣ(ω = 0) and Γ(µ) the scattering rate Γ(µ) = −ImΣ(ω = 0), σ 0 ∝ e 2 / ̵ h the universal conductivity unit and E 0 ∝ t 0 /ea is a characteristic electric field scale determined by inter-atomic hopping energy t 0 and by the lattice constant a [37].
It is worth to stress again that Eqs.
(3)-(4) are tied together since they must descend in a compelling way from a common approximation for the self-energy. Heretofore, although many approaches for the self-energy have been discussed for the linear response, the third-order response has been analyzed only in the simplistic case of a phenomenological constant scattering rate Γ(µ) = Γ. Since such phenomenological self-energy does not depend on the applied external field, the third-order response function reduces to the first "square" diagram of Fig. 1a dropping all the vertex renormalization processes, i.e. replacing the filled circles with empty ones (= bare electron-photon coupling). A similar scheme can as well be employed for the linear response. Under these ultrasimplified conditions, one can see that the linear and third-order dc transport depend uniquely on the semiclassical parameter x = µ/Γ(µ), i.e. f 1 (x; y) = f 1 (x), f 3 (x; y) = f 3 (x). An analytical expression for the functions f 1 (x), f 3 (x) is obtained in the SM [45]. In particular, in the Boltzmann regime one gets results f 1 (∞) ≈ 2µ/πΓ, f 3 (∞) ≈ −3πΓ/32µ, implying that nonlinear effects lead to a reduction of the dc conductivity in the Boltzmann regime. A similar analysis is performed in the quantum regime, giving [45] f 1 (0) = 8/π 2 , f 3 (0) = 2/5, meaning that nonlinear effects should yield an enhancement of the dc conductivity in the quantum regime.
The above predictions, based on the phenomenological model of a constant scattering rate Γ, are challenged when many-body effects are computed in a compelling conserving scheme. In Fig. 2 we show the characteristic dc transport [46] function f 3 as a function of the semiclassical parameter x = µ/Γ(µ), from the extreme quantum limit (x ≪ 1) to the Boltzmann regime (x ≫ 1). From the computational point of view, since the presence of a finite cut-off energy scale W , the quantum limit can be conveniently investigated by fixing U and varying µ, whereas the Boltzmann regime is more eas- ily spanned by fixing µ and varying U . Let's discuss first the Boltzmann regime (Fig. 2b). We notice that a compelling many-body analysis recovers qualitatively (but with an increase in the magnitude of factor 10) the predictions of the phenomenological constant-Γ model On the other hand, in the quantum regime x ≪ x * (Fig. 2a), f 3 shows a significant dependence on µ 0 , signalizing that the nonlinear dc transport properties are no more governed uniquely by the semiclassical parameter x but that the detailed value of µ 0 (or conversely, of U ) starts playing a relevant role. Some striking things are worth being pointed out: (i) counterintuitively, the third-order contribution to the dc transport appears to be magnified approaching the clean limit U → 0; (ii) there are two isosbestic points (i.e. x-points where f 3 does not depend on U ) coinciding in a very good approximation with the zeroes of the f 3 function; (iii) whereas the phenomenological constant-Γ modelling predicts a a well-determined positive sign of the nonlinear dc correction in the quantum regime (implying an increase of the total conductivity), the sign of the nonlinear terms of the full conserving many-body theory in the quantum regime is not univocally determined, presenting a positive region in the crossover range and a negative sign in the extreme quantum limit.
We can rationalize points (i)-(ii) by assuming that in the quantum regime the nonlinear characteristic dc transport function f 3 (x, U ) can be factorized as [45]: (with γ > 0), where the strength U of the interaction rules governs the intensity of the third-order dc transport, while the semiclassical parameter x seems to dictate the sign of the third-order correction. To assess the meaningfulness of such description we plot in Fig. 3a in a log-log scale the absolute value of the function f 3 (x, U ) versus U for a representative case x = 0.03 in the quantum regime. We find a perfect agreement with a scaling behavior f 3 (x, U ) ∝ 1/U γ with γ slightly smaller than 2 (γ = 1.82) signalizing that in the clean limit the third-order dc transport is expected to be dominant with respect to the linear one. As detailed in [45], a similar analysis is valid in the whole quantum regime.
In order to gain a full understanding of these novel features, we analyze separately in Fig. 3a,b the relevance of each family of diagrams contributing to the total third-order conductivity as depicted in Fig. 1a. We can thus realize that the contribution of the conventional "square" diagram (which is the only one present in the non-interacting case and for the phenomenological damping model), is essentially marginal, as well as the contribution of the last "bubble" associated with the renormalized three-photon vertex. The dominant role is instead played by the "triangle" diagrams containing the renormalized two-photon vertex. A quantitative analysis shows that each family of diagrams obeys Eq. (5) with an approximately integer exponent (i.e. γ ■ ≈ 2, γ ▲ ≈ 2, and γ • ≈ 1 for the square, triangle and the bubble diagrams). The dominance of the triangle diagrams results thus in an exponent very close to 2 (γ ≈ γ ▲ ). The selfconsistent BS renormalization of the TPV (Fig. 1c) is a crucial ingredient in such novel scenario. This can be assessed in Fig. 3a,b where one can see that, once neglected the BS renormalization, the contribution of the triangle diagrams results to be of the same order (even smaller) of that of the conventional square diagram. The dominant role of the TPV renormalization appears even more evident by investigating the scaling of the characteristic third-order dc transport function f 3 (x, U ) versus U . As depicted in Fig. 3a, once replaced the BS renormalized TPV (Fig. 1c) with the "bare" one ( Fig. 1e), the third-order dc conductivity scales as 1/U (γ △ ≈ 1), with an additional sign change change, as shown in Fig. 3b. This means that the BS renormalization of the TPV gives rise in the quantum regime to an additional dependence ∼ 1/U that diverges in the clean limit. The n-photon vertex matrix structure, which readsΛ n = (−evσ y ) n Λ n , plays a crucial role in the relevance of the BS renormalization effect. The impressive effect is peculiar of the TPV renormalization Λ 2 = Λ and does not appear in the BS renormalization of the one-three- This different impact can be traced down to the different structure in the Pauli space. As detailed in [45], we get indeed In the quantum regime, one can thus show that in the dc limit U X 1 ≈ U X 3 ∝ U , whereas U X 2 ≈ 1 + O(U ), so that resulting in an effective divergence of Λ 2 /Λ (0) 2 in the dc limit at zero temperature and in the clean limit (U → 0).
The impact of the two-photon renormalization can be understood in more details by investigating the two- 10 dc | factor with E = 1mV/nm versus chemical potential µ and relaxation rate Γ(µ) in the full conserving models. The sign of σ where and where z 1 and z 2 are the electronic frequencies in the complex plane. Particularly enlightening is the analysis of the retarded-retarded (RR) channel. In the dc limit (ω → 0) Q RR 2 ( , + ω) is simply given by where ω is the photon energy and is the electronic energy from the Fermi surface. For low-energy excitations = 0 we have thus Q RR 2 The Boltzmann regime is achieved as µ 0 ≫ 2U S(0). In the clean limit U → 0 we get thus S(0) = µ 0 and Q RR 2 The quantum regime is on the other hand characterized by µ 0 ≪ 2U S(0), and we get Q RR 2 = 1/2U , leading thus to a huge enhancement (divergence) in the clean limit. A similar behavior Q 2 ∝ 1/U appears also in the retarded-advance (RA) channel, although is a more delicate way. In the dc limit, we find indeed the leading term in Q RA 2 ( , + ω) ≈ −i2Γ( )/ω which shows a divergence as a function of the photon energy ω → 0. Once plugged this behavior in the response function associated with the "triangle" diagrams containing the renormalized BS two-phonon vertex, χ ren. triangles (ω) ∼ Γ( ) χ unren. triangles (ω)/ω the divergence is ω implies that χ unren. triangles (ω) must be expanded to a higher order in ω, involving the second derivative S ′′ (0) = −1/[4U 2 Γ(0)] and higher orders. As a net result, the divergence 1/ω in Q RA 2 ( , + ω) is reflected in a consequent divergence −1/U in the response function, similarly as for the RR channel, but with a negative sign. The balance between RR and RA terms determines the change of sign of the third-order dc conductivity as a function of x n the quantum regime. We must stress that the huge enhancement of the third-order dc transport is governed by the dominant role of the twophonon vertex renormalization. Since Λ (0) 2 scales as U and 1/[1 − U X 2 ] scales as 1/U , such enhancement can be regarded as TPV self-generation (Λ 2 ≠ 0) which survives in the weak-coupling (clean) limit U → 0.
The net result on the dc transport is summarized in Fig. 3c where we plot the sign and the magnitude of the third-order conductivity for a given electric field E = 1 mV/nm normalized to the linear order conductivity, g = |σ dc |, in the physical space of the effective chemical potential µ and scattering rate Γ(µ), as they can be obtained directly in an experimental way. The Boltzmann regime corresponds thus to the right-lower corner whereas the extreme quantum regime (x → 0) is recovered in the left-upper corner. As noticed before, at odds with the predictions of the phenomenological model, we find that the third-order conductivity σ (3) dc is negative not only in the Boltzmann regime, but also in the quantum regime. Note that the zeroes of the third-order dc conductivity, see x 1 and x 2 in Fig. 2, appear in this plot as straight dashed lines. This is a consequence of the factorizable expression for the characteristic third-order dc transport function as shown in Eq. (5). We mark with tiny dotted in this plot the regions where the third-order terms start to be relevant g ≈ 0.1 and where they become of the same order than the linear dc term g ≈ 1 [47].
Note that in the Boltzmann regime, by increasing µ → ∞, one should correspondingly need infinitesimally small Γ → 0 in order to detect third-order corrections, while in the quantum regime µ → 0 third-order corrections appear to be relevant up to large values of Γ dictated only by the applied electric field (and by the ultra-violet cut-off of the Dirac model). An alternative and maybe more direct way to assess the relevance of the nonlinear conductivity is to evaluate in graphene [37] as a paradigmatic 2D Dirac material the critical electric field E max above which thirdorder corrections to the dc transport become of the same order of the linear term, |σ dc . In the phenomenological constant-Γ model we obtain E max = αE 0 with α ≈ 1.47 × µΓ/t 2 0 in the Boltzmann regime and α ≈ 1.42 × Γ 2 /t 2 0 in the quantum limit. These values can be compared with the estimates for the full conserving theory that gives α ≈ 0.32 × µΓ(µ)/t 2 0 in the Boltzmann regime and α ≈ 1.
[48] a roughly constant value Γ ≈ 15 meV was estimated in the wide range µ ∼ 0−200 meV. With these values the phenomenological model would estimate a critical field E max ≈ 10.8 mV/nm for µ ≈ 200 meV and E max ≈ 74.7 µV/nm for µ = 0 with a quantum-Boltzmann crossover at µ * ≈ 150 meV, whereas the full conserving theory predicts E max ≈ 2.35 mV/nm for µ ≈ 200 meV and E max ≈ 0.051 mV/nm for µ = 0, in a more observable range.
In conclusion, in this Letter, we have presented a fully conserving theory of nonlinear transport response in 2D Dirac materials. Our results show that the previous analyses in literature, based on phenomenological scattering models, can be qualitatively (but not quantitatively) reliable in the Boltzmann regime but they completely fail in the quantum regime. We have shown that, in a wide region of the phase diagram, close to the neutral point, the nonlinear dc transport response is dominated by novel physical processes where the two-photon vertex, absent in the bare Dirac Hamiltonian, plays a relevant role. It should be furthermore stressed that our results, focused on the dc limit, implies that the current knowledge about the nonlinear optical response in the terahertz regime should be deeply revised. Our work opens new scenarios for a deep understanding of the electromagnetic response of Dirac systems, whose relevance ranges from condensed matter to high-energy physics. York, 1984).
[37] For numerical calculations we consider a two-dimensional Dirac material with realistic parameters for graphene. In particular we consider a honeycomb lattice with interatom distance b = 1.42Å, with lattice parameter a = 3b = 2.46Å and nearest-neighbor hopping t 0 ≈ 3 eV, corresponding to a Dirac velocity v = 3t 0 a/2 ̵ h ≈ 10 6 m/s. The size of the unit-cell reads thus S cell = 3a 2 /2. We consider double spin and valley degeneracy N s = N v = 2. In order to preserve the size of the Brilloiun zone  dc is negative should not be regarded as an onset of negative total dc conductivity, rather as a sign that the expansion at the third-order in E starts to be a poor approximation and higher order terms in powers of E must be included in the analysis. 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 [41,42]. Here we show how such approach provides a consistent framework for evaluating self-energy and susceptibilities. In our work we employ dimensional regularization that endures gauge invariance. As a benchmark example, and for the sake of simplicity, 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)⟩ = 0 and the effective scattering potential reads (S1) Note that the average ⟨. . . ⟩ imp is meant over all the impurity configurations and γ imp = n imp V 2 imp 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 readŝ Note that S = N cell S cell is the system area with S cell being the unit-cell area. Due to the isotropic impurity scattering the self-energy spinor structure is trivial asΣ(ik n ) = Σ(ik n )Î and therefore the Green's function can be explicitly written as followsĜ where z = ik n . Accordingly, we find where S(z) = z + µ 0 − Σ(z). In arbitrary D dimensions, we have 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 [42] (S6) Therefore, we find Now we set D = d − where d = 2 is the physical dimension and → 0. Note that Γ E ( /2) ≈ 2/ for → 0 and      Note that we use the prescription lim →0 1 ≡ ln[W ] where W is the ultra-violet energy cut-off. Eventually, we obtain the following self-consistent formula for the self-energy (S10) After solving the above relation self-consistently, the real, ∆ = Re[Σ(ω = 0)], and imaginary, Γ = −Im[Σ(ω = 0)], parts of the self-energy at the Fermi surface are depicted in Fig. S1. For the undoped regime and at the Dirac point, The above procedure of dimensional regularization is employed in similar way in the evaluation of other momentum integrals in this study.
(S12) We use a shorthand notation 1 ≡ (r 1 , t 1 ) for the space-time coordinate. Using Dyson recursive relation, the full field-dependent and interacting Green's function is given in terms a field-dependent self-energy,Σ, and a bare Green's function,Ĝ 0 , as followsĜ The inverse of bare Green's function readŝ We assume an external gauge field along y axis, A(1) = A(1)ŷ. The one-photon current vertex is given in terms of variational derivative of bare Green's function versus the gauge field: The thermodynamic physical current, J(1; A) = J(1; A)ŷ, in Dirac systems reads Note that 1 ′+ ≡ (r 1 ′ , t 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 ]. 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's function [43,44] followsĜ where ⟨. . . ⟩ stands for the thermodynamical average and T is for the time-ordering operation.

A. Conserving linear response theory in Dirac systems
We define the linear susceptibility as follows Using Eq. (S17), the linear susceptibility reads Using Eq. (S14), we find Since the self-energy depends on the external potential only through its dependance on the Green's function, we can write down We define Bethe-Salpeter kernel as bellowΞ Therefore, using Eq. (S21) the self-consistent Bethe-Salpater relation for the one-photon vertex function followŝ (1,5)Λ 1 (5,6; 2)Ĝ(6, 1 ′ ) (S26) Note that the dressed one-photon vertex function are defined aŝ For the sake of simplicity, we extend the definition of space-time parameter to include also spinor indexes as1 = (r1, t1, s1) and, for instance, we can drop "ˆ" symbol inΛ 1Ĝ instead use Λ 1 G since the spinor multiplication is taken into account when we replace ∫ d1 ∑ s1 → ∫ d1. Moreover we use shorthand notation for C = ∫ d1A(. . . ,1)B(1, . . . ) as C = AB. We use this compact notation from now on Note that Tr[. . . ] in the above formula stands for the sum over un-contracted spinor index. In the compact notation 1 and 2 are space-time symbols.

B. Conserving third order response theory in Dirac systems
Third-order response function is given by Using GG −1 = 1 and Eq. (S21), we evaluate the third derivative of the Green's function versus vector potential component, Using Eq. (S27), we have Once more we use GG −1 = 1 and perform second derivative of the Green's function which leads The two-photon vertex function is defined as follows where P(2; 3; 4) stands for all six permutations among 2,3, and 4 space-time coordinates. Therefore, the third-order response function reads Tr Λ The two-photon vertex function is defined as follows The second derivative of the bare Green's function vanishes in Dirac systems which implies the absence of bare twophoton vertex. However, an interaction induced two-photon vertex function is obtained owing to the field-dependent self-energy. Since the self-energy depends on the external field only through the dependence on the Green's function, we have Using Eq. (S35), we obtain

Interaction induced three-photon vertex
The three-photon vertex function is defined as follows Λ 3 (2, 3, 4) = 1 2 The third derivative of the bare Green's function vanishes in Dirac systems which implies the absence of bare threephoton vertex. However, an interaction induced three-photon vertex function is obtained owing to the field-dependent self-energy. Since the self-energy depends on the external field only through the dependence on the Green's function, we have Using Eq. (S37), we obtain To summarise, the formal derivation given in the current Section guides us in constructing a conserving diagrammatic theory for the third-order response function in Dirac systems. Lengthy mathematical relations for different 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 the main text. Quantitative evaluation of these diagrams are explicitly discussed with great details in the next Section.

S3. ANALYTICAL EXPRESSIONS FOR NONLINEAR CONDUCTIVITY
We present here the analytical expression of the third-order optical response function which in the Matsubara space can be formally written as: THG (m 1 , m 2 , m 3 ) = 1 β n P (n, n + m 1 , n + m 1 + m 2 , n + m 1 + m 2 + m 3 ), where n stands for a short notation n = iω n (ω n being a fermionic frequency), and m = iω m (ω m being a bosonic frequency). For the third-harmonic generation we have m 1 = m 2 = m 3 = m and After a straightforward algebra we perform the Matsubara summation and analytic continuation as iω m → ̵ hω + i0 + , we obtain (see Section S6) Note that P RRRR ( 0 , 1 , 2 , 3 ) with j = + jω + iη j means that all frequency arguments are in the retarded channel (i.e. (η j → 0 + ) while P ARRR implies that the first argument is in the advanced channel, η 0 → −0 + , but the other are retarded, η j≠0 → 0 + . Therefore, we have P AAAA = (P RRRR ) * . By knowing the response function, the third-harmonic optical conductivity reads

A. Different diagram contributions
Using Baym-Kadanoff analysis summarised in the previous section, we construct diagrams for the third-order response 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 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: 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. Note that Ω 1 (z 0 , z 1 , z 2 , z 3 ) is the bare square diagram in the absence of vertex renormalization, and the one-photon vertexΛ 1 (z i , z j ) = (−evσ y )Λ (0) 1 = 1 is the bare one-photon vertex function and Q 1 (z i , z j ) is the one-photon Bethe-Salpeter renormalization factor which is given in the following subsection. In similar way we can write the contributions of the two triangles diagrams, namely where Note thatΛ 2 (z 0 , z 1 , z 2 ) = (−evσ y ) 2 Λ (0) 2 (z 0 , z 1 , z 2 )Q 2 (z 0 , z 2 ) is the two-photon vertex function where Λ (0) 2 (z 0 , z 1 , z 2 ) is the unrenormalized two-photon vertex function and Q 2 (z 0 , z 2 ) is the two-photon Bethe-Salpeter renormalization factor (see subsections below). We have also the further triangle diagram: where Finally we make the bubble term explicit: is the three-photon vertex function with Λ (0) 3 (z 0 , z 1 , z 2 , z 3 ) being the unrenormalized three-photon vertex function and Q 3 (z 0 , z 3 ) is the three-photon Bethe-Salpeter renormalizationfactor (see subsections below). The explicit expressions of Λ z 1 , z 2 , z 3 ), and Q n=1,2,3 (z i , z j ), Ω 1 (z 0 , z 1 , z 2 , z 3 ), Ω 2 (z 0 , z 2 , z 3 ), and Ω 3 (z 0 , z 1 , z 3 ) 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 readŝ Λ 1 (p, p + q; n, n + m) −Λ (0) 1 (p, p + q; n, n + m) = γ imp kĜ (k, n)Λ 1 (k, k + q; n, n + m)Ĝ(k + q, n + m), (S58) where m ≡ iq m and n ≡ ik n = ip n stand for the bosonic and fermionic Matsubara frequencies, respectively. Note that in the integrand we have shifted the dummy momentum k as k + p → k and therefore we can see that vertex correction does not depends on the fermion momentum p. For the optical (or dipole) approximation we have q = 0. Therefore, the Bethe-Salpeter relation for the one-photon vertex function readŝ Λ 1 (n, n + m) =Λ (0) 1 + γ imp kĜ (k, n)Λ 1 (n, n + m)Ĝ(k, n + m) . (S59) Note that we haveΛ (0) 1 = δĜ −1 0 /δA| A→0 = −evσ y whereĜ 0 stands for the non-interacting Green's function. We assume the following ansatz for the vertex function Λ 1 (n, n + m) = aÎ + 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 = 0 and eventually the following result for the vertex functionΛ 1 = (−evσ y )Λ 1 where Note that Λ (0) 1 = 1 and Using dimensional regularization, we find the following formula for the X 1 function: The above equation can be straightforwardly generalized in the generic complex space by replacing n → iω n → z 1 , n + m → iω n + iω m → z 2 . Eq. (S61) defines the one-photon vertex renormalization factor: Obviously X 1 (z 1 , z 2 ) = X 1 (z 2 , z 1 ).

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 the fermionic momentum, p. Moreover, for the optical limit we can neglect the photon momentum q. The self-consistent Bethe-Salpeter relation for the two-photon vertex function is depicted in Fig. 1c of the main text and it readŝ Λ 2 (n, n + m, n + 2m) =Λ (0) 2 (n, n + m, n + 2m) + γ imp kĜ (k, n)Λ 2 (n, n + m, n + 2m)Ĝ(k, n + 2m) . (S65) In the non-interacting Dirac system the "bare" two-photon vertex function is zero,Λ 2 ∝ δ 2Ĝ−1 0 /δA 2 | A→0 = 0, due to the linear momentum dependence of the Hamiltonian. However, due to interaction the unrenormalized two-photon vertexΛ (0) 2 is finite given by the following relation (see Fig. 1e of the main text) Λ (0) 2 (n, n + m, n + 2m) = −γ imp kĜ (k, n)Λ y (n, n + m)Ĝ(k, n + m)Λ y (n + m, n + 2m)Ĝ(k, n + 2m) . (S66) From now on we adopt the short-hand notation z j = n + jm with j = 0, 1, 2, 3. We findΛ in which Q 1 (z i , z j ) is the one-photon renormalization factor defined in the previous subsection, and where By performing the momentum integration using the dimensional regularization, we obtain In a compact form we can write By solving the the self-consistent Bethe-Salpeter relation for the two-photon vertex given in Eq. (S65), we obtain in which Q 2 (z 0 , z 2 ) is the two-photon Bethe-Salpeter renormalization factor and where We explicitly obtain Note that X 2 (z, z ′ ) = X 2 (z ′ , z). Using the self-energy relation Eq. (S10) we have Therefore, the two-photon renormalizationfactor reads (Eq. (6) of the main text) It is useful to evaluate of Q RR 2 ( , ) which is given in Eq. (7) of the main text. Using the above relation in the retarded-retarded (RR) channel and employing the self-energy relation Eq. (S10), we obtain Therefore, we have Consequently, we arrive at Using Eq. (S10), we obtain Therefore, we find the following relation which is given in Eq. (7) of the main text.

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 in Fig. 1f of the main text. Accordingly we findΛ Here Q 1 (z i , z j ), Q 2 (z i , z j ) are the Bethe-Salpeter one-and two-photon renormalization functions, respectively. The explicit expression for Ω 1 function is given by Similarly, one can obtain Finally, the Bethe-Salpeter renormalization of the three-photon vertex function gives: where Note that X 3 (z 1 , z 2 ) = X 1 (z 1 , z 2 ).

S4. LINEAR CONDUCTIVITY
Linear response function is obtained after performing a Matsubara summation as follows Note that β = 1/k B T where k B is the Boltzmann constant and T stands for the electronic temperature. The Pfunction is analytical in the complex plain except two branch cuts at and − m where span over whole real axes. After performing the summation and an analytical continuation as m → ̵ hω + i0 + , we find [49] where n F (x) = 1/(e βx + 1) is the Fermi-Dirac distribution function. Note that "R" and "A" superscripts stand for the retarded and advanced, respectively. Accordingly, we have P RR ( , + ̵ hω) = [P AA ( , + ̵ hω)] * = P ( +i0 + , + ̵ hω+i0 + ) and P AR ( , + ̵ hω) = P ( − i0 + , + ̵ hω + i0 + ). The linear optical conductivity reads For the 2D Dirac model, the P -function reads 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 = 2 and N v = 2. The linear dc-conductivity follows where σ 0 = e 2 /4 ̵ h. The retarded self-energy can be decomposed into its real and imaginary parts Σ( ) = ∆( ) − iΓ( ) with ∆( ) and Γ( ) > 0 being odd and even real functions, respectively. Using this notation we find X RR 1 (0, 0) = −1, X AR 1 (0, 0) = w(x) with x = µ/Γ(µ) in which the renormalized chemical potential is given by µ = µ 0 − ∆(µ), and w(x) reads (S99) Therefore, we find ; U (S100) where . (S101) The functional dependence of f 1 (x; U ) on x and U is illustrated in Fig. S2. have (S102) The asymptotic form of f 1 (x) for small and large x follows

S5. NONLINEAR DC CONDUCTIVITY
A. Analytical derivation f 3 (x) function in the constant-Γ model Nonlinear dc conductivity is given by where the imaginary part of nonlinear (third-harmonic) response function follows

S6. ANALYTICAL CONTINUATION FOR THE THIRD ORDER RESPONSE FUNCTION
The summation of the fermionic Matsubara frequency n is performed by the contour integration technique.

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 σ  In Fig. S5, we illustrate the universal scaling of f 3 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 which support 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 completely different from that of the full quantum theory which is given in Fig. 3c of the main. text. We can see only one sign-change in the constant-Γ model in contrast to that of full quantum theory which gives two sing-changes. Unlike the full quantum theory, the constant-Γ model predicts a positive nonlinear correction in the quantum regime.  dc is written on the plot where the sign-switch border is highlighted by a dashed red line. Green and blue dotted lines stand for the contour lines with g = 1 and g = 0.1, respectively.