{'fraction_non_alphanumeric': 0.038901746932563964, 'fraction_numerical': 0.008306648489267937, 'mean_word_length': 4.474140923290524, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Photogalvanic effect in Weyl semimetals

We theoretically study the impact of impurities on the photogalvanic effect (PGE) in Weyl semimetals with weakly tilted Weyl cones. Our calculations are based on a two-nodes model with an inversion symmetry breaking offset and we employ a kinetic equation approach in which both optical transitions as well as particle-hole excitations near the Fermi energy can be taken into account. We focus on the parameter regime with a single photoactive node and control the calculation in small impurity concentration. Internode scattering is treated generically and therefore our results allow to continuously interpolate between the cases of short range and long range impurities. We find that the time evolution of the circular PGE may be nonmonotonic for intermediate internode scattering. Furthermore, we show that the tilt vector introduces three additional linearly independent components to the steady state photocurrent. Amongst them, the photocurrent in direction of the tilt takes a particular role inasmuch it requires elastic internode scattering or inelastic intranode scattering to be relaxed. It may therefore be dominant. The tilt also generates skew scattering which leads to a current component perpendicular to both the incident light and the tilt. We extensively discuss our findings and comment on the possible experimental implications.


I. INTRODUCTION
In the present days, Weyl semimetals [1] enjoy significant scientific interest which is generated by a multitude of reasons. From the theoretical viewpoint, this class of materials is special since it represents a solid state realization of Weyl's theory of chiral relativistic fermions [2] -an elegant theory which regrettably seems to have lost its connection to fundamental particle physics after the discovery of neutrino oscillations [3]. In the context of transport theories, the unique features related to the chiral anomaly [4] are at the basis [5] of the experimentally observed giant magnetoresistance [6]. From the viewpoint of applications, Weyl semimetals are attractive as their topologically protected band touching associated with a Berry curvature singularity may allow to explore and exploit novel regimes and phenomena of semiconductor physics. As we understand now, the appearance of Weyl nodes in systems lacking inversion or time reversal symmetry is far from being exceptional and the observation of Weyl physics has been already reported in various materials [7]. One particularly exciting phenomenon occurring in Weyl semimetals is the photogalvanic effect (PGE) [8][9][10][11], i.e. the generation of current due to the exposure to light. In contrast to most ordinary semiconductors, Weyl semimetals are susceptible to lowest frequencies and allow for novel technologies in the infrared. This consequence of the protected gapless spectrum comes along with the theoretical prediction of a quantized circular PGE [12] which is awaiting its experimental verification.
The theoretical foundations of PGE go back to the seventies (for review see Refs. [13][14][15]) while modern derivations rely on the nonlinear susceptibility framework [16], Floquet theory [17] and the Keldysh quantum kinetic equation approach [18]. In this paper, we concentrate on the dc response. It is common to distinguish two contributions: the injection current and the shift current. The former of the two represents the current generated by photoelectrons which are excited by optical (vertical) transitions with an anisotropic transition rate. In the steady state, this contribution diverges linearly when the current relaxation time becomes infinite. Differently stated, in the absence of a relaxation mechanism, the injection current grows linearly in time. In contrast, the shift current [19,20] is due to the real space displacement that electrons undergo upon an optical transition. This contribution is always finite and thus subleading for long relaxation times.
Therefore, in the present paper we keep only the dc injection current. We focus on the parameter regime with a single photoactive node which materializes the predicted topological quantization [12]. We investigate its relaxation due to impurities taking into account both intra and inter valley scattering in a model of two Weyl nodes. We scrutinize the time evolution of the current in response to a sudden illumination as well as its steady state characteristics. Furthermore we study the effect of a weak tilt in the Weyl spectrum. As we show, the interplay of impurity scattering and tilt leads to skew scattering and to the appearance of current components transversal to the direction of illumination.
It is also worthwhile to emphasize, that the current injection mechanism for the nonlinear PGE response at large frequencies is strongly different from the nonlinear intraband response at small frequencies considered in Refs. [21][22][23]. We explicitly show that the latter exactly vanishes in our model for the case of absent tilt and intervalley scattering and that it is generically subdominant for frequencies exceeding the elastic scattering rate.
This article is structured as follows. In Sec. II we introduce the model under consideration and the assumptions for the calculation. In Sec. III we present a sketch of the technical derivation as well as the main results, which are subsequently discussed in Sec. IV. We conclude with a summary and an outlook about the experimental relevance of our findings, Sec. V. The Appendix contains details on the model (Appendix A), the elastic scattering rates (Appendix B), the relaxation of the photocarriers (Appendix C), the intraband response (Appendix D), and a microscopic tight binding model (Appendix E).

II. MODEL AND ASSUMPTION
The terms entering the free Hamiltonian H 0 (p, x ) = H kin (p) + H dis (x ) are as follows. The kinetic term is a generalized low energy version of lattice models presented in Refs. [12,24] and Appendix E (1a) Here, Pauli matrices σ = (σ 1 , σ 2 , σ 3 ) act in spin space, while Pauli matrices κ = (κ 1 , κ 2 , κ 3 ) act in node space. The terms proportional to velocities v, u with 0 ≤ u v represent Weyl dispersion and the tilt in directiont respectively, the energy scale γ is the offset between the two nodes. The spectrum ξζ (p) = ζ(ξv|p − ζb| + ut · [p − ζb] + γ) of the kinetic Hamiltonian is characterized by quantum numbers ζ = ±1 (node), ξ = ±1 (band) and p (momentum) and plotted in Fig. 1. This kinetic term is supplemented by a disorder potential where V (x ) = Rj V(x − R j ) is a sum over impurity potentials V(x ) which are centered at uniformly distributed positions R j in R 3 . Finally there is a homogeneous electric field eE (t) = e ± E ± e ±iωt = −∇Φ(x , t) enclosed in our model Hamiltonian by means of H = H 0 (p, x ) + Φ(x , t) (e is the elementary charge).
By construction, the model defined by Eqs.
(1) is not invariant under time reversal, space inversion or rotation symmetries, see App. A for more details. We note that the tilt direction could also be assumed to be equal in both nodes [25]. In contrast, we chose a model with opposite tilt in opposite cones, so that the tilt preserves inversion symmetry.
Our calculation of the photocurrent in Weyl semimetals relies on the following simplifying assumptions. First, we neglected any spatial dependence of the electric field which is justified for v/c → 0, where c is the speed of light. In doing so, we omit the photon drag effect.
The assumption of a homogeneous electric field treatment within a 3D bulk Hamiltonian also relies on the second assumption of a long Thomas-Fermi screening length. In this limit, the contribution of surface states can be expected to be subdominant. Third, in neglecting all other sources of scattering we assume impurities to be the dominant source of photocurrent relaxation. While this statement is generally expected to be true at sufficiently low temperatures where electron-phonon and electron-electron scattering rates are small, we will explicitly show the limitations of this picture. Fourth, the current relaxation is calculated using a semiclassical kinetic equation approach which is justified if the mean free path of photoelectrons and photoholes exceeds their wavelength. This implies that the frequency of light is much larger than the elastic scattering rate. Formally, we implement this assumption by keeping only the leading O(n −1 imp ) terms in small impurity density. For the evaluation of quantum transition probabilities we truncate the T-matrix at next to leading (i.e. second) order in powers of weak impurity potentials. We also assume that the impurity potential is short ranged as compared to the wavelength of photocarriers (but it may be longranged as compared to 1/|b|). Fifth, our calculation at the lowest temperatures is exponentially accurate if the energy difference between Fermi energy and the energy of photocarriers exceeds the temperature. Finally, all presented formulae are valid to linear order in u/v 1, unless stated otherwise.

III. SKETCH OF THE CALCULATION AND RESULTS
In this section we present a sketch of the derivation of the photocurrent for a tilted, disordered Weyl semimetal. For the sake of a clearer presentation, we here restrict ourselves to the case of photocarrier generation at the ζ = +1 cone only (see Fig. 1) and of absent intervalley scattering and relegate details and the more general case of finite intervalley scattering to Appendices B and C.
Our calculation is based on the Boltzmann kinetic equation [26,27] describing the time evolution of the distribution function f (p, t) (2) Here,ṗ = eE (t) and the collision integral contains two contributions. First, there is a term describing the excitation rates of photocarriers Here, we assumed states with energy ξζ (p) [ −ξζ (p)] to be empty [filled]. In the considered parameter range, ξ = 1 (ξ = −1) for photoelectrons (photoholes) and ζ = 1, see Fig. 1. Clearly, terms proportional to the Berry curvature Ω ξ = −ξp/2p 2 change sign when the chirality of the photoactive node is reversed. Note that St inj contains two anisotropic terms proportional to the direct productp ⊗p and to Ω ξ as well as an isotropic term proportional to 1. The second part of the collision integral describes the relaxation of photocarrieres by impurities We introduced the shorthand notation p = d 3 p/(2π) 3 . Since impurity scattering is fully elastic, it may only relax anisotropic terms in the photoexcitation rates. The intervalley scattering discussed below and in Appendix C may also relax the valley imbalance of the isotropic contribution to the photoexcitation rates.
In the limit u/v = 0, time reversal and rotational symmetries dictate that the scattering probability w p p contains only terms proportional to 1 and p · p and thus w p p = w pp . This statement is fulfilled also when intervalley scattering is present, see Appendix B. However, the symmetry of the scattering probability is lost in the presence of a finite tilt u/v = 0. Here we present the dominant contributions of symmetric and antisymmetric scattering probability: Here, V p is the Fourier transform of the impurity potential and the density of states (DOS) is κ  Fig. 1 we assumed ω = 3γ for this plot and furthermore we assumed Γ sk,+ = Γ+/3.
The skew scattering probability obtained in Eq. (4b) relies on the third moment of the distribution function of the disorder potential. We note that, by means of both semiclassical [28] and diffractive [29,30] mechanisms, skew scattering may also occur for Gaussian disorder models. However, skew scattering probabilities due to Gaussian disorder contain an additional factor of n imp and are thus subleading in the perturbation scheme employed here.
To obtain the photocurrent in the steady state, the solution for the distribution function can be obtained by equating the full collision integral to zero. We obtain the following nonequilibrium corrections to the Fermi-Dirac distribution function in band ξ where we introduced We also introduced the intravalley scattering rates for momentum relaxation and skew scattering which both implicitly depend on the energy of the photocarriers by means of the density of states. In Eq. (6a) the scalar part of the distribution function was omitted (it will be recovered shortly). Furthermore, Eq. (6a) was derived using the premise that only O(u/v) corrections proportional to odd powers ofp enter the final expression for the photocurrent.
We insert this expression into the definition of the current Here the sum over ξ reflects the two types of photocarriers that the injection term generates: electrons in band ξ = +1 at energy e = γ + ω/2 and holes in band ξ = −1 at energy h = γ − ω/2. As a side remark, we note that anomalous velocity terms δv ∼ Ω ξ × E are unimportant for the photocarriers, since their nonequilibrium distribution function Eq. (6a) is quadratic in electric field.
The final steady state current can be written as We present the microscopic values of the effective scattering rates τ CPGE , τ skew , τt, τ ⊥ in Eqs. (9b)-(9i), below. In Appendix C, we generalize the calculation presented here to a finite intervalley scattering described by the param- The latter interpolates between longrange impurities (κ = 0, no intervalley scattering) and short range impurities (κ = 1, strong intervalley scattering). The relaxation times in the limit of weak intervalley scattering κ → 0 behave as In contrast, for κ → 1 we obtain General formulae which interpolate between these two limits are presented in Eq. (C23) of Appendix C and are plotted in Fig. 2. Note that all relaxation rates Γ ζ ( ) and Γ sk,ζ ( ) introduced in these formulae are implicitly energy dependent. However, in the photoactive cone, the particle-hole symmetry about the Weyl node implies for both photoelectrons and photoholes equal scattering rates Γ + ( e ) = Γ + ( h ), Γ sk,+ ( e ) = Γ sk,+ ( h ). The rate τt is not defined in a noninteracting model without intervalley scattering as it is determined by the relaxation time of the isotropic part in Eq. (3a). This can only be achieved by inelastic scattering or by means of intervalley scattering. We discuss this issue and other physical implication of the result presented in Eqs. (9), Eq. (C23) and Fig. 2 in Sec. IV.
In Appendix C we also present a derivation of the time evolution after a sudden illumination in the limiting case u/v = 0. Technically, this amounts to adding the solution of the homogenoeous Boltzmann equation (Liouville equation), Eq. (2), to the particular steady state solution, Eq. (6a). In terms of the final result Eq. (9a) this amounts to the following replacement τ CPGE → T (t) with where we assumed the light to be switched on at time t = 0. For a more general formula of T (t) which interpolates between the two limits κ = 0 and κ = 1, we refer the reader to Eq. (C26) of Appendix C as well as to Fig. 3. Again, the discussion of this result is relegated to Sec. IV.
Finally, the Appendix D contains a calculation of the intraband rectified response due to particle hole excitations near the Fermi surface. There we show, that this contribution vanishes in the absence of tilt and intervalley scattering and that it is generically suppressed by the parameter Γ + /ω 1.

IV. DISCUSSION
In this section we analyze and discuss the physical content of the major results of this work, Eqs. (9) and (10) as well as Figs. 2 and 3. Γ + T(t) Color online] Time evolution of the circular photogalvanic effect after a sudden illumination at time t = 0 for κ = 0, 1/3, 2/3, 1 from top to bottom. In this plot ω = 3γ was chosen, i.e. the case represented in Fig. 1.
We begin the discussion considering the limit of absent tilt u/v = 0. In this case, the dc photocurrent is proportional to −i(E + × E − ) and is thus crucially dependent on circularly polarized light. This is the injection mechanism: vertical transitions inject photocarriers moving predominantly in direction −i(E + × E − ), see Eq. (3a). This leads to the linearly increasing time-dependent response, Fig. 3 at small t 1/Γ + . The slope is topologically quantized [12] in the present case of a single photoactive Weyl node. Disorder exponentially relaxes the linear increase and a steady state is formed, see Fig. 3 at large t 1/Γ + . It is worthwhile to notice that 1/Γ + is the "transport" mean free time, which is a factor of 3 larger than the "quantum" mean free time responsible for the level broadening. While these features are generic to any disorder model, intervalley scattering has the following three additional effects. First, it increases the scattering rate by increasing the density of available final states. This implies faster relaxation and a smaller steady state value. Second, since the relaxation rates at energies e (photoelectrons) and h (photoholes) are different in the ζ = −1 node, intervalley scattering introduces two different relaxation rates for the circular PGE. Third, for intermediate κ and provided Γ − < Γ + intervalley scattering can lead to a peculiar nonmonotonic time dependence, see Eq. (C26) and the κ = 1/3 and κ = 2/3 curves in Fig. 3. For the parameters chosen in our plot the nonmonotonicity stems from the photoholes living at energy h . Those scatter from the photoactive ζ = +1 node to the ζ = −1 node where the current is opposite and the decay much slower in view of a smaller DOS.

