CURRENT RECTIFICATION BASED ON NONCENTROSYMMETRIC QUANTUM MATERIALS

Abstract

Rectification is a process that converts electromagnetic fields into direct current (DC). Such a process underlies a wide range of technologies, including wireless communication, wireless charging, energy harvesting, and infrared detection. Existing rectifiers are mostly based on semiconductor diodes, with limited applicability to small voltages or high frequency inputs. Here, we present an alternative approach to current rectification that uses the electronic properties of quantum crystals without semiconductor junctions. We identify a new mechanism for rectification from skew scattering due to the 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. These effects can be used in high-frequency rectifiers by rational material design and quantum wavefunction engineering.

Claims

1. A current rectifier comprising: an antenna to receive an alternating current (AC) electric field; a noncentrosymmetric material, in electrical communication with the antenna, to generate a direct current (DC) current via a second-order nonlinear response by the noncentrosymmetric material to the AC electric field received by the antenna; and a pair of electrodes, in electrical communication with the noncentrosymmetric material, to conduct the DC current.

2. The current rectifier of claim 1, wherein the AC electric field oscillates at a frequency from about 0.1 THz to about 10 THz.

3. The current rectifier of claim 1, wherein the AC electric field oscillates in a first direction in a plane of the noncentrosymmetric material and the DC current flows from the noncentrosymmetric material in a second direction orthogonal to the first direction and in the plane of the noncentrosymmetric material.

4. The current rectifier of claim 1, wherein the AC electric field oscillates in a first direction in a plane of the noncentrosymmetric material and the DC current flows from the noncentrosymmetric material in the first direction.

5. The current rectifier of claim 1, wherein the noncentrosymmetric material lacks symmetry in a direction perpendicular to a plane of the noncentrosymmetric material.

6. The current rectifier of claim 1, wherein the noncentrosymmetric material has a band gap of about 26 meV at room temperature.

7. The current rectifier of claim 1, wherein the noncentrosymmetric material is configured to generate the DC current from the AC electric field at room temperature.

8. The current rectifier of claim 1, wherein the noncentrosymmetric material comprises a graphene monolayer disposed on a dielectric layer that breaks an inversion symmetry of the graphene monolayer.

9. The current rectifier of claim 1, wherein the noncentrosymmetric material comprises a graphene bilayer and further comprising: a first insulator layer disposed on a first side of the graphene bilayer; a first gate electrode, on a far side of the first insulator layer from the graphene bilayer, to apply a first gate voltage to the graphene bilayer; a second insulator layer disposed on a second side of the graphene bilayer opposite the first side of the graphene bilayer; and a second gate electrode, on a far side of the second insulator layer from the graphene bilayer, to apply a second gate voltage to the graphene bilayer, the first gate voltage and the second gate voltage breaking an inversion symmetry of the graphene bilayer.

10. The current rectifier of claim 1, wherein the noncentrosymmetric material comprises a transition metal dichalcogenide.

11. A method of rectifying current, the method comprising: receiving an alternating current (AC) electric field with an antenna; and coupling the AC electric field from the antenna into a noncentrosymmetric material that generates a direct current (DC) current via a second-order nonlinear response to the AC electric field received by the antenna.

12. The method of claim 11, wherein the AC electric field oscillates at a frequency from about 0.1 THz to about 10 THz.

13. The method of claim 11, wherein the AC electric field oscillates in a first direction in a plane of the noncentrosymmetric material and the DC current flows from the noncentrosymmetric material in a second direction orthogonal to the first direction and in the plane of the noncentrosymmetric material.

14. The method of claim 11, wherein the AC electric field oscillates in a first direction in a plane of the noncentrosymmetric material and the DC current flows from the noncentrosymmetric material in the first direction.

15. The method of claim 11, wherein the noncentrosymmetric material lacks symmetry in a direction perpendicular to a plane of the noncentrosymmetric material.

16. The method of claim 11, wherein the noncentrosymmetric material has a band gap of about 26 meV at room temperature.

17. The method of claim 11, wherein the noncentrosymmetric material generates the DC current from the AC electric field at room temperature.

18. The method of claim 11, wherein the noncentrosymmetric material comprises a graphene monolayer disposed on a dielectric layer that breaks a symmetry of the graphene monolayer.

19. The method of claim 11, wherein the noncentrosymmetric material comprises a graphene bilayer and further comprising: applying an electric field to the graphene bilayer to break a symmetry of the graphene bilayer.

20. The method of claim 11, wherein the noncentrosymmetric material comprises a transition metal dichalcogenide.

Description

BRIEF DESCRIPTIONS OF THE DRAWINGS

[0014] The skilled artisan will understand that the drawings primarily are for illustrative purposes and are not intended to limit the scope of the inventive subject matter described herein. The drawings are not necessarily to scale; in some instances, various aspects of the inventive subject matter disclosed herein may be shown exaggerated or enlarged in the drawings to facilitate an understanding of different features. In the drawings, like reference characters generally refer to like features (e.g., functionally similar and/or structurally similar elements).

[0015] FIGS. 1A and 1B show a high-frequency (e.g., THz) current rectifier with a noncentrosymmetric quantum material that converts an oscillating electric field into a transverse DC current.

[0016] FIG. 2 shows a high-frequency (e.g., THz) current rectifier with a noncentrosymmetric quantum material that converts an oscillating electric field into a longitudinal DC current.

[0017] FIG. 3A is a side view a noncentrosymmetric quantum material formed of a monolayer graphene/dielectric heterostructure where the monolayer graphene's inversion symmetry is broken by a dielectric layer.

[0018] FIG. 3B is a top view of the monolayer graphene/dielectric heterostructure of FIG. 3A.

[0019] FIG. 4 is a side view of a noncentrosymmetric quantum material formed of a bilayer graphene stack whose inversion symmetry is broken by an applied bias electric field.

[0020] FIG. 5A illustrates transverse (χ.sub.yxx) skew scattering on a honeycomb (e.g., graphene) lattice.

[0021] FIG. 5B illustrates the electric field and rectified current for transverse (χ.sub.yxx) skew scattering on a honeycomb lattice.

[0022] FIG. 5C illustrates Fermi surface displacement at frequency ω (green) for transverse (χ.sub.yxx) skew scattering on a honeycomb lattice.

[0023] FIG. 5D illustrates stationary Fermi surface displacement for transverse (χ.sub.yxx) skew scattering on a honeycomb lattice.

[0024] FIG. 6A illustrates longitudinal (χ.sub.yyy) skew scattering on a honeycomb (e.g., graphene) lattice.

[0025] FIG. 6B illustrates the electric field and rectified current for longitudinal (χ.sub.yyy) skew scattering on a honeycomb lattice.

[0026] FIG. 6C illustrates Fermi surface displacement at frequency ω (green) for longitudinal (χ.sub.yyy) skew scattering on a honeycomb lattice.

[0027] FIG. 6D illustrates stationary Fermi surface displacement for longitudinal (χ.sub.yyy) skew scattering on a honeycomb lattice.

[0028] FIG. 7A is a plot of the second-order conductivity ω for a graphene monolayer heterostructure.

[0029] FIG. 7B is a plot of the reduced voltage responsivity η.sub.V for a graphene monolayer heterostructure.

[0030] FIG. 7C is a plot of the second-order conductivity ω for a graphene bilayer heterostructure.

[0031] FIG. 7D is a plot of the reduced voltage responsivity η.sub.V for a graphene bilayer heterostructure.

[0032] FIG. 8A illustrates the lattice structure and crystal axes of monolayer graphene.

[0033] FIG. 8B illustrates the lattice structure and crystal axes of bilayer graphene.

[0034] FIG. 8C shows the Brillouin zone and energy contours for the monolayer graphene model of FIG. 8A.

[0035] FIG. 9A is a plot of the Fermi energy ε.sub.F for monolayer graphene (upper trace) and bilayer graphene (lower trace).

[0036] FIG. 9B is a plot of the Fermi velocity v.sub.F for monolayer graphene (flat trace) and bilayer graphene (sloped trace).

[0037] FIG. 9C is a plot of the linear DC conductivity σ at zero temperature for monolayer graphene (sublinearly increasing trace) and bilayer graphene (linearly increasing trace).

[0038] FIG. 10A is a plot of the linear conductivity σ of monolayer graphene versus frequency for carrier densities ranging from 0.5×10.sup.12 cm.sup.−2 (bottom trace) to 5×10.sup.12 cm.sup.−2 (top trace) and a temperature of zero degrees Kelvin.

[0039] FIG. 10B is a plot of the second-order conductivity χ of monolayer graphene versus frequency for carrier densities ranging from 0.5×10.sup.12 cm.sup.−2 to 5×10.sup.12 cm.sup.−2 and a temperature of zero.

[0040] FIG. 10C is a plot of the reduced current responsivity η.sub.I of monolayer graphene versus frequency for carrier densities ranging from 0.5×10.sup.12 cm.sup.−2 to 5×10.sup.12 cm.sup.−2 and a temperature of zero.

[0041] FIG. 10D is a plot of the reduced voltage responsivity η.sub.V of monolayer graphene versus frequency for carrier densities ranging from 0.5×10.sup.12 cm.sup.−2 to 5×10.sup.12 cm.sup.−2 and a temperature of zero.

[0042] FIG. 10E is a plot of the linear conductivity σ of bilayer graphene versus frequency for carrier densities ranging from 0.5×10.sup.12 cm.sup.−2 (bottom trace) to 5×10.sup.12 cm.sup.−2 (top trace) and a temperature of zero degrees Kelvin.

[0043] FIG. 10F is a plot of the second-order conductivity χ of bilayer graphene versus frequency for carrier densities ranging from 0.5×10.sup.12 cm.sup.−2 to 5×10.sup.12 cm.sup.−2 and a temperature of zero.

[0044] FIG. 10G is a plot of the reduced current responsivity η.sub.I of bilayer graphene versus frequency for carrier densities ranging from 0.5×10.sup.12 cm.sup.−2 to 5×10.sup.12 cm.sup.−2 and a temperature of zero.

[0045] FIG. 10H is a plot of the reduced voltage responsivity η.sub.V of bilayer graphene versus frequency for carrier densities ranging from 0.5×10.sup.12 cm.sup.−2 to 5×10.sup.12 cm.sup.−2 and a temperature of zero.

