# Diverging shift current responses in the gapless limit of two-dimensional
systems

Hiroki Yoshida  Department of Physics, Tokyo Institute of Technology, 2-12-1
Ookayama, Meguro-ku, Tokyo 152-8551, Japan  Shuichi Murakami  Department of
Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo
152-8551, Japan

###### Abstract

The shift current responses of two-dimensional systems in the gapless limit
are investigated. As the energy gap becomes smaller, the interband transition
probability becomes larger at the band edges. We found the divergence of the
shift current response in a manner proportional to the inverse of the gap size
when the system becomes a two-dimensional Weyl semimetal in the gapless limit.
This behavior is different from that for electric polarization, which has a
jump across the gapless case. It means that the quantitative relationship
between electric polarizations and shift currents is broken in this gapless
limit.

##  I introduction

The bulk photovoltaic effect (BPVE) is a DC current response to the uniformly
irradiated light in non-centrosymmetric systems [1, 2, 3, 4, 5, 6, 7, 8]. The
BPVE has been actively investigated as it is a candidate for an alternative
way to convert the light energy to the electric power. Since currents
generated by the application of an external field oscillate in phase with that
field, this BPVE is inherently a non-linear optical effect. Leading
contributions to the BPVE are the second-order optical responses caused by the
interband excitations of Bloch electrons. In the clean limit, generated
second-order currents can be divided into two parts called injection currents
[9, 10] and shift currents [11, 12]. Injection currents and shift currents are
caused by the change of the velocity and the position of electrons by the
interband transition, respectively [5].

Injection and shift currents have been attracting much attention because they
are related to the quantum geometry [13, 5, 10, 14, 15, 16, 17, 18, 19]. It is
known that the injection current has a direct connection with the Berry
curvature [5] and even quantization of the response occurs in Weyl semimetals
[20]. On the other hand, the shift current response is much more complicated
and there are no direct connections between shift currents and geometric
quantities such as Berry curvatures and Berry connections. Under some
assumptions, it is found that the shift current is directly related with the
Berry connection. As the integration of the Berry connection over the
Brillouin zone is the electric polarization [21, 22, 23], shift current can be
described as a difference of the electric polarizations between conduction and
valence bands [7]. There are many previous studies that have investigated
shift currents in topological materials [9, 24, 25, 26, 27], and quantization
of shift currents by multi-gap topology has also been predicted [28].

In this paper, we investigate the shift current response of two-dimensional
systems when their band gap approaches zero. In our previous studies, we found
that the electric polarization of a two-dimensional system has a jump when the
system changes across a two-dimensional Weyl semimetal phase [29, 30].
Surprisingly, we found that this jump is solely determined by the Weyl dipole,
described by the relative positions and monopole charges of Weyl points in the
reciprocal space. Since the shift current is related with the electric
polarization, we naively expect a jump or other singular behaviors of the
shift current response across the two-dimensional Weyl semimetal phase. As a
result, in the present paper we found that the shift current at the band edge
diverges, in inverse proportion to the band gap in the limit. In the gapless
limit, contributions from band edges solely determine the shift current, and
when there are one conduction and one valence bands near the band edge, the
coefficient of the divergence can be written using the first- and second-order
derivatives of the two-band Hamiltonian with respect to the Bloch wave vector.
These results indicate that the diverging behavior is a universal property of
two-dimensional systems.

This paper is structured as follows. In Sec. II, we first develop a general
theory of the shift current response at the band edge of two-dimensional
systems in the gapless limit. We mainly study the case when the systems become
Weyl semimetals in this limit. We also study the cases with higher-order
energy dispersions at the end of the section. Section III deals with the two-
dimensional tight-binding model of electrons on a hexagonal lattice and
numerically confirm the diverging behavior of the shift current. In Sec. IV,
we examine our results based on the relationship between the electric
polarization and the shift current. Finally, Sec. V is devoted to the
conclusion.

##  II general theory of shift current divergences in two-dimensional gapless
limits

We first review the theory of second-order optical responses. After
introducing the shift current, we investigate the shift current response in a
two-band system in the limit where the energy gap closes and the bands form a
Dirac cone. Then, we comment on general multi-band systems with energy
dispersions of higher orders in the wavevector. Throughout this paper, we
focus on spinless and non-interacting two-dimensional systems at zero
temperature for simplicity.

###  II.1 Review of the shift current

The dc current density induced by the second-order optical response can be
written as

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

where $\omega$ is the frequency of the incident light, $j_{\mathrm{dc}}^{a}$
is the dc-current in the $a$ direction, $\sigma_{\mathrm{dc}}^{abc}$ is the
conductivity tensor, $\mathbf{E}(\omega)$ is the electric field of the
incident light, and $a,b,c$ are the Cartesian coordinates. In the clean limit,
the third-rank tensor $\sigma_{\mathrm{dc}}^{abc}$ can be divided into the
injection $\sigma_{\mathrm{inj}}^{abc}$ and the shift currents
$\sigma_{\mathrm{shift}}^{abc}$ as

$\displaystyle\sigma_{\mathrm{dc}}^{abc}=\sigma_{\mathrm{inj}}^{abc}+\sigma_{\mathrm{shift}}^{abc}.$
(2)

