# Emerging nonlinear Hall effect in Kane-Mele two-dimensional topological
insulators

Rajesh K. Malla  Avadh Saxena  Theoretical Division and Center for Nonlinear
Studies, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
Wilton J. M. Kort-Kamp  Theoretical Division, Los Alamos National Laboratory,
Los Alamos, New Mexico 87545, USA

###### Abstract

The recent observations of nonlinear Hall effect in time-reversal symmetry
protected systems and on the surface of three-dimensional topological
insulators due to an in-plane magnetic field have attracted immense
experimental and theoretical investigations in two-dimensional transition
metal dichalcogenides and Weyl semimetals. The origin of this type of second
order effect has been attributed to the emergence of a Berry curvature dipole,
which requires a low-symmetry environment. Here, we propose a mechanism for
generating such a second order nonlinear Hall effect in Kane-Mele two-
dimensional topological insulators due to spatial and time reversal symmetry
breaking in the presence of Zeeman and Rashba couplings. By actively tuning
the energy gaps with external electromagnetic fields we also demonstrate that
the nonlinear Hall effect shows remarkable signatures of topological phase
transitions existing in the considered two-dimensional systems.

Two-dimensional materials with strong intrinsic spin-orbit coupling (SOC,
$\lambda_{SO}$) have attracted intense investigations due to the possibility
to control and manipulate both electric and spin currents for application in
various optoelectronic and spintronic devices. One of the intriguing phenomena
emerging from the SOC is the quantum spin Hall effect, which describes a
nonzero transverse spin current in systems with time-reversal (TRS) and
spatial-inversion (SIS) symmetries Kane and Mele [2005a, b]. It has been
experimentally verified in various semiconducting materials, including HgTe
Bernevig et al. [2006]; König et al. [2007] and InAs/GaSb quantum wells Liu et
al. [2008]; Du et al. [2015], as well as in WTe2 transition metal
dichalcogenides (TMD) Qian et al. [2014]; Tang et al. [2017]. The inclusion of
Rashba spin-orbit coupling ($\lambda_{R}$) breaks SIS in these systems and
destroys spin angular momentum conservation while the gapless edge states and
the spin topological invariant (i.e., spin Chern number) remain unchanged
provided $\lambda_{R}<\lambda_{SO}$ Sheng et al. [2006]. Moreover, since the
robustness of the spin Chern number is fully determined by the bulk band gap,
it has been shown that the quantum spin Hall phase persists even when TRS is
broken Prodan [2009]. This further allows one to extend spin manipulation to
TRS-broken quantum spin Hall insulator (QSHI) states, e.g., via exchange
interaction within the Kane-Mele Hamiltonian Yang et al. [2011]. Nevertheless,
studies on the optical properties and electronic transport in these cases have
been largely focused on the linear response regime Stille et al. [2012];
Tabert and Nicol [2013]; Ledwith et al. [2018]; Kort-Kamp [2017]; Rodriguez-
Lopez et al. [2018]; Ezawa [2013]; Rodriguez-Lopez et al. [2017]; Tabert and
Nicol [2014]; Farias et al. [2018]; Malla and Raikh [2018a]; Wu et al. [2018,
2020].

Second order nonlinear electromagnetic phenomena, including second harmonic
generation (SHG), sum- and difference frequency generation, optical
rectification, and the Pockels effect are very sensitive to the symmetry of a
crystal. In a purely semiclassical scenario, periodic systems protected by
either TRS or SIS do not show any second order electric current Hamamoto et
al. [2017]. However, it has been predicted that a second order Hall current
can emerge in noncentrosymmetric TRS invariant systems due to a nonzero Berry
curvature dipole moment Moore and Orenstein [2010]; Sodemann and Fu [2015].
This effect has been theoretically investigated in various TRS invariant
systems You et al. [2018]; Zhang et al. [2018]; Son et al. [2019]; Xu et al.
[2018]; Du et al. [2018]; Facio et al. [2018]; Battilomo et al. [2019]; Lahiri
et al. [2021] and was experimentally verified in bilayer WTe2 Ma et al.
[2019]. In TRS protected systems the nonlinear Hall effect emerges due to the
presence of a non-zero Berry curvature dipole, which requires a low-symmetry
environment. A non-zero dipole has been achieved as a result of either the
tilt of a massive Dirac cone Sodemann and Fu [2015]; Nandy and Sodemann
[2019]; Du et al. [2019], engineering of a finite strain Son et al. [2019], or
a hexagonal warping effect Battilomo et al. [2019]. Recently, disorder induced
nonlinear Hall effect was demonstrated in a TRS protected system Du et al.
[2021a]; Xiao et al. [2019]. On the other hand, a nonlinear planar Hall effect
was experimentally realized on the surface of a three-dimensional nonmagnetic
topological insulator, Bi2Se3, with broken TRS due to an in-plane magnetic
field He et al. [2019]. The spin-momentum locking near the surface together
with the hexagonal warping is responsible for such a nonlinear response Zhang
and Vignale [2018]. In 3D materials, such as Weyl and Dirac semimetals, the
origin of nonlinear planar Hall effect has been attributed to chiral anomaly
Li et al. [2021]; Nandy et al. [2021], which has been extensively studied in
the case of linear planar Hall effect Zheng et al. [2020]; Burkov [2017];
Nandy et al. [2017]; Chen et al. [2018]; He et al. [2014]; Liang et al.
[2015]; Zhang et al. [2016]; Li et al. [2016]; Xiong et al. [2015];
Hirschberger et al. [2016]. Recently, a signature of second order charge-to-
spin conversion due to an unclear mechanism was reported in transition metal
dichalcogenide (TMD)-graphene heterostructures supporting the Rashba-Edelstein
and spin Hall effects Ghiasi et al. [2019].

Here, we propose a mechanism for generating second order charge, spin, and
valley Hall currents in Kane-Mele type 2D topological materials emerging from
the “cross-talk” between Rashba spin-orbit coupling and in-plane magnetic
field-induced Zeeman effect (Fig. 1). We show that the Rashba-Zeeman interplay
creates the necessary asymmetries in space and time to generate such second
order nonlinear responses, and examine how it affects the energy band
structure, the Berry curvature, and the ensuing electronic conductivities. The
Rashba-Zeeman interplay induced transport has been previously studied in
different quantum systems, including a 1D chain Malla and Raikh [2018b], a
spinful realization of Aubry-Andre model Malla and Raikh [2018c], and in a 2D
electron gas Malla and Raikh [2019]. We argue that two-dimensional systems
such as, for instance, siliceneVogt et al. [2012], germanene Dávila et al.
[2014], stanene Zhu et al. [2015], and plumbene Yuhara et al. [2019],
collectively called graphene family materials, as well as antiferromagnetic
manganese chalcogenophosphates (MnPX3, X = S, Se) Li et al. [2013] and
perovskites Liang et al. [2013], are potential platforms to verify this
mechanism. Recently a new class of families known as Jacutingaite has been
shown to host the Kane-Mele type QSH phase Marrazzo et al. [2018, 2019, 2020];
Cucchi et al. [2020]. By tailoring the energy gaps in the Kane-Mele
Hamiltonian of these systems via external interactions, e.g., using static
electric fields or non-resonant circularly polarized light, one can steer the
monolayers through various topological phase transitions. We demonstrate that
only the Berry curvature induced second order Hall current distinguishes
between various topological phases with different Chern numbers, and find a
strong enhancement of Hall current near the phase transitions.

Figure 1: (a) Schematic diagram of the lattice structure of hexagonal
honeycomb monolayer. (b) Representation of Rashba and Zeeman coupling. For
illustration purposes we consider that the sublattices A and B have a buckled
structure, characteristic of two-dimensional topological insulators such as
silicene, germanane, stanene, and plumbene. (c) (Left) Energy band structure
near the $K$ valley in the absence of Rashba and Zeeman coupling. (Center)
Individual effects of each of these interactions in the energy-momentum
dispersion. (Right) Illustration of the effect of the Rashba-Zeeman interplay
in the energy band structure. (d, e) Spin-distributions of the two conduction
bands in momentum space are shown for the cases considered in (c).

##  I Second order nonlinear dynamics in Dirac-like 2D systems

We consider nonlinear effects in light-matter interactions between an ac
electric field and an arbitrary monolayer with unperturbed Hamiltonian
${\hat{H}}_{0}$. The equation of motion for the density matrix operator
$\hat{\rho}(t)$ reads

$i\hbar\ \frac{\partial\hat{\rho}(t)}{\partial
t}=\left[{\hat{H}}_{0}+\hat{H}_{i}(t),\hat{\rho}(t)\right]-\Gamma({\hat{\rho}(t)}-{\hat{\rho}}^{(0)}),$
(1)

where ${\hat{H}}_{i}(t)=-e\phi({\bm{r}},t)$ is the interaction Hamiltonian,
and
$\phi({\bm{r}},t)=\int_{-\infty}^{\infty}d\omega\sum_{\bm{q}}\phi_{\bm{q}\omega}e^{i({\bm{q}}\cdot{\bm{r}}-\omega
t)}$ is the scalar electromagnetic potential, which relates to the electric
field as $\bm{{\cal E}}=-\nabla\phi$. Also,
${\hat{\rho}}^{(0)}=\sum_{lk}f_{lk}|lk\rangle\langle lk|$ is the equilibrium
density matrix. Here, $f_{lk}$ describes the Fermi-Dirac distribution (FD) of
electrons with momentum $\hbar{\bm{k}}$ in the $l$-th band of the unperturbed
Hamiltonian, and $\Gamma$ is the phenomenological decay rate. We solve Eq. (1)
within the framework of perturbation theory by expanding the time-dependent
density matrix in powers of the electric field amplitude
$\hat{\rho}(t)=\sum_{n}{\hat{\rho}}^{(n)}(t)\propto{\cal E}(t)^{n}$, and
compute the $n$-th order two-dimensional current density at ${\bm{r}}_{0}$ as
${\bm{j}}^{(n)}({\bm{r}}_{0},t)=-(e/2)\sum_{ll^{\prime}{\bm{k}}{\bm{k}}^{\prime}}\langle
l^{\prime}{\bm{k}}^{\prime}|\\{{\hat{\bm{v}}\delta({\bm{r}}_{0}-{\bm{r}})\\}_{+}|l{\bm{k}}\rangle\langle
l{\bm{k}}|\hat{\rho}^{(n)}(t)|l^{\prime}{\bm{k}}^{\prime}\rangle}$, where
$\\{{\hat{a}},{\hat{b}}\\}_{+}={\hat{a}}{\hat{b}}+{\hat{b}}{\hat{a}}$. The
velocity operator is defined as
${\hat{\bm{v}}}=-i[{\hat{\bm{r}}},{\hat{H}}_{0}]$. We consider the local limit
${\bm{q}}\rightarrow 0$ by neglecting spatial dispersion since we assume a
normally incident electromagnetic wave, and derive expressions for the optical
conductivity tensor up to second order in the perturbation expansion by
enforcing that the current be proportional to the second power of the incident
field, ${\bm{j}}_{\alpha}=\sigma_{\alpha\beta\gamma}{\cal E}_{\beta}{\cal
E}_{\gamma}$, where $\sigma_{\alpha\beta\gamma}$ is the the optical
conductivity tensor describing a second order dynamics.