[0046] FIG. 11A is a linear-scale plot of the linear conductivity σ of monolayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0047] FIG. 11B is a linear-scale plot of the second-order conductivity χ of monolayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0048] FIG. 11C is a linear-scale plot of the reduced current responsivity η.sub.I of monolayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0049] FIG. 11D is a linear-scale plot of the reduced voltage responsivity η.sub.V of monolayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0050] FIG. 11E is a linear-scale plot of the linear conductivity σ of bilayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0051] FIG. 11F is a linear-scale plot of the second-order conductivity χ of bilayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0052] FIG. 11G is a linear-scale plot of the reduced current responsivity η.sub.I of bilayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0053] FIG. 11H is a linear-scale plot of the reduced voltage responsivity η.sub.V of bilayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0054] FIG. 12A is a logarithmic-scale plot of the linear conductivity σ of monolayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0055] FIG. 12B is a logarithmic-scale plot of the second-order conductivity χ of monolayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0056] FIG. 12C is a logarithmic-scale plot of the reduced current responsivity η.sub.I of monolayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0057] FIG. 12D is a logarithmic-scale plot of the reduced voltage responsivity η.sub.V of monolayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0058] FIG. 12E is a logarithmic-scale plot of the linear conductivity σ of bilayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0059] FIG. 12F is a logarithmic-scale plot of the second-order conductivity χ of bilayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0060] FIG. 12G is a logarithmic-scale plot of the reduced current responsivity η.sub.I of bilayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

[0061] FIG. 12H is a logarithmic-scale plot of the reduced voltage responsivity η.sub.V of bilayer graphene versus frequency and temperature for a carrier density n=1×10.sup.12 cm.sup.−2.

DETAILED DESCRIPTION

[0062] FIGS. 1A and 1B illustrate a current rectifier 100 that can rectify radio-frequency and terahertz electric fields (e.g., alternating current (AC) electric fields oscillating at 0.1 THz to 10 Hz). The current rectifier 100 includes an antenna 100 that is coupled to a thin patch of noncentrosymmetric material 120, which in turn is coupled to a pair of electrodes 130a and 130b (collectively, electrodes 130). The electrodes 130 connect the rectifier 100 to a load 140, which may have an impedance of ohms to kiloohms. In principle, there is no high-frequency coupling to these electrodes 130, but they may be DC-coupled (e.g., using a low-pass filter) to remove high-frequency electric field components and pass DC current.

[0063] In this example, the patch of noncentrosymmetric material 120 is rectangular, with two halves of the antenna 110 are coupled to one pair of opposing edges and the electrodes 130 coupled to the other pair of opposing edges. The antenna 110, noncentrosymmetric material 120, and electrodes 130 can be deposited on any suitable dielectric substrate (e.g., a silicon oxide substrate). Notably, this current rectifier 100 does not include or form a diode.

[0064] The antenna 110 is shown in FIGS. 1A and 1B as a bowtie antenna but could take any other suitable shape. It is made of metal or another conductive material and is configured to detect linearly polarized incident electric fields. For example, the bowtie antenna 110 has a long axis stretching horizontally in FIG. 1B and is configured to detect a horizontally polarized electric field, as indicated by the double-headed arrow representing the electric field oscillation direction in FIG. 1B. If the incident electric field is not linearly polarized along the bowtie antenna's long axis-for example, if the incident electric field is unpolarized or elliptically polarized-then the bowtie antenna 110 couples the horizontal projection of the incident electric field into the noncentrosymmetric material 120.

[0065] The noncentrosymmetric material 120 can include a two-dimensional material, such as monolayer graphene, bilayer graphene, multi-layer graphene, or a transition metal dichalcogenide, or a three-dimensional material, such as a bulk transition metal dichalcogenide or a Weyl semimetal without inversion (e.g., tantalum arsenide (TaAs)). It is formed in a roughly planar rectangular patch whose dimensions can be about 1-10 μm wide to about 1-10 μm long to about 10-100 nm thick. The patch may be formed in other shapes (e.g., another polygon) so long as the noncentrosymmetric material 120 is coupled to the antenna 110 and electrodes 130. As explained in greater detail below, this patch may include multiple layers, including dielectric layers and/or additional electrodes that break the inversion symmetry of a two-dimensional material.

[0066] The noncentrosymmetric material 120 lacks inversion symmetry along an axis parallel to the DC current. This axis is generally, but not always, perpendicular to the noncentrosymmetric material's largest surface (i.e., along the propagation direction of a normally incident electric field). For the transverse response illustrated in FIG. 1B, the noncentrosymmetric material 120 lacks inversion symmetry along an axis perpendicular to the electric field oscillation direction.

[0067] In operation, the antenna 110 couples an incident electric field E.sub.ω into the noncentrosymmetric material 120, which exhibits a second-order nonlinear response to the incident electric field E.sub.ω, as explained in greater detail below. More specifically, skew scattering by the lattice of the noncentrosymmetric material 120 produces a net DC current J.sub.DC. In the rectifier 100 shown in FIGS. 1A and 1B, this scattering is transverse—that is, the DC current J.sub.DC propagates in a direction that is perpendicular to the oscillation direction (indicated by the double-headed arrow in FIG. 1B) of the electric field E.sub.ω. The current direction is set by the polarization of the electric field E.sub.ω and the crystal asymmetry. The DC current J.sub.DC flows within the layer of noncentrosymmetric material 120 to the load 140 via the electrodes 130.

[0068] FIG. 2 shows a current rectifier 200 that generates a “longitudinally” propagating rectified current that J.sub.DC from a high-frequency (e.g., 0.1-10 THz) electric field E.sub.ω. Like the current rectifier 100 in FIGS. 1A and 1B, this current rectifier 200 includes an antenna 210; a thin, rectangular patch of noncentrosymmetric material 220; and electrodes 230a and 230b (collectively, electrodes 230) on a substrate 202. The electrodes 230 are electrically connected to a load 240, which may also be on the same substrate 202 as the current rectifier 200. In this case, however, the antenna elements and electrodes 230 are coupled to the same sides of the noncentrosymmetric material 220.

[0069] In operation, the antenna 210 couples the electric field E.sub.ω into the noncentrosymmetric material 220 such that the electric field E.sub.ω oscillates in the direction indicated by the double-headed arrow in FIG. 2. Again, the noncentrosymmetric material 220 generates a net DC current J.sub.DC in response to the electric field E.sub.ω via a second-order skew scattering response. In this case, however, the DC current J.sub.DC propagates in a direction parallel to the electric field E.sub.ω oscillation direction. The DC current J.sub.DC flows within the plane of noncentrosymmetric material 220 to the load 240 via the electrodes 230.

[0070] FIGS. 3A, 3B, and 4 show different examples of noncentrosymmetric heterostructures suitable for use as the noncentrosymmetric materials in the current rectifiers shown in FIGS. 1A, 1B, and 2. More specifically, FIGS. 3A and 3B show side and top views, respectively, of a noncentrosymmetric heterostructure 320 formed of monolayer graphene 322 on a layer of dielectric boron nitride (BN) 324 on a substrate 302. This heterostructure 320 may be up to about 100 nm thick. By itself, the monolayer graphene 322 has inversion symmetry. The lattice mismatch between the monolayer graphene 322 and the BN layer 324 breaks the monolayer graphene's inversion symmetry, causing the heterostructure 320 to be noncentrosymmetric. (A layer of silicon carbide (SiC) can be used instead of the BN layer 324 to break the monolayer graphene's inversion symmetry.)

[0071] FIG. 4 is a side view of a noncentrosymmetric heterostructure 420 with bilayer graphene 422 sandwiched between a pair of BN layers 424a and 424b (collectively, BN layers 424), which are themselves sandwiched between a pair of electrodes 426a and 426b (collectively, electrodes 426). Each BN layer 424a, 424b may be about 10-30 nm thick (e.g., 20 nm) and can have different thicknesses (i.e., the upper BN layer 424a can be thicker or thinner than the lower BN layer 424b). Together, the bilayer graphene 422, BN layers 424, and electrodes 426 form a double capacitive structure on a substrate 402.

[0072] Like monolayer graphene, the bilayer graphene 422 has inversion symmetry. This inversion symmetry is broken by applying bias fields V.sub.BG and V.sub.TG to the upper and lower layers of the bilayer graphene 422 with the electrodes 426. These bias fields are orthogonal to the surface of the substrate 402, to the polarization of the incident (AC) electric field, and to the rectified current. The bias fields generate a displacement field with the graphene bilayer 422 that are roughly proportional to (V.sub.BG/d.sub.1−V.sub.TG/d.sub.2), where d.sub.1 is the thickness of the upper BN layer 424a and d.sub.2 is the thickness of the lower BN layer 424b. The displacement field amplitude can range from about −1.5 V/nm to about +1.5 V/nm for applied voltages on the order of 2 V. The applied voltages and BN layer thicknesses are also selected based on the desired carrier density within the bilayer graphene 422. This carrier density is proportional to (V.sub.BG/d.sub.1+V.sub.TG/d.sub.2) and may be about ±10.sup.13 cm.sup.−2.

[0073] FIGS. 5A-5D and 6A-6D illustrate the second-order response by skew scattering on a honeycomb lattice for the longitudinal (χ.sub.yxx) and transverse (χ.sub.yyy) cases, respectively. FIGS. 5A and 6A show the schematics of skew scattering. When a self-rotating wave packet (large circles) is scattered by an inversion-symmetric potential (small circle), its motion is deflected like the Magnus effect. Two wave packets moving in the opposite directions produces zero net current. However, if the two self-rotate in different directions, skew scattering produces net current in the perpendicular direction.

[0074] FIGS. 5B and 6B show the electric field and rectified current on a honeycomb lattice. There are two schemes for measuring a strong rectification effect: the AC oscillating field can be along either the zigzag direction (FIG. 5B) or armchair direction (FIG. 6B) and the DC rectification signal is collected along the armchair direction.

[0075] FIGS. 5C and 6C shows Fermi surface displacement at frequency ω (solid outline). The oscillating electric field E(w) forces the Fermi surface to swing back and forth from its equilibrium position (dashed outline). Owing to the Fermi surface anisotropy, each valley yields finite velocity V along the k.sub.x direction after time averaging. This velocity is canceled with the two valleys, and there is no DC current generated as a linear response.

[0076] FIGS. 5D and 6D show the stationary Fermi surface displacement. The electric field and skew scattering produce the stationary Fermi surface displacement (solid outline) from the equilibrium state (dashed outline) as a second-order response. Finite rectified current is observed when the contributions from the two valleys do not cancel.