The injection current is caused by the change of velocities of electrons upon
optical excitations and its response tensor $\sigma_{\mathrm{inj}}^{abc}$ is
linearly proportional to the relaxation time. The shift current is caused by
the change of the positions of electrons by the excitation and its response
tensor $\sigma_{\mathrm{shift}}^{abc}$ is independent of the relaxation time.
In this paper, we focus on non-magnetic systems with time-reversal symmetry.
In this case, depending on the polarization of the irradiated light, we can
distinguish the injection and shift currents. When we irradiate a circularly
polarized light to non-magnetic systems, the second-order response is limited
to the injection current, while for the linearly polarized light, we only get
the shift current. Therefore, we here focus on a shift current response to
linearly polarized light [9, 31, 32, 33]. Then, for the linearly polarized
light, the tensor $\sigma_{\mathrm{shift}}^{abb}$ can be explicitly written as
[5, 16, 1, 3, 34]

$\displaystyle\sigma_{\mathrm{shift}}^{abb}$ $\displaystyle=-\frac{2\pi
e^{3}}{\hbar^{2}}\int_{\mathbf{k}}\sum_{n,m}f_{nm}R^{ab}_{mn}\absolutevalue{r_{nm}^{b}}^{2}\delta(\omega_{mn}-\omega),$
(3)

where $-e\ (e>0)$ is the charge of an electron,
$\int_{\mathbf{k}}\coloneqq\int_{\mathrm{BZ}}\frac{d^{2}k}{(2\pi)^{2}}$ is the
integration over the two-dimensional Brillouin zone, $f_{n}\coloneqq
f_{F}(\varepsilon_{n})$ is the Fermi distribution function of the $n$-th band
with energy $\varepsilon_{n}=\hbar\omega_{n}$, $f_{nm}\coloneqq f_{n}-f_{m}$
and $\hbar\omega_{nm}=\hbar\omega_{n}-\hbar\omega_{m}$ are the differences of
the Fermi distribution functions and energies of the $n$-th and $m$-th bands,
respectively, and $r_{nm}^{b}=\bra{n}i\partial_{k_{b}}\ket{m}$ is the
$(n,\,m)$-element of the position operator. The quantity
$R^{ab}_{mn}=-\partial_{k_{a}}\mathrm{arg}(r_{mn}^{b})+r_{mm}^{a}-r_{nn}^{a}$
is called the shift vector. The above formula of the quantity $R^{ab}_{mn}$
involves derivatives of energy eigenvectors, and thus is not convenient for
our study. Thus we use the formula for the shift current without the use of
derivatives of energy eigenvectors in Ref. [35]. An equality for the off-
diagonal elements of the position operator
$r_{nm}^{b}=iv^{b}_{nm}/\omega_{mn}$ holds, where $v_{nm}^{b}=\bra{n}\partial
H/\partial k_{b}\ket{m}/\hbar$ is the element of the velocity operator. Then,
using these relations, $I_{mn}^{abb}\coloneqq
R^{ab}_{mn}\absolutevalue{r_{nm}^{b}}^{2}$ can be expressed in terms of the
velocity operators as

$\displaystyle I_{mn}^{abb}$
$\displaystyle=\frac{1}{\varepsilon_{nm}^{2}}\imaginary\left[\frac{v_{nm}^{a}\quantity(v_{mm}^{b}-v_{nn}^{b})+\quantity(a\leftrightarrow
b)}{\varepsilon_{nm}}v_{mn}^{b}\right.$
$\displaystyle\quad\left.+w_{nm}^{ab}v_{mn}^{b}-\sum_{p\neq
n,m}\quantity(\frac{v_{np}^{b}v_{pm}^{a}}{\varepsilon_{pm}}-\frac{v_{np}^{a}v_{pm}^{b}}{\varepsilon_{np}})v_{mn}^{b}\right],$
(4)

where $\varepsilon_{mn}=\hbar\omega_{mn}$ and
$w_{nm}^{ab}=\bra{n}\partialderivative{H}{k_{a}}{k_{b}}\ket{m}$. The first two
terms are direct excitation processes and the third term represents three-band
processes, which vanishes for two-band models. Below, we consider a zero-
temperature limit where $f_{nm}=-1$ if the $n$-th band is occupied and the
$m$-th band is empty, $f_{nm}=0$ if both the $n$-th and $m$-th bands are
occupied or empty.

###  II.2 Shift current responses in the gapless limit

For simplicity, we first confine our discussions to two-band gapped systems
and calculate $I_{12}^{abb}$, where the band indices $m=1$ and $2$ denote the
occupied and the unoccupied bands, respectively. Here we consider a system
described by the Bloch Hamiltonian where $g(\mathbf{k})$ is a complex function
of $\mathbf{k}=\quantity(k_{x},k_{y})$ and $M$ is a real parameter. Then, the
energies of this system are

$\displaystyle\varepsilon_{1}$ $\displaystyle=-\varepsilon,$ (5)
$\displaystyle\varepsilon_{2}$ $\displaystyle=+\varepsilon,$ (6)
$\displaystyle\varepsilon$ $\displaystyle=\sqrt{\absolutevalue{g}^{2}+M^{2}}.$
(7)

We assume that the complex function $g(\mathbf{k})$ becomes zero somewhere in
the Brillouin zone. Then, the gap is given by $2\absolutevalue{M}$, and this
system becomes gapless if and only if $g(\mathbf{k})=0$ and $M=0$. By a
straightforward calculation of Eq. (4), $I_{12}^{abb}$ in the integrand can be
written as

$\displaystyle I_{12}^{abb}$
$\displaystyle=-\frac{M}{4\varepsilon^{3}}\imaginary\quantity[\partialderivative{g^{*}}{k_{b}}\partialderivative{g}{k_{b}}{k_{a}}-\frac{1}{2\varepsilon^{2}}\partialderivative{\absolutevalue{g}^{2}}{k_{b}}\partialderivative{g^{*}}{k_{b}}\partialderivative{g}{k_{a}}].$
(8)

