Photogalvanic response in multi-Weyl semimetals

We investigate the dependence of the photogalvanic response of a multi-Weyl semimetal on its topological charge, tilt, and chemical potential. We derive analytical expressions for the shift and injection conductivities for tilted charge-$n$ Weyl points $(n=1,2,3)$ using a low energy two-band effective Hamiltonian. For double-Weyl semimetals, we also compute the response from two-band and four-band tight-binding models with broken time-reversal symmetry to study the effect of band bending and the contributions from higher bands. We find a significant deviation in the responses obtained from the effective low-energy continuum model and more realistic four-band continuum and tight-binding models. We analyze several different limits of these models. We describe the nature of the deviations and provide estimates of their dependence on the frequency and other model parameters. Our analysis provides a simple explanation for the first-principle calculation based frequency dependence of the injection current in SrSi$_2$. Additionally, we find interesting parameter regimes where the frequency dependence of the non-linear optical response can be directly used to probe the type-I/type-II nature of the Weyl cone. We obtain analytical results for the charge-4 Weyl semimetal by reducing the original problem involving a triple $k$-space integral to one with only a double integral. This simplification allows us to extract all relevant information about the nature of its second-order dc response and the precise condition for observing circular photogalvanic effect quantization. The semi-analytical approach presented here can also be extended to a systematic study of second harmonic generation and first-order optical conductivity in charge-4 Weyl semimetals.