The generation of second order optical current involves two virtual electronic
transitions between the energy bands and each transition can take place either
within the same band $l=l^{\prime}$ (intraband) or across two different bands
$l\neq l^{\prime}$ (interband). The second order optical conductivity
$\sigma_{\alpha\beta\gamma}$ can be separated into two terms
${\tilde{\sigma}}_{\alpha\beta\gamma}^{(1)}$ and
${\tilde{\sigma}}_{\alpha\beta\gamma}^{(2)}$ that are either proportional to
the FD functions or their derivatives, respectively (see appendix). The
expressions for ${\tilde{\sigma}}_{\alpha\beta\gamma}^{(1)}$ and
${\tilde{\sigma}}_{\alpha\beta\gamma}^{(2)}$ are

${\tilde{\sigma}}_{\alpha\beta\gamma}^{(1)}(\omega_{1},\omega_{2})=\frac{e^{3}\hbar^{2}}{S}\sum\limits_{l\neq
l^{\prime}\neq
l^{\prime\prime}{\bm{k}}}\frac{\langle\hat{v}_{\alpha}\rangle_{l^{\prime}l}}{\left(E_{l^{\prime}l}+\hbar(\omega_{1}+\omega_{2}+i\Gamma)\right)}\\\
\times\frac{1}{E_{ll^{\prime\prime}}E_{l^{\prime\prime}l^{\prime}}}\Big{[}\frac{\langle\hat{v}_{\beta}\rangle_{ll^{\prime\prime}}\langle\hat{v}_{\gamma}\rangle_{l^{\prime\prime}l^{\prime}}f_{l^{\prime}l^{\prime\prime}}}{E_{l^{\prime}l^{\prime\prime}}+\hbar(\omega_{1}+i\Gamma)}-\frac{\langle\hat{v}_{\gamma}\rangle_{ll^{\prime\prime}}\langle\hat{v}_{\beta}\rangle_{l^{\prime\prime}l^{\prime}}f_{l^{\prime\prime}l}}{E_{l^{\prime\prime}l}+\hbar(\omega_{1}+i\Gamma)}\Big{]}\\\
+\frac{e^{3}\hbar}{S}\sum\limits_{l\neq
l^{\prime}{\bm{k}}}\frac{\langle\hat{v}_{\alpha}\rangle_{l^{\prime}l}f_{l^{\prime}l}}{\left(E_{l^{\prime}l}+\hbar(\omega_{1}+\omega_{2}+i\Gamma)\right)}\\\
\times\frac{\partial}{\partial
k_{\gamma}}\left(\frac{\langle\hat{v}_{\beta}\rangle_{ll^{\prime}}}{E_{ll^{\prime}}\left(E_{l^{\prime}l}+\hbar(\omega_{1}+i\Gamma)\right)}\right),$
(2)
${\tilde{\sigma}}_{\alpha\beta\gamma}^{(2)}(\omega_{1},\omega_{2})={\tilde{\sigma}}_{\alpha\beta\gamma}^{(2,GV)}+{\tilde{\sigma}}_{\alpha\beta\gamma}^{(2,BD)}+{\tilde{\sigma}}_{\alpha\beta\gamma}^{(2,BDC)}=\\\
\frac{e^{3}}{S}\sum\limits_{l{\bm{k}}}\frac{\langle{\hat{v}}_{\alpha}\rangle_{ll}\frac{\partial^{2}f_{l}}{\partial
k_{\beta}\partial
k_{\gamma}}}{\hbar(\omega_{1}+\omega_{2}+i\Gamma)(\hbar(\omega_{1}+i\Gamma))}\\\
+\frac{e^{3}\hbar}{S}\sum\limits_{l\neq
l^{\prime}{\bm{k}}}\frac{\langle\hat{v}_{\alpha}\rangle_{l^{\prime}l}}{\left(E_{l^{\prime}l}+\hbar(\omega_{1}+\omega_{2}+i\Gamma)\right)}\\\
\times\Bigg{[}\frac{\langle\hat{v}_{\gamma}\rangle_{ll^{\prime}}\frac{\partial
f_{l^{\prime}l}}{\partial
k_{\beta}}}{E_{ll^{\prime}}\hbar(\omega_{1}+i\Gamma)}+\frac{\langle\hat{v}_{\beta}\rangle_{ll^{\prime}}\frac{\partial
f_{l^{\prime}l}}{\partial
k_{\gamma}}}{E_{ll^{\prime}}\left(E_{l^{\prime}l}+\hbar(\omega_{1}+i\Gamma)\right)}\Bigg{]},$
(3)

where $S$ is the surface area and
$\sigma_{\alpha\beta\gamma}={\tilde{\sigma}}_{\alpha\beta\gamma}^{(1)}+{\tilde{\sigma}}_{\alpha\beta\gamma}^{(2)}$
is the total conductivity, $\omega_{1}$ and $\omega_{2}$ are frequencies of
the incident electric fields, and $E_{ll^{\prime}}=E_{l}-E_{l^{\prime}}$ with
$E_{l}$ being the energy of the $l$th band at momentum ${\bm{k}}$.

Equation (2) becomes important for larger frequencies and shows spectral
resonances whenever a bandgap coincides with energy of either one photon
$\hbar\omega_{1}$ or two photons $\hbar(\omega_{1}+\omega_{2})$. We restrict
our discussion to near static regime where
${\tilde{\sigma}}_{\alpha\beta\gamma}^{(1)}(\omega_{1},\omega_{2})$ becomes
negligible. The first term ${\tilde{\sigma}}_{\alpha\beta\gamma}^{(2,GV)}$ in
Eq. (3) originates from the charge carrier’s group velocity (GV) $\partial
E_{l}({\bm{k}})/\partial k_{\alpha}$, and a nonzero contribution from it
requires breaking the symmetry of the energy bands $E_{l}({\bm{k}})\neq
E_{l}(-{\bm{k}})$. The second term
${\tilde{\sigma}}_{\alpha\beta\gamma}^{(2,BD)}$ emerges from the nonzero
dipole moment of the Berry curvature in the semiclassical (low frequency)
limit $\mu\gg\hbar\omega$. The Berry dipole (BD) requires asymmetry of the
Berry curvature and is given bySodemann and Fu [2015]

$D_{\alpha
z}=\sum\limits_{l{\bm{k}}}\frac{\partial\Omega_{z,l}({\bm{k}})}{\partial
k_{\alpha}}f_{lk},$ (4)

where $\alpha$ is the dipole direction and
$\Omega_{z,l}({\bm{k}})=2\leavevmode\nobreak\ \text{Im}[\sum_{l^{\prime}\neq
l}(\langle{\hat{v}}_{y}\rangle_{ll^{\prime}}\langle{\hat{v}}_{x}\rangle_{l^{\prime}l})/E_{l^{\prime}l}^{2}]$
is the Berry curvature of $l$th band of a monolayer lying in the $xy$-plane.
Note that $D_{\alpha z}$ is nonzero provided
$\partial\Omega_{z,l}({\bm{k}})/\partial
k_{\alpha}\neq-\partial\Omega_{z,l}(-{\bm{k}})/\partial k_{\alpha}$, which
occurs in 2D materials with a single mirror symmetry line, resulting in a
dipole moment of the Berry curvature in the monolayer’s plane and orthogonal
to the symmetry line Sodemann and Fu [2015]. The third term
${\tilde{\sigma}}_{\alpha\beta\gamma}^{(2,BDC)}$ in Eq. (3) is very similar to
the BD term, however, the denominator includes a third power in energy, and
does not appear within the standard semiclassical approach. This effect has
been previously investigated in the dynamics of Bloch electrons under uniform
electromagnetic fields and magnetically parity-violating systems Gao et al.
[2014]; Watanabe and Yanase [2021] and is typically referred to as intrinsic
Fermi surface effect, since it includes a derivative of the FD function. In
the small frequency regime this term is significantly smaller
($\sim\hbar\omega/\mu$) compared with the BD contribution, and therefore we
call it “Berry dipole correction” (BDC). Finally, we comment that in
centrosymmetric systems the functions within the summations in Eqs. (2) and
(3) have odd parity in momentum, leading to a vanishing second order current.

We comment that our approach can also be extended to both spin and valley
currents. Since, in general the spin will no longer be a good quantum number,
we can use a spin-polarized velocity operator to filter the spin current along
a particular spin polarization. In order to get the spin current density,
denoted by ${\bm{j}}_{\nu}^{(n)}({\bm{r}_{0},t})$, the standard velocity
operator can be replaced by a spin-polarized velocity operator
${\hat{v}}_{\alpha}^{(\nu)}=\\{{\hat{v}}_{\alpha},\hat{s}_{\nu}\\}/4$ in the
expression for electric current density ${\bm{j}}^{(n)}({\bm{r}_{0},t})$,
where $\nu$ is the direction of the spin polarization and $\hat{s}$ is the
spin Pauli matrixHamamoto et al. [2017]. In the following, we will keep our
discussion to spin polarization along $z$-direction, and denote the spin
conductivity as $\sigma_{s,\alpha\beta\gamma}^{(2)}$. The valley conductivity
is obtained by summing the contributions arising from decoupled valleys
weighted by the valley index $\eta$, in a similar fashion as the spin
conductivity is the sum of the conductivities from all Dirac cones weighted by
the spin index $s$,
$\sigma_{\eta,\alpha\beta\gamma}=\sum_{\eta}\eta\sigma_{\alpha\beta\gamma}^{(\eta)}$,
where $\sigma_{\alpha\beta\gamma}^{(\eta)}$ is the conductivity contribution
from the valley $\eta$. Finally, we will also denote the spin-polarized and
valley-polarized Berry curvature as
$\Omega_{z,l}^{(s)},\Omega_{z,l}^{(\eta)}$, and the spin-polarized and valley-
polarized Berry curvature dipole as $D_{\alpha z}^{(s)},D_{\alpha
z}^{(\eta)}$, respectively.