We can see that for all $\mathbf{k}$ such that $g(\mathbf{k})\neq 0$,
$I_{12}^{abb}(\mathbf{k})$ approaches $0$ in the gapless limit $M\to 0$ due to
the overall factor $M$ in Eq. (8). On the other hand, at the band edge
$\mathbf{k}_{0}$ such that $g(\mathbf{k}=\mathbf{k}_{0})=0$, $\varepsilon$
becomes $\absolutevalue{M}$, $I_{12}^{abb}(\mathbf{k})$ has an overall factor
$M^{-2}$ and the second term in the square bracket in Eq. (8) becomes zero. As
a result, the shift current in the $M\to 0$ limit is solely determined by the
value of the first term in the square bracket in Eq. (8) evaluated at the band
edge $\mathbf{k}_{0}$.

To calculate $\sigma_{\mathrm{shift}}^{abb}$, we need to consider the joint
density of states (JDOS) part
$\int_{\bf{k}}\delta\quantity(\omega_{mn}-\omega)$. If the leading term of
$g(\mathbf{k})$ around the band edge is proportional to
$\mathbf{k}-\mathbf{k}_{0}$, the integration of the delta function in Eq. (3)
reduces to a term proportional to $\absolutevalue{M}$ as we see below. In
order to carry out the integration with respect to the reciprocal vector, we
need to change the variable of the delta function from $\omega$ to
$\mathbf{k}$. Here, we follow the discussion in Ref. [35] and integrate around
the band edge using a polar coordinate. Around the band edge $\mathbf{k}_{0}$
such that $g(\mathbf{k}_{0})=0$, we can expand the function as
$g(\mathbf{k})\approx\nabla_{\mathbf{k}}g\cdot(\mathbf{k}-\mathbf{k}_{0})$,
where the derivative is with respect to the reciprocal vectors. Then, we can
approximate the energy $\varepsilon$ up to parabolic terms using
$\delta\mathbf{k}\coloneqq\mathbf{k}-\mathbf{k}_{0}$ as

$\displaystyle\varepsilon\approx\absolutevalue{M}+\frac{1}{2\absolutevalue{M}}\absolutevalue{\partial_{k_{x}}g(\mathbf{k}_{0})\cdot\delta
k_{x}+\partial_{k_{y}}g(\mathbf{k}_{0})\cdot\delta k_{y}}^{2}.$ (9)

We introduce a polar coordinate around $\mathbf{k}_{0}$ by variables $k$ and
$\theta$ (see Appendix A and set $n=1$). Using these new variables, the JDOS
at can be calculated as

$\displaystyle\int_{\mathbf{k}}h(\mathbf{k})\delta(2\varepsilon-\omega)$
$\displaystyle\approx\int\frac{\mathrm{d}k\mathrm{d}\theta}{(2\pi)^{2}}\cdot\frac{h(\mathbf{k})\cdot
2\absolutevalue{M}k}{\absolutevalue{\imaginary[\partial_{k_{x}}g(\mathbf{k}_{0})(\partial_{k_{y}}g(\mathbf{k}_{0}))^{*}]}}$
$\displaystyle\qquad\times\frac{\delta\quantity(k-\sqrt{(\omega-2\absolutevalue{M})\absolutevalue{M}})}{4k}$
$\displaystyle=\frac{1}{2(2\pi)^{2}}\frac{\absolutevalue{M}}{\absolutevalue{\imaginary[\partial_{k_{x}}g(\mathbf{k}_{0})(\partial_{k_{y}}g(\mathbf{k}_{0}))^{*}]}}$
$\displaystyle\qquad\times\int\mathrm{d}\theta
h(k,\theta)\evaluated{}_{k=\sqrt{(\omega-2\absolutevalue{M})\absolutevalue{M}}},$
(10)

where $h(\mathbf{k})$ is an arbitrary function. As shown above,
$h(\mathbf{k})$ for the calculation of the shift conductivity
$\sigma_{\mathrm{shift}}^{abb}$ is $I_{12}^{abb}$, which is non-zero only at
$k=0$. Then, combining Eqs. (8) and (10), in the gapless limit $M\to 0$, we
get

$\displaystyle\sigma_{\mathrm{shift}}^{abb}(\omega)\evaluated{}_{\omega=2M}\sim\frac{e^{3}}{4\hbar^{2}}\sum_{g(\mathbf{k})=0}\frac{\imaginary\quantity[\partialderivative{g^{*}}{k_{b}}\partialderivative{g}{k_{b}}{k_{a}}]}{\absolutevalue{\imaginary[\partial_{k_{x}}g(\partial_{k_{y}}g)^{*}]}}\cdot\frac{1}{2M}.$
(11)

On the other hand, contributions from points other than the band edge
$(\omega>2M)$, the shift conductivity tensor becomes zero in the $M\to 0$
limit. This is one of our main results. We can see that the shift current
response for the light corresponding to the band edge is inversely
proportional to the gap size and is divergent in the long-wavelength limit, or
equivalently, the gapless limit. The coefficient of $1/M$ includes the first-
and second-order derivatives of $g(\mathbf{k})$. These are generally non-zero
in realistic systems.

This divergence is different from the low-frequency divergence proportional to
$1/\omega^{4-d}$, where $d$ is the dimension of the Weyl semimetals discussed
in Ref. [9]. The low-frequency divergence discussed in Ref. [9] is originated
from the linear dispersion and it is present only if the Dirac cone is tilted
and the system is in a state called type-II Weyl (Dirac) semimetal phase. The
$1/\omega$ divergence discussed in this paper is present without the tilting
of the Dirac cone.