Results

Second-Order Response

[0077] 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

j.sub.a=χ.sub.abcE*.sub.bE.sub.c,   (1)

where χ.sub.abc is the rank-three tensor for the second-order conductivity, E.sub.a=E.sub.a(ω) is the (complex) amplitude of the external electric field of frequency ω, and the indices a, b, c label the coordinates. χ.sub.abc satisfies χ.sub.abc=χ*.sub.acb. Since the current j.sub.a is odd under inversion while the electric field E.sub.a is even, finite χ.sub.abc is possible only in noncentrosymmetric media. The second-order response can be decomposed into two parts as j.sub.a=χ.sub.abc′E*.sub.bE.sub.c+iχ.sub.abc″E*.sub.bE.sub.c, with χ.sub.abc′=χ.sub.acb′ and ω.sub.abc″=−χ.sub.acb″. The former describes the response to a linearly-polarized field and the latter to a circularly-polarized field. 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 χ′.sub.abc requires the presence of crystal anisotropy.

[0078] Nonlinear responses are studied extensively in the optical frequency regime. 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 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, here we consider intra-band nonlinear response at frequencies below the inter-band transition threshold.

Semiclassical Boltzmann Analysis

[0079] We analyze the second-order electrical transport using the semiclassical Boltzmann theory. For a homogeneous system, the Boltzmann equation is given by (ℏ=1):

[00001]$\begin{array}{cc}\frac{\partial f}{\partial t}+\stackrel{.}{k}.Math.\frac{\partial f}{\partial k}=-\mathrm{��}\left[f\right],& \left(2\right)\end{array}$

where ƒ(k, ε) is the distribution function at energy ε and [ƒ] denotes the collision integral. The time derivative of the wavevector is equal to the force felt by an electron {dot over (k)}=−eE under the external electric field E. The distribution function ƒ 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 [ƒ], whose explicit form is shown below in a general form. It includes the scattering rate w.sub.k,k, the probability of the transition per unit time from a Bloch state with a wavevector k to that with k′. We consider scattering by impurities. The impurity density should be sufficiently lower than the electron density for the semiclassical Boltzmann analysis to be valid; see below for further discussions.