##  II Low-energy Kane-Mele Hamiltonian in 2D topological insulators

We apply the above formalism to a Kane-Mele two-dimensional topological
insulator interacting with a monochromatic ac electric field $\bm{{\cal
E}}(t)$ impinging normally on the monolayer plane (Fig. 1). The low-energy
unperturbed Hamiltonian of the system reads

$\displaystyle{\hat{H}}_{0}^{\eta s}$ $\displaystyle=$ $\displaystyle\hbar
v_{F}(\eta
k_{x}\tau_{x}+k_{y}\tau_{y})+\Delta_{\eta}^{s}\tau_{z}+\frac{\lambda_{R}}{2}(\eta\tau_{x}s_{y}-\tau_{y}s_{x})$
(5) $\displaystyle+$ $\displaystyle\Delta_{B}(\cos\theta s_{x}+\sin\theta
s_{y}),$ (6)

where the particle momentum is measured from the Dirac points $K$ and
$K^{\prime}$ and $v_{F}$ is the Fermi velocity. Pauli matrices $\tau_{i}$
correspond to sublattice pseudospins ($A$, $B$), while $s_{i}$ describe the
electron spins. The spin and valley dependent Dirac gap is
$\Delta_{\eta}^{s}=\eta s\lambda_{SO}$, where the strength of the SOC
$\lambda_{SO}$ depends on the buckling of the honeycomb lattice as well as the
atomic size. For instance, it is found to be $3.9$ meV for silicene, $20$ meV
for germanene, $300$ meV for stanene, and $0.4$ eV for plumbeneEzawa [2015];
Rodriguez-Lopez et al. [2017]; Yuhara et al. [2019]. The third term in Eq. (5)
represents the Rashba coupling $\lambda_{R}$ (Fig. 1b), which originates from
the nearest neighbor hopping in the tight-binding model and requires external
parameters such as static electric fieldBychkov and Rashba [1984]; Tse et al.
[2011], metal-atom adsorptionShikin et al. [2008]; Weeks et al. [2011], or
presence of a substrateRashba [2009], which breaks SIS in the normal
direction. Finally, we consider a Zeeman interaction via the last term in Eq.
(5) by applying an in-plane magnetic field ${\bf B}$
$(\Delta_{B}=g_{L}\mu_{B}|B|/2)$ arbitrarily oriented along the plane of the
sample with an angle $\theta$ measured from the $x$-axis (Fig. 1b). Here
$g_{L}$ is the Lande g factor and $\mu_{B}$ is the Bohr magneton. Equation (5)
without the Rashba-Zeeman terms, accounts for the centrosymmetric part of the
full tight-binding KM Hamiltonian. It is protected by both TRS, ${\cal
T}{\hat{H}_{\eta}}({\bm{k}}){\cal T}^{-1}={\hat{H}}_{\eta}(-{\bm{k}})$, and
SIS, ${\cal P}{\hat{H}_{\eta}}({\bm{k}}){\cal
P}^{-1}={\hat{H}}_{-\eta}(-{\bm{k}})$, where the symmetry operators are given
as ${\cal T}=i\gamma_{0}\tau_{1}s_{2}{\cal K}$ and ${\cal
P}=\gamma_{1}\tau_{1}s_{0}$, and $\gamma$, $\tau$, and $s$ are Pauli matrices
associated with the valley, sublattice, and spin degrees of freedom, and
${\cal K}$ is the complex conjugation operator. In this case, the energy band-
structure is degenerate in both valley and spin (Fig. 1c) and the expectation
value of the spin operator is orthogonal to the monolayer plane and locked to
the direction of motion of the charge carriers (Fig. 1d).

The Rashba and Zeeman terms break respectively SIS and TRS and, when
simultaneously applied to the monolayer, lift the the degeneracy of the energy
bands $E(-{\bm{k}})\neq E({\bm{k}})$. In the presence of a nonzero
$\lambda_{R}$, the energies describing the conduction ($E_{+}$) and valence
($E_{-}$) bands read
$E_{+(-)}=\pm\lambda_{R}/2+(-)[\hbar^{2}v_{F}^{2}k^{2})+\left(\lambda_{SO}\mp(\lambda_{R}/2)\right)^{2}]^{1/2}$.
For momentum $k=0$, the energies of the valence bands are shifted by
$-\lambda_{SO}\pm\lambda_{R}$ while the conduction bands remain unaffected,
which leads to new energy gaps equal to $2|\lambda_{SO}\pm\lambda_{R}|$ see
Fig. 1c. The Zeeman field affects the energy band structure as
$E=\pm[\left(\hbar v_{F}k\pm\Delta_{B}\right)^{2}+\lambda_{SO}^{2}]^{1/2}$,
but does not change the magnitude of the energy gaps. On the other hand, the
minimum energy is shifted radially in momentum space and lies along the
contour $\hbar v_{F}k=|\Delta_{B}|$, as shown in Fig. 1c. The corresponding
effects for these interaction terms in the spin distribution in momentum space
are shown in Fig. 1d,e. While the Rashba term results in a spin-momentum
locking mechanism that leads the in-plane components of the spin to be
orthogonal to the momentum, the Zeeman term tends to orient the spin parallel
to the direction of the applied in-plane magnetic field.

The independent violation of SIS or TRS via Rashba or Zeeman interaction does
not break the even parity of the energy band, i. e.,
$E(-{\bm{k}})=E({\bm{k}})$, and therefore Rashba and Zeeman couplings cannot
individually lead to any energy band dependent second order electric current.
Rather, both contributions are simultaneously necessary to enable a spin-to-
charge conversion mechanism that in the nonlinear regime results in an
emerging nonlinear Hall currentDu et al. [2021b]; Ortix . Assuming that the
direction of the magnetic field is along the $x$-axis ($\theta=0$), one can
show that the eigenenergies of the unperturbed Hamiltonian satisfy a quartic
order equation, which includes a term proportional to
$\lambda_{R}\Delta_{B}k_{y}$. Because this term is odd in $k_{y}$, it results
in a band structure that does not preserve momentum inversion symmetry. Note
that nonzero Zeeman and Rashba couplings are needed for this term to be
relevant, further confirming that both SIS and TRS need to be broken to enable
new effects in the KM model Hamiltonian. We mention that the energy dispersion
along $k_{x}$, however, remains symmetric, $E_{l}(k_{x})=E_{l}(-k_{x})$. Note
also that when both interactions are turned on, the spin distribution due to
each valley is distorted in the direction orthogonal ($k_{y}$) to the magnetic
field (Fig. 1d). This asymmetry involving spins with opposite momentum is
crucial to create unbalanced spin currents, which then results in charge
transport when the sample is driven by an electromagnetic waveHe et al.
[2019].

##  III Rashba-Zeeman enabled electronic, spin, and valley transport

The interplay between the Rashba coupling $\lambda_{R}$ and Zeeman energy
$\Delta_{B}$ has a strong effect on all the three terms of the conductivity
tensor in Eq. (3), and we analyze each of them in detail. We assume that the
Zeeman coupling is due to a magnetic field applied along the $x$-direction,
and we consider a monochromatic incident field of frequency $\omega$ and focus
on second harmonic generation processes, i.e., $\omega_{1}=\omega_{2}=\omega$.

###  III.1 Group velocity contribution

Since the interplay between Rashba and Zeeman interactions does not affect the
energy and spin distributions in the $x$ direction, all GV conductivity terms
with an odd number of ‘$x$’s must vanish. For concreteness, in the following
we will discuss the non-zero $\\{\alpha\beta\gamma\\}=\\{yyy\\},\\{yxx\\}$
conductivity tensor components. The former (latter) describes the longitudinal
(Hall) current generated in the system when the ac electromagnetic wave is
applied orthogonal (parallel) to the static in-plane magnetic field. In both
cases, the currents are perpendicular to the Zeeman field direction. The
conductivity tensors ${\tilde{\sigma}}_{yyy}^{(2,GV)}$ and
${\tilde{\sigma}}_{yxx}^{(2,GV)}$ are plotted as a function of chemical
potential for a fixed frequency and $\Delta_{B}=0.5\lambda_{SO}$, see Fig. 2.
Note that both conductivities show non-monotonic behavior with a maximum near
$\mu_{0}\approx 1.11\lambda_{SO}$. The magnitude of this maximum increases for
higher values of the Rashba coupling because this corresponds to larger
asymmetries of the energy dispersion. The position of the maximum conductivity
can be explained by closely looking at the energy dispersion of the conduction
bands. When the chemical potential is below $\mu_{0}$, only the lowest
conduction band contributes to the conductivity tensor, while above $\mu_{0}$,
a second contribution appears from the highest conduction band. Each
conduction band, however, gives contributions of opposite sign and the
cumulative effect results in a decrease of the total conductivity for
$\mu>\mu_{0}$. For values of the chemical potential much larger than the spin-
orbit coupling one notices that the conductivities approach zero. This is a
consequence of the fact that more energetic charge carriers, which dominate
the electronic transport in this regime, are not sensitive to the broken SIS
and TRS that occur at lower energies near the bandgap.

The physical origin of the GV contribution to a second order conductivity
arises from the Rashba-Zeeman induced distortion of Fermi contours. Rashba
coupling alone introduces spin-momentum locking in the system and the spins at
opposite wave vectors are oriented in opposite directions (Fig. 1c, second
panel) and are equally populated. This introduces a second order spin current
in the system, i.e., the spins oriented perpendicular to the electric field
give a longitudinal spin current and spins oriented parallel (or anti-
parallel) to the field generate a transverse spin current. The charge current,
however, vanishes since the spins with opposite orientations move in opposite
directions. Note that the spin current here refers to the spins with
orientation along the plane of the sample. When a Zeeman field is introduced,
in addition to the Rashba coupling, it distorts the Fermi contours and creates
population imbalance of opposite spins at opposite wave vectors (Fig. 1c, last
panel). This results in an asymmetry in the number of charges moving in the
opposite directions, and leads to a second order nonzero charge current along
both longitudinal and transverse directions.