Figure 1: The schematic illustration of the multi-band system. The unoccupied
$m$-th band and occupied $n$-th band form a Weyl point in the limit $M\to 0$.
The $p$-th band $(p\neq m,n)$ is not degenerate with these two bands.

Next, we consider general multi-band systems. For systems with more than two
bands, the third term in Eq. (4) is nonzero. Here we analyze the $M$
dependency of this term. We consider three bands indexed as $m$, $n$, and $p$,
where the $m$-th band is unoccupied and the $n$-th band is occupied and they
form a Weyl point in the Weyl semimetal limit. The $p$-th band $(p\neq
m,\,p\neq n)$ is not degenerate with these two bands (Fig. 1). First,
$\varepsilon_{pm}$ approaches a finite value even in the Weyl semimetal limit
because only the two bands $m$ and $n$ are degenerate in this limit. Likewise,
terms of the form of $v_{np}^{a}$ are expected to approach a finite value in
the limit $M\to 0$. Next we discuss the terms that involve derivatives of
eigenstates, such as
$v_{np}^{a}=\bra{n}\ket{\partial_{k_{a}}p}/\varepsilon_{pn}$, $p\neq n$. As
discussed above, $\varepsilon_{pn}$ approaches a finite value as $M\to 0$. The
amplitude of $\bra{n}\ket{\partial_{k_{a}}p}$ is proportional to the
transition probability between two states $\ket{p}$ and $\ket{n}$. As the
energies $\varepsilon_{p}$ and $\varepsilon_{n}$ are far apart, we expect the
mixing of these two states does not behave singularly in the $M\to 0$ limit.
Therefore, the terms containing the sum over $p$ in Eq. (4) are of the order
of $M^{0}$. Thus even in the presence of these terms, the leading divergence
in the two-dimensional Weyl semimetal limit is still $1/M$.

To end this section, we consider a higher-order band degeneracy such as
quadratic band degeneracy as a comparison with a Weyl point. Suppose that
$g(\mathbf{k})=(c_{1}\delta k_{x}+c_{2}\delta
k_{y})^{n}+\order{\absolutevalue{\delta\mathbf{k}}^{n+1}}\quad(c_{1},c_{2}\in\mathbb{C},\,n\in\mathbb{N})$
around the gap minimum. Even in this case, contributions from
$\mathbf{k}$-points where $g(\mathbf{k})\neq 0$ goes to zero in the $M\to 0$
limit as the integrand in Eq. (4) has an overall factor $M$. For
$\mathbf{k}$-points around $g(\mathbf{k})=0$, we can still carry out the
integration of the JDOS by using a polar coordinate. Then, by similar
arguments with the case of $n=1$ (for the details of the calculation, see
Appendix A), we get

$\displaystyle\sigma_{\mathrm{shift}}(\omega)\sim\frac{\mathrm{sgn}(M)}{\absolutevalue{M}^{2-\frac{1}{n}}\quantity(\omega-2\absolutevalue{M})^{1-\frac{1}{n}}}\quad(\omega\gtrapprox
2\absolutevalue{M}),$ (12)

for $\omega$ slightly larger than $2\absolutevalue{M}$. When $n>1$, the shift
current response is divergent at the frequency of the band edge even with a
nonzero value of $M$. The divergent term
$\absolutevalue{\omega-2\absolutevalue{M}}^{\frac{1}{n}-1}$ comes from the
JDOS and we can see that the $1/M$ divergence discussed above is unique to the
Weyl point with linear energy dispersion ($n=1$), where the JDOS is not
divergent at the band edge.

##  III tight-binding model calculation

In this section, we numerically calculate
$\sigma_{\mathrm{shift}}^{abb}(\omega)$ for a two-dimensional tight-binding
model of electrons on an anisotropic hexagonal lattice (Fig. 2) and confirm
the divergent behavior of the shift current response in the gapless limit.

Figure 2: The tight-binding model on a two-dimensional hexagonal lattice.
Anisotropic nearest-neighbor hoppings are parametrized by $t_{1}$ for the
solid lines and $t_{2}$ for the dashed lines. The on-site potentials at the
sites in the sublattices A and B are $+M$ and $-M$, respectively. The lattice
constant is given by $a$. This system has an inversion symmetry if $M=0$ and
has a $C_{3}$-rotational symmetry if $t_{1}=t_{2}$.

The tight-binding Hamiltonian of the system is given by

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=\sum_{\langle
i,j\rangle}t_{ij}\hat{c}_{i}^{\dagger}\hat{c_{j}}+\sum_{i}M_{i}\hat{c}_{i}^{\dagger}\hat{c}_{i},$
(13)

where the first sum is over all the pairs of neighboring sites. The parameter
$t_{ij}\ (>0)$ describes the nearest-neighbor hoppings between neighboring
sites $i$ and $j$. For the hoppings indicated by the solid lines in Fig. 2, we
set $t_{ij}=t_{1}$ and for those by the dashed lines, we set $t_{ij}=t_{2}$.
When $t_{1}\neq t_{2}$, this system lacks $C_{3}$-rotational symmetry. The
real parameter $M_{i}$ is the on-site potential and takes the values $+M$ and
$-M$ on the sublattices A and B, respectively. When $M\neq 0$, this term
breaks inversion symmetry. Together with the $C_{3}$-rotational symmetry
breaking, this lack of inversion symmetry allows non-zero shift current
generation. After the Fourier transformation, the Bloch Hamiltonian of this
system is given in the form of Eq. (LABEL:eq:2level_Hamiltonian) with

