## 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 wave function engineering.

## 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 marked change in the model for energy production, consumption, and transportation. Because of the explosive growth of wireless technologies and portable devices, there is now increasing effort toward 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 direct current (DC). 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. These diodes face two fundamental limitations (*1*, *2*). First, rectification requires a threshold input voltage *V*_{T} = *k*_{B}*T*/*e*, known as thermal voltage (about 26 meV 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 harvestors of electromagnetic field in the microwave and terahertz 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 (*3*, *4*), 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 χ(ω) is an intrinsic material property that characterizes the DC current generated by an external electric field oscillating at frequency ω. Rectification in a single homogeneous 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 (*5*, *6*) such as photocurrent and second-harmonic generation are extensively studied (*7*, *8*), much less is known about second-order response at radio, microwave, and infrared frequencies (*9*, *10*, *11*). In this direction, recent works have predicted intraband photocurrent (*12*, *13*) and second-order nonlinear Hall effect due to “Berry curvature dipole” (*14*) 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 (*14*). Very recently, this effect was observed in low-frequency (around 100 Hz) transport measurements on the two-dimensional (2D) transition metal dichalcogenide (TMD) 1T′-WTe_{2} (*15*, *16*). The second-order Hall conductivity of bilayer WTe_{2} is remarkably large (*15*), in agreement with its large Berry curvature dipole from the tilted Dirac dispersion (*17*). 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 (*18*–*26*). 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 wave functions via quantum-mechanical scattering rates. We identify a previously unknown contribution to rectification from skew scattering with nonmagnetic impurities. Skew scattering arises from the intrinsic chirality of Bloch wave functions 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. On the basis of 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 *E* of frequency ω. The second-order response also involves a 2ω 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* _{abc}* is the rank-three tensor for the second-order conductivity;

*E*=

_{a}*E*(ω) is the (complex) amplitude of the external electric field of frequency ω; and the indices

_{a}*a*,

*b*, and

*c*label the coordinates. χ

*satisfies*

_{abc}*j*is odd under inversion, while the electric field

_{a}*E*is even, finite χ

_{a}*is possible only in noncentrosymmetric media. The second-order response can be decomposed into two parts as*

_{abc}*′ = χ*

_{abc}*′ and χ*

_{acb}*″ = − χ*

_{abc}*″. The former describes the response to a linearly polarized field and the latter to a circularly polarized field (*

_{acb}*27*). Note that in an isotropic medium, second-order response to a linearly polarized field must vanish. In other words, a nonzero

Nonlinear responses are extensively studied in the optical frequency regime (*5*). 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 bandgap 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 (*27*). For a homogeneous system, the Boltzmann equation is given by (ħ = 1)*f*(** k**, ϵ) is the distribution function at energy ϵ and C[

*f*] denotes the collision integral. The time derivative of the wave vector is equal to the force felt by an electron

**. The distribution function**

*E**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 C[

*f*], whose explicit form is shown later in a general form. It includes the scattering rate

*w*

_{k′k}, the probability of the transition per unit time from a Bloch state with a wave vector

**to that with**

*k***′. 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.**

*k*To obtain a nonzero second-order conductivity χ* _{abc}*, we need to capture the effect of inversion breaking. For time-reversal invariant systems, the energy dispersion ϵ

*and the band velocity*

_{k}

*v*_{0}(

**) = ∇**

*k*_{k}ϵ

_{k}satisfy the conditions ϵ

_{k}= ϵ

_{−k}and

*v*_{0}(

**) = −**

*k*

*v*_{0}(−

**), which hold with or without inversion symmetry: That is, the absence of inversion is encoded in the wave function but not in the energy dispersion.**

*k*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 wave function ψ* _{k}*. For and only for noncentrosymmetric crystals, the transition rate from

**to**

*k***′ and from −**

*k***to −**

*k***′ can be different:**

*k**w*

_{k′k}≠

*w*