We refrain from showing the GV contribution to the spin and valley
conductivities since those are zero for the case under consideration. Indeed,
the $z$-polarized spin conductivity
${\tilde{\sigma}}_{s,\alpha\beta\gamma}^{(2,GV)}$ is nonzero only if the up-
spin current and down-spin current are different. Because the Rashba-Zeeman
interplay only creates asymmetry among the planar components of the spins, one
can show that the second order current due to the up spins cancels the
contribution from down spins, therefore resulting in
${\tilde{\sigma}}_{s,\alpha\beta\gamma}^{(2,GV)}=0$. Regarding valley
transport, although the asymmetric group velocity of the charge carriers
individually leads to nonzero valley conductivities
${\tilde{\sigma}}_{\eta,\alpha\beta\gamma}^{(2,GV)}$ at $K$ and $K^{\prime}$
Dirac points, their opposite value results in a perfect cancellation of the
total valley conductivity.

Figure 2: Group velocity contribution to the second order intraband
longitudinal $\tilde{\sigma}_{yyy}^{(2),GV}(2\omega)$ and (inset) transverse
$\tilde{\sigma}_{yxx}^{(2),GV}(2\omega)$ electronic conductivities for
$\Delta_{B}=0.5\lambda_{SO}$. The frequency dependence is captured in the
normalization parameter
${\sigma_{0}}=e^{3}v_{F}/\left[S\hbar(2\omega+i\Gamma)(\omega+i\Gamma)\right]$.

###  III.2 Berry dipole contribution

Next, we address the effect of Rashba-Zeeman interplay on charge, spin-
polarized and valley-polarized Berry curvature profiles. The BD contribution
to the electronic, spin, and valley transport has an anomalous origin and may
generate a second order transverse current when the incident ac field has a
component perpendicular to the Zeeman magnetic field. For the Hamiltonian in
consideration in Eq. (5), we observe that only the second order spin-transport
survives while the charge and valley transports vanish. Indeed, the Berry
curvatures for up and down spins are the negative of each other and thus the
total BD per Dirac cone vanishes, resulting in a zero second order
conductivity. Similarly, the valley-polarized BD at $K$ and $K^{\prime}$ are
exactly zero, which prevents valley transport to take place. On the other
hand, the spin-polarized Berry curvature takes the spin weight into account
and the contributions from up and down spins add up. This results in a finite
spin-polarized Berry dipole profile at a particular Dirac cone and a nonzero
spin transport.

Figure 3: (a) Momentum distribution of the spin-polarized Berry curvature
$\Omega_{z,l}^{(s,c)}({\bm{k}})$ and (b) corresponding derivative
$\partial\Omega_{z,l}^{(s,c)}({\bm{k}})/\partial k_{y}$ for $\Delta_{B}=0$,
$\lambda_{R}=0$. Index ‘c’ refers to contributions from the two conduction
bands. (c), (d) Same as (a) and (b) for $\Delta_{B}=0.5\lambda_{SO}$,
$\lambda_{R}=0.3\lambda_{SO}$.

In Fig. 3 we plot the momentum dependence of spin-polarized Berry curvature
profile
$\Omega_{z,l}^{(s,c)}({\bm{k}})=\sum_{l=3,4}\Omega_{z,l}^{(s)}({\bm{k}})$ and
its corresponding derivative $\partial\Omega_{z,l}^{(s,c)}({\bm{k}})/\partial
k_{y}$ for the conduction bands to demonstrate how they are affected by the
Zeeman-Rashba coupling. When both $\Delta_{B}$ and $\lambda_{R}$ are zero
(panel 3(a)) we notice that spin-polarized Berry curvature is symmetric and
its corresponding derivative is an odd function of $k_{y}$ (panel 3(b)),
therefore the Berry dipole moment vanishes. When the Rashba-Zeeman interplay
is turned on, however, it creates an asymmetry in the Berry curvature as shown
in Fig. 3c, where we observe that the center of the distribution is shifted
along $k_{y}$-direction. As a result, the corresponding derivative in Fig. 3d
is no longer an odd function, and this guarantees that the integral in
momentum space appearing in Eq. (4) is nonzero, and thus the BD moment is
nonzero. As discussed before, the third contribution of (3),
${\tilde{\sigma}}_{s,xyy}^{(2,BDC)}$, follows the same principles as BD
contribution and acts as a small correction to the BD contribution of the
conductivity. In Fig. 4, the contributions due to the BD to the spin-polarized
conductivity $\tilde{\sigma}_{s,\alpha\beta\gamma}^{(2)}$ are shown for a
fixed value of the Zeeman interaction ($\Delta_{B}=0.5\lambda_{SO}$) and for
three values of Rashba coupling, as in Fig. 2. The BD-enabled spin-
conductivity shows a non-monotonic behavior similar to that seen in the GV
contribution term, and at large chemical potential the conductivity saturates
near zero. We also observe a peculiar feature for chemical potential near
$\mu_{0}$, which arises due to the numerical derivative taken near the
crossing of two conduction bands.

Figure 4: Conductivities $\tilde{\sigma}_{s,xyy}^{(2),BD}(2\omega)$ for
$\Delta_{B}=0.5\lambda_{SO}$. The frequency dependence is captured in the
parameter ${\sigma_{1}}=e^{3}\hbar v_{F}/\left[S(\omega+i\Gamma)\right]$.

###  III.3 Interplay between Berry dipole moment and tailored topological
phase transitions

The Dirac half-gaps $\Delta_{\eta}^{s}$ in the Hamiltonian in Eq. (5) can be
externally tailored to host various topological phases. For example, for a 2D
Kane-Mele topological insulator with two staggered sub-lattices a static
electric field can be used to generate on-site potentials that will modify the
gap as $\Delta_{\eta}^{s}=\eta s\lambda_{SO}-\lambda_{E}$, where $\lambda_{E}$
is the coupling with the field. Analogously, a high-frequency non-resonant
circularly polarized laser can induce Floquet topological states in the
monolayer, providing an additional knob for tailoring the energy band
structure via asymmetric valley coupling, i.e., the Dirac gap change as
$\Delta_{\eta}^{s}=\eta s\lambda_{SO}-\lambda_{E}-\eta\lambda_{L}$, where
$\lambda_{L}$ represents the laser induced interaction. The topological phases
accessible via these driving fields can be characterized by topological
invariants, namely the Chern ${\cal C}=\sum_{\eta,s}{\cal C}_{s}^{\eta}$, spin
Chern ${\cal C}_{s}=\sum_{\eta,s}s\leavevmode\nobreak\ {\cal C}_{s}^{\eta}/2$,
valley Chern ${\cal C}_{\eta}=\sum_{\eta,s}\eta\leavevmode\nobreak\ {\cal
C}_{s}^{\eta}$, and spin-valley Chern $C_{s\eta}=\sum_{\eta,s}\eta
s\leavevmode\nobreak\ {\cal C}_{s}^{\eta}/2$ numbers. Here ${\cal
C}_{s}^{\eta}=\eta\leavevmode\nobreak\ \text{sign}[\Delta_{s}^{\eta}]/2$. By
varying parameters $\lambda_{E}$ and $\lambda_{L}$, one can map a 2D
topological phase diagram, which hosts electronic phases, including quantum
spin Hall insulator (QSHI), band insulator (BI), anomalous quantum Hall
insulator (AQHI), and polarized spin quantum Hall insulator (PS-QHI). These
phases are shown in Fig. 5a and correspond to systems with all Dirac gaps
open. Topological phase transitions take place at the borderlines between
these phases (black solid curves) which mark the closing of at least one Dirac
cone, while the two points
$(\lambda_{E}/\lambda_{SO},\lambda_{L}/\lambda_{SO})=(1,0),(0,1)$ represent
closing of two mass gaps. As we discuss next, the asymmetry introduced by
modifying the Dirac gaps with external agents enables BD-induced charge, spin,
and valley transport conductivities to be nonzero.

For concreteness, here we will focus on the conductivity tensor
${\tilde{\sigma}}_{xyy}^{(2,BD)}$, since it originates from the Berry
curvature dipole moment and includes signatures of various topological
properties of the system. In Fig. 5b, we show the phase diagram of
${\tilde{\sigma}}_{xyy}^{(2,BD)}$ for fixed values of chemical potential
$\mu=5\lambda_{SO}$, Zeeman interaction $\Delta_{B}=0.5\lambda_{SO}$, and
Rashba coupling $\lambda_{R}=0.3\lambda_{SO}$. Away from the phase boundaries
the conductivity shows nonlinear dependence on the magnitude of the Dirac gaps
and is clearly seen in the PS-QHI phase, where the red and blue colors
coexist, even though the sign of all the Dirac gaps remains the same. Near the
topological phase transitions the conductivity is dominated by the smallest
mass gap. The change of color near the phase boundaries reflects a change in
Chern number, and a decrease in Chern number is associated with a color change
from red to blue. This color change also marks a change in the Berry dipole
direction. This behavior is consistent across all phase boundaries and can be
better understood by comparing with the phase diagram in Fig. 5a. Although the
phase diagram of the conductivity demonstrates qualitative signatures of
topological phase transitions and a decrease in Chern number, the magnitude of
this change is difficult to estimate.

A quantitative analysis of the second order transport properties across
various phase transitions is shown in Fig. 5c,d. There we plot the
conductivities corresponding to charge-, spin-, and valley-currents along
three paths in the phase diagram, namely: (A) QSHI$\rightarrow$AQHI (B)
QSHI$\rightarrow$BI (C) AQHI$\rightarrow$PS-QHI$\rightarrow$BI, as shown by
marked arrows in Fig. 5a. Along path (A), the Chern numbers $\\{{\cal C},{\cal
C}_{s},{\cal C}_{\eta}\\}$ change as $\\{0,1,0\\}\rightarrow\\{-2,0,0\\}$, as
we move from the QSHI to AQHI phase. Both charge and spin-polarized
conductivities show similar behavior near the phase transition point
$(\lambda_{E}/\lambda_{SO},\lambda_{L}/\lambda_{SO})=(0,1)$, with a resonance-
antiresonance shaped behavior, indicating the change of Berry dipole direction
across the phase transition. The sign of Berry curvature near a particular
Dirac cone is dependent on the sign of that Dirac gap. When one moves across a
phase boundary the Dirac gap closes and reopens and this change in sign is
reflected in the Berry curvature and an increase/decrease in the Chern number.
This sign change is subsequently transferred to the Berry dipole. Furthermore,
we observe that near the phase transitions the conductivities show strong
enhancement in magnitude due to the fact that the inter-subband gap becomes
small and the Berry dipole increases as $1/|\Delta_{\eta}^{s}|^{2}$. Note that
the valley conductivity vanishes along this path. The similarity between
charge and spin conductivity is due to the fact that both charge and spin
Chern numbers decrease in magnitude. Furthermore, we observe that the
magnitude of the peaks across a phase transition for the charge conductivity
is nearly two times that of the spin conductivity. This quantitative
difference in the conductivities demonstrates that the charge- and spin-
dipoles strongly respond to the topological phase transition and this response
inherits the quantitative information about the change in Chern numbers.

