High-frequency rectification via chiral Bloch electrons

TL;DR

This paper introduces a novel rectification mechanism based on chiral Bloch electrons in quantum crystals, enabling high-frequency rectifiers without semiconductor junctions, with potential for advanced electronic applications.

Contribution

It uncovers a new rectification mechanism from skew scattering in chiral quantum materials, expanding the design possibilities for high-frequency electronic devices.

Findings

Large, tunable rectification effects in graphene multilayers.

Identification of skew scattering as a rectification mechanism.

Potential for high-frequency rectifiers through material engineering.

Abstract

Rectification is a process that converts electromagnetic fields into a direct current. Such a process underlies a wide range of technologies such as wireless communication, wireless charging, energy harvesting, and infrared detection. Existing rectifiers are mostly based on semiconductor diodes, with limited applicability to small voltage or high frequency inputs. Here, we present an alternative approach to current rectification that uses the intrinsic electronic properties of quantum crystals without using semiconductor junctions. We identify a previously unknown mechanism for rectification from skew scattering due to the inherent chirality of itinerant electrons in time-reversal-invariant but inversion-breaking materials. Our calculations reveal large, tunable rectification effects in graphene multilayers and transition metal dichalcogenides. Our work demonstrates the possibility of realizing high-frequency rectifiers by rational material design and quantum wavefunction engineering.

Equations

$$j_a={\chi }_{abc}{E}_b^{\ast }{E}_c,$$

$$\frac{\partial f}{\partial t}+\dot{\mathbf{k}}\cdot \frac{\partial f}{\partial \mathbf{k}}=-\mathcal{C}[f],$$

$$w_{{\mathbf{k}}^{\prime }\mathbf{k}}^{(A)}=-(2\pi {)}^2{\int }_{\mathbf{q}}\mathrm{Im}⟨V_{{\mathbf{k}}^{\prime }\mathbf{q}}{V}_{\mathbf{q}\mathbf{k}}{V}_{\mathbf{k}{\mathbf{k}}^{\prime }}⟩\delta ({ϵ}_{\mathbf{k}}-{ϵ}_{{\mathbf{k}}^{\prime }})\delta ({ϵ}_{{\mathbf{k}}^{\prime }}-{ϵ}_{\mathbf{q}}),$$

$${\chi }_{abc}=$$

$$\displaystyle+\tau^{\prime}_{2}\tau_{2}\tau_{1}(\omega)\int_{\bm{k}}v_{0,a}\int_{\bm{k}^{\prime}}w_{\bm{k}\bm{k}^{\prime}}^{(A)}\partial_{k^{\prime}_{b}}\partial_{k^{\prime}_{c}}f_{0}(\bm{k}^{\prime})\bigg{]}$$

$$+(b\leftrightarrow c{)}^{\ast },$$

$$j_2\sim en{v}_F\cdot {\left(\frac{eE\Delta t}{p_F}\right)}^2\cdot \frac{\tau }{\tilde{\tau }}$$

$$\begin{array}{c}{\mathfrak{R}}_I=\frac{1}{L}\frac{\chi }{{\sigma }_{\omega }}=\frac{{\eta }_I}{L}\frac{\tau }{\tilde{\tau }},{\eta }_I\sim \frac{e\tau }{p_F},\\ {\mathfrak{R}}_V=\frac{1}{L}\frac{\chi }{{\sigma }_{\omega }{\sigma }_0}=\frac{{\eta }_V}{L}\frac{\tau }{\tilde{\tau }},{\eta }_V\sim \frac{1}{en{v}_F},\end{array}$$

$${\eta }_P=\frac{{\chi }^2{E}^2}{{\sigma }_0{\sigma }_{\omega }}\propto {\left(\frac{eE\Delta t}{p_F}\right)}^2\cdot {\left(\frac{\tau }{\tilde{\tau }}\right)}^2.$$

$$H_s(\mathbf{k})$$

$$\mathcal{C}[f]={\int }_{{\mathbf{k}}^{\prime }}\left[w_{{\mathbf{k}}^{\prime }\mathbf{k}}f(\mathbf{k},{ϵ}_{\mathbf{k}})-w_{\mathbf{k}{\mathbf{k}}^{\prime }}f({\mathbf{k}}^{\prime },{ϵ}_{\mathbf{k}}+\delta {ϵ}_{{\mathbf{k}}^{\prime }\mathbf{k}})\right].$$

$$f^{\text{scatt}}$$

$$f^{\text{adist}}$$

$$\mathbf{v}$$

$$H=H_0+V+U,$$

$$w_{{\mathbf{k}}^{\prime }\mathbf{k}}=2\pi |T_{{\mathbf{k}}^{\prime }\mathbf{k}}|\delta ({ϵ}_{{\mathbf{k}}^{\prime }}-{ϵ}_{\mathbf{k}}),$$

$$\mathcal{P}:{\Omega }_a(\mathbf{k})\to {\Omega }_a(-\mathbf{k}),\mathcal{T}:{\Omega }_a(\mathbf{k})\to -{\Omega }_a(-\mathbf{k}).$$

$$\delta {ϵ}_{{\mathbf{k}}^{\prime }\mathbf{k}}=-\mathbf{F}\cdot \delta {\mathbf{r}}_{{\mathbf{k}}^{\prime }\mathbf{k}}.$$

$$\frac{\partial f}{\partial t}+\dot{\mathbf{k}}\cdot \frac{\partial f}{\partial \mathbf{k}}=-\mathcal{C}[f].$$

$$\mathcal{C}[f]={\int }_{{\mathbf{k}}^{\prime }}[w_{{\mathbf{k}}^{\prime }\mathbf{k}}f(\mathbf{k},ϵ)-w_{\mathbf{k}{\mathbf{k}}^{\prime }}f({\mathbf{k}}^{\prime },ϵ+\delta {ϵ}_{{\mathbf{k}}^{\prime }\mathbf{k}})],$$

$$w_{{\mathbf{k}}^{\prime }\mathbf{k}}=2\pi |T_{{\mathbf{k}}^{\prime }\mathbf{k}}|\delta ({ϵ}_{{\mathbf{k}}^{\prime }}-{ϵ}_{\mathbf{k}}),$$

$$w_{\mathbf{k}{\mathbf{k}}^{\prime }}=w_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(S)}+w_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(A)},$$

$$w_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(S)}=\frac{1}{2}(w_{\mathbf{k}{\mathbf{k}}^{\prime }}+w_{{\mathbf{k}}^{\prime }\mathbf{k}}),w_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(A)}=\frac{1}{2}(w_{\mathbf{k}{\mathbf{k}}^{\prime }}-w_{{\mathbf{k}}^{\prime }\mathbf{k}}).$$

$$w_{{\mathbf{k}}^{\prime }\mathbf{k}}^{(S)}=w_{-\mathbf{k},-{\mathbf{k}}^{\prime }}^{(S)}=w_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(S)},$$

$$w_{{\mathbf{k}}^{\prime }\mathbf{k}}^{(A)}=-w_{-\mathbf{k},-{\mathbf{k}}^{\prime }}^{(A)}=-w_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(A)}.$$

$${\int }_{{\mathbf{k}}^{\prime }}{w}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(A)}=0.$$

$$w_{{\mathbf{k}}^{\prime }\mathbf{k}}^{(S)}=2\pi ⟨|V_{{\mathbf{k}}^{\prime }\mathbf{k}}{|}^2⟩\delta ({ϵ}_{{\mathbf{k}}^{\prime }}-{ϵ}_{\mathbf{k}}),$$

$$w_{{\mathbf{k}}^{\prime }\mathbf{k}}^{(A)}=-(2\pi {)}^2{\int }_{\mathbf{q}}\mathrm{Im}⟨V_{{\mathbf{k}}^{\prime }\mathbf{q}}{V}_{\mathbf{q}\mathbf{k}}{V}_{\mathbf{k}{\mathbf{k}}^{\prime }}⟩\delta ({ϵ}_{\mathbf{k}}-{ϵ}_{{\mathbf{k}}^{\prime }})\delta ({ϵ}_{{\mathbf{k}}^{\prime }}-{ϵ}_{\mathbf{q}}),$$

$$f=f_0+f^{\text{scatt}}+f^{\text{adist}}.$$

$$f^{\text{scatt}}(\mathbf{k})=\sum _{n\ge 1}{f}_n^{\text{scatt}}(\mathbf{k}),f_n^{\text{scatt}}(\mathbf{k})=\sum _{m\ge 0}{f}_n^{(m)}(\mathbf{k}),$$

$$f^{\text{adist}}(\mathbf{k})=\sum _{n\ge 1}{f}_n^{\text{adist}}(\mathbf{k}),f_n^{\text{adist}}(\mathbf{k})=\sum _{m\ge 0}{g}_n^{(m)}(\mathbf{k}).$$

$$\frac{\partial f_1^{(0)}}{\partial t}-e\mathbf{E}\cdot \frac{\partial f_0}{\partial \mathbf{k}}=-{\int }_{{\mathbf{k}}^{\prime }}{w}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(S)}[f_1^{(0)}(\mathbf{k})-f_1^{(0)}({\mathbf{k}}^{\prime })],$$

