Difference frequency generation in topological semimetals

F. de Juan, 2 Y. Zhang, T. Morimoto, Y. Sun, Joel E. Moore, and A. G. Grushin Donostia International Physics Center, P. Manuel de Lardizabal 4, 20018 Donostia-San Sebastian, Spain IKERBASQUE, Basque Foundation for Science, Maria Diaz de Haro 3, 48013 Bilbao, Spain Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA Department of Applied Physics, The University of Tokyo, Tokyo, 113-8656, Japan Max Planck Institute for Chemical Physics of Solids, 01187 Dresden, Germany Department of Physics, University of California, Berkeley, CA 94720, USA Institut Neél, CNRS and Université Grenoble Alpes, Grenoble, France (Dated: July 8, 2019)

One of the most striking predictions that follow from the protected point-like band crossings in topological semimetals is the first quantized observable defined for a metal. It is a circular photocurrent that grows in time at a universal rate given by fundamental constants only in chiral topological semimetals [1,2], which lack mirror symmetries [3][4][5]. This steadily growing current, known as injection current, originates from resonant, interband transitions and is proportional to the intensity of light. Its universal rate of growth, given by the trace of the circular photogalvanic tensor β ab , cannot be directly extracted from a steady-state experiment because at times longer than the scattering time t τ , the current saturates to a τ -dependent value, which makes the measurement of the quantized response challenging.
To measure the intrinsic quantized current rate it is rather desirable to use time-dependent electric fields in the form of light pulses of duration shorter than τ , and measure the emitted THz fields. However, typical experiments are performed in the opposite regime [6][7][8], where assumptions on the nature of scattering are required to extract the injection current [9].
In this work, we propose an alternative to access the quantized current rate, which is to measure difference frequency generation (DFG), where two monochromatic light beams of frequencies ω ± ∆ω/2 produce a slowly oscillating current of frequency ∆ω. In the limit ∆ω ω, DFG is formally equivalent to a photogalvanic effect, but if in addition we demand that ∆ω τ −1 , the response is intrinsic and τ independent, exposing the universal quantum.
To show this, we calculate the circular DFG response for all metals in the regime ω ∆ω τ −1 , using the length gauge formalism [10][11][12][13][14]. We find that the interband contribution to circular DFG oscillates exactly out of phase with respect to the incoming light, and is given by the intrinsic injection rate β ab , becoming topologically quantized and independent of material parameters in chiral topological metals. In addition, we find that there is a free-carrier contribution to circular DFG [15] for any metal which oscillates in phase with the incoming light and displays strong ω dependence due to an extra term beyond the semiclassical Berry-dipole [16]. Due to this additional term, which reduces to a universal function for linear band crossings at time-reversal invariant momenta, the total circular DFG in metals can have any phase shift compared to the incoming light. We present predictions for both interband and free-carrier contributions to the circular DFG using several models for topological semimetals as well as ab-initio calculations. Finally, we briefly discuss the subtleties of the free-carrier DC photocurrents in the opposite limit ∆ω τ −1 [17,18]. Our work contributes to elucidate how topological properties, responsible for unusual non-linear effects in topological semimetals [19][20][21][22][23][24][25][26][27][28], determine non-monochromatic responses.
Difference frequency generation in metals -DFG is the current response obtained when shining two monochromatic beams E a i (t) = Re[E a i e −iωit ] with i = 1, 2 with average frequency ω = (ω 1 + ω 2 )/2 and difference ∆ω = ω 1 −ω 2 . For concreteness we consider equal polarizations E a 1 = E a 2 = E a , and assume time-reversal symmetry arXiv:1907.02537v1 [cond-mat.mes-hall] 4 Jul 2019 throughout the manuscript (we consider the general case in the Supplemental Material [29]) . In the mentioned regime ω ∆ω τ −1 , the generated DFG current can be expanded in perturbation theory in ∆ω/ω as Defining γ ab = γ ab 1 +γ ab 2 , the explicit expressions for these tensors are where we assume ω > 0. In these equations C = e 3 / 2 , k = d 3 k/(2π) 3 , f n is the Fermi function which we take at zero temperature, f nm = f n − f m , the Bloch eigenstates are H |n = ω n |n , ω nm = ω n − ω m , ξ a nn = i n|∂ ka n is the diagonal Berry connection, and r a nm = i n|∂ ka m is the interband position matrix element, which is zero unless n = m. We represent derivatives with a comma, as in f n,a = ∂ ka f n , while a semicolon denotes the generalized derivative, as in r b nm;a = ∂ ka r b nm − i(ξ a nn − ξ a mm )r b nm . The Berry curvature is Ω a = abc ξ c nn,b . The tensors β ab and σ abc in Eq. (1) are of interband origin and are the DFG analogs of the monochromatic photogalvanic effects known as injection and shift currents, respectively. The tensor γ ab is a free-carrier, intraband contribution only present in metals. The focus of this work lies on β ab and γ ab which we label as circular DFG contributions since, unlike σ abc , they lead to currents which change sign when the helicity of circularly polarized light is reversed.
The interband contribution to the circular DFG, β ab , has the same functional form as the interband circular photogalvanic effect. As anticipated, this implies that the trace of β ab , and thus the corresponding DFG contribution Eq. (1), will be quantized in terms of fundamental constants in chiral topological semimetals [1,2].
The free-carrier contribution to the circular DFG, γ ab , is closely related to the optical Hall conductivity σ ab Hall (ω) [30]. As we show in [29], γ ad dbc can be obtained from σ bc Hall (ω) by replacing f n → f n,a /ω in its integrand, which explicitly reveals its metallic origin [15].
Within γ ab the free-carrier part γ ab 2 given by Eq. (5), is the well studied Berry-dipole term of semiclassical origin and a Fermi surface property [16,[31][32][33][34][35][36]. The novel contribution we report is γ ab 1 in Eq. (4), which is also written as Fermi surface integral. Note that its integrand depends on all bands so γ ab 1 is not a Fermi surface property. Unlike the 1/ω dependence of the semiclassical contribution, γ ab 1 displays a strong ω dependence away from zero frequency. At every frequency where a new interband transition becomes active, the energy denominator in Eq. (4) vanishes and the response diverges. Interestingly, γ ab 1 can also be written in terms of Kramers-Kronig transforms of resonant contributions, a useful property to calculate γ ab 1 using ab-initio calculations (see Supplemental Material [29]).
DFG of a single Weyl node -As the simplest example, we first present the circular DFG tensors for a general Weyl semimetal node of the form written in terms of the Fermi velocity matrix v ij , and its tilt u i in the crystallographic coordinate system. We present the detailed computation of the tensors β ab , γ ab 1 , γ ab 2 in [29] and expressed their final form in terms of trace and trace-free parts as where , and similar expressions hold for γ ab 1 and γ ab 2 . In terms of the universal photogalvanic constant β 0 ≡ iC/4π = iπe 3 /h 2 the trace parts are and the trace-free parts are where with a = ( 2µ ω − 1)/ũ, g 2 (ω) = arctanh 4ũµω 4µ 2 +(ũ 2 −1)ω 2 and g 3 (ω) = arctanh The trace (β) and traceless (β F ) parts are plotted in Fig. 1 a) and b) for zero and finite tilt. As shown in Ref. [1], β displays a quantized plateau once the resonant manifold of optical transitions becomes closed, because it is then determined by the monopole charge of the node. In the tilted case,ũ = 0 this manifold is open for 2µ/(1 +ũ) < ω < 2µ/(1 −ũ), and β becomes quantized for ω > 2µ/(1 −ũ). β F is finite only for 2µ/(1 +ũ) < ω < 2µ/(1 −ũ), and vanishes in the zero tilt limit. Off-diagonal components of β ab therefore require finite tilt [37][38][39].
The free-carrier trace (γ = γ 1 + γ 2 ) and traceless (γ F = γ 1,F +γ 2,F ) parts are shown in Fig. 1 c) and d). In the zero tilt limit, γ 1 ∝ ω/(4µ 2 − ω 2 ). At finite tilt, the divergence at 2µ splits into two logarithmic divergences at ω ± = 2µ/(1 ±ũ), where γ 1,F also displays strong singularities. We note that γ ab 1 is well behaved when ω → 0. The singularities at ω → 0 originate from γ 2 and γ 2,F . While the trace of the semiclassical free-carrier contribution γ 2 gives a universal 1/ω divergence that is independent of the chemical potential and tilt [35], γ 2,F displays a non-universal 1/ω divergence, which depends on the tilt. Since in any material the number of left and right chiralities must be equal, the sum over Weyl nodes will cancel all the trace contributions γ 2 in pairs. The total free-carrier current, however, may have a non-universal 1/ω pole due to γ 2,F because of its tilt dependence and the tensor structure in Eq. (7).
DFG quantization in chiral topological semimetals -For the trace β to be quantized, a topological semimetal is required to have left and right chirality nodes at different energies [1,2], which is only allowed in a chiral lattice structure, which lacks mirror symmetries.
Consider first two Weyl nodes of opposite chiralities away from time reversal invariant points located at energies µ L and µ R measured from the chemical potential and tiltsũ L andũ R (e.g. SrSi 2 without spin-orbit coupling [4]). In the presence of time-reversal symmetry two more symmetry related nodes exist, which contribute in exactly the same way and simply double the result we present. In Fig. 2 a) and b) we show the circular DFG trace parts β and γ. The quantized plateau seen for β is realized in the range 2µ L /(1 −ũ L ) < ω < 2µ R /(1 +ũ R ), which determines how large the tilts can be before quantization is lost. Note also that, upon summing over the two chiralities, the total trace-part γ presents no 1/ω pole, as discussed above.
Our second example are multifold fermions. They are low energy excitations close to degeneracy nodes where three, four, or six bands meet, and which requires additional crystalline symmetries to remain degenerate [40][41][42][43][44]. They have a definite chirality and larger monopole charge compared to Weyl nodes. In space group 198, to linear order in momentum and neglecting spin-orbit coupling, a threefold fermion exists at the Γ point with a Hamiltonian H = v F k · S, where S is a vector of three spin-1 matrices. At the R point, a fourfold fermion exists composed of two Weyl nodes of equal chiralities, H = −v F k · σ ⊗ 1, separated in energy from the threefold at Γ. This situation is realized in CoSi [45][46][47], RhSi [47], and AlPt [48], which motivate our example.
In Fig. 2 c) and d) we show the circular DFG trace parts β and γ corresponding to a threefold and a fourfold nodes at energies µ 3f and µ 4f (see the Supplemental Material [29] for analytic expressions). Due to the monopole charge carried by the multifold fermions [2], β displays a quantized plateau at β 0 /2, two times that of a Weyl node. The free-carrier part γ displays divergences at every energy where a new band becomes resonant with the Fermi level. As for the chiral Weyl case, upon summing over chiralities the low-energy divergent semiclassical part is absent. The remaining freecarrier contribution is This is a universal function for linear nodal points occurring at time-reversal invariant momenta.
Similarly to multifods, Kramers Weyl [5] (e.g. occurring in elemental Te or TlTe 2 O 6 ) present Weyl nodes at time-reversal invariant momenta separated in energy. The full DFG tensor can be calculated analytically, even with quadratic corrections, and it is detailed in the Supplemental Material [29] since it is conceptually similar to our previous examples.
DFG in TaAs and RhSi from first principles -In the presence of mirror symmetries the trace of γ ab and β ab vanish, with the consequent loss of quantization. In a