Figure 5: (a) Topological phase diagram of the KM Hamiltonian without the
Rashba and Zeeman terms. (b) The phase diagram of
$\tilde{\sigma}_{xyy}^{(2,BD)}$ as a function of $\lambda_{E}$ and
$\lambda_{L}$ in the presence of Rashba-Zeeman couplings. The conducitvity is
normalized with respect to the corresponding maximum value. (c, d)
Conductivities $\tilde{\sigma}_{xyy}^{(2,BD)}$ (black),
$\tilde{\sigma}_{s,xyy}^{(2,BD)}$ (blue), and
$\tilde{\sigma}_{\eta,xyy}^{(2,BD)}$ (red) (c) along paths A
($\lambda_{P}=\lambda_{E}$), B ($\lambda_{P}=\lambda_{L}$), and (d) C, as
shown in panel (a). Note that, for path B, the red color represents the charge
conductivity and the black color represents the valley conductivity. The
conductivities are evaluated for fixed values of $\Delta_{B}=0.5\lambda_{SO}$,
$\lambda_{R}=0.3\lambda_{SO}$, and $\mu=5\lambda_{SO}$. The conductivities in
(c, d) are normalized by $\sigma_{1}S$.

The fact that the conductivities and BD behavior inherit topology signatures
near phase transitions can be seen across all the phase boundaries. For
example, along path (B) $\lambda_{L}=0$ and the spin conductivity behavior
remains the same as path (A), while the behavior of charge and valley
conductivities are switched. This is indeed in line with the fact that the
change in charge Chern number across path (A) is equivalent to the change in
valley Chern number across path (B). In order to further strengthen our
argument, we also consider a more interesting path (C), where the system
undergoes two phase transitions. The charge Chern number increases across both
boundaries, while the valley Chern number decreases across both boundaries.
This is shown in Fig. 5d, where the charge conductivity is the negative of the
valley conductivity, and each conductivity shows the same behavior across both
the boundaries. The spin Chern number, on the other hand, increases across the
first boundary and then decreases across the second boundary. The
corresponding conductivity behavior near the first boundary is opposite of the
behavior at the second boundary. The above observations and arguments
demonstrate clear signatures of topology in the conductivity near the phase
transitions and could be extended to a three-dimensional topological phase
diagram by adding another parameter to the mass term $\Delta_{\eta}^{s}$Malla
and Kort-Kamp [2021].

Next, we compare ${\tilde{\sigma}}_{xyy}^{(2,BD)}$ evaluated at two different
points in the phase diagram,
$(\lambda_{E},\lambda_{L})/\lambda_{SO}=(0.3,1.5)$ (AQHI) and
$(\lambda_{E},\lambda_{L})/\lambda_{SO}=(1.5,0.3)$ (BI). These two points have
the same band structure, however the signs of two of the mass gaps are
different, and therefore they belong to different topological phases. To
understand how the difference between the two phases affects the second order
electronic conductivity of the system, we plot contributions from $K$ and
$K^{\prime}$ valleys separately (and the total conductivity) as a function of
$\mu$ in Fig. 6a (AQHI) and Fig. 6b (BI). Note that, the contributions from
$K$ valley are identical in both the phases, however that associated to the
$K^{\prime}$ valley for the AQHI phase is exactly the negative of the
contribution in the BI phase. The sharp changes in the conductivity occur when
$\mu$ coincides with one of the Dirac half-gaps $\Delta_{\eta}^{s}$. For small
$\mu$, the conductivity only comes from the smallest gap, and the sign of that
Dirac gap becomes equally important. Finally, in Fig. 6c we plot the
magnitudes of the conductivities $|{\tilde{\sigma}}_{xyy}^{(2,BD)}|$ and
$|{\tilde{\sigma}}_{yxx}^{(2,BD)}|$ for the point considered in Fig 6a as a
function of the angle $\theta$ of the magnetic field as defined in Eq. (5).
The electric current is maximum when the electric field and the magnetic field
are perpendicular to each other. The transverse nature of current suggests
that the current flows parallel to the magnetic field (perpendicular to the
electric field). Note that, the magnitude of the conductivities corresponding
to points $(\lambda_{E},\lambda_{L})/\lambda_{SO}=(0.3,1.5)$ and
$(\lambda_{E},\lambda_{L})/\lambda_{SO}=(1.5,0.3)$ are exactly equal, since
they represent identical band structures.

Figure 6: (a, b) The $\mu$ dependence of the charge conductivities
${\tilde{\sigma}}_{xyy}^{(2,BD)}$ for
$(\lambda_{E},\lambda_{L})/\lambda_{SO}=(0.3,1.5)$ (AQHI) and
$(\lambda_{E},\lambda_{L})/\lambda_{SO}=(1.5,0.3)$ (BI), respectively. (c) The
polar plot of $|{\tilde{\sigma}}_{xyy}^{(2,BD)}|$ (blue) and
$|{\tilde{\sigma}}_{yxx}^{(2,BD)}|$ (red) are shown for the point in the AQHI
phase as considered in (a). The angle of the magnetic field is shown in
degree. The conductivities are normalized by $\sigma_{1}S$,
$\mu=5\lambda_{SO}$.

##  IV Discussions and Conclusion

In summary, we have showed that massive Dirac fermions in centrosymmetric KM
model can generate a second order transverse current when an out-of-plane
symmetry broken Rashba coupling coexists with an in-plane magnetic field. We
comment that the full KM tight-binding model also includes non-centrosymmetric
part which in the low-energy regime can be included via the trigonal warping
term near the Dirac conesBattilomo et al. [2019]. A second order current due
to this warping term has been demonstrated in various honeycomb
systemsMargulis et al. [2013]; Golub and Tarasenko [2014]. Comparing with our
study, we note that the asymmetry in energy spectra created by Rashba Zeeman
interplay is proportional to $k_{\alpha}^{1/2}$, while the asymmetry created
by the warping term is proportional to $k_{\alpha}^{3/2}$. So in a full KM
model the Rashba-Zeeman induced second order term should dominate the
optoelectronic response near the Dirac point, while the warping induced second
order current dominates away from it. We also mention that a second kind of
Rashba coupling ($\lambda_{R2}$) of intrinsic origin can appear in a typical
KM material. Nevertheless, the strength of this coupling is proportional to
the next-nearest-neighbor hopping parameter, therefore being often very small
compared to $\lambda_{SO}$. Hence, we neglected it from our analysis.

In the quasi-static limit, we have demonstrated that the intraband current
includes contributions from the asymmetry in the energy dispersion as well as
from the asymmetry created in the Berry curvature. The GV contribution
generates a second order current along longitudinal as well as the transverse
direction, and remains normal to the Zeeman field. The BD contribution,
however, only produces transverse current and is nonzero only when a component
of the electric field is perpendicular to the Zeeman field. Since the
transverse contributions from the GV and BD terms are perpendicular to each
other, both contributions can be experimentally probed separately. We briefly
comment on the case of large frequency limit
$\hbar\omega\sim\Delta_{\eta}^{s}$, where the interband conductivity
${\tilde{\sigma}}_{\alpha\beta\gamma}^{(1)}$ dominates. When $\mu$ lies in the
middle of the band gap, interband transitions become the only contributing
factor, but the second order nonlinear response in this case is zero in the
absence of Rashba-Zeeman interplay Malla and Kort-Kamp [2021]. When both
Rashba and Zeeman couplings are present the asymmetry in the energy bands
generates a resonance-induced second harmonic electric current perpendicular
to the in-plane magnetic field. The material response then shows two spectral
resonances corresponding to frequencies $\hbar\omega=2|\Delta_{\eta}^{s}|$ and
$\hbar\omega=|\Delta_{\eta}^{s}|$. The resonance corresponding to $2\omega$
dominates since it resonates with a smaller gap. The Berry curvature
contribution does not affect the interband current at the origin of the phase
diagram, but becomes significant away from it.

We also demonstrated that the topology of the low-energy Hamiltonian can be
tailored via external parameters and the asymmetry between the Dirac gaps for
different spins and valleys can play a significant role in the second order
transverse current. We also show that the the Berry curvature induced dipole
moment strength increases for Dirac cones with the smallest gap and the dipole
changes sign across phase transitions. The recent synthesis of two-dimensional
buckled members of graphene family, including silicene, germanene, stanene,
and plumbene, as well as topologically non-trivial antiferromagnetic manganese
chalcogenophosphates and perovskites provide a perfect platform to
experimentally investigate second order Hall transport effects in Kane-Mele
two-dimensional topological insulators. The ability to control each Dirac gap
in these monolayers makes them a unique all-in-one material for opto-
electronics as well as spintronics based applications.

##  Appendix A Method

We have followed the method prescribed in Ref. Mikhailov, 2016, and included
the electric field in the Hamiltonian as a scalar potential
$\phi({\bm{r}},t)$. The equation of motion (1) is solved perturbatively in the
energy eigen-basis $|\lambda\rangle=|l{\bm{k}}\rangle$ and then the expression
for the electric current is
${\bm{j}}^{(n)}({\bm{r}}_{0},t)=-e/2\sum_{ll^{\prime}{\bm{k}}{\bm{k}}^{\prime}}\langle
l^{\prime}{\bm{k}}^{\prime}|\\{\hat{{\bm{v}}}\delta({\bm{r}}_{0}-{\bm{r}})\\}_{+}|l{\bm{k}}\rangle\langle
l{\bm{k}}|\hat{\rho}^{(n)}(t)|l^{\prime}{\bm{k}}^{\prime}\rangle$ is used to
evaluate the second order current. The first order solution leads to the Kubo
formula for the linear conductivity. The second-order density matrix can be
obtained from the first order density matrix and is expressed as