I. INTRODUCTION
The quantum geometry (QG) of Bloch wavefunctions can significantly influence the electronic properties and response functions of a material [1].The quantum anomalous Hall effect in the absence of a magnetic field is the seminal example of such QG effects, originating in this case from the Chern number [2].More generally, the anomalous contribution from band topology can overcome limitations in non-topological systems on many physical properties like superfluid weight [3], exciton stability [4], transport coefficients [5], and optical responses [6,7].The bulk photovoltaic effect (BPVE) is one such effect where quantum geometry contributions have been shown to be of immense importance [8,9].
The BPVE is a second-order optical response where a DC current is produced in response to an AC electric field.It has been shown that in many noncentrosymmetric materials, the non-trivial structure of Bloch wavefunctions engenders a BPVE without creating any macroscopic electric field or carrier concentration gradient in the sample [10].This allows one to overcome the Shockley-Queisser limit [11] present in traditional p-n junctions.
Based on the mechanism of generation, the bulk photovoltaic effects can be divided into shift and injection currents [8].The shift current results from the realspace shift in the electron wavepacket due to interband photoexcitation [12], and the injection current is * raj.a@northeastern.educaused by change in electron velocity upon inter-band transition [8].The properties of these responses are determined by the polarization of light, and presence of time-reversal and space-inversion symmetries [13].The shift current response occurs for linearly-polarized light even when time-reversal symmetry is present.On the other hand, the injection current requires circularly polarized light and is also known as the circular photogalvanic effect (CPGE).However, when time-reversal symmetry is broken, both the shift current and injection current can occur for circularly and linearly polarized light, respectively [14].
These mechanisms for a BPVE are intimately related to the quantum geometry of the electronic wavefunction [15,16], and thus are proving to be reliable tools to probe and utilize the band topology [6,7].The bulk photovoltaic effects in Weyl semimetals have attracted enormous research interest as they provide a mechanism to generate photocurrents in the infrared and THz regime [17,18].It was shown in the seminal work [19], that the CPGE contribution from a Weyl node would exhibit quantization proportional to the charge of the Weyl node.Following these theory works, CPGE was measured in many different Weyl semimetals including TaAs [20], RhSi [21], and TaIrTe 4 [22] which showed interesting helicity-dependent behavior arising from the chirality of Weyl nodes.These experimental works also highlighted the importance of using more realistic models: the tilt and higher bands were shown to play an important role in determining the CPGE response [23].
Materials which can host Weyl nodes with Berry monopole charge higher than two are not known but triple-WSM can be possibly obtained from cubic Dirac semimetals [34] by applying a magnetic field or by Floquet engineering [35].Another interesting feature of these multi-Weyl semimetals is that the dispersion around Weyl node is no longer linear in all directions but becomes quadratic (cubic) for two directions for charge two (three).This leads to a strong anisotropy in the velocity matrix and also modifies the density of states which has been known to affect the transport coefficients [36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51] and linear optical responses [52][53][54][55][56][57] of multi-Weyl semimetals.These unusual properties of multi-Weyl semimetals are also believed to significantly influence the second-order optical responses, such as the BPVE and second-harmonic generation.A deeper understanding of how different properties of these MWSMs affect the shift current and injection current can possibly lead to a mechanism to probe the topological charge of Weyl semimetals.
Most theoretical works on the BPVE employ effective two-band low-energy Hamiltonians.These models have proven quite useful for general predictions like the quantization of the injection current conductivity, but the experimental signatures are often complicated by the discrepancy between effective low-energy models and real electronic band structure where the band curvatures and higher energy bands start to play an important role.As a result, the predictions of the continuum model usually agree only in a small energy window.This necessitates the need to analyze the role of different model parameters and understand the frequency behavior of Weyl semimetals in different regimes away from this small energy window.
In our work, we first provide a complete analytical solution to the two-band charge-n low energy Hamiltonian along with an analysis of its important features, including CPGE quantization.These analytical expressions elucidate the role of tilt and non-linear dispersion on different components of the shift and injection current conductivities in multi-Weyl semimetals.We also numerically evaluate the response in tight-binding models and observe a significant deviation in some components of second-order conductivity which highlight the importance of band curvature.
For multi-Weyl semimetals, the validity of two-band models becomes further restricted.Double Weyl nodes are obtained when two charge-1 Weyl nodes are pinned to a high-symmetry point and two of the four bands are gapped out by some symmetry allowed perturbations.As a result, even if the effective two-band picture is valid for each charge-1 Weyl node in a given energy range, it might not be valid for the double-Weyl node if the perturbation is not strong enough to push the other two bands out of that energy window.This type of scenario occurs in the charge-2 WSM SrSi 2 [29] where the two charge-1 Weyl nodes are gapped out by a spin-orbit coupling resulting in a very small gap between the bands hosting a double-Weyl node and the higher energy bands.
Inspired by the band structure of SrSi 2 , we also consider a four-band continuum model and find a significant deviation from the two-band continuum model.We find that the CPGE quantization is destroyed and instead a very different behavior is observed at small frequencies.In the particular case of the fourband model, we find two opposite limits in the parameter space with good and poor agreement.We notice that the agreement is better when the perturbation induced gap is large.Our analysis provides a simple-explanation for the results from first-principle calculations in Ref. [58] where quantization is observed only above a certain cutoff frequency.We attribute this discrepancy to the contribution from higher bands.
Finally, we also investigate the charge-4 case by using a two-band effective low energy model and a tightbinding model.We derive semi-analytical expressions for different components of the shift and injection current conductivities.Most importantly, we obtain the analytical limits for the frequency window where CPGE quantization can be observed.
Our paper is organized as follows.In Sec.II, we provide a brief introduction to the shift and injection current conductivities along with the symmetry requirements to observe their effects.In Sec.III, we derive expressions for different components of these second-order conductivity tensors by considering an effective two-band low-energy Hamiltonian for a Weyl node with arbitrary charge n.We also include a finite tilt in the z-direction in our analysis and systematically study how tilt affects these different components at different chemical potentials and frequencies.In Sec.IV, we focus on double-Weyl semimetals and consider two different models.First, we compare different conductivities for a two-band tightbinding model and an effective low-energy Hamiltonian.Next, we consider a four-band model inspired by the SrSi 2 band structure around its double Weyl node and study the second-order conductivities in different limits.In Sec.V, we derive the joint density of states (JDOS), and the shift and injection current conductivity expressions for a charge-4 model.In Sec.VI, we discuss the implications of our results.

II. PHOTOGALVANIC RESPONSE
In materials lacking inversion symmetry, the photogalvanic effect (PGE) refers to the generation of directed photocurrent as a second-order response to an external time-varying electromagnetic field.For light of frequency ω (and wavelength much larger than the sample size so the electric field has uniform amplitude), the second-order dc response is given by, where the second-order conductivity σ abc (ω) can be divided into a shift current conductivity, σ abc shift and an injection current conductivity, σ abc inj .These two quantities are given by, where, n, m label the energy bands, Numerical calculation of these quantities by direct evaluation of wavefunction derivatives can be difficult as it would require fixing a smooth gauge for the wavefunctions at each point.However, it is possible to circumvent this problem completely by making use of r b nm = −iv b nm /ω nm = −i n| ∂ ∂k b H |m /ω nm , and the sum rule [8,13,59], where, The condition on the summation l =n,m is understood as ω l = ω n , ω m [13].
The consequences of time-reversal symmetry can be seen directly by analyzing the integrand in Eq. (2) and Eq.(3) under a time-reversal operation.Time-reversal symmetry enforces the real part of the integrand to be odd in k space, and hence makes σ abc shift real and σ abc inj imaginary [60].In other words, when time-reversal symmetry is preserved, the shift current conductivity is non-zero only for linearly-polarized light and the injection current requires circularly polarized light.However, no such restrictions are present once the timereversal symmetry is broken.