[0080] To obtain a nonzero second-order conductivity χ.sub.abc, we should capture the effect of inversion breaking. For time-reversal-invariant systems, the energy dispersion ε.sub.k and the band velocity v.sub.0(k)=∇.sub.kε.sub.k satisfy the conditions ε.sub.k=ε.sub.−k and v.sub.0(k)=−v.sub.0(−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.

[0081] 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 ψ.sub.k. For and only for noncentrosymmetric crystals, the transition rate from k to k′ and from −k to −k′ can be different: w.sub.k,k≠w.sub.−k′,″k. A scattering process with such asymmetry is referred to as skew scattering. Skew scattering is a source of the second-order DC response, i.e., rectification. The strength of skew scattering is characterized by w.sub.k′k.sup.(A)=(w.sub.k′,k″w.sub.−k′,k)/2, whereas the symmetric component w.sub.k′k.sup.(S)=(w.sub.k′k+w.sub.−k′,−k)/2 does not rely on inversion breaking. The presence of time-reversal symmetry guarantees w.sub.k′k=w.sub.−k,−k′, which allows us to write w.sub.k′k.sup.(S,A)=(w.sub.k′k±w.sub.kk′)/2.

[0082] The transition rate w.sub.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 w.sub.k′k.sup.(S)=2π|V.sub.k′k|.sup.2δ(ε.sub.k′″ε.sub.k), where V.sub.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 k⇄k′ and an indirect process k.fwdarw.k′ via an intermediate state q:

w.sub.k′k.sup.(A)=−(2π).sup.2∫.sub.qlmV.sub.k′qV.sub.qkV.sub.k′δ(ε.sub.k−ε.sub.k′)δ(ε.sub.k′−ε.sub.q),   (3)

which satisfies w.sub.k′k.sup.(A)=−w.sub.kk′.sup.(A). Here we use the notation ∫.sub.q=∫d.sup.dq/(2π).sup.d (d is the spatial dimension).

[0083] Skew scattering provides an extrinsic contribution to the anomalous Hall effect and the spin Hall effect. These are linear response phenomena in systems with broken time-reversal and spin-rotational symmetries, respectively. Skew scattering exists and contributes to the 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(r) is inversion asymmetric, such as a dipole or octupole potential. Even for a symmetric single-impurity potential V(r)=V(−r), such as a delta function potential, skew scattering can also appear. In such cases, skew scattering is due to the chirality of the electron wavefunction in a noncentrosymmetric crystal. In the simplest case of a delta function impurity potential, we have V.sub.k′k∝u.sub.k′|u.sub.k where |u.sub.k is the Bloch wavefunction within a unit cell. Then the expression for w.sup.(A) given in Eq. (3) involves the wavefunction overlap at different k points on the Fermi surface, similar to the Berry phase formula, and also skew scattering is known to be related to Berry curvature in a certain limit.

[0084] While there is no net chirality in nonmagnetic systems, the Bloch wavefunction |u.sub.k of an electron at finite momentum can 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 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 and the viscosity results from scattering with impurities; see FIGS. 5A and 6A.

[0085] To obtain the second-order conductivity, the distribution function ƒ(k) should be calculated to second order in the external electric field E(ω). Moreover, since skew scattering w.sup.(A) is parametrically smaller than w.sup.(S), we expand the distribution function up to first order in w.sup.(A). Solving the Boltzmann equation self-consistently at each order of E and w.sup.(A) yields the second-order DC current response j=−e∫.sub.kv.sub.0(k)ƒ(k) due to skew scattering in time-reversal-invariant systems. The second-order conductivity χ.sub.abc can be formally expressed in terms of various scattering times τ.sub.n, τ′.sub.n:

χ.sub.abc=−e.sup.3[τ′.sub.2τ′.sub.1(ω)τ.sub.1(ω)∫.sub.kv.sub.0,a∂.sub.k.sub.b∫.sub.k, w.sub.kk′.sup.(A)∂.sub.k′.sub.cƒ.sub.0(k′)+τ′.sub.2τ′.sub.2τ.sub.1(ω)∫.sub.kv.sub.0,a∫.sub.k, w.sub.kk′.sup.(A)∂.sub.k′.sub.b∂.sub.k′.sub.cƒ.sub.0(k′)]+(b⇄c)*,    (4)

where τ.sub.n(ω) is a shorthand for τ.sub.n(ω)=τ.sub.n/(1−iωτ.sub.n) and likewise for τ′.sub.n(ω). The scattering times τ.sub.n, τ′.sub.n with n=1, 2 are associated with the dominant symmetric scattering w.sup.(S); see below for details. They determine the steady-state distribution function up to second order in the electric field, when skew scattering is neglected. A detailed discussion of a self-consistent solution of the Boltzmann equation and Joule heating appears below.

[0086] A rough order-of-magnitude estimate of x can be obtained from Eq. (4) by using two scattering times τ and {tilde over (τ)} for the symmetric and skew scattering, respectively. Then, the second-order DC current of a metal or degenerate semiconductor under a linearly-polarized electric field is given by

[00002]$\begin{array}{cc}{j}_2\text{∼}en{v}_F.Math.{\left(\frac{eE\Delta t}{p_F}\right)}^2.Math.\frac{\tau }{\tilde{\tau }}& \left(5\right)\end{array}$

with t=τ for ωτ<<1 and t=1/ω for ωτ>>1. Here v.sub.F is the Fermi velocity, p.sub.F=ℏk.sub.F is the Fermi momentum, and n is the electron density.

[0087] This estimate provides a heuristic understanding of the second-order response. In the low-frequency limit ω.fwdarw.0, the ratio of the second-order current j.sub.2 and the linear response current j.sub.1=σ.sub.ωE is a product of two dimensionless quantities eEτ/p.sub.F and τ/{tilde over (τ)}. 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 τ/{tilde over (τ)} characterizes the strength of skew scattering which is responsible for second-order electrical response, relative to the symmetric scattering.

[0088] The short-circuit current responsivity .sub.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.sup.2, the current responsivity is given by .sub.I≡j.sub.2L/(j.sub.1EL.sup.2). Similarly, the voltage responsivity .sub.V is defined for the generated DC voltage as .sub.V≡(j.sub.2L/σ.sub.0)/(j.sub.1EL.sup.2). Both .sub.I and .sub.V are independent of the magnitude of the external electric field. Using the linear response conductivity σ.sub.ω˜env.sub.F. (eλ.sub.Fτ/ℏ)Re[(1−iωτ).sup.−1] obtained from the same Boltzmann equation gives the short-circuit current responsivity and voltage responsivity:

[00003]$\begin{array}{cc}{\mathrm{.Math.}}_I=\frac{1}{L}\frac{\chi }{{\sigma }_{\omega }}=\frac{{\eta }_I}{L}\frac{\tau }{\tilde{\tau }},{\eta }_I\text{∼}\frac{e\tau }{p_F},{\mathrm{.Math.}}_V=\frac{1}{L}\frac{\chi }{{\sigma }_{\omega }{\sigma }_0}=\frac{{\eta }_V}{L}\frac{\tau }{\tilde{\tau }},{\eta }_V\text{∼}\frac{1}{{\mathrm{env}}_F},& \left(6\right)\end{array}$

[0089] Here we define the reduced current and voltage responsivities η.sub.I and η.sub.V, respectively, which are independent of the sample size and the ratio τ/{tilde over (τ)}. The approximate relations for η.sub.I and η.sub.V hold in both low- and high-frequency limits. We also consider the ratio of the generated DC power and input power η.sub.P≡(j.sub.2.sup.2/σ.sub.0)/(σ.sub.107 E.sup.2), which characterizes the power conversion efficiency (assuming the load resistance and internal resistance are comparable). We find

[00004]$\begin{array}{cc}{\eta }_P=\frac{{\chi }^2{E}^2}{{\sigma }_0{\sigma }_{\omega }}\propto {\left(\frac{eE\Delta t}{p_F}\right)}^2.Math.{\left(\frac{\tau }{\tilde{\tau }}\right)}^2.& \left(7\right)\end{array}$

These figures of merit show that a material with low carrier density n or small Fermi momentum p.sub.F is desirable for efficient rectification. Moreover, a longer mean free time leads to a higher current responsivity. 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

[0090] Pristine graphene has high mobility and low carrier density, but it is centrosymmetric. Nevertheless, as a van der Waals material, we can assemble multilayer stacks to break inversion. Realizations of noncentrosymmetric structures include monolayer graphene on a hexagonal boron nitride (hBN) substrate, electrically-biased bilayer graphene, and multilayer graphene, such as ABA-stacked trilayer graphene. Anisotropy of energy dispersions, also used for second-order response, naturally arises from trigonal crystal structures.

[0091] 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.Math.p Hamiltonian considered here is

[00005]$\begin{array}{cc}{H}_s\left(k\right)=\left(\begin{array}{cc}\Delta & {\mathrm{svk}}_{-s}-\lambda k_s^2\\ {\mathrm{svk}}_s-\lambda k_{-s}^2& -\Delta \end{array}\right),& \left(8\right)\end{array}$

where s=±1 denotes two valleys at K and K′, respectively, and we define k.sub.±=k.sub.x±ik.sub.y. This Hamiltonian has a band gap of 2 in the energy spectrum and describes two graphene systems: (1) monolayer graphene on hBN, where the spatially varying atomic registries break the carbon sublattice symmetry and opens up gaps at Dirac points, and (2) bilayer graphene with an out-of-plane electric field that breaks the layer symmetry and opens up gaps. The difference between the monolayer and bilayer cases lies in the relative strength of v and λ. For the monolayer (bilayer) case, λk.sub.F.sup.2<<vk.sub.F(vk.sub.F<<k.sub.F.sup.2) is responsible for the trigonal warping of the linear (quadratic) energy dispersion.

[0092] The Bloch wavefunction of H.sub.S(k) has two components associated with sublattice (layer) degrees of freedom in monolayer (bilayer) graphene. When the band gap is present (≠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′ due to time reversal symmetry. Because of the three-fold rotation symmetry of graphene, Berry curvature dipole vanishes, leaving skew scattering as the only mechanism for rectification.

[0093] We calculate the second-order conductivity from Eq. (4), and find nonvanishing elements χ.sub.xxy=χ.sub.xyx=χ.sub.yxx=−χ.sub.yyy≡χ, consistent with the point group symmetry. The relation among the electric field, induced current, and underlying crystalline lattice is depicted in FIGS. 5B and 6B, along with the schematic pictures of the Fermi surface displacement in FIGS. 5C, 5D, 6C, and 6D. The function χ has the form χ=e.sup.3v.sub.F(τ.sup.3/{tilde over (τ)})ζ(ε.sub.F, ω), with the dimensionless function ζ given below. The ratio τ/{tilde over (τ)} is proportional to the gap and the Fermi surface trigonal warping. The gap describes the effect of inversion breaking on the electronic structure, while the Fermi surface trigonal warping is responsible for Fermi surface asymmetry within a valley E.sub.s(k)≠E.sub.s(−k).

[0094] FIGS. 7A-7D show the second-order conductivity χ and the reduced voltage responsivity η.sub.V for both monolayer and bilayer gapped graphene as functions of frequency v=ω/(2π). FIGS. 7A and 7C show the second-order conductivity χ for monolayer graphene and bilayer graphene, respectively. FIGS. 7B and 7D show the reduced voltage responsivity η.sub.V for monolayer graphene and bilayer graphene, respectively. The carrier density n is changed from 0.5×10.sup.12 cm.sup.−2 to 5×10.sup.12 cm.sup.−2. We use the values v=0.94×10.sup.6 m/s, =15 meV, and τ=1.13 ps for the monolayer case; =50 meV and τ=0.96 ps for the bilayer case, from transport, infrared spectroscopy, and scanning tunneling microscopy/spectroscopy measurements. For the bilayer case λ=(2m).sup.−1 is determined by the effective mass m≈0.033m.sub.e (m.sub.e: electron mass) and v≈1×10.sup.5 m/s.

[0095] 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λ.sub.F.sup.2=constant. The second-order conductivity then behaves as χ˜e.sup.3v.sub.Fτ.sup.2/ℏ.sup.2.Math.(τ/{tilde over (τ)}) for ωx<<1 and χ˜e.sup.3v.sub.Fl(ℏω).sup.2.Math.(τ/{tilde over (τ)}) for ωτ>>1.

[0096] 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 ω.sup.−2, the responsivity saturates in the high frequency limit. With a sample size of 10 μm, carrier density n=0.5×10.sup.12 cm.sup.−2 , and ratio τ/{tilde over (τ)}=0.01, reasonable estimates of mean free time and impurity concentration give a saturated voltage responsivity .sub.V˜30 V/W for monolayer graphene and 40 V/W for bilayer graphene (see below for the frequency dependence). As expected from the qualitative estimate Eq. (6), high responsivity is realized with low carrier density. The semiclassical approach is valid when the mean free path is much longer than the Fermi wavelength.

[0097] Both monolayer and bilayer graphene show broadband response. The saturation of responsivities lasts up to the onset frequency of an inter-band transition (not included in our study). For example, at the lowest density of n=0.5×10.sup.12 cm.sup.−2 shown in FIGS. 7A-7D, the Fermi energy for the monolayer case is ε.sub.F˜80 meV, so that due to the Pauli blocking the inter-band transition does not occur at frequencies below 2ε.sub.F˜40 THz. The inter-band threshold frequency also increases with the band gap, which is tunable in bilayer graphene by the displacement field.

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

[0099] 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, but with much larger band gaps. Hence, they should exhibit large intra-band second-order response as well. Topological insulator surface states also show this effect (see below). In addition to 2D systems, three-dimensional bulk materials without inversion are also worth considering. As discussed, they should have large responsivities in low carrier densities. From this respect, inversion-breaking Weyl semimetals could also be used in place of graphene for current rectification.

Discussion

[0100] So far, we have considered the skew scattering contribution to the second-order DC response by neglecting the effects from the Berry curvature Ω directly. It contributes to the distribution function through the collision integral and the electron's velocity as anomalous and side-jump velocities. Hence the second-order conductivity is also modified. The Berry curvature is odd (even) under time reversal (inversion), and thus is allowed to exist in time-reversal-invariant noncentrosymmetric materials. The effect of Berry curvature to linear transport properties has been well studied in the context of the anomalous Hall effect. Finite Berry curvature enforces a coordinate shift δr.sub.kk′ after scattering, giving a displacement of the center of an electron wave packet. It is accompanied by a potential energy shift δε.sub.k′k in the presence of the external field, which results in the collision integral

[ƒ]=∫.sub.k, [w.sub.k′kƒ(k,ε.sub.k)−w.sub.kk′ƒ(k′,ε.sub.k+δε.sub.k′k)].   (9)

[0101] One may decompose the distribution function in nonequilibrium as ƒ=ƒ.sub.0+ƒ.sup.scatt+ƒ.sup.adist, where ƒ.sub.0 is the equilibrium distribution function, ƒ.sup.scatt describes the effect of scattering to the distribution (with δε.sub.k′k=0), and ƒ.sup.adist is the anomalous distribution due to finite δε.sub.k′k. Because of the different origins, these distributions have distinct dependence on the scattering time τ. For low frequencies with ωτ<<1, the τ dependence of each term can be estimated by a simple power counting. Noting that the scattering rate w amounts to τ.sup.−1, we can expand the distribution function in powers of the electric field E, with coefficients depending on τ:

ƒ.sup.scatt=Σ.sub.n≥1ƒ.sub.n.sup.scatt, ƒ.sub.n.sup.scatt∝τ.sup.nE.sup.n,

ƒ.sup.adist=Σ.sub.n≥1ƒ.sub.n.sup.adist, ƒ.sub.n.sup.adist∝τ.sup.n−1E.sup.n,   (10),

in the weak impurity limit.

[0102] The electric current density is obtained by j=−e∫.sub.kv(k)ƒ(k), where v(k) is the electron's velocity,

[00006]$\begin{array}{cc}v=\frac{\partial {ϵ}_k}{\partial k}-\stackrel{.}{k}\times \Omega +{\int }_{k^{\prime }}{w}_{k^{\prime }k}\delta r_{k^{\prime }k}\equiv v_0+v_{\mathrm{av}}+v_{\mathrm{sj}}.& \left(11\right)\end{array}$

[0103] Besides the velocity given by the band dispersion v.sub.0, v contains two additional terms: The anomalous velocity v.sub.av is associated with the Berry curvature and the last term v.sub.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 (ƒ.sup.scatt+ƒ.sup.adist) and the electron velocity (v.sub.0+v.sub.av+v.sub.sj).

[0104] The dominant contribution in the weak impurity limit (τ.fwdarw.∞) emerges from the skew scattering with v.sub.0ƒ.sub.2.sup.scatt, resulting in χ∝τ.sup.2. The Berry curvature dipole contribution corresponds to v.sub.vaƒ.sub.1.sup.scatt, which amounts to χ∝τ. It shows the frequency dependence as χ′∝e.sup.3τD, ψ″≠0 for low frequencies (ωτ<<1) and χ′∝e.sup.3D/(ω.sup.2τ), χ″∝e.sup.3D/ω for high frequencies (ωτ>>1), where D is the Berry curvature dipole. The Berry curvature dipole is a quantity with dimensions of length in a two-dimensional system; it is dimensionless in a three-dimensional system. The Berry curvature dipole and skew scattering have the same frequency dependence although they are distinct in the x dependence.

[0105] Inversion breaking results in finite skew scattering w.sup.(A). Along with the accompanying anisotropic Fermi surface, the second-order electronic response induced by skew scattering should persist in any noncentrosymmetric material. 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.

[0106] In summary, the analysis above is a systematic study of second-order electrical response due to an intra-band process. Skew scattering makes the dominant contribution in the weak impurity limit. Skew scattering is responsible for current rectification and opens a new way to low-power energy harvesters and terahertz detection.

Semiclassical Boltzmann Theory

[0107] We describe in detail the analysis of the semiclassical Boltzmann transport theory for the calculation of second-order response in noncentrosymmetric materials. The following analysis is not meant to be limiting and to some extent is parallel with that for the anomalous Hall effect. However, a fundamental difference is found in symmetry; for a second-order response, inversion is broken whereas time-reversal symmetry is broken for the anomalous Hall effect.

[0108] The total Hamiltonian is given by

H=H.sub.0+V+U,

where H.sub.0 is the Hamiltonian for electrons in the clean limit with the eigenvalue ε.sub.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: ∇.sub.rU=−F. In the following, we consider an external electric field of frequency ω, namely, F=−eE(ω).

[0109] A formally rigorous analysis involves 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 may have 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 k can be designated. However, the uncertainty principle imposes the limitation rk≥1. We desire position and momentum resolutions set by the mean free path and the Fermi wavevector k.sub.F; i.e., k.sub.F>>1 is assumed for the semiclassical analysis. In the following, we deal with electron scattering by impurities. Since and k.sub.F depend on the impurity and electron densities, respectively, the condition is satisfied when the impurity density is sufficiently lower than the electron density.

[0110] We deal with spatially homogeneous systems and hence the electron distribution function ƒ(k, ε) does not depend on the spatial position, but on the wavevector k and the energy ε. We assume that the frequency of the external field ω is lower than the inter-band spacing, so that inter-band transitions are suppressed and negligible. We calculate the distribution function in the presence of the external field by the semiclassical Boltzmann equation (ℏ=1)

[00007]$\begin{array}{cc}\frac{\partial f}{\partial t}+\stackrel{.}{k}.Math.\frac{\partial f}{\partial k}=-\mathrm{��}\left[f\right].& \left(12\right)\end{array}$

The time derivative of the momentum is equal to the force felt by an electron wave packet: {dot over (k)}=F. The collision integral C[ƒ] is given by

[ƒ]=∫.sub.k′[w.sub.k′kƒ(k, ε)−w.sub.kk′ƒ(k′,ε+δε.sub.k′k)],

where w.sub.k′k is the scattering rate from a state with momentum k to one with k′. In the following, we consider elastic scattering. Then, the scattering rate is defined by

w.sub.k′k=2π|T.sub.k′k|δ(ε.sub.k′−ε.sub.k),

with the scattering T matrix T.sub.k′k. It is given by T.sub.k′k=k′|V|ψ.sub.k, where |k is the eigenstate of H.sub.0, the Hamiltonian in the clean limit, and |ψ.sub.k is the eigenstate of H.sub.0+V, including the scattering. |ψ.sub.k can be obtained as the solution to the Lippman-Schwinger equation |ψ.sub.k=|k+(ε.sub.k−H.sub.0+iδ).sup.−1V|ψ.sub.k.

[0111] The energy shift δÅ.sub.k′k occurs after a scattering owing to the Berry curvature Ω. The Berry curvature is defined from the Berry connection A as

A.sub.a=iu.sub.k∂.sub.k.sub.a|u.sub.k, Ω.sub.a=ε.sub.abc∂.sub.k.sub.bA.sub.c

where |ψ.sub.k is the lattice periodic part of the Bloch function at k and ε.sub.abc is the Levi-Civita symbol. The Berry curvature transforms under inversion and time reversal as

:Ω.sub.a(k).fwdarw.Ω.sub.a(−k), :Ω.sub.a(k).fwdarw.−ψ.sub.a(−k).

[0112] It is deduced that the Berry curvature vanishes everywhere in the Brillouin zone when symmetry exists. Finite Berry curvature causes the coordinate shift at a scattering process, denoted by δr.sub.k′k, which describes the displacement of the center of a wave packet by a scattering k.fwdarw.k′. The energy shift δ∈.sub.k′k is given by using the coordinate shift δr.sub.k′k as

δ∈.sub.k′k=−F.Math.δr.sub.k′k.

For a weak scattering with a small momentum change, the coordinate shift can be approximated as δr.sub.k′k≈(k′−k)×Ω(k). The coordinate shift is independent of the scattering time τ even though it is related to scattering; it describes the displacement after a single impurity scattering.

[0113] We decompose the distribution function ƒ as

ƒ=ƒf.sub.0+ƒ.sup.scatt+ƒ.sup.adist.

The first term ƒ.sub.0 is the distribution function in equilibrium, i.e., the Fermi-Dirac distribution. The second term ƒ.sup.scatt describes the scattering contribution without the Berry curvature and the last term ƒ.sup.adist is the anomalous distribution due to the Berry curvature and the energy shift δ∈.sub.k′k.

Scattering Rate in Noncentrosymmetric Media

[0114] In noncentrosymmetric media, the probabilities of a scattering process k.fwdarw.k′ and the inverted process −k.fwdarw.−k′ can be different. For time-reversal systems, this imbalance is captured by decomposing the scattering rate into two parts:

w.sub.kk′=w.sub.kk′.sup.(S)+w.sub.kk′.sup.(A),

where the symmetric and antisymmetric parts w.sup.(S) and w.sup.(A), respectively, are defined by

w.sub.kk′.sup.(S)=1/2(w.sub.kk′+w.sub.k′k), w.sub.kk′.sup.(A)=1/2(w.sub.kk′+w.sub.k′k).

When time reversal is preserved, the scattering rate satisfies the reversibility w.sub.kk′=w.sub.−k′,−k, which leads to the relations

w.sub.k′k.sup.(S)=w.sub.−k,−′.sup.(S)=w.sub.kk′.sup.(S),

w.sub.k′k.sup.(A)=w.sub.−k,−′.sup.(A)=−w.sub.kk′.sup.(A).

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

∫.sub.k′w.sub.kk′.sup.(A)=0.   (13)

When the elastic scattering is due to impurities, w.sup.(S) and w.sup.(A) are obtained to the lowest order in the impurity scattering potential by

w.sub.k′k.sup.(S)=2π|V.sub.k′k|.sup.2δ(ε.sub.k′−ε.sub.k),

w.sub.k′k.sup.(A)=−(2π).sup.2∫.sub.qlmV.sub.k′qV.sub.qkV.sub.kk′δ(ε.sub.k−ε.sub.k′)δ(ε.sub.k′−ε.sub.q),   (14, 15)

where denotes the average over the impurity distribution and V.sub.k′k is the matrix element for the single impurity scattering. 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.

[0115] V. I. Belinicher et al., “The photogalvanic effect in media lacking a center of symmetry,” Sov. Phys. Usp. 23, 199-223 (1980), 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.sub.k′k.sup.(A), similarly to our analysis. A difference can be found in the origin of finite w.sub.k′k. The asymmetry may be imposed on scattering potentials for impurity scattering, ionization, photoexcitation, and recombination. 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.sub.k′k.sup.(A), but even an isotropic scattering potential can generate w.sub.k′k.sup.(A) through a wavefunction.

Formal Solutions

[0116] The semiclassical Boltzmann equation (12) should be solved self-consistently to obtain the distribution function ƒ(k). Here, we expand the distribution function ƒ=ƒ.sub.0+ƒ.sup.scatt+ƒ.sup.adist with respect to the electric field E(ω) and the asymmetric part of the scattering rate w.sup.(A):

[00008]$f^{\mathrm{scatt}}\left(k\right)=\underset{n\ge 1}{.Math.}{f}_n^{scatt}\left(k\right),f_n^{scatt}\left(k\right)=\underset{m\ge 0}{.Math.}{f}_n^{\left(m\right)}\left(k\right),f^{\mathrm{adist}}\left(k\right)=\underset{n\ge 1}{.Math.}{f}_n^{adist}\left(k\right),f_n^{adist}\left(k\right)=\underset{m\ge 0}{.Math.}{g}_n^{\left(m\right)}\left(k\right).$

The subscript n and the superscript m correspond to the orders of the electric field E and the antisymmetric scattering rate w.sup.(A), respectively. As we have seen in Eqs. (14) and (15), w.sup.(A) is smaller than w.sup.(S), so that we treat the former as a perturbation. Then, the Boltzmann equation is decomposed for each ƒ.sub.n.sup.(m) or g.sub.n.sup.(m) to become

[00009]$\frac{\partial f_1^{\left(0\right)}}{\partial t}-\mathrm{eE}.Math.\frac{\partial f_0}{\partial k}=-{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}\left[f_1^{\left(0\right)}\left(k\right)-f_1^{\left(0\right)}\left(k^{\prime }\right)\right],\frac{\partial f_1^{\left(1\right)}}{\partial t}=-{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}\left[f_1^{\left(1\right)}\left(k\right)-f_1^{\left(1\right)}\left(k^{\prime }\right)\right]+{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(A\right)}{f}_1^{\left(0\right)}\left(k^{\prime }\right),\frac{\partial f_2^{\left(0\right)}}{\partial t}-\mathrm{eE}.Math.\frac{\partial f_1^{\left(0\right)}}{\partial k}=-{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}\left[f_2^{\left(0\right)}\left(k\right)-f_2^{\left(0\right)}\left(k^{\prime }\right)\right],\frac{\partial f_2^{\left(1\right)}}{\partial t}-\mathrm{eE}.Math.\frac{\partial f_1^{\left(1\right)}}{\partial k}=-{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}\left[f_2^{\left(1\right)}\left(k\right)-f_2^{\left(1\right)}\left(k^{\prime }\right)\right]+{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(A\right)}{f}_2^{\left(0\right)}\left(k^{\prime }\right),\frac{\partial g_1^{\left(0\right)}}{\partial t}=-{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}\left[g_1^{\left(0\right)}\left(k\right)-g_1^{\left(0\right)}\left(k^{\prime }\right)\right]+{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}{f}_0^{\prime }\left(k^{\prime }\right)\left(\mathrm{eE}.Math.\delta r_{k\prime k}\right),\frac{\partial g_1^{\left(1\right)}}{\partial t}=-{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}\left[g_1^{\left(1\right)}\left(k\right)-g_1^{\left(1\right)}\left(k^{\prime }\right)\right]+{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(A\right)}{g}_1^{\left(0\right)}\left(k^{\prime }\right)+{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(A\right)}{f}_0^{\prime }\left(k^{\prime }\right)\left(\mathrm{eE}.Math.\delta r_{k\prime k}\right),\frac{\partial g_2^{\left(0\right)}}{\partial t}-\mathrm{eE}.Math.\frac{\partial g_1^{\left(0\right)}}{\partial k}=-{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}\left[g_2^{\left(0\right)}\left(k\right)-g_2^{\left(0\right)}\left(k^{\prime }\right)\right]+{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}\left[f_1^{\prime \left(0\right)}\left(k^{\prime }\right)+g_1^{\prime \left(1\right)}\left(k^{\prime }\right)\right]\left(\mathrm{eE}.Math.\delta r_{k^{\prime }k}\right)+\frac{1}{2}{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}{f}_0^{″}\left(k^{\prime }\right){\left(\mathrm{eE}.Math.\delta r_{k\prime k}\right)}^2,\frac{\partial g_2^{\left(1\right)}}{\partial t}-\mathrm{eE}.Math.\frac{\partial g_1^{\left(1\right)}}{\partial k}=-{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}\left[g_2^{\left(1\right)}\left(k\right)-g_2^{\left(1\right)}\left(k^{\prime }\right)\right]+{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(S\right)}\left[f_1^{\prime \left(1\right)}\left(k^{\prime }\right)+g_1^{\prime \left(1\right)}\right]\left(\mathrm{eE}.Math.\delta r_{k\prime k}\right)+{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(A\right)}{g}_2^{\left(0\right)}\left(k^{\prime }\right)+{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(A\right)}\left[f_1^{\prime \left(0\right)}\left(k^{\prime }\right)+g_1^{\prime \left(0\right)}\right]\left(\mathrm{eE}.Math.\delta r_{k\prime k}\right)+\frac{1}{2}{\int }_{k\prime }{w}_{\mathrm{kk}\prime }^{\left(A\right)}{f}_0^{″}\left(k^{\prime }\right){\left(\mathrm{eE}.Math.\delta r_{k\prime k}\right)}^2.$

A prime symbol ′ added to the distribution functions stands for the derivative with respect to the energy: ƒ′(k, ε)=∂.sub.εƒ(k, ε).

[0117] Each equation should be solved self-consistently; however, this is usually not easy. Here, we solve a part of the collision integral involving w.sup.(S) self-consistently, while the other part with w.sup.(A) is treated as a perturbation. The collision integrals with w.sup.(S) define the scattering times τ.sub.n.sup.(m) and τ.sub.n.sup.(m) as follows:

[00010]$\begin{array}{cc}{\int }_{k^{\prime }}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}\left[f_n^{\left(m\right)}\left(k\right)-f_n^{\left(m\right)}\left(k^{\prime }\right)\right]=\frac{1}{{\tau }_n^{\left(m\right)}}{f}_n^{\left(m\right)}\left(k\right),& \left(16\right)\\ {\int }_{k^{\prime }}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}\left[g_n^{\left(m\right)}\left(k\right)-g_n^{\left(m\right)}\left(k^{\prime }\right)\right]=\frac{1}{t_n^{\prime \left(m\right)}}{g}_n^{\left(m\right)}\left(k\right).& \end{array}$