$\displaystyle\langle\lambda|\hat{\rho}^{(2)}|\lambda^{\prime}\rangle_{t}=e^{2}\int_{-\infty}^{\infty}d\omega_{1}d\omega_{2}e^{-i(\omega_{1}+\omega_{2})t}$
$\displaystyle\hskip
56.9055pt\times\sum\limits_{\lambda^{\prime\prime}}\frac{\langle\lambda|\phi_{\omega_{2}}({\bm{r}_{2}})|\lambda^{\prime\prime}\rangle\langle\lambda^{\prime\prime}|\phi_{\omega_{1}}({\bm{r}_{1}})|\lambda^{\prime}\rangle}{E_{\lambda^{\prime}}-E_{\lambda}+\hbar(\omega_{1}+\omega_{2}+i\Gamma)}$
$\displaystyle\times\Bigg{[}\frac{f_{\lambda^{\prime}}-f_{\lambda^{\prime\prime}}}{E_{\lambda^{\prime}}-E_{\lambda^{\prime\prime}}+\hbar(\omega_{1}+i\Gamma)}-\frac{f_{\lambda^{\prime\prime}}-f_{\lambda}}{E_{\lambda^{\prime\prime}}-E_{\lambda}+\hbar(\omega_{2}+i\Gamma)}\Bigg{]},$
(7)

and the corresponding expression for the second order current is given by

${\bm{j}}^{(2)}({\bm{r}}_{0},t)=-\frac{e^{3}}{2S}\int_{-\infty}^{\infty}d\omega_{1}d\omega_{2}\sum_{\tilde{q},{\bm{q}}_{1},{\bm{q}}_{2}}\phi_{{\bm{q}}_{1}\omega_{1}}\phi_{{\bm{q}}_{2}\omega_{2}}\\\
\times\sum_{\lambda,\lambda^{\prime}}\frac{\langle\lambda^{\prime}|\\{\hat{{\bm{v}}},e^{-i\tilde{{\bm{q}}}.{\bm{r}}_{0}}\\}_{+}|\lambda\rangle
e^{i\tilde{{\bm{q}}}.{\bm{r}}_{0}-i(\omega_{1}+\omega_{2})t}}{E_{\lambda^{\prime}}-E_{\lambda}+\hbar(\omega_{1}+\omega_{2}+i\Gamma)}\\\
\times\langle\lambda|e^{i{\bm{q}}_{2}.{\bm{r}}_{2}}|\lambda^{\prime\prime}\rangle\langle\lambda^{\prime\prime}|e^{i{\bm{q}}_{1}.{\bm{r}}_{1}}|\lambda^{\prime}\rangle\\\
\times\Bigg{[}\frac{f_{\lambda^{\prime}}-f_{\lambda^{\prime\prime}}}{E_{\lambda^{\prime}}-E_{\lambda^{\prime\prime}}+\hbar(\omega_{1}+i\Gamma)}-\frac{f_{\lambda^{\prime\prime}}-f_{\lambda}}{E_{\lambda^{\prime\prime}}-E_{\lambda}+\hbar(\omega_{2}+i\Gamma)}\Bigg{]}.$
(8)

Note that, the current does not change if we switch $(q_{2},\omega_{2},r_{2})$
with $(q_{1},\omega_{1},r_{1})$ in the second term of the last line in (8).

Now using the matrix element
$\langle\lambda|e^{i{\bm{q}}_{1}.{\bm{r}}_{1}}|\lambda^{\prime}\rangle=\delta_{l,l^{\prime}}-q_{1}M_{l,l^{\prime}}$,
we can write the conductivity in a form which shows all types of virtual
transitions involved in the second order process, and the expression reads

${\bm{j}}^{(2)}({\bm{r}}_{0},t)=-\frac{e^{3}}{2S}\int_{-\infty}^{\infty}d\omega_{1}d\omega_{2}\sum_{{\bm{q}}_{1},{\bm{q}}_{2}}\phi_{{\bm{q}}_{1}\omega_{1}}\phi_{{\bm{q}}_{2}\omega_{2}}e^{i({\bm{q}}_{1}+{\bm{q}}_{2}).{\bm{r}}_{0}-i(\omega_{1}+\omega_{2})t}\sum_{l,l^{\prime},l^{\prime\prime},{\bm{k}}}\frac{\langle
l^{\prime},{\bm{k}}|\\{\hat{{\bm{v}}},e^{-i({\bm{q}}_{1}+{\bm{q}}_{2}).{\bm{r}}_{0}}\\}_{+}|l,\bm{k+q_{1}+q_{2}}\rangle}{E_{\lambda^{\prime}}-E_{\lambda}+\hbar(\omega_{1}+\omega_{2}+i\Gamma)}\\\
\times\Bigg{[}(\delta_{l,l^{\prime\prime}}-q_{2,\gamma}M_{l,l^{\prime\prime},\bm{k+q_{1}}})(\delta_{l^{\prime\prime},l^{\prime}}-q_{1,\beta}M_{l^{\prime\prime},l^{\prime},{\bm{k}}})\frac{f_{l^{\prime},{\bm{k}}}-f_{l^{\prime\prime},\bm{k+q_{1}}}}{E_{l^{\prime},{\bm{k}}}-E_{l^{\prime\prime},\bm{k+q_{1}}}+\hbar(\omega_{1}+i\Gamma)}-\\\
(\delta_{l,l^{\prime\prime}}-q_{1,\beta}M_{l,l^{\prime\prime},\bm{k+q_{2}}})(\delta_{l^{\prime\prime},l^{\prime}}-q_{2,\gamma}M_{l^{\prime\prime},l^{\prime},{\bm{k}}})\frac{f_{l^{\prime\prime},\bm{k+q_{2}}}-f_{l,\bm{k+q_{1}+q_{2}}}}{E_{l^{\prime\prime},\bm{k+q_{2}}}-E_{l,\bm{k+q_{1}+q_{2}}}+\hbar(\omega_{1}+i\Gamma)}\Bigg{]}.$
(9)

The delta-functions in Eq. (9) represent virtual intraband transitions, while
the terms with matrix elements $M_{l,l^{\prime},{k_{\alpha}}}=2\hbar\langle
v_{\alpha}\rangle_{l,l^{\prime}}/(E_{l}-E_{l^{\prime}})$ represent virtual
interband transitions. We can separate terms depending on the type of virtual
transitions, and we consider three cases, (i) $l=l^{\prime}=l^{\prime\prime}$
(two intraband transitions) (ii) $l^{\prime\prime}=l\neq l^{\prime}$ or
$l^{\prime\prime}=l^{\prime}\neq l$ (one intraband and one interband
transitions) (iii) $l\neq l^{\prime}\neq l^{\prime\prime}$ (two interband
transitions).

Case i: When virtual transitions occur within the same band only the delta-
functions survive in Eq. (9), and we simplify the expression inside the square
bracket in (9) to

$\frac{f_{l,{\bm{k}}}-f_{l,\bm{k+q_{1}}}}{E_{l,{\bm{k}}}-E_{l,\bm{k+q_{1}}}+\hbar(\omega_{1}+i\Gamma)}\\\
-\frac{f_{l,\bm{k+q_{1}}}-f_{l,\bm{k+q_{1}+q_{2}}}}{E_{l,\bm{k+q_{1}}}-E_{l,\bm{k+q_{1}+q_{2}}}+\hbar(\omega_{1}+i\Gamma)},$
(10)

using the identities $f_{l,k}-f_{l,k+q_{1}}\approx-q_{1\beta}\frac{\partial
f_{lk}}{\partial k_{\beta}}$ and $\frac{\partial f_{l,k}}{\partial
k_{\alpha}}-\frac{\partial f_{l,k+q_{1}}}{\partial k_{\alpha}}\approx-
q_{1\beta}\frac{\partial^{2}f_{lk}}{\partial k_{\alpha}\partial k_{\beta}}$
and arrive at a much simpler expression
$\frac{q_{1\beta}q_{2\gamma}}{\hbar\omega_{1}+i\hbar\Gamma}\frac{\partial^{2}f_{lk}}{\partial
k_{\gamma}\partial k_{\beta}}$. From here it is straightforward to find the GV
contribution in Eq. (3).

Case ii: In the mixed case when the intraband transition and the interband
transition coexist, the term in brackets in Eq. (9) for $l=l^{\prime\prime}$
becomes,

$\frac{(-q_{1,\beta}M_{l,l^{\prime},{\bm{k}}})(f_{l^{\prime},{\bm{k}}}-f_{l,\bm{k+q_{1}}})}{E_{l^{\prime},{\bm{k}}}-E_{l,\bm{k+q_{1}}}+\hbar(\omega_{1}+i\Gamma)}\\\
-\frac{(-q_{2,\gamma}M_{l,l^{\prime},{\bm{k}}})(f_{l,\bm{k+q_{1}}}-f_{l,\bm{k+q_{1}+q_{2}}})}{\hbar(\omega_{1}+i\Gamma)}$
(11)

and for $l^{\prime}=l^{\prime\prime}$,

$(-q_{2,\Gamma}M_{l,l^{\prime},\bm{k+q_{1}}})\frac{f_{l^{\prime},{\bm{k}}}-f_{l^{\prime},\bm{k+q_{1}}}}{\hbar(\omega_{1}+i\Gamma)}\\\
-(-q_{1,\beta}M_{l,l^{\prime},\bm{k+q_{2}}})\frac{f_{l^{\prime},\bm{k+q_{2}}}-f_{l,\bm{k+q_{1}+q_{2}}}}{E_{l^{\prime},\bm{k+q_{2}}}-E_{l,\bm{k+q_{1}+q_{2}}}+\hbar(\omega_{1}+i\Gamma)}.$
(12)

We can combine the two terms to arrive at the following expression

$\frac{q_{1\beta}q_{2\Gamma}M_{l,l^{\prime},k}}{\hbar\omega_{1}+i\hbar\Gamma}\Big{[}\frac{\partial
f_{l^{\prime},k}-f_{l,k}}{\partial k_{\beta}}\Big{]}\\\
+q_{1,\beta}q_{2,\Gamma}\frac{\partial}{\partial
k_{\Gamma}}\left[\frac{M_{l,l^{\prime},k}(f_{l^{\prime},\bm{k}}-f_{l,\bm{k}})}{E_{l^{\prime},\bm{k}}-E_{l,\bm{k}}+\hbar(\omega_{1}+i\Gamma)}\right].$
(13)

The first term gives the BD contribution in Eq. (3). The second term can be
expressed as a sum of two terms where one term contains the derivatives of the
FD functions and the other term includes the remaining contributions including
the FD functions. Note that, the term with the derivative of FD functions
leads to the BDC term in (3), and the other term becomes the second part of
Eq. (2).

Case iii: When both the virtual transitions are of interband type, the
conductivity Eq. (9) needs no further simplification and we can simply put
$q_{1}$ and $q_{2}$ to be zero and arrive at the first term of Eq. (2).

##  Appendix B Derivation of BD