III. RESULTS FOR THE CHARGE-n LOW-ENERGY WEYL HAMILTONIAN
We begin with a low-energy effective Hamiltonian for a two-band charge-n Weyl point where ζ = ±1, u z and u t are, respectively, the effective velocity and tilt along ẑ.Here, kx,y = k x,y /k 0 , and µ is the chemical potential.The values k 0 , ε 0 are materialdependent parameters with units of momentum and energy, respectively.We will assume ε 0 > 0 and set k 0 = 1.The chirality of this Weyl point is χ = sgn(u z ζ).
The energy eigenvalues are given by, It should be noted that although all our derivations will hold for n being any positive integer, it makes physical sense to only take n = 1, 2, 3 due to symmetry restrictions in actual lattice systems [61,62].Two-band charge-4 Weyl points are allowed but have different low energy Hamiltonian [61,63] and are discussed in a later section.
In order to use Eq.(2), and Eq.(3) to find the shift and injection conductivity tensors, we note that the delta and Fermi-Dirac distribution functions restrict the domain of integration.In our calculations, we assume temperature, T = 0 K which simplifies the Fermi-Dirac distribution to f (E) = 1 − Θ(E) where Θ is the Heaviside function.The delta function forces the integration to be performed over the surface 2ε 0 ( k2 x + k2 y ) n + u 2 z k 2 z /ε 2 0 − ω = 0, while the theta function further selects out a portion of this surface.By making suitable substitutions, this surface can be transformed into a sphere which makes it easier to perform the integral analytically (see Appendix A) for arbitrary charge n.
After accounting for the finite tilt of the Weyl cone, the Pauli blocking condition restricts the integration region on this sphere to region S as shown in Fig. 1 with θ 1 and θ 2 given by: for p = 1, 2 where ϕ p = 1 W sgn ut uz 2µ ω + (−1) p , and W = |u t /u z | is an important quantity which determines if the WSM is type-I (W < 1) or type-II (W > 1).The behavior of θ 1 , θ 2 is mainly determined by the amount of tilt (W ) and doping (µ), and is crucial to understanding the basic features of the response.
For zero doping, θ 2 = −θ 1 = π/2 for type-I and θ 2 = −θ 1 = arcsin(1/W ) for type-II WSM.It should also be noted that the angles lose dependence on chirality in this case.It is important to note that these results contain an implicit ω dependence.In the transformed coordinates, where the integration surface is a sphere, these angles are measured from the x-axis in the xz-plane and determine which part of that surface is not Pauli-blocked (region S in Fig. 1).First, we evaluate the join density of states using the expression where the factor f nm accounts for Pauli-blocking effects.
In the absence of the tilt, we obtain the expected ω 2/n dependence for a charge-n Weyl node.However, at finite tilt and finite chemical potential, this ω 2/n dependence is modulated by the angular factor of θ2 θ1 cos 2/n−1 θ dθ.Next, we calculate different components of shift and injection current tensors.The resulting expressions are given in Table I.We notice that all conductivity tensors are directly proportional to the charge of the Weyl point, except for σ zzz inj .It should be noted that analytical results for the n = 1 and untilted n = 2 cases have been given in Refs.[13,59,64,65] and Ref. [66], respectively.Here, we have extended the analytical results to arbitrary chiral charge-n with finite tilt.
Let's first analyze the shift current conductivity results.As shown in Table I, there are two kinds of non-zero components: (i) purely imaginary which is responsible for a second-order dc photocurrent from circular polarization, and (ii) purely real which leads to a photogalvanic effect from linearly polarized light.For the shift current, the circular polarization components always vanish at zero doping since θ 1 = −θ 2 for µ = 0 from Eq. (7).Similarly, when time-reversal symmetry (TRS) is preserved, the circular polarization current from a time-reversed pair of nodes would also vanish as u z → −u z under time-reversal.
On the other hand, the linear polarization component σ xyz shift shows a very interesting behavior and can even provide estimates of tilt and chemical potential.We note that among all the non-zero conductivity tensors, σ xyz shift alone changes sign with frequency and can be used to estimate µ.For type-I and type-II with W < 2, this sign change occurs at ω = 2|µ| which can be understood from Eq. ( 7) which indicates that while one of the angles is zero the other becomes ±π/2 leading to σ xyz shift = 0.The latter stays at ±π/2 for small variation in ω, while the former changes sign going through ω = 2|µ|, causing σ xyz shift to do the same (as it has a sin θ cos 2 θ dependence).The W ≥ 2 case is not so straightforward but after some work we find that the sign change occurs at 2|µ| Our results show that the tilt parameter plays an important role for all shift current components.When the tilt vanishes, all the shift current conductivity components also vanish.This can be easily understood from the behavior of θ p from Eq.( 7) in the limit W → 0. For ω < 2|µ|, θ 1 = θ 2 = ±π/2 which simply means that the entire ω 21 = ω surface is Pauli-blocked.However, when ω > 2|µ|, the entire surface becomes Pauli-unblocked (as captured by θ 2 = −θ 1 = π/2) which again leads to a vanishing shift conductivity.Now, we turn our attention to injection current conductivity components-some of which are known to exhibit quantization proportional to the Berry charge of the Weyl node.Here again, there are two kind of components: (i) purely imaginary which leads to CPGE, and (ii) purely real which leads to a photogalvanic effect from linearly polarized light.When time-reversal symmetry is preserved, the contribution from timereversed Weyl node pairs is such that the real components vanish and only CPGE survives, as expected.Also, all the real components of the injection current conductivity would disappear at zero doping and also at zero tilt.As a result in order to get an injection current for linearly polarized light not only time-reversal must be broken but doping and tilt should be finite as well.