$$\frac{\partial f_1^{(1)}}{\partial t}=-{\int }_{{\mathbf{k}}^{\prime }}{w}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(S)}[f_1^{(1)}(\mathbf{k})-f_1^{(1)}({\mathbf{k}}^{\prime })]+{\int }_{{\mathbf{k}}^{\prime }}{w}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(A)}{f}_1^{(0)}({\mathbf{k}}^{\prime }),$$

$$\frac{\partial f_2^{(0)}}{\partial t}-e\mathbf{E}\cdot \frac{\partial f_1^{(0)}}{\partial \mathbf{k}}=-{\int }_{{\mathbf{k}}^{\prime }}{w}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(S)}[f_2^{(0)}(\mathbf{k})-f_2^{(0)}({\mathbf{k}}^{\prime })],$$

Full text

Introduction

An important goal of basic research on quantum materials is to breed new quantum technologies that can address the increasingly complex energy challenges. The past decade has seen a dramatic change in the model for energy production, consumption, and transportation. Due to the explosive growth of wireless technologies and portable devices, there is now increasing effort towards developing microscaled devices that are able to harvest ambient energy into usable electrical energy. The physical process central to harvesting electromagnetic energy is rectification, which refers to the conversion from an oscillating electromagnetic field to a DC current. Existing rectifiers operating at radio frequency are mostly based on electrical circuits with diodes, where the built-in electric field in the semiconductor junction sets the direction of the DC current. Such diodes face two fundamental limitations rectifier ; energy . First, rectification requires a threshold input voltage $V_{T}=k_{B}T/e$, known as thermal voltage (about $26\,\text{mV}$ at room temperature). Second, the responsivity is limited by the transition time in diodes (typically order of nanoseconds) and drops at high frequencies. On the other hand, because of the fast-developing microwatts and nanowatts electronics and next generation wireless networks, energy harvesters of electromagnetic field in the microwave and terahertz (THz) frequency range are in great demand.

High-frequency rectifiers can also be used in sensor and detector technology for the infrared, far infrared and submillimeter bands detector1 ; detector2 , which has wide-ranging applications in medicine, biology, climatology, meteorology, telecommunication, astronomy, etc. However, there is a so-called terahertz gap (0.1 to 10 THz) between the operating frequencies of electrical diodes and photodiodes. At frequencies within this range, efficient detection technology remains to be developed.

Instead of using semiconductor junctions, rectification can be realized as the nonlinear electrical or optical response of noncentrosymmetric crystals (Fig. 1). In particular, the second-order nonlinearity $\chi(\omega)$ is an intrinsic material property that characterizes the DC current generated by an external electric field oscillating at frequency $\omega$. Rectification in a single homogenous material is not limited by the thermal voltage threshold or the transition time innate to a semiconductor junction. Moreover, nonlinear electrical or optical response of metals and degenerate semiconductors is much faster than photothermal effects used in bolometers for thermal radiation detection.

However, second-order nonlinearity of most materials is small. While nonlinear optical properties of quantum materials Boyd ; Chemla such as photocurrent and second-harmonic generation are extensively studied Weyl2 ; Weyl3 , much less is known about second-order response at radio, microwave, and infrared frequencies Genkin ; Soref ; Leuthold . In this direction, recent works have predicted intraband photocurrent Moore ; Spivak and second-order nonlinear Hall effect due to “Berry curvature dipole” Sodemann in nonmagnetic materials at zero magnetic field. In particular, the nonlinear Hall effect is predicted to be prominent in materials with titled Dirac or Weyl cones, which are a source of large Berry curvature dipole Sodemann . Very recently, this effect was observed for the first time in low-frequency (around 100 Hz) transport measurements on the two-dimensional (2D) transition metal dichalcogenide (TMD) 1T*′*-WTe2 dipole ; Mak . The second-order Hall conductivity of bilayer WTe2 is remarkably large dipole , in agreement with its large Berry curvature dipole from the tilted Dirac dispersion Qian . This and other types of nonlinear response from intraband processes are also being explored in topological insulator surface states, Weyl semimetals, Rashba systems and heavy fermion materials Ganichev ; Morimoto ; Ishizuka ; Koenig ; nonreciprocal ; Zhang ; Rostami ; Facio ; Dzsaber . Despite the recent progress, mechanisms for nonlinear electric and terahertz response of noncentrosymmetric crystals remain to be thoroughly studied.

In this work, we present a systematic theoretical study of intraband second-order response in time-reversal-invariant noncentrosymmetric materials, using a semiclassical Boltzmann equation and taking account of electron Bloch wavefunctions via quantum-mechanical scattering rates. We identify a new contribution to rectification from skew scattering with nonmagnetic impurities. Importantly, such skew scattering arises from the intrinsic chirality of Bloch wavefunctions in momentum space and does not require spin-orbit coupling. We show that the contributions to rectification from skew scattering and from Berry curvature dipole scale differently with the impurity concentration, and skew scattering predominates at low impurity concentration. Moreover, skew scattering is allowed in all noncentrosymmetric crystals, whereas Berry curvature dipole requires more strict symmetry conditions. Based on this new mechanism, we predict large and highly tunable rectification effects in graphene multilayers and heterostructures, as well as 2H-TMD monolayers.

Results

Second-order response

We study the DC current in a homogeneous material generated at second order by an external electric field $\bm{E}$ of frequency $\omega$. The second-order response also involves a $2\omega$ component, corresponding to the second-harmonic generation, which we do not focus on in this work. We write down the second-order DC response as

[TABLE]

where $\chi_{abc}$ is the rank-three tensor for the second-order conductivity, $E_{a}=E_{a}(\omega)$ is the (complex) amplitude of the external electric field of frequency $\omega$, and the indices $a,b,c$ label the coordinates. $\chi_{abc}$ satisfies $\chi_{abc}=\chi_{acb}^{*}$. Since the current $j_{a}$ is odd under inversion while the electric field $E_{a}$ is even, finite $\chi_{abc}$ is possible only in noncentrosymmetric media. The second-order response can be decomposed into two parts as $j_{a}=\chi_{abc}^{\prime}E_{b}^{*}E_{c}+i\chi_{abc}^{\prime\prime}E_{b}^{*}E_{c}$, with $\chi_{abc}^{\prime}=\chi_{acb}^{\prime}$ and $\chi_{abc}^{\prime\prime}=-\chi_{acb}^{\prime\prime}$. The former describes the response to a linearly-polarized field and the latter to a circularly-polarized field Belinicher2 . It is important to note that in an isotropic medium second-order response to a linearly-polarized field must vanish. In other words, a nonzero $\chi^{\prime}_{abc}$ requires the presence of crystal anisotropy.

Nonlinear responses are extensively studied in the optical frequency regime Boyd . The classical approach to nonlinear optics considers electrons bound to a nucleus by an anharmonic potential, so that the restoring force to the displacement of electrons becomes nonlinear. In the quantum-mechanical theory, the nonlinear optical response at frequencies larger than the band gap is usually dominated by electric dipole transitions between different bands. On the other hand, with energy harvesting and infrared detection in mind, in this work we consider intraband nonlinear response at frequencies below the interband transition threshold.

Semiclassical Boltzmann analysis

We analyze the second-order electrical transport using the semiclassical Boltzmann theory Belinicher2 . For a homogeneous system, the Boltzmann equation is given by $(\hbar=1)$

[TABLE]

where $f(\bm{k},\epsilon)$ is the distribution function at energy $\epsilon$ and $\mathcal{C}[f]$ denotes the collision integral. The time derivative of the wavevector is equal to the force felt by an electron $\dot{\bm{k}}=-e\bm{E}$ under the external electric field $\bm{E}$. The distribution function $f$ deviates from the equilibrium Fermi–Dirac distribution when electrons are accelerated by the external field. A nonequilibrium steady state is obtained by the balance of the acceleration and the relaxation due to scattering, described by the collision integral $\mathcal{C}[f]$, whose explicit form is shown later in a general form. It includes the scattering rate $w_{\bm{k}^{\prime}\bm{k}}$, the probability of the transition per unit time from a Bloch state with a wavevector $\bm{k}$ to that with $\bm{k}^{\prime}$. We will consider scattering by impurities. The impurity density should be sufficiently lower than the electron density for the semiclassical Boltzmann analysis to be valid; see Materials and Methods for further discussions.

To obtain a nonzero second-order conductivity $\chi_{abc}$, we need to capture the effect of inversion breaking. For time-reversal-invariant systems, the energy dispersion $\epsilon_{\bm{k}}$ and the band velocity $\bm{v}_{0}(\bm{k})=\nabla_{\bm{k}}\epsilon_{\bm{k}}$ satisfy the conditions $\epsilon_{\bm{k}}=\epsilon_{-\bm{k}}$ and $\bm{v}_{0}(\bm{k})=-\bm{v}_{0}(-\bm{k})$, which hold with or without inversion symmetry: that is, the absence of inversion is encoded in the wavefunction but not in the energy dispersion.