B. Finite u/v
We now turn to the steady state solution at finite u/v. The finite tilt allows for the presence of additional contributions to the current, see Figs. 2 and 4.
First, skew scattering leads to a term −i(E + × E − ) ×t. Just as the circular PGE, it is only present for circular polarization of light. As we already mentioned, the momentum of excited photo electrons predominantly points in direction −i(E + × E − ). The time-reversal symmetry breaking tilt introduces a finite anisotropy leading to a preferred direction of momentum relaxation. The resultant imbalance between electrons moving towards −i(E + ×E − )×t as compared the opposing direction generates the skew scattering term proportional to τ skew . It is in accordance with intuitive expectation that opposite tilts in opposite nodes imply antagonistic skew scattering contributions from opposite Weyl nodes. It is however surprising to find, that the O(u/v) skew scattering contribution exactly vanishes at strongest intervalley scattering κ = 1, see the blue dot-dashed curve in Fig. 2.
Apart from skew scattering, the finite tilt in the spectrum also introduces terms which do not rely on circular polarization of light and which may be present also for linearly polarized light. We first focus on the term proportional tot(E + · E − ). This term stems from the (u/v)t term in the definition of the current, Eq. (8). Therefore, it stems from the scalar (isotropic) part of the distribution function of photocarriers, or differently said, from the scalar term proportional to 1 in Eq. (3a). As we already mentioned, without intervalley scattering impurities cannot relax such a term. This is why the photocurrent in directiont, represented by a red arrow in Fig. 4, is much stronger than all other contributions. The presence of intervalley scattering κ > 0 introduces a finite τt, which however diverges as 1/κ for small κ, see Eq. (9d). This is represented in the red dashed curve in Fig. 2 which, for the sake of a better presentation, had to be downscaled by an extra factor of 20. In the chosen parameter regime, τt is particularly large, because Γ − is relatively small at h . Even though inelastic scattering is beyond the scope of the present paper, we mention that in practice, the steady state value of the current in t direction is determined by the smaller of τ inel (inelastic scattering time) and τt.
Finally, the steady state result, Eq. (9a) contains a term proportional to Re[E + ⊗ E − ]t. This term stems from the tensor contribution to the photocarrier excitation rate, i.e. thep ⊗p term in Eq. (3a) and is absent in the absence of intervalley scattering.