If τ.sub.n.sup.(m) and τ′.sub.n.sup.(m) do not satisfy the self-consistency, one may regard the equations above as the definitions for the relaxation time approximation. The notation of the scattering times τ.sub.n and τ′.sub.n used in the main text; they correspond to τ.sub.n≡τ.sub.n.sup.(0) and τ′.sub.n≡τ.sub.n.sup.(1).

[0118] We assume τ and τ′ are independent of momentum k in solving the equations. Then, we obtain the formal solutions for the distribution functions ƒ.sub.n.sup.(m) and g.sub.n.sup.(m). The second-order response contains different frequencies, 0 and ±2ω. Since now we focus on rectification, i.e., zero-frequency response, we write down only the second-order solutions with ω=0. Now, we list below the formal solutions up to the second order in the electric field (n≤2) and to the first order in w.sup.(A) (m≤1):

[00011]$\begin{array}{cc}{f}_1^{\left(0\right)}\left(k,\omega \right)={\tau }_{1\omega }^{\left(0\right)}{\mathrm{eE}}_a{\partial }_{k_a}{f}_0\left(k\right),& \left(17\right)\\ f_1^{\left(1\right)}\left(k,\omega \right)={\tau }_{1\omega }^{\left(1\right)}{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{\tau }_{1\omega }^{\left(0\right)}{\mathrm{eE}}_a{\partial }_{k_a^{\prime }}{f}_0\left(k^{\prime }\right),& \\ f_2^{\left(0\right)}\left(k,0\right)={\tau }_2^{\left(0\right)}{\mathrm{eE}}_a^{*}{\partial }_{k_a}{f}_1^{\left(0\right)}\left(k,\omega \right)+c.c.={\tau }_2^{\left(0\right)}{\mathrm{eE}}_a^{*}{\partial }_{k_a}{\tau }_{1\omega }^{\left(0\right)}{\mathrm{eE}}_b{\partial }_{k_b}{f}_0\left(k\right)+c.c.,& \\ f_2^{\left(1\right)}\left(k,0\right)={\tau }_2^{\left(1\right)}{\mathrm{eE}}_a^{*}{\partial }_{k_a}{f}_1^{\left(1\right)}\left(k,\omega \right)+c.c.+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{f}_2^{\left(0\right)}\left(k^{\prime },0\right)={\tau }_2^{\left(1\right)}{\mathrm{eE}}_a^{*}{\partial }_{k_a}{\tau }_{1\omega }^{\left(1\right)}{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{\tau }_{1\omega }^{\left(0\right)}{\mathrm{eE}}_b{\partial }_{k_b^{\prime }}{f}_0\left(k^{\prime }\right)+{\tau }_2^{\left(1\right)}{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{\tau }_2^{\left(0\right)}{\mathrm{eE}}_a^{*}{\partial }_{k_a^{\prime }}{\tau }_{1\omega }^{\left(0\right)}{\mathrm{eE}}_b{\partial }_{k_b^{\prime }}{f}_0\left(k^{\prime }\right)+c.c.,& \\ g_1^{\left(0\right)}\left(k,\omega \right)={\tau }_{1\omega }^{\prime \left(0\right)}{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}{f}_0^{\prime }\left(k^{\prime }\right){\mathrm{eE}}_a\delta r_{a;k^{\prime }k},& \left(18\right)\\ g_1^{\left(1\right)}\left(k,\omega \right)={\tau }_{1\omega }^{\prime \left(1\right)}\left[{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{\tau }_{1\omega }^{\prime \left(0\right)}{\int }_{k^{″}}^{}{w}_{k^{\prime }{k}^{″}}^{\left(S\right)}{f}_0^{\prime }\left(k^{″}\right){\mathrm{eE}}_a\delta r_{a;k^{″}{k}^{\prime }}+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{f}_0^{\prime }\left(k^{\prime }\right){\mathrm{eE}}_a\delta r_{a;k^{\prime }k}\right],& \\ g_2^{\left(0\right)}\left(k,0\right)={\tau }_2^{\prime \left(0\right)}\left[{\mathrm{eE}}_a^{*}{\partial }_{k_a}{g}_1^{\left(0\right)}\left(k,\omega \right)+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}\left[f_1^{\prime \left(0\right)}{k}^{\prime },\omega \right)+g_1^{\prime \left(0\right)}\left(k^{\prime },\omega \right)\right]{\mathrm{eE}}_a^{*}\delta r_{a;k^{\prime }k}+\frac{1}{2}{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}{f}_0^{″}\left(k^{\prime }\right){\mathrm{eE}}_a^{*}\delta r_{a;k^{\prime }k}{\mathrm{eE}}_b\delta r_{b;k^{\prime }k}]+c.c.={\tau }_2^{\prime \left(0\right)}\left[{\mathrm{eE}}_a^{*}{\partial }_{k_a}{\tau }_{1\omega }^{\prime \left(0\right)}{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}{f}_0^{\prime }\left(k^{\prime }\right){\mathrm{eE}}_{b;k^{\prime }k}+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}{\mathrm{eE}}_a^{*}\delta r_{a;k^{\prime }k}{\partial }_{{ϵ}_{k^{\prime }}}{\tau }_{1\omega }^{\left(0\right)}{\mathrm{eE}}_b{\partial }_{k_b^{\prime }}{f}_0\left(k^{\prime }\right)+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}{\mathrm{eE}}_a^{*}\delta r_{a;k^{\prime }k}{\partial }_{{ϵ}_{k^{\prime }}}{\tau }_{1\omega }^{\prime \left(0\right)}{\int }_{k^{″}}^{}{w}_{k^{\prime }{k}^{″}}^{\left(S\right)}{f}_0^{\prime }\left(k^{″}\right){\mathrm{eE}}_b\delta r_{b;k^{″}{k}^{\prime }}+\frac{1}{2}{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}{f}_0^{″}\left(k^{\prime }\right){\mathrm{eE}}_a^{*}\delta r_{a;k^{\prime }k}{\mathrm{eE}}_b\delta r_{b;k^{\prime }k}\right]+c.c.,& \\ g_2^{\left(1\right)}\left(k,0\right)={\tau }_2^{\prime \left(1\right)}\left[{\mathrm{eE}}_a^{*}{\partial }_{k_a}{g}_1^{\left(1\right)}\left(k,\omega \right)+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}\left[f_1^{\prime \left(1\right)}\left(k^{\prime },\omega \right)+g_1^{\prime \left(1\right)}\left(k^{\prime },\omega \right)\right]{\mathrm{eE}}_a^{*}\delta r_{a;k^{\prime }k}+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{g}_2^{\left(0\right)}\left(k^{\prime },0\right)+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}\left[f_1^{\prime \left(0\right)}\left(k^{\prime },\omega \right)+g_1^{\prime \left(0\right)}\left(k^{\prime },\omega \right)\right]{\mathrm{eE}}_a^{*}\delta r_{a;k^{\prime }k}+\frac{1}{2}{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{f}_0^{″}\left(k^{\prime }\right){\mathrm{eE}}_a^{*}\delta r_{a;k^{\prime }k}{\mathrm{eE}}_b\delta r_{b;k^{\prime }k}\right]+c.c.={\tau }_2^{\prime \left(1\right)}\left\{{\mathrm{eE}}_a^{*}{\partial }_{k_a}{\tau }_{1\omega }^{\prime \left(1\right)}\left[{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{\tau }_{1\omega }^{\prime \left(0\right)}{\int }_{k^{″}}^{}{w}_{k^{\prime }{k}^{″}}^{\left(S\right)}{f}_0^{\prime }\left(k^{″}\right){\mathrm{eE}}_b\delta r_{b;k^{″}{k}^{\prime }}+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{f}_0^{\prime }\left(k^{\prime }\right){\mathrm{eE}}_b\delta r_{b;k^{\prime }k}\right]+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(S\right)}{\mathrm{eE}}_a^{*}\delta r_{a;k^{\prime }k}{\partial }_{{ϵ}_{k^{\prime }}}\left[{\tau }_{1\omega }^{\left(1\right)}{\int }_{k^{″}}^{}{w}_{k^{\prime }{k}^{″}}^{\left(A\right)}{\tau }_{1\omega }^{\left(0\right)}{\mathrm{eE}}_b{\partial }_{k_b^{″}}{f}_0\left(k^{″}\right)+{\tau }_{1\omega }^{\prime \left(1\right)}{\int }_{k^{″}}^{}{w}_{k^{\prime }{k}^{″}}^{\left(A\right)}{\tau }_{1\omega }^{\prime \left(0\right)}{\int }_{k^{\mathrm{\prime \prime \prime }}}^{}{w}_{k^{″}{k}^{\mathrm{\prime \prime \prime }}}^{\left(S\right)}{f}_0^{\prime }{k}^{\mathrm{\prime \prime \prime }}\right)\hspace{1em}{\mathrm{eE}}_b\delta r_{b;k^{\mathrm{\prime \prime \prime }}{k}^{″}}+{\tau }_{1\omega }^{\prime \left(1\right)}{\int }_{k^{″}}^{}{w}_{k^{\prime }{k}^{″}}^{\left(A\right)}{f}_0^{\prime }\left(k^{″}\right){\mathrm{eE}}_b\delta r_{b;k^{″}{k}^{\prime }}]+{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{\tau }_2^{\prime \left(0\right)}\left[{\mathrm{eE}}_a^{*}{\partial }_{k_a^{\prime }}{\tau }_{1\omega }^{\prime \left(0\right)}{\int }_{k^{″}}^{}{w}_{k^{\prime }{k}^{″}}^{\left(S\right)}{f}_0^{\prime }\left(k^{″}\right){\mathrm{eE}}_b\delta r_{b;k^{″}{k}^{\prime }}+{\int }_{k^{″}}^{}{w}_{k^{\prime }{k}^{″}}^{\left(S\right)}{\mathrm{eE}}_a^{*}\delta r_{a;k^{″}{k}^{\prime }}{\partial }_{{ϵ}_{k^{″}}}{\tau }_{1\omega }^{\left(0\right)}{\mathrm{eE}}_b{\partial }_{k_b^{″}}{f}_0\left(k^{″}\right)+{\int }_{k^{″}}^{}{w}_{k^{\prime }{k}^{″}}^{\left(S\right)}{\mathrm{eE}}_a^{*}\delta r_{a;k^{″}{k}^{\prime }}{\partial }_{{ϵ}_{k^{″}}}{\tau }_{1\omega }^{\prime \left(0\right)}{\int }_{k^{\mathrm{\prime \prime \prime }}}{w}_{k^{″}{k}^{\mathrm{\prime \prime \prime }}}^{\left(S\right)}{f}_0^{\prime }\left(k^{\mathrm{\prime \prime \prime }}\right){\mathrm{eE}}_b\delta r_{b;k^{\mathrm{\prime \prime \prime }}{k}^{″}}+\frac{1}{2}{\int }_{k^{″}}{w}_{k^{\prime }{k}^{″}}^{\left(S\right)}{f}_0^{″}\left(k^{″}\right){\mathrm{eE}}_a^{*}\delta r_{a;k^{″}{k}^{\prime }}{\mathrm{eE}}_b\delta r_{b;k^{″}{k}^{\prime }}\right]+{\int }_{k^{\prime }}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{\mathrm{eE}}_a^{*}\delta r_{a;k^{\prime }k}\left[{\tau }_{1\omega }^{\left(0\right)}{\mathrm{eE}}_b{\partial }_{k_a}{f}_0^{\prime }\left(k\right)+{\partial }_{{ϵ}_{k^{\prime }}}{\tau }_{1\omega }^{\prime \left(0\right)}{\int }_{k^{″}}^{}{w}_{k^{\prime }{k}^{″}}^{\left(S\right)}{f}_0^{\prime }\left(k^{″}\right){\mathrm{eE}}_b\delta r_{b;k^{″}{k}^{\prime }}\right]+\frac{1}{2}{\int }_{k^{\prime }}^{}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{f}_0^{″}\left(k^{\prime }\right){\mathrm{eE}}_a^{*}\delta r_{b;k^{\prime }k}{\mathrm{eE}}_b\delta r_{b;k^{\prime }k}\right\}+c.c.& \end{array}$

