Intrinsic Fermi Surface Contribution to the Bulk Photovoltaic Effect

We study the Fermi surface contribution to the nonlinear DC photocurrent at quadratic order in a spatially uniform optical field in the ultra-clean limit. In addition to shift and injection current, we find that polarized light incident on a metallic system generates an intrinsic contribution to the bulk photovoltaic effect deriving from photoinduced electronic transitions on the Fermi surface. In {velocity} gauge, this contribution originates in both the coherent band off-diagonal and diagonal parts of the density matrix, describing respectively, the coherent wave function evolution and the carrier dynamics of an excited population. We derive a formula for the intrinsic Fermi surface contribution for a time-reversal invariant chiral Weyl semimetal illuminated with circularly-polarized light. At low frequency, this response is proportional to the frequency of the driving field, with its sign determined by the topological charge of the Weyl nodes and with its magnitude being comparable to the recently discovered quantized circular photogalvanic effect. Our work presents a complete derivation for all contributions to nonlinear DC photocurrent and classify them according to the polarization of light in the presence and absence of TRS.

Introduction-Interest in developing new platforms for efficient solar energy conversion has drawn attention to the photovoltaic properties of new materials and the physics of their light-matter interactions. The bulk photovoltaic effect (BPVE), sometimes also referred to as the photogalvanic effect (PGE), has attracted much attention, as it can directly convert light to a DC current [1]. The BPVE is a second-order nonlinear response that can be decomposed into terms that are symmetric and antisymmetric in the polarization states of the light, corresponding to a linear photogalvanic effect (LPGE) and a circular photogalvanic effect (CPGE), respectively [1]. The ballistic current is an important mechanism of the BPVE and emerges in both LPGE and CPGE [2][3][4][5], where an asymmetric distribution of photo-excited charge carriers on the conduction band are induced by electron-phonon or electron-electron scattering. Apart from the ballistic current, the CPGE is usually described as an injection current. Owing to the phase lag between orthogonal components of a circularly-polarized beam, an asymmetry in the excited state population at time-reversed momenta k and −k can be induced [6,7]. In a two-band model, the carrier generation rate can be related to the Berry curvature, which relates the trace of the CPGE tensor to the quantized topological charge of degenerate points in the band structure for a Weyl semimetal [8]. In LPGE, other than ballistic current, shift current can be described as a coherent response associated with the real-space shift of an electron induced by a dipole-mediated vertical interband transition [9]. Recent photo-Hall measurements of both LPGE and CPGE in an applied magnetic field have successfully separated the shift and ballistic contributions to the electric charge current [10,11]. First principles studies [12][13][14] of the shift current have led to the prediction and discovery of many new photovoltaic materials [15][16][17][18][19][20][21].
Light-matter interactions in these materials are described by coupling electrons in the material to the electromagnetic potentials A µ . In velocity gauge, the external electric potential is taken to vanish, and the external electromagnetic vector potential can be incorporated into the electronic Hamiltonian through a minimal coupling procedure that augments the electronic momentum operatorp →p + eA(r, t)/ [22,23]. For a spatially uniform, but time varying electric field and a vanishing magnetic field, a time-dependent gauge transformation on the electronic wave function can bring the effective Hamiltonian back to its original unperturbed form with the addition of an electric-field-induced perturbation δH = eE(t) ·r [24,25]. This form of the Hamiltonian is called length gauge and has been used to derive expressions for the contributions to nonlinear electric currents in insulators and semimetals [6,26,27].
Using time-dependent perturbation theory in velocity gauge, we derive formulas for the nonlinear photocurrents induced at quadratic order in an external driving field. We find a geometric contribution to the secondorder current related to the curvature of the electronic Bloch bands that must vanish for insulating materials, but can be nonzero for metals and semimetals. This contribution to the current derives from both off-diagonal and time-independent band diagonal contributions in the density matrix, and it is proportional to k-space derivatives of the unperturbed Fermi occupation factors multiplied by the band-resolved Berry curvature or quantum metric tensor [28,29], depending on the polarization of the light. For a time-reversal (TR) invariant system, the current is only non-vanishing under circularly-polarized (CP) illumination. By breaking time-reversal symmetry (TRS), both circularly-and linearly-polarized light can induce a current. At zero temperature, the current derives from dipole-allowed vertical inter-band transitions arXiv:2011.06542v3 [cond-mat.mtrl-sci] 4 Dec 2021 for electrons whose crystal momenta are near the Fermi surface. Thus we name this new photocurrent the "intrinsic Fermi surface contribution", defining a new type of CPGE and LPGE response.
As an example, we calculate this contribution to the CPGE in a three-dimensional time-reversal symmetric, but inversion symmetry broken Weyl semimetal. We show that for an isotropic Weyl node, this response is proportional to the Weyl node's charge weighted by a factor that depends on the energy of the Bloch bands near the Fermi surface. In the small ω limit, the photovoltaic response for an isolated Weyl node is purely quantized and would therefore directly measure the topological charge of the Weyl point. When mirror symmetries are broken in a chiral Weyl semimetal, point nodes with opposite topological charges are offset in energy, allowing a nonzero DC charge current to flow. This can be distinguished from the extrinsic quantized circular photogalvanic effect deriving from the injection contribution to the current, which is proportional to a scattering time τ [8]. General theory-We study the Hamiltonian for a particle with charge −e coupled to a time-dependent vector potential A(t) H(r,p, t) = (p + eA(t)) 2 2m e − e V (r) (1) where V (r) is the crystal potential, m e is the mass of the electron, andp andr are the electronic momentum and position operators. We work in velocity gauge, where the response to the electric field is made by the identification Here, we focus on the single-particle Hamiltonian,Ĥ(r,p, t) and do not include the contributions deriving from many-body interactions such as excitonic effects in our model. The nonlinear DC charge current is determined by calculating the trace of the product of the velocity operatorv(t) = i[Ĥ(r,p, t),r]/ and density matrixρ(t). To isolate the nonlinear response, we expand bothv(t) andρ(t) in a power series of A(t) up to quadratic order in the external driving field (denoted by the superscripts on the operators below) and calculate the current as (See Appendix for details). To solve forρ(t), we employ the von Neumann equation which describes the time evolution of this quantum operator [30]: The HamiltonianĤ(r,p, t) can be divided into two pieces.Ĥ 0 (r,p) = p 2 2me +V (r) describes the unperturbed Hamiltonian before application of the external field and H (r,p, t) = eA(t)·p me + e 2 A(t)·A(t) 2me describes the interaction between electrons in the material and the external field. With this substitution, we can solve equation 3 for the density matrix order by order in the vector potential A(t). Here, for simplicity Eq. 1 starts from a nonrelativistic quadratic Hamiltonian. In the Appendix, we show an equivalent derivation for a general Bloch Hamil-tonianĤ(k).
We focus on external driving fields with few nonzero Fourier components and write A(t) = ω A(ω )e i(ω −iη)t with A(ω ) = A * (−ω ). The limit as η → 0 denotes an adiabatic turning on of the electromagnetic field. At second order in the electromagnetic vector potential, the current couples to two different electromagnetic fields with unique Fourier components at two general frequencies ω 1 and ω 2 . The DC limit is found by expanding the current in powers of ω 1 + ω 2 to extract divergent and finite contributions to the current as ω 1 → −ω 2 and η → 0. Fermi surface contribution-We parse the nonlinear DC photocurrent into its diagonal j dia ∼ ev nn (k)ρ nn (k) and off-diagonal j off ∼ n =m ev nm (k)ρ mn (k) parts. Here O nm (k) = Ψ n (k)|Ô |Ψ m (k) , where |Ψ n (k) are the Bloch eigenstates of the unperturbed Hamiltonian H 0 (r,p). j dia describes the current generated from the dynamics of excited carrier populations: terms in j dia are proportional to both the velocity, v nn (k), and population density, ρ nn (k), of Bloch electrons in bands n with crystal momentum k. As shown by Eqs. (S20)-(S23) in the Appendix, j dia can be broken into three pieces but belonging to two types: (I) the divergent terms j dia1 proportional to 1/η e 2ηt that diverge as η → 0 and (II) j dia3 terms that are finite as η → 0. However, the generation rate ∂ t j dia1 remains finite. This is called injection current j inj , induced by resonant excitations between Bloch electrons in different bands with the same crystal momentum whose energy differs by ω.
The non-divergent diagonal contribution, j dia3 , can be written as Here, µ is the electron chemical potential. We note that j dia3 is rarely discussed in previous studies, but is an intrinsic contribution to the induced photocurrent: different from j dia1 and j dia2 , it is not proportional to 1/η. The contribution j off describes a current generated from the coherence between electronic Bloch states in different bands. At second order in the perturbing field, a third band state l is inserted as an intermediate transition state between the coherent pair n and m. As shown in appendix (Eqs. (S24)(S25)), the three-band processes n − → l − → m can be converted to an effective n ↔ m inter-band transition with an energy factor 1 εn(k)−εm(k)+ ω −i η . This can be decomposed into a resonant j off1 with δ function and an off-resonant j off2 with principal part, representing two types of couplings. j off1 is associated with the shift current j shift , where the electron coordinate in real space is shifted along with the resonant excitation from band n to m. For j off2 , after removing terms that vanish identically, it can be written as: We note that for insulators with a minimum band gap energy, E gap , terms contributing to both j off2 and j dia3 involve field-induced dipole mediated transitions between Bloch states that lead to nonzero currents even for light with frequency | ω| < E gap . Though j off2 and j dia3 derive from different lightinduced transition processes between Bloch states, the sum of the two contributions to the current can be simplified into a single term proportional to k-space derivatives of the unperturbed Fermi occupation factors f T n (k, µ): At zero temperature, k-space derivatives of the Fermi occupation functions are proportional to a delta function that is nonzero only for crystal momenta along the Fermi surface. This important simplification clarifies the absence of DC current generation for insulators perturbed by light with frequency | ω| < E gap . In experiment, only the total current is measured, and these sub-gap currents arising from j off2 and j dia3 cancel, leaving only contributions to the current for materials with non-vanishing Fermi surfaces. Therefore, this contribution is nonzero only for metallic systems; we denote this contribution to the current in equation 6 as j Fermi . At finite temperature, the sharp Fermi surface will smear out, allowing states with nonzero Fermi occupations above and below the Fermi surface to contribute to generation of current. Interestingly, one recent study discovered that the subbandgap photocurrents can be used to probe the magnitude of finite lifetimes [31], in contrast to the "clean limit" and independence of τ here. For a general Bloch HamiltonianĤ 0 (k), in addition to j Fermi there is another Fermi surface contribution to the current: Derivation of this contribution is given in the general Bloch Hamiltonian part in Appendix. As shown in appendix, the sum of two Fermi surface contributions (Eq. 6 and Eq.7) is equivalent to the sum of the Drude contribution, Berry curvature dipole contribution and the free-carrier contribution derived in length gauge. The Drude and Berry curvature dipole contributions can be decomposed from j Fermi2 but they only represent partial metallic contributions [32,33]. We note j Fermi2 vanishes if the energy dispersion relation ofĤ 0 (k) is linear or quadratic in the crystal momentum k. By lifting TRS, v nm (k) and v nm (−k) are no longer negative complex conjugates. BPVE phenomena are enriched in magnetic systems, as j inj , j shift and j Fermi are non-vanishing for driving fields with either LP and CP components. We summarize and classify our results accordingly in Table I. Example: single Weyl point-We calculate this new Fermi surface contribution j Fermi for a single isotropic Weyl node in a minimal two-band model: Here v 0 is the Fermi velocity, and σ α are the Pauli matrices. The system is illuminated by CP light described by the vector field A(t) = A 0 (cos(ωt)x + sin(ωt)ŷ) that propagates in theẑ-direction. With linearly dispersive bands, j Fermi2 is zero. At zero temperature, we write j Fermi in a different form: ) is the band-resolved Berry curvature. In this way, the DC charge current can be related with the integration of a weighted Berry curvature dipole over the Brillouin zone [27,34,35]. With A(t) introduced above, j Fermi only has aẑ component, wheren z (k) is the unit vector normal to the Fermi surface, Ω z (k) is theẑ-component of the Berry curvature, and µ is the chemical potential. For this simple model,  I. The classification of all types of contributions to the nonlinear DC photocurrent. "Intrinsic/extrinsic" means whether the current is independent/dependent on τ . "Diagonal/off-diagonal" refers to the part of ρ from which the current arises. The last four columns give the correspondence between each current and light polarization with and without TRS. Each contribution is gauge invariant, independent of the choice of phase of the Bloch wavefunctions. |ε n (k) − ε m (k)| = 2|µ| and is constant across the Fermi surface. For arbitrary CP light illumination, the current can be written as where χ ij (ω) is a purely imaginary photovoltaic tensor with the property where Q = 1 2π F S dS · Ω(k) is the charge of the Weyl point. In the limit of | ω| |µ|, the current is proportional to the Weyl node's charge, and the trace of χ(ω) is proportional to the Chern number of the Fermi surface which encloses this topological degeneracy. The charge is a purely topological aspect of the band structure and is robust even when the isotropic symmetry of the Weyl node is broken and the Fermi surface degenerates into an elliptic surface. Consistent with the injection current, the sign of the photocurrent is also dictated by the charge of the Weyl point, and can be used to detect the chirality of the band singularity [36]. We note that for a tilted Weyl cone the energy difference between bands is not constant along the Fermi surface and the Tr[χ ij (ω)] cannot be reduced to the simple form given in equation 11.
Unlike the quantized circular photogalvanic effect [8], here the current does not originate from coherence of optically coupled band states, but instead originates from electrons with crystal momentum along the Fermi surface. In addition, the effect is intrinsic and only depends on the light frequency ω and not on an extrinsic scattering time τ . This distinguishes itself from the injection current in that due to its coherent nature it does not require a finite τ to populate carriers [37]. For metallic systems illuminated by mid-infrared light, j Fermi is significant and comparable to j inj . Unlike the difference frequency generation scenario studied in other work [32], j Fermi can be excited by one monochromatic, polarized light. Example: multiple Weyl nodes-We now demonstrate calculation of j Fermi for a time-reversal symmetric, but inversion and mirror broken Weyl semimetal. For systems with TRS, Weyl points in the band structure must come in pairs. The two Weyl points in a pair have the same energy and topological charge, but are located at points in the Brillouin zone with opposite crystal momentum. The topology of the Brillouin zone demands that the sum of the charges of all Weyl points in the Brillouin zone must vanish: n Q n = 0 forcing the number of Weyl points in a time-reversal symmetric system to be a multiple of four. If mirror symmetry is broken, Weyl points of opposite sign need not occur at the same energy [8], allowing the energy differences between bands, ∆(k) = n (k) − m (k), near each Weyl point to be inequivalent (see Figure 1). We have seen that for a single Weyl point and for light illumination | ω| |µ|, the trace of χ(ω) is simply proportional to the integral of the Berry curvature across the Fermi surface. For a system of isotropic Weyl points illuminated by light with arbitrary ω, ∆(k) → ∆ and we may write the trace of χ(ω) as Here Q n is the charge of node n and ∆ n is the energy difference between Bloch states with crystal momentum along the Fermi surface near Weyl node n. With mirror symmetry, ∆ n are the same for differently-charged Weyl points. No charge current flows, but only a chiral current represents a charge pumping between Weyl nodes with different chiralities, which is not observable and is similar to the chiral anomaly [38][39][40]. If ∆ n are different, the trace of χ(ω) will be nonzero. Unlike the quantized CPGE induced by injection current where Pauli blocking can forbid the photocurrent when upper state is filled [8], the chemical potential can sit either above or below the Weyl node. Each node's contribution to the current does not change its sign whether the Fermi surface is an electron or hole pocket. In the limit | ω| |∆ n |, we can expand Tr[χ ij (ω)] in powers of 2 ω 2 /(∆ n ) 2 , We see that the leading-order term in the expansion is linear in ω, with slope determined by the ratio of the charges of Weyl nodes to the energy differences of the Bloch bands near the Weyl nodes. The relationships between ω and ωTr[χ(ω)], and ω and Tr[χ(ω)] for a minimal four-Weyl node system are plotted in Fig. 2. Breaking the isotropy of the Weyl nodes takes ∆ n → ∆ n (k), and the leading contribution to the trace of χ(ω) will no longer be directly proportional to the charges of the Weyl nodes. However, the first-order non-vanishing contribution to the current will maintain a linear relationship to the frequency of light.
Candidate materials-Noncentrosymmetric metals characterized by the coexistence of metallicity and ferroelectric distortions provide a promising platform for observing this novel addition to the nonlinear current. Experiments have demonstrated that metallic LiOsO 3 and Cd 2 Re 2 O 7 experience a centrosymmetric to noncentrosymmetric phase transition at 140 K and at 200 K, respectively [41,42]; while the engineering of interfaces in ANiO 3 /LaAlO 3 heterostructures provides another scheme for achieving other interesting noncentrosymmetric metals [43]. In addition, recent studies on few-layer topological semimetal WTe 2 have demonstrated a switchable ferroelectric polarization that could also provide a platform to observe large j Fermi under illumination by CP light [44]. can also be derived in length gauge [32]. We also acknowledge Hiruku Watanabe and Yoichi Yanese for informing us of the update of [33]. We finally acknowledge Daniel Kaplan for useful discussion at the finite lifetime limit [31].