One consequence of inversion asymmetry is manifested in a scattering process, as the transition rate depends not only on the dispersion but also on the Bloch wavefunction $\psi_{\bm{k}}$. For and only for noncentrosymmetric crystals, the transition rate from $\bm{k}$ to $\bm{k}^{\prime}$ and from $-\bm{k}$ to $-\bm{k}^{\prime}$ can be different: $w_{\bm{k}^{\prime}\bm{k}}\neq w_{-\bm{k}^{\prime},-\bm{k}}$. A scattering process with such asymmetry is referred to as skew scattering. We will see that skew scattering is a source of the second-order DC response, i.e., rectification. The strength of skew scattering is characterized by $w_{\bm{k}^{\prime}\bm{k}}^{(A)}=(w_{\bm{k}^{\prime}\bm{k}}-w_{-\bm{k}^{\prime},-\bm{k}})/2$, whereas the symmetric component $w_{\bm{k}^{\prime}\bm{k}}^{(S)}=(w_{\bm{k}^{\prime}\bm{k}}+w_{-\bm{k}^{\prime},-\bm{k}})/2$ does not rely on inversion breaking. The presence of time-reversal symmetry guarantees $w_{\bm{k}^{\prime}\bm{k}}=w_{-\bm{k},-\bm{k}^{\prime}}$, which allows us to write $w_{\bm{k}^{\prime}\bm{k}}^{(S,A)}=(w_{\bm{k}^{\prime}\bm{k}}\pm w_{\bm{k}\bm{k}^{\prime}})/2$.

The transition rate $w_{\bm{k}^{\prime}\bm{k}}$ is usually calculated to the lowest order in scattering potential by Fermi’s golden rule. For elastic impurity scattering, it is given by $w_{\bm{k}^{\prime}\bm{k}}^{(S)}=2\pi\langle|V_{\bm{k}^{\prime}\bm{k}}|^{2}\rangle\delta(\epsilon_{\bm{k}^{\prime}}-\epsilon_{\bm{k}})$, where $V_{\bm{k}^{\prime}\bm{k}}$ is the matrix element of a single scatterer and $\langle\ \rangle$ denotes the impurity average. However, this formula is symmetric under the exchange $\bm{k}\leftrightarrow\bm{k}^{\prime}$ and does not capture skew scattering. The latter arises at the next-leading order from the interference of a direct transition $\bm{k}\to\bm{k}^{\prime}$ and an indirect process $\bm{k}\to\bm{k}^{\prime}$ via an intermediate state $\bm{q}$ Luttinger1 :

[TABLE]

which indeed satisfies $w_{\bm{k}^{\prime}\bm{k}}^{(A)}=-w_{\bm{k}\bm{k}^{\prime}}^{(A)}$. Here we use the notation $\int_{\bm{q}}=\int d^{d}q/(2\pi)^{d}$ ($d$ is the spatial dimension).

It is well known that skew scattering provides an extrinsic contribution to the anomalous Hall effect AHE and the spin Hall effect SHE . These are linear response phenomena in systems with broken time-reversal and spin-rotational symmetries, respectively. We emphasize that skew scattering exists and contributes to second-order response even in materials with both time-reversal and spin-rotational symmetry. With regard to its microscopic origin, skew scattering can be caused by nonmagnetic impurities when the single-impurity potential $V(\bm{r})$ is inversion asymmetric, such as a dipole or octupole potential Belinicher2 . Here we show that even for a symmetric single-impurity potential $V(\bm{r})=V(-\bm{r})$ such as a delta function potential, skew scattering can also appear. In such case, skew scattering is due to the inherent chirality of the electron wavefunction in a noncentrosymmetric crystal. In the simplest case of a delta function impurity potential, we have $V_{\bm{k}^{\prime}\bm{k}}\propto\langle u_{\bm{k}^{\prime}}|u_{\bm{k}}\rangle$ where $|u_{\bm{k}}\rangle$ is the Bloch wavefunction within a unit cell. Then the expression for $w^{(A)}$ given in Eq. (3) involves the wavefunction overlap at different $\bm{k}$ points on the Fermi surface, similar to the Berry phase formula Vanderbilt , and also skew scattering is known to be related to Berry curvature in a certain limit Nagaosa . Importantly, while there is no net chirality in nonmagnetic systems, the Bloch wavefunction $|u_{\bm{k}}\rangle$ of an electron at finite momentum is allowed to be complex and carry $\bm{k}$-dependent orbital magnetization, spin polarization, or Berry curvature. In other words, mobile electrons in inversion-breaking nonmagnetic crystals can exhibit momentum-dependent chirality that is opposite at $\bm{k}$ and $-\bm{k}$. Then, skew scattering occurs due to electron chirality in a similar way as the classical Magnus effect, where a spinning object is deflected when moving in a viscous medium. The spinning motion is associated with the self-rotation of a wave packet Chang and the viscosity results from scattering with impurities; see Fig. 2A.

To obtain the second-order conductivity, the distribution function $f(\bm{k})$ should be calculated to second order in the external electric field $\bm{E}(\omega)$. Moreover, since skew scattering $w^{(A)}$ is parametrically smaller than $w^{(S)}$, we expand the distribution function up to first order in $w^{(A)}$. By solving the Boltzmann equation self-consistently at each order of $\bm{E}$ and $w^{(A)}$, we obtain the second-order DC current response $\bm{j}=-e\int_{\bm{k}}\bm{v}_{0}(\bm{k})f(\bm{k})$ due to skew scattering in time-reversal-invariant systems. The second-order conductivity $\chi_{abc}$ can be formally expressed in terms of various scattering times $\tau_{n},\tau^{\prime}_{n}$:

[TABLE]

where $\tau_{n}(\omega)$ is a shorthand for $\tau_{n}(\omega)=\tau_{n}/(1-i\omega\tau_{n})$ and likewise for $\tau^{\prime}_{n}(\omega)$. The scattering times $\tau_{n},\tau^{\prime}_{n}$ with $n=1,2$ are associated with the dominant symmetric scattering $w^{(S)}$. They determine the steady-state distribution function up to second order in the electric field, when skew scattering is neglected. Detailed descriptions of a self-consistent solution of the Boltzmann equation and discussions about Joule heating can be found in Supplementary Materials.

A rough order-of-magnitude estimate of $\chi$ is obtained from Eq. (Semiclassical Boltzmann analysis) by using two scattering times $\tau$ and $\tilde{\tau}$ for the symmetric and skew scattering, respectively, leading to the second-order DC current of a metal or degenerate semiconductor under a linearly-polarized electric field:

[TABLE]

with $\Delta t=\tau$ for $\omega\tau\ll 1$ and $\Delta t=1/\omega$ for $\omega\tau\gg 1$. Here $v_{F}$ is the Fermi velocity, $p_{F}=\hbar k_{F}$ is the Fermi momentum, and $n$ is the electron density.

This estimate provides a heuristic understanding of the second-order response. In the low-frequency limit $\omega\rightarrow 0$, the ratio of the second-order current $j_{2}$ and the linear response current $j_{1}=\sigma_{\omega}E$ is a product of two dimensionless quantities $eE\tau/p_{F}$ and $\tau/\tilde{\tau}$. The first term is the change of electron’s momentum under the external field during the mean free time, divided by the Fermi momentum. The second term $\tau/\tilde{\tau}$ characterizes the strength of skew scattering which is responsible for second-order electrical response, relative to the symmetric scattering.

The short-circuit current responsivity $\mathfrak{R}_{I}$ is a metric for a rectifier, defined as the ratio of the generated DC current to the power dissipation. With a sample dimension $L^{2}$, the current responsivity is given by $\mathfrak{R}_{I}\equiv j_{2}L/(j_{1}EL^{2})$. Similarly, the voltage responsivity $\mathfrak{R}_{V}$ is defined for the generated DC voltage as $\mathfrak{R}_{V}\equiv(j_{2}L/\sigma_{0})/(j_{1}EL^{2})$. Both $\mathfrak{R}_{I}$ and $\mathfrak{R}_{V}$ are independent of the magnitude of the external electric field. Using the linear response conductivity $\sigma_{\omega}\sim env_{F}\cdot(e\lambda_{F}\tau/\hbar)\operatorname{Re}[(1-i\omega\tau)^{-1}]$ obtained from the same Boltzmann equation, we have

[TABLE]

Here we define the reduced current and voltage responsivities $\eta_{I}$ and $\eta_{V}$, respectively, which are independent of the sample size and the ratio $\tau/\tilde{\tau}$. The approximate relations for $\eta_{I}$ and $\eta_{V}$ hold in both low- and high-frequency limits. We also consider the ratio of the generated DC power and input power $\eta_{P}\equiv(j_{2}^{2}/\sigma_{0})/(\sigma_{\omega}E^{2})$, which characterizes the power conversion efficiency (assuming the load resistance and internal resistance are comparable). We find

[TABLE]

From these figure of merits, it is clear that a material with low carrier density $n$ or small Fermi momentum $p_{F}$ is desirable for efficient rectification. Moreover, to achieve high current responsivity a long mean free time is preferred. In this respect, high-mobility semiconductors and semimetals are promising as rectifiers for infrared detection. In the following, we analyze graphene systems and estimate their efficiency of rectification.

Graphene multilayers and transition metal dichalcogenides