τ.sub.nω.sup.(m) and τ′.sub.mω.sup.(m) are the shorthand notations for

[00012]${\tau }_{n\omega }^{\left(m\right)}=\frac{{\tau }_n^{\left(m\right)}}{1-i\omega {\tau }_n^{\left(m\right)}},{\tau }_{n\omega }^{\prime \left(m\right)}=\frac{{\tau }_n^{\prime \left(m\right)}}{1-i\omega {\tau }_n^{\prime \left(m\right)}}.$

Those formal solutions find the scattering time dependence of the distribution functions for low frequencies (ωτ<<1):

ƒ.sub.n.sup.scatt∝τ.sup.nE.sup.n,

ƒ.sub.n.sup.adist∝(τ.sup.0E.sup.n, τ.sup.1E.sup.n, . . . , τ.sup.n−1E.sup.n).

Higher-order corrections for the anomalous distribution (n≥2) have a nested structure of ƒ.sub.n−1.sup.scatt and ƒ.sub.n−1.sup.adist. In the weak disorder limit of τ.fwdarw.∞, we obtain

ƒ.sub.n.sup.adist∝τ.sup.n−1E.sup.n.

We can define another scattering time f from the antisymmetric part w.sup.(A). It is roughly speaking given by

[00013]$\begin{array}{cc}{\int }_{k^{\prime }}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}f\left(k^{\prime }\right)\text{∼}\frac{1}{\tilde{\tau }}f\left(k\right).& \left(19\right)\end{array}$