NONLINEAR DC PHOTOCURRENT: INFINITE BANDS LIMIT
We start from a quadratic Hamiltonian (S1) At infinite bands limit, the commutation relation holds [1,2]: We define the velocity operator asv For Hamiltonian (S1),v =p me . Taking Bloch state |χ n (k)⟩ as basis, velocity matrix elements v nm (k) are diagonal of k: and can be computed as where u n (k), u m (k) are corresponding periodic functions of χ n (k) and χ m (k), respectively. When n ̸ = m, v nm (k) can be replaced with the position matrix element r nm (k): and r nm (k) is defined through: The first term contains a k derivative to the δ function, which is not well defined mathematically. It vanishes unless n = m. For n = m, we integrate it over k = k ′ : where ϵ − → 0. In the equation above we replace r with R to denote the band-resolved Berry connection. Below we show explicitly how to use commutation relation and position matrix elements to derive nonlinear DC photocurrent at infinite bands limit.
Under velocity gauge, through minimal coupling, the Hamiltonian becomes: From (S3), the velocity operator becomes:v and the corrected velocity matrix where the superscript "A" denotes the minimal coupling. We note in (S11), the correction to velocity matrix element is only at the first order of vector potential A: all higher order corrections vanish due to commutation relation (S2) [1].
To describe the interaction between the external field and electrons of the system, the perturbation is written as is dropped as it is proportional to the identity and will not contribute to the current at quadratic order in the electromagnetic field. For time-dependent perturbation, Von Neumann equation describes the time evolution of the density matrix ρ(t): where subscript "I" denotes the interaction picture. By integrating equation (S13), we can solve ρ I (t) order by order: Here we only present the first order ρ (1) (t) and the second order ρ (2) (t): The macroscopic current is computed as: By separating v(t) and ρ(t) into different orders of A according to (S11) and (S14), the nonlinear photocurrent is: Since v A mnk (t) only has the first order correction of A on diagonal terms, while according to (S15) the correction on ρ (1) nmk (t) is only on off-diagonal terms, the nonlinear photocurrent only comes from v mn (k)ρ (2) nmk (t). In (S16) we require ω 1 + ω 2 = 0 as the condition for DC response, thus ρ 2 nmk (t) is reduced to: As described in the main manuscript, we separate j(t) into diagonal contribution First we deal with j dia . Using identity: j dia can be broken into three pieces. The contribution from the first resonant term is: This is generally referred to as injection current (j inj ). For systems with time-reversal symmetry (TRS), v nn (k) = −v nn (−k), injection current is non-vanishing only with CP light. The second principal term vanishes trivially by switching indices i ↔ j, n ↔ m, ω ′ ↔ −ω ′ . For the third term, using the identity it can be simplified as: Similar to injection current, with time-reversal symmetry (TRS), this term is nonzero only with CP light illumination. Later we show how this term can be combined with j off and constitute the j Fermi . We now switch to j off . By replacing v mn (k) with r mn (k) according to (S6), j off is written as: The first two lines are demonstrated to be zero by rotating l − → n, m − → l and using commutation relation [r i , P j ] = iℏδ ij . Using (S7) and (S8), we get: Similar to diagonal parts, j off can be separated as a resonant contribution related with δ function: and the remaining nonresonant contribution is related with principal part: With TRS, by switching n ↔ m and k ↔ −k, j off2 vanishes with LP light illumination, and resonant j off1 remains. This is generally referred to as shift current (j shift ). With CP light, resonant j off1 vanishes, but nonresonant j off2 is nonvanishing. By switching n ↔ m, i ↔ j, and ω ′ ↔ −ω ′ , the term R mm (k) − R nn (k) v i nm (k)v j mn (k) vanishes trivially. j off2 is reduced as: Above we have inserted a factor of 1 2 due to symmetrization by switching n ↔ m, i ↔ j, and ω ′ ↔ − ω ′ . Inspecting (S23) and (S28), ∇k are applied on different parts of the same equation. Adding them together is equivalent to applying ∇k derivative on (f T n (k, µ) − f T m (k, µ)): In . The current is contributed by excitations from both occupied states to Fermi surface and from Fermi surface to unoccupied states (j Fermi ). By switching n ↔ m, we rewrite (S29) as:

GAUGE INVARIANCE OF NONLINEAR DC PHOTOCURRENT
We add arbitrary phases to the periodic function part of Bloch states u n (k) − → e iϕn(k) u n (k), u m (k) − → e iϕm(k) u m (k) to see whether expressions of photocurrent will change. Since velocity operatorv k does not enclose any phase information, velocity matrix element v nm (k) will transform following: According to Eqs. S21 and S30, j inj and j Fermi are thus gauge invariant. Below we prove the gauge invariance of j shift . By adding arbitrary phases, we have the following transformation: Therefore, j shift is also gauge invariant.

NONLINEAR DC PHOTOCURRENT: TIME-REVERSAL SYMMETRY BROKEN
When TRS is broken in magnetic systems, the relation v nm (k) = −v * nm (−k) will no longer hold. Here we reinspect the expressions for each contribution.
We look at the diagonal contribution j dia = e V n,k v nn (k)ρ nn (k) first. As shown in Eqs. S20-S23, it can be parsed into the resonant injection current j dia1 (j inj ) and nonresonant j dia2 and j dia3 . For the resonant part, without v nn (k) ̸ = −v nn (−k), j inj can also survive with the illumination of LP light. This is referred to as "magnetic injection current" and has been computed recently [3,4]. j dia2 vanishes trivially and the cancellation does not rely on TRS. The remaining j dia3 is associated with the nonresonant part of j off = e V n̸ =m,k v nm (k)ρ mn (k) below and constitutes the j Fermi . For the off-diagonal contribution j off = e V n̸ =m,k v mn (k)ρ nm (k), as shown by Eq. S26, j shift represents the resonant part. When v nm (k) = −v * nm (−k) does not hold, j shift is nonvanishing with both LP and CP light. This is in contrast with nonmagnetic systems where j shift is nonvanishing only with LP light. For the nonresonant part of j off (Eq. S27), the term R mm (k) − R nn (k) v i nm (k)v j mn (k) vanishes trivially independent of TRS. With TRS breaking, the term v i nm (k)i∇ k v j mn (k) is not cancelled between k and −k. Combined with j dia3 , it arrives in the form of Eq. S30 (j Fermi ) and is nonvanishing with both LP and CP light.
As a summary, by breaking TRS, the nonlinear DC photocurrent will have j inj (Eq. S21), j shift (Eq. S26), and j Fermi (Eq. S30) with illumination of both CP and LP light.

