Difference frequency generation in topological semimetals

TL;DR

This paper derives a universal expression for difference frequency generation (DFG) in topological semimetals, revealing quantized responses under circular polarization and proposing experimental detection methods for topological photocurrents in metals.

Contribution

It provides the first general derivation of DFG in metals, showing quantization in topological semimetals and identifying additional free-carrier contributions beyond semiclassical theory.

Findings

DFG in chiral topological semimetals is quantized in units of e^3/h^2.

Quantized DFG response is independent of material parameters when Δω ≫ τ^{-1}.

Additional free-carrier contributions to DFG oscillate with frequency and phase, beyond semiclassical predictions.

Abstract

When two lasers are applied to a non-centrosymmetric material, light can be generated at the difference of the incoming frequencies $\Delta\omega$, a phenomenon known as difference frequency generation (DFG), well characterized in semiconductors. In this work, we derive a general expression for DFG in metals, which we use to show that the DFG in chiral topological semimetals under circular polarized light is quantized in units of $e^3/h^2$ and independent of material parameters, including the scattering time $\tau$, when $\Delta\omega \gg \tau^{-1}$. In this regime, DFG provides a simpler alternative to measure a quantized response in metals compared to previous proposals based on single frequency experiments. Our general derivation unmasks, in addition, a free-carrier contribution to the circular DFG beyond the semiclassical one. This contribution can be written as a Fermi surface integral, features strong frequency dependence, and oscillates with a $\pi/2$ shift with respect to the quantized contribution. We make predictions for the circular DFG of chiral and non-chiral materials using generic effective models, and ab-initio calculations for TaAs and RhSi. Our work provides a complete picture of the DFG in the length gauge approach, in the clean, non-interacting limit, and highlights a plausible experiment to measure topologically quantizated photocurrents in metals.

Equations

$$J^a(t)=$$

$$+2\cos (\Delta \omega t){\sigma }^{abc}(\omega )(E^b{E}^{c\ast }+E^c{E}^{b\ast }),$$

$${\beta }^{ab}(\omega )$$

$${\sigma }^{abc}(\omega )$$

$${\gamma }_1^{ab}(\omega )$$

$${\gamma }_2^{ab}(\omega )$$

$$H={\sigma }_i{v}_{ij}{k}_j+u_i{k}_i,$$

$${\beta }^{ab}=\chi \left(\frac{{\delta }^{ab}}{3}\beta +\left[\frac{{\delta }^{ab}}{3}-v_{ac}{v}_{bd}^{-1}\frac{{\tilde{u}}_c{\tilde{u}}_d}{{\tilde{u}}^2}\right]{\beta }_F\right),$$

$$\beta$$

$${\gamma }_1$$

$${\beta }_F$$

$${\gamma }_{1,F}$$

$$\displaystyle\Big{.}+3(4\mu^{2}+(1-\tilde{u}^{2})\omega^{2})g_{3}(\omega)\Big{]},$$

$${\gamma }_{2,F}$$

$$J(t)$$

$$P_{\mathrm{inter}}^a(t)$$

$$J_{\mathrm{intra}}^a(t)$$

$${\partial }_t{c}_{mn}+i{\omega }_{mn}{c}_{mn}=\frac{e{E}^b(t)}{\hslash }\left[-c_{mn;b}+i\sum _p{r}_{mp}^b{c}_{pn}-c_{mp}{r}_{pn}^b\right],$$

$$c_{mn}^{(1)}(t)$$

$$c_{mn}^{(2)}(t)$$

$$-\frac{f_{n,bc}{\delta }_{nm}}{{\omega }_{\beta }{\omega }_{\Sigma }}+i\frac{r_{mn}^c{f}_{nm,b}}{{\omega }_{\beta }({\omega }_{\Sigma }-{\omega }_{mn})}+i\frac{r_{mn}^b{f}_{nm,c}}{({\omega }_{\beta }-{\omega }_{mn})({\omega }_{\Sigma }-{\omega }_{mn})}],$$

$$J^a(t)$$

$${\sigma }^{ab}(\omega )$$

$${\sigma }_{\mathrm{abs}}^{ab}=\frac{1}{2}\left[{\sigma }^{ab}+({\sigma }^{ba}{)}^{\ast }\right]$$

$${\sigma }_{\mathrm{rea}}^{ab}=\frac{1}{2}\left[{\sigma }^{ab}-({\sigma }^{ba}{)}^{\ast }\right]$$

$${\partial }_t{P}_{\mathrm{inter}}^a(t)$$

$$J_{\mathrm{intra}}^a(t)$$

$${\chi }_{\mathrm{intra}}^{abc}({\omega }_{\beta },{\omega }_{\gamma })=\frac{1}{-i{\omega }_{\Sigma }}{\eta }^{abc}({\omega }_{\beta },{\omega }_{\gamma })+{\sigma }^{abc}({\omega }_{\beta },{\omega }_{\gamma })+{\rho }_2^{abc}({\omega }_{\beta },{\omega }_{\gamma }).$$

$${\eta }^{abc}({\omega }_{\beta },{\omega }_{\gamma })$$

$${\sigma }^{abc}({\omega }_{\beta },{\omega }_{\gamma })$$

$${\Gamma }^{abc}(\omega )$$

$${\Lambda }^{abc}(\omega )$$

$${\rho }_2^{abc}({\omega }_{\beta },{\omega }_{\gamma })$$

$$J_{\mathrm{intra}}^a(t)$$

$$=2\sum _{{\omega }_{\Sigma }=\pm \Delta \omega }{\chi }_{intra}^{abc,(2)}(\omega -\frac{{\omega }_{\Sigma }}{2},-\omega -\frac{{\omega }_{\Sigma }}{2})e^{-i{\omega }_{\Sigma }t}{E}^b{E}^{c\ast },$$

$${\chi }_{intra}^{abc,(2)}(\omega -\frac{{\omega }_{\Sigma }}{2},-\omega -\frac{{\omega }_{\Sigma }}{2})$$

$$+{\sigma }^{abc}(\omega ,-\omega )+{\rho }_2^{abc}(\omega ,-\omega )+O({\omega }_{\Sigma }).$$

$${\eta }^{abc}(\omega ,-\omega )$$

$${\partial }_{{\omega }_{\Sigma }}{\eta }^{abc}(\omega -\frac{{\omega }_{\Sigma }}{2},-\omega -\frac{{\omega }_{\Sigma }}{2}){|}_{{\omega }_{\Sigma }=0}$$

$${\sigma }^{abc}(\omega ,-\omega )$$

$${\rho }_2^{abc}(\omega ,-\omega )$$

Full text

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 de Juan et al. (2017); Flicker et al. (2018), which lack mirror symmetries Hirayama et al. (2015); Huang et al. (2016); Chang et al. (2018). 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 $\beta^{ab}$, cannot be directly extracted from a steady-state experiment because at times longer than the scattering time $t\gg\tau$, the current saturates to a $\tau$-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 $\tau$, and measure the emitted THz fields. However, typical experiments are performed in the opposite regime Laman et al. (2005); Bieler et al. (2006); Bas et al. (2016), where assumptions on the nature of scattering are required to extract the injection current Rees et al. (2019).

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 $\omega\pm\Delta\omega/2$ produce a slowly oscillating current of frequency $\Delta\omega$. In the limit $\Delta\omega\ll\omega$, DFG is formally equivalent to a photogalvanic effect, but if in addition we demand that $\Delta\omega\gg\tau^{-1}$, the response is intrinsic and $\tau$ independent, exposing the universal quantum.

To show this, we calculate the circular DFG response for all metals in the regime $\omega\gg\Delta\omega\gg\tau^{-1}$, using the length gauge formalism Sipe and Ghahramani (1993); Aversa and Sipe (1995); Sipe and Shkrebtii (2000); Nastos and Sipe (2006, 2010). 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 $\beta^{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 Genkin and Mednis (1968) for any metal which oscillates in phase with the incoming light and displays strong $\omega$ dependence due to an extra term beyond the semiclassical Berry-dipole Sodemann and Fu (2015). 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 $\Delta\omega\ll\tau^{-1}$ Belinicher et al. (1986, 1982). Our work contributes to elucidate how topological properties, responsible for unusual non-linear effects in topological semimetals Yang et al. (2017); Zhang et al. (2018a); Ma et al. (2017); Osterhoudt et al. (2019); Wu et al. (2017); Patankar et al. (2018); Ji et al. (2018); Hosur (2011); Morimoto and Nagaosa (2016); Parker et al. (2018), determine non-monochromatic responses.

Difference frequency generation in metals - DFG is the current response obtained when shining two monochromatic beams $\mathcal{E}_{i}^{a}(t)={\rm Re}[E^{a}_{i}e^{-i\omega_{i}t}]$ with $i=1,2$ with average frequency $\omega=(\omega_{1}+\omega_{2})/2$ and difference $\Delta\omega=\omega_{1}-\omega_{2}$. For concreteness we consider equal polarizations $E^{a}_{1}=E^{a}_{2}=E^{a}$, and assume time-reversal symmetry throughout the manuscript (we consider the general case in the Supplemental Material Sup (2019)) . In the mentioned regime $\omega\gg\Delta\omega\gg\tau^{-1}$, the generated DFG current can be expanded in perturbation theory in $\Delta\omega/\omega$ as

[TABLE]

Defining $\gamma^{ab}=\gamma_{1}^{ab}+\gamma_{2}^{ab}$, the explicit expressions for these tensors are

[TABLE]

where we assume $\omega>0$. In these equations $C=e^{3}/\hbar^{2}$, $\int_{k}=\int d^{3}k/(2\pi)^{3}$, $f_{n}$ is the Fermi function which we take at zero temperature, $f_{nm}=f_{n}-f_{m}$, the Bloch eigenstates are $H\left|n\right>=\omega_{n}\left|n\right>$, $\omega_{nm}=\omega_{n}-\omega_{m}$, $\xi^{a}_{nn}=i\left<n|\partial_{k_{a}}n\right>$ is the diagonal Berry connection, and $r_{nm}^{a}=i\left<n|\partial_{k_{a}}m\right>$ is the interband position matrix element, which is zero unless $n\neq m$. We represent derivatives with a comma, as in $f_{n,a}=\partial_{k_{a}}f_{n}$, while a semicolon denotes the generalized derivative, as in $r^{b}_{nm;a}=\partial_{k_{a}}r_{nm}^{b}-i(\xi^{a}_{nn}-\xi^{a}_{mm})r_{nm}^{b}$. The Berry curvature is $\Omega^{a}=\epsilon^{abc}\xi^{c}_{nn,b}$.

The tensors $\beta^{ab}$ and $\sigma^{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 $\gamma^{ab}$ is a free-carrier, intraband contribution only present in metals. The focus of this work lies on $\beta^{ab}$ and $\gamma^{ab}$ which we label as circular DFG contributions since, unlike $\sigma^{abc}$, they lead to currents which change sign when the helicity of circularly polarized light is reversed.

The interband contribution to the circular DFG, $\beta^{ab}$, has the same functional form as the interband circular photogalvanic effect. As anticipated, this implies that the trace of $\beta^{ab}$, and thus the corresponding DFG contribution Eq. (1), will be quantized in terms of fundamental constants in chiral topological semimetals de Juan et al. (2017); Flicker et al. (2018).

The free-carrier contribution to the circular DFG, $\gamma^{ab}$, is closely related to the optical Hall conductivity $\sigma_{\rm Hall}^{ab}(\omega)$  Steiner et al. (2017). As we show in Sup (2019), $\gamma^{ad}\epsilon^{dbc}$ can be obtained from $\sigma^{bc}_{\rm Hall}(\omega)$ by replacing $f_{n}\to f_{n,a}/\omega$ in its integrand, which explicitly reveals its metallic origin  Genkin and Mednis (1968).

Within $\gamma^{ab}$ the free-carrier part $\gamma_{2}^{ab}$ given by Eq. (5), is the well studied Berry-dipole term of semiclassical origin and a Fermi surface property Deyo et al. (2009); Moore and Orenstein (2010); Sodemann and Fu (2015); Morimoto et al. (2016); Ishizuka et al. (2017); Rostami and Polini (2017); Zhang et al. (2018b). The novel contribution we report is $\gamma_{1}^{ab}$ in Eq. (4), which is also written as Fermi surface integral. Note that its integrand depends on all bands so $\gamma_{1}^{ab}$ is not a Fermi surface property. Unlike the $1/\omega$ dependence of the semiclassical contribution, $\gamma_{1}^{ab}$ displays a strong $\omega$ 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, $\gamma_{1}^{ab}$ can also be written in terms of Kramers-Kronig transforms of resonant contributions, a useful property to calculate $\gamma_{1}^{ab}$ using ab-initio calculations (see Supplemental Material Sup (2019)).

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

[TABLE]

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 $\beta^{ab}$, $\gamma_{1}^{ab}$, $\gamma_{2}^{ab}$ in Sup (2019) and expressed their final form in terms of trace and trace-free parts as

[TABLE]

where $\chi=\tfrac{{\rm det}v}{|{\rm det}v|}=\pm 1$ is the chirality of the node, $v^{-1}_{ij}\equiv(v^{-1})_{ij}$ and $\tilde{u}_{i}=v^{-1}_{ij}u_{j}$, $\tilde{u}=\sqrt{\tilde{u}_{i}^{2}}$, and similar expressions hold for $\gamma_{1}^{ab}$ and $\gamma_{2}^{ab}$. In terms of the universal photogalvanic constant $\beta_{0}\equiv iC/4\pi=i\pi e^{3}/h^{2}$ the trace parts are

[TABLE]

and the trace-free parts are

[TABLE]

where $g_{1}(\omega)=a\Theta(1+a)\Theta(1-a)+\Theta(a-1)-\Theta(-a-1)$ with $a=(\tfrac{2\mu}{\omega}-1)/\tilde{u}$, $g_{2}(\omega)={\rm arctanh}\left(\frac{4\tilde{u}\mu\omega}{4\mu^{2}+(\tilde{u}^{2}-1)\omega^{2}}\right)$ and $g_{3}(\omega)={\rm arctanh}\left(\frac{2\tilde{u}\omega^{2}}{-4\mu^{2}+(\tilde{u}^{2}+1)\omega^{2}}\right)$.

The trace ($\beta$) and traceless ($\beta_{F}$) parts are plotted in Fig. 1 a) and b) for zero and finite tilt. As shown in Ref. de Juan et al. (2017), $\beta$ 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, $\tilde{u}\neq 0$ this manifold is open for $2\mu/(1+\tilde{u})<\omega<2\mu/(1-\tilde{u})$, and $\beta$ becomes quantized for $\omega>2\mu/(1-\tilde{u})$. $\beta_{F}$ is finite only for $2\mu/(1+\tilde{u})<\omega<2\mu/(1-\tilde{u})$, and vanishes in the zero tilt limit. Off-diagonal components of $\beta^{ab}$ therefore require finite tilt Chan et al. (2017); König et al. (2017); Golub et al. (2017).

The free-carrier trace ($\gamma=\gamma_{1}+\gamma_{2}$) and traceless ($\gamma_{F}=\gamma_{1,F}+\gamma_{2,F}$) parts are shown in Fig. 1 c) and d). In the zero tilt limit, $\gamma_{1}\propto\omega/(4\mu^{2}-\omega^{2})$. At finite tilt, the divergence at $2\mu$ splits into two logarithmic divergences at $\omega_{\pm}=2\mu/(1\pm\tilde{u})$, where $\gamma_{1,F}$ also displays strong singularities. We note that $\gamma_{1}^{ab}$ is well behaved when $\omega\to 0$. The singularities at $\omega\to 0$ originate from $\gamma_{2}$ and $\gamma_{2,F}$. While the trace of the semiclassical free-carrier contribution $\gamma_{2}$ gives a universal $1/\omega$ divergence that is independent of the chemical potential and tilt Rostami and Polini (2017), $\gamma_{2,F}$ displays a non-universal $1/\omega$ 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 $\gamma_{2}$ in pairs. The total free-carrier current, however, may have a non-universal $1/\omega$ pole due to $\gamma_{2,F}$ because of its tilt dependence and the tensor structure in Eq. (7).

DFG quantization in chiral topological semimetals - For the trace $\beta$ to be quantized, a topological semimetal is required to have left and right chirality nodes at different energies de Juan et al. (2017); Flicker et al. (2018), 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 $\mu_{L}$ and $\mu_{R}$ measured from the chemical potential and tilts $\tilde{u}_{L}$ and $\tilde{u}_{R}$ (e.g. SrSi2 without spin-orbit coupling Huang et al. (2016)). 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 $\beta$ and $\gamma$. The quantized plateau seen for $\beta$ is realized in the range $2\mu_{L}/(1-\tilde{u}_{L})<\omega<2\mu_{R}/(1+\tilde{u}_{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 $\gamma$ presents no $1/\omega$ 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 Mañes (2012); Bradlyn et al. (2016); Chang et al. (2017); Tang et al. (2017); Bouhon and Black-Schaffer (2017). 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 $\Gamma$ point with a Hamiltonian $H=v_{F}\mathbf{k}\cdot\mathbf{S}$, where $\mathbf{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}\mathbf{k}\cdot\mathbf{\sigma}\otimes 1$, separated in energy from the threefold at $\Gamma$. This situation is realized in CoSi Takane et al. (2019); Rao et al. (2019); Sanchez et al. (2019), RhSi Sanchez et al. (2019), and AlPt Schröter et al. (2019), which motivate our example.

In Fig. 2 c) and d) we show the circular DFG trace parts $\beta$ and $\gamma$ corresponding to a threefold and a fourfold nodes at energies $\mu_{3f}$ and $\mu_{4f}$ (see the Supplemental Material Sup (2019) for analytic expressions). Due to the monopole charge carried by the multifold fermions Flicker et al. (2018), $\beta$ displays a quantized plateau at $\beta_{0}/2$, two times that of a Weyl node. The free-carrier part $\gamma$ 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 free-carrier contribution is $\gamma_{1}=(\beta_{0}/\pi)\sum_{i}\chi_{i}\omega/(\mu_{i}^{2}-\omega^{2})$ where $\mu_{i}=\mu_{3f},2\mu_{4f}$ and $\chi_{i}=\pm$. This is a universal function for linear nodal points occurring at time-reversal invariant momenta.

Similarly to multifods, Kramers Weyl Chang et al. (2018) (e.g. occurring in elemental Te or TlTe2O6) 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 Sup (2019) 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 $\gamma^{ab}$ and $\beta^{ab}$ vanish, with the consequent loss of quantization. In a system with $C_{4v}$ symmetry like TaAs Weng et al. (2015); Huang et al. (2015); Lv et al. (2015a); Xu et al. (2015); Yang et al. (2015), the only allowed component is antisymmetric $\beta^{xy}=-\beta^{yx}$, which can only be due to the traceless part in Eq. (7). The same requirement holds for $\gamma^{ab}$. Moreover, from Eq. (7) both tilt and non-trivial $v_{ij}$ are required to produce a finite result in TaAs. If $v_{ij}=\delta_{ij}$ this component must also be zero Golub et al. (2017); Leppenen et al. (2019), 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 Sup (2019) for details). The results are shown in Fig. 3 a) and b). We observe that $\gamma^{yx}$ shows the $1/\omega$ divergence, and several sign changes close to $2\mu_{W_{i}}$, where $\mu_{W_{i}}$ is the energy of the two types of Weyl nodes W1 and W2 present in TaAs Weng et al. (2015); Huang et al. (2015); Buckeridge et al. (2016). Consistently, $\beta^{yx}$ exhibits characteristic peaks around $2\mu_{W_{i}}$. 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. Flicker et al. (2018) we find that $\beta$ 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 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 $\gamma_{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 $\omega\gg\Delta\omega\gg\tau^{-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, $\beta^{ab}$ and $\gamma^{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 $\beta^{ab}$ can be achieved if the hierarchy of scales $\omega\gg\Delta\omega\gg\tau^{-1}$ is met. Interestingly evaluating if this condition is met is possible without prior knowledge of $\tau$: when $\Delta\omega\gg\tau^{-1}$ the interband circular DFG signal oscillates exactly out of phase with the incoming light, while in the dirty limit $\Delta\omega\ll\tau^{-1}$ it oscillates in phase. Similarly, when $\Delta\omega\gg\tau^{-1}$ the free-carrier circular DFG $\gamma^{ab}$ oscillates in-phase and displays a characteristic singular $\omega$ 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 $\hbar\omega\sim 0.7$ eV Rees et al. (2019), so setting $\Delta\omega=100$ meV may allow scattering rates as large as $\tau^{-1}\sim 20$ meV or $\tau\sim 200$ fs. It should be noted that the quantized plateau is corrected in multifold materials by multi-band corrections Flicker et al. (2018), 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 $\Delta\omega=0$ with finite $\tau$) 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 Belinicher et al. (1986); Belinicher and Sturman (1980); Belinicher et al. (1982). In the dirty limit ($\Delta\omega\ll\tau^{-1}$) there are $\tau$-independent disorder-induced contributions which can cancel the apparently intrinsic contributions Belinicher et al. (1986), in a situation reminiscent of the anomalous Hall effect. In view of other disorder-induced photocurrent calculations König et al. (2017); Isobe et al. (2018); Golub and Ivchenko (2018); Nandy and Sodemann (2019), the extent to which this cancellation happens remains to be understood. Due to our assumption that $\Delta\omega\gg\tau^{-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. Belinicher et al. (1986). 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 Matsyshyn and Sodemann (2019). A. G. G. acknowledges financial support from the Marie Curie programme under EC Grant agreement No. 653846. J. E. M. was supported by the Quantum Materials program at LBNL, funded by the US Department of Energy under Contract No. DE-AC02-05CH11231.

Appendix A Appendix A: Difference frequency generation in the length gauge

The length gauge approach to the computation of non-linear response functions has been developed in detail for semiconductors in Refs. Sipe and Ghahramani (1993); Aversa and Sipe (1995); Sipe and Shkrebtii (2000); Nastos and Sipe (2006, 2010). 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\vec{E}\cdot\vec{r}$ with $e<0$, and the current expectation value can be expressed in terms of the matrix elements of the density matrix $c_{nm}=\left<a^{\dagger}_{n}(k)a_{m}(k)\right>$ as the sum of interband and intraband pieces

[TABLE]

where $\int_{k}=\int d^{3}k/(2\pi)^{3}$ and $D^{a}_{nm}=v_{nn}^{a}\delta_{nm}-(e/\hbar)E^{b}(t)[r^{b}_{nm;a}+\delta_{nm}\epsilon^{dab}\Omega^{d}_{n}]$ is a generalized velocity which contains interband and intraband anomalous velocities, and $\Omega^{a}_{n}=\epsilon^{abc}\xi^{c}_{nn,b}$ is the Berry curvature. The intraband anomalous velocity, omitted in Ref. Sipe and Shkrebtii (2000) but included in Ref. Aversa and Sipe (1995), 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\hbar\partial_{t}c=[H,c]$, which for the matrix elements implies

[TABLE]

To obtain the non-linear response functions, $c_{mn}$ are expanded iteratively to the desired order in the electric field Sipe and Shkrebtii (2000), starting with $c_{nm}^{(0)}=f_{n}\delta_{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}_{\beta}e^{-i\omega_{\beta}t}$ where $\beta$ runs over different monochromatic components and is summed over the terms up to second order are

[TABLE]

with $\mathcal{B}^{b}_{mn}=\frac{f_{nm}r_{mn}^{b}}{\omega_{\beta}-\omega_{mn}}$ and $\omega_{\Sigma}=\omega_{\beta}+\omega_{\gamma}$. The frequencies $\omega_{\beta}$ should be understood as $\omega_{\beta}+i\epsilon$ where the limit $\epsilon\rightarrow 0$ is taken after the thermodynamic limit. These coefficients were calculated in Ref. Hipolito et al. (2016), 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 $\epsilon\rightarrow 0$, are valid only when the frequency denominators are strictly finite. In the absence of scattering mechanisms, taking the limits $\omega_{\beta}\rightarrow 0$ or $\omega_{\Sigma}\rightarrow 0$ will lead to divergences in $c_{nn}^{(1)}$ and $c_{nn}^{(2)}$ 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 Culcer et al. (2017). As discussed previously Belinicher et al. (1986); Belinicher and Sturman (1980); Belinicher et al. (1982); Golub and Ivchenko (2018), 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 $\omega_{nm}\rightarrow 0$.

In light of this discussion, the only intrinsic response that can be modeled without a scattering mechanism is the response in the limit $\omega_{\Sigma}\gg 1/\tau$. 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 $\omega\rightarrow 0$ with $\omega\gg 1/\tau$. In the photocurrent context, the limit $\omega_{\Sigma}\gg 1/\tau$ 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 $\tau$.

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

[TABLE]

Where integration by parts and $v_{nn}^{a}=\omega_{n,a}$ were used. Note $\sigma^{ab}(\omega)=(\sigma^{ab}(-\omega))^{*}$. For further reference we can separate this into the absorptive (or dissipative) part $\sigma^{ab}_{\rm abs}(\omega)$ and reactive (or non-dissipative) part $\sigma^{ab}_{\rm rea}(\omega)$ as

[TABLE]

Note both parts have real and imaginary parts. In the presence of time reversal symmetry, $\sigma^{ab}_{\rm abs}$ is real, $\sigma^{ab}_{\rm rea}$ is imaginary, and both are symmetric under $a\leftrightarrow b$. The equivalent expression in the velocity gauge and how to transform it into this one are discussed in Sipe and Ghahramani (1993), see Eq. 1.9.

Next we consider the second order contributions to the current, which take the general form

[TABLE]

$\chi^{abc}_{\rm inter}$ has the insulator contribution derived in Ref. Sipe and Shkrebtii (2000) plus those obtained from the terms containing $f_{n,a}$ in Eq. (18). Due to the extra factor $\omega_{\Sigma}$, both contributions to $\partial_{t}P_{\rm inter}$, which give rise to an optical rectification current, are subdominant when $\omega_{\Sigma}$ is small and will be neglected compared to those in $J_{\rm intra}^{a}$ in this work.

The dominant contributions to frequency difference generation come from $\chi^{abc}_{\rm intra}$, which is given by

[TABLE]

The first two terms in this expression are the standard ones obtained for semiconductors Sipe and Shkrebtii (2000). They can be expressed in terms of single frequency tensors

[TABLE]

which take the form

[TABLE]

where $C=e^{3}/\hbar^{2}$. Using $(r^{c}_{nm})^{*}=r^{c}_{mn}$, $(r^{c}_{nm;a})^{*}=r^{c}_{mn;a}$, we see these single frequency tensors satisfy $\Gamma^{abc}(-\omega)=\left[\Gamma^{abc}(\omega)\right]^{*}$ and $\Lambda^{abc}(-\omega)=\left[\Lambda^{abc}(\omega)\right]^{*}$ so both tensors satisfy the standard Kramers-Kronig relations.

Finally, there is also the new Fermi surface contribution

[TABLE]

In the last term we have used $\frac{1}{\omega_{\beta}}+\frac{1}{\omega_{\gamma}}=\frac{\omega_{\Sigma}}{\omega_{\beta}\omega_{\gamma}}$. $\rho^{abc}_{2}$ is new compared to the insulating case and contains two terms: the first reproduces the the semiclassical contributions Deyo et al. (2009); Moore and Orenstein (2010); Sodemann and Fu (2015); Morimoto et al. (2016); Ishizuka et al. (2017); Rostami and Polini (2017); Zhang et al. (2018b) originating in the diagonal anomalous velocity. The second term is purely classical Sodemann and Fu (2015).

A.1 Frequency difference generation at small $\Delta\omega$

We now consider the second order response to two monochromatic beams of frequencies $\omega_{1}$ and $\omega_{2}$, $E^{b}(t)=(E^{b}_{1}e^{-i\omega_{1}t}+E^{b*}_{1}e^{i\omega_{1}t}+E^{b}_{2}e^{-i\omega_{2}t}+E^{b*}_{2}e^{i\omega_{2}t})$. Our aim is to take slightly detuned frequencies $\omega_{1}=\omega-\tfrac{\Delta\omega}{2}$ and $\omega_{2}=\omega+\tfrac{\Delta\omega}{2}$ and look at the frequency difference generation terms which oscillate with a low frequency $\Delta\omega=\omega_{2}-\omega_{1}$ (i.e. when $\omega_{\Sigma}=\pm\Delta\omega$). Keeping only these terms in the sum in Eq. (24) we have

[TABLE]

where we have used the symmetry of $\chi^{abc,(2)}(\omega_{\beta},\omega_{\gamma})=\chi^{acb,(2)}(\omega_{\gamma},\omega_{\beta})$ which originates from the exchange of dummy indices, and in the last step we have considered the case where the two beams have the same polarization $E_{1}^{b}=E_{2}^{b}$ for simplicity.

We now consider the limit where $|\omega_{\Sigma}|\ll\omega$. The dominant term in this limit comes from the $1/\omega_{\Sigma}$ pole accompanying $\eta^{abc}$ in Eq. (25). Since we are also interested in the next leading order terms which are independent of $\omega_{\Sigma}$, we need to expand $\eta^{abc}$ to first order in $\omega_{\Sigma}$. We therefore have

[TABLE]

When evaluating these terms explicitly, they will naturally separate into resonant (or absorptive) terms, proportional to delta functions, and non-resonant (or reactive) terms proportional to principal value integrals. In the absence of time reversal symmetry, both resonant and non-resonant parts are in general complex. Using that $\Gamma^{acb}(-\omega)=\left[\Gamma^{acb}(\omega)\right]^{*}$ and $\Lambda^{acb}(-\omega)=\left[\Lambda^{acb}(\omega)\right]^{*}$, the four terms in Eq. (32) can be computed as

[TABLE]

where the different absorptive and reactive parts take the form

[TABLE]

Note the covariant derivatives in Eq. (40) turned into regular ones because the connections drop out. Also note that by construction $\Lambda_{\rm abs}^{abc}$ and $\Lambda^{abc}_{\rm rea}$ denote only the resonant and non-resonant parts of Eq. (29) which give rise to a real current, which are either real and symmetric or imaginary and antisymmetric. Also note in Eq. (36) we used $\epsilon^{abc}\Omega^{c}=\xi_{a,b}-\xi_{b,a}$, in last term we used $v^{a}_{nn}=\omega_{n,a}$, and we integrated by parts in both terms.

Putting all terms together we have

[TABLE]

While this expression appears to contain reactive terms, we next show that these vanish. For this we follow Nastos and Sipe (2010) and use $\partial_{\omega}\left(\tfrac{1}{\omega_{nm}-\omega}\right)=-\partial_{\omega_{mn}}\left(\tfrac{1}{\omega_{nm}-\omega}\right)=-\frac{1}{\Delta^{a}_{mn}}\partial_{k_{a}}\left(\tfrac{1}{\omega_{nm}-\omega}\right)$ and integrate by parts in Eq. (34) to show that

[TABLE]

where the first term is a new Fermi surface contribution

[TABLE]

Finally, substituting Eq. (42) into Eq. (41) we get

[TABLE]

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

[TABLE]

where we defined the total free-carrier contribution $\rho^{abc}(\omega)=\rho_{1}^{abc}(\omega)+\rho_{2}^{abc}(\omega)$. In the main text we drop the “intra” label since $J_{\rm intra}$ is the leading contribution to the total current in the limit of small $\Delta\omega$.

It is worth noting that the total free-carrier contribution takes the form

[TABLE]

which, remarkably, is obtained from the reactive linear conductivity $\sigma^{ab}_{\rm rea}$ in Eq. (22) making the replacement $f_{n}\rightarrow f_{n,a}/\omega$. 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, $\sigma^{ab}_{\rm rea}=\sigma^{ab}_{\rm Hall}$ is the optical Hall conductivity, and only two terms survive

[TABLE]

which will lead to the contributions discussed in the main text, $\gamma^{ab}_{1}$ and $\gamma^{ab}_{2}$, when transformed into two index tensors as discussed below.

Appendix B Kramers-Kronig relations

Here we show the explicit Kramers-Kronig relations. For any complex function which is analytic in the upper half plane $\chi(\omega)$ and satisfies $\chi(-\omega)=\chi(\omega)^{*}$ we have the usual relations

[TABLE]

Recalling that $\sigma^{ab}(-\omega)=[\sigma^{ab}(\omega)]^{*}$, $\Gamma^{abc}(-\omega)=[\Gamma^{abc}(\omega)]^{*}$ and $\Lambda^{abc}(-\omega)=[\Lambda^{abc}(\omega)]^{*}$, 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 $\Gamma^{abc}$ and $\Lambda^{abc}$ we have

[TABLE]

One application of these identities is the numerical computation of $\rho_{1}^{abc}$, since its definition in Eq. (43) might be unstable due to the presence of the energy denominator. Solving for $\rho_{1}^{abc}$ in Eq. (42) we also find

[TABLE]

which allows to compute $\rho_{1}^{abc}$ via Kramers-Kronig relations. We have found this to be more stable for ab-initio calculations.

B.1 Explicit real and imaginary parts

Following Ref. Sipe and Shkrebtii (2000), 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, $\sum_{nm}=\sum_{n>m}+\sum_{n<m}$ (the terms with $n=m$ are excluded since $f_{nn}=0$). One may rewrite all expressions with sums of the type $\sum_{n>m}$ only by relabeling dummy indices, and then interchanging them, and using that $f_{nm}=-f_{mn}$, $\omega_{nm}=-\omega_{mn}$, $r^{a}_{nm}(k)=(r_{mn}^{a}(k))^{*}$. The linear response conductivity reads

[TABLE]

where we define the sum and difference of energy denominators

[TABLE]

with the $i\epsilon$ factors made explicit. These functions should be interpreted as $F_{\pm}=\frac{1}{\omega_{mn}-\omega}\pm\frac{1}{\omega_{mn}+\omega}+i\pi[-\delta(\omega_{mn}-\omega)\pm\delta(\omega_{mn}+\omega)]$ and satisfy

[TABLE]

For the non-linear response we find

[TABLE]

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 $\Gamma_{\rm abs}^{abc}(\omega)$ give rise to an injection current (which grows linearly in time) while those in $\Lambda_{\rm abs}^{abc}(\omega)$ and $\rho_{1,2}^{abc}(\omega)$ 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 $\omega_{nm}(k)=\omega_{nm}(-k)$, $r^{a}_{nm}(-k)=r_{mn}^{a}(k)$, $\Delta^{a}_{nm}(-k)=-\Delta^{a}_{nm}(k)$ and $r_{nm;a}^{c}(k)=-r_{mn;a}^{c}(-k)$. The last one follows since the operator $D_{a}=\partial_{k_{a}}+i(\xi^{a}_{nn}-\xi^{a}_{mm})$ gives an extra minus sign. Splitting integrals into $\int_{k}F(k)=1/2[\int_{k}F(k)+\int_{-k}F(-k)]=1/2[\int_{k}(F(k)+F(-k))]$. we realize that

[TABLE]

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 $\Lambda_{\rm abs}^{abc}(\omega)$ is purely linear and is identified as the shift current.

Also note that in the limit of $\omega=0$, ${\rm Re}F_{-}=0$ but ${\rm Re}F_{+}=1/\omega_{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. Gao et al. (2014) by deriving the semiclassical equations of motion to second order in the electric field.

Appendix C Mapping to two index tensors

To connect to Eq. (1) in the main text we rewrite Eq. (45) in terms of two rank tensors:

[TABLE]

where we used that in the presence of time reversal symmetry $\Gamma^{abc}_{\rm abs},\rho_{1}^{abc}$ and $\rho_{2}^{abc}$ are antisymmetric and can be written in terms of a two index tensor, such that $\beta^{ad}=\tfrac{1}{2}\epsilon^{dbc}\Gamma_{\rm abs}^{abc}$, $\gamma^{ad}_{1}=\tfrac{1}{2}\epsilon^{dbc}\rho_{1}^{abc}$ and $\gamma^{ad}_{2}=\tfrac{1}{2}\epsilon^{dbc}\rho_{2}^{abc}$. Using the definitions Eqs. (58) to (61), assuming $\omega>0$, and using that ${\rm Re}[F_{-}(\omega_{mn},\omega)]=2\omega/(\omega_{nm}^{2}-\omega^{2})$ and that ${\rm Im}[F_{+}(\omega_{mn},\omega)]=\pi[-\delta(\omega_{mn}-\omega)+\delta(\omega_{mn}+\omega)]$ we arrive at the equations in the main text which we repeat here for completeness

[TABLE]

The analog of Eq. (54) for the two-index tensor $\gamma_{1}^{ab}$ assuming time-reversal symmetry is

[TABLE]

where

[TABLE]

Appendix D Appendix B: Calculation for tilted Weyl node

Here we present the computation of the circular photogalvanic tensors $\beta^{ab}$, $\gamma_{1}^{ab}$, $\gamma_{2}^{ab}$ for a general two band model of the form

[TABLE]

where $v_{ij}$ is a general Fermi velocity matrix that can always be made symmetric by an appropriate choice of the $\sigma_{i}$, $u_{i}$ is the tilt vector and $\mu>0$ is the chemical potential. Some of our results generalize those presented in Refs.Chan et al. (2017); Yang et al. (2017); Carbotte (2016); Steiner et al. (2017).

We work in the coordinate system given by crystallographic axes, where both $v_{ij}$ and $u_{i}$ will in general be arbitrary, as for a Weyl node at a generic point in the Brillouin zone there are no symmetry constraints on them. Symmetries will constrain the total tensors after adding the contributions of all nodes related by symmetry. This model has energies $\omega_{1}=\vec{u}\vec{k}-\sqrt{v_{ij}k_{j}v_{il}k_{l}}$ and $\omega_{2}=\vec{u}\vec{k}+\sqrt{v_{ij}k_{j}v_{il}k_{l}}$ and the Fermi functions are $f_{1}=1$ and $f_{2}=\Theta(\mu-\omega_{2})$.

The expressions for the CPGE tensors in the two-band, clean limit read

[TABLE]

The calculation can be simplified by defining a new coordinate $k^{\prime}_{i}=v_{ij}k_{j}$, in terms of which we have

[TABLE]

Using these relations, we can perform the change of variables in the integral with $\int_{k}=\tfrac{1}{|{\rm det}v|}\int_{k^{{}^{\prime}}}$, obtaining

[TABLE]

where $\chi=\tfrac{{\rm det}v}{|{\rm det}v|}=\pm 1$ is the chirality of the node and

[TABLE]

where all quantities are now computed from the Hamiltonian $H=\sigma_{i}k_{i}^{\prime}+\tilde{u}_{i}k^{\prime}_{i}$, with $\tilde{u}_{i}=(v^{-1})_{ij}u_{j}$.

The tensor structure of $\beta^{\prime ij}$, $\gamma_{1}^{\prime ij}$ and $\gamma_{2}^{\prime ij}$ can only contain $\tilde{u}_{i}\tilde{u}_{j}/\tilde{u}^{2}$ and $\delta^{ij}$, and by separating into trace and traceless prefactors we have

[TABLE]

Knowing the tensor structure, we can compute the coefficients by choosing spherical coordinates where the $z$ axis is aligned with the tilt, $\tilde{u}_{i}=(0,0,\tilde{u})$, so that $\beta=\beta^{\prime ab}\delta^{ab}$ and $\beta_{F}$ can be computed as $\beta_{F}=(\beta-3\beta^{\prime zz})/2$. Note that

[TABLE]

where $v^{a}=\hat{k}^{a}+\tilde{u}^{a}$ and $\epsilon_{k}=k+\tilde{u}k\cos\theta$.

The computation of $\beta^{\prime ab}$ can be done as follows

[TABLE]

To evaluate this integral define $a=(\tfrac{2\mu}{\omega}-1)/\tilde{u}$ and $g_{1}(\omega)=a\Theta(1+a)\Theta(1-a)+\Theta(a-1)-\Theta(-a-1)$. Changing variables to $x=\cos\theta$ we have

[TABLE]

so $\beta_{F}=-\tfrac{i\pi C}{8}\tfrac{(-g_{1}(\omega)+g^{3}_{1}(\omega))}{(2\pi)^{2}}$. Note that when $\tilde{u}\to 0$ then $g_{1}(\omega)\to-1$ for $\omega>2\mu$ and we get $\beta=-\beta_{0}/4$ with $\beta_{0}=\tfrac{i\pi e^{3}}{h^{2}}$. Together with the factor 4 in Eq. (63) the trace of the interband DFG is quantized to $\beta_{0}$.

Along similar lines the computation of $\gamma_{1}^{ab}$ reduces to

[TABLE]

The integrals can be evaluated in Mathematica resulting in the expressions provided in the main text:

[TABLE]

where $g_{2}(\omega)={\rm arctanh}\left(\frac{4\tilde{u}\mu\omega}{4\mu^{2}+(\tilde{u}^{2}-1)\omega^{2}}\right)$ and $g_{3}(\omega)={\rm arctanh}\left(\frac{2\tilde{u}\omega^{2}}{-4\mu^{2}+(\tilde{u}^{2}+1)\omega^{2}}\right)$. Using $\gamma_{1,F}=(\gamma_{1}-3\gamma^{\prime zz}_{1})/2$ we obtain

[TABLE]

which is given in the main text. It is also possible to check that these expressions indeed coincide with those obtained with the Kramers-Kronig transformations, given by Eq. (68).

Finally, the computation of the semiclassical contribution $\gamma_{2}^{ab}$ is as follows

[TABLE]

where $k(\mu,\theta)=\mu/(1+\tilde{u}\cos\theta)$. Using $\Omega^{a}=\tfrac{k^{a}}{2k^{3}}$ we get

[TABLE]

where we have assumed $\tilde{u}<v_{F}$. Decomposing intro trace and traceless parts with Eq. (7), we only need to compute

[TABLE]

Using $\gamma_{2,F}=(\gamma_{2}-3\gamma^{\prime zz}_{2})/2$ we obtain

[TABLE]

Appendix E Weyl node pinned at a TRIM

In this section we consider a Weyl node near the Fermi level and pinned to a time-reversal invariant momentum (TRIM) by time-reversal symmetry. These nodes appear in chiral structures Chang et al. (2018) and their partners of opposite chirality are generally far away in energy and hence Pauli blocked. A single Weyl node can therefore account for the response of the whole crystal, if quadratic corrections to the dispersion are included. Due to time-reversal symmetry, $u_{i}=0$ and the only quadratic term allowed is proportional to the identity

[TABLE]

For simplicity we now restrict the example to chiral materials with point group $D_{3}$ such as TlTe2O6 Chang et al. (2018), where the only non-zero parameters are $v_{xx}=v_{yy}\equiv v_{\parallel}$, $v_{zz}\equiv v_{\perp}$ and $A_{xx}=A_{yy}\equiv A_{\parallel}$, $A_{zz}\equiv v_{\perp}$. This model is also applicable to the conduction band of elemental Te around the H point (not a TRIM) because the spin-dependent quadratic terms are negligible Ivchenko and Pikus (1975).

Here we computate the traces $\beta$ and $\gamma$ for the Hamiltonian in Eq. (99). The energies after the change of variables are given by $\epsilon_{k^{\prime}}=\pm k^{\prime}+ak^{\prime 2}+bk^{\prime 2}_{z}$ with $a=A_{\parallel}/v_{\parallel}^{2}$ and $b=A_{\perp}/v_{\perp}^{2}-A_{\parallel}/v_{\parallel}^{2}$. We assume that $a>0$ and $b>0$. The only modification from the quadratic term occurs in the Fermi functions, which are $f_{1}=\Theta(-k^{\prime}+ak^{\prime 2}+bk^{\prime 2}_{z}-\mu)$ and $f_{2}=\Theta(k^{\prime}+ak^{\prime 2}+bk^{\prime 2}_{z}-\mu)$.

[TABLE]

Defining $c_{\pm}=(\pm 2/\omega+a-4\mu/\omega^{2})/b$ and the functions $f_{\pm}=\Theta(c_{\pm}+1)\Theta(-c_{\pm})c_{\pm}+\Theta(c_{\pm})$ and changing variables to $x=\cos\theta$ we have

[TABLE]

To obtain $\gamma_{1}$ we use the Kramers-Kronig formula (68)

[TABLE]

and change variables to $y_{\pm}=|f_{\pm}(\omega)|^{1/2}$, $dy_{\pm}=\partial_{\omega^{\prime}}|f_{\pm}(\omega^{\prime})|^{1/2}d\omega^{\prime}$ to get

[TABLE]

where $\omega^{\prime 2}(y_{\pm})=\left(\frac{1\pm\sqrt{1+4\mu(a+by_{\pm})}}{a+by_{\pm}}\right)^{2}$. The integrals give

[TABLE]

where $A=\sqrt{1+4\mu a}$, $B=\sqrt{1+4\mu(a+b)}$, $g_{\pm}=\tfrac{4\mu}{\omega}\pm 1$.

We note that for $H_{\mathrm{TRIM}}$ $\gamma_{2}^{ab}$ was computed ab-initio for Te in Ref. Tsirkin et al. (2018) and the injection piece was computed in Ma et al. (2017).

E.1 Threefold and fourfold fermion

For a symmetric threefold fermion where the lowest bands (1,2) are completely filled, the matrix elements restrict the transitions to be zero unless $m=3$ and $n=2$. Evaluating the integrals as in the previous sections with $\Delta^{a}_{23}=-\partial_{k_{a}}\omega_{3}=-v_{F}k^{a}/k$, $\Omega^{b}_{3}=k^{b}/k^{3}$, $\Omega^{b}_{2}=0$ and Eq. (83) we get

[TABLE]

where $\beta_{0}=iC/(4\pi)=i\pi e^{3}/h$ is the universal circular photogalvanic quantum. A symmetric fourfold fermion is composed by two Weyls with the same chirality. The analytic results are given by

[TABLE]

Fig. 2 in the main text is computed by summing the contributions of a threefold with chemical potential $\mu_{3f}$ and a fourfold with chemical potential $\mu_{4f}$ of opposite chiralities.

Appendix F Details of the ab-initio calculation

TaAs has space group Lv et al. (2015b) $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 $\beta^{xy}=-\beta^{yx}$, and similar for $\gamma^{ab}$ Sturman and Fridkin (1992). 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) Koepernik and Eschrig (1999) within the generalized gradient approximation (GGA) Perdew et al. (1996). 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\times 200\times 200$ to $960\times 960\times 960$ Zhang et al. (2018b). Satisfactory convergence was achieved for a k-grid of size $400\times 400\times 400$. To compute $\gamma^{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. $\gamma^{ab}_{1}$ is obtained via Eq. (68) with the use of Kramers-Kronig transformations.