JDOS
point.For the latter, the response window is proportional to |µ|.
For type-II case, asymptotically while the remaining components asymptotically approach zero.
The quantization condition for CPGE for the injection current conductivity can be easily obtained as the trace of the CPGE tensor, which gives the perfect quantized value equal −n sgn(u z ζ) only when θ 2 = −θ 1 = π/2.The contribution of the factor 1 2 (sin θ 2 − sin θ 1 ) is easy to understand when interpreted as the fraction of the solid angle available for integration, which leads to a reduced value of quantized response when either |θ 1 |, |θ 2 | < π/2.For µ = 0, type-I WSM gives perfect quantization where as in type-II WSM, the quantization value is reduced by a factor of 1/W .When µ = 0, type-I WSMs show perfect quantization above a certain frequency cutoff, i.e for ω ≥ 2|µ|/(1 − W ) where as type-II WSMs show a reduced quantization for ω ≥ 2|µ|/(W − 1), as shown in Fig. 2(i).Note that in the case of type-II, while individual terms in the CPGE trace only approach their respective quantized values asymptotically, the trace itself is fully quantized for ω > 2|µ|/(W − 1).This feature is captured in Fig. 2(g,h,i).
When TRS is broken, injection current can also be generated by linearly polarized light, and the non-zero components for this case depend on the tilt direction.For tilt along the z-axis, non-zero linear photogalvanic effect (LPGE) injection current conductivities include σ yzy , σ yyz , σ xzx , σ xxz , σ zxx , σ zyy , and σ zzz .The last one is the only component among all shift and injection current conductivities which can allow for a current in the direction of linear polarization if it coincides with the direction of the tilt.Thus, a measurement of σ aaa inj can provide a simple way to determine the direction of the tilt in charge 2 and 3 which have linearly dispersing bands in only one direction.