NON-LINEAR PHOTOCURRENTS FOR GENERAL BLOCH HAMILTONIANS
Here we give an equivalent derivation of the above contributions to the non-linear DC photocurrent for a general Bloch Hamiltonian. We show that in addition to the shift, injection, and j F ermi currents there is another contribution to the non-linear photocurrent that is only non-vanishing for systems with a Fermi surface. We denote this contribution J Fermi2 .

TIME EVOLUTION OF THE QUANTUM DENSITY MATRIX
Here we solve for the nonlinear response of charge currents to quadratic order in a perturbing electric field by first solving the quantum density matrix,ρ(t), to quadratic order in the external electromagnetic vector potential. To do so we employ the von Neumann equation which describes the time evolution of this quantum operator [5].
The Hamiltonian,Ĥ(r,p, t), can be divided into two parts:Ĥ 0 (r,p) andĤ ′ (r,p, t). The time independent part, H 0 (r,p), describes noninteracting electrons in a crystal potential, and the time dependent part describes the interaction between electrons in the material and the external perturbing electric field. Due to the translation symmetry of the crystal it is advantageous to write the operators in the basis of the Bloch eigenstates of H 0 (r,p): O nm (k, k ′ ) = ⟨u n (k)|Ô |u m (k ′ )⟩, where |u n (k)⟩ is the periodic part of the Bloch wave function in band n, with crystal momentum k, and energy ε n (k). Here we will study spatially homogeneous electric fields, E(r, t) → E(t), that only induce electronic transitions between Bloch states with the same crystal momentum k and thus, define The Bloch HamiltonianĤ 0 (r,p) can be written in the basis of localized atomic orbitals { φ Rj } whose matrix elements we denote as Here R i indexes the position of orbital i in the unit cell at R. Due to translation symmetry H 0 . These orbitals are chosen such that they satisfy where τ i is the intracellular position of orbital i. We can then write the operator aŝ whereĉ † φ (R i ) andĉ φ (R i ) are electronic creation and annihilation operators for the i orbital in the unit cell at R.
To make connection with our Bloch states, we can build Bloch like function out of these localized atomic orbitals by taking the linear combination We can define our Bloch Hamiltonian and its components by taking the inner product of the Hamiltonian with respect to these states HereR = R i − R j and note that H 0 ij (R) is periodic under translation by a lattice vector. We can diagonalize These coefficients u nk (i) make up the components of the periodic part of the Bloch vectors |u n (k)⟩ written in the orbital basis.