V. SUMMARY AND OUTLOOK
In summary, we have derived and analyzed the disorder induced relaxation of the photogalvanic effect (PGE) in a simple model for a Weyl semimetal. We took into account intra and intervalley scattering as well as leading order  Fig. 2). Here, γ is the direction of incident light. The orange −i(E+ × E−) term is the "quantized" response for circularly polarized light [12]. The blue vector −i(E+ × E−) ×t is directed into the plane spanned by the directions of light propagation and tilt t. Just as the green Re[E+ ⊗ E−]t and redt(E+ · E−) vectors it is a consequence of finite tilt. In the absence of intervalley scattering, the red contribution proportional tot can not be relaxed by impurities and may therefore be dominant, see The observation of decisive qualitative and quantitative consequences of a weak spectral tilt corroborates the findings of earlier studies on different observables such as the linear conductivity tensor [31][32][33], the polarization function [25] as well as transport through Weyl tunnel junctions [34]. Further implications on characteristic features of Weyl fermions, e.g. the natural optical activity [35], will be the subject of a separate studies [36].
We conclude our paper with an outlook on the experimental relevance of our findings. First of all, we comment on the minimal two Weyl node model that we are considering. Such materials are not discovered so far, but they were suggested theoretically in heterostructures [37] and alloys [38] of Chern insulator materials as well as, most recently, in certain magnetic Heusler compounds [39]. At the same time our simplified model may be applied to present day materials when pairs of Weyl nodes are well separated in momentum space.
Concerning optical experiments, we are aware of only two published results on the PGE in Weyl semimetals. The work of Ref. [9] is devoted to the second harmonic generation, i.e. a different observable as compared to the focus of this paper. In the second experiment [11] the chirality of the Weyl fermions, which is given by the sign of the Berry monopole charge, was inferred from the photocurrent response to mid-infrared light. Two key points are important: (i) Experiment reveals a current component perpendicular to the direction of light which strongly depends on the polarization (linear, circular positive, circular negative) and which is claimed to be proportional to the chirality of the photoactive nodes. As we show, skew scattering also induces a term of the very same tensor structure. (ii) Polarization independent photocurrent is observed in a direction perpendicular to both light and to the polarization dependent current. As we show, such terms may also be induced by the tilt. In addition, unpublished THz spectroscopy data on the dc PGE in TaAs presented at the APS march meeting 2017 [40] indicates that the signals of radiated electric field in direction perpendicular or parallel to the polar axis show very different behavior and that for the latter case the signal is nearly independent of the polarization of incident light and is relatively long lived. The preliminary interpretation focused on anisotropic scattering and is thus related to the skew scattering mechanism investigated here. Another route for interpretation could be contributions similar to the termst(E + · E − ) considered in this paper. We repeat that these are indeed long lived and independent on the polarization of incident light. At the same time, we note that the tilt of the Weyl fermions in TaAs is substantial (v/u > 2) and, therefore, our perturbative calculations should not be expected to quantitatively describe those experiments since our time-reversal breaking two node toy model has limited resemblance with the time reversal conserving 24 Weyl node material TaAs. However, we are convinced that the present results provide a proof of principle for various types of photocurrent contributions which will trigger future investigations.