_{−k′,−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*

_{k′k}=

*w*

_{−k,−k′}, which allows us to write

The transition rate *w*_{k′k} is usually calculated to the lowest order in scattering potential by Fermi’s golden rule. For elastic impurity scattering, it is given by *V*_{k′k} is the matrix element of a single scatterer and 〈 〉 denotes the impurity average. However, this formula is symmetric under the exchange ** k ↔ k**′ and does not capture skew scattering. The latter arises at the next-leading order from the interference of a direct transition

**′ and an indirect process**

*k → k***′ via an intermediate state**

*k → k***(**

*q**28*)

_{q}= ∫

*d*/(2π)

^{d}q*(*

^{d}*d*is the spatial dimension).

It is well known that skew scattering provides an extrinsic contribution to the anomalous Hall effect (*29*) and the spin Hall effect (*30*). 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 symmetries. With regard to its microscopic origin, skew scattering can be caused by nonmagnetic impurities when the single-impurity potential *V*(** r**) is inversion asymmetric, such as a dipole or octupole potential (

*27*). Here, we show that even for an inversion-symmetric single-impurity potential

*V*(

**) =**

*r**V*( −

**) such as a delta function potential, skew scattering can also appear. In this case, skew scattering is due to the inherent chirality of the electron wave function in a noncentrosymmetric crystal. In the simplest case of a delta function impurity potential, we have**

*r**V*

_{k′k}∝ 〈

*u*

_{k′}∣

*u*

_{k}〉 where∣

*u*

_{k}〉 is the Bloch wave function within a unit cell. Then, the expression for

*w*

^{(A)}given in Eq. 3 involves the wave function overlap at different

**points on the Fermi surface, similar to the Berry phase formula (**

*k**31*), and also, skew scattering is known to be related to Berry curvature in a certain limit (

*32*). While there is no net chirality in nonmagnetic systems, the Bloch wave function ∣

*u*〉 of an electron at finite momentum is allowed to be complex and carry

_{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**

*k***and −**

*k***. Then, skew scattering occurs because of 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 (**

*k**33*), and the viscosity results from scattering with impurities; see Fig. 2A.

To obtain the second-order conductivity, we should calculate the distribution function *f*(** k**) to second order in the external electric field

**(ω). Moreover, since skew scattering**

*E**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

**and**

*E**w*

^{(A)}, we obtain the second-order DC current response

**= −**

*j**e*∫

_{k}

*v*_{0}(

**)**

*k**f*(

**) due to skew scattering in time-reversal invariant systems. The second-order conductivity χ**

*k**can be formally expressed in terms of various scattering times*

_{abc}*(ω) is a shorthand for τ*

_{n}*(ω) = τ*

_{n}*/(1 −*

_{n}*i*ωτ

*) and likewise for*

_{n}*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 the Supplementary Materials.

A rough order-of-magnitude estimate of χ is obtained from Eq. 4 by using two scattering times τ and *t* = τ for ωτ ≪ 1, and Δ*t* = 1/ω for ωτ ≫ 1. Here, *v*_{F} is the Fermi velocity, *p*_{F} = **ħ***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 ω → 0, the ratio of the second-order current *j*_{2} and the linear response current *j*_{1} = σ_{ω}*E* is a product of two dimensionless quantities *eE*τ/*p*_{F} and

The short-circuit current responsivity *L*^{2}, the current responsivity is given by _{ω} ∼ *env*_{F} · (*e*λ_{F}τ/ħ)Re [(1 − *i*ωτ)^{−1}] (with the Fermi wavelength λ_{F} = 2π/k_{F}) obtained from the same Boltzmann equation, we have

Here, we define the reduced current and voltage responsivities η* _{I}* and η

*, respectively, which are independent of the sample size and the ratio*

_{V}*and η*

_{I}*hold in both low- and high-frequency limits. We also consider the ratio of the generated DC power and input power*

_{V}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 TMDs

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 (*34*–*38*), electrically biased bilayer graphene (*39*–*41*), and multilayer graphene such as ABA-stacked trilayer graphene (*42*, *43*). 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 wave functions and contributes to the second-order conductivity in graphene heterostructures with gaps induced by inversion symmetry breaking. The *k* · *p* Hamiltonian that we consider here is*s* = ±1 denotes two valleys at *K* and *K*^{′}, respectively, and we define *k*_{±} = *k _{x}* ±

*ik*. This Hamiltonian has a bandgap of 2Δ 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 (

_{y}*34*–

*38*); and (ii) bilayer graphene with an out-of-plane electric field that breaks the layer symmetry and opens up gaps (

*39*–

*41*). The difference between the monolayer and bilayer cases lies in the relative strength of

*v*and λ. For the monolayer (bilayer) case,

The Bloch wave function of *H _{s}*(

**) has two components associated with sublattice (layer) degrees of freedom in monolayer (bilayer) graphene. When the bandgap is present (Δ ≠ 0), the wave function 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**K*and

*K*′ because of time reversal symmetry. Note that because of the threefold rotation symmetry of graphene, Berry curvature dipole vanishes (

*14*), leaving skew scattering as the only mechanism for rectification.

We calculate the second-order conductivity from Eq. 4 and find nonvanishing elements χ* _{xxy}* = χ

*= χ*

_{xyx}*= − χ*

_{yxx}*≡ χ, 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 χ has the form*

_{yyy}*E*

_{s}(

**) ≠**

*k**E*

_{s}(−

**). Both ingredients are necessary for finite second-order response.**

*k*The second-order conductivity χ and the reduced voltage responsivity η* _{V}* for both monolayer and bilayer gapped graphene are shown in Fig. 3 as the functions of frequency ν = ω/(2π). The frequency dependence of the second-order conductivity is qualitatively understood with Eq. 5. We can further simplify the relation in two dimensions for

The responsivity is affected by the linear conductivity σ as well as the second-order conductivity χ. 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 σ and χ decrease as ω^{−2}, the responsivity saturates in the high-frequency limit. With a sample size of 10 μm, the carrier density *n* = 0.5 × 10^{12} cm^{−2}, and the ratio

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 × 10^{12} cm^{−2} shown in Fig. 3, the Fermi energy for the monolayer case is ϵ_{F} ≈ 80 meV, so that due to the Pauli blocking, the interband transition does not occur at frequencies below 2ϵ_{F} ∼ 40 THz. The interband threshold frequency also increases with the bandgap, 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 uses 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 (*44*) but with much larger bandgaps (*45*, *46*). Hence, they should exhibit a large intraband second-order response as well. Topological insulator surface states also show the effect (see the 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 **Ω** 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 (*29*); see Materials and Methods and the Supplementary Materials for details. Finite Berry curvature enforces a coordinate shift δ*r*_{kk′} after scattering, giving a displacement of the center of an electron wave packet. It is accompanied by a potential energy shift δϵ_{k′k} in the presence of the external field, which results in the collision integral

One may decompose the distribution function in nonequilibrium as *f* = *f*_{0} + *f*^{scatt} + *f*^{adist}, where *f*_{0} is the equilibrium distribution function, *f*^{scatt} describes the effect of scattering to the distribution (with δϵ_{k′k} = 0), and *f*^{adist} is the anomalous distribution due to finite δϵ_{k′k}. Because of the the different origins, they have distinct dependence on the scattering time τ. For low frequencies with ωτ ≪ 1, τ dependence of each term is easily estimated by a simple power counting. Noting that the scattering rate *w* amounts to τ^{−1}, we can expand the distribution function in powers of the electric field *E*, with coefficients depending on τ

The electric current density is obtained by ** j** = −

*e*∫

_{k}**(**

*v***)**

*k**f*(

**), where**

*k***(**

*v***) is the electron’s velocity**

*k*Besides the velocity given by the band dispersion *v*_{0}, ** v** contains two additional terms: The anomalous velocity

*v*_{av}is associated with the Berry curvature, and the last term

*v*_{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 χ from combinations in the product of the out-of-equilibrium distribution function (

*f*

^{scatt}+

*f*

^{adist}) and the electron velocity (

*v*_{0}+

*v*_{av}+

*v*_{sj}).

The dominant contribution in the weak impurity limit (τ → ∞) emerges from the skew scattering with ^{2}. The Berry curvature dipole contribution found in (*14*) corresponds to *e*^{3}τ*D*, χ″ ≈ 0 for low frequencies (ωτ ≪ 1) and χ′ ∝ *e*^{3}*D*/(ω^{2}τ), χ″ ∝ *e*^{3}*D*/ω for high frequencies (ωτ ≫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 τ 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 a threefold rotation axis, do not show nonlinear response from a Berry curvature dipole (*14*).

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 use 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 (*27*, *29*), 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*H*_{0} is the Hamiltonian for electrons in the clean limit with the eigenvalue ϵ_{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: ∇*r**U* = −** F**. In the following, we consider an external electric field of frequency ω, namely,

**= −**

*F**e*

**(ω).**

*E*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 ** r** and momentum

**can be designated. However, the uncertainty principle imposes the limitation Δ**

*k**r*Δ

*k*≳ 1. We require position and momentum resolutions set by the mean free path 𝓁 and the Fermi wave vector

*k*

_{F}; i.e.,

*k*

_{F}𝓁 ≫ 1 is necessary for the semiclassical analysis. In the following, we deal with electron scattering by impurities. Since 𝓁 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*_{k′k} and the energy shift δϵ_{k′k}. The scattering rate is given by*T* matrix *T*_{k′k}. The energy shift appears for a scattering process in the presence the Berry curvature **Ω**. It is defined from the Berry connection *A _{a}* =

*i*〈

*u*

_{k}∣∂

_{ka}

*u*

_{k}〉 as Ω

*= ε*

_{a}*, where*

_{abc}∂_{kb}A_{c}**|**

*u*

_{k}〉 is the lattice periodic part of the Bloch function at

**and ε**

*k**is the Levi-Civita symbol. The Berry curvature transforms under inversion P and time reversal T as*

_{abc}It is deduced that the Berry curvature vanishes everywhere in the Brillouin zone when PT symmetry exists. Finite Berry curvature causes the coordinate shift δ*r*_{k′k}, which describes the displacement of the center of a wave packet after a scattering process ** k** →

**′. For a weak scattering with a small momentum change, the coordinate shift can be approximated as δ**

*k*

*r*_{k′k}≈ (

**′ −**

*k***) ×**

*k***Ω**(

**). The energy shift δϵ**

*k*_{k′k}is given by using the coordinate shift δ

*r*_{k′k}as

Note that the coordinate shift is independent of the scattering time τ, although 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.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/13/eaay2497/DC1

Section S1. Semiclassical Boltzmann theory

Section S2. Current response

Section S3. Graphene-based models with trigonal lattice structures

Section S4. Estimate of parameters

Section S5. Surface state of a topological insulator

Fig. S1. Monolayer and bilayer graphene models.

Fig. S2. Carrier density dependence of the material properties.

Fig. S3. Frequency dependence of the response.

Fig. S4. Frequency and temperature dependence of the response.

Fig. S5. Frequency and temperature dependence of the response with a logarithmic scale for the frequency axis.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank Q. Ma, M. Soljačić, B. Halperin, and L. Liu for helpful discussions.

**Funding:**The work was supported by the U.S. Army Research Laboratory and the U.S. 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). The authors 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.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).