$\displaystyle
g(\mathbf{k})=2t_{1}\cos\quantity(\frac{\sqrt{3}a}{2}k_{x})e^{i\frac{a}{2}k_{y}}+t_{2}e^{-iak_{y}}.$
(14)

This $g(\mathbf{k})$ becomes zero at
$(k_{x},k_{y})=\quantity(2\cos^{-1}\quantity(t_{2}/(2t_{1}))/(\sqrt{3}a),2\pi/(3a))$
and $\quantity(2\cos^{-1}\quantity(-t_{2}/(2t_{1}))/(\sqrt{3}a),4\pi/(3a))$ if
$\absolutevalue{t_{2}/t_{1}}\leq 2$, where we set $0\leq\cos^{-1}x\leq\pi$. At
these points, the bands have an energy gap of $2\absolutevalue{M}$. When
$M=0$, two bands are degenerate at these $\mathbf{k}$-points with linear
dispersion around them and this system is a two-dimensional Weyl semimetal. By
substituting this function $g(\mathbf{k})$ into the general formula in Eq.
(8), we can calculate the shift current conductivity. For the numerical
calculation, we approximated the delta function with a Lorenzian as
$\delta(\omega)\approx 1/\pi\cdot\epsilon/(\epsilon^{2}+\omega^{2})$, where
$\epsilon$ is a positive infinitesimal.

Figure 3: The shift conductivity $\sigma_{\mathrm{shift}}^{yxx}(\omega)$ for
the tight-binding model (13) for $M/t_{1}=0.05,0.15,0.25,0.385,0.5,1$.
Vertical lines indicate the band edges $\omega=2M$. We set $t_{2}/t_{1}=0.5$.
We set the Fermi energy $E_{F}$ to be $E_{F}=0$, which means that the bands
are half-filled.

The calculated $\sigma^{yxx}_{\mathrm{shift}}$ is plotted as a function of
$\omega$ for several values of $M$ in Fig. 3. Solid curves with different
colors indicate the calculated $\sigma^{yxx}_{\mathrm{shift}}$ of this system.
The dashed curve indicates the expected $1/M$ divergence in Eq. (11). The
sizes of the band gap $2\absolutevalue{M}$ are indicated by vertical lines for
each $M$ in the same color as $\sigma_{\mathrm{shift}}^{yxx}$. We can see that
the shift current response is not divergent at band edges
$\omega=2\absolutevalue{M}$ shown by the vertical lines for each value of $M$.
On the other hand, when we make $M$ small, the shift current responses of the
system at band edges diverges as $1/M$ and the curve shows a good agreement
with the expected divergence from Eq. (11).

Similar results can be obtained for $\sigma^{yyy}_{\mathrm{shift}}$. Note that
since the system has a mirror symmetry with respect to the $y$-axis,
$\sigma^{xaa}_{\mathrm{shift}}=0$ for $a=x,y$, and no currents are generated
in the $x$-direction.

##  IV Relationship with electric polarization

The shift current is caused by the shift of positions of electrons along with
the optical excitation and hence, related to the electric polarization of the
material [7]. In this section, we show that a simple relationship between
electric polarizations and the shift current is lost in the two-dimensional
Weyl semimetal limit in general.

Figure 4: A Weyl point in the reciprocal space. Along the blue line passing
the Weyl point in the direction of a unit vector $\mathbf{e}_{a}$, that has an
azimuthal angle $\varphi$ and a polar angle $\theta$, the off-diagonal element
of the position operator $\mathbf{r}_{12}=i\bra{1}\nabla_{\mathbf{k}}\ket{2}$
has no radial components around the Weyl point.

Consider an effective Hamiltonian of a three-dimensional Weyl point

$\displaystyle H_{\mathrm{eff}}=v\mathbf{k}\cdot\boldsymbol{\sigma},$ (15)

where $v$ is a positive constant and
$\boldsymbol{\sigma}=\quantity(\sigma_{x},\,\sigma_{y},\,\sigma_{z})$ are
Pauli matrices. Then, eigenstates of this Hamiltonian can be written as

$\displaystyle\ket{n}$
$\displaystyle=\frac{1}{\sqrt{2\absolutevalue{\mathbf{k}}}}\matrixquantity(\eta\sqrt{\absolutevalue{\mathbf{k}}+\eta
k_{z}}\\\ e^{i\beta}\sqrt{\absolutevalue{\mathbf{k}}-\eta k_{z}}),$ (16)

where $\eta=(-1)^{n}\ (n=1,2)$ and $\beta=\tan^{-1}\quantity(k_{y}/k_{x})$ is
an azimuthal angle of $\mathbf{k}$ in a spherical coordinate. The
corresponding energies are $\varepsilon_{1}=-v\absolutevalue{\mathbf{k}}$ and
$\varepsilon_{2}=+v\absolutevalue{\mathbf{k}}$.

We first check the behavior of the off-diagonal element of the velocity
operator $\mathbf{v}_{12}$ in a direction parallel to an arbitrary unit vector
$\mathbf{e}_{a}$. Then, under the gauge choice of Eq. (16),
$v_{12}^{a}\coloneqq\mathbf{v}_{12}\cdot\mathbf{e}_{a}$ is