Deformed coordinate system
In view of the tilt, the Fermi surfaces are ellipsoidal. Here we introduce a coordinate transformation which deforms the ellipse into a circle at the expense of certain Jacobians.
We first use that at a given energy the equienergy surface of a band characterized by (ξ, ζ) is an ellipse in the p t − p ⊥ plane (p t =t · p, p ⊥ = p 2 − p 2 t ) with center at (p t , It is useful to choose the following energy dependent parametrization of the three vector p .
The integration measure transforms as and, since a delta function in energy space can be written as we can write Note that, when k is a unit vector, we can writê We will use the following parametrization of f ξ,ζ (p) where the transformation p → ( ,k) depends on the cone ξ, ζ under consideration. We repeat that the band index ξ is uniquely determined by and ζ.

Appendix B: Elastic scattering rates
The scattering probability for scattering from the state l to the state l is where are the matrix elements of the retarded transition matrix. We next present restrictions on the scattering probability imposed by the symmetries discussed in the main text.

Symmetry considerations
As we mentioned in the main text, our model is not invariant under any of time reversal, rotation or inversion symmetries. We find We used the representationsT = σ y K (K is the complex conjugation);Î = κ 1 I (I : The symmetry operations on the wave function are as followsT Here, again we usedp 1 + ip 2 = (1 −p 2 3 )e iφ . For the rotational symmetry it is more useful to investigate its action on a projector onto an eigenstatê We exploit those transformations in the analysis of the scattering probabilities in the limiting case u/v = 0. For the present model, rotational symmetry plays a major role and we find that, for u = 0 and after average over disorder configuration, which implies that w l l is a scalar or pseudo scalar under transformations and can be expanded by the basis functions 1,p ·p ,p · b,p · b, (p ×p ) · b. We can use the fact that after disorder average, the average translational inveriance implies that b enters only in the form of |V b | 2 under the assumption of sufficiently short ranged impurities. Therefore, rotational invariance implies that w l l = [f ζ,ξ;ζ ,ξ + g ζ,ξ;ζ ,ξ p ·p ]δ( l − l ).
In particular, w l l is symmetric under the exchange of momenta.
We can further use the implication of TR symmetry and find so that and thus the absence of skew scattering: w l l = w ll in the limit u = 0.

Scattering probabilities a. Born approximation
The scattering probability in Born approximation can be readily evaluated from Eq. (B1) For the skew scattering we expand T l l l |Ĥ dis +Ĥ dis G R 0Ĥdis |l and obtain [41] w (3a) l l For our specific model this leads to We use that This concludes the derivation of Eq. (4b) of the main text.

Scattering rates
It will be useful to expand the distribution function for the band (ξ, ζ) by means of representations of the rotation group. We write Furthermore, the kinetic equation and the collision integral shall be written as a vector f = (f ζ=+1 , f ζ=−1 ) T in valley space and denote it by an underbar. Note that, at a given energy there are two bands involved, one at each valley.
We can write the collision integral entering the equation for electrons in pocket ξ, ζ The notation St[{f }] reflects that in general the collision integral depends on the whole set of distribution functions, i.e. in the present case on the distribution functions of both nodes. All energies are taken at given ξ,ζ which is either e or h . We employ the valley space vector notation for the scalar part of the distribution function to O(u/v) Here, the matrix ξ = diag(ξ + , ξ − ). As we shall see later, for the vector part we don't need u/v corrections which contain an even power of vectorsk. Under this premise (reflected by the sign . =), only the nontilted contribution survives Next, we turn our attention to the matrix distribution function. Here we obtain The sum over repeated indices is to be understood. We will find later that the solution for the tensorial part of the distribution function f ijk is first order in u/v. Therefore, we only need the zeroth order in u/v from the collision integral, i.e.
For the skew scattering the leading order is (u/v) 1 and we only need to consider the contributions of the scalar part (vanishes), of the vector part (nonvanishing) and of the matrix part (vanishes) of the distribution function.
After angular integration the skew-scattering collision integral for electrons in band (ξ, ζ) becomes We define the band dependent skew scattering rate Γ (sk) ζ = π 2 3 n imp V 3 0 ν 2 ξζ . We further use that Γ + /Γ − = ν + /ν − to reexpress the ratios of DOSs. Then Note that Γ (sk) Appendix C: Relaxation of photocarriers We use the notation ζ 0 = sgn( F ) for photoactive node and we further write ξ 0 for the band index in that node at a given energy e or h .

Injection term
The injection term, i.e. the collision term in the Boltzmann equation generating photoelectrons/photoholes in band (ξ 0 , ζ 0 ) in vector notation is given by the Fermi's golden-rule expression For the product of transition matrix elements we calculate symmetric (round brackets) and antisymmetric (square brackets) components separately This leads to Eq. (3a) of the main text. After the change of coordinates fromp tok the injection term in vector notation can be written as where we introduced (recall E(t) = ± E ± e ±iωt ) as well as As we mentioned previously, we are only interested in the dc part of these expressions since ω/Γ ζ0 1 implies that second harmonic generation is only weakly affected by disorder.

Static solution to O[(u/v) 0 ]
We expand f = f (0) + (u/v)f (1) and keep only the zeroth order in (u/v) in this section. Upon equating the full collision integral St inj +St (2) +St (3a) to zero we obtain The solutions for matrix and vector part immediately followk We insert these solutions back into the Eq. (C6) and then solve for the scalar part The symbol " . =" means that this equation is valid only if f (0) is multiplied by (1, −1) from the left. This is, because the scalar symmetric part (f ζ=+ + f ζ=− ) corresponds to a zero mode of the collision integral: The inverse of Γ −1 s means (symbolically in the subspace of the antisymmetric part (f ζ=+ − f ζ=− )) As we shall show below, for the calculation of the current we only need terms containing an odd number ofks in the linear corrections to the distribution function. We denote this assumption by a " . =" in the following equation for the first order terms After some algebra, this leads to This solution is inserted back into Eq. (C10) and the solution leads to

Static current response
The current is defined as Since we found that distribution function is proportional to δ(ω − 2vp) it is useful to switch to spherical coordinates in p space (i.e. we transform back from thê k representation). From this equation we readlily see that in the terms of first order in (u/v) of f (p) we need only terms which contain odd powers ofp.
We will write Here we introduced as well as We used that both Γ t and ξ are diagonal matrices.

Dynamic solution to O[(u/v) 0 ]
In order to obtain the time dependence of the photocurrent after a sudden illumination, we add the solution of the homogeneous Boltzmann equation to our steady state solution, Eq. (C7b). Under the assumption of thermal equilibrium at t = 0 we obtain.