Pristine graphene has high mobility and low carrier density, but it is centrosymmetric. Nevertheless, as a van der Waals material, we can easily assemble multilayer stacks to break inversion. Realizations of noncentrosymmetric structures include monolayer graphene on a hexagonal boron nitride (hBN) substrate Louie ; LeRoy ; Pablo1 ; Dean1 ; Wallece , electrically-biased bilayer graphene McCann ; graphene_review1 ; graphene_review2 , and multilayer graphene such as ABA-stacked trilayer graphene Wu ; Young . Anisotropy of energy dispersions, also required for second-order response, naturally arises from trigonal crystal structures.

We now show how skew scattering arises from the chirality of quantum wavefunctions and contributes to the second-order conductivity in graphene heterostructures with gaps induced by inversion symmetry breaking. The $k\cdot p$ Hamiltonian we consider here is

[TABLE]

where $s=\pm 1$ denotes two valleys at $K$ and $K^{\prime}$, respectively, and we define $k_{\pm}=k_{x}\pm ik_{y}$. This Hamiltonian has a band gap of $2\Delta$ in the energy spectrum and describes two graphene systems: (i) monolayer graphene on hBN, where the spatially varying atomic registries break the carbon sublattice symmetry and opens up gaps at Dirac points Louie ; LeRoy ; Pablo1 ; Dean1 ; Wallece and (ii) bilayer graphene with an out-of-plane electric field that breaks the layer symmetry and opens up gaps McCann ; graphene_review1 ; graphene_review2 . The difference between the monolayer and bilayer cases lies in the relative strength of $v$ and $\lambda$. For the monolayer (bilayer) case, $\lambda k_{F}^{2}\ll vk_{F}$ ($vk_{F}\ll\lambda k_{F}^{2}$) is responsible for the trigonal warping of the linear (quadratic) energy dispersion.

The Bloch wavefunction of $H_{s}(\bm{k})$ has two components associated with sublattice (layer) degrees of freedom in monolayer (bilayer) graphene. When the band gap is present ($\Delta\neq 0$), the wavefunction in each valley is chiral and carries finite Berry curvature, leading to skew scattering from nonmagnetic impurities. The chirality and Berry curvature are opposite for valley $K$ and $K^{\prime}$ due to time reversal symmetry. Note that because of the three-fold rotation symmetry of graphene, Berry curvature dipole vanishes Sodemann , leaving skew scattering as the only mechanism for rectification.

We calculate the second-order conductivity from Eq. (Semiclassical Boltzmann analysis), and find nonvanishing elements $\chi_{xxy}=\chi_{xyx}=\chi_{yxx}=-\chi_{yyy}\equiv\chi$, consistent with the point group symmetry. The relation among the electric field, induced current, and underlying crystalline lattice is depicted in Fig. 2B, along with the schematic pictures of the Fermi surface displacement (Fig. 2, C and D). The function $\chi$ has the form $\chi=e^{3}v_{F}(\tau^{3}/\tilde{\tau})\zeta(\epsilon_{F},\omega)$, with the dimensionless function $\zeta$ given in the Supplementary Materials. The ratio $\tau/\tilde{\tau}$ is proportional to the gap $\Delta$ and the Fermi surface trigonal warping. The former describes the effect of inversion breaking on the electronic structure, while the latter is responsible for Fermi surface asymmetry within a valley $E_{s}(\bm{k})\neq E_{s}(-\bm{k})$. Both ingredients are necessary for finite second-order response.

The second-order conductivity $\chi$ and the reduced voltage responsivity $\eta_{V}$ for both monolayer and bilayer gapped graphene are shown in Fig. 3 as the functions of frequency $\nu=\omega/(2\pi)$. The frequency dependence of the second-order conductivity is qualitatively understood with Eq. (5). We can further simplify the relation in two dimensions for $n\lambda_{F}^{2}=\text{const.}$ The second-order conductivity then behaves as $\chi\sim e^{3}v_{F}\tau^{2}/\hbar^{2}\cdot(\tau/\tilde{\tau})$ for $\omega\tau\ll 1$ and $\chi\sim e^{3}v_{F}/(\hbar\omega)^{2}\cdot(\tau/\tilde{\tau})$ for $\omega\tau\gg 1$.

The responsivity is affected by the linear conductivity $\sigma$ as well as the second-order conductivity $\chi$. At high frequencies with the carrier density fixed, we observe a decrease of the rectified current while the energy dissipation by the linear response also decreases. Since both $\sigma$ and $\chi$ decrease as $\omega^{-2}$, the responsivity saturates in the high frequency limit. With a sample size of $10\,\mu\text{m}$, the carrier density $n=0.5\times 10^{12}\,\text{cm}^{-2}$, and the ratio $\tau/\tilde{\tau}=0.01$, from reasonable estimates of mean free time and impurity concentration we find the saturated voltage responsivity $\mathfrak{R}_{V}\sim 30\,\text{V/W}$ for monolayer and $40\,\text{V/W}$ for bilayer (see also the Supplementary Materials for the frequency dependence). As expected from the qualitative estimate Eq. (6), high responsivity is realized with low carrier density. We note that the semiclassical approach is valid when the mean free path is much longer than the Fermi wavelength.

Both monolayer and bilayer graphene show broadband response. The saturation of responsivities lasts up to the onset frequency of an interband transition (not included in our study). For example, at the lowest density of $n=0.5\times 10^{12}\,\text{cm}^{-2}$ shown in Fig. 3, the Fermi energy for the monolayer case is $\epsilon_{F}\approx 80\,\text{meV}$, so that due to the Pauli blocking the interband transition does not occur at frequencies below $2\epsilon_{F}\sim 40\,\text{THz}$. The interband threshold frequency also increases with the band gap, which is tunable in bilayer graphene by the displacement field.

Finite temperature affects the intraband second-order response and responsivity through thermal smearing of the distribution function and the change of scattering times. Here thermal smearing does not change our result much as the Fermi energy is much higher than room temperature. Since it utilizes the material’s intrinsic nonlinearity, rectification from the intraband process considered here does not suffer from the noise associated with thermally excited electron-hole pairs in photodiodes.

Our analysis also applies to 2H-TMD monolayers, where transition metal and chalcogen atoms form trigonal crystal structures. Similar to gapped graphene, their band structures can also be described as massive Dirac fermions with trigonal warping TMD_warping , but with much larger band gaps TMD_review ; TMD_review2 . Hence they should exhibit large intraband second-order response as well. Topological insulator surface states also shows the effect (see Supplementary Materials). In addition to 2D systems, 3D bulk materials without inversion are also worth considering. As discussed, we expect large responsivities in low carrier densities. From this respect, inversion-breaking Weyl semimetals are promising candidates.

Discussion

So far we have considered the skew scattering contribution to the second-order DC response. Besides, the Berry curvature $\bm{\Omega}$ contributes to the distribution function through the collision integral and the electron’s velocity as the anomalous and side-jump velocities, which modifies the second-order conductivity. The Berry curvature is odd (even) under time reversal (inversion), and thus it is allowed to exist in time-reversal-invariant noncentrosymmetric materials. A similar Boltzmann transport analysis for linear response is reported in the context of the anomalous Hall effect without time-reversal symmetry AHE ; see Materials and Methods, and Supplementary Materials for detail. Finite Berry curvature enforces a coordinate shift $\delta\bm{r}_{\bm{k}\bm{k}^{\prime}}$ after scattering, giving a displacement of the center of an electron wave packet. It is accompanied by a potential energy shift $\delta\epsilon_{\bm{k}^{\prime}\bm{k}}$ in the presence of the external field, which results in the collision integral

[TABLE]

One may decompose the distribution function in nonequilibrium as $f=f_{0}+f^{\text{scatt}}+f^{\text{adist}}$, where $f_{0}$ is the equilibrium distribution function, $f^{\text{scatt}}$ describes the effect of scattering to the distribution (with $\delta\epsilon_{\bm{k}^{\prime}\bm{k}}=0$), and $f^{\text{adist}}$ is the anomalous distribution due to finite $\delta\epsilon_{\bm{k}^{\prime}\bm{k}}$. Because of the the different origins, they have distinct dependence on the scattering time $\tau$. For low frequencies with $\omega\tau\ll 1$, $\tau$ dependence of each term is easily estimated by a simple power counting. Noting that the scattering rate $w$ amounts to $\tau^{-1}$, we can expand the distribution function in powers of the electric field $E$, with coefficients depending on $\tau$:

[TABLE]

in the weak impurity limit.

The electric current density is obtained by $\bm{j}=-e\int_{\bm{k}}\bm{v}(\bm{k})f(\bm{k})$, where $\bm{v}(\bm{k})$ is the electron’s velocity,

[TABLE]

Besides the velocity given by the band dispersion $\bm{v}_{0}$, $\bm{v}$ contains two additional terms: The anomalous velocity $\bm{v}_{\text{av}}$ is associated with the Berry curvature and the last term $\bm{v}_{\text{sj}}$ describes the side-jump contribution, which arises from a combined effect of scattering and Berry curvature. Now, we can see several contributions to the second-order conductivity $\chi$ from combinations in the product of the out-of-equilibrium distribution function $(f^{\text{scatt}}+f^{\text{adist}})$ and the electron velocity $(\bm{v}_{0}+\bm{v}_{\text{av}}+\bm{v}_{\text{sj}})$.