IV. CHARGE-2 WEYL SEMIMETALS
For concreteness, we numerically calculate the conductivity tensors for the following two-band tightbinding model for a charge-2 WSM with broken inversion, time-reversal and mirror symmetries, which has nodes at (k x , k y , k z ) = (0, 0, ± acos(M − 2)) for 1 ≤ M ≤ 3. The low-energy Hamiltonian near the nodes is given by, where, Chirality of this node is given by χ = sgn(u z ).The bands disperse as t(u t ± u z )k z when k x , k y = 0 and ±t(k 2 x + k 2 y ) when k z = 0. Based on the possible range of band inversion strength and k-space node separation given in Ref. [30], we take t = 1eV, g = 0.1, M = 2.958 and µ = −0.0287t.The chosen parameters give W = 0.33 and put the lower energy node near zero energy.
The second-order dc conductivity results obtained for Eq.( 11) are shown in orange in Fig. 3. To compare these results against those obtained simply by treating each node separately based on Eq.( 12), we have included the blue curve which represents the sum total of  I).The orange and green curves correspond to type-I (W1 = 0.334) and type-II (W2 = 2.334) case, respectively.We have chosen W2 − W1 = 2 just to keep the plots neat.We have taken uz = 0.287, ε0 = 1, µ = −0.03.Note that except for the JDOS, remaining plots will show similar behavior for charges 1 and 3.
contributions from individual nodes using expressions from Table I.This is reasonable if the contribution from at least one of the nodes is constant over the energy range under consideration, as is the case here.Looking at Fig. 3, it is clear that higher order terms present in the tight binding model lead to significant deviations from the low-energy predictions of Table I. Surprisingly, we find that the injection conductivities σ xxz , σ zxx and σ zzz (d-f), develop a plateau up to an energy of about 0.07t.Additionally, a shift in the response energy window to the left by about 0.02t is seen for all the conductivities (b-h) and the CPGE quantization (i).A shift in the quantization window has also been seen for T -symmetric charge-2 Weyl system [58] and is believed to arise from higher order terms.
We probe the origin of these deviations by explicitly adding higher order terms to Eq. ( 12) (see Appendix B for details).Specifically, we find that including the second order terms ( 1 2 (k ), not only matches the energy shift, but also captures all other features of the tight-binding results.Most importantly, we find the plateaus to come from the node situated close to zero energy (in our case, it is the node with negative u z ) and their heights to be −sgn(uz) 64 , sgn(uz) 64 and sgn(uz)(2−M )
We should point out that the σ xxz , σ zxx plateaus can be obtained by just including the 1 2 (k 2 x +k 2 y )σ z term where as the one for σ zzz can be explained with the ( 1 2 M − 1)k 2 z )σ z term alone.We believe that the plateaus should be present in any charge-2 WSM where these higher order terms show up.Note that despite the energy shift, σ xyz Lastly, we also find perfect CPGE quantization up to  11) obtained numerically.The blue curve is sum of contributions from the two nodes based on analytical results for the low energy model Eq.( 12) while the green curve is obtained by including higher order terms 12) (see Appendix B for details).The energy separation between the nodes is |2guz| and µ = µ − guz.We have taken t = 1, M = 2.958, g = 0.1, µ = −0.0287.(a) JDOS for each node for the low energy model with (solid) and without (dashed) higher order correction terms.Pink and purple correspond to the nodes at −µ − guz and −µ + guz, respectively.
an energy of about 0.07t as shown in Fig 3(i).Behavior of the JDOS for each node with (solid) and without (dashed) higher order corrections is shown in Fig. 3(a).With the higher order terms, the JDOS for the node at energy −µ + gu z becomes non-zero after about 0.071t and explains why the quantization ceases earlier than the predicted value 0.086t.
Beyond two-band models, the CPGE quantization is no longer guaranteed to hold.This has been explored in Ref. [19] for a charge-1 WSM.In order to better understand the contributions coming from higher bands and the extent to which they destroy quantization in charge-2 WSM, we study the following four-band tightbinding model (taking inspiration from Ref. [29]) with broken time-reversal symmetry, where τ and σ are Pauli matrices acting on the orbital and spin space, respectively.
where u z = ± (3 − M )(M − 1), and u t = g(M − 2).We begin with |∆| large compared to |gu z | and gradually decrease it to below |gu z |.The two bands which touch, disperse as t(u t + u z )k z , t(u t − u z )k z when k x , k y = 0 and x +k 2 y ) 2|∆−guz| when k z = 0. We use the same g, M, t values from before.We note that unlike the two-band case Eq.( 12), the quadratic dispersion now has a dependence on ∆.For ∆ = 0.5 (recall gu z = 0.0287), the dispersion becomes almost the same for the two cases and provides a good starting point for comparison.Also, since the gap between the highest occupied and the lowest unoccupied bands is ∼ |∆|, the effect of higher bands should be more prominent for smaller values of ∆.
Results obtained for Eq.( 13) are shown in Fig. 4. We find large deviation from perfect quantization for small gaps as seen in Fig. 4(f).However, for ∆ |gu z | we do see almost perfect quantization.Also, the plateaus seen earlier in σ xxz inj , σ zxx inj , σ zzz inj continue to show up when ∆ is at least a few times larger than |gu z | as shown in (a-c).Their heights become dependent on ∆ and are empirically found to be about ∆/32, −∆/32, and ∆(M − 2)/32, respectively.