The expression for BD follows from the second term in Eq. (3) and reads

${\tilde{\sigma}}_{\alpha\beta\gamma}^{(2,BD)}=\frac{e^{3}\hbar}{S}\sum\limits_{l\neq
l^{\prime}{\bm{k}}}\frac{\langle\hat{v}_{\alpha}\rangle_{l^{\prime}l}}{\left[E_{l^{\prime}l}+\hbar(\omega_{1}+\omega_{2}+i\Gamma)\right]}\\\
\times\frac{\langle\hat{v}_{\gamma}\rangle_{ll^{\prime}}\frac{\partial
f_{l^{\prime}l}}{\partial
k_{\beta}}}{E_{ll^{\prime}}\hbar(\omega_{1}+i\Gamma)}.$ (14)

In the semiclassical limit,
$E_{l^{\prime}l}\gg\hbar\omega_{1},\hbar\omega_{2}$, we can neglect the
frequencies in the denominator
$E_{l^{\prime}l}+\hbar(\omega_{1}+\omega_{2}+i\Gamma)$, and then Eq. (14)
simplifies to

${\tilde{\sigma}}_{\alpha\beta\gamma}^{(2,BD)}=\frac{e^{3}\hbar}{S\hbar(\omega_{1}+i\Gamma)}\sum\limits_{l\neq
l^{\prime}{\bm{k}}}\frac{\langle\hat{v}_{\alpha}\rangle_{l^{\prime}l}\langle\hat{v}_{\gamma}\rangle_{ll^{\prime}}}{E_{l^{\prime}l}^{2}}\frac{\partial
f_{l^{\prime}l}}{\partial k_{\beta}}.$ (15)

By using the fact that $\partial f_{l^{\prime}l}/\partial k_{\beta}=-\partial
f_{ll^{\prime}}/\partial k_{\beta}$, and then separating the terms
proportional to $\partial f_{l}/\partial k_{\beta}$ and $\partial
f_{l^{\prime}}/\partial k_{\beta}$, we obtain

${\tilde{\sigma}}_{\alpha\beta\gamma}^{(2,BD)}=\frac{e^{3}\hbar}{S\hbar(\omega_{1}+i\Gamma)}\sum\limits_{l^{\prime}{\bm{k}}}\sum\limits_{l\neq
l^{\prime}}\frac{\left(\langle\hat{v}_{\alpha}\rangle_{l^{\prime}l}\langle\hat{v}_{\gamma}\rangle_{ll^{\prime}}\right)-\left(\langle\hat{v}_{\alpha}\rangle_{l^{\prime}l}\langle\hat{v}_{\gamma}\rangle_{ll^{\prime}}\right)^{*}}{E_{l^{\prime}l}^{2}}\frac{\partial
f_{l^{\prime}}}{\partial k_{\beta}}.$ (16)

Equation (16) can be expressed as

${\tilde{\sigma}}_{\alpha\beta\gamma}^{(2,BD)}=\frac{e^{3}\hbar}{S\hbar(\omega_{1}+i\Gamma)}\sum\limits_{l^{\prime}{\bm{k}}}\Omega_{\delta,l}({\bm{k}})\frac{\partial
f_{l^{\prime}}}{\partial k_{\beta}},$ (17)

where $\delta=\alpha\times\gamma$, is the direction of Berry curvature. For a
2D system, the Berry curvature vector always points in the $z$ direction and
$\alpha$ and $\gamma$ must be linearly independent for a non-zero Berry
curvature. Using integration by parts the summation in Eq. (17) can be
separated into two parts, where the first part is proportional to the FD
functions at ${\bm{k}}\rightarrow\infty$ and the second part resembling Eq.
(4). The FD function vanishes at ${\bm{k}}\rightarrow\infty$ and Eq. (17) can
be expressed as

$D_{\alpha
z}=\sum\limits_{l{\bm{k}}}\frac{\partial\Omega_{z,l}({\bm{k}})}{\partial
k_{\alpha}}f_{lk}.$ (18)

## Acknowledgements

We acknowledge the Laboratory Directed Research and Development program of Los
Alamos National Laboratory under project number 20190574ECR. R.K.M. also
thanks the Center for Nonlinear Studies at LANL for financial support under
project 20190495CR.