The dominant contribution in the weak impurity limit $(\tau\to\infty)$ emerges from the skew scattering with $\bm{v}_{0}f^{\text{scatt}}_{2}$, resulting in $\chi\propto\tau^{2}$. The Berry curvature dipole contribution found in Ref. Sodemann corresponds to $\bm{v}_{\text{av}}f^{\text{scatt}}_{1}$, which amounts to $\chi\propto\tau$. It shows the frequency dependence as $\chi^{\prime}\propto e^{3}\tau D$, $\chi^{\prime\prime}\approx 0$ for low frequencies $(\omega\tau\ll 1)$ and $\chi^{\prime}\propto e^{3}D/(\omega^{2}\tau)$, $\chi^{\prime\prime}\propto e^{3}D/\omega$ for high frequencies $(\omega\tau\gg 1)$, where $D$ is the Berry curvature dipole. It is a quantity with the dimension of length in two dimensions. The Berry curvature dipole and skew scattering have the same frequency dependence although they are distinct in the $\tau$ dependence.

Inversion breaking results in finite skew scattering $w^{(A)}$. Along with the accompanying anisotropic Fermi surface, the second-order electronic response induced by skew scattering persists in any noncentrosymmetric materials. This is in marked contrast to the Berry curvature dipole mechanism, which imposes more symmetry constraints in addition to inversion breaking. Some symmetry classes without inversion, e.g., those containing three-fold rotation axis, do not show nonlinear response from a Berry curvature dipole Sodemann .

In summary, we have performed a systematic study of second-order electrical response due to intraband process and have identified the skew scattering mechanism as the dominant contribution in the weak impurity limit. This mechanism is responsible for current rectification and opens a new way to low-power energy harvesters and terahertz detection.

Materials and Methods

We employ the semiclassical Boltzmann transport theory to analyze the second-order response in noncentrosymmetric materials. The analysis is to some extent in parallel with that for the anomalous Hall effect Belinicher2 ; AHE , but there is a fundamental difference in symmetry: Inversion should be broken for finite second-order response whereas time-reversal symmetry is broken for the anomalous Hall effect.

The total Hamiltonian is given by

[TABLE]

where $H_{0}$ is the Hamiltonian for electrons in the clean limit with the eigenvalue $\epsilon_{\bm{k}}$, $V$ describes electron scattering, and $U$ includes the external electric field, which drives the system into nonequilibrium. In a semiclassical description, $U$ results in the force acting on an electron’s wave packet: $\nabla_{\bm{r}}U=-\bm{F}$. In the following, we consider an external electric field of frequency $\omega$, namely, $\bm{F}=-e\bm{E}(\omega)$.

A formally rigorous analysis requires a quantum-mechanical calculation of density matrices, which usually requires considerable effort. The semiclassical analysis is simpler and offers intuitive understanding; the downside is that it has limitations owing to the uncertainty principle. The semiclassical description relies on wave packets of Bloch states. A wave packet is localized both in the real space and the momentum space so that the mean position $\bm{r}$ and momentum $\bm{k}$ can be designated. However, the uncertainty principle imposes the limitation $\Delta r\Delta k\gtrsim 1$. We require position and momentum resolutions set by the mean free path $\ell$ and the Fermi wavevector $k_{F}$; i.e., $k_{F}\ell\gg 1$ is necessary for the semiclassical analysis. In the following, we deal with electron scattering by impurities. Since $\ell$ and $k_{F}$ depends on the impurity and electron densities, respectively, the condition is satisfied when the impurity density is sufficiently lower than the electron density.

The semiclassical Boltzmann equation for the distribution function $f$ includes the scattering rate $w_{\bm{k}^{\prime}\bm{k}}$ and the energy shift $\delta\epsilon_{\bm{k}^{\prime}\bm{k}}$. The scattering rate is given by

[TABLE]

with the scattering $T$ matrix $T_{\bm{k}^{\prime}\bm{k}}$. The energy shift appears for a scattering process in the presence the Berry curvature $\bm{\Omega}$. It is defined from the Berry connection $A_{a}=i\Braket{u_{\bm{k}}}{\partial_{k_{a}}u_{\bm{k}}}$ as $\Omega_{a}=\varepsilon_{abc}\partial_{k_{b}}A_{c}$, where $\ket{u_{\bm{k}}}$ is the lattice periodic part of the Bloch function at $\bm{k}$ and $\varepsilon_{abc}$ is the Levi-Civita symbol. The Berry curvature transforms under inversion $\mathcal{P}$ and time reversal $\mathcal{T}$ as

[TABLE]

It is deduced that the Berry curvature vanishes everywhere in the Brillouin zone when $\mathcal{PT}$ symmetry exists. Finite Berry curvature causes the coordinate shift $\delta\bm{r}_{\bm{k}^{\prime}\bm{k}}$, which describes the displacement of the center of a wave packet after a scattering process $\bm{k}\to\bm{k}^{\prime}$. For a weak scattering with a small momentum change, the coordinate shift can be approximated as $\delta\bm{r}_{\bm{k}^{\prime}\bm{k}}\approx(\bm{k}^{\prime}-\bm{k})\times\bm{\Omega}(\bm{k})$. The energy shift $\delta\epsilon_{\bm{k}^{\prime}\bm{k}}$ is given by using the coordinate shift $\delta\bm{r}_{\bm{k}^{\prime}\bm{k}}$ as

[TABLE]

It is worth noting that the coordinate shift is independent of the scattering time $\tau$ even though it is related to scattering as it describes the displacement after a single impurity scattering. Further details of the analysis are presented in the Supplementary Materials.

Acknowledgments

We thank Q. Ma, M. Soljačić, B. Halperin, and L. Liu for helpful discussions. Funding: The work was supported by the US Army Research Laboratory and the US Army Research Office through the Institute for Soldier Nanotechnologies. L.F was supported in part by a Simons Investigator Award from the Simons Foundation. Author contributions: H.I. and L.F. conceived the idea. H.I. performed the theoretical calculations. All authors discussed the results and wrote the manuscript. Competing interests: The authors are part of inventors on a patent application related to this work filed by the MIT (application no. 62/779,025). They declare that they have no other competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

S1 Semiclassical Boltzmann theory

We describe in detail the analysis of the semiclassical Boltzmann transport theory for the calculation of second-order response in noncentrosymmetric materials. We deal with spatially homogeneous systems and hence the electron distribution function $f(\bm{k},\epsilon)$ does not depend on the spatial position, but on the wavevector $\bm{k}$ and the energy $\epsilon$. We assume that the frequency of the external field $\omega$ is lower than the interband spacing, so that interband transitions are suppressed and negligible.

We calculate the distribution function in the presence of the external field by the semiclassical Boltzmann equation $(\hbar=1)$ Luttinger1 ; Sinitsyn_review

[TABLE]

The time derivative of the momentum is equal to the force felt by an electron wave packet: $\dot{\bm{k}}=\bm{F}$. We have an external electric field $\bm{E}(\omega)$, so that the force is $\bm{F}=-e\bm{E}(\omega)$. The collision integral $C[f]$ is given by Luttinger1 ; Sinitsyn2

[TABLE]

where $w_{\bm{k}^{\prime}\bm{k}}$ is the scattering rate from a state with momentum $\bm{k}$ to one with $\bm{k}^{\prime}$. In the following, we consider elastic scattering. Then, the scattering rate becomes

[TABLE]

with the scattering $T$ matrix $T_{\bm{k}^{\prime}\bm{k}}$. It is given by $T_{\bm{k}^{\prime}\bm{k}}=\braket{\bm{k}^{\prime}}{V}{\psi_{\bm{k}}}$, where $\ket{\bm{k}}$ is the eigenstate of $H_{0}$, the Hamiltonian in the clean limit, and $\ket{\psi_{\bm{k}}}$ is the eigenstate of $H_{0}+V$, including the scattering $V$. $\ket{\psi_{\bm{k}}}$ can be obtained as the solution to the Lippman–Schwinger equation $\ket{\psi_{\bm{k}}}=\ket{\bm{k}}+(\epsilon_{\bm{k}}-H_{0}+i\delta)^{-1}V\ket{\psi_{\bm{k}}}$. The energy shift is obtained as $\delta\epsilon_{\bm{k}^{\prime}\bm{k}}=-\bm{F}\cdot\delta\bm{r}_{\bm{k}^{\prime}\bm{k}}$ with the coordinate shift $\delta\bm{r}_{\bm{k}^{\prime}\bm{k}}\approx(\bm{k}^{\prime}-\bm{k})\times\bm{\Omega}(\bm{k})$ for a weak scattering process with a small momentum change Sinitsyn1 .

S1.1 Scattering rate in noncentrosymmetric media

In noncentrosymmetric media, the probability of a scattering process $\bm{k}\to\bm{k}^{\prime}$ and the inverted process $-\bm{k}\to-\bm{k}^{\prime}$ can be different. For time-reversal systems, this imbalance is captured by decomposing the scattering rate into two parts:

[TABLE]

where the symmetric and antisymmetric parts $w^{(S)}$ and $w^{(A)}$, respectively, is defined by

[TABLE]

When time reversal is preserved, the scattering rate satisfies the reversibility $w_{\bm{k}\bm{k}^{\prime}}=w_{-\bm{k}^{\prime},-\bm{k}}$, which leads to the relations

[TABLE]