$\displaystyle v_{12}^{a}$
$\displaystyle=\frac{v}{\hbar}\left\\{\partialderivative{k_{x}}{k_{a}}\quantity(\cos\alpha\cos\beta-i\sin\beta)\right.$
$\displaystyle\quad+\partialderivative{k_{y}}{k_{a}}\quantity(\cos\alpha\sin\beta+i\cos\beta)$
$\displaystyle\quad\left.-\partialderivative{k_{z}}{k_{a}}\sin\alpha\right\\},$
(17)

where $\alpha=\tan^{-1}\quantity(\sqrt{k_{x}^{2}+k_{y}^{2}}/k_{z})$ is a polar
angle of $\mathbf{k}$. As shown in Fig. 4, let $\theta$ and $\varphi$ denote
the polar and azimuthal angles of the vector $\mathbf{e}_{a}$ in the spherical
coordinate. Then, using $\partial k_{x}/\partial
k_{a}=\cos\varphi\sin\theta,\,\partial k_{y}/\partial
k_{a}=\sin\varphi\sin\theta,$ and $\partial k_{z}/\partial k_{a}=\cos\theta$,
the quantity $\absolutevalue{v_{12}^{a}}^{2}$ can be written using $\varphi$,
$\theta$, and $\beta=\tan^{-1}\quantity(k_{y}/k_{x})$ as

$\displaystyle\absolutevalue{v_{12}^{a}}^{2}$
$\displaystyle=\frac{v^{2}}{\hbar^{2}}\left\\{\quantity(\cos\alpha\sin\theta\cos\quantity(\beta-\varphi)-\sin\alpha\cos\theta)^{2}\right.$
$\displaystyle\qquad\left.+\sin^{2}\theta\sin^{2}\quantity(\beta-\varphi)\right\\}.$
(18)

This quantity becomes zero if and only if
$\quantity(\varphi,\,\theta)=\quantity(\beta,\,\alpha)$ or
$\quantity(\beta+\pi,\,\pi-\alpha)$, meaning that
$\absolutevalue{v_{12}^{a}}=0$ on a line in $\mathbf{k}$-space crossing the
Weyl point and parallel to the direction of $\mathbf{e}_{a}$ (blue line in
Fig. 4). This result is independent of the choice of gauges as
$\absolutevalue{v_{12}^{a}}$ is a gauge-independent quantity. Since
$r_{12}^{a}=iv_{12}^{a}/(2\varepsilon)$, $r_{12}^{a}$ also becomes zero on the
line crossing the Weyl point and parallel to the $\mathbf{e}_{a}$.

Now, we discuss whether the shift current and the electric polarization behave
similarly. In the Ref. [7], it is shown that if $r_{12}^{a}(\mathbf{k})$ can
be approximated to be constant as a function of $\mathbf{k}$; then
$\absolutevalue{r_{12}^{a}}^{2}$ in Eq. (3) can be taken out of the integral,
leading to the quantitative relationship between the shift current
conductivities and the difference of electric polarizations as

$\displaystyle\int\mathrm{d}\omega\sigma_{\mathrm{shift}}^{aaa}\approx-\frac{\pi
e^{2}\absolutevalue{r^{a}_{12}}^{2}}{\hbar^{2}}\quantity(P^{a}_{2}-P^{a}_{1}).$
(19)

This relationship is consistent with the understanding of the shift current as
a shift of the positions of electrons via optical excitations. In our setting,
however, the element $\absolutevalue{r_{12}^{a}}$ goes to zero around the Weyl
point in the two-dimensional Weyl semimetal limit. Hence, the assumption of
the argument above is not valid in this case and the quantitative relationship
between the shift conductivity and the electric polarization is lost in the
two-dimensional Weyl semimetal limit.

##  V conclusion

In this paper, we investigated the shift current response of two-dimensional
systems when the system approaches a gapless limit and found a divergence of
the shift conductivity when the system in the gapless limit is a Weyl
semimetal. This work is motivated by the authors’ previous findings that by a
change of the system across a two-dimensional Weyl semimetal phase, the
electric polarization has a discrete jump. By noting the quantitative
relationships between the shift conductivity and electric polarizations in the
previous works, we expected a jump of the shift conductivity across a two-
dimensional Weyl semimetal phase.

The calculation of the shift conductivity around the band edge can be divided
into the geometric part and the joint density of states. For the direct
excitation of electrons from one band to another band, the geometric part is
proportional to $M^{-2}$, where $M$ parameterizes the size of the band gap. As
the joint density of state part is proportional to $M$ in the Weyl semimetal
limit, the contribution from this direct excitation to the shift conductivity
is proportional to $M^{-1}$ and diverges as we take the Weyl semimetal limit.
We found that the coefficient of this divergence is given by the first- and
second-order derivatives of the Hamiltonian with respect to the Bloch
wavevector and confirmed its validity by numerical calculations. For multi-
band systems, indirect excitation processes via transitions to another band
are involved, making the divergence coefficient more complex but still we can
show that the contribution to the shift conductivity also diverges as
$M^{-1}$. Taken together, the shift conductivity in the two-dimensional Weyl
semimetal limit is divergent, in inverse proportion to the size of the gap.

This divergence crucially depends on the dimensionality and the dispersion of
the energy degeneracy and hence, is unique to the two-dimensional Weyl
semimetal limit. The unexpected finding of divergence rather than a jump in
shift conductivity was due to the fact that the quantitative relationship
between electric polarization and shift conductivity is lost in the Weyl
semimetallic limit, and the intuitive picture of the shift current as a change
in the position of the center of charge can no longer be applied.

###### Acknowledgements.