V. CHARGE-4 WEYL SEMIMETALS
Having looked at the charge-2 case in some detail, we move on to investigate the behaviour of the injection conductivity and JDOS for charge-4 WSMs.
The existence of CPGE quantization in such systems has been discussed in earlier studies [63].In our study, we want to develop a full understanding of how model parameters and doping affect these responses.In order to do that, we take the following two-band Hamiltonian based on Ref. [63], which has nodes of opposite chirality at (0, 0, 0) and (π, π, π).The low-energy Hamiltonian near Γ-point is given by, where µ = µ + 6c 1 .The chirality of this Weyl point is given by χ = sgn(c 3 ).We derive all our results using Eq. ( 16) with c 1 > 0 (the opposite case is an easy generalization, which we discuss later).Its eigenvalues  16) with c1, ω > 0 (when c1 < 0, replace c1 → |c1|, µ → −µ and use the c1 > 0 results).Integrals given here are to be evaluated over the region satisfying Here, we have defined 2π iτ e 3 / 2 abc σ abc are given by The presence of the sixth order term k 2 x k 2 y k 2 z above does not allow us to fully evaluate Eq.(2), Eq.(3), and Eq.( 8) analytically.However, it is possible to integrate out k y (one can pick any one out of the three k coordinates) and get rid of the delta function in exchange for a new constraint (see Appendix E for details).
The biggest advantage of going from a triple to a double integral is that the new constraint now defines a closed area compared to a closed surface before, which makes it much easier to analyze.The expressions for JDOS and injection conductivities (nonzero components) thus obtained are given in Table II.Note that the shift conductivities are zero.Although these integrals appear complicated, they are easy to evaluate numerically.A key result of our analysis is the precise location of the energy window and condition under which trace of the CPGE tensor is quantized for different amounts of doping.
When µ = 0, quantization is seen only for c 1 /|c 2 | < 1 starting at a frequency of as shown in Fig. 5 (a).The situation for 1 < c 1 /|c 2 | < 2 and 2 < c 1 /|c 2 | is also shown in Fig. 5 (b) and (c), respectively.While the trace is non-zero for any finite frequency in the former case, it turns non-zero only after where ω p is the unique real positive root of (ω − 2|µ|) 3 − 54 this happens after ω p .The three cases are shown in Fig. 5 (d), (e) and (f), respectively.
When µ > 0, we are presented with a wider range of possibilities for observing quantization.We find that, irrespective of the c 1 /|c 2 | value, the trace becomes nonzero after min ω p , 2µ When this condition is met, perfect quantization is seen but only for a finite window of energy ω satisfying max ω p , 2µ which is shown in Fig. 5(h).The situation when no quantization is possible for µ > 0 is shown in Fig. 5(i).It is clear that while |c 1 /c 2 | plays a crucial role, c 1 , c 2 , c 3 and µ intricately determine the behavior of the CPGE trace and its quantization.Note that the plots in Fig. 5 have been obtained by numerically evaluating the integrals in Table II.We have included the JDOS plot at zero doping in Fig. 6(a).As shown in the figure, the JDOS has a √ ω behavior going towards zero frequency.Note that the JDOS result from Table I also predicts a √ ω dependence if we take n = 4.This seems more like a coincidence as that model still has a linearly dispersing band along k z , very different from the C-4 model in Eq.( 16).
We would like to point out that so far the results presented assume c 1 > 0. It turns out that we can continue to use the same results for a charge-4 node with 0.00 0.05 0.10 0.15 0.20 0.25 0.00 0.05 0.10 0.15 0.20 0.25  16) with c1 > 0, obtained from the numerical evaluation of integrals in Table II for different combinations of µ and c1/|c2| (respective values shown in the inset).We have taken c1 = 0.0665, c3 = 0.4.(a)-(c) µ = 0, (d)-(f) µ < 0. In both these cases, perfect quantization is seen only when c1/|c2| < 1. (g)-(i) µ > 0, perfect quantization is guaranteed for c1/|c2| < 1 however, unlike previous two cases, it can also be seen for c1/|c2| > 1 as long as c 1 < 0 (and a chemical potential µ) by treating it as a |c 1 | node with chemical potential −µ.With this small but important extension, our analysis covers all the possible cases.
In Fig. 6(b), we chose c 1 /|c 2 | < 1 and therefore expect both nodes to show perfect quantization.Since µ = 0, the two nodes will show quantization starting at different frequencies which results in an overall finite quantization window.The dashed curves represent the sum of contributions from the two nodes based on the low energy result from Table II (as remarked earlier, this makes sense here because at no point in the energy range under consideration do the contributions from both nodes become non-constant simultaneously).
In Fig. 6(c), we chose 1 < c 1 /|c 2 | and µ = −0.3.This gives µ Γ = 0.1 and µ R = −0.7.Since the node closer to zero energy falls under the µ > 0 category when using the low-energy results, we can choose c 1 , c 2 , c 3 such that it shows quantization for a finite window (blue curve) or no quantization at all (red curve).For the former, we also ensure that the contribution from other node starts only after the end of the quantization window.The dashed curves show contribution from the Γ-node in both cases.In both Fig. 6(b) and (c), we find excellent  15) with µ = −0.3,0.5, respectively.Note that µΓ = µ + 6c1, µR = µ − 6c1.The corresponding dashed green curve is obtained by evaluating expressions from Table II for each node separately and then adding the results.We have used the same c1, c2, c3 from before.(c) Results with c1/c2 > 1, µ = −0.3(same c1, c3 as before).The dashed curves show contribution from the node closer to E = 0 (obtained using results from Table II) for each case and are discontinued after contribution from the other node becomes non-zero.
agreement between the tight binding and low energy results, showing that the higher order terms are not at play in this parameter range and can be neglected.

