# Photogalvanic response in multi-Weyl semimetals

Arpit Raj<EMAIL_ADDRESS>Department of Physics, Northeastern
University, Boston, MA 02115, USA  Swati Chaudhary  Department of Physics, The
University of Texas at Austin, Austin, Texas 78712, USA  Department of
Physics, Northeastern University, Boston, MA 02115, USA  Department of
Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Gregory A. Fiete  Department of Physics, Northeastern University, Boston, MA
02115, USA  Department of Physics, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA

###### Abstract

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 SrSi2. 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 non-centrosymmetric 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 real-space shift in the electron wavepacket due to inter-band
photoexcitation [12], and the injection current is caused 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 [7, 6]. 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].

In recent years, many different kinds of Weyl semimetals have been discovered
[24, 25, 26, 27]. In certain materials, it has been shown that the Weyl node
carrying a charge higher than $n=1$ can be stabilized by crystal symmetries
[28]. These semimetals, also known as multi-Weyl semimetals (MWSMs) have been
proposed in SrSi2 [29, 30], Cu2Se, and RhAs3 [31] which can host Weyl nodes
with charge $n=2$. It was shown in Ref. [32, 33] that such double-Weyl nodes
can also be engineered in Luttinger semimetals like $\alpha$-Sn by applying
strain and magnetic fields or via Floquet engineering.

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 SrSi2 [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 SrSi2, 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 four-band 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 tight-binding 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 tight-binding model and an effective low-energy
Hamiltonian. Next, we consider a four-band model inspired by the SrSi2 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 $\omega$
(and wavelength much larger than the sample size so the electric field has
uniform amplitude), the second-order dc response is given by,

$\displaystyle j_{dc}^{a}$
$\displaystyle=\sigma^{abc}(\omega)E_{b}(\omega)E_{c}(-\omega),$ (1)

where the second-order conductivity $\sigma^{abc}(\omega)$ can be divided into
a shift current conductivity, $\sigma_{\text{shift}}^{abc}$ and an injection
current conductivity, $\sigma_{\text{inj}}^{abc}$. These two quantities are
given by,

$\displaystyle\begin{split}\sigma_{\text{shift}}^{abc}&=\frac{-i\pi
e^{3}}{\hbar^{2}}\int_{\mathbf{k}}\sum_{n>m}f_{nm}\Big{(}r_{nm}^{b}r_{mn;a}^{c}-r_{mn}^{c}r_{nm;a}^{b}\Big{)}\\\
&\hskip 128.0374pt\times\delta(\omega_{nm}-\omega),\end{split}$ (2)
$\displaystyle\sigma_{\text{inj}}^{abc}$ $\displaystyle=\tau\frac{2\pi
e^{3}}{\hbar^{2}}\int_{\mathbf{k}}\sum_{n>m}f_{nm}\Delta_{nm}^{a}r_{nm}^{b}r_{mn}^{c}\delta(\omega_{nm}-\omega),$
(3)

where, $n,m$ label the energy bands,
$\int_{\mathbf{k}}=\int\differential[3]{k}/(2\pi)^{3}$,
$\omega_{nm}=\omega_{n}-\omega_{m}$ is the energy difference between bands $n$
and $m$, $f_{nm}=f_{n}-f_{m}$ where $f$ is the Fermi-Dirac distribution
function, $\Delta_{nm}^{a}=v^{a}_{nn}-v^{a}_{mm}$ with $v^{a}_{nn}$ being the
velocity matrix elements, and $\tau$ is the relaxation time. The interband
Berry connection is given by
$r_{nm}^{b}=\bra{n}i\partialderivative{k_{b}}\ket{m}\text{ for }n\neq m$ and
zero otherwise, with its generalized derivative defined as
$r_{nm;a}^{b}=\partialderivative{r_{nm}^{b}}{k_{a}}-i(\xi_{nn}^{a}-\xi_{mm}^{a})r_{nm}^{b}$,
where $\xi_{nn}^{a}=\bra{n}i\partialderivative{k_{a}}\ket{n}$ is the intraband
Berry connection.

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_{nm}^{b}=-iv_{nm}^{b}/\omega_{nm}=-i\bra{n}\partialderivative{k_{b}}\mathcal{H}\ket{m}/\omega_{nm}$,
and the sum rule [8, 59, 13],

$\displaystyle\begin{split}r_{nm;a}^{b}&=\frac{i}{\omega_{nm}}\Bigg{[}\frac{\Delta_{nm}^{b}v_{nm}^{a}+\Delta_{nm}^{a}v_{nm}^{b}}{\omega_{nm}}-w_{nm}^{ba}\\\
&\qquad+\sum_{l\neq
n,m}\bigg{(}\frac{v_{nl}^{b}v_{lm}^{a}}{\omega_{lm}}-\frac{v_{nl}^{a}v_{lm}^{b}}{\omega_{nl}}\bigg{)}\Bigg{]},\quad
n\neq m\end{split}$ (4)

where,
$w_{nm}^{ba}=\bra{n}\partialderivative{k_{b}}\partialderivative{k_{a}}\mathcal{H}\ket{m}$.
The condition on the summation $\sum_{l\neq n,m}$ is understood as
$\omega_{l}\neq\omega_{n},\omega_{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
$\vec{k}$ space, and hence makes $\sigma^{abc}_{\text{shift}}$ real and
$\sigma^{abc}_{\text{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 time-
reversal 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 $\zeta=\pm 1$, $u_{z}$ and $u_{t}$ are, respectively, the effective
velocity and tilt along $\hat{\mathbf{z}}$. Here,
$\tilde{k}_{x,y}=k_{x,y}/k_{0}$, and $\mu$ is the chemical potential. The
values $k_{0},\varepsilon_{0}$ are material-dependent parameters with units of
momentum and energy, respectively. We will assume $\varepsilon_{0}>0$ and set
$k_{0}=1$. The chirality of this Weyl point is $\chi=\text{sgn}(u_{z}\zeta)$.
The energy eigenvalues are given by,

$\displaystyle E_{n,\pm}$
$\displaystyle=u_{t}k_{z}-\mu\pm\varepsilon_{0}\sqrt{(\tilde{k}_{x}^{2}+\tilde{k}_{y}^{2})^{n}+u_{z}^{2}k_{z}^{2}/\varepsilon_{0}^{2}}.$
(5)

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-\Theta(E)$ where $\Theta$ is the Heaviside function. The delta
function forces the integration to be performed over the surface
$2\varepsilon_{0}\sqrt{(\tilde{k}_{x}^{2}+\tilde{k}_{y}^{2})^{n}+u_{z}^{2}k_{z}^{2}/\varepsilon_{0}^{2}}-\omega=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 $\theta_{1}$ and $\theta_{2}$ given by:

$\displaystyle\theta_{p}$ $\displaystyle=\begin{dcases}-\pi/2,\text{ if
}\varphi_{p}<-1\\\ \arcsin(\varphi_{p}),\text{ if }-1\leq\varphi_{p}\leq 1\\\
+\pi/2,\text{ if }1<\varphi_{p}\end{dcases},$ (6)

for $p=1,2$ where
$\varphi_{p}=\frac{1}{W}\left(\text{sgn}\left(\frac{u_{t}}{u_{z}}\right)\frac{2\mu}{\omega}+(-1)^{p}\right)$,
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 $\theta_{1},\theta_{2}$ is
mainly determined by the amount of tilt $(W)$ and doping $(\mu)$, and is
crucial to understanding the basic features of the response.

For zero doping, $\theta_{2}=-\theta_{1}=\pi/2$ for type-I and
$\theta_{2}=-\theta_{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 $\omega$ 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).

Figure 1: The surface defined by $\delta(\omega_{21}-\omega)$ in the
transformed coordinates (see Appendix A). The factor $f_{21}$ restricts the
integral in Eq.(2), Eq.(3), Eq.(7) to the Pauli-unblocked region S (shown in
brown).

First, we evaluate the join density of states using the expression

$\displaystyle\text{JDOS}(\omega)$
$\displaystyle=\int_{\mathbf{k}}\sum_{n>m}f_{mn}\delta(\omega_{nm}-\omega),$
(7)

where the factor $f_{nm}$ accounts for Pauli-blocking effects. In the absence
of the tilt, we obtain the expected $\omega^{2/n}$ dependence for a charge-$n$
Weyl node. However, at finite tilt and finite chemical potential, this
$\omega^{2/n}$ dependence is modulated by the angular factor of
$\int_{\theta_{1}}^{\theta_{2}}\cos^{2/n-1}{\theta}\differential{\theta}$.

Table 1: Results for the low-energy charge-$n$ Hamiltonian in Eq.(LABEL:eq:Hn) with $\varepsilon_{0},\omega>0$. Note that $\chi=\text{sgn}(u_{z}\zeta)$. $\mathrm{JDOS}$ |  $\frac{k_{0}^{2}}{8\pi^{2}|u_{z}|}\frac{1}{n}\left(\frac{\omega}{2\varepsilon_{0}}\right)^{2/n}\int_{\theta_{1}}^{\theta_{2}}\cos^{2/n-1}{\theta}\differential{\theta}$
---|---
Shift conductivity |
$\sigma^{xzx}=-\sigma^{xxz}=\sigma^{yzy}=-\sigma^{yyz}$ |  $n\frac{i\mathrm{sgn}(u_{z})e^{3}k_{0}^{2}}{32\pi\hbar^{2}}\big{(}\sin^{2}{\theta_{2}}-\sin^{2}{\theta_{1}}\big{)}\frac{1}{\omega}$
$\sigma^{xyz}=\sigma^{xzy}=-\sigma^{yzx}=-\sigma^{yxz}$ |  $n\frac{\mathrm{sgn}(u_{z}\zeta)e^{3}k_{0}^{2}}{32\pi\hbar^{2}}\big{(}\sin{\theta_{2}}\cos^{2}{\theta_{2}}-\sin{\theta_{1}}\cos^{2}{\theta_{1}}\big{)}\frac{1}{\omega}$
Injection conductivity |
$\sigma^{zxy}=-\sigma^{zyx}$ |  $n\frac{i\tau\mathrm{sgn}(u_{z}\zeta)e^{3}k_{0}^{2}}{24\pi\hbar^{2}}\big{(}\sin^{3}{\theta_{1}}-\sin^{3}{\theta_{2}}\big{)}$
$\sigma^{yzy}=\sigma^{yyz}=\sigma^{xzx}=\sigma^{xxz}$ |  $n\frac{\tau\mathrm{sgn}(u_{z})e^{3}k_{0}^{2}}{64\pi\hbar^{2}}\big{(}\cos^{4}{\theta_{1}}-\cos^{4}{\theta_{2}}\big{)}$
$\sigma^{xyz}=-\sigma^{xzy}=\sigma^{yzx}=-\sigma^{yxz}$ |  $n\frac{i\tau\mathrm{sgn}(u_{z}\zeta)e^{3}k_{0}^{2}}{48\pi\hbar^{2}}\big{(}3\sin{\theta_{1}}-\sin^{3}{\theta_{1}}-3\sin{\theta_{2}}+\sin^{3}{\theta_{2}}\big{)}$
$\sigma^{zxx}=\sigma^{zyy}$ |  $n\frac{\tau\mathrm{sgn}(u_{z})e^{3}k_{0}^{2}}{256\pi\hbar^{2}}\big{(}-6\cos(2\theta_{1})+\cos^{2}(2\theta_{1})+6\cos(2\theta_{2})-\cos^{2}(2\theta_{2})\big{)}$
$\sigma^{zzz}$ |  $\frac{\tau u_{z}^{2}\mathrm{sgn}(u_{z})e^{3}k_{0}^{2}}{2^{2+2/n}\varepsilon_{0}^{2/n}(1+n)\pi\hbar^{2}\omega^{2-2/n}}\big{(}\cos^{2+2/n}{\theta_{2}}-\cos^{2+2/n}{\theta_{1}}\big{)}$

Next, we calculate different components of shift and injection current
tensors. The resulting expressions are given in Table 1. We notice that all
conductivity tensors are directly proportional to the charge of the Weyl
point, except for $\sigma^{zzz}_{\text{inj}}$. It should be noted that
analytical results for the $n=1$ and untilted $n=2$ cases have been given in
Refs. [59, 64, 65, 13] 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
1, 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 $\theta_{1}=-\theta_{2}$ for $\mu=0$ from Eq. (6).
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}\rightarrow-u_{z}$ under time-reversal.

On the other hand, the linear polarization component
$\sigma^{xyz}_{\text{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, $\sigma^{xyz}_{\text{shift}}$ alone changes
sign with frequency and can be used to estimate $\mu$. For type-I and type-II
with $W<2$, this sign change occurs at $\omega=2|\mu|$ which can be understood
from Eq. (6) which indicates that while one of the angles is zero the other
becomes $\pm\pi/2$ leading to $\sigma^{xyz}_{\text{shift}}=0$. The latter
stays at $\pm\pi/2$ for small variation in $\omega$, while the former changes
sign going through $\omega=2|\mu|$, causing $\sigma^{xyz}_{\text{shift}}$ to
do the same (as it has a $\sin{\theta}\cos^{2}{\theta}$ dependence).

The $W\geq 2$ case is not so straightforward but after some work we find that
the sign change occurs at $2|\mu|\sqrt{\tfrac{3}{W^{2}-1}}$ (see Appendix D
for details). Note that
$\tfrac{2|\mu|}{1+W}<2|\mu|\sqrt{\tfrac{3}{W^{2}-1}}\leq 2|\mu|$ with the
equality holding at $W=2$, as one would expect. Interestingly, for $\mu=0$,
both components $\sigma^{xyz}_{\text{shift}},\sigma^{xxz}_{\text{shift}}$ show
a $1/\omega$ divergence for a type-II WSM. Additionally, for type-I, all shift
current conductivities are non-zero (shown in Fig. 2(b,c)) only in a finite
frequency window determined by the tilt parameter $W$ and doping.

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
$\theta_{p}$ from Eq.(6) in the limit $W\rightarrow 0$. For $\omega<2|\mu|$,
$\theta_{1}=\theta_{2}=\pm\pi/2$ which simply means that the entire
$\omega_{21}=\omega$ surface is Pauli-blocked. However, when $\omega>2|\mu|$,
the entire surface becomes Pauli-unblocked (as captured by
$\theta_{2}=-\theta_{1}=\pi/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 time-reversed 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.

For finite doping, the conductivities become non-zero after $2|\mu|/(1+W)$ for
both type-I and type-II WSMs (note that the $1/\omega$ divergence gets cut off
in case of type-II). For type-I,
$\sigma^{xyz}_{\text{inj}},\sigma^{yzx}_{\text{inj}},\sigma^{zxy}_{\text{inj}}$
reach their quantized value of $-n\,\text{sgn}(u_{z}\zeta)/12\pi$ after
$2|\mu|/(1-W)$ where as other components become zero beyond this point. For
the latter, the response window is proportional to $|\mu|$. For type-II case,
$\sigma^{xyz}_{\text{inj}},\sigma^{yzx}_{\text{inj}},\sigma^{zxy}_{\text{inj}}$
approach their respective quantized values
$-n\,\text{sgn}(u_{z}\zeta)\frac{3W^{2}-1}{12\pi W^{3}}$,
$-n\,\text{sgn}(u_{z}\zeta)\frac{3W^{2}-1}{12\pi W^{3}}$,
$\frac{-n\,\text{sgn}(u_{z}\zeta)}{6\pi W^{3}}$ 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,

$\displaystyle\frac{2\pi}{i\tau
e^{3}k_{0}^{2}/\hbar^{2}}\epsilon_{abc}\sigma^{abc}_{\text{inj}}$
$\displaystyle=-n\,\mathrm{sgn}(u_{z}\zeta)\left(\frac{\sin{\theta_{2}}-\sin{\theta_{1}}}{2}\right),$
(18)

which gives the perfect quantized value equal $-n\,\mathrm{sgn}(u_{z}\zeta)$
only when $\theta_{2}=-\theta_{1}=\pi/2$. The contribution of the factor
$\tfrac{1}{2}(\sin{\theta_{2}}-\sin{\theta_{1}})$ is easy to understand when
interpreted as the fraction of the solid angle available for integration,

$\displaystyle\begin{split}\frac{1}{4\pi}\int_{S}\differential{\Omega}&=\frac{1}{4\pi}\Big{(}4\pi-\big{(}2\pi(1-\sin{\theta_{2}})\\\
&\hskip 71.13188pt+2\pi(1+\sin{\theta_{1}})\big{)}\Big{)},\\\ \end{split}$
(19)

which leads to a reduced value of quantized response when either
$|\theta_{1}|,|\theta_{2}|<\pi/2$. For $\mu=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 $\mu\neq 0$, type-I WSMs show perfect quantization above
a certain frequency cutoff, i.e for $\omega\geq 2|\mu|/(1-W)$ where as type-II
WSMs show a reduced quantization for $\omega\geq 2|\mu|/(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 $\omega>2|\mu|/(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
$\sigma^{yzy},\sigma^{yyz},\sigma^{xzx},\sigma^{xxz},\sigma^{zxx},\sigma^{zyy},$
and $\sigma^{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 $\sigma^{aaa}_{\text{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.

Figure 2: Plots showing (a) JDOS, (b)-(c) shift conductivity, (d)-(h)
injection conductivity, (i) CPGE quantization for a single charge-2 Weyl point
(obtained using the expressions in Table 1). The orange and green curves
correspond to type-I ($W_{1}=0.334$) and type-II ($W_{2}=2.334$) case,
respectively. We have chosen $W_{2}-W_{1}=2$ just to keep the plots neat. We
have taken $u_{z}=0.287,\varepsilon_{0}=1,\mu=-0.03$. Note that except for the
JDOS, remaining plots will show similar behavior for charges 1 and 3.

##  IV Charge-2 Weyl semimetals

For concreteness, we numerically calculate the conductivity tensors for the
following two-band tight-binding model for a charge-2 WSM with broken
inversion, time-reversal and mirror symmetries,

$\displaystyle\mathcal{H}^{2b}$ $\displaystyle=$ $\displaystyle
t\big{(}2(\cos(k_{y})-\cos(k_{x}))\sigma_{x}+2\sin(k_{x})\sin(k_{y})\sigma_{y}$
(20) $\displaystyle\quad+(M-\cos(k_{x})-\cos(k_{y})-\cos(k_{z}))\sigma_{z}$
$\displaystyle\quad+g\sin{k_{z}}\sigma_{0}\big{)}-\mu\sigma_{0},$

which has nodes at $(k_{x},k_{y},k_{z})=(0,0,\pm\acos(M-2))$ for $1\leq M\leq
3$. The low-energy Hamiltonian near the nodes is given by,

$\displaystyle\begin{split}\mathcal{H}^{2b}_{\pm}&=t\big{(}(k_{x}^{2}-k_{y}^{2})\sigma_{x}+2k_{x}k_{y}\sigma_{y}+u_{z}k_{z}\sigma_{z}\\\
&\quad+\left(u_{t}k_{z}-(\mu/t-gu_{z})\right)\sigma_{0}\big{)},\end{split}$
(21)

where, $u_{z}=\pm\sqrt{(3-M)(M-1)}$, and $u_{t}=g(M-2)$. Chirality of this
node is given by $\chi=\text{sgn}(u_{z})$. The bands disperse as $t(u_{t}\pm
u_{z})k_{z}$ when $k_{x},k_{y}=0$ and $\pm t(k_{x}^{2}+k_{y}^{2})$ 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=1\text{eV},\,g=0.1,\,M=2.958$ and $\mu=-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.(20) are shown in
orange in Fig. 3. To compare these results against those obtained simply by
treating each node separately based on Eq.(21), we have included the blue
curve which represents the sum total of contributions from individual nodes
using expressions from Table 1. 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.

Figure 3: Injection and shift conductivities for type-I charge-2 WSM. The
orange curve represents the response for the two-band tight binding model
Eq.(20) obtained numerically. The blue curve is sum of contributions from the
two nodes based on analytical results for the low energy model Eq.(21) while
the green curve is obtained by including higher order terms
$(\tfrac{1}{2}(k_{x}^{2}+k_{y}^{2})+(\tfrac{1}{2}M-1)k_{z}^{2})\sigma_{z}-\tfrac{1}{2}gu_{z}k_{z}^{2}\sigma_{0}$
in Eq.(21) (see Appendix B for details). The energy separation between the
nodes is $|2gu_{z}|$ and $\widetilde{\mu}=\mu-gu_{z}$. We have taken
$t=1,M=2.958,g=0.1,\mu=-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 $-\mu-gu_{z}$ and $-\mu+gu_{z}$,
respectively.

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 1. Surprisingly, we find that the injection conductivities
$\sigma^{xxz},\sigma^{zxx}$ and $\sigma^{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 $\mathcal{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. (21) (see Appendix B for details). Specifically, we find that
including the second order terms
$(\tfrac{1}{2}(k_{x}^{2}+k_{y}^{2})+(\tfrac{1}{2}M-1)k_{z}^{2})\sigma_{z}-\tfrac{1}{2}gu_{z}k_{z}^{2}\sigma_{0}$
in Eq.(21), 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
$\frac{-\text{sgn}(u_{z})}{64}$, $\frac{\text{sgn}(u_{z})}{64}$ and
$\frac{\text{sgn}(u_{z})(2-M)}{64}$, respectively.

We should point out that the $\sigma^{xxz},\sigma^{zxx}$ plateaus can be
obtained by just including the $\tfrac{1}{2}(k_{x}^{2}+k_{y}^{2})\sigma_{z}$
term where as the one for $\sigma^{zzz}$ can be explained with the
$(\tfrac{1}{2}M-1)k_{z}^{2})\sigma_{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, $\sigma^{xyz}_{\text{shift}}$
crosses zero at $2|\mu-gu_{z}|$ as seen from Fig 3(b). This may not hold for
arbitrarily chosen $\mu$. 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.

Lastly, we also find perfect CPGE quantization up to 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 $-\mu+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
tight-binding model (taking inspiration from Ref. [29]) with broken time-
reversal symmetry,

$\displaystyle\begin{split}\mathcal{H}^{4b}&=t\big{(}\sin(k_{x})\tau_{x}+\sin(k_{y})\tau_{y}\\\
&\quad+(M-\cos(k_{x})-\cos(k_{y})-\cos(k_{z}))\tau_{z}\\\
&\quad+\Delta\left(\tau_{x}\sigma_{x}+\tau_{y}\sigma_{y}\right)+g\sin(k_{z})\tau_{z}\sigma_{z}\big{)}-\mu,\end{split}$
(22)

where $\tau$ and $\sigma$ are Pauli matrices acting on the orbital and spin
space, respectively.

Figure 4: (a)-(e) Injection conductivities and (f) trace of the CPGE tensor
for charge-2 WSM obtained using the four-band model Eq.(22). The dashed brown
curve shows the corresponding result for the two-band model (green curve from
Fig. 3) which is close to the $\Delta=0.5$ curve, as expected. We see
significant deviations from perfect CPGE quantization for
$\Delta\lesssim|gu_{z}|$. The
$\sigma^{xxz}_{\text{inj}},\sigma^{zxx}_{\text{inj}},\sigma^{zzz}_{\text{inj}}$
plateaus continue to show up for $\Delta\gg|gu_{z}|$ with heights of about
$\Delta/32$, $-\Delta/32$ and $\Delta(M-2)/32$, respectively.

The $k_{z}$ dependent $\tau_{z}\sigma_{z}$ term produces tilt while
$\Delta(\tau_{x}\sigma_{x}+\tau_{y}\sigma_{y})$ gives rise to the quadratic
band dispersion along $k_{x},k_{y}$. The low-energy Hamiltonian near nodes at
$(0,0,\pm\acos(M-2))$ is,

$\displaystyle\begin{split}\mathcal{H}_{\pm}^{4b}&=t\big{(}k_{x}\tau_{x}+k_{y}\tau_{y}+u_{z}k_{z}\tau_{z}+\Delta\left(\tau_{x}\sigma_{x}+\tau_{y}\sigma_{y}\right)\\\
&\quad+(u_{t}k_{z}+gu_{z})\tau_{z}\sigma_{z}\big{)}-\mu,\end{split}$ (23)

where $u_{z}=\pm\sqrt{(3-M)(M-1)}$, and $u_{t}=g(M-2)$.

We begin with $|\Delta|$ 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
$\frac{|t|(k_{x}^{2}+k_{y}^{2})}{2|\Delta+gu_{z}|}$,
$\frac{-|t|(k_{x}^{2}+k_{y}^{2})}{2|\Delta-gu_{z}|}$ when $k_{z}=0$. We use
the same $g,M,t$ values from before. We note that unlike the two-band case
Eq.(21), the quadratic dispersion now has a dependence on $\Delta$. For
$\Delta=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
$\sim|\Delta|$, the effect of higher bands should be more prominent for
smaller values of $\Delta$.

Results obtained for Eq.(22) are shown in Fig. 4. We find large deviation from
perfect quantization for small gaps as seen in Fig. 4(f). However, for
$\Delta\gg|gu_{z}|$ we do see almost perfect quantization. Also, the plateaus
seen earlier in
$\sigma^{xxz}_{\text{inj}},\sigma^{zxx}_{\text{inj}},\sigma^{zzz}_{\text{inj}}$
continue to show up when $\Delta$ is at least a few times larger than
$|gu_{z}|$ as shown in (a-c). Their heights become dependent on $\Delta$ and
are empirically found to be about $\Delta/32$, $-\Delta/32$, and
$\Delta(M-2)/32$, respectively.

##  V Charge-4 Weyl semimetals

Table 2: Results for JDOS and injection conductivity for the Hamiltonian in Eq.(31) with $c_{1},\omega>0$ (when $c_{1}<0$, replace $c_{1}\rightarrow|c_{1}|$, $\mu\rightarrow-\mu$ and use the $c_{1}>0$ results). Integrals given here are to be evaluated over the region satisfying $x>0,z>0,\left(4c_{2}^{2}\left(x^{2}-xz+z^{2}\right)-c_{1}^{2}\right)\left(-c_{3}^{2}x^{2}\omega z-8c_{2}^{2}c_{1}(x-2z)^{2}+8c_{1}^{3}\right)>0$. Here, we have defined $F_{\frac{2\mu}{\omega}}(x,z)=\Theta\left(-x-z+1+\frac{2\mu}{\omega}\right)\Theta\left(x+z+1-\frac{2\mu}{\omega}\right)\Theta\left(1+\frac{2\mu}{\omega}\right)$. JDOS | $\int_{0}^{1+\frac{2\mu}{\omega}}\differential{z}\int_{0}^{1+\frac{2\mu}{\omega}}\differential{x}\frac{c_{1}\sqrt{\omega}F_{\frac{2\mu}{\omega}}(x,z)}{4\pi^{3}\sqrt{z}\sqrt{\left(4c_{2}^{2}\left(x^{2}-xz+z^{2}\right)-c_{1}^{2}\right)\left(-c_{3}^{2}x^{2}\omega z-8c_{2}^{2}c_{1}(x-2z)^{2}+8c_{1}^{3}\right)}}$
---|---
Injection conductivity |
$\sigma^{xyz}=-\sigma^{xzy}$ | $\int_{0}^{1+\frac{2\mu}{\omega}}\differential{z}\int_{0}^{1+\frac{2\mu}{\omega}}\differential{x}\frac{-i\tau e^{3}\sqrt{3}c_{2}^{2}c_{3}\sqrt{\omega}\left(4c_{2}^{2}c_{3}^{2}x\omega z^{2}\left(x^{2}-z^{2}\right)+c_{3}^{2}c_{1}^{2}x\omega z^{2}-96c_{2}^{4}c_{1}xz^{2}(x-2z)-16c_{2}^{2}c_{1}^{3}(x+z)^{2}+4c_{1}^{5}\right)F_{\frac{2\mu}{\omega}}(x,z)}{2\pi^{2}\hbar^{2}c_{1}^{3}\sqrt{z}\left(24c_{1}c_{2}^{2}-c_{3}^{2}\omega z\right)\sqrt{\left(4c_{2}^{2}\left(x^{2}-xz+z^{2}\right)-c_{1}^{2}\right)\left(-c_{3}^{2}x^{2}\omega z-8c_{2}^{2}c_{1}(x-2z)^{2}+8c_{1}^{3}\right)}}$
$\sigma^{yzx}=-\sigma^{yxz}$
$\sigma^{zxy}=-\sigma^{zyx}$ | $\int_{0}^{1+\frac{2\mu}{\omega}}\differential{z}\int_{0}^{1+\frac{2\mu}{\omega}}\differential{x}\frac{-i\tau e^{3}\sqrt{3}c_{2}^{2}c_{3}\sqrt{\omega}xz^{2}\left(4c_{3}^{2}c_{2}^{2}\omega\left(z^{2}-x^{2}\right)+96c_{1}c_{2}^{4}(x-2z)+c_{1}^{2}c_{3}^{2}\omega\right)F_{\frac{2\mu}{\omega}}(x,z)}{\pi^{2}\hbar^{2}c_{1}^{3}\sqrt{z}\left(24c_{1}c_{2}^{2}-c_{3}^{2}\omega z\right)\sqrt{\left(4c_{2}^{2}\left(x^{2}-xz+z^{2}\right)-c_{1}^{2}\right)\left(-c_{3}^{2}x^{2}\omega z-8c_{2}^{2}c_{1}(x-2z)^{2}+8c_{1}^{3}\right)}}$
$\frac{2\pi}{i\tau e^{3}/\hbar^{2}}\epsilon_{abc}\sigma^{abc}$ | $\int_{0}^{1+\frac{2\mu}{\omega}}\differential{z}\int_{0}^{1+\frac{2\mu}{\omega}}\differential{x}\frac{8\sqrt{3}c_{2}^{2}c_{3}\sqrt{\omega}\left(-2c_{1}^{3}+8c_{1}c_{2}^{2}(x+z)^{2}-c_{3}^{2}x\omega z^{2}\right)F_{\frac{2\mu}{\omega}}(x,z)}{\pi c_{1}\sqrt{z}\left(24c_{1}c_{2}^{2}-c_{3}^{2}\omega z\right)\sqrt{\left(4c_{2}^{2}\left(x^{2}-xz+z^{2}\right)-c_{1}^{2}\right)\left(-c_{3}^{2}x^{2}\omega z-8c_{2}^{2}c_{1}(x-2z)^{2}+8c_{1}^{3}\right)}}$

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],

$\displaystyle\begin{split}\mathcal{H}_{4}&=-2c_{1}(\cos(k_{x})+\cos(k_{y})+\cos(k_{z}))\sigma_{0}\\\
&\quad+2c_{2}\big{(}\sqrt{3}(\cos(k_{y})-\cos(k_{x}))\sigma_{x}\\\
&\quad-(\cos(k_{x})+\cos(k_{y})-2\cos(k_{z}))\sigma_{z}\big{)}\\\
&\quad+c_{3}\sin(k_{x})\sin(k_{y})\sin(k_{z})\sigma_{y}-\widetilde{\mu}\sigma_{0},\end{split}$
(30)

which has nodes of opposite chirality at $(0,0,0)$ and $(\pi,\pi,\pi)$. The
low-energy Hamiltonian near $\Gamma$-point is given by,

$\displaystyle\begin{split}\mathcal{H}_{4}^{\Gamma}&=c_{1}\left(k_{x}^{2}+k_{y}^{2}+k_{z}^{2}\right)\sigma_{0}+c_{2}\Big{(}\sqrt{3}\left(k_{x}^{2}-k_{y}^{2}\right)\sigma_{x}\\\
&\quad+\left(k_{x}^{2}+k_{y}^{2}-2k_{z}^{2}\right)\sigma_{z}\Big{)}+c_{3}k_{x}k_{y}k_{z}\sigma_{y}-\mu\sigma_{0},\end{split}$
(31)

where $\mu=\widetilde{\mu}+6c_{1}$. The chirality of this Weyl point is given
by $\chi=\text{sgn}(c_{3})$. We derive all our results using Eq. (31) with
$c_{1}>0$ (the opposite case is an easy generalization, which we discuss
later). Its eigenvalues are given by

$\displaystyle\begin{split}E^{\Gamma}_{4,\pm}&=c_{1}\left(k_{x}^{2}+k_{y}^{2}+k_{z}^{2}\right)\pm
2|c_{2}|\bigg{(}k_{x}^{4}+k_{y}^{4}+k_{z}^{4}\\\ &\quad-
k_{x}^{2}k_{y}^{2}-k_{y}^{2}k_{z}^{2}-k_{z}^{2}k_{x}^{2}+\left(\frac{c_{3}}{2c_{2}}\right)^{2}k_{x}^{2}k_{y}^{2}k_{z}^{2}\bigg{)}^{\frac{1}{2}}-\mu.\end{split}$
(32)

The presence of the sixth order term $k_{x}^{2}k_{y}^{2}k_{z}^{2}$ above does
not allow us to fully evaluate Eq.(2), Eq.(3), and Eq.(7) 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 (non-zero components) thus obtained are given in Table 2. 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 $\mu=0$, quantization is seen only for $c_{1}/|c_{2}|<1$ starting at a
frequency of $\frac{54c_{1}^{3}}{c_{3}^{2}}$ 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
$\frac{54c_{1}^{3}}{c_{3}^{2}}$ for the latter.

For $\mu<0$, a perfect quantization is again only possible for
$c_{1}/|c_{2}|<1$. When this is the case, the trace becomes non-zero after
$\text{min}\left(\omega_{p},\frac{2|\mu|}{1-\frac{c_{1}}{2|c_{2}|}}\right)$
and reaches $\pm 4$ at
$\text{max}\left(\omega_{p},\frac{2|\mu|}{1-\frac{c_{1}}{|c_{2}|}}\right)$
where $\omega_{p}$ is the unique real positive root of
$(\omega-2|\mu|)^{3}-54\frac{c_{1}^{3}}{c_{3}^{2}}\omega^{2}=0$. For
$1<c_{1}/|c_{2}|<2$, the trace becomes non-zero after
$\text{min}\left(\omega_{p},\frac{2|\mu|}{1-\frac{c_{1}}{2|c_{2}|}}\right)$
while for $2<c_{1}/|c_{2}|$, this happens after $\omega_{p}$. The three cases
are shown in Fig. 5 (d), (e) and (f), respectively.

Figure 5: Trace of the CPGE tensor for a single charge-4 Weyl point Eq.(31)
with $c_{1}>0$, obtained from the numerical evaluation of integrals in Table 2
for different combinations of $\mu$ and $c_{1}/|c_{2}|$ (respective values
shown in the inset). We have taken $c_{1}=0.0665,c_{3}=0.4$. (a)-(c) $\mu=0$,
(d)-(f) $\mu<0$. In both these cases, perfect quantization is seen only when
$c_{1}/|c_{2}|<1$. (g)-(i) $\mu>0$, perfect quantization is guaranteed for
$c_{1}/|c_{2}|<1$ however, unlike previous two cases, it can also be seen for
$c_{1}/|c_{2}|>1$ as long as
$\frac{2\mu}{c_{1}/|c_{2}|-1}>\text{max}\left(\omega_{p},\frac{2\mu}{1+c_{1}/2|c_{2}|}\right)$,
as shown in (h).

When $\mu>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 non-zero after
$\text{min}\left(\omega_{p},\frac{2\mu}{1+\frac{c_{1}}{|c_{2}|}}\right)$ where
$\omega_{p}$ is now the unique real positive root of the cubic equation
$(\omega-2\mu)^{3}+54\frac{c_{1}^{3}}{c_{3}^{2}}\omega^{2}=0$. For
$c_{1}/|c_{2}|<1$, it goes on to reach a saturation value of $\pm 4$ at
$\text{max}\left(\omega_{p},\frac{2\mu}{1+\frac{c_{1}}{2|c_{2}|}}\right)$, as
shown in Fig. 5(g). Interestingly for $\mu>0$, quantization becomes possible
even for $1<c_{1}/|c_{2}|$ provided
$\text{max}\left(\omega_{p},\frac{2\mu}{1+\frac{c_{1}}{2|c_{2}|}}\right)<\frac{2\mu}{\frac{c_{1}}{|c_{2}|}-1}$.
When this condition is met, perfect quantization is seen but only for a finite
window of energy $\omega$ satisfying
$\text{max}\left(\omega_{p},\frac{2\mu}{1+\frac{c_{1}}{2|c_{2}|}}\right)<\omega<\frac{2\mu}{\frac{c_{1}}{|c_{2}|}-1}$
which is shown in Fig. 5(h). The situation when no quantization is possible
for $\mu>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 $\mu$ 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 2. We have
included the JDOS plot at zero doping in Fig. 6(a). As shown in the figure,
the JDOS has a $\sqrt{\omega}$ behavior going towards zero frequency. Note
that the JDOS result from Table 1 also predicts a $\sqrt{\omega}$ 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.(31).

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 $c_{1}<0$ (and a chemical potential $\mu$) by treating it as a $|c_{1}|$
node with chemical potential $-\mu$. With this small but important extension,
our analysis covers all the possible cases.

Figure 6: (a) JDOS for a single charge-4 Weyl point with
$\mu=0,c_{1}=0.0665,c_{2}=0.4668,c_{3}=4$ ($c_{1}/c_{2}<1$). (b) The red and
blue curves capture the CPGE quantization for the two-band tight binding model
Eq.(30) with $\mu=-0.3,0.5$, respectively. Note that
$\mu_{\Gamma}=\mu+6c_{1},\mu_{R}=\mu-6c_{1}$. The corresponding dashed green
curve is obtained by evaluating expressions from Table 2 for each node
separately and then adding the results. We have used the same
$c_{1},c_{2},c_{3}$ from before. (c) Results with $c_{1}/c_{2}>1$, $\mu=-0.3$
(same $c_{1},c_{3}$ as before). The dashed curves show contribution from the
node closer to $E=0$ (obtained using results from Table 2) for each case and
are discontinued after contribution from the other node becomes non-zero.

For completeness, we also compute the CPGE trace using the full tight-binding
model Eq.(30), and the results are shown in Fig. 6(b),(c). Note that in this
model, when going from $\Gamma$ to $R$, we find
$(c_{1},c_{2},c_{3})\rightarrow(-c_{1},-c_{2},-c_{3})$. Since $c_{1}$ turns
negative, when using results from Table 2, we treat the $(c_{1}<0,\mu_{R})$
$R$-node as a $(|c_{1}|,-\mu_{R})$ node. Also, since $c_{3}$ flips sign too,
$\chi_{\Gamma}=-\chi_{R}$, as expected.

In Fig. 6(b), we chose $c_{1}/|c_{2}|<1$ and therefore expect both nodes to
show perfect quantization. Since $\mu\neq 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 2
(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 $\mu=-0.3$. This gives
$\mu_{\Gamma}=0.1$ and $\mu_{R}=-0.7$. Since the node closer to zero energy
falls under the $\mu>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
$\Gamma$-node in both cases. In both Fig. 6(b) and (c), we find excellent
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 ($\mu$). 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 time-reversal broken
systems for the charge-2 case.

Beyond the extremely important 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 low-energy 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_{1}^{3}/c_{3}^{2}$ and $\mu$ 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}\delta(\omega_{21}-\omega)$ term in it at $T=0$ K.

##  VII Acknowledgement

S.C. would like to acknowledge the funding from the National Science
Foundation through the Center for Dynamics and Control of Materials: an NSF
MRSEC under Cooperative Agreement No. DMR- 1720595\. G.A.F. acknowledges
additional support from NSF DMR-2114825

##  Appendix A Analytical expression for shift and injection conductivity
tensors

We work with the low-energy effective Hamiltonian,

with eigenvalues

$\displaystyle E_{n,\pm}$
$\displaystyle=u_{t}k_{z}-\mu\pm\varepsilon_{0}\sqrt{(\tilde{k}_{x}^{2}+\tilde{k}_{y}^{2})^{n}+u_{z}^{2}k_{z}^{2}/\varepsilon_{0}^{2}}.$
(33)

The domain for the integrals in Eq.(2), Eq.(3) is determined by $f_{21}$ and
$\delta(\omega_{21}-\omega)=\delta(2\varepsilon_{0}\sqrt{(\tilde{k}_{x}^{2}+\tilde{k}_{y}^{2})^{n}+u_{z}^{2}k_{z}^{2}/\varepsilon_{0}^{2}}-\omega)$.
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:
$\int\differential{k_{x}}\int\differential{k_{y}}=\int_{+}\differential{k_{x}}\int_{+}\differential{k_{y}}+\int_{+}\differential{k_{x}}\int_{-}\differential{k_{y}}+\int_{-}\differential{k_{x}}\int_{+}\differential{k_{y}}+\int_{-}\differential{k_{x}}\int_{-}\differential{k_{y}}$,
and combine them into a single integral over the first quadrant by making
substitutions $k_{x}=\pm k_{0}\sqrt{x}$ and $k_{y}=\pm k_{0}\sqrt{y}$
depending on the sign (both $x,y>0$). We also put
$k_{z}=\frac{\varepsilon_{0}}{u_{z}}z$. By making
$x\rightarrow\frac{x-y}{\sqrt{2}}$ and $y\rightarrow\frac{x+y}{\sqrt{2}}$, we
rotate the $x-y$ axis counterclockwise by $\pi/4$, and scale $z$ by
$z\rightarrow 2^{n/4}z$. Finally, we let $x\rightarrow x^{2/n}$ to get
$\delta(2^{1+n/4}\varepsilon_{0}\sqrt{x^{2}+z^{2}}-\omega)$. These
transformations also change the integration measure
$\int_{\mathbf{k}}\rightarrow\int\tfrac{2^{1/2+n/4}k_{0}^{2}\varepsilon_{0}x^{2/n-1}\differential{x}\differential{z}\differential{y}}{16\pi^{3}|u_{z}|n\sqrt{x^{2}-y^{2}}}$,
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{\theta},z=r\sin{\theta}$ to obtain
$\delta(r-\omega/(2^{1+n/4}\varepsilon_{0}))/(2^{1+n/4}\varepsilon_{0})$.

Since we are taking the temperature to be zero,
$f_{21}=\Theta(E_{1})-\Theta(E_{2})$ where, $\Theta$ is the Heaviside step
function. Because of the condition put by the delta function, we have

$\displaystyle f_{21}$
$\displaystyle=\Theta\left(u_{t}k_{z}-\mu-\frac{\omega}{2}\right)-\Theta\left(u_{t}k_{z}-\mu+\frac{\omega}{2}\right).$
(34)

Since we have assumed $\omega>0$, the only non-zero value for $f_{21}$ is $-1$
when $\mu-\omega/2<u_{z}k_{z}<\mu+\omega/2$. Using coordinate transformations
from before, this condition becomes

$\displaystyle\frac{2\mu}{\omega}-1<\frac{u_{t}}{u_{z}}\sin{\theta}<\frac{2\mu}{\omega}+1,$
(35)
$\displaystyle\frac{\text{sgn}\left(\frac{u_{t}}{u_{z}}\right)\frac{2\mu}{\omega}-1}{W}<\sin{\theta}<\frac{\text{sgn}\left(\frac{u_{t}}{u_{z}}\right)\frac{2\mu}{\omega}+1}{W},$
(36)

where $W=|u_{t}/u_{z}|$. The definitions for $\theta_{1},\theta_{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
$\mathcal{H}_{n}$, and combine them to obtain analytical expressions for the
JDOS, shift and injection conductivities.

##  Appendix B Higher order terms for charge-2 WSM

Based on the higher order terms appearing in the expansion of Eq.(20) near its
nodes, we look at the effect of including
$(\tfrac{1}{2}(k_{x}^{2}+k_{y}^{2})+u_{m}k_{z}^{2})\sigma_{z}-\tfrac{1}{2}gu_{z}k_{z}^{2}\sigma_{0}$
in Eq.(21) where $u_{m}=\tfrac{1}{2}M-1$. As before, we use a series of
transformations to simplify the Dirac delta constraint. Key steps are as
follows (with $\varepsilon_{0}=1$),

  1. 1.

$k_{x}\rightarrow\pm\sqrt{x}$, $k_{y}\rightarrow\pm\sqrt{y}$,
$k_{z}\rightarrow\frac{u_{z}}{u_{m}}z$.

  2. 2.

$x\rightarrow\frac{1}{\sqrt{2}}(x-y)$, $y\rightarrow\frac{1}{\sqrt{2}}(x+y)$,
$z\rightarrow\frac{1}{2}(\sqrt{z}-1)$.

  3. 3.

$z\rightarrow\frac{\sqrt{40}u_{m}}{u_{z}^{2}}z+1$, integrate out $y$ (from
$-x,x$).

  4. 4.

$x\rightarrow\frac{x}{2\sqrt{5+\sqrt{5}}}-\frac{z}{2\sqrt{5-\sqrt{5}}}$,
$z\rightarrow\frac{x}{2\sqrt{5+\sqrt{5}}}+\frac{z}{2\sqrt{5-\sqrt{5}}}$.

  5. 5.

$x\rightarrow\omega\cos{\theta}$, $z\rightarrow\omega\sin{\theta}$, integrate
from $\theta_{1},\theta_{2}$.

Analytical expressions for JDOS, shift and injection conductivity tensors can
be obtained as before. We still have $f_{21}=-1$, however, the condition that
determines $\theta_{1},\theta_{2}$ becomes,

$\displaystyle\begin{split}2\mu-\omega<\frac{u_{t}u_{z}}{u_{m}}\left(\sqrt{\frac{\omega\sqrt{5}u_{m}\sin(\theta+\beta)}{u_{z}^{2}}+1}-1\right)\\\
-\frac{gu_{z}^{3}}{4u_{m}^{2}}\left(\sqrt{\frac{\omega\sqrt{5}u_{m}\sin(\theta+\beta)}{u_{z}^{2}}+1}-1\right)^{2}<2\mu+\omega,\end{split}$
(37)

where $\beta=\arctan(\varphi-1)$,
$-\frac{\pi}{2}-\arctan(\varphi)\leq\theta\leq\frac{\pi}{2}-\arctan(\varphi)$,
and $\varphi$ is the golden ratio. Allowed values of $\theta$ can be found by
solving this inequality numerically. When solutions turn out to be disjoint
intervals, each interval defines its own $\theta_{1},\theta_{2}$. The
analytical expression is evaluated for each interval and then summed.

##  Appendix C Tilt and Zeeman terms for charge-2 WSM

We can also include additional terms of the form
$A(\tilde{k}_{x}^{2}+\tilde{k}_{y}^{2})\sigma_{0}$ and $B\sigma_{z}$ into
Eq.(LABEL:eq:Hn). These correspond to second-order tilt and Zeeman terms,
respectively. The $B$ term only shifts the origin along $k_{z}$, modifying the
$k_{z}\rightarrow z$ transformation to
$k_{z}=\frac{\varepsilon_{0}z-B}{u_{z}}$. It not difficult to see that these
terms only affect $f_{21}$,

$\displaystyle\begin{split}f_{21}&=\Theta\left(\frac{A}{\varepsilon_{0}}\frac{\omega\cos{\theta}}{2}+\frac{u_{t}}{u_{z}}\frac{\omega\sin{\theta}}{2}-\widetilde{\mu}-\frac{\omega}{2}\right)\\\
&\qquad-\Theta\left(\frac{A}{\varepsilon_{0}}\frac{\omega\cos{\theta}}{2}+\frac{u_{t}}{u_{z}}\frac{\omega\sin{\theta}}{2}-\widetilde{\mu}+\frac{\omega}{2}\right),\end{split}$
(38)

where $\widetilde{\mu}=\mu+Bu_{t}/u_{z}$. Since $\omega>0$, we have
$f_{21}=-1$ when,

$\displaystyle\frac{2\widetilde{\mu}}{\omega}-1<\frac{A}{\varepsilon_{0}}\cos{\theta}+\frac{u_{t}}{u_{z}}\sin{\theta}<\frac{2\widetilde{\mu}}{\omega}+1.$
(39)

By defining
$\widetilde{W}=\sqrt{A^{2}/\varepsilon_{0}^{2}+u_{t}^{2}/u_{z}^{2}}$,
$\sin{\phi}=\frac{A/\varepsilon_{0}}{\widetilde{W}}$, and
$\cos{\phi}=\frac{|u_{t}/u_{z}|}{\widetilde{W}}$, we obtain with
$-\pi/2\leq\theta,\phi\leq\pi/2$. Also, we define
$\alpha=\text{sgn}\left(\tfrac{u_{t}}{u_{z}}\right)\phi$ and
$\widetilde{\varphi}_{p}=\frac{1}{\widetilde{W}}\left(\text{sgn}\left(\frac{u_{t}}{u_{z}}\right)\frac{2\widetilde{\mu}}{\omega}+(-1)^{p}\right)$
with $p=1,2$. The inequality becomes
$\widetilde{\varphi}_{1}<\sin(\theta+\alpha)<\widetilde{\varphi}_{2}$ which
can be solved for the minimum ($\widetilde{\theta_{1}}$) and maximum
($\widetilde{\theta_{2}}$) allowed values of $\theta$. These can be obtained
as the left and right end points of the intervals,

$\displaystyle(\widetilde{\theta_{1}},\widetilde{\theta_{2}})$
$\displaystyle=\begin{dcases}(\pi/2,\pi/2),\quad
1<\widetilde{\varphi}_{1},1<\widetilde{\varphi}_{2}\\\
(-\pi/2,-\pi/2),\quad\widetilde{\varphi}_{1}<-1,\widetilde{\varphi}_{2}<-1\\\
I_{0},\quad\widetilde{\varphi}_{1}<-1,1<\widetilde{\varphi}_{2}\\\ I_{0}\cap
I_{2},\quad 0<\widetilde{\varphi}_{1}<1,1<\widetilde{\varphi}_{2}\\\ I_{0}\cap
I_{1},\quad-1<\widetilde{\varphi}_{1}<0,1<\widetilde{\varphi}_{2}\\\ I_{0}\cap
I_{3},\quad\widetilde{\varphi}_{1}<-1,-1<\widetilde{\varphi}_{2}<0\\\
I_{0}\cap I_{4},\quad\widetilde{\varphi}_{1}<-1,0<\widetilde{\varphi}_{2}<1\\\
I_{0}\cap I_{1}\cap
I_{4},\quad-1<\widetilde{\varphi}_{1}<0<\widetilde{\varphi}_{2}<1\\\ I_{0}\cap
I_{1}\cap I_{3},\quad-1<\widetilde{\varphi}_{1}<\widetilde{\varphi}_{2}<0\\\
I_{0}\cap I_{2}\cap I_{4},\quad
0<\widetilde{\varphi}_{1}<\widetilde{\varphi}_{2}<1\\\ \end{dcases}$ (41)

where,

$\displaystyle I_{0}$ $\displaystyle=(-\pi/2,\pi/2),$ $\displaystyle I_{1}$
$\displaystyle=(-\pi-\alpha,-\pi-\alpha-\arcsin{\widetilde{\varphi}_{1}})\cup(-\alpha+\arcsin{\widetilde{\varphi}_{1}},\pi-\alpha),$
$\displaystyle I_{2}$
$\displaystyle=(-\alpha+\arcsin{\widetilde{\varphi}_{1}},\pi-\alpha-\arcsin{\widetilde{\varphi}_{1}}),$
$\displaystyle I_{3}$
$\displaystyle=(-\pi-\alpha-\arcsin{\widetilde{\varphi}_{2}},-\alpha+\arcsin{\widetilde{\varphi}_{2}}),$
$\displaystyle I_{4}$
$\displaystyle=(-\pi-\alpha,-\alpha+\arcsin{\widetilde{\varphi}_{2}})\cup(\pi-\alpha-\arcsin{\widetilde{\varphi}_{2}},\pi-\alpha).$

##  Appendix D Sign changing of $\sigma^{xyz}_{\mathrm{shift}}$ for $W>2$

For $W<2$, we found that the sign change occurred at $\omega=2|\mu|$ when one
angle was $\pm\pi/2$ and the other zero. However, when $W>2$,
$(\theta_{2},\theta_{1})$ cannot take either $(0,-\pi/2)$, or $(\pi/2,0)$.
Finding the point of sign change now requires us to seek other solutions of
$\sin{\theta_{2}}\cos^{2}{\theta_{2}}-\sin{\theta_{1}}\cos^{2}{\theta_{1}}=0$.
Converting cosine into sine we get
$(\sin{\theta_{2}}-\sin{\theta_{1}})(\sin^{2}{\theta_{2}}+\sin{\theta_{2}}\sin{\theta_{1}}+\sin^{2}{\theta_{1}}-1)=0$.
Let’s look for solutions other than $\theta_{2}=\theta_{1}$, $(0,-\pi/2)$, and
$(\pi/2,0)$. We can solve for $(\sin{\theta_{2}},\sin{\theta_{1}})$ to get,

$\displaystyle(\sin{\theta_{2}},\sin{\theta_{1}})$
$\displaystyle=\begin{dcases}\left(\tfrac{-x-\sqrt{4-3x^{2}}}{2},x\right),\,-1<x<\tfrac{-1}{\sqrt{3}}\\\
\left(x,\tfrac{-x-\sqrt{4-3x^{2}}}{2}\right),\,\tfrac{-1}{\sqrt{3}}<x<0\\\
\left(\tfrac{-x+\sqrt{4-3x^{2}}}{2},x\right),\,0<x<\tfrac{1}{\sqrt{3}}\\\
\left(x,\tfrac{-x+\sqrt{4-3x^{2}}}{2}\right),\,\tfrac{1}{\sqrt{3}}<x<1\end{dcases}$
(42)

Using definitions of $\theta_{1},\theta_{2}$, we solve for $\omega$ by
eliminating $x$ to obtain $\omega=2|\mu|\sqrt{\frac{3}{W^{2}-1}}$. Note that
for $W=2$, this gives $\omega=2|\mu|$ as expected.

##  Appendix E Analytical results for charge-4 WSM

The delta function constraint $\delta(\omega-\omega_{21})$ translates to

$\displaystyle\begin{split}&4|c_{2}|\bigg{(}k_{x}^{4}+k_{y}^{4}+k_{z}^{4}-k_{x}^{2}k_{y}^{2}-k_{y}^{2}k_{z}^{2}-k_{z}^{2}k_{x}^{2}\\\
&\quad\hskip
85.35826pt+\left(\frac{c_{3}}{2c_{2}}\right)^{2}k_{x}^{2}k_{y}^{2}k_{z}^{2}\bigg{)}^{\frac{1}{2}}=\omega.\end{split}$
(43)

To simplify this, we use the following transformations,

  1. 1.

$k_{x}\rightarrow\pm\sqrt{x}$, $k_{y}\rightarrow\pm\sqrt{y}$,
$k_{z}\rightarrow\pm\sqrt{z}$ (reduce the integral to $x,y,z>0$ octant).

  2. 2.

$x\rightarrow\frac{1}{\sqrt{2}}(x-y)$, $y\rightarrow\frac{1}{\sqrt{2}}(x+y)$.

  3. 3.

Integrate out $y$ (from $-x,x$). To do this, we need to find roots of the
equation
$\omega=\left(2c_{3}^{2}z\left(x^{2}-y^{2}\right)+8c_{2}^{2}\left(\left(x-\sqrt{2}z\right)^{2}+3y^{2}\right)\right)^{\frac{1}{2}}$.
The condition for existence of real roots satisfying $-x\leq y\leq x$ is given
by (post step 4 substitution).

$\displaystyle\begin{split}&\left(4c_{2}^{2}\left(x^{2}-xz+z^{2}\right)-c_{1}^{2}\right)\big{(}-c_{3}^{2}x^{2}\omega
z\\\ &\hskip
71.13188pt-8c_{2}^{2}c_{1}(x-2z)^{2}+8c_{1}^{3}\big{)}>0.\end{split}$ (44)

  4. 4.

$x\rightarrow\frac{\omega}{2\sqrt{2}c_{1}}x$,
$z\rightarrow\frac{\omega}{2c_{1}}z$ .

Using these transformations along with the eigenvalues and normalized
eigenfunctions of Eq.(31), we simply Eq.(2), Eq.(3), Eq.(7) to obtain the
expressions shown in Table 2 (the shift conductivities are zero).

Figure 7: The green area is the region defined by Eq.(44) (we have used
$c_{1}=0.0665,c_{2}=0.4668,c_{3}=4$, $\omega=0.4$, but the features shown here
are quite general). The blue, orange and magenta curves correspond to
$-c_{3}^{2}x^{2}\omega z-8c_{2}^{2}c_{1}(x-2z)^{2}+8c_{1}^{3}=0$ for
$\omega=0.1,0.4$ and $4.0$, respectively. The green curve represents the
ellipse $4c_{2}^{2}\left(x^{2}-xz+z^{2}\right)-c_{1}^{2}=0$. The blue dots are
placed at $\left(\frac{2^{1/3}c_{1}}{|c_{3}|^{2/3}\omega^{1/3}},\frac{2\times
2^{1/3}c_{1}}{|c_{3}|^{2/3}\omega^{1/3}}\right)$.

Behavior of $\frac{2\pi}{i\tau e^{3}/\hbar^{2}}\epsilon_{abc}\sigma^{abc}$ is
determined by the interplay between conditions set by
$\Theta\left(-x-z+1+\frac{2\mu}{\omega}\right)$,
$\Theta\left(x+z+1-\frac{2\mu}{\omega}\right)$ and Eq.(44). Since $c_{1}>0$ by
choice and $c_{2},c_{3}$ appear only as their squares, the analysis of region
defined by Eq.(44) becomes quite general. To understand this, let us focus on
the curves $4c_{2}^{2}\left(x^{2}-xz+z^{2}\right)-c_{1}^{2}=0$ and
$-c_{3}^{2}x^{2}\omega z-8c_{2}^{2}c_{1}(x-2z)^{2}+8c_{1}^{3}=0$ for $x,z>0$.
They intersect the $x$-axis at $x=c_{1}/2|c_{2}|$ and $x=c_{1}/|c_{2}|$,
respectively, but cross the $z$-axis together at $z=c_{1}/2|c_{2}|$
(intercepts are independent of $\omega$). Tangents to these curves with slope
$-1$ are important. For the ellipse this happens at
$(z=c_{1}/2|c_{2}|,x=c_{1}/2|c_{2}|)$, the tangent has equation
$x+z=c_{1}/|c_{2}|$. For the second curve we have several cases. For
$\omega<48\frac{|c_{2}|^{3}}{c_{3}^{2}}$, there is only one such tangent at
($z=\frac{2^{1/3}c_{1}}{|c_{3}|^{2/3}\omega^{1/3}},x=\frac{2\times
2^{1/3}c_{1}}{|c_{3}|^{2/3}\omega^{1/3}}$), described by
$x+z=\frac{54^{1/3}c_{1}}{|c_{3}|^{2/3}\omega^{1/3}}$. For larger $\omega$,
there is another tangent with slope $-1$, but its presence is of no
consequence to our analysis. The important thing to note is that
$x+z=\frac{54^{1/3}c_{1}}{|c_{3}|^{2/3}\omega^{1/3}}$ is completely sandwiched
between $x+z=c_{1}/|c_{2}|$ and $x+z=c_{1}/2|c_{2}|$ for
$54\frac{|c_{2}|^{3}}{c_{3}^{2}}<\omega<8\times
54\frac{|c_{2}|^{3}}{c_{3}^{2}}$. These features are illustrated in Fig. 7.
With these key observations in mind, we now analyze the $\mu=0$, $\mu<0$, and
$\mu>0$ cases separately.

For $\mu=0$, the theta function constraints reduce to $x+z<1$. When
$c_{1}/|c_{2}|<1$, the CPGE trace is non-zero for any finite $\omega$ and
becomes $\pm 4$ after $\frac{54c_{1}^{3}}{c_{3}^{2}}$. When $c_{1}/|c_{2}|>1$,
some portion of Eq.(44) is always left out and we do not get perfect
quantization. For $1<\frac{c_{1}}{|c_{2}|}<2$, the trace is non-zero for any
finite $\omega$ where as for $\frac{c_{1}}{|c_{2}|}>2$, this happens only
after $\frac{54c_{1}^{3}}{c_{3}^{2}}$.

For $\mu<0$, the condition set by
$\Theta\left(x+z+1-\frac{2\mu}{\omega}\right)$ is always satisfied where as
$\Theta\left(-x-z+1+\frac{2\mu}{\omega}\right)$ requires
$x+z<1-\frac{2|\mu|}{\omega}$. An important thing to note here is that the
term $1-\frac{2|\mu|}{\omega}\in(-\infty,1)$. We are only interested when it
lies in $(0,1)$ which happens for $\omega>2|\mu|$. Since it can ever only
reach 1, full overlap with region Eq.(44) 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-\frac{2|\mu|}{\omega}$ and Eq.(44) is
determined by solutions to equations
$1-\frac{2|\mu|}{\omega}=\frac{c_{1}}{|c_{2}|}$,
$1-\frac{2|\mu|}{\omega}=\frac{c_{1}}{2|c_{2}|}$, and
$1-\frac{2|\mu|}{\omega}=\frac{54^{1/3}c_{1}}{|c_{3}|^{2/3}\omega^{1/3}}$. The
last equation can be rewritten as
$(\omega-2|\mu|)^{3}-54\frac{c_{1}^{3}}{c_{3}^{2}}\omega^{2}=0$. This cubic
equation never has three real roots. Since the product of its roots is
$8|\mu|^{3}>0$, the only real root, $\omega_{p}$, is always positive. Note
that $\omega_{p}>2|\mu|$. The CPGE trace becomes non-zero after
$\text{min}\left(\omega_{p},\frac{2|\mu|}{1-\frac{c_{1}}{2|c_{2}|}}\right)$,
and saturates to $\pm 4$ after
$\text{max}\left(\omega_{p},\frac{2|\mu|}{1-\frac{c_{1}}{|c_{2}|}}\right)$.
When $1<c_{1}/|c_{2}|<2$, the trace is non-zero after
$\text{min}\left(\omega_{p},\frac{2|\mu|}{1-\frac{c_{1}}{2|c_{2}|}}\right)$
where as for $c_{1}/|c_{2}|>2$, this happens after $\omega_{p}$ (it never
reaches $\pm 4$ in either case).

For $\mu>0$, the possibilities become even more interesting.
$\Theta\left(-x-z+1+\frac{2\mu}{\omega}\right)\Theta\left(x+z+1-\frac{2\mu}{\omega}\right)$
sets bounds on the integration region, requiring
$\frac{2\mu}{\omega}-1<x+z<\frac{2\mu}{\omega}+1$. The term
$1+\frac{2\mu}{\omega}\in(1,\infty)$, which means that if $c_{1}/|c_{2}|>1$, a
portion of Eq.(44) will necessarily be left out for
$\omega>\frac{2\mu}{\frac{c_{1}}{|c_{2}|}-1}$ (perfect quantization still
possible for smaller energies). Now, the solutions to equations
$\frac{2\mu}{\omega}-1=\frac{c_{1}}{|c_{2}|}$,
$\frac{2\mu}{\omega}-1=\frac{c_{1}}{2|c_{2}|}$, and
$\frac{2\mu}{\omega}-1=\frac{54^{1/3}c_{1}}{|c_{3}|^{2/3}\omega^{1/3}}$ become
crucial in determining the amount of region Eq.(44) available for integration.
The last equation can be rewritten as the cubic equation
$(\omega-2\mu)^{3}+54\frac{c_{1}^{3}}{c_{3}^{2}}\omega^{2}=0$. The product of
its roots is $8\mu^{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 $\mu<4\frac{c_{1}^{3}}{c_{3}^{2}}$ where
as for all roots to be positive, $\mu>9\frac{c_{1}^{3}}{c_{3}^{2}}$. Thus, we
always get exactly one positive root, $\omega_{p}$. Note that
$\omega_{p}<2\mu$ in this case. The CPGE trace becomes non-zero after
$\text{min}\left(\omega_{p},\frac{2\mu}{1+\frac{c_{1}}{|c_{2}|}}\right)$. For
$c_{1}/|c_{2}|<1$, it goes on to reach a saturation value of $\pm 4$ after
$\text{max}\left(\omega_{p},\frac{2\mu}{1+\frac{c_{1}}{2|c_{2}|}}\right)$.
When $c_{1}/|c_{2}|>1$, we see quantization for
$\text{max}\left(\omega_{p},\frac{2\mu}{1+\frac{c_{1}}{2|c_{2}|}}\right)<\omega<\frac{2\mu}{\frac{c_{1}}{|c_{2}|}-1}$.
Perfect quantization is not possible when
$\frac{2\mu}{\frac{c_{1}}{|c_{2}|}-1}<\text{max}\left(\omega_{p},\frac{2\mu}{1+\frac{c_{1}}{2|c_{2}|}}\right)$.

## References

  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).
  * Haldane [1988] F. D. M. Haldane, Model for a quantum hall effect without landau levels: Condensed-matter realization of the ”parity anomaly”, Phys. Rev. Lett. 61, 2015 (1988).
  * Xie  _et al._ [2020] F. Xie, Z. Song, B. Lian, and B. A. Bernevig, Topology-bounded superfluid weight in twisted bilayer graphene, Phys. Rev. Lett. 124, 167002 (2020).
  * Hu  _et al._ [2022] X. Hu, T. Hyart, D. I. Pikulin, and E. Rossi, Quantum-metric-enabled exciton condensate in double twisted bilayer graphene, Phys. Rev. B 105, L140506 (2022).
  * Xiao  _et al._ [2006] D. Xiao, Y. Yao, Z. Fang, and Q. Niu, Berry-phase effect in anomalous thermoelectric transport, Phys. Rev. Lett. 97, 026603 (2006).
  * Ma  _et al._ [2021] Q. Ma, A. G. Grushin, and K. S. Burch, Topology and geometry under the nonlinear electromagnetic spotlight, Nature materials 20, 1601 (2021).
  * Orenstein  _et al._ [2021] J. Orenstein, J. Moore, T. Morimoto, D. Torchinsky, J. Harter, and D. Hsieh, Topology and symmetry of quantum materials via nonlinear optical responses, Annual Review of Condensed Matter Physics 12, 247 (2021).
  * Sipe and Shkrebtii [2000] J. E. Sipe and A. I. Shkrebtii, Second-order optical response in semiconductors, Phys. Rev. B 61, 5337 (2000).
  * Dai and Rappe [2023] Z. Dai and A. M. Rappe, Recent progress in the theory of bulk photovoltaic effect, Chemical Physics Reviews 4, 011303 (2023).
  * Young and Rappe [2012] S. M. Young and A. M. Rappe, First principles calculation of the shift current photovoltaic effect in ferroelectrics, Phys. Rev. Lett. 109, 116601 (2012).
  * Shockley and Queisser [1961] W. Shockley and H. J. Queisser, Detailed balance limit of efficiency of p-n junction solar cells, Journal of applied physics 32, 510 (1961).
  * Fregoso  _et al._ [2017] B. M. Fregoso, T. Morimoto, and J. E. Moore, Quantitative relationship between polarization differences and the zone-averaged shift photocurrent, Phys. Rev. B 96, 075421 (2017).
  * Ahn  _et al._ [2020] J. Ahn, G.-Y. Guo, and N. Nagaosa, Low-frequency divergence and quantum geometry of the bulk photovoltaic effect in topological semimetals, Phys. Rev. X 10, 041041 (2020).
  * Holder  _et al._ [2020] T. Holder, D. Kaplan, and B. Yan, Consequences of time-reversal-symmetry breaking in the light-matter interaction: Berry curvature, quantum metric, and diabatic motion, Phys. Rev. Research 2, 033100 (2020).
  * Ahn  _et al._ [2022] J. Ahn, G.-Y. Guo, N. Nagaosa, and A. Vishwanath, Riemannian geometry of resonant optical responses, Nature Physics 18, 290 (2022).
  * Morimoto and Nagaosa [2016] T. Morimoto and N. Nagaosa, Topological nature of nonlinear optical effects in solids, Science advances 2, e1501524 (2016).
  * Chan  _et al._ [2017] C.-K. Chan, N. H. Lindner, G. Refael, and P. A. Lee, Photocurrents in weyl semimetals, Phys. Rev. B 95, 041104 (2017).
  * Nagaosa  _et al._ [2020] N. Nagaosa, T. Morimoto, and Y. Tokura, Transport, magnetic and optical properties of weyl materials, Nature Reviews Materials 5, 621 (2020).
  * de Juan  _et al._ [2017] F. de Juan, A. G. Grushin, T. Morimoto, and J. E. Moore, Quantized circular photogalvanic effect in weyl semimetals, Nature communications 8, 1 (2017).
  * Ma  _et al._ [2017] Q. Ma, S.-Y. Xu, C.-K. Chan, C.-L. Zhang, G. Chang, Y. Lin, W. Xie, T. Palacios, H. Lin, S. Jia, _et al._ , Direct optical detection of weyl fermion chirality in a topological semimetal, Nature Physics 13, 842 (2017).
  * Rees  _et al._ [2020] D. Rees, K. Manna, B. Lu, T. Morimoto, H. Borrmann, C. Felser, J. Moore, D. H. Torchinsky, and J. Orenstein, Helicity-dependent photocurrents in the chiral weyl semimetal rhsi, Science advances 6, eaba0509 (2020).
  * Ma  _et al._ [2019] J. Ma, Q. Gu, Y. Liu, J. Lai, P. Yu, X. Zhuo, Z. Liu, J.-H. Chen, J. Feng, and D. Sun, Nonlinear photoresponse of type-ii weyl semimetals, Nature materials 18, 476 (2019).
  * Sadhukhan and Nag [2021a] B. Sadhukhan and T. Nag, Role of time reversal symmetry and tilting in circular photogalvanic responses, Phys. Rev. B 103, 144308 (2021a).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * Burkov [2018] A. Burkov, Weyl metals, Annual Review of Condensed Matter Physics 9, 359 (2018).
  * Chang  _et al._ [2017] G. Chang, S.-Y. Xu, B. J. Wieder, D. S. Sanchez, S.-M. Huang, I. Belopolski, T.-R. Chang, S. Zhang, A. Bansil, H. Lin, and M. Z. Hasan, Unconventional chiral fermions and large topological fermi arcs in rhsi, Phys. Rev. Lett. 119, 206401 (2017).
  * Chang  _et al._ [2018] G. Chang, B. J. Wieder, F. Schindler, D. S. Sanchez, I. Belopolski, S.-M. Huang, B. Singh, D. Wu, T.-R. Chang, T. Neupert, _et al._ , Topological quantum properties of chiral crystals, Nature materials 17, 978 (2018).
  * Fang  _et al._ [2012] C. Fang, M. J. Gilbert, X. Dai, and B. A. Bernevig, Multi-weyl topological semimetals stabilized by point group symmetry, Phys. Rev. Lett. 108, 266802 (2012).
  * Huang  _et al._ [2016] S.-M. Huang, S.-Y. Xu, I. Belopolski, C.-C. Lee, G. Chang, T.-R. Chang, B. Wang, N. Alidoust, G. Bian, M. Neupane, _et al._ , New type of weyl semimetal with quadratic double weyl fermions, Proceedings of the National Academy of Sciences 113, 1180 (2016).
  * Singh  _et al._ [2018] B. Singh, G. Chang, T.-R. Chang, S.-M. Huang, C. Su, M.-C. Lin, H. Lin, and A. Bansil, Tunable double-weyl fermion semimetal state in the srsi2 materials class, Scientific reports 8, 10540 (2018).
  * Zhu  _et al._ [2018] Z. Zhu, Y. Liu, Z.-M. Yu, S.-S. Wang, Y. X. Zhao, Y. Feng, X.-L. Sheng, and S. A. Yang, Quadratic contact point semimetal: Theory and material realization, Phys. Rev. B 98, 125104 (2018).
  * Zhang  _et al._ [2018a] D. Zhang, H. Wang, J. Ruan, G. Yao, and H. Zhang, Engineering topological phases in the luttinger semimetal $\alpha$-sn, Phys. Rev. B 97, 195139 (2018a).
  * Ghorashi  _et al._ [2018] S. A. A. Ghorashi, P. Hosur, and C.-S. Ting, Irradiated three-dimensional luttinger semimetal: A factory for engineering weyl semimetals, Phys. Rev. B 97, 205402 (2018).
  * Liu and Zunger [2017] Q. Liu and A. Zunger, Predicted realization of cubic dirac fermion in quasi-one-dimensional transition-metal monochalcogenides, Physical Review X 7, 021019 (2017).
  * Hübener  _et al._ [2017] H. Hübener, M. A. Sentef, U. De Giovannini, A. F. Kemper, and A. Rubio, Creating stable floquet–weyl semimetals by laser-driving of 3d dirac materials, Nature communications 8, 13940 (2017).
  * Chen and Fiete [2016] Q. Chen and G. A. Fiete, Thermoelectric transport in double-weyl semimetals, Phys. Rev. B 93, 155125 (2016).
  * Fu and Wang [2022a] L. X. Fu and C. M. Wang, Thermoelectric transport of multi-weyl semimetals in the quantum limit, Phys. Rev. B 105, 035201 (2022a).
  * Menon and Basu [2020] A. Menon and B. Basu, Anomalous hall transport in tilted multi-weyl semimetals, Journal of Physics: Condensed Matter 33, 045602 (2020).
  * Sinha and Sengupta [2019] D. Sinha and K. Sengupta, Transport across junctions of a weyl and a multi-weyl semimetal, Physical Review B 99, 075153 (2019).
  * Kulikov  _et al._ [2020] K. Kulikov, D. Sinha, Y. M. Shukrinov, and K. Sengupta, Josephson junctions of weyl and multi-weyl semimetals, Phys. Rev. B 101, 075110 (2020).
  * Park  _et al._ [2017] S. Park, S. Woo, E. J. Mele, and H. Min, Semiclassical boltzmann transport theory for multi-weyl semimetals, Phys. Rev. B 95, 161113 (2017).
  * Dantas  _et al._ [2018] R. Dantas, F. Pena-Benitez, B. Roy, and P. Surówka, Magnetotransport in multi-weyl semimetals: A kinetic theory approach, Journal of High Energy Physics 2018, 1 (2018).
  * Nag and Nandy [2020] T. Nag and S. Nandy, Magneto-transport phenomena of type-i multi-weyl semimetals in co-planar setups, Journal of Physics: Condensed Matter 33, 075504 (2020).
  * Huang  _et al._ [2017] Z.-M. Huang, J. Zhou, and S.-Q. Shen, Topological responses from chiral anomaly in multi-weyl semimetals, Phys. Rev. B 96, 085201 (2017).
  * Fu and Wang [2022b] L. Fu and C. Wang, Thermoelectric transport of multi-weyl semimetals in the quantum limit, Physical Review B 105, 035201 (2022b).
  * Gorbar  _et al._ [2017] E. Gorbar, V. Miransky, I. Shovkovy, and P. Sukhachov, Anomalous thermoelectric phenomena in lattice models of multi-weyl semimetals, Physical Review B 96, 155138 (2017).
  * Roy and Narayan [2022] S. Roy and A. Narayan, Non-linear hall effect in multi-weyl semimetals, Journal of Physics: Condensed Matter 34, 385301 (2022).
  * Zeng  _et al._ [2023] C. Zeng, S. Nandy, P. Liu, S. Tewari, and Y. Yao, Quantum oscillations as a robust fingerprint of chiral anomaly in nonlinear response in weyl semimetals, Phys. Rev. B 107, L081107 (2023).
  * Nag and Kennes [2022] T. Nag and D. M. Kennes, Distinct signatures of particle-hole symmetry breaking in transport coefficients for generic multi-weyl semimetals, Phys. Rev. B 105, 214307 (2022).
  * Bouhlal  _et al._ [2022] A. Bouhlal, A. Abbout, A. Jellal, H. Bahlouli, and M. Vogl, Tunneling phase diagrams in anisotropic multi-weyl semimetals, Annalen der Physik 534, 2200267 (2022).
  * Bouhlal  _et al._ [2021] A. Bouhlal, A. Jellal, H. Bahlouli, and M. Vogl, Tunneling in an anisotropic cubic dirac semi-metal, Annals of Physics 432, 168563 (2021).
  * Ahn  _et al._ [2017] S. Ahn, E. Mele, and H. Min, Optical conductivity of multi-weyl semimetals, Physical Review B 95, 161112 (2017).
  * Sun and Wang [2017] Y. Sun and A.-M. Wang, Magneto-optical conductivity of double weyl semimetals, Phys. Rev. B 96, 085147 (2017).
  * Mukherjee and Carbotte [2018] S. P. Mukherjee and J. P. Carbotte, Doping and tilting on optics in noncentrosymmetric multi-weyl semimetals, Phys. Rev. B 97, 045150 (2018).
  * Yadav  _et al._ [2023] S. Yadav, S. Sekh, and I. Mandal, Magneto-optical conductivity in the type-i and type-ii phases of weyl/multi-weyl semimetals, Physica B: Condensed Matter , 414765 (2023).
  * Ghosh  _et al._ [2022] S. Ghosh, A. Sahoo, and S. Nandy, Theoretical investigations on kerr and faraday rotations in topological multi-weyl semimetals, arXiv:2209.11217 (2022).
  * Das  _et al._ [2022] S. Das, K. Das, and A. Agarwal, Nonlinear magnetoconductivity in weyl and multi-weyl semimetals in quantizing magnetic field, Phys. Rev. B 105, 235408 (2022).
  * Sadhukhan and Nag [2021b] B. Sadhukhan and T. Nag, Electronic structure and unconventional nonlinear response in double weyl semimetal $\mathrm{Sr}{\mathrm{si}}_{2}$, Phys. Rev. B 104, 245122 (2021b).
  * Cook  _et al._ [2017] A. M. Cook, B. M. Fregoso, F. De Juan, S. Coh, and J. E. Moore, Design principles for shift current photovoltaics, Nature communications 8, 1 (2017).
  * Zhang  _et al._ [2018b] Y. Zhang, H. Ishizuka, J. van den Brink, C. Felser, B. Yan, and N. Nagaosa, Photogalvanic effect in weyl semimetals from first principles, Phys. Rev. B 97, 241118 (2018b).
  * Zhang  _et al._ [2020] T. Zhang, R. Takahashi, C. Fang, and S. Murakami, Twofold quadruple weyl nodes in chiral cubic crystals, Physical Review B 102, 125148 (2020).
  * Yu  _et al._ [2022] Z.-M. Yu, Z. Zhang, G.-B. Liu, W. Wu, X.-P. Li, R.-W. Zhang, S. A. Yang, and Y. Yao, Encyclopedia of emergent particles in three-dimensional crystals, Science Bulletin 67, 375 (2022).
  * Cui  _et al._ [2021] C. Cui, X.-P. Li, D.-S. Ma, Z.-M. Yu, and Y. Yao, Charge-four weyl point: Minimum lattice model and chirality-dependent properties, Phys. Rev. B 104, 075115 (2021).
  * de Juan  _et al._ [2020] F. de Juan, Y. Zhang, T. Morimoto, Y. Sun, J. E. Moore, and A. G. Grushin, Difference frequency generation in topological semimetals, Physical Review Research 2, 012017 (2020).
  * Yang  _et al._ [2017] X. Yang, K. Burch, and Y. Ran, Divergent bulk photovoltaic effect in weyl semimetals, arXiv preprint arXiv:1712.09363 (2017).
  * Mandal [2020] I. Mandal, Effect of interactions on the quantization of the chiral photocurrent for double-weyl semimetals, Symmetry 12, 919 (2020).