This work is partly supported by Japan Society for the Promotion of Science
(JSPS) KAKENHI Grant Numbers JP24KJ1109, JP22K18687, JP22H00108, and
JP24H02231. HY is also supported by JST SPRING, Grant Number JPMJSP2106.

##  Appendix A Shift conductivity for energy degeneracies with higher order
dispersions

Here, we calculate the shift conductivity of a two-dimensional system in a
gapless limit with higher-order energy dispersion in $\mathbf{k}$ discussed at
the end of Sec. II. If an effective Hamiltonian is given in the same form as
Eq. (LABEL:eq:2level_Hamiltonian), we only need to consider contributions from
band edges. Suppose that the off-diagonal component of an effective
Hamiltonian $g(\mathbf{k})$ is written as

$\displaystyle g(\mathbf{k})=\quantity(c_{1}\delta k_{x}+c_{2}\delta
k_{y})^{n}+\mathcal{O}(\absolutevalue{\delta\mathbf{k}}^{n+1}),$ (20)

where $c_{1}$ and $c_{2}$ are complex constants and $\delta k_{x}$ and $\delta
k_{y}$ are the small displacements from the band edge. Then, the energy is
approximated as

where $\lambda_{1,2}$ are positive numbers satisfying
$\lambda_{1}+\lambda_{2}=\absolutevalue{c_{1}}^{2}+\absolutevalue{c_{2}}^{2}$
and $\lambda_{1}\lambda_{2}=\quantity(\imaginary(c_{1}c_{2}^{*}))^{2}$. The
original coordinate $(\delta k_{x},\,\delta k_{y})$ is transformed by an
orthogonal matrix $O$ as $(\delta k_{x}^{\prime},\,\delta
k_{y}^{\prime})^{\mathrm{T}}=O(\delta k_{x},\,\delta k_{y})^{\mathrm{T}}$.
Then, we introduce a polar coordinate $(k,\,\theta)$ around the band edge as
$\sqrt{\lambda_{1}}\delta k_{x}^{\prime}=k\cos\theta$ and
$\sqrt{\lambda_{2}}\delta k_{y}^{\prime}=k\sin\theta$. Using these variables,
we get

$\displaystyle\varepsilon\approx\absolutevalue{M}+\frac{1}{2\absolutevalue{M}}k^{2n}.$
(21)

In general, we can calculate the integration of our interest around the band
edge as

$\displaystyle\phantom{=}\int\frac{\mathrm{d}\delta k_{x}\mathrm{d}\delta
k_{y}}{(2\pi)^{2}}h(\mathbf{k})\delta\quantity(2\varepsilon-\omega)$
$\displaystyle=\int\frac{\mathrm{d}k\mathrm{d}\theta}{(2\pi)^{2}}h(\mathbf{k})\absolutevalue{\frac{\partial\quantity(\delta
k_{x},\delta
k_{y})}{\partial\quantity(k,\theta)}}\delta\quantity(\frac{1}{\absolutevalue{M}}k^{2n}-\quantity(\omega-2\absolutevalue{M}))$
$\displaystyle=\int\frac{\mathrm{d}k\mathrm{d}\theta}{(2\pi)^{2}}\frac{h(\mathbf{k})\cdot
k}{\sqrt{\lambda_{1}\lambda_{2}}}\frac{\delta\quantity(k-\quantity[(\omega-2\absolutevalue{M})\absolutevalue{M}]^{\frac{1}{2n}})}{\absolutevalue{\frac{2nk^{2n-1}}{\absolutevalue{M}}}}$
$\displaystyle=\frac{1}{2n(2\pi)^{2}\absolutevalue{\imaginary(c_{1}c_{2}^{*})}}\cdot\frac{\absolutevalue{M}^{\frac{1}{n}}}{\quantity(\omega-2\absolutevalue{M})^{1-\frac{1}{n}}}$
$\displaystyle\hskip 60.0pt\times\int\mathrm{d}\theta
h\quantity(k,\,\theta)\evaluated{}_{k=\quantity[(\omega-2\absolutevalue{M})\absolutevalue{M}]^{\frac{1}{2n}}}.$
(22)

For the calculation of $\sigma_{\mathrm{shift}}^{abb}$, $h(\mathbf{k})$ is
non-zero only at band edge in the limit $M\to 0$ and thus, the $\theta$
integral yields $2\pi$. Combining this result with a $1/M^{2}$ dependence of
other integrands of $\sigma_{\mathrm{shift}}^{abb}$, the $M$ dependence of
$\sigma_{\mathrm{shift}}^{abb}$ is

$\displaystyle\sigma_{\mathrm{shift}}^{abb}(\omega)\sim\frac{\mathrm{sgn}(M)}{\absolutevalue{M}^{2-\frac{1}{n}}\quantity(\omega-2\absolutevalue{M})^{1-\frac{1}{n}}}.$
(23)