Those equalities do not hold in general when time reversal is broken since states at $\bm{k}$ and $-\bm{k}$ usually have different energies. The optical theorem for elastic scattering guarantees the relation

[TABLE]

When the elastic scattering is due to impurities, $w^{(S)}$ and $w^{(A)}$ are obtained to the lowest order in the impurity scattering potential by

[TABLE]

where $\langle\ \rangle$ denotes the average over the impurity distribution and $V_{\bm{k}^{\prime}\bm{k}}$ is the matrix element for the single impurity scattering. We can see that the symmetric part of the scattering rate is obtained at the lowest order of the Born approximation, whereas the antisymmetric part is found at the next leading order. The detailed balance is broken when the antisymmetric part is finite.

Reference Belinicher2 reviews the photogalvanic effect in noncentrosymmetric materials, considering the asymmetry of scattering. The photovoltaic effect is studied with the Boltzmann equation and the antisymmetric component of the scattering rate $w_{\bm{k}^{\prime}\bm{k}}^{(A)}$, similarly to our analysis. A difference can be found in the origin of finite $w_{\bm{k}^{\prime}\bm{k}}$. In Ref. Belinicher2 , the asymmetry is imposed on scattering potentials for impurity scattering, ionization, photoexcitation, and recombination. We note that the asymmetry of the scattering rate can be finite due to an asymmetric scattering potential and also a wavefunction of a noncentrosymmetric medium. For the latter case, a scattering potential does not need to be inversion asymmetric to have finite $w_{\bm{k}^{\prime}\bm{k}}^{(A)}$, but even an isotropic scattering potential can generate $w_{\bm{k}^{\prime}\bm{k}}^{(A)}$ through a wavefunction.

S1.2 Formal solutions

The semiclassical Boltzmann equation (S1) should be solved self-consistently to obtain the distribution function $f(\bm{k})$. A fully self-consistent solution is generally difficult to obtain; here we decompose the distribution function $f$ as

[TABLE]

The first term $f_{0}$ is the distribution function in equilibrium, i.e., the Fermi–Dirac distribution. The second term $f^{\text{scatt}}$ describes the scattering contribution without the Berry curvature and the last term $f^{\text{adist}}$ is the anomalous distribution due to the Berry curvature and the energy shift $\delta\epsilon_{\bm{k}^{\prime}\bm{k}}$. We further expand $f^{\text{scatt}}$ and $f^{\text{adist}}$ with respect to the electric field $\bm{E}(\omega)$ and the asymmetric part of the scattering rate $w^{(A)}$:

[TABLE]

The subscript $n$ and the superscript $m$ correspond to the orders of the electric field $E$ and the antisymmetric scattering rate $w^{(A)}$, respectively. As we have seen in Eqs. (S9) and (S10), $w^{(A)}$ is smaller than $w^{(S)}$, so that we treat the former as a perturbation. Then, the Boltzmann equation is decomposed for each $f_{n}^{(m)}$ or $g_{n}^{(m)}$ to become

[TABLE]

A prime symbol ′ added to the distribution functions stands for the derivative with respect to the energy: $f^{\prime}(\bm{k},\epsilon)=\partial_{\epsilon}f(\bm{k},\epsilon)$.

Each equation has to be solved self-consistently; however it is usually not easy. Here, we solve a part of the collision integral involving $w^{(S)}$ self-consistently, while the other part with $w^{(A)}$ is treated as a perturbation. The collision integrals with $w^{(S)}$ define the scattering times $\tau_{n}^{(m)}$ and $\tau_{n}^{(m)}$ as follows:

[TABLE]

If $\tau_{n}^{(m)}$ and $\tau_{n}^{\prime(m)}$ do not satisfy the self-consistency, one may regard the equations above as the definitions for the relaxation time approximation. We also note the notation of the scattering times $\tau_{n}$ and $\tau^{\prime}_{n}$ used in the main text; they correspond to $\tau_{n}\equiv\tau_{n}^{(0)}$ and $\tau^{\prime}_{n}\equiv\tau_{n}^{(1)}$.

We assume $\tau$ and $\tau^{\prime}$ are independent of momentum $\bm{k}$ in solving the equations. Then, we obtain the formal solutions for the distribution functions $f_{n}^{(m)}$ and $g_{n}^{(m)}$. The second-order response contains different frequencies, [math] and $\pm 2\omega$. Since now we focus on rectification, i.e., zero-frequency response, we only write down the second-order solutions with $\omega=0$. Now, we list below the formal solutions up to the second order in the electric field $(n\leq 2)$ and to the first order in $w^{(A)}$ $(m\leq 1)$:

[TABLE]

$\tau_{n\omega}^{(m)}$ and $\tau_{n\omega}^{\prime(m)}$ are the shorthand notations for

[TABLE]

Those formal solutions find the scattering time dependence of the distribution functions for low frequencies $(\omega\tau\ll 1)$:

[TABLE]

Higher-order corrections for the anomalous distribution $(n\geq 2)$ has a nested structure of $f_{n-1}^{\text{scatt}}$ and $f_{n-1}^{\text{adist}}$. In the weak disorder limit of $\tau\to\infty$, we obtain

[TABLE]

We can define another scattering time $\tilde{\tau}$ from the antisymmetric part $w^{(A)}$. It is roughly speaking given by

[TABLE]

We note that the definition of $\tilde{\tau}$ is accompanied with the distribution function because the optical theorem for elastic scattering concludes Eq. (S8).

S2 Current response

From the distribution function, the electric current response is obtained by

[TABLE]

$\bm{v}(\bm{k})$ is the electron’s group velocity, given by

[TABLE]

The group velocity from the energy band dispersion $\bm{v}_{0}$ is independent of the Berry curvature, while the Berry curvature induces the anomalous velocity $\bm{v}_{\text{av}}$ Karplus ; Luttinger3 ; Adams and the side-jump velocity $\bm{v}_{\text{sj}}$ Berger :

[TABLE]

In terms of the electric field $E$, $\bm{v}_{\text{av}}$ is linear in $E$ whereas $\bm{v}_{0}$ and $\bm{v}_{\text{sj}}$ do not depend on $E$. Lastly, only the side-jump velocity depends on the scattering time: $\bm{v}_{\text{sj}}\propto\tau^{-1}$. Since there is no current in equilibrium, the electric current is described by AHE ; Sinitsyn_review

[TABLE]

which allows us to decompose the current into contributions of different origins. We note that the contribution from the anomalous distribution $f^{\text{adist}}$ could be recognized as a part of the side-jump effect since it is originated both from the Berry curvature and scattering.

Before calculating the second-order conductivity, we note the problem about the energy dissipation or Joule heating because of $\bm{j}\cdot\bm{E}\neq 0$; see also e.g., Ref. Luttinger1 for discussions. This concludes that a stationary state solution cannot be obtained for a closed system with energy conservation. To circumvent this issue, we neglect parts of the distribution function which are isotropically coupled to the electric field. To be more explicit, the distribution function at the second order in the electric field involves a term proportional to $|\bm{E}|^{2}f^{\prime\prime}_{0}$. Technically, such an excess term arises when we substitute a distribution function $f_{n}$ $(n\geq 1)$ into the collision integral $\mathcal{C}[f]$, which involves only elastic scattering. Since this term is isotropic, it cannot be relaxed by elastic scattering, which conserves energy, so that a stationary state solution does not exist. Inelastic scattering resolves the problem as it does not conserve energy and excess energy is dissipated as heat. Note that an inelastic scattering time is typically much longer than an elastic scattering time. In calculating the second-order conductivity, we simply subtract and neglect such isotropic terms although they potentially change the temperature for a closed system. Importantly, they do not contribute to current from the band velocity since $\int_{\bm{k}}\bm{v}_{0}f^{\prime\prime}_{0}=\int_{\bm{k}}\nabla_{\bm{k}}f^{\prime}_{0}=0$.

We can see the six possible combinations in Eq. (S40). For the second-order current response, the scattering time dependence of each contribution for low frequencies $\omega\tau\ll 1$ is

[TABLE]

In the weak disorder limit of $\tau\to\infty$, the first term, the skew scattering contribution, $\bm{v}_{0}f^{\text{scatt}}$ predominates, which we focus on in this work. This is equivalent to neglect all Berry curvature related effects. The second-order current response arising from this term is

[TABLE]

We define the second-order conductivity by the relation

[TABLE]

and we obtain from Eqs. (S26) and (S27)

[TABLE]

$\chi_{abc}^{(0)}$ vanishes identically in time-reversal-invariant systems because of $\epsilon_{\bm{k}}=\epsilon_{-\bm{k}}$; it follows from $\bm{v}_{0}(\bm{k})=-\bm{v}_{0}(-\bm{k})$, and we see that the integrand is odd in $\bm{k}$ for $\chi_{abc}^{(0)}$. Therefore, the second-order contribution from the skew-scattering contribution arises from $\chi_{abc}^{(1)}$, i.e., $\chi_{abc}=\chi_{abc}^{(1)}$.

To calculate the second-order response, it is useful to express Eqs. (S46) and (S47) in the following forms:

[TABLE]