## References

  * Kane and Mele [2005a] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005a), URL https://link.aps.org/doi/10.1103/PhysRevLett.95.146802.
  * Kane and Mele [2005b] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005b), URL https://link.aps.org/doi/10.1103/PhysRevLett.95.226801.
  * Bernevig et al. [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006), ISSN 0036-8075, URL https://science.sciencemag.org/content/314/5806/1757.
  * König et al. [2007] M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007), ISSN 0036-8075, URL https://science.sciencemag.org/content/318/5851/766.
  * Liu et al. [2008] C. Liu, T. L. Hughes, X.-L. Qi, K. Wang, and S.-C. Zhang, Phys. Rev. Lett. 100, 236601 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.100.236601.
  * Du et al. [2015] L. Du, I. Knez, G. Sullivan, and R.-R. Du, Phys. Rev. Lett. 114, 096802 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.096802.
  * Qian et al. [2014] X. Qian, J. Liu, L. Fu, and J. Li, Science 346, 1344 (2014), ISSN 0036-8075, URL https://science.sciencemag.org/content/346/6215/1344.
  * Tang et al. [2017] S. Tang, C. Zhang, D. Wong, Z. Pedramrazi, H.-Z. Tsai, C. Jia, B. Moritz, M. Claassen, H. Ryu, S. Kahn, et al., Nature Physics 13, 683 (2017), URL https://doi.org/10.1038/nphys4174.
  * Sheng et al. [2006] D. N. Sheng, Z. Y. Weng, L. Sheng, and F. D. M. Haldane, Phys. Rev. Lett. 97, 036808 (2006), URL https://link.aps.org/doi/10.1103/PhysRevLett.97.036808.
  * Prodan [2009] E. Prodan, Phys. Rev. B 80, 125327 (2009), URL https://link.aps.org/doi/10.1103/PhysRevB.80.125327.
  * Yang et al. [2011] Y. Yang, Z. Xu, L. Sheng, B. Wang, D. Y. Xing, and D. N. Sheng, Phys. Rev. Lett. 107, 066602 (2011), URL https://link.aps.org/doi/10.1103/PhysRevLett.107.066602.
  * Stille et al. [2012] L. Stille, C. J. Tabert, and E. J. Nicol, Phys. Rev. B 86, 195405 (2012), URL https://link.aps.org/doi/10.1103/PhysRevB.86.195405.
  * Tabert and Nicol [2013] C. J. Tabert and E. J. Nicol, Phys. Rev. Lett. 110, 197402 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.110.197402.
  * Ledwith et al. [2018] P. Ledwith, W. J. M. Kort-Kamp, and D. A. R. Dalvit, Phys. Rev. B 97, 165426 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.165426.
  * Kort-Kamp [2017] W. J. M. Kort-Kamp, Phys. Rev. Lett. 119, 147401 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.147401.
  * Rodriguez-Lopez et al. [2018] P. Rodriguez-Lopez, W. J. M. Kort-Kamp, D. A. R. Dalvit, and L. M. Woods, Phys. Rev. Materials 2, 014003 (2018), URL https://link.aps.org/doi/10.1103/PhysRevMaterials.2.014003.
  * Ezawa [2013] M. Ezawa, Phys. Rev. Lett. 110, 026603 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.110.026603.
  * Rodriguez-Lopez et al. [2017] P. Rodriguez-Lopez, W. J. M. Kort-Kamp, D. A. R. Dalvit, and L. M. Woods, Nature Communications 8, 14699 (2017), URL https://doi.org/10.1038/ncomms14699.
  * Tabert and Nicol [2014] C. J. Tabert and E. J. Nicol, Phys. Rev. B 89, 195410 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.89.195410.
  * Farias et al. [2018] M. B. Farias, W. J. M. Kort-Kamp, and D. A. R. Dalvit, Phys. Rev. B 97, 161407 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.161407.
  * Malla and Raikh [2018a] R. K. Malla and M. E. Raikh, Phys. Rev. B 98, 115437 (2018a), URL https://link.aps.org/doi/10.1103/PhysRevB.98.115437.
  * Wu et al. [2018] W. Wu, W. Zhang, S. Chen, X. Ling, W. Shu, H. Luo, S. Wen, and X. Yin, Opt. Express 26, 23705 (2018), URL http://www.opticsexpress.org/abstract.cfm?URI=oe-26-18-23705.
  * Wu et al. [2020] W. Wu, S. Chen, W. Xu, Z. Liu, R. Lou, L. Shen, H. Luo, S. Wen, and X. Yin, Photon. Res. 8, B47 (2020), URL http://www.osapublishing.org/prj/abstract.cfm?URI=prj-8-12-B47.
  * Hamamoto et al. [2017] K. Hamamoto, M. Ezawa, K. W. Kim, T. Morimoto, and N. Nagaosa, Phys. Rev. B 95, 224430 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.95.224430.
  * Moore and Orenstein [2010] J. E. Moore and J. Orenstein, Phys. Rev. Lett. 105, 026805 (2010), URL https://link.aps.org/doi/10.1103/PhysRevLett.105.026805.
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.216806.
  * You et al. [2018] J.-S. You, S. Fang, S.-Y. Xu, E. Kaxiras, and T. Low, Phys. Rev. B 98, 121109 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.121109.
  * Zhang et al. [2018] Y. Zhang, J. van den Brink, C. Felser, and B. Yan, 2D Materials 5, 044001 (2018), URL https://doi.org/10.1088/2053-1583/aad1ae.
  * Son et al. [2019] J. Son, K.-H. Kim, Y. H. Ahn, H.-W. Lee, and J. Lee, Phys. Rev. Lett. 123, 036806 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.036806.
  * Xu et al. [2018] S.-Y. Xu, Q. Ma, H. Shen, V. Fatemi, S. Wu, T.-R. Chang, G. Chang, A. M. M. Valdivia, C.-K. Chan, Q. D. Gibson, et al., Nature Physics 14, 900 (2018), URL https://doi.org/10.1038/s41567-018-0189-6.
  * Du et al. [2018] Z. Z. Du, C. M. Wang, H.-Z. Lu, and X. C. Xie, Phys. Rev. Lett. 121, 266601 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.266601.
  * Facio et al. [2018] J. I. Facio, D. Efremov, K. Koepernik, J.-S. You, I. Sodemann, and J. van den Brink, Phys. Rev. Lett. 121, 246403 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.246403.
  * Battilomo et al. [2019] R. Battilomo, N. Scopigno, and C. Ortix, Phys. Rev. Lett. 123, 196403 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.196403.
  * Lahiri et al. [2021] S. Lahiri, T. Bhore, K. Das, and A. Agarwal, _Berry curvature induced nonlinear magnetoresistivity in two dimensional systems_ (2021), eprint 2109.08979.
  * Ma et al. [2019] Q. Ma, S.-Y. Xu, H. Shen, D. MacNeill, V. Fatemi, T.-R. Chang, A. M. Mier Valdivia, S. Wu, Z. Du, C.-H. Hsu, et al., Nature 565, 337 (2019), URL https://doi.org/10.1038/s41586-018-0807-6.
  * Nandy and Sodemann [2019] S. Nandy and I. Sodemann, Phys. Rev. B 100, 195117 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.195117.
  * Du et al. [2019] Z. Z. Du, C. M. Wang, S. Li, H.-Z. Lu, and X. C. Xie, Nature Communications 10, 3047 (2019), URL https://doi.org/10.1038/s41467-019-10941-3.
  * Du et al. [2021a] Z. Z. Du, C. M. Wang, H.-P. Sun, H.-Z. Lu, and X. C. Xie, Nature Communications 12, 5038 (2021a), URL https://doi.org/10.1038/s41467-021-25273-4.
  * Xiao et al. [2019] C. Xiao, Z. Z. Du, and Q. Niu, Phys. Rev. B 100, 165422 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.165422.
  * He et al. [2019] P. He, S. S.-L. Zhang, D. Zhu, S. Shi, O. G. Heinonen, G. Vignale, and H. Yang, Phys. Rev. Lett. 123, 016801 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.016801.
  * Zhang and Vignale [2018] S. S.-L. Zhang and G. Vignale, in _Spintronics XI_ (International Society for Optics and Photonics, 2018), vol. 10732, p. 1073215\.
  * Li et al. [2021] R.-H. Li, O. G. Heinonen, A. A. Burkov, and S. S.-L. Zhang, Phys. Rev. B 103, 045105 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.045105.
  * Nandy et al. [2021] S. Nandy, C. Zeng, and S. Tewari, _Chiral anomaly induced nonlinear hall effect in multi-weyl semimetals_ (2021), eprint 2104.14969.
  * Zheng et al. [2020] S.-H. Zheng, H.-J. Duan, J.-K. Wang, J.-Y. Li, M.-X. Deng, and R.-Q. Wang, Phys. Rev. B 101, 041408 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.101.041408.
  * Burkov [2017] A. A. Burkov, Phys. Rev. B 96, 041110 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.041110.
  * Nandy et al. [2017] S. Nandy, G. Sharma, A. Taraphder, and S. Tewari, Phys. Rev. Lett. 119, 176804 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.176804.
  * Chen et al. [2018] F. C. Chen, X. Luo, J. Yan, Y. Sun, H. Y. Lv, W. J. Lu, C. Y. Xi, P. Tong, Z. G. Sheng, X. B. Zhu, et al., Phys. Rev. B 98, 041114 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.041114.
  * He et al. [2014] L. P. He, X. C. Hong, J. K. Dong, J. Pan, Z. Zhang, J. Zhang, and S. Y. Li, Phys. Rev. Lett. 113, 246402 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.113.246402.
  * Liang et al. [2015] T. Liang, Q. Gibson, M. N. Ali, M. Liu, R. J. Cava, and N. P. Ong, Nature Materials 14, 280 (2015), URL https://doi.org/10.1038/nmat4143.
  * Zhang et al. [2016] C.-L. Zhang, S.-Y. Xu, I. Belopolski, Z. Yuan, Z. Lin, B. Tong, G. Bian, N. Alidoust, C.-C. Lee, S.-M. Huang, et al., Nature Communications 7, 10735 (2016), URL https://doi.org/10.1038/ncomms10735.
  * Li et al. [2016] Q. Li, D. E. Kharzeev, C. Zhang, Y. Huang, I. Pletikosić, A. V. Fedorov, R. D. Zhong, J. A. Schneeloch, G. D. Gu, and T. Valla, Nature Physics 12, 550 (2016), URL https://doi.org/10.1038/nphys3648.
  * Xiong et al. [2015] J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan, M. Hirschberger, W. Wang, R. J. Cava, and N. P. Ong, Science 350, 413 (2015), ISSN 0036-8075, eprint https://science.sciencemag.org/content/350/6259/413.full.pdf, URL https://science.sciencemag.org/content/350/6259/413.
  * Hirschberger et al. [2016] M. Hirschberger, S. Kushwaha, Z. Wang, Q. Gibson, S. Liang, C. A. Belvin, B. A. Bernevig, R. J. Cava, and N. P. Ong, Nature Materials 15, 1161 (2016), URL https://doi.org/10.1038/nmat4684.
  * Ghiasi et al. [2019] T. S. Ghiasi, A. A. Kaverzin, P. J. Blah, and B. J. van Wees, Nano Letters 19, 5959 (2019), URL https://doi.org/10.1021/acs.nanolett.9b01611.
  * Malla and Raikh [2018b] R. K. Malla and M. E. Raikh, Phys. Rev. B 98, 144202 (2018b), URL https://link.aps.org/doi/10.1103/PhysRevB.98.144202.
  * Malla and Raikh [2018c] R. K. Malla and M. E. Raikh, Phys. Rev. B 97, 214209 (2018c), URL https://link.aps.org/doi/10.1103/PhysRevB.97.214209.
  * Malla and Raikh [2019] R. K. Malla and M. E. Raikh, Phys. Rev. B 99, 205426 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.99.205426.
  * Vogt et al. [2012] P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet, and G. Le Lay, Phys. Rev. Lett. 108, 155501 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.108.155501.
  * Dávila et al. [2014] M. E. Dávila, L. Xian, S. Cahangirov, A. Rubio, and G. L. Lay, New Journal of Physics 16, 095002 (2014), URL https://doi.org/10.1088/1367-2630/16/9/095002.
  * Zhu et al. [2015] F.-f. Zhu, W.-j. Chen, Y. Xu, C.-l. Gao, D.-d. Guan, C.-h. Liu, D. Qian, S.-C. Zhang, and J.-f. Jia, Nature Materials 14, 1020 (2015), URL https://doi.org/10.1038/nmat4384.
  * Yuhara et al. [2019] J. Yuhara, B. He, N. Matsunami, M. Nakatake, and G. Le Lay, Advanced Materials 31, 1901017 (2019), eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.201901017, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.201901017.
  * Li et al. [2013] X. Li, T. Cao, Q. Niu, J. Shi, and J. Feng, Proceedings of the National Academy of Sciences 110, 3738 (2013), ISSN 0027-8424, URL https://www.pnas.org/content/110/10/3738.
  * Liang et al. [2013] Q.-F. Liang, L.-H. Wu, and X. Hu, New Journal of Physics 15, 063031 (2013), URL https://doi.org/10.1088/1367-2630/15/6/063031.
  * Marrazzo et al. [2018] A. Marrazzo, M. Gibertini, D. Campi, N. Mounet, and N. Marzari, Phys. Rev. Lett. 120, 117701 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.117701.
  * Marrazzo et al. [2019] A. Marrazzo, M. Gibertini, D. Campi, N. Mounet, and N. Marzari, Nano Letters 19, 8431 (2019), URL https://doi.org/10.1021/acs.nanolett.9b02689.
  * Marrazzo et al. [2020] A. Marrazzo, N. Marzari, and M. Gibertini, Phys. Rev. Research 2, 012063 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.012063.
  * Cucchi et al. [2020] I. Cucchi, A. Marrazzo, E. Cappelli, S. Riccò, F. Y. Bruno, S. Lisi, M. Hoesch, T. K. Kim, C. Cacho, C. Besnard, et al., Phys. Rev. Lett. 124, 106402 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.106402.
  * Gao et al. [2014] Y. Gao, S. A. Yang, and Q. Niu, Phys. Rev. Lett. 112, 166601 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.112.166601.
  * Watanabe and Yanase [2021] H. Watanabe and Y. Yanase, Phys. Rev. X 11, 011001 (2021), URL https://link.aps.org/doi/10.1103/PhysRevX.11.011001.
  * Ezawa [2015] M. Ezawa, Journal of the Physical Society of Japan 84, 121003 (2015), eprint https://doi.org/10.7566/JPSJ.84.121003, URL https://doi.org/10.7566/JPSJ.84.121003.
  * Bychkov and Rashba [1984] Y. A. Bychkov and É. I. Rashba, Soviet Journal of Experimental and Theoretical Physics Letters 39, 78 (1984).
  * Tse et al. [2011] W.-K. Tse, Z. Qiao, Y. Yao, A. H. MacDonald, and Q. Niu, Phys. Rev. B 83, 155447 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.83.155447.
  * Shikin et al. [2008] A. M. Shikin, A. Varykhalov, G. V. Prudnikova, D. Usachov, V. K. Adamchuk, Y. Yamada, J. D. Riley, and O. Rader, Phys. Rev. Lett. 100, 057601 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.100.057601.
  * Weeks et al. [2011] C. Weeks, J. Hu, J. Alicea, M. Franz, and R. Wu, Phys. Rev. X 1, 021001 (2011), URL https://link.aps.org/doi/10.1103/PhysRevX.1.021001.
  * Rashba [2009] E. I. Rashba, Phys. Rev. B 79, 161409 (2009), URL https://link.aps.org/doi/10.1103/PhysRevB.79.161409.
  * Du et al. [2021b] Z. Z. Du, H.-Z. Lu, and X. C. Xie, Nature Reviews Physics (2021b), URL https://doi.org/10.1038/s42254-021-00359-6.
  * [77] C. Ortix, Advanced Quantum Technologies n/a, 2100056 (????), eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/qute.202100056, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/qute.202100056.
  * Malla and Kort-Kamp [2021] R. K. Malla and W. J. M. Kort-Kamp, Scientific Reports 11, 9734 (2021), URL https://doi.org/10.1038/s41598-021-89219-y.
  * Margulis et al. [2013] V. A. Margulis, E. Muryumin, and E. Gaiduk, Journal of Physics: Condensed Matter 25, 195302 (2013).
  * Golub and Tarasenko [2014] L. E. Golub and S. A. Tarasenko, Phys. Rev. B 90, 201402 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.90.201402.
  * Mikhailov [2016] S. A. Mikhailov, Phys. Rev. B 93, 085403 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.085403.