Light matter interaction between electrons in the material and an external electromagnetic field can be taken into account via the Peierls substitution [6]
Here we will work in temporal gauge such that the external electromagnetic scalar potential ϕ ext (r, t) can be taken to vanish. For spatially homogeneous electromagnetic vector fields, A(t), we can expand the Hamiltonian in the orbital basis asĤ The external field has modified the matrix elements of the Hamiltonian such that Using equation S41 it can be shown that in the Bloch basis We see that the introduction of spatially homogeneous electromagnetic vector potential has simply boosted the momentum of the matrix elements of the Bloch Hamiltonian.
To extractĤ ′ (r,p, t) one only needs to expandĤ(r,p, t) in a power series of the electromagnetic vector potential A(t). By substitution into the von Neumann equation we can begin to solve for the density matrix in powers of the external electromagnetic field with the identification E(t) = −∂ t A(t) made possible by use of the temporal gauge.
The von Neumann equation is a first order linear differential equation in time. To solve the equation an initial condition on the density matrix must be specified. Before application of the external electromagnetic field the material is in its ground state and can be describe by f T n (k, µ) the Fermi occupation function describing the probability for an electron to be in the Bloch state in band n at momentum k given the temperature T and chemical potential µ of the system.
Here we will mainly be concerned with materials perturbed by light with a few nonzero Fourier components ω. We thus Fourier transform equation S36 and write it in the Bloch basis where ρ nm (k, ω) = dtρ nm (k, t)e −iωt and is the Fourier transform of the perturbed Bloch Hamiltonian written in the basis of eigenstates of H 0 ij (k). In doing so our differential equation has become an algebraic equation one can readily solve in powers of E(t).