[0119] The definition of {tilde over (τ)} is accompanied with the distribution function because the optical theorem for elastic scattering concludes Eq. (13).

Current Response

[0120] From the distribution function, the electric current response is obtained by

j=−e∫.sub.kv(k)ƒ(k).

v(k) is the electron's group velocity, given by

[00014]$\begin{array}{cc}\begin{array}{c}v=\frac{\partial {ϵ}_k}{\partial k}-\stackrel{.}{k}\times \Omega +{\int }_{k^{\prime }}^{}{w}_{k^{\prime }k}\delta r_{k^{\prime }k}\\ \equiv v_0+v_{av}+v_{sj}.\end{array}& \end{array}$

[0121] The group velocity from the energy band dispersion v.sub.0 is independent of the Berry curvature, while the anomalous velocity v.sub.av=−{dot over (k)}×Ω=eE×Ω and the side-jump velocity v.sub.sj=∫.sub.k′w.sub.k′kδr.sub.k′k are induced by the Berry curvature. In terms of the electric field E, v.sub.av is linear in E whereas v.sub.0 and v.sub.sj do not depend on E. Lastly, only the side-jump velocity depends on the scattering time: v.sub.sj∝τ.sup.−1. Since there is no current in equilibrium, the electric current is described by

j=−e∫.sub.k(v.sub.0+v.sub.av+v.sub.sj)(ƒ.sup.scatt+ƒ.sup.adist),   (20)

which allows us to decompose the current into contributions of different origins. The contribution from the anomalous distribution ƒ.sup.adist could be recognized as a part of the side-jump effect since it is originated both from the Berry curvature and scattering.

[0122] Before calculating the second-order conductivity, consider the problem of energy dissipation or Joule heating because of j.Math.E≠0. Previous works indicates 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 |E|.sup.2ƒ″.sub.0. Technically, such an excess term arises when we substitute a distribution function f.sub.n (n≥1) into the collision integral [ƒ], 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. 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. The isotropic terms do not contribute to current from the band velocity since ∫.sub.kv.sub.0ƒ″.sub.0=∫.sub.k∇.sub.kƒ′.sub.0=0.

[0123] We can see the six possible combinations in Eq. (20). For the second-order current response, the scattering time dependence of each contribution for low frequencies ωτ<<1 is

v.sub.0ƒ.sub.2.sup.scatt∝τ.sup.2, v.sub.avƒ.sub.1.sup.scatt∝τ.sup.1, v.sub.sjƒ.sub.2.sup.scatt∝τ.sup.1,

v.sub.0ƒ.sub.2.sup.adist∝(τ.sup.1,τ.sup.0), v.sub.avƒ.sub.1.sup.adist∝τ.sup.0, v.sub.sjƒ.sub.2.sup.adist∝(τ.sup.0, τ.sup.−1).

In the weak disorder limit of τ.fwdarw.∞, the first term, the skew scattering contribution, v.sub.0ƒ.sup.scatt predominates, which we focus on here. This is equivalent to neglecting effects related to Berry curvature. The second-order current response arising from this term is

j.sub.2=−e∫.sub.kv.sub.0(k)ƒ.sub.2.sup.scatt(k)=−e∫.sub.kv.sub.0(k)=−e∫.sub.kv.sub.0(k)[ƒ.sub.2.sup.(0)(k)+ƒ.sub.2.sup.(1)(k)].

We define the second-order conductivity by the relation

j.sub.2,a=χ.sub.abcE*.sub.bE.sub.c,

and we obtain from Eqs. (17) and (18)

χ.sub.abc=χ.sub.abc.sup.(0)+χ.sub.abc.sup.(1),

χ.sub.abc.sup.(0)=−e.sup.3τ.sub.2.sup.(0)τ.sub.1ω.sup.(0)∫.sub.kv.sub.0,a∂.sub.k.sub.b∂.sub.k.sub.cƒ.sub.0(k)+(b⇄c)*,

χ.sub.abc.sup.(1)=−e.sup.3[τ.sub.2.sup.(1)τ.sub.1ω.sup.(1)τ.sub.1ω.sup.(0)∫.sub.kv.sub.0,a∂.sub.k.sub.b∫.sub.k′w.sub.s,kk′.sup.(A)∂.sub.k′.sub.cƒ.sub.0(k′)+τ.sub.2.sup.(1)τ.sub.2.sup.(0)τ.sub.1ω.sup.(0)∫.sub.kv.sub.0,a∫.sub.k′w.sub.s,kk′.sup.(A)∂.sub.k′.sub.b∂.sub.k′.sub.cƒ.sub.0(k′)]+(b⇄c)*.   (21, 22)

χ.sub.abc.sup.(0) vanishes identically in time-reversal-invariant systems because of ε.sub.k=ε.sub.−k; it follows from v.sub.0(k)=−v.sub.0(−k), and we see that the integrand is odd in k for χ.sub.abc.sup.(0). Therefore, the second-order contribution from the skew-scattering contribution arises from χ.sub.abc.sup.(1), i.e., χ.sub.abc=χ.sub.abc.sup.(1).

[0124] To calculate the second-order response, it is useful to express Eqs. (21) and (22) in the following forms:

[00015]$\begin{array}{c}{\chi }_{\alpha bc}^{\left(0\right)}=-e^3{\tau }_2^{\left(0\right)}{\tau }_{1\omega }^{\left(0\right)}{\int }_k\left[-f_0^{\prime }\left(k\right)\right]v_{0,a}\left({\partial }_{k_b}{\partial }_{k_c}{ϵ}_k\right)+{\left(b.Math.c\right)}^{*}\\ =-2e^3\mathrm{Re}\left[{\tau }_2^{\left(0\right)}{\tau }_{1\omega }^{\left(0\right)}\right]{\int }_k\left[-f_0^{\prime }\left(k\right)\right]v_{0,a}\left({\partial }_{k_b}{\partial }_{k_c}{ϵ}_k\right),\end{array}$${\chi }_{\alpha bc}^{\left(1\right)}=-e^3\left[{\tau }_2^{\left(1\right)}{\tau }_{1\omega }^{\left(1\right)}{\tau }_{1\omega }^{\left(0\right)}{\int }_k\left[-f_0^{\prime }\left(k\right)\right]\left({\partial }_{k_a}{\partial }_{k_b}{ϵ}_k\right){\int }_{k^{\prime }}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{v}_{0,c}-{\tau }_2^{\left(1\right)}{\tau }_2^{\left(0\right)}{\tau }_{1\omega }^{\left(0\right)}{\int }_k\left[-f_0^{\prime }\left(k\right)\right]v_{0,a}{\int }_{k^{\prime }}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}\left({\partial }_{k_b}{\partial }_{k_c}{ϵ}_k\right)+{\tau }_2^{\left(1\right)}{\tau }_2^{\left(0\right)}{\tau }_{1\omega }^{\left(0\right)}{\int }_k\left[-f_0^{\prime }\left(k\right)\right]{\partial }_{k_a}{\int }_{k^{\prime }}{w}_{{\mathrm{kk}}^{\prime }}^{\left(A\right)}{v}_{0,b}{v}_{0,c}\right]+{\left(b.Math.c\right)}^{*}.$

[0125] We also show the expression for χ.sub.abc.sup.(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 ε.sub.F is replaced by the chemical potential μ.sub.chem with the carrier density fixed. More explicitly, from the following relation

−ƒ′.sub.0(k)=∫dε[−ƒ′.sub.0(ε;μ.sub.chem)]δ(ε.sub.k−ε.sub.F),

the second-order conductivity at finite temperature is obtained by

χ.sub.abc(μ.sub.chem, T)=∫dε[−ƒ′.sub.0(ε;μ.sub.chem)]χ.sub.abc(ε,T=0).

Linear Response

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

j.sub.1,a(ω)=σ.sub.abE.sub.b(ω)=−e∫.sub.kv.sub.0,aƒ.sub.1.sup.scatt,

leading to the linear conductivity

σ.sub.ab=σ.sub.ab.sup.(0)+σ.sub.ab.sup.(1),

σ.sub.ab.sup.(0)=−e.sup.2τ.sub.1ω.sup.(0)∫.sub.kv.sub.0,a∝.sub.k.sub.bƒ.sub.0(k)=e.sup.2τ.sub.1ω.sup.(0)∫.sub.kv.sub.0,av.sub.0,b[−ƒ′.sub.0(k)],

σ.sub.ab.sup.(1)=−e.sup.2τ.sub.1ω.sup.(1)τ.sub.1ω.sup.(0)∫.sub.kv.sub.0,a∫.sub.k′w.sub.kk′.sup.(A)∝.sub.k.sub.b′∫.sub.0(k′).   (23)

Graphene-Based Models with Trigonal Lattice Structures

[0127] FIGS. 8A-8C illustrate explicit models for evaluating the current response. We consider two graphene-based models, utilizing monolayer and bilayer graphene. FIG. 8A shows the lattice structure of the monolayer graphene along with the crystal axes. FIG. 8B shows the lattice structure of bilayer graphene. The shaded circles depict distinct sublattice sites with different potential energies, and the filled and empty circles correspond to atoms on different layers. Both lattices have reflection planes parallel to the y axis and its symmetry partners under the three-fold rotation. FIG. 8A shows Brillouin zone and energy contours for the monolayer model. With small doping, there are trigonally-warped Fermi surfaces around K and K′ points, facing the opposite directions. The bilayer model exhibits the similar Fermi surfaces as the monolayer ones.

[0128] The crystalline lattices for the monolayer and bilayer models belong to a crystallographic point group D.sub.3h and C.sub.3v, respectively. Inversion is equivalent to the in-plane two-fold rotation in 2D models, which is absent in the present models. Around K and K′ points,

[00016]$K=\left(\frac{4\pi }{3\sqrt{3}},0\right)\mathrm{and}{K}^{\prime }=\left(-\frac{4\pi }{3\sqrt{3}},0\right),$

the k.Math.p Hamiltonian to second order is

[00017]$\begin{array}{cc}\begin{array}{c}{H}_s\left(k\right)=\left(\begin{array}{cc}\Delta & sv{k}_{-s}-\lambda k_s^2\\ sv{k}_s-\lambda k_{-s}^2& -\Delta \end{array}\right)\\ =\left[\mathrm{sv}{k}_x-\lambda \left(k_x^2-k_y^2\right)\right]{\sigma }_x+\left(v{k}_y+2s\lambda k_x{k}_y\right){\sigma }_y+\Delta {\sigma }_Z,\end{array}& \left(24\right)\end{array}$

where k.sub.x and k.sub.y are measured from K or K′ point, labeled by s=±1, and k.sub.s is defined by k.sub.s=k.sub.x+isk.sub.y. The energy dispersion of the positive branch is given by

ε.sub.k=√{square root over (v.sup.2k.sup.2+λ.sup.2k.sup.4+Δ.sup.2−2svλk.sup.3 cos 3θ)},

and the corresponding normalized wavefunction is

[00018]$.Math.k,s.Math.=\frac{1}{2\sqrt{{ϵ}_k\left({ϵ}_k+\Delta \right)}}\left\{{ϵ}_k+\Delta \left[sv{k}_x-\lambda \left(k_x^2-k_x^2\right)\right]+i\left(v{k}_y+2s\lambda k_x{k}_y\right)\right\}$

The Berry curvature is also calculated as

[00019]${\Omega }_z=-\frac{s\Delta \left(v^2-4{\lambda }^2{k}^4\right)}{2{ϵ}_k^3}.$

[0129] The Hamiltonian Eq. (24) is available to both monolayer and bilayer models, which are shown in FIGS. 8A and 8B, respectively. The Pauli matrix σ.sub.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 λ is dominant, which is explained below. Δ describes the potential energy difference at the two sublattice sites or layers, which opens a gap to the energy spectrum. We consider the two cases separately below.

Monolayer Graphene

[0130] For monolayer graphene, the terms with v dominate those with λ(vk.sub.F>>λk-hd F.sup.2), where λ characterizes the trigonal warping of the Fermi surface. s=+1 corresponds to K valley and s=−1 to K′ 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 λ are given by v=3ta/2 and λ=3ta.sup.2/8, respectively.

[0131] We assume that the trigonal warping is small and obtain the energy dispersion ε.sub.k to first order in λ as

[00020]${ϵ}_k=\sqrt{v^2{k}^2+{\Delta }^2}-sv\lambda \frac{k^3}{\sqrt{v^2{k}^2+{\Delta }^2}}\cos 3{\theta }_k+O\left({\lambda }^2\right),$

where the Fermi wavevector k.sub.F is given by

[00021]$k_F=\frac{\sqrt{{ϵ}_F^2-{\Delta }^2}}{v}.$

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

[00022]$D_0=\frac{{ϵ}_F}{2\pi v^2},$

and θ.sub.k denotes the polar angle of k. The Fermi surface can be parameterized by θ.sub.k as

[00023]$k\left({\theta }_k\right)=k_F+s\lambda \frac{k_F^2}{v}\cos 3{\theta }_k.$

[0132] We use the wavefunction on the Fermi surface to evaluate the scattering rate. The normalized wavefunction on the Fermi surface is given by

[00024]$s\lambda \frac{k_F^2}{v}\cos 3{\theta }_k.{\left(\sqrt{\frac{{ϵ}_F+\Delta }{2{ϵ}_F}},\frac{\left(s\cos {\theta }_k+i\sin {\theta }_k\right)\left(v{k}_F+i\lambda k_F^2\sin 3{\theta }_k\right)}{\sqrt{2{ϵ}_F\left({ϵ}_F+\Delta \right)}}\right)}^T$

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

[00025]${\Omega }_z=-\frac{s{v}^2\Delta }{2{ϵ}_F^3}+O\left({\lambda }^2\right).$

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

[00026]$V\left(r\right)=V_i\underset{j}{.Math.}\delta \left(r-r_j\right),$

but restrict scattering within each valley. V.sub.i is the impurity potential strength, n.sub.i is the impurity density, and r.sub.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 λ as

[00027]$\begin{array}{cc}{w}_{s,{\mathrm{kk}}^{\prime }}^{\left(S\right)}=\frac{\pi n_i{V}_i^2}{2{{ϵ}_F^2\left({ϵ}_F+\Delta \right)}^2}\left[{.Math.{\left({ϵ}_F+\Delta \right)}^2+v^2{k}_F^2{e}^{i\left(\theta -{\theta }^{\prime }\right)}.Math.}^2-2\mathrm{s\lambda }v{k_F^3\left({ϵ}_F+\Delta \right)}^2\sin \left(\theta -{\theta }^{\prime }\right)\left(\sin 3\theta -\sin 3{\theta }^{\prime }\right)\right]\delta \left({ϵ}_k-{ϵ}_{k^{\prime }}\right),w_{s,{\mathrm{kk}}^{\prime }}^{\left(A\right)}=-\frac{\pi n_i{V}_i^3\Delta k_F^2}{2{ϵ}_F^2}\left\{s\sin \left(\theta -{\theta }^{\prime }\right)+\lambda \frac{k_F}{\nu }\left[\sin \left(\theta -4{\theta }^{\prime }\right)+\sin \left(4\theta -{\theta }^{\prime }\right)-\sin \left(\theta +2{\theta }^{\prime }\right)+\sin \left(2\theta +{\theta }^{\prime }\right)\right]\right\}\delta ({ϵ}_k-{ϵ}_{k^{\prime }}),& \left(25\right)\end{array}$