We also show the expression for $\chi_{abc}^{(0)}$ for reference although it vanishes. These expressions show that results at finite temperature can be obtained by considering the thermal broadening of the Fermi–Dirac distribution, where the Fermi energy $\epsilon_{F}$ is replaced by the chemical potential $\mu_{\text{chem}}$ with the carrier density fixed. More explicitly, from the following relation

[TABLE]

the second-order conductivity at finite temperature is obtained by

[TABLE]

S2.1 Linear response

We also calculate the linear response to obtain responsivities (defined below). From the Boltzmann equation, the linear current response is given by

[TABLE]

leading to the linear conductivity

[TABLE]

S3 Graphene-based models with trigonal lattice structures

Now we evaluate the current response for explicit models. We consider two graphene-based models, utilizing monolayer and bilayer graphene. The crystalline lattices belong to a crystallographic point group $D_{3h}$ and $C_{3v}$, respectively. We note that inversion is equivalent to the in-plane two-fold rotation in 2D models, which is absent in the present models. Around $K$ and $K^{\prime}$ points, $\bm{K}=\left(\frac{4\pi}{3\sqrt{3}},0\right)$ and $\bm{K}^{\prime}=\left(-\frac{4\pi}{3\sqrt{3}},0\right)$, the $k\cdot p$ Hamiltonian to second order is

[TABLE]

where $k_{x}$ and $k_{y}$ are measured from $K$ or $K^{\prime}$ point, labeled by $s=\pm 1$, and $k_{s}$ is defined by $k_{s}=k_{x}+isk_{y}$. The energy dispersion of the positive branch is given by

[TABLE]

and the corresponding normalized wavefunction is

[TABLE]

The Berry curvature is also calculated as

[TABLE]

The Hamiltonian Eq. (S56) is available to both monolayer and bilayer models (Fig. S1). The Pauli matrix $\sigma_{a}$ describes the two sublattices for monolayer graphene and two layers for AB (Bernal) stacked bilayer graphene. Another difference is found in which of $v$ and $\lambda$ is dominant, which is explained later. $\Delta$ describes the potential energy difference at the two sublattice sites or layers, which opens a gap to the energy spectrum. In the following, we consider the two cases separately.

S3.1 Monolayer graphene

For monolayer graphene, the terms with $v$ dominate those with $\lambda$ $(vk_{F}\gg\lambda k_{F}^{2})$, where $\lambda$ characterizes the trigonal warping of the Fermi surface. $s=+1$ corresponds to $K$ valley and $s=-1$ to $K^{\prime}$ valley. The Dirac mass can be induced for example by placing a graphene sheet on a hexagonal boron nitride (hBN) substrate. For a tight-binding model on a honeycomb lattice with the nearest-neighbor hopping amplitude $t$ and the distance between sites $a$, the velocity $v$ and the trigonal warping $\lambda$ are given by $v=3ta/2$ and $\lambda=3ta^{2}/8$, respectively graphene_review1 ; graphene_review2 .

We assume that the trigonal warping is small and obtain the energy dispersion $\epsilon_{\bm{k}}$ to first order in $\lambda$ as

[TABLE]

where the Fermi wavevector $k_{F}$ is given by

[TABLE]

The density of states (DOS) at the Fermi energy is

[TABLE]

and $\theta_{\bm{k}}$ denotes the polar angle of $\bm{k}$. The Fermi surface can be parameterized by $\theta_{\bm{k}}$ as

[TABLE]

To evaluate the scattering rate, we need the wavefunction on the Fermi surface. The normalized wavefunction on the Fermi surface is given by

[TABLE]

Also, it is useful to see the Berry curvature on the Fermi surface, which is given by

[TABLE]

The second-order conductivity is considered in the presence of an elastic impurity scattering. Here, we assume short-ranged random scalar impurities

[TABLE]

but restrict scattering within each valley. $V_{i}$ is the impurity potential strength, $n_{i}$ is the impurity density, and $\bm{r}_{j}$ denotes an impurity position, where the summation is taken over all impurity sites. Then, we obtain the scattering rates for the valley $s$ to first order in $\lambda$ as

[TABLE]

with $\theta=\theta_{\bm{k}}$ and $\theta^{\prime}=\theta_{\bm{k}^{\prime}}$ parametrizing the Fermi surface.

In order to obtain the distribution function, we need to determine the scattering times that satisfy the relation Eq. (S22). For the present Hamiltonian, the scattering times are obtained as

[TABLE]

These solutions are obtained to $O(\lambda^{0})$, which are enough to evaluate conductivities to order $O(\lambda^{1})$ (see the next paragraph). Also, there is no directional dependence in the scattering times at order $\lambda^{0}$. For clarity, we define

[TABLE]

and the dimensionless quantity

[TABLE]

which are defined to satisfy $\tau_{1}^{(0)}=\tau_{1}^{(1)}=\tau_{2}^{(1)}=\tau$ and $\tau_{1}^{(1)}=\gamma\tau$.

We calculate the second-order contribution from each valley separately, namely $\chi_{abc}=\chi_{abc}^{(0)}+\chi_{abc}^{(1)}$ with

[TABLE]

The contributions from each valley $s$ are obtained from Eqs. (S46) and (S47) at zero temperature as

[TABLE]

where the factor of two from spin is multiplied. We note that $D_{3h}$ symmetry concludes that the imaginary part of the second-order DC conductivity vanishes Belinicher2 and that the system does not have second-order response to a circularly-polarized field. We define another scattering time $\tilde{\tau}$, related to $w^{(A)}$ [see Eq. (S35)], as

[TABLE]

This quantity measures the amounts of inversion breaking $(m)$ and the Fermi surface warping $(\lambda)$. The ratio $\tau/\tilde{\tau}$ becomes

[TABLE]

Roughly speaking, the ratio is determined by the product of the following dimensionless quantities: inversion breaking $\Delta/\epsilon_{F}$, the Fermi surface warping $\lambda k_{F}^{2}/(vk_{F})$, and the impurity strength $V_{i}k_{F}^{2}/\epsilon_{F}$. We can see from Eq. (S68) that the antisymmetric scattering becomes finite when inversion is broken by $\Delta\neq 0$ but it does not require finite Fermi surface warping $\lambda$. However, the second-order response vanishes when the Fermi surface warping is absent. The ratio $\tau/\tilde{\tau}$ measures the amount of the contribution to the second-order response from skew scattering, but not the magnitude of skew scattering itself.

$\chi^{(0)}$ is originated from the absence of time-reversal or inversion symmetry within each valley, resulting from the Fermi surface anisotropy. However, time reversal is preserved with the two valleys, and thus contributions from different valleys cancel. This effect of $\chi^{(0)}$ can be seen in Fig. 2(c), where the two valleys generates the opposite static velocity $\bm{V}$. In contrast, $\chi^{(1)}$ adds up with the two valleys, since this contribution captures the opposite chirality of the wavefunctions around each valley, which is not visible from the Fermi surface shape. Finally, we obtain the nonvanishing elements of the second-order conductivity

[TABLE]

with the average Fermi velocity on the Fermi surface

[TABLE]

and the dimensionless function

[TABLE]

Here the factor of four in Eq. (S78) corresponds to spin and valley degrees of freedom. The finite elements of $\chi_{abc}$ are consistent with $D_{3h}$ symmetry of the lattice.

S3.1.1 Linear conductivity and responsivity

To estimate the magnitude of the second-order response and to calculate responsivities, we also calculate the linear conductivity, particularly $\sigma_{s,ab}^{(0)}$. From Eq. (S54), we obtain

[TABLE]

for each valley with the spin degrees of freedom included. The linear conductivity is diagonal, and hence we write

[TABLE]

Since the linear conductivity can be written as $\sigma=ne\mu$ with the electron density $n$ and the mobility $\mu$, in low frequencies we find the relation

[TABLE]

We note that the electron density is obtained by $n=k_{F}^{2}/\pi$.

The short-circuit current responsivity $\mathfrak{R}_{I}$ is defined as the ratio of the generated DC current $j_{2}$ to the power of the incident electric field absorbed by the sample. When a sample has a dimension of $L^{2}$ and the incident field is uniform on the sample, the current responsivity is given by

[TABLE]

We note that the dimension $L$ corresponds to either width or length of a sample, depending on whether the second-order response of interest is longitudinal or transverse to the incident electric field. We also define the voltage responsivity $\mathfrak{R}_{V}$ as the ratio of the voltage generated by the second-order response to the incident power:

[TABLE]

Both current responsivity $\mathfrak{R}_{I}$ and voltage responsivity $\mathfrak{R}_{V}$ depend on the sample dimension $L$ and the ratio $\tau/\tilde{\tau}$. To remove those factors, we define the reduced current responsivity $\eta_{I}$ and the reduced voltage responsivity $\eta_{V}$ as follows:

[TABLE]

For the present model, $\eta_{I}$ and $\eta_{V}$ are obtained as

[TABLE]

S3.2 Bilayer graphene

Bernal stacked bilayer graphene has $C_{3v}$ symmetry when the two layer have different potential energies. Because of the symmetry, we can use the same Hamiltonian Eq. (S56); however, $v$ is treated as a perturbation instead of $\lambda$ for the bilayer case $(\lambda k_{F}^{2}\gg vk_{F})$, as we can see from a tight-binging Hamiltonian graphene_review1 ; graphene_review2 . Treating $v$ as a perturbation, we obtain the energy dispersion $\epsilon_{\bm{k}}$ to first order in $v$ as