VI. CONCLUSION AND DISCUSSION
In summary, we have presented a comprehensive and unified study of the second-order dc response in tilted multi-Weyl systems with a focus on the roles played by tilt (W ) and doping (µ).For charges n=1, 2 and 3, we have derived analytical expressions for shift and injection conductivity using a low energy continuum model and then compared its predictions against more realistic two-and four-band tight binding models of timereversal broken the charge-2 case.
extremely CPGE quantization, we also report other features of the photogalvanic response arising mainly from the finite tilt and band curvatures.We systematically investigated the role of tilt, band curvatures, and higher bands in deciding the shift and injection current conductivities of multi-Weyl semimetals.We find that in TRS broken multi-Weyl semimetals, finite tilt can lead to non-zero injection current from linearly polarized light which not only provides a probe for the tilt direction but can possibly also provide a way to engineer the injection current by using strain or some other mechanism which controls the tilt of Weyl nodes.
We have also provided the first complete analysis of the photogalvanic response in charge-4 WSM based on a lowenergy two-band model, covering all possibilities arising from different combinations of model parameters and the chemical potential.Although C-4 WSMs do not have a tilt in the usual sense (like the other three charges which have linear dispersing band in at least one direction), the ratio |c 1 /c 2 | plays a similar role, and together with c 3 1 /c 2 3 and µ determines nature of the response.Within the confines of the low-energy model, our results help point out exactly when CPGE quantization can be seen in C-4 WSMs.
We believe that the new approach we have taken here to study C-4 would find applications in studying many other optical responses as well.For example, it can easily be extended to the study of SHG and first-order conductivity for the low energy two band model.In principle it should work with any quantity that requires evaluating a k-space integral with f 21 δ(ω 21 − ω) term in it at T = 0 K.We work with the low-energy effective Hamiltonian, with eigenvalues The domain for the integrals in Eq.(2), Eq.( 3) is determined by f 21 and δ(ω Let us focus on the delta function first.To simplify things, we split the integral in k x − k y plane over the four quadrants: dk x dk y = + dk x + dk y + + dk x − dk y + − dk x + dk y + − dk x − dk y , and combine them into a single integral over the first quadrant by making substitutions k x = ±k 0 √ x and k y = ±k 0 √ y depending on the sign (both x, y > 0).We also put k z = ε0 uz z.By making x → x−y √ 2 and y → x+y √ 2 , we rotate the x − y axis counterclockwise by π/4, and scale z by z → 2 n/4 z.Finally, we let x → x 2/n to get δ(2 1+n/4 ε 0 √ x 2 + z 2 − ω).These transformations also change the integration measure k → , with x > 0 and −x < y < x (we do the y integral first with these limits).The delta function defines a circle in xz-plane which lets us use x = r cos θ, z = r sin θ to obtain δ(r − ω/(2 1+n/4 ε 0 ))/(2 1+n/4 ε 0 ).
Since we are taking the temperature to be zero, f 21 = Θ(E 1 ) − Θ(E 2 ) where, Θ is the Heaviside step function.Because of the condition put by the delta function, we have Since we have assumed ω > 0, the only non-zero value for f 21 is −1 when µ − ω/2 < u z k z < µ + ω/2.Using coordinate transformations from before, this condition becomes where W = |u t /u z |.The definitions for θ 1 , θ 2 given in main text follow from this.With this, we can easily compute other ingredients of the integral from the eigenvalues and normalized eigenfunctions of H n , and combine them to obtain analytical expressions for the JDOS, shift and injection conductivities.
We can also include additional terms of the form A( k2 x + k2 y )σ 0 and Bσ z into Eq.( 5).These correspond to second-order tilt and Zeeman terms, respectively.The B term only shifts the origin along k z , modifying the k z → z transformation to k z = ε0z−B uz .It not difficult to see that these terms only affect f 21 , where µ = µ + Bu t /u z .Since ω > 0, we have with −π/2 ≤ θ, φ ≤ π/2.Also, we define α = sgn ut uz φ and ϕ p = 1 W sgn ut uz 2 µ ω + (−1) p with p = 1, 2. The inequality becomes ϕ 1 < sin(θ + α) < ϕ 2 which can be solved for the minimum ( θ 1 ) and maximum ( θ 2 ) allowed values of θ.These can be obtained as the left and right end points of the intervals, where, For µ = 0, the theta function constraints reduce to x + z < 1.When c 1 /|c 2 | < 1, the CPGE trace is nonzero for any finite ω and becomes ±4 after .When c 1 /|c 2 | > 1, some portion of Eq.(E2) is always left out and we do not get perfect quantization.For 1 < c1 |c2| < 2, the trace is non-zero for any finite ω where as for c1 |c2| > 2, this happens only after ω .An important thing to note here is that the term 1 − 2|µ| ω ∈ (−∞, 1).We are only interested when it lies in (0, 1) which happens for ω > 2|µ|.Since it can ever only reach 1, full overlap with region Eq.(E2) is possible if c 1 /|c 2 | < 1, the condition to get perfect quantization.When this condition is met, the amount of overlap between x + z < 1 − 2|µ| For µ > 0, the possibilities become even more interesting.Θ −x − z + 1 + 2µ ω Θ x + z + 1 − 2µ ω sets bounds on the integration region, requiring 2µ ω − 1 < x + z < 2µ ω + 1.The term 1 + 2µ ω ∈ (1, ∞), which means that if c 1 /|c 2 | > 1, a portion of Eq.(E2) will necessarily be left out for ω > 2µ ω 2 = 0.The product of its roots is 8µ 3 > 0 which means that when two roots are complex (conjugate pair), the real root must be positive.However, when all roots are real, there are two possibilities − one positive two negative roots or three positive roots.It turns out the latter case is not possible because the condition for all roots being real is µ < 4 where as for all roots to be positive, µ > 9
shift crosses zero at 2|µ−gu z | as seen from Fig 3(b).This may not hold for arbitrarily chosen µ.In our case, we have carefully put one node close to zero energy which makes the other node almost entirely dictate the behavior of the response.

1+ c 1
|c 2 | where ω p is now the unique real positive root of the cubic equation (ω − 2µ) 0. For c 1 /|c 2 | < 1, it goes on to reach a saturation value of ±4 at max ω p , in Fig. 5(g).Interestingly for µ > 0, quantization becomes possible even for 1

TABLE I .
Results for the low-energy charge-n Hamiltonian in Eq.(

TABLE II .
Results for JDOS and injection conductivity for the Hamiltonian in Eq.(