FIG. 2. a) and b):
Diagonal parts of the circular DFG, β (blue) and γ (green) for a chiral Weyl semimetal with two Weyl nodes at energies µL and µR, with µL/µR = 0.25, and tiltsũL = 0.2ũR = 0.25. c) and d): Same quantities for the linear approximation of a multifold material in space group 198, with a threefold at µ 3f and a double Weyl fourfold fermion at µ 4f with µ 3f /µ 4f = 0. 25. system with C 4v symmetry like TaAs [49][50][51][52][53], the only allowed component is antisymmetric β xy = −β yx , which can only be due to the traceless part in Eq. (7). The same requirement holds for γ ab . Moreover, from Eq. (7) both tilt and non-trivial v ij are required to produce a finite result in TaAs. If v ij = δ ij this component must also be zero [39,54], regardless of the tilt, since Eq. (7) becomes symmetric. In order to provide a qualitative prediction, we have calculated the different parts of the circular DFG for TaAs using density functional theory (DFT) (see [29] for details). The results are shown in Fig. 3 a)

and b).
We observe that γ yx shows the 1/ω divergence, and several sign changes close to 2µ Wi , where µ Wi is the energy of the two types of Weyl nodes W1 and W2 present in TaAs [49,50,55]. Consistently, β yx exhibits characteristic peaks around 2µ Wi . These features follow qualitatively those expected for the traceless circular DFG components, shown in Fig. 3 b) and d).
Unlike TaAs, materials in the cubic space group 198, such as RhSi, CoSi or AlPt, are chiral and lack mirror symmetry. Taking RhSi as an example, and using the cubic nature of the space group we have calculated the diagonal components of the DFG using DFT (see Fig. 1 c) and d)). Consistent with Ref. [2] we find that β displays a flat region between the activation frequencies of the threefold and fourfold fermions. However the response is not exactly quantized: it is corrected by the presence of First principles prediction for circular DFG. a) and b) show βxy and γxy for TaAs with spin-orbit coupling (SOC). The activation frequencies of the W2 and W1 pairs of nodes at 25 and 45 meV (dashed vertical lines) mark the peaks in a) and b), in qualitative agreement with Fig. 1 b) and d). In c) and d) we show the trace β (with and without SOC) and γxx with SOC for RhSi. The activation frequencies of the threefold at Γ and fourfold fermion at R, 0.1 and 0.625 eV respectively, are marked by dashed vertical lines. The horizontal dashed line in b) marks the effective monopole charge of 4 corresponding to 4β0 = 0.318e 3 / 2 . For b) and d) the total free-carrier contribution γ (green-solid) is composed of the semiclassical (γ2, red-dashed) and the novel free-carrier contribution (γ1, orange-dotted), which qualitatively agree with Fig. 2 d).
additional bands deviating from exact quantization even in the absence of spin-orbit coupling. These corrections are expected to decrease upon decreasing the chemical potential. Lastly, the free-carrier contribution γ xx shows peaks and sign changes close to the activation frequencies of the threefold and fourfold fermions, consistent with Fig. 2 d).
Discussion -In this work, we have described the circular DFG response in metals in the clean limit ω ∆ω τ −1 , where it is an intrinsic response of the band structure. We have found that it is composed of two contributions that oscillate out of phase, β ab and γ ab : the trace of the former is quantized in chiral topological metals, and the latter is a free-carrier contribution that vanishes in the absence of a Fermi surface.
The measurement of the quantized trace of β ab can be achieved if the hierarchy of scales ω ∆ω τ −1 is met. Interestingly evaluating if this condition is met is possible without prior knowledge of τ : when ∆ω τ −1 the interband circular DFG signal oscillates exactly out of phase with the incoming light, while in the dirty limit ∆ω τ −1 it oscillates in phase. Similarly, when ∆ω τ −1 the free-carrier circular DFG γ ab oscillates in-phase and displays a characteristic singular ω dependence. The observation of an out-of-phase smooth component together with an in-phase component with singularities is therefore a strong signal that the clean limit has been achieved.
In the currently available RhSi, the plateau extends up to ω ∼ 0.7 eV [9], so setting ∆ω = 100 meV may allow scattering rates as large as τ −1 ∼ 20 meV or τ ∼ 200 fs. It should be noted that the quantized plateau is corrected in multifold materials by multi-band corrections [2], as confirmed by our DFT calculations, so it is desirable to search for new chiral Weyl semimetals with simpler band structures to reach exact quantization.
The relation of the free-carrier DFG response with the more often measured DC photogalvanic responses (where ∆ω = 0 with finite τ ) is quite subtle. Both general arguments and explicit calculations have been used to argue that the free-carrier part of the DC photocurrent vanishes in the absence of a specific mechanism for dissipation [17,18,56]. In the dirty limit (∆ω τ −1 ) there are τ -independent disorder-induced contributions which can cancel the apparently intrinsic contributions [17], in a situation reminiscent of the anomalous Hall effect. In view of other disorder-induced photocurrent calculations [38,[57][58][59], the extent to which this cancellation happens remains to be understood. Due to our assumption that ∆ω τ −1 the above subtleties do not affect the DFG regime pertinent to this work.
In conclusion, our work shows how difference frequency generation provides a disorder-independent route to separately measure photocurrent quantization, and a novel Fermi surface contribution beyond semiclassics. More broadly, it suggests that non-monochromatic responses encompass a rich set of experimental avenues to expose topology.
Acknowledgements -We are grateful to L. Wu for enlightening discussions, and to L. Golub and M. Teixido for providing us with copies of Ref. [17]. We acknowledge O. Matsyshyn and I. Sodemann for valuable discussions, and for sharing their manuscript on a study employing a similar formalism prior to publication [60]. The length gauge approach to the computation of non-linear response functions has been developed in detail for semiconductors in Refs. [10][11][12][13][14]. The same approach can be used for metals in the non-interacting limit with a small number of modifications. In this formalism, the coupling of the bare Hamiltonian H 0 to the electric field takes the form H = H 0 − e E · r with e < 0, and the current expectation value can be expressed in terms of the matrix elements of the density matrix c nm = a † n (k)a m (k) as the sum of interband and intraband pieces is a generalized velocity which contains interband and intraband anomalous velocities, and Ω a n = abc ξ c nn,b is the Berry curvature. The intraband anomalous velocity, omitted in Ref. [12] but included in Ref. [11], is immaterial for insulators but important for metals, in particular because it gives rise to the semiclassical photogalvanic effects.
The density matrix satisfies the dynamical equation i ∂ t c = [H, c], which for the matrix elements implies To obtain the non-linear response functions, c mn are expanded iteratively to the desired order in the electric field [12], starting with c (0) nm = f n δ nm where f n is the zero temperature Fermi function for band n. For insulators, f n = 1 for occupied bands and f n = 0 otherwise, but in the presence of a Fermi surface the c nm acquire extra pieces because f n depends on k. Expanding in the electric field E b (t) = E b β e −iω β t where β runs over different monochromatic components and is summed over the terms up to second order are with B b mn = fnmr b mn ω β −ωmn and ω Σ = ω β + ω γ . The frequencies ω β should be understood as ω β + i where the limit → 0 is taken after the thermodynamic limit. These coefficients were calculated in Ref. [61], yet we find some sign differences, e.g. the absence of a relative sign between the two terms in Eq. (17).
It is important to note that the diagonal (n = m) parts of these expressions, understood in the limit → 0, are valid only when the frequency denominators are strictly finite. In the absence of scattering mechanisms, taking the limits ω β → 0 or ω Σ → 0 will lead to divergences in c (1) nn and c (2) nn respectively. The dynamical equations for c nn are simply not defined in this limit.
These divergences are physical: they represent the response of an electron system in the idealized case of no relaxation. If we wish to know the response of the system when one of these frequencies is zero, scattering must be included rigorously in the calculation with the quantum kinetic equation [62]. As discussed previously [17,18,56,58], this is of paramount importance in particular for the Fermi surface free-carrier contributions to the photocurrent, because certain disorder-induced mechanisms give rise to photocurrents that turn out to be independent of the disorder details. These contributions are not captured in the phenomenological approach that simply introduces a constant relaxation time in the dynamical equation. The same reasoning applies for the off diagonal c nm terms if the chermical potential is near a band degeneracy where ω nm → 0.
In light of this discussion, the only intrinsic response that can be modeled without a scattering mechanism is the response in the limit ω Σ 1/τ . This situation is familiar from the context of the anomalous Hall effect, which also has intrinsic and disorder induced contributions, and the intrinsic one can be measured as the AC Hall conductivity in the limit ω → 0 with ω 1/τ . In the photocurrent context, the limit ω Σ 1/τ represents the frequency difference response to two slightly detuned monochromatic lasers. The response functions in this limit can also be used to compute the transient photocurrent response to pulsed light when the duration of the pulse is much shorter than τ .
With these caveats in mind, we now proceed to compute the response functions. As a warmup, we first review the linear response, where the presence of a Fermi surface gives to the Drude divergence. Collecting the terms to linear order, the full first order response is given by Where integration by parts and v a nn = ω n,a were used. Note σ ab (ω) = (σ ab (−ω)) * . For further reference we can separate this into the absorptive (or dissipative) part σ ab abs (ω) and reactive (or non-dissipative) part σ ab rea (ω) as Note both parts have real and imaginary parts. In the presence of time reversal symmetry, σ ab abs is real, σ ab rea is imaginary, and both are symmetric under a ↔ b. The equivalent expression in the velocity gauge and how to transform it into this one are discussed in [10], see Eq. 1.9.
Next we consider the second order contributions to the current, which take the general form χ abc inter has the insulator contribution derived in Ref. [12] plus those obtained from the terms containing f n,a in Eq. (18). Due to the extra factor ω Σ , both contributions to ∂ t P inter , which give rise to an optical rectification current, are subdominant when ω Σ is small and will be neglected compared to those in J a intra in this work. The dominant contributions to frequency difference generation come from χ abc intra , which is given by The first two terms in this expression are the standard ones obtained for semiconductors [12]. They can be expressed in terms of single frequency tensors which take the form where C = e 3 / 2 . Using (r c nm ) * = r c mn , (r c nm;a ) * = r c mn;a , we see these single frequency tensors satisfy Γ abc (−ω) = Γ abc (ω) * and Λ abc (−ω) = Λ abc (ω) * so both tensors satisfy the standard Kramers-Kronig relations.
Finally, there is also the new Fermi surface contribution a real current, which are either real and symmetric or imaginary and antisymmetric. Also note in Eq. (36) we used abc Ω c = ξ a,b − ξ b,a , in last term we used v a nn = ω n,a , and we integrated by parts in both terms. Putting all terms together we have While this expression appears to contain reactive terms, we next show that these vanish. For this we follow [14] and and integrate by parts in Eq. (34) to show that where the first term is a new Fermi surface contribution Finally, substituting Eq. (42) into Eq. (41) we get Indeed this response contains no reactive part (even in the absence of time-reversal symmetry). The final expression for the time dependent DFG current (31) takes the form where we defined the total free-carrier contribution ρ abc (ω) = ρ abc 1 (ω) + ρ abc 2 (ω). In the main text we drop the "intra" label since J intra is the leading contribution to the total current in the limit of small ∆ω.
It is worth noting that the total free-carrier contribution takes the form which, remarkably, is obtained from the reactive linear conductivity σ ab rea in Eq. (22) making the replacement f n → f n,a /ω. Therefore, for every contribution to the reactive conductivity there is an intraband photocurrent. And because of the extra k derivative in f n,a , the time-reversal even parts of the photocurrent correspond to the time reversal odd parts of the reactive linear conductivity and viceversa. In the presence of time-reversal invariance, σ ab rea = σ ab Hall is the optical Hall conductivity, and only two terms survive which will lead to the contributions discussed in the main text, γ ab 1 and γ ab 2 , when transformed into two index tensors as discussed below.

Kramers-Kronig relations
Here we show the explicit Kramers-Kronig relations. For any complex function which is analytic in the upper half plane χ(ω) and satisfies χ(−ω) = χ(ω) * we have the usual relations Recalling that σ ab (−ω) = [σ ab (ω)] * , Γ abc (−ω) = [Γ abc (ω)] * and Λ abc (−ω) = [Λ abc (ω)] * , the Kramers-Kronig relations can be used for any of these three tensors directly. However, it should be noted that in the absence of time-reversal symmetry real and imaginary parts do not map to absorptive and reactive parts. Since one is usually interested in computing the absorptive parts numerically because integrating delta functions is stable, and then reproducing the reactive parts with Kramers-Kronig, it is useful to spell out these relations solving for the reactive components. For the non-linear response tensors Γ abc and Λ abc we have One application of these identities is the numerical computation of ρ abc 1 , since its definition in Eq. (43) might be unstable due to the presence of the energy denominator. Solving for ρ abc 1 in Eq. (42) we also find which allows to compute ρ abc 1 via Kramers-Kronig relations. We have found this to be more stable for ab-initio calculations.
Explicit real and imaginary parts Following Ref. [12], all these results can be rewritten in a way that naturally separates real vs imaginary and symmetric vs antisymmetric parts. To do so, we note that the sum over states counts every pair twice, nm = n>m + n<m (the terms with n = m are excluded since f nn = 0). One may rewrite all expressions with sums of the type n>m only by relabeling dummy indices, and then interchanging them, and using that f nm = −f mn , ω nm = −ω mn , r a nm (k) = (r a mn (k)) * . The linear response conductivity reads where we define the sum and difference of energy denominators with the i factors made explicit. These functions should be interpreted as For the non-linear response we find Λ abc abs (ω) = − C 2 k n>m f nm Im(r c nm;a r b mn + r b nm;a r c mn )ImF − + iRe(r c nm;a r b mn − r b nm;a r c mn )ImF + , ρ abc 2 (ω) = − C 2ω k n f n,a iΩ d n dbc + P ω n,bc ω .
Recalling that (r c nm r b mn ) * = r b nm r c mn we see that each of these four terms contains a symmetric, real part which contributes to the linear photogalvanic effect (LGPE) plus an antisymmetric imaginary part which contributes to the circular photogalvanic effect (CPGE). The contributions in Γ abc abs (ω) give rise to an injection current (which grows linearly in time) while those in Λ abc abs (ω) and ρ abc 1,2 (ω) give rise to a current that is constant in time. Finally, we can now check how the different pieces transform under the time reversal operation. This symmetry imposes ω nm (k) = ω nm (−k), r a nm (−k) = r a mn (k), ∆ a nm (−k) = −∆ a nm (k) and r c nm;a (k) = −r c mn;a (−k). The last one follows since the operator D a = ∂ ka + i(ξ a nn − ξ a mm ) gives an extra minus sign. Splitting integrals into k F (k) = 1/2[ k F (k) + −k F (−k)] = 1/2[ k (F (k) + F (−k))]. we realize that ReΓ abc abs (ω) = ImΛ abc abs (ω) = Reρ abc 1 (ω) = Reρ abc 2 (ω) = 0, while the rest of contributions remain the same. With time-reversal symmetry, the injection current and Fermi surface contributions are therefore purely circular, while the current coming from Λ abc abs (ω) is purely linear and is identified as the shift current. Also note that in the limit of ω = 0, ReF − = 0 but ReF + = 1/ω nm is finite, so this term gives an additional semiclassical contribution (beyond the Berry curvature dipole and the Drude peak derivative). This contribution was found in Ref. [63] by deriving the semiclassical equations of motion to second order in the electric field.

Details of the ab-initio calculation
TaAs has space group [67] I4 1 md (#109), with point group 4mm (C 4v ). The two index circular photogalvanic tensors have the same symmetry as the gyrotropic tensor, which in C 4v only has components β xy = −β yx , and similar for γ ab [68]. To calculate the free-carrier current and injection current in TaAs, we obtain the density-functional theory (DFT) Bloch wave functions from the Full-Potential Local-Orbital program (FPLO) [69] within the generalized gradient approximation (GGA) [70]. By projecting the Bloch wave functions onto Wannier functions, we obtain a tight-binding Hamiltonian with 32 bands, which we use for efficient evaluation of the photocurrent. To implement the circular DFG integrals the Brillouin zone was sampled by k-grids from 200 × 200 × 200 to 960 × 960 × 960 [36]. Satisfactory convergence was achieved for a k-grid of size 400 × 400 × 400. To compute γ ab 1 from ab-initio, it is better to go back to the Kramers-Kronig formula because the Fermi surface formula requires the derivative of Fermi-Dirac distribution function on momentum space, which is numerically very unstable and hard to converge. γ ab 1 is obtained via Eq. (68) with the use of Kramers-Kronig transformations.