## References

  * Kraut and von Baltz [1979] W. Kraut and R. von Baltz, Anomalous bulk photovoltaic effect in ferroelectrics: A quadratic response theory, Phys. Rev. B 19, 1548 (1979).
  * Belinicher and Sturman [1980] V. I. Belinicher and B. I. Sturman, The photogalvanic effect in media lacking a center of symmetry, Sov. Phys. Usp. 23, 199 (1980).
  * von Baltz and Kraut [1981] R. von Baltz and W. Kraut, Theory of the bulk photovoltaic effect in pure crystals, Phys. Rev. B 23, 5590 (1981).
  * Aversa and Sipe [1995] C. Aversa and J. E. Sipe, Nonlinear optical susceptibilities of semiconductors: Results with a length-gauge analysis, Phys. Rev. B 52, 14636 (1995).
  * Sipe and Shkrebtii [2000] J. E. Sipe and A. I. Shkrebtii, Second-order optical response in semiconductors, Phys. Rev. B 61, 5337 (2000).
  * Fridkin [2001] V. M. Fridkin, Bulk photovoltaic effect in noncentrosymmetric crystals, Crystallogr. Rep. 46, 654 (2001).
  * 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).
  * Pusch  _et al._ [2023] A. Pusch, U. Römer, D. Culcer, and N. J. Ekins-Daukes, Energy Conversion Efficiency of the Bulk Photovoltaic Effect, PRX Energy 2, 013006 (2023).
  * 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).
  * Ahn  _et al._ [2022] J. Ahn, G.-Y. Guo, N. Nagaosa, and A. Vishwanath, Riemannian geometry of resonant optical responses, Nat. Phys. 18, 290 (2022).
  * 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).
  * Young  _et al._ [2012] S. M. Young, F. Zheng, and A. M. Rappe, First-Principles Calculation of the Bulk Photovoltaic Effect in Bismuth Ferrite, Phys. Rev. Lett. 109, 236601 (2012).
  * Belinicher  _et al._ [1982] V. I. Belinicher, E. L. Ivchenko, and B. I. Sturman, Kinetic theory of the displacement photovoltaic effect in piezoelectrics, Sov. Phys. JETP 56, 359 (1982), [Zh. Eksp. Teor. Fiz. 83, 649].
  * Hsu  _et al._ [2023] H.-C. Hsu, J.-S. You, J. Ahn, and G.-Y. Guo, Nonlinear photoconductivities and quantum geometry of chiral multifold fermions, Phys. Rev. B 107, 155434 (2023).
  * Morimoto and Nagaosa [2016a] T. Morimoto and N. Nagaosa, Topological aspects of nonlinear excitonic processes in noncentrosymmetric crystals, Phys. Rev. B 94, 035117 (2016a).
  * Morimoto and Nagaosa [2016b] T. Morimoto and N. Nagaosa, Topological nature of nonlinear optical effects in solids, Sci. Adv. 2, e1501524 (2016b).
  * Resta [2011] R. Resta, The insulating state of matter: a geometrical theory, Eur. Phys. J. B 79, 121 (2011).
  * Ma  _et al._ [2021] Q. Ma, A. G. Grushin, and K. S. Burch, Topology and geometry under the nonlinear electromagnetic spotlight, Nat. Mater. 20, 1601 (2021).
  * Morimoto  _et al._ [2023] T. Morimoto, S. Kitamura, and N. Nagaosa, Geometric Aspects of Nonlinear and Nonequilibrium Phenomena, J. Phys. Soc. Jpn. 92, 072001 (2023).
  * de Juan  _et al._ [2017] F. de Juan, A. G. Grushin, T. Morimoto, and J. E. Moore, Quantized circular photogalvanic effect in Weyl semimetals, Nat. Commun. 8, 15995 (2017).
  * Resta [1992] R. Resta, Theory of the electric polarization in crystals, Ferroelectrics 136, 51 (1992).
  * King-Smith and Vanderbilt [1993] R. D. King-Smith and D. Vanderbilt, Theory of polarization of crystalline solids, Phys. Rev. B 47, 1651 (1993).
  * Vanderbilt and King-Smith [1993] D. Vanderbilt and R. D. King-Smith, Electric polarization as a bulk quantity and its relation to surface charge, Phys. Rev. B 48, 4442 (1993).
  * 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).
  * Osterhoudt  _et al._ [2019] G. B. Osterhoudt, L. K. Diebel, M. J. Gray, X. Yang, J. Stanco, X. Huang, B. Shen, N. Ni, P. J. W. Moll, Y. Ran, and K. S. Burch, Colossal mid-infrared bulk photovoltaic effect in a type-I Weyl semimetal, Nat. Mater. 18, 471 (2019).
  * Kim  _et al._ [2017] K. W. Kim, T. Morimoto, and N. Nagaosa, Shift charge and spin photocurrents in Dirac surface states of topological insulator, Phys. Rev. B 95, 035134 (2017).
  * Tan and Rappe [2016] L. Z. Tan and A. M. Rappe, Enhancement of the Bulk Photovoltaic Effect in Topological Insulators, Phys. Rev. Lett. 116, 237402 (2016).
  * Jankowski and Slager [2024] W. J. Jankowski and R.-J. Slager, Quantized shift response in multi-gap topological phases (2024), arXiv:2402.13245.
  * Yoshida  _et al._ [2023a] H. Yoshida, T. Zhang, and S. Murakami, Polarization jumps by breaking symmetries of two-dimensional Weyl semimetals, Phys. Rev. B 107, 035122 (2023a).
  * Yoshida  _et al._ [2023b] H. Yoshida, T. Zhang, and S. Murakami, Polarization jumps across topological phase transitions in two-dimensional systems, Phys. Rev. B 108, 075160 (2023b).
  * Watanabe and Yanase [2021] H. Watanabe and Y. Yanase, Chiral Photocurrent in Parity-Violating Magnet and Enhanced Response in Topological Antiferromagnet, Phys. Rev. X 11, 011001 (2021).
  * 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. Res. 2, 033100 (2020).
  * 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).
  * Barik and Sau [2020] T. Barik and J. D. Sau, Nonequilibrium nature of nonlinear optical response: Application to the bulk photovoltaic effect, Phys. Rev. B 101, 045201 (2020).
  * 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, Nat. Commun. 8, 14176 (2017).