QUANTUM CURRENT OPERATOR
Here we wish to determine the induced charge current in a material order by order in the perturbing external electric field. This can be calculated with knowledge ofρ(t) and the quantum current operatorĵ(t) via We consider perturbation from spatially homogeneous electric fields that gives rise to spatially homogeneous currents. We thus can identify the current operator with the product of the electric charge e and the velocity operatorv(t) whose general form isv(t) = i[Ĥ(r,p, t),r]/ℏ. Due to the translation symmetry of the crystal it is convenient to express the velocity operator in the basis of our Bloch like functions χ k i . We take the inner product of the operator with respect to two such states and find It is important to recognize that the representation of the velocity operator on Bloch states only takes this form in the gauge choice on χ k j shown above. With this gauge choice we choose the Bloch state to have the property ψ n(k+G) = |ψ nk ⟩ which implies that u n(k+G) (j) = e −iG·τj u nk (j). This is indeed Here we will be considering spatially homogeneous electric fields with few nonzero frequency components ω, and thus focus on the Fourier transform of the current Again due to the translation symmetry of the crystal, it is convenient to compute the trace of these operators in the Bloch basis. With knowledge of the representation ofv(ω) andρ(ω) in the Bloch basis, equation S52 can be solve order by order in the electric field to compute the induced current.

CHARGE CURRENTS FIRST ORDER IN THE ELECTRIC FIELD
To solve for charge currents linear in the electric field we first expand the quantum density matrix and velocity operator in a power series with respect to A ext (t) where p indexes the order to which ρ There will be two contributions to the induced current at this order in the electric field. We write the total first order current j (1) (ω) as the sum of these two contributions j  B (ω). The first derives from the first order contribution ρ (1) nm (k, ω) to the quantum density matrix traced against the zeroth order contribution to the velocity operator v (0) (ω). The second comes from the zeroth order contribution to the density operator ρ (0) nm (k, ω) traced against the first order contribution to the velocity operatorv (1) (ω).
Using equation S50 and equation S46 the zeroth and first order terms to the velocity operators written in the Bloch basis are v respectively. The zeroth order contribution to the density matrix nm (k, ω) describes the occupation of the electrons in the ground state before the external electric field is applied to the material, which in the Bloch basis is simply the Fermi occupation functions We see that the first contribution to equation S54 can be simply written as To derive the contribution to the current proportional to v (0) nm (k, ω) we must solve the von Neumann equation to obtain ρ (1) nm (k, ω). Equating first order contributions in equation S47 we have Where we have used equation S46 to expand the Bloch Hamiltonian in powers of A(ω). Using equation S55, we can substitute ρ (0) nm (k, ω) into the above and solve for ρ (1) nm (k, ω).
Where we have defined the unperturbed velocity operator v nm (k) = 1/ℏ ⟨u n (k)| ∇ kĤ0 (k) |u m (k)⟩. Tracinĝ ρ (1) (k, ω) againstv (0) (k, ω) gives the second contribution to j (1) (ω). By expressing the electromagnetic vector potential A(ω) in terms of the electric field E(ω) we find the total contribution to j (1) (ω) can then be written as Integration by parts in the first term of equation S59 will lead to a term proportional to the Fermi surface and other terms proportional to matrix elements of the dipole operator R nm (k) = ⟨u n (k)| i∇ |u m (k)⟩. Note that the representation of ther on the period part of the Bloch functions is i∇. These extra terms, when combined with the second term in equation S59 resolve the divergence in the term, as ω → 0, stemming from A(ω) = iE(ω)/ω. This integration by parts and subsequent cancellation leads to an expression matching the canonical equation found when computing the first order current j (1) (ω) in length gauge [2,7].

CHARGE CURRENTS SECOND ORDER IN THE ELECTRIC FIELD
At second order in the electric field, the current j (2) (ω) has three unique contributions deriving from different orders in the power series expansions of the velocity operator and quantum density matrix. We can write them as The second order contribution to the velocity operator written in the Bloch basis is nm (k) and v (2) nm (k, ω), and using equations S58 and S55 the first two contributions j A (ω) and j To derive the last contribution to the second order current j C (ω) we first solve forρ (2) (ω ′ ). We do this in the Bloch basis by collecting terms that are second order in the electric field in equation S47.
Later we will find it useful to breakup the contributions to ρ (2) nm (k, ω) into two pieces. The first piece ρ (2),0 nm (k, ω) is proportional to ρ (0) nm (k, ω) and whose contribution to the current we denote j  nm (k, ω) and whose contribution to the current we denote j (2) C1 (ω). The total contribution to j (2) C (ω) can thus be written as

DC CURRENTS SECOND ORDER IN THE ELECTRIC FIELD
Here we focus on the zero frequency contribution to the second order current j (2) (0). We will show how the contributions to j (2) (0) can be combined in this DC limit such that they reproduce the established shift and injection current contributions to this nonlinear DC current. In doing so we will show that the principle parts of the resonant pieces of the propagators in the expressions for j (2) (0) vanish leaving contributions that are constrained to resonant transitions in the Brillouin zone and transitions on the Fermi surface.
We start by analyzing j C1 (0). We break this contribution into pieces proportional to the diagonal elements of ρ (2),1 nm (k, ω) we denote j (2) C1A (0) and pieces proportional to the off diagonal elements of ρ (2),1 nm (k, ω) we denote j (2) C1B (0). To regulate the fermionic propagators we move the resonances off the real axes by taking ℏω ′ → ℏω ′ +iη for all frequencies. At the end of the calculation we will take the limit as η → 0. Using equation S66 we find j (2) We now investigate the denominator on the right hand side of equation S68. In the limit as η → 0 we have The second contribution to the denominator can be shown to vanish by relabeling the dummy indices in the sum of equation S68 (i ↔ j, ω ′ ↔ −ω ′ , and n ↔ m) leaving two contributions to j C1A (0). Now we investigate j Switching n ↔ l and m ↔ l, and using the relationship between off diagonal elements of the velocity operator and off diagonal elements of the dipole operator: v nm (k)/(ε n (k) − ε m (k) + iη) = iR nm (k)/ℏ for m ̸ = n and η → 0, we get j (2) Using l |u l (k)⟩ ⟨u l (k)| = 1 and ⟨u n (k)| ∂ ki |u m (k)⟩ = −(∂ ki ⟨u n (k)|) |u m (k)⟩ deriving from the orthogonality of the Bloch states, we get We arrive at j (2) We now add j B (0). Using equation S65 and S73 we find We can now break up the denominator into its resonant and principle parts leading to two contributions to j B (0). We now will combine j (2) C1B (0), and j (2) B (0). We use the identity in the principle part of j C1A (0) and then integrate by parts. This leads to a term proportional to k-space derivatives of the Fermi occupation factors and a contribution equal to the negate of the principle part of equation S74. We note that terms in equation S74 proportional to matrix elements of the dipole operator have no principle contribution as can be shown by the relabeling of dummy indices (i ↔ j, ω ↔ −ω, and n ↔ m).
This leads to the following three contributions to the current arising from the j C1B (0), and j B (0) contributions to j (2) (0). We can write them as (S76) The first term, j inj (0), is the injection current. It is divergent for η → 0, but can be regulated with a relaxation time τ (k) which amounts to adding a term in the equation of motion for the density matrix (equation S47) of the form This term encodes the effects of other interactions in the Hamiltonian not considered above that are linear in the deviations of the density matrix from its equilibrium value, ρ nm (k, ω). Its inclusion into the equation of motion of the density matrix ultimately amounts to substitution of η for 1/τ nm (k). The injection current, j inj (0), can be shown to vanish for time reversal symmetric systems perturbed by linear polarized light, by noting the constraint time reversal symmetry puts on the matrix elements of the velocity operator: v nm (k) = −v mn (−k).
The second term in equation S76, j shift (0), is the shift current. It is an intrinsic contribution to j (2) (0) arising from k-space derivatives of matrix elements of the velocity matrix. These k-space derivatives are accompanied with diagonal elements of the dipole operator such that the current remains invariant under a gauge transformation of the Bloch states of the form |u n (k)⟩ → e iϕn(k) |u n (k)⟩ for all n. It can be shown to vanish for time reversal symmetric materials exposed to external electric fields with circular polarization, but is nonzero for electric fields with linear polarization.
Both the injection and shift current derive from electronic transition between Bloch electrons with energy differences ε n (k) − ε m (k) = −ℏω ′ independent of crystal momentum k. Connection with the canonical expressions for these currents in length gauge can be made using this constraint upon replacement of the electromagnetic vector potential with the electric field and the use of the relationship between off diagonal elements of the velocity and dipole operators.
Lastly the third term in equation S76, j Fermi (0), is a contribution to the current proportional to k-space derivatives of the equilibrium Fermi occupation factors f T n (k, µ). At zero temperature ∇ k f T =0 n (k, µ) → ℏv nn (k)δ(ε n (k) − µ) such that the sum over the Brillouin zone is constrained to crystal momentum k along the Fermi surface. For time reversal symmetric systems this contribution to the current can be shown to vanish for external electric fields with linear polarization, but can be nonzero for circular polarized light. The last contributions to j (2) (0) are from j (2) A (0) and j (2) C0 (0). From equation S66 we find j (2) We can perform an integration by parts of j A (0) and add the expression to j C0 (0). Using the relationship between matrix elements of the velocity operator and matrix elements of the dipole operator and the identity ⟨u n (k)| ∂ ki |u m (k)⟩ = −(∂ ki ⟨u n (k)|) |u m (k)⟩ we find the final contribution to the current j Fermi2 (0) = j (2) C0 (0) can be written as At zero temperature, like j Fermi2 (0), this term j Fermi2 (0) has k-space contributions constrained to crystal momentum along the Fermi surface. Due to the second-order derivative of Hamiltonian in k space, j Fermi2 (0) is nonzero for systems whose energy dispersion relation H(k) is higher than second order in k. As shown in Eq. S79, j Fermi2 (0) is independent of TRS but only depends upon the polarization of the light. Since CP light is antisymmetric ( (0) is nonvanishing only with LP light. It is also straightforward to see j Fermi2 (0) is gauge invariant. Our derivation is consistent with formulas derived using Feynman diagrammatic approach [8,9].

EQUIVALENCE BETWEEN VELOCITY GAUGE AND LENGTH GAUGE
The equivalence between velocity gauge and length gauge can be proved if we rewrite our j Fermi (0) (Eq. S76) and j Fermi2 (0) (Eq. S79), which is composed of j (2) A (0) (Eq. S64) and j (2) C0 (0) (Eq. S78). According to the useful relation of a general operatorÔ(k) [1], A (0) can be rewritten as: We denote the first term in Eq. S81 as j Here and below ε nm (k) = ε n (k) − ε m (k) and f nm (k) = f n (k) − f m (k) are defined as the band energy and band occupation difference, respectively. The last step is obtained by switching band indices n ↔ m. In Eq. S82, the first term proportional to ∂ ki ∇ k f n (k) is referred to as the classical Drude contribution derived in length gauge [10], which describes the occupation change perturbed by the electric field at second order. This can be captured by semiclassical transport theory even in the single band case. Using Eq. S6, by converting v mn (k) to R mn (k), the second term in j (2) A1 (0) is written as: In Eq. S83, we break j A1−b (0) into a piece with photon energy ℏω ′ . This is the Berry curvature dipole contribution derived in length gauge [10][11][12]. To prove that, we rewrite this term as: Using Eq. 13 in [7]: and the definition of Berry curvature Ω(k): Eq. S84 can be rewritten as: which is the same as three-index tensor derived in length gauge [12]. The remaining term in j A1−b (0) (we denote it as j (2) A1−b2 (0)) can be added to j Fermi (0). j Fermi (0) is equivalent to: j Fermi (0) = k,n,m,i,j,ω ′ =±ω − e 3 2V ℏ v i nm (k)v j mn (k) ε nm (k) + ℏω ′ ∇ k f nm (k)A i (ω ′ )A j (−ω ′ ) = k,n,m,i,j,ω ′ =±ω e 3 2V ℏ 3 ε nm (k)ε mn (k)R i nm (k)R j mn (k) ε nm (k) + ℏω ′ ∇ k f nm (k)A i (ω ′ )A j (−ω ′ ).
For isotropic Weyl semimetal with two bands, the current is: where ε F is the Fermi energy relative to the Weyl node, k F is the magnitude of the Fermi vector,n(k) is the normal vector to the Fermi surface, and Ω(k) = (Ω x 12 (k), Ω y 12 (k), Ω z 12 (k)), where 1, 2 refer to the low and high energy bands, respectively. For a CP light, the current can also be defined through j i = χ ij (ω)[E(ω) × E * (ω)] j , where χ ij is an imaginary tensor. For CP light adopted above, [E(ω) × E * (ω)] ⃗ e z = 2E x (ω)E y (−ω) = 2ω 2 A x (ω)A y (−ω) ⃗ e z = iE 2 0 /2 ⃗ e z = iω 2 A 2 0 /2 ⃗ e z , and j F ermi,z = χ zz [E(ω) × E * (−ω)] z . j F ermi,x and j F ermi,y can be established in an analogous approach, and the total j F ermi is related with the trace over χ ij : where Q n = 1 2π F S dS · Ω(k) is the charge of the Weyl point n. For time-reversal invariant system, weyl points at K and −K have the same charge so won't annihilate. The overall contribution is: