Divergent bulk photovoltaic effect in Weyl semimetals
Xu Yang Kenneth Burch Ying Ran
cond-mat.mes-hall



Department of Physics,
Boston College, Chestnut Hill, MA 02467, USA


Department of Physics,
Boston College, Chestnut Hill, MA 02467, USA


Department of Physics,
Boston College, Chestnut Hill, MA 02467, USA


 Weyl semimetals (WSM) have been discovered in time-reversal symmetric materials, featuring monopoles of Berry's curvature in momentum space. WSM have been distinguished between Type-I and II where the velocity tilting of the cone in the later ensures a finite area Fermi surface. To date it has not been clear whether the two types results in any qualitatively new phenomena. Here we focus on the shift-current response (σ_shift(ω)), a second order optical effect generating photocurrents. Surprisingly we find that up to an order unity constant, σ_shift(ω)∼e^3/h^21/ω in Type-II WSM, diverging in the low frequency ω→ 0 limit. This is in stark contrast to the vanishing behavior (σ_shift(ω)∝ω) in Type-I WSM. In addition, in both Type-I and Type-II WSM, a nonzero chemical potential μ relative to nodes leads to a large peak of shift-current response with a width ∼ |μ|/ħ and a height ∼e^3/h1/|μ|, the latter diverging in the low doping limit. We show that the origin of these divergences is the singular Berry's connections and the Pauli-blocking mechanism. Similar results hold for the real part of the second harmonic generation, a closely related nonlinear optical response.


Divergent bulk photovoltaic effect in Weyl semimetals
    Ying Ran
    December 30, 2023
=====================================================



Introduction— Recently, Weyl semimetals (WSM) have been discovered in many materials with strong spin-orbit coupling<cit.>. The topology of the electronic band structures gives rise to linear band touching points — Weyl nodes in momentum space, which are monopoles of the Berry's connection and related to the chiral anomaly in the context of high energy physics<cit.>. These topological semimetals have been shown to host various exotic properties such as surface Fermi arcs<cit.>, semi-quantized anomalous Hall effect<cit.>, angle-dependent negative magnetoresistance<cit.>, novel nonlinear optical effects<cit.>, and a recent observation of a Colossal Bulk PhotoVoltaic Effect (BPVE)<cit.>. The non-linear optical response has received increasing attention as a means to probe the Berry curvature of materials in general. This suggests non-linear optical effects can be used to distinguish between materials with different Fermi surface topologies, a question particularly relevant to WSM. Indeed, shortly after the discovery of the first Type-I WSM material in TaAs,<cit.> it was realized the tilt of velocity of the cone can be severe as to result in finite Fermi surfaces at all doping levels.<cit.> Nonetheless a clear distinguishing experimental consequence between these Type-II and their Type-I counterparts has yet to emerge. 

To linear order, the electrons near a Weyl node can be described by the effective two-by-two Hamiltonian: H^Weyl= (∑_ak^a  α^a-μ) σ_0+∑_a,b k^a β_abσ_b, where μ is the chemical potential, k⃗ is the momentum relative to the Weyl node, a/b=x,y,z, σ_b are the Pauli matrices and σ_0 is the identity matrix. After choosing the frame of the principal axes and a proper basis, the off-diagonal elements in β can be eliminated:

    H^Weyl=(ħk⃗·v⃗_t-μ) σ_0+ħ∑_a (k^a v^a  σ_a).

Here apart from the generally anisotropic velocity v⃗, another velocity v⃗_t describes the tilting of the bands and breaks the degeneracy between the bands at k⃗ and -k⃗, with the chemical potential (μ) usually nonzero.

It is known that either inversion symmetry or time-reversal symmetry needs to be broken in order to realize the Weyl nodes, and the experimentally confirmed Weyl semimetals have been overwhelmingly non-centrosymmetric and time-reversal symmetric, including the (Nb,Ta)(As,P) and W_1-xMo_xTe_2 series<cit.>. When the tilting velocity is large enough, which is realized in materials like W_1-xMo_xTe_2, the system can becomes a Type-II WSM<cit.>, since along some directions the two energy bands share the same sign of the velocity. Such a system must have finite size Fermi surfaces  (Fig.<ref>(d))[ Consequently it may be more accurate to call it a Type-II Weyl metal instead of a semimetal. However to be consistent with existing literature we follow the existing naming convention.]. The type of a Weyl node is determined by a dimensionless number W:

    W≡√((v_t^x/v^x)^2+(v_t^y/v^y)^2+(v_t^z/v^z)^2),

and a Type-I(Type-II) Weyl node is realized if W<1 (W>1). 

As we show later, the bulk photovoltaic effect<cit.> with its direct connection to Berry curvature,<cit.> offers a method to distinguish these two types. It has also attracted significant interest due to its potential applications in renewable energy generation<cit.> and fast photo-detectors<cit.>. The intrinsic contributions to the BPVE can be expressed as the second order nonlinear photocurrent response:

    𝐣^a=σ_2^abc(ω) E^b(ω)E^c(-ω),

where, a/b/c=x,y,z, 𝐣 is the DC electric current density, E⃗(t)=[E⃗(ω)e^-iω t+E⃗(-ω)e^iω t] is the electric field of the light. In order to have a nonzero σ_2, inversion symmetry needs to be broken, which happens to be also a condition to realize WSM. Generally speaking, σ_2^abc has both intraband and interband contributions. However, for time-reversal symmetric materials with linear polarized light, it turns out that this photocurrent response only has interband contributions, which has been coined the shift-current(σ_2^abc=σ_shift^abc) as it results from a change in the center of mass of the electrons upon optical excitation. Perturbation theory within the single-particle framework<cit.> gives:

    σ_shift^abc(ω)=2π e^3/ħ^2∫d^3k⃗/(2π)^3∑_n,mI_mn^abc[f_nm·δ(ω_mn-ω)],

where n,m label energy bands, ħω_mn(k⃗)≡ E_m(k⃗)-E_n(k⃗) is the energy difference between the two bands. f_nm(k⃗)=f_n(k⃗)-f_m(k⃗) is the difference of the Fermi-Dirac function between the two bands. The gauge invariant quantity I_mn^abc(k⃗)≡1/2·[r_mn^b r^c_nm;a+r_mn^c r_nm;a^b], where r^a_mn(k⃗)≡ i⟨ u_m(k⃗)|∂_k^a|u_n(k⃗)⟩ is nothing but the non-Abelian Berry's connection (with r^a_mn;b(k⃗) its generalized derivative: r^a_mn;b≡∂ r_mn^a/∂ k^b-i[A_m^b(k⃗)-A_n^b(k⃗)]r_mn^a(k⃗), and A_n^b(k⃗)≡ i⟨ u_n(k⃗)|∂_k^b|u_n(k⃗)⟩ the usual intraband Berry's connection).

It has been pointed out that the shift current response is related to the topology of the band structure<cit.>. Indeed the quantities responsible for σ_shift directly involves the Abelian and non-Abelian Berry's phases. Because Weyl nodes are monopoles of the Abelian Berry's connection, these quantities are expected to be diverging and singular near Weyl nodes, which motivates us to carefully study the resulting nonlinear optical effects. 

Eq.(<ref>) has a familiar form of the Fermi's golden rule. Indeed, considering the case of linearly polarized light along the b-direction, one has:

    I^abb_mn=|r_nm^b|^2 R_nm,b^a,

where the gauge invariant real space displacement R⃗_nm,b(k⃗) is the so-called “shift-vector”, defined as:

    R_nm,b^a(k⃗)≡ -∂[r_nm^b]/∂ k^a+A^a_n(k⃗)-A^a_m(k⃗),

If one interprets Eq.(<ref>,<ref>) as the Fermi's golden rule, |r_nm^b|^2 is just the matrix-element factor for the optical absorption. Thus, eR⃗_nm,b(k⃗) should be viewed as the dipole moment induced the photo-excited particle-hole pair, giving rise to a rate of change in polarization, i.e., DC photocurrent. 


Motivated by photogalvanic applications<cit.>, previously the shift-current response has been mainly discussed in the context of insulators<cit.>, where the low temperature σ_shift vanishes when ħω is below the band gap.  It is convenient to introduce the optical joint density of states including the factor f_nm responsible for the Pauli-blocking effect:

    JDOS(ω)≡∫d^3k⃗/(2π)^3∑_m,n[f_nm·δ(ω_mn-ω)]

For most materials, I^abc_mn(k⃗) is a smooth function of momentum, and according to Eq.(<ref>) the response σ_shift is essentially proportional to JDOS. It has been proposed that engineering JDOS in semiconductors may be a route to optimize the BPVE<cit.>. Although the JDOS∝ω^2 is small due to the linear dispersion in WSM, recent second harmonic<cit.> and photocurrent experiments<cit.> show that these materials host large nonlinear optical effects at least in the infrared regimes. In this paper we show that the WSM actually feature divergent nonlinear optical responses in the low frequency regime due to the singular Berry's phases near the Weyl nodes.






Results: Our main results are summarized in Fig. <ref> and Fig. <ref>. For simplicity let us start with the zero doping case μ=0. In the absence of extrinsic scattering processes, simple dimensional analysis shows that up to fundamental constants e^3/h, σ_shift(ω) is proportional to the inverse of an energy scale, which could involve either an intrinsic energy scale t of the material, the temperature k_B T, or the photon energy ħω. 

At zero doping and zero temperature, the linear Weyl equation Eq.(<ref>) does not contain an energy scale, therefore σ_shift(ω) must be either ∼e^3/h^21/ω or vanishing in this linear approximation. The former divergent case is exactly realized in Type-II WSM (Fig.<ref>(f)). In Type-I WSM, however, we will show that zero doping and zero temperature leads to σ_shift(ω)∼e^3/hħω/t^2 (Fig.<ref>(c)), where t is an intrinsic energy scale due to the band-bending (i.e., deviation from linear dispersion) —typically ∼0.1-0.2eV.

Next, we discuss the reason for the drastic difference between Type-I and Type-II WSM, as well as the effect of doping. A careful evaluation of σ_shift in Eq.(<ref>) based on the linear Weyl equation Eq.(<ref>) shows that, at T=0, a single Weyl node contributes (At T≠0, a general and more complicated analytic result is shown in supplemental material Eq.(<ref>,<ref>)):

    σ^abc_shift   =πχ/8W^2·∑_d[ϵ_dcav_t^bv_t^d/(v^d)^2+(b↔ c)]·σ^linear_node(ω)
              
    σ^linear_node   ≡ (cos[θ_1]sin[θ_1]^2-cos[θ_2]sin[θ_2]^2)·e^3/h^2ω,
    θ_1    ≡[2μ/ħω+1/W], θ_2≡[2μ/ħω-1/W].

Up to an order unity constant σ^abc_shift(ω) is determined by  σ^linear_node(ω). Here W is defined in Eq.(<ref>), ϵ_abc is the Levi-Civita antisymmetric tensor, and we have chosen the convention ω>0. χ≡[v^xv^yv^z]=±1 is the chirality of the Weyl node (i.e., monopole charge). Note that we have chosen the frame x,y,z to be right-handed. We have also defined a function [s]≡0 if s⩾1, [s]≡π if s⩽-1, and [s]≡[s] if -1<s<1. 



Interestingly, due to the chirality factor in this result, we know that this 1/ω divergent term is absent in Dirac semimetals, because in those systems each Dirac node can be viewed as two Weyl nodes with opposite chiralities sharing the same set of μ, v⃗_t and opposite v⃗. It can also be immediately seen from Eq. (<ref>) that the 1/ω divergence is absent for  σ^aaa_shift due to the Levi-Civita antisymmetric tensor.

The physical meaning of the two angles θ_1,θ_2 is the following. For the moment it is convenient to re-scale the momentum in x,y,z directions so that v^x=v^y=v^z=v. The energy conservation δ(ω_mn-ω) constrains our consideration on a sphere in the re-scaled momentum space, whose radius equals 2ω/v. The re-scaled v⃗_t can be used to define a special axis to set up a spherical coordinate system. Generally speaking, Pauli blocking takes over in certain solid angle regions. Namely the factor f_nm further constrains the sphere into a region between the polar angles θ_1 and θ_2 (See Fig.<ref>(e) for an illustration). 

In Type-I WSM (W<1) at zero doping, the whole sphere contributes, consistent with θ_1=0 and θ_2=π (Fig.<ref>(b)). But this is exactly a situation when the θ-dependent factor in Eq.(<ref>), and thus the 1/ω term in σ_shift, vanishes. (Note that, in this case, a careful analysis including band-bending effects shows that even the constant order vanishes, leaving the next order σ_shift∝ω in the low frequency limit. See supplemental information.) But in Type-II WSM (W>1) at zero doping, only the part of the sphere between θ_1,2=
[±1/W] contributes (Fig.<ref>(e)), leading to the 1/ω divergent response.

This results from the strongly angle dependent diverging Berry connection in WSM, despite the integrand I_mn^abc in Eq.(<ref>) scaling as 1/ω^3. The full angular average over the 4π solid angle would annihilate the 1/ω term in σ_shift. However, Pauli-blocking could take over in certain angular regions, removing the net cancellation and retaining the divergent term.

In particular, at T=0, in Type-I WSM (W<1) with finite μ, this 1/ω term survives only over a frequency range: 2|μ|/1+W⩽ħω⩽2|μ|/1-W. (Thus in the extreme case when v_t=0 and W=0, this term vanishes.) In fact, it is straightforward to show that, σ^linear_node would also change sign exactly at ħω=2|μ|. Altogether this leads to a large peak in the σ_shift(ω) with a width ∼ |μ|/ħ and a height ∼e^3/h·1/|μ|, featuring a characteristic sign-changing line-shape (Fig.<ref> top). Similarly a finite μ truncates the 1/ω divergence in Type-II WSM and leads to a sign-changing large peak when ħω∼μ (Fig.<ref> bottom).

As shown in Fig.<ref>, the line-shapes of these peaks are quasi-universal: they only depend on the dimensionless number W. Different values of doping μ only re-scale the peak's frequency range by a μ factor and its height by a 1/μ factor. This means that one may use μ to control the frequency range of the peak to engineer tunable frequency-sensitive photo-electric devices. 


Temperature and impurities— At finite temperatures, the Fermi-Dirac distribution would smear out and truncate the divergences when ħω≲ k_B T (which is confirmed in our tight-binding model calculations Fig.<ref>). But at low temperatures k_B T≪ħω these divergences are not significantly modified. (See supplemental material Eq.(<ref>,<ref>) for a general analytic form of the σ_shift within linear-dispersion approximation.) However even at zero temperature, impurities give rise to scattering, while finite temperature enables other scattering mechanisms (e.g. electron-phonon). These scatterings, which can be phenomenologically characterized by a scattering time τ, have been ignored so far. Namely, even at low temperatures k_B T≪ħω our result Eq.(<ref>) holds only in the long scattering time ωτ≫1 limit with the divergences truncated when ω≲1/τ. For instance, previous experiments report that τ in the Type-I WSM TaAs is of the order of a pico-second<cit.>. This suggests that our predicted striking response can be observed in the Terahertz or higher frequency regimes.

A conceptually interesting situation occurs when a finite temperature is introduced in the μ=0 Type-I WSM, where divergence is absent at T=0 due to angular cancellation discussed before. However when T≠ 0, the thermally excited particle-hole pairs partially play the role of doping, and there is no reason for a full angular cancellation. In this case, simple dimensional analysis leads to striking results: we expect a temperature-induced peak of the intrinsic σ_shift when ħω∼ k_BT, whose height ∼e^3/hk_B T diverging in the low T limit. (see supplemental material Eq.(<ref>) and Fig.<ref>) This is observed in our tight-binding model calculations (see Fig.<ref>(a)).

The size of the effects —It is interesting to estimate the size of effects in the divergent regimes; e.g., when ħω∼μ for Type-I WSM or ħω≳μ in Type-II WSM, in the presence of a temperature k_BT≪ħω. Previously shift-current responses have not be much studied in the low frequency (e.g. Terahertz) regimes. Note that even in a generic multiband metal, σ_shift is expected to vanish at zero temperature in the low frequency regimes. This is simply because the energy conservation δ(ω_mn-ω) and f_nm in JDOS constrain both the valence band and the conduction band at the Fermi level when ω→ 0, which would not occur due to band-repulsions.

Plugging in ω=1THz, the estimated size of σ_shift∼e^3/h^21/ω  in the divergent regimes is ∼ 0.01 A/V^2, several orders of magnitudes larger than known reported values in visible or infrared regimes<cit.>







Tight-binding model— To concretely illustrate the predicted responses we compute the σ_shift tensor using a minimal time-reversal symmetric 4-band tight-binding model featuring 4 Weyl nodes:<cit.>

    H^TB=t[(2.5 - cos k^x - cos k^y - cos k^z) τ_x +sin k^y τ_y 
       + 0.5 cos k^y s_xτ_y + sin k^z s_xτ_z + (ξcos k^x - μ_0)]-μ,

where s_a(spin) and τ_a(orbital) are two sets of Pauli matrices. t is an overall energy scale which can be 0.1-0.5eV, to be broadly consistent with the relevant energy scales in existing Weyl materials. Symmetries in this model include: time-reversal is_yK (K is conjugation), x→ -x mirror s_x, and y→ -y mirror s_yτ_x. The four Weyl nodes are located at (k^x=± 0.920,k^y=± 0.464,k^z=0), with an energy μ_0t and μ_0=0.606ξ. A positive parameter ξ controls the tilting velocity v⃗_t: for ξ<1(ξ>1), Type-I(Type-II) WSM is realized. We choose ξ=0.5 (ξ=1.5), corresponding to W=0.539(W=1.616), as the representative for Type-I (Type-II) WSM.


We numerically compute σ^abc_shift, using the full formula Eq.(<ref>) without resorting to the linear approximation. The results of a particular component σ^zxx_shift are plotted in Fig.<ref> for various doping levels and temperatures, which are fully consistent with previous discussions.


Second harmonic generation— Another directly related second order nonlinear optical response is the second harmonic generation (SHG), in which light at frequency ω drives current at frequency 2ω<cit.>. Defining SHG response tensor 𝐣^a(2ω)=σ_SHG^abc(ω,2ω)E^b(ω)E^c(ω), it is known that the real part [σ_SHG] is given by the interband contribution, and within the two-band approximation and linear-node approximation, we have

    [σ^abc_SHG(ω,2ω)]=-3/2σ^abc_shift(ω)
 
where σ^abc_shift(ω) given by Eq. (<ref>) (see supplemental information). [Previously another identity is known within general two-band approximation: [σ^aaa_SHG(ω,2ω)]=-σ^aaa_shift(ω)+1/2σ^aaa_shift(2ω) <cit.>, which naively is inconsistent with Eq.(<ref>). However within linear-node approximation we have shown that σ^aaa_shift(ω)=0, which is consistent with both identities.] [Interestingly, the intraband contribution of SHG is responsible for [σ_SHG], which has been pointed out to be a Fermi surface Berry's curvature effect and contains a 1/ω term<cit.>]. For example, similar to σ_shift, the zero temperature [σ_SHG] features a similar 1/ω divergence in Type-II WSM when μ=0, and large peak behaviors in both Type-I and Type-II WSM when μ≠0. And even at μ=0, a finite temperature induces a large peak in [σ_SHG] in Type-I WSM with a width ∼ k_B T and a height ∼e^3/h k_B T.



Possible applications— In this paper we report the 1/ω diverging DC photovoltaic effect in the low frequency regime in WSM, due to the combination of the diverging Berry's curvature and Pauli-blocking effect. Fast Terahertz photon detection has been a long standing challenge. The reported large photovoltaic effect may be useful for this purpose. In addition, the doping-induced large peak regimes of σ_shift in WSM, whose frequency ranges are controlled by μ, may be useful as a tunable frequency-sensitive probe for far-infrared or Terahertz photons, i.e., a spectrum analyzer. 

XY and YR acknowledge support from the National Science Foundation under Grant No. DMR-1151440. KSB is grateful for support from the National Science Foundation, Grant No. DMR-1709987.




apsrev





Supplemental Material:
 Divergent bulk photovoltaic effect in Weyl semimetals



Xu Yang, Kenneth Burch, and Ying Ran





§ SHIFT CURRENT IN TYPE-I WEYL SEMI-METAL

In this section we shall prove that for a generic type-I Weyl semi-metal with Fermi level at the Weyl nodes, the leading term in the shift current tensor will be proportional to ω when ω→ 0.

The low energy physics can be captured by the following generic 2-band Hamiltonian with chemical potential μ set to zero

    H=f_0σ_0+∑_i=x,y,zf_iσ_i,

where f_0,f_i are functions of k and f_i=0 when k=0. The eigenvalues are E_c=f_0+ϵ, E_v=f_0-ϵ with ϵ=√(∑_i f_i^2). Since the tilting will not affect the shift current tensor, we will set f_0 to be zero.

The shift current tensor for the 2-band model is obtained by doing the following integral<cit.>

    σ_shift^abc(ω)=2π e^3/ħ^2∫d^3k⃗/(2π)^3I_cv^abc[f_cv·δ(ω_cv-ω)],

where I^abc_cv has the following explicitly gauge-invariant expression

    I^abc_cv=∑_i,j,m[1/8ϵ^3(f_mf_i,bf_j,ac-f_i,bf_j,af_mϵ_,c/ϵ)ϵ_ijm+(b↔ c)],

where ϵ_ijm is the Levi-Civita symbol.

Now let's prove that terms proportional to 1/ω (denoted as σ^(-1)(ω) below) and terms independent of ω (denoted as σ^(0)(ω) below) in σ^abc vanish in the low frequency limit ω→ 0. A simple order of estimate with k∼ω when ω→ 0 tells us that only the k-linear terms and k-quadratic terms in f_i will contribute to σ^(-1)(ω) and σ^(0)(ω). Therefore it suffices to consider a linear node plus some quadratic corrections. It is always possible to choose the following form of f_i by an affine transformation which does not affect the integral (ħ v_F is set to be 1 throughout this section)

    f_i=k_i+α_ijlk_jk_l,

where α_ijl is a rank-3 tensor symmetric respect to the interchange j↔ l and repeated indices are summed over.

In the spherical coordinate system, we have

    k_x=r sin[θ]cos[ϕ],k_y=r sin[θ]sin[ϕ],k_z=r cos[θ].

And the unit vector k̂=(sin[θ]cos[ϕ],sin[θ]sin[ϕ],cos[θ]) will be used below to simplify the notation.

When ω is small, k_i∼ r∼ω and we can expand ϵ in powers of r

    ϵ=r+r^2f(θ,ϕ)+𝒪(r^3),

where f(θ,ϕ) is a function of α_ijl,θ,ϕ. It is easy to obtain the explicit form of f(θ,ϕ), but for our purpose we only need the fact that f(θ,ϕ)=-f(π-θ,π+ϕ) since it comes from the angular dependence of terms involving 3 k_i's.

In the same spirit we can expand I^abc in powers of r

    I^abc=1/8ϵ^3[I^(0)_abc+I^(1)_abc+𝒪(r^2)], with 
       I^(0)=-ϵ_bamk_mk_c/ϵ^2+(b↔ c)=-ϵ_bamr^2k̂_mk̂_c/ϵ^2+(b↔ c),
       I^(1)=2ϵ_bjmα_jack_m+h^abc(k_x,k_y,k_z)/ϵ^2+(b↔ c)
    =2rϵ_bjmα_jack̂_m+r^3h^abc(k̂_x,k̂_y,k̂_z)/ϵ^2+(b↔ c)

where h^abc(k_x,k_y,k_z) is a homogeneous polynomial of degree 3, whose explicit form, for our purpose, is not important.


The integral then becomes

    σ^abc(ω)=const*∫ dΩ(∫ r^2drI^abc(r,θ,ϕ)δ(2ϵ(r,θ,ϕ)-ω)),

where dΩ=sin[θ]dθ dϕ.



When ω→ 0, we know that for any fixed θ,ϕ there is only one solution r(θ,ϕ) to the equation 2ϵ(r(θ,ϕ),θ,ϕ)=ω, from which we can solve r as a function of θ and ϕ. In fact, we can expand r(θ,ϕ) in powers of ω with f(θ,ϕ) defined before

    r(θ,ϕ)=ω/2-ω^2/4f(θ,ϕ)+𝒪(ω^3).



The integral can then be written as

    σ^abc(ω)=const*∫ dΩ1/ω^3 r(θ,ϕ)^2/∂_rϵ|_r=r(θ,ϕ)[I^(0)_abc+I^(1)_abc+𝒪(ω^2)]
       =const*∫ dΩ1/ω^3[ r(θ,ϕ)^2-2r(θ,ϕ)^3f(θ,ϕ)+𝒪(ω^4)][I^(0)_abc+I^(1)_abc+𝒪(ω^2)]
       =const*∫ dΩ1/ω[ 1-2ω f(θ,ϕ)+𝒪(ω^2)][I^(0)_abc+I^(1)_abc+𝒪(ω^2)]


Terms that are of order 1/ω comes from the integral over I^(0)_abc with r(θ,ϕ) set to ω/2 in Eq.(<ref>)

    σ^(-1)(ω)=const*∫ dΩ1/ω(I^(0)_abc)=const*∫ dΩ1/ω[-ϵ_bam(k̂)_m(k̂)_c+(b↔ c)]=const*(-ϵ_bac-ϵ_cab)=0,

where we have used the fact that ∫ (k̂)_m(k̂)_ndΩ=4π/3δ_m,n.

Terms that are independent of ω comes from the following 3 integrals:

    σ^(0)_I=const*∫ dΩ1/ω(I^(0)_abc)=const*∫ dΩ f(θ,ϕ)[-ϵ_bamk̂_mk̂_c+(b ↔ c)],


    σ^(0)_II=const*∫ dΩ1/ω(I^(1)_abc)=const*∫ dΩ [ϵ_bjmα_jack̂_m+h^abc(k̂_x,k̂_y,k̂_z)/2+(b↔ c)]

and 

    σ^(0)_III=const*∫ dΩ f(θ,ϕ)(I^(0)_abc)=const*∫ dΩ f(θ,ϕ)[-ϵ_bam(k̂)_m(k̂)_c+(b↔ c)].



It is easy to see that under (θ,ϕ)→ (π-θ,π+ϕ), the integrand of Eq.(<ref>),(<ref>),(<ref>) all change sign, therefore they are all equal to zero.

In conclusion, we have analytically shown that in the low frequency limit ω→ 0, term that diverges as 1/ω and term that is independent of ω in the shift current tensor σ^abc vanish. Therefore the leading term in σ^abc will be proportional to ω.




§ ANALYTICAL FORMULA FOR THE SHIFT CURRENT IN WEYL SEMI-METAL WITH TILTING AND DOPING IN LOW-FREQUENCY LIMIT

In this section we will obtain an analytical formula of the shift-current tensor for a Weyl node with both tilting and nonzero chemical potential within the linear approximation at zero temperature. This formula captures the physics of type-I and type-II Weyl semimetals with or without doping in a unified fashion.

Let's consider the following generic Hamiltonian

    H=(ħk⃗·v⃗^t-μ)σ_0+∑_i=x,y,zħ k_iv_iσ_i.


We will choose the following parameterization v_ik_i=rO_ijk̂_j, where O_ij is an orthogonal matrix and k̂=(sin[θ]cos[ϕ],sin[θ]sin[ϕ],cos[θ]). The orghogonal matrix O is chosen such that

    (v^t_x/v_x,v^t_y/v_y,v^t_z/v_z)· O=(0,0,W),
 where W=√((v^t_x/v_x)^2+(v^t_y/v_y)^2+(v^t_z/v_z)^2). In fact the third column of O is fully determined
    (O_13,O_23,O_33)=1/W(v^t_x/v_x,v^t_y/v_y,v^t_z/v_z).


In this new coordinate system, we have

    σ_shift^abc(ω)=2π e^3/ħ^2∫d^3k⃗/(2π)^3I_12^abc[f_21·δ(ω_12-ω)]= e^3sign[v_xv_yv_z]/4π^2ħ^2v_xv_yv_z∫ r^2dr dΩ I_12^abcf_21·δ(2r-ω).

where I^abc_12 is given by

    -v_av_bv_c/8r^3[O_mjk̂_jO_clk̂_lϵ_bam+(b↔ c)]

and f_21 is given by

    f_21=1/e^β (ħ Wrcos[θ]-ħ r-μ)+1-1/e^β (ħ Wrcos[θ]+ħ r-μ)+1,

which is a function of r and θ.

Therefore we have

    σ_shift^abc(ω)=- e^3sign(v_xv_yv_z)v_av_bv_c/8ω h^2v_xv_yv_z∫sin[θ]dθ dϕ f_21(ω/2,θ) [O_mjk̂_jO_clk̂_lϵ_bam+(b↔ c)].


The angular integration can be easily done: 

    ∫ dΩ f_21(ω/2,θ)k̂_mk̂_n=∫ dΩ f_21(ω/2,θ)(k̂_x^2 δ_m,n+(k̂_z^2-k̂_x^2)δ_m,nδ_m,3)=c·δ_m,n-π· g(ϖ,z,W)δ_m,nδ_m,3,

where c is a constant, z=e^βμ is the fugacity and ϖ=βħω/2. The isotropic part c·δ_m,n does not contribute to the shift current tensor as in the case of type-I Weyl semimetal without doping. And the anisotropic part can be evaluated by first integrating over ϕ and then integrating over x=cos[θ] to yield

    g(ϖ,z,W)=-∫^1_-1(3x^2-1)f_21(ω/2,x)dx=2/Wϖln[1+z^-1e^(1-W)ϖ/1+z^-1e^(1+W)ϖ·1+z^-1e^(W-1)ϖ/1+z^-1e^-(1+W)ϖ]
       -∑_n=2,36/W^nϖ^n[Li_n(-z^-1e^(1-W)ϖ)+Li_n(-z^-1e^(1+W)ϖ)-Li_n(-z^-1e^-(1+W)ϖ)-Li_n(-z^-1e^(-1+W)ϖ)],

where Li_n(x) is the polylogarithm of order n.

Therefore after inserting Eq.(<ref>) and Eq.(<ref>) into Eq.(<ref>) we have the following result

    σ^abc_shift(ω)   =πsign(v_xv_yv_z)/8W^2·[ϵ_dcav^t_bv^t_d/(v_d)^2+(b↔ c)]·σ^linear_node(ω)
       σ^linear_node(ω)≡[g(ϖ,z,W)·e^3/h^2ω]





Therefore up to an order 1 constant, the shift-current response tensor is determined by the function σ^linear_node. Let's discuss two simple limits of Eq. (<ref>). First, when μ=0, the function g only depends on the dimensionless variable ϖ, which gives us the scaling form 
    σ_shift(ω)|_β=b·σ_shift(bω)|_β/b
 
We plot this scaling behavior at finite temperatures and μ=0 for a few representative values of W in Fig.<ref>. In particular for Type-I WSM at μ=0, T≠ 0 induces a peak of σ_shift whose width∝ k_BT and height∝ 1/k_B T. This behavior is also observed in Fig.3(a) in our tight-binding model calculations.

Next, when β→∞, we can obtain a simpler form of g(ϖ,z,W). In fact, an easier way is to replace the Fermi-Dirac distribution function by the Heaviside step function: 
    f_21=Θ(r-k⃗·v⃗^t+μ/ħ)-Θ(-r-k⃗·v⃗^t+μ/ħ)=Θ(r(1-W cos[θ])+μ/ħ)-Θ(-r(1+Wcos[θ])+μ/ħ).


The Heaviside step function together with the δ function constraint the integration region of θ, therefore we will introduce the following two θ angles to characterize the upper and lower limit of the integration

    θ_1 ≡[2μ/ħω+1/W], θ_2≡[2μ/ħω-1/W]

where the function [s] is defined in the following way: [s]≡0 if s⩾1, [s]≡π if s⩽-1, and [s]≡[s] if -1<s<1.

After integration we have

    σ^linear_node(ϖ,z,W)|_T→ 0=(cos[θ_1]sin[θ_1]^2-cos[θ_2]sin[θ_2]^2)e^3/h^2ω,

which gives us the shift-current tensor at zero-temperature:

    σ_shift^abc(ω)=πsign(v_xv_yv_z)/8W^2·[ϵ_dcav^t_bv^t_d/(v_d)^2+(b↔ c)]· (cos[θ_1]sin[θ_1]^2-cos[θ_2]sin[θ_2]^2)e^3/h^2ω







§ ANALYTICAL FORMULA FOR THE SECOND-HARMONIC-GENERATION IN WEYL SEMI-METAL WITH TILTING AND DOPING IN LOW-FREQUENCY LIMIT

The second-harmonic-generation (SHG) response tensor is defined via 𝐣^a(2ω)=σ_SHG^abc(ω,2ω)E^b(ω)E^c(ω).

From the standard time-dependent perturbation theory<cit.>, we have the following expression for the real part of the SHG response tensor

    Re[σ_SHG^abc(ω,2ω)]=iπ e^3/2ħ^2ω^2∑_m,n,p∫_BZdk^3/(2π)^3[(v^a_mnw^bc_nm+2v^a_mn{v^b_npv^c_pm}/ω_mp+ω_np)f_mnδ(2ω-ω_nm)
       +(w^ab_mnv^c_nm+w^ac_mnv^b_nm)f_mnδ(ω-ω_nm)+v^a_mn{v^b_npv^c_pm}/ω_pm+ω_pn(f_mpδ(ω-ω_pm)-f_npδ(ω-ω_np))],

where v^i_mn≡1/ħ⟨m|∂_k_iH|n⟩, w^ij_mn≡1/ħ^2⟨m|∂_k_i∂_k_jH|n⟩ and {v^b_npv^c_pm}=v^b_npv^c_pm+v^c_npv^b_pm. Note that since we are dealing with generic tight-binding models, a careful derivation following Ref.<cit.> yields the extra w^ij_mn terms which are absent in the literatures listed above.

For a 2-band model within the linear approximation, Eq. (<ref>) can be simplified to be

    Re[σ_SHG^abc(ω,2ω)]=iπ e^3/2ħ^2∫_BZdk^3/(2π)^3v^a_vc{v^b_cvΔ^c_cv}/ω_cv^3f_cv[8δ(2ω-ω_cv)-δ(ω-ω_cv)],

where Δ^i_nm≡ v^i_nn-v^i_mm and v/c corresponds to valence/conduction bands respectively.

In the linear approximation, we have the following identity

    r^b_mnr^c_nm;a+r^c_mnr^b_nm;a=-v^a_nm{v^b_mnΔ^c_mn}/ω_mn^3-Δ^a_mn{v^b_mnv^c_nm}/ω_mn^3.


Therefore we can rewrite Eq. (<ref>) as follows

    Re[σ_SHG^abc(ω,2ω)]=iπ e^3/2ħ^2∫_BZdk^3/(2π)^3f_vc(r^b_cvr^c_vc;a+r^c_cvr^b_vc;a)[8δ(2ω-ω_cv)-δ(ω-ω_cv)],

where we have used the fact that the integration over Δ^a_cv{v^b_cvv^c_vc}/ω_cv^3 vanishes due to time-reversal symmetry.

Eq. (<ref>) assumes a very similar form to the expression of shift-current. In fact, it's easy to see that within the linear approximation we have
Re[σ_SHG^abc(ω,2ω)]=-4σ^abc_shift(2ω)+1/2σ^abc_shift(ω)=-3/2σ^abc_shift(ω), where σ^abc_shift(ω) is given by Eq. (<ref>). As a consequence, all the discussions of leading-order terms in 
    Re[σ_SHG^abc(ω,2ω)]
 with the presence of tilting and doping at non-zero temperature naturally follow that of σ^abc_shift(ω) in the last section.