[TABLE]

The DOS at the Fermi energy $D_{0}$ is given by

[TABLE]

The normalized wavefunction on the Fermi surface to order $v$ is

[TABLE]

where $k_{F}$ is defined by

[TABLE]

The Berry curvature on the Fermi surface is

[TABLE]

Using the wavefunction, we obtain the symmetric and antisymmetric scattering rates $w^{(S)}$ and $w^{(A)}$, respectively:

[TABLE]

For the bilayer case, the scattering times become

[TABLE]

We define $\tau$ as

[TABLE]

resulting in $\tau_{1}^{(0)}=\tau$ and $\tau_{1}^{(1)}=\tau_{2}^{(0)}=\tau_{2}^{(1)}=\gamma^{-1}\tau$, with the same $\gamma$ for the monolayer case.

The contributions to the second-order conductivity from valley $s$ are obtained from Eqs. (S46) and (S47) at zero temperature as

[TABLE]

For the bilayer model, the ratio $\tau/\tilde{\tau}$ is given by

[TABLE]

Similarly to the monolayer case, $\chi_{s,abc}^{(0)}$ from each valley cancels but $\chi_{s,abc}^{(1)}$ adds up, leading to the second-order conductivity

[TABLE]

where for the bilayer model the average velocity on the Fermi surface $v_{F}$ is

[TABLE]

and the dimensionless function $\zeta(\omega)$ becomes

[TABLE]

S3.2.1 Linear conductivity and responsivity

We also evaluate the linear conductivity $\sigma_{s,ab}^{(0)}$ for each valley $s$ with the spin degrees of freedom included [Eq. (S54)]:

[TABLE]

Similarly to the monolayer case, we define

[TABLE]

For low frequencies, the mobility is calculated from the relation $\sigma=ne\mu$ by

[TABLE]

where the electron density is $n=k_{F}^{2}/\pi$.

For the bilayer case, the reduced current responsivity becomes

[TABLE]

and the reduced voltage responsivity is

[TABLE]

S4 Estimate of parameters

Parameters for a realistic evaluation are taken and estimated from experiments on monolayer and bilayer graphene. For monolayer graphene on an hBN substrate, the imbalance of the potential energy within a sublattice is induced by the substrate and the lattice mismatch of graphene and hBN creates a superlattice with a periodicity of about 13 nm. In addition to the original (primary) Dirac points, the moiré superlattice produces new (secondary) superlattice Dirac points, located at a different energy Louie . Those two have different velocities and energy gaps: at ambient pressure, $v=0.94\times 10^{6}\,\text{m/s}$ and $2\Delta=30\,\text{meV}$ for the original Dirac point, and $v=0.5\times 10^{6}\,\text{m/s}$ and $2\Delta=20\,\text{meV}$ for the superlattice Dirac points LeRoy ; Dean1 . The size of the gap depends on the superlattice periodicity and it could be enhanced by interaction Song . We note that even at the original Dirac point, the velocity is slightly modified. We consider the original Dirac points instead of the superlattice Dirac points. The latter have the valley degeneracy of six and anisotropic Dirac cones whereas the valley degeneracy of the former is two.

The scattering time is estimated from a transport measurement with a high-quality sample Wang . At 300 K, they observed the mobility $\mu\approx 90,000\,\text{cm}^{2}/(\text{V\,s})$ with the carrier density $n=1\times 10^{12}\,\text{cm}^{-2}$, which amounts to the scattering time $\tau\approx 1.13\,\text{ps}$. The corresponding frequency is $(2\pi\tau)^{-1}\approx 141\,\text{GHz}$ and the mean free path is $\ell\sim v_{F}\tau\approx 1.6\,\mu\text{m}$. We note that the Fermi wavelength at $n=0.5\times 10^{12}\,\text{cm}^{-2}$ is $\lambda_{F}\approx 0.05\,\mu\text{m}$, which is still shorter than the mean free path: $\lambda_{F}\ll\ell$. To estimate the impurity strength, we need to know the impurity concentration. From transport and scanning tunneling microscopy/spectroscopy Dean2 ; Xue , we approximate the impurity density as $n_{i}\approx 1\times 10^{9}\,\text{cm}^{-2}$, which results in $V_{i}\approx 2.77\times 10^{-13}\,\text{eV\,cm}^{2}$ ($nV_{i}\approx 277\,\text{meV}$ at $n=1\times 10^{12}\,\text{cm}^{-2}$). Since the velocity is only slightly modified from the value without a superlattice structure, we employ the value $\lambda/v=a/4$ obtained from a tight-binding model, where $a$ is the carbon atom spacing, given by $a=1.42\,\text{\AA}$ graphene_review1 ; graphene_review2 . Those values yields the ratio $\tau/\tilde{\tau}\approx 0.003$ at $n=1\times 10^{12}\,\text{cm}^{-2}$. In the present calculations, we do not consider the temperature dependence of the parameters or the origin of scattering and inhomogeneity of samples for simplicity, although the scattering time becomes longer in lower temperatures.

For bilayer graphene, high-quality bilayer graphene grown by chemical vapor deposition (CVD) and detached on hBN realizes the mobility $\mu=30,000\,\text{cm}^{2}/\text{(V\,s)}$ at 300 K Schmitz . An out-of-plane electric field (displacement field) $D$ opens a band gap of about $2\Delta=100\,\text{meV}$ with $D\approx 1.1\,\text{V/nm}$ Lu . The coefficient $\lambda=(2m)^{-1}$ is determined by the effective mass $m\approx 0.033m_{e}$ ($m_{e}$: electron mass) and $v\approx 1\times 10^{5}\,\text{m/s}$ bilayer1 ; bilayer2 . Evaluating at $n=1\times 10^{12}\,\text{cm}^{-2}$, we find the scattering time $\tau\approx 0.96\,\text{ps}$, corresponding to the frequency $(2\pi\tau)^{-1}\approx 166\,\text{GHz}$, and the mean free path $\ell\sim v_{F}\tau\approx 0.35\,\mu\text{m}$. We assume the impurity concentration $n_{i}\approx 1\times 10^{9}\,\text{cm}^{-2}$, considering its high mobility, to obtain the impurity strength $V_{i}\approx 1.06\times 10^{-13}\,\text{eV\,cm}^{2}$ ($nV_{i}\approx 106\,\text{meV}$ at $n=1\times 10^{12}\,\text{cm}^{-2}$). At the carrier density $n=1\times 10^{12}\,\text{cm}^{-2}$, we have the ratio $\tau/\tilde{\tau}\approx 0.018$.

The carrier density dependence of $\epsilon_{F}$, $v_{F}$, and $\sigma(0)$ is shown in Fig. S2, and the frequency dependence of $\sigma$, $\chi$, $\eta_{I}$, and $\eta_{I}$ is depicted in Fig. S3. The frequency and temperature dependence of the response are given in Figs. S4 and S5. We note that the scattering times used for the calculations are those at 300 K. The scattering time becomes longer in lower temperature, which makes the linear conductivity $\sigma$, the second-order conductivity $\chi$, and the reduced current responsivity $\eta_{I}$ larger, whereas the low- and high-frequency values of the reduced voltage responsivity $\eta_{V}$ barely change.

S5 Surface state of a topological insulator

A surface state of a topological insulator also lacks inversion and thus hosts the rectification effect Ganichev . Here we take an example of the surface state of the topological insulator Bi2Te3 TIsurface ; Ganichev , where the Hamiltonian is

[TABLE]

leading to the energy spectrum

[TABLE]

The origin $\bm{k}=0$ is located at $\Gamma$ point and the $k_{x}$ axis is parallel to $\Gamma K$ line. Unlike the previous two cases based on graphene, it has a single Dirac spectrum centered at $\Gamma$ point. Furthermore, the topology of the wavefunction in the bulk prohibits a mass gap of the surface states. In this case, a surface breaks inversion by itself and the symmetry is $C_{3v}$.

For the present Hamiltonian, short-ranged impurities without momentum dependence do not produce finite $w^{(A)}$, nor does the second-order response. This issue can be circumvented instead by considering Coulomb impurities, where the impurity potential is inversely proportional to the distance from the impurity position: $V(\bm{r})\propto\sum_{j}|\bm{r}-\bm{r}_{j}|^{-1}$. Assuming that the hexagonal warping of the Fermi surface is small, we have the matrix element of the Coulomb impurity potential on the Fermi surface as

[TABLE]

Then, we obtain the scattering rates

[TABLE]

with the numerical factor $c_{1}\approx 0.93$ obtained from the numerical integration over momenta.

Since there is a single Fermi surface and time-reversal symmetry relates the Fermi surface with itself, $\chi^{(0)}$ vanishes with the single Fermi surface. A finite second-order conductivity can be found in $\chi^{(1)}$ as follows:

[TABLE]

Here, we assume for simplicity the constant scattering time $\tau\approx(2\pi D_{0}n_{i}V_{c}^{2})^{-1}$ obtained from $w^{(S)}$, and the scattering time originated from $w^{(A)}$ is defined by

[TABLE]

The two numerical factors $c_{1}\approx 0.87$ and $c_{2}\approx 0.14$ are evaluated by numerical integrations.