Strain-Fluctuation-Induced Near-Quantization of Valley Hall Conductivity in Graphene Systems

TL;DR

This paper presents a theory explaining near-quantized valley Hall conductivity in graphene caused by strain-induced gauge potentials, highlighting the role of gauge disorder in anomalous Hall phenomena.

Contribution

The paper introduces a theoretical framework for valley Hall effect in graphene considering strain fluctuation-induced gauge potentials as the main disorder source.

Findings

Near-quantized valley Hall conductivity in the band transport regime.

Reproduction of temperature and gap dependence of the valley Hall effect.

Identification of gauge disorder as a key factor in anomalous Hall physics.

Abstract

We develop a theory of the valley Hall effect in high-quality graphene samples, in which strain fluctuation-induced random gauge potentials have been suggested as the dominant source of disorder. We find a near-quantized value of valley Hall conductivity in the band transport regime, which originates from an enhanced side jump of a Dirac electron when it scatters off the gauge potential. By assuming a small residue charge density our theory reproduces qualitatively the temperature- and gap-dependence of the observed valley Hall effect at the charge neutral point. Our study suggests that the valley Hall effect in graphene systems represents a new paradigm for the anomalous Hall physics where gauge disorder plays an important role.

Equations

$$H_0=\hslash v\mathbf{k}\cdot \mathbf{\sigma }+\Delta {\sigma }_z,$$

$$|u_{\mathbf{k}}^c⟩=(\genfrac{}{}{0ex}{}{\cos \frac{{\theta }_{\mathbf{k}}}{2}}{\sin \frac{{\theta }_{\mathbf{k}}}{2}{e}^{i{\varphi }_{\mathbf{k}}}}),|u_{\mathbf{k}}^v⟩=\left(\genfrac{}{}{0ex}{}{\sin \frac{{\theta }_{\mathbf{k}}}{2}{e}^{-i{\varphi }_{\mathbf{k}}}}{-\cos \frac{{\theta }_{\mathbf{k}}}{2}}\right),$$

$${\Omega }_c(\mathbf{k})=-{\Omega }_v(\mathbf{k})=-\frac{{\hslash }^2{v}^2\Delta }{2({\hslash }^2{v}^2{k}^2+{\Delta }^2{)}^{3/2}}.$$

$${\sigma }_H^{\text{int}}=4\frac{e^2}{\hslash }\sum _{\mathbf{k}}{\Omega }_c(\mathbf{k})\Theta ({\epsilon }_F-{\epsilon }_{\mathbf{k}})=\frac{2e^2}{h}(1-\cos {\theta }_F),$$

$$V_{\text{imp}}(\mathbf{r})=d_x(\mathbf{r}){\sigma }_x+d_y(\mathbf{r}){\sigma }_y.$$

$$d_{\pm }(\mathbf{q})=F(\mathbf{q})q^2{e}^{\mp 2i{\phi }_{\mathbf{q}}},$$

$$\delta {\mathbf{r}}_{{\mathbf{k}}^{\prime }\mathbf{k}}^c={\mathbf{A}}_{{\mathbf{k}}^{\prime }}^c-{\mathbf{A}}_{\mathbf{k}}^c-{\hat{\mathbf{D}}}_{{\mathbf{k}}^{\prime },\mathbf{k}}\arg (V_{{\mathbf{k}}^{\prime }\mathbf{k}}^c),$$

$$V_{{\mathbf{k}}^{\prime }\mathbf{k}}^c=\frac{1}{\mathcal{S}}\int d\mathbf{r}{e}^{-i\mathbf{q}\cdot \mathbf{r}}⟨u_{{\mathbf{k}}^{\prime }}^c|V_{\text{imp}}(\mathbf{r})|u_{\mathbf{k}}^c⟩,$$

$$⟨{\hat{\mathbf{D}}}_{{\mathbf{k}}^{\prime },\mathbf{k}}\arg (V_{{\mathbf{k}}^{\prime }\mathbf{k}}^c){⟩}_{\text{dis}}=\mathrm{Im}\frac{⟨V_{\mathbf{k}{\mathbf{k}}^{\prime }}^c{\hat{\mathbf{D}}}_{{\mathbf{k}}^{\prime },\mathbf{k}}{V}_{{\mathbf{k}}^{\prime }\mathbf{k}}^c{⟩}_{\text{dis}}}{⟨|V_{{\mathbf{k}}^{\prime }\mathbf{k}}^c{|}^2{⟩}_{\text{dis}}},$$

$$V_{{\mathbf{k}}^{\prime }\mathbf{k}}^c=d_{-}(\mathbf{q})⟨u_{{\mathbf{k}}^{\prime }}^c|{\sigma }_{+}|u_{\mathbf{k}}^c⟩+d_{+}(\mathbf{q})⟨u_{{\mathbf{k}}^{\prime }}^c|{\sigma }_{-}|u_{\mathbf{k}}^c⟩,$$

$$\begin{array}{rl}⟨|V_{{\mathbf{k}}^{\prime }\mathbf{k}}^c{|}^2{⟩}_{\text{dis}}& =\frac{1}{2}⟨d_{-}(\mathbf{q})d_{+}(\mathbf{q}){⟩}_{\text{dis}}{\sin }^2{\theta }_{\mathbf{k}},\end{array}$$

$$\begin{array}{rl}\delta {\mathbf{r}}_{{\mathbf{k}}^{\prime }\mathbf{k}}^c& =\frac{{\mathbf{\Omega }}_c(\mathbf{k})\times (\mathbf{k}-{\mathbf{k}}^{\prime })}{{\sin }^2{\theta }_{\mathbf{k}}}.\end{array}$$

$$\begin{array}{rl}\delta {\mathbf{r}}_{{\mathbf{k}}^{\prime }\mathbf{k}}^v& =\frac{{\mathbf{\Omega }}_v(\mathbf{k})\times (\mathbf{k}-{\mathbf{k}}^{\prime })}{{\sin }^2{\theta }_{\mathbf{k}}}.\end{array}$$

$$\begin{array}{rl}\delta {\mathbf{r}}_{{\mathbf{k}}^{\prime }\mathbf{k}}^n& =-\frac{{\mathbf{\Omega }}_n(\mathbf{k})}{|⟨u_{{\mathbf{k}}^{\prime }}^n|u_{\mathbf{k}}^n⟩{|}^2}\times (\mathbf{k}-{\mathbf{k}}^{\prime }),\end{array}$$

$$\begin{array}{rl}{\omega }_{{\mathbf{k}}^{\prime }\mathbf{k}}& =\frac{2\pi }{\hslash }⟨|V_{{\mathbf{k}}^{\prime }\mathbf{k}}^c{|}^2{⟩}_{\text{dis}}\delta ({ϵ}_{c,\mathbf{k}}-{ϵ}_{c,{\mathbf{k}}^{\prime }})\end{array}$$

$$\begin{array}{rl}& \frac{1}{{\tau }^{tr}}=\frac{2\pi }{\hslash }\sum _{{\mathbf{k}}^{\prime }}⟨|V_{{\mathbf{k}}^{\prime }\mathbf{k}}^c{|}^2{⟩}_{\text{dis}}(1-\cos ({\varphi }_{\mathbf{k}}-{\varphi }_{{\mathbf{k}}^{\prime }}))\delta ({ϵ}_F-{ϵ}_{c,{\mathbf{k}}^{\prime }}).\end{array}$$

$$\begin{array}{rl}& v_x^{sj}(\mathbf{k})=\sum _{{\mathbf{k}}^{\prime }}{\omega }_{{\mathbf{k}}^{\prime }\mathbf{k}}(\delta {\mathbf{r}}_{{\mathbf{k}}^{\prime }\mathbf{k}}^c{)}_x=\frac{\cos {\theta }_{\mathbf{k}}}{2k_F}\frac{1}{{\tau }^{tr}}\sin {\varphi }_{\mathbf{k}},\end{array}$$

$$g_{\mathbf{k}}=-\frac{\partial n_e}{\partial {ϵ}_{c,\mathbf{k}}}e{E}_y{v}_y^c(\mathbf{k}){\tau }^{tr},$$

$$\begin{array}{rl}{\sigma }_{xy}^{direct}& =e\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\frac{g_{\mathbf{k}}}{E_y}{v}_x^{sj}(\mathbf{k})=\frac{e^2}{4h}\cos {\theta }_F.\end{array}$$

$$\sum _{{\mathbf{k}}^{\prime }}{\omega }_{{\mathbf{k}}^{\prime }\mathbf{k}}(g_{\mathbf{k}}^{adist}-g_{{\mathbf{k}}^{\prime }}^{adist}+(-\frac{\partial n_e}{\partial {ϵ}_{c,\mathbf{k}}})e{E}_y(\delta {\mathbf{r}}_{{\mathbf{k}}^{\prime }\mathbf{k}}^c{)}_y)=0$$

$$\begin{array}{r}{\gamma }_{\mathbf{k}}=-(\frac{\partial n_e}{\partial {ϵ}_{c,\mathbf{k}}})e{E}_y\frac{\cos {\theta }_{\mathbf{k}}}{2k^2}.\end{array}$$

$$\begin{array}{rl}{\sigma }_{xy}^{adist}& =e\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\frac{g_{\mathbf{k}}^{adist}}{E_y}{v}_x^c(\mathbf{k})=\frac{e^2}{4h}\cos {\theta }_F.\end{array}$$

$$\begin{array}{rl}{\sigma }_{xy}^{sj}& =4({\sigma }_{xy}^{direct}+{\sigma }_{xy}^{adist})=\frac{2e^2}{h}\cos {\theta }_F,\end{array}$$

$$H_{\text{BLG}}=(\hslash v{)}^2[(k_x^2-k_y^2){\sigma }_x+2k_x{k}_y{\sigma }_y]+\Delta {\sigma }_z.$$

$$\delta n_0=\sum _{\mathbf{k}}[n_e({ϵ}_{\mathbf{k}},\mu ,T)-n_h({ϵ}_{\mathbf{k}},\mu ,T)],$$

$$\begin{array}{rl}& {\sigma }_{xy}^{direct}={\sigma }_{xy}^{adist}=\frac{e^2}{4h}\frac{\Delta }{k_{B}T}\\ & \times {\int }_{\Delta /k_{B}T}^{\infty }dx[\frac{e^{x-\frac{\mu }{k_{B}T}}}{(1+e^{x-\frac{\mu }{k_{B}T}}{)}^2}\frac{1}{x}+\frac{e^{x+\frac{\mu }{k_{B}T}}}{(1+e^{x+\frac{\mu }{k_{B}T}}{)}^2}\frac{1}{x}],\\ & {\sigma }_{xy}^{int}=\frac{e^2}{2h}\frac{\Delta }{k_{B}T}\\ & \times {\int }_{\Delta /k_{B}T}^{\infty }dx[\frac{1}{1+e^{x-\frac{\mu }{k_{B}T}}}\frac{1}{x^2}+\frac{1}{1+e^{x+\frac{\mu }{k_{B}T}}}\frac{1}{x^2}],\end{array}$$

$$\begin{array}{rl}& {\sigma }_H^v=4[{\sigma }_{xy}^{direct}+{\sigma }_{xy}^{adist}+{\sigma }_{xy}^{int}]\\ & =\frac{2e^2}{h}[\frac{1}{1+e^{(\Delta -\mu )/k_{B}T}}+\frac{1}{1+e^{(\Delta +\mu )/k_{B}T}}].\end{array}$$

$${\sigma }_H^v=\frac{2e^2}{h}[1-\frac{(4+{\sin }^2{\theta }_F)\cos {\theta }_F}{\frac{4}{w_0}\frac{{\tau }_{sc}}{{\tau }_a}+4-3{\sin }^2{\theta }_F}].$$

$$\begin{array}{rl}\mathbf{d}(\mathbf{q})& =g_1\frac{1}{|\mathbf{q}{|}^4}\mathcal{F}(\mathbf{q})(q_x^2-q_y^2,-2q_x{q}_y),\\ V(\mathbf{q})& =g_2\frac{1}{|\mathbf{q}{|}^2}\mathcal{F}(\mathbf{q}),\end{array}$$

$$\mathcal{F}(\mathbf{q})$$

Full text

I Introduction

Charge carriers in graphene can be described by the two-dimensional (2D) Dirac equation, which exhibit a slew of interesting electronic properties Castro Neto et al. (2009). One of the consequences is that strains behave as pseudo-magnetic fields for Dirac electrons, and carefully designed lattice deformation pattern can result in the formation of Landau levels Guinea et al. (2010). Even without the engineered strains, random strain fluctuations are inevitable in 2D materials Meyer et al. (2007); Fasolino et al. (2007); Vozmediano et al. (2010); Amorim et al. (2016). They appear in the form of either out-of-plane corrugations due to thermal ripples or in-plane displacement from the interaction with substrates. Strain fluctuations then act like random magnetic (gauge) field, and can significantly affect transport behaviors. Their effect has been extensively explored in the longitudinal transport phenomena Ludwig et al. (1994); Couto et al. (2014); Engels et al. (2014); Ochoa et al. (2013); Vicent et al. (2017), such as weak localization Couto et al. (2014); Engels et al. (2014) and spin relaxation phenomena Ochoa et al. (2013); Vicent et al. (2017).

In this work, we investigate the role of random strain fluctuations on a particular type of transverse transport phenomena—the valley Hall effect of 2D Dirac electrons Xiao et al. (2007). The valley Hall effect can be regarded as two opposite copies of the anomalous Hall effect of a pair of gapped Dirac points related by time-reversal symmetry Nagaosa et al. (2010). That is, carriers in the two valleys will flow in the opposite transverse direction upon the application of a longitudinal electric field. Our motivation is twofold. First, in gapped monolayer and bilayer graphene systems, the valley Hall effect has been observed experimentally Gorbachev et al. (2014); Sui et al. (2015); Shimazaki et al. (2015). Thus, detailed experimental study of Hall-type transport in the presence of strain fluctuations is feasible. Secondly, all recent valley Hall measurements in graphene are carried out in high-quality devices in which strain fluctuations are the dominant source of disorder Couto et al. (2014); Engels et al. (2014). However, no existing theories have discussed its effect on the Hall transport. We will show that, the strain fluctuations are essential to understand the valley Hall effect of 2D Dirac electrons, and provide a new insight to recent debates Lensky et al. (2015); Kirczenow (2015); Li et al. (2011); Zhu et al. (2017); Song and Vignale ; Brown et al. (2018) on the observed nonlocal signals at the charge neutral point Gorbachev et al. (2014); Sui et al. (2015); Shimazaki et al. (2015).

Our main results are summarized below. Focusing on the band transport regime, we find that the valley Hall conductivity exhibits a singular behavior in the presence of strain fluctuation-induced long-range gauge disorder: as the Fermi level sweeps across the band edge, it jumps from zero to a nearly quantized value, $2e^{2}/h$ for monolayer graphene and $4e^{2}/h$ for bilayer graphene. The origin of this singular behavior is traced back to an enhanced side jump of a Dirac electron when it scatters off the gauge potential Yang et al. (2011a, b). Furthermore, at the charge neutral point, by assuming a small residue charge density we calculate the temperature- and gap-dependence of the valley Hall conductivity, which qualitatively agrees with the experiment Sui et al. (2015). We also find that strain-induced long-range scalar potential can reduce the valley Hall conductivity from its quantized value. Our study suggests that the valley Hall effect in graphene systems represents a new paradigm for the anomalous Hall physics where gauge disorder plays an important role. Our theory can also be applied to other 2D valley Hall materials such as transition metal dichalcogenides as well Wu et al. (2019); Hung et al. (2019).

The paper is organized as follows. In Sec. II we present the intrinsic valley Hall effect. In Sec. III we study the random strain-induced side jump and the resulting valley Hall effect. Temperature dependence and rigorous numerical analysis are shown in Sec. IV and Sec. V, respectively. In Sec. VI, effect of long-range scalar potential is investigated. Finally, discussion and conclusion are made in Sec. VII. Technical details are relegated to the appendixes.

II Intrinsic valley Hall effect

We begin with the following effective Hamiltonian

[TABLE]

which describes the low-energy electron dynamics in one of the Dirac valleys in gapped graphene. Here $v$ is the velocity, $\bm{k}=(k_{x},k_{y})$ is the two-dimensional wave vector, $2\Delta$ is the band gap opened by inversion symmetry breaking, and $\bm{\sigma}$ represents the sublattice indices. The Hamiltonian for the other valley can be obtained by performing a time-reversal operation on $H_{0}$. The energy dispersion is given by $\varepsilon_{c,v}=\pm\varepsilon_{\bm{k}}=\pm(\hbar^{2}v^{2}k^{2}+\Delta^{2})^{1/2}$ with the corresponding eigenstates

[TABLE]

where the superscript $c$ and $v$ label the conduction and valence bands, respectively, and the angular variables $\theta_{\bm{k}}$ and $\phi_{\bm{k}}$ are defined as $\theta_{\bm{k}}\equiv\cos^{-1}(\Delta/\varepsilon_{\bm{k}})$ and $\phi_{\bm{k}}\equiv\tan^{-1}(k_{y}/k_{x})$. The opening of the band gap gives rise to nonzero Berry curvature Xiao et al. (2007), defined by $\bm{\Omega}_{n}(\bm{k})=\Omega_{n}(\bm{k})\hat{z}=i\langle\nabla_{\bm{k}}u_{\bm{k}}^{n}|\times|\nabla_{\bm{k}}u_{\bm{k}}^{n}\rangle$ Xiao et al. (2010). For the two-band model given in Eq. (1) we have

[TABLE]

The Berry curvature in the other valley has opposite sign, as required by time-reversal symmetry.

Assuming weak inter-valley scattering, we can decouple the valley Hall effect into two copies of the anomalous Hall effect for each valley species. However, this decoupling must be treated with care. In the anomalous Hall effect, there is an intrinsic contribution to the Hall conductivity, given by the summation of the Berry curvature over all occupied Bloch states Nagaosa et al. (2010). It can be divided into two parts. One comes from fully occupied bands. This part manifests as chiral edge states at the Fermi energy, and gives rise to a quantized contribution to the Hall conductivity. On the other hand, the valley Hall systems considered here are topologically trivial without protected edge states. Therefore no electronic transport is possible when the Fermi energy is inside the band gap 111Specifically, if one integrates the Berry curvature over the valence band of the gapped Dirac Hamiltonian, one would obtain a Hall conductivity of $\pm e^{2}/h$ for the two valleys (counting the spin degeneracy). Naively taking their difference seems to suggest that the valley Hall conductivity should be $2e^{2}/h$ when the Fermi energy is in the band gap. However, since there is no topological edge states, this result is inconsistent with the general notion that electronic transport should be a Fermi surface property and fully occupied topologically trivial bands should not contribute to electronic transport. The nonzero result is an artifact due to treating the two Dirac valleys separately, which cannot capture the global topology of the bands.. Consequently, we shall drop this contribution in the calculation of the valley Hall conductivity. This leaves us with the contribution from partially occupied bands. Haldane has argued that this contribution can be written as the Berry phase of quasiparticles moving on the Fermi surface, and thus it can be regarded as a Fermi surface property Haldane (2004). For an electron-doped sample, the Fermi surface contribution is

[TABLE]

where the factor of 4 counts the spin and valley degeneracy, $\varepsilon_{F}$ is the Fermi energy, and $\theta_{F}\equiv\theta_{k_{F}}$. When $\varepsilon_{F}$ approaches the band edge, $\sigma_{H}^{\text{int}}$ vanishes. If the sample is hole doped, then one should consider the Fermi surface contribution of holes, which is opposite to that of electrons. Obviously, the intrinsic contribution alone cannot explain the observed valley Hall effect around the charge neutral point Gorbachev et al. (2014); Sui et al. (2015); Shimazaki et al. (2015).

III Random strain-induced side jump

To remedy this situation, we consider the effect of random strain-induced gauge disorder. In graphene systems, strain can be induced either by out-of-plane corrugations or by in-plane displacements of the carbon atoms. Both will generate a random gauge potential for Dirac electrons Vozmediano et al. (2010),

[TABLE]

The point group symmetry of graphene requires that $\bm{d}(\bm{r})\propto(u_{xx}-u_{yy},-2u_{xy})$, where $u_{\alpha\beta}$ is the strain tensor defined in terms of the deformation field $\bm{u}(\bm{r})$, $u_{\alpha\beta}=(\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha})/2$. The Fourier transform of $\bm{d}(\bm{r})$ has the form Couto et al. (2014)

[TABLE]

where $d_{\pm}(\bm{q})=d_{x}(\bm{q})\pm id_{y}(\bm{q})$, and $F(\bm{q})$ is a prefactor depending on the details of the strain field. The appearance of the phase angle $2\varphi_{\bm{q}}$, defined by $\varphi_{\bm{q}}\equiv\tan^{-1}(q_{y}/q_{x})$, is due to the fact that $u_{\alpha\beta}$ is a second-order derivative of either the height field (out-of-plane corrugations) or random potentials from substrates (in-plane displacements) Couto et al. (2014); Vozmediano et al. (2010). Both modes are long-wavelength elastic modes, as indicated by the $\bm{q}$-dependence of $F(\bm{q})$.

To reveal the effect of long-range gauge disorder on the Hall conductivity, we invoke a recently developed semiclassical Boltzmann theory Sinitsyn et al. (2006, 2007). Such theory has been widely used to investigate Hall-type transport under scalar disorder in various systems, whereas our work for the first time generalizes it to the long-range gauge disorder. Without loss of generality, we consider electrons scattering in the conduction band. The key quantity here is a sudden coordinate shift Berger (1970) experienced by an electron wave packet as it scatters off an impurity (see Fig. 1 (a)), given by Sinitsyn et al. (2006)

[TABLE]

where $\bm{A}^{c}_{\bm{k}}=\langle u_{\bm{k}}^{c}|i\bm{\nabla}_{\bm{k}}u_{\bm{k}}^{c}\rangle=-\sin^{2}\frac{\theta_{\bm{k}}}{2}(\hat{z}\times\bm{k})/k^{2}$ is the Berry connection of the conduction band and $\hat{\bm{D}}_{\bm{k}^{\prime},\bm{k}}=\bm{\nabla}_{\bm{k}^{\prime}}+\bm{\nabla}_{\bm{k}}$. The information of the impurity is encoded in the quantity

[TABLE]

where $\cal S$ is the system area and $\bm{q}\equiv\bm{k}^{\prime}-\bm{k}$ is the momentum transfer of electrons. The overall effect of the disorder is obtained by taking the disorder average, under which the last term of Eq. (7) becomes

[TABLE]

where $\langle\dots\rangle_{\text{dis}}$ stands for disorder or thermal average.

To proceed further, it is convenient to write $V^{c}_{\bm{k}^{\prime}\bm{k}}$ using the chiral basis

[TABLE]

where $\sigma_{\pm}=(\sigma_{x}\pm i\sigma_{y})/2$. The correlation $\langle|V_{\bm{k}^{\prime}\bm{k}}^{c}|^{2}\rangle_{\text{dis}}$ then breaks up into terms with opposite chirality $\langle d_{\mp}(\bm{q})d_{\pm}(-\bm{q})\rangle_{\text{dis}}$, and terms with the same chirality $\langle d_{\pm}(\bm{q})d_{\pm}(-\bm{q})\rangle_{\text{dis}}$. The former contain no phase factor, whereas the latter, according to Eq. (6), is proportional to $\exp(\mp 4i\varphi_{\bm{q}})$. Note that such term should be small after angular average, and we can neglect the term. The accuracy of this approximation is demonstrated by rigorous numerical analysis in Appendix  B. By such approximation, the correlation reduces to

[TABLE]

where we have used the condition of elastic scattering, i.e., $\theta_{\bm{k}}=\theta_{\bm{k}^{\prime}}$.

This approximation also applies to the correlation $\langle V_{\bm{k}\bm{k}^{\prime}}^{c}\hat{\bm{D}}_{\bm{k}^{\prime},\bm{k}}V_{\bm{k}^{\prime}\bm{k}}^{c}\rangle_{\text{dis}}$. Making use of the fact that $\hat{\bm{D}}_{\bm{k}^{\prime},\bm{k}}d_{\pm}(\bm{q})=(\bm{\nabla}_{\bm{k}}+\bm{\nabla}_{\bm{k}^{\prime}})d_{\pm}(\bm{k}^{\prime}-\bm{k})=0$, we obtain

[TABLE]

Similar calculation can be applied to electrons from valence band:

[TABLE]

We can see that in the coordinate shift $\delta\bm{r}^{c}_{\bm{k}^{\prime}\bm{k}}$ or $\delta\bm{r}^{v}_{\bm{k}^{\prime}\bm{k}}$, the strain-related perfactor $F(\bm{q})$ drops out completely, thus this expression is generally applicable for both out-of-plane and in-plane modes of strain fluctuations. Physically, $\delta\bm{r}_{\bm{k}^{\prime}\bm{k}}$ describes a coordinate shift transverse to the momentum change $\bm{k}-\bm{k}^{\prime}$, leading to a Hall-like current.

It is useful to compare with the short-range scalar disorder Sinitsyn et al. (2007). In that case, the coordinate shift is given by

[TABLE]

where $n=c/v$ refers to conduction (valence) band. In the denominator $|\langle u^{n}_{\bm{k}^{\prime}}|u^{n}_{\bm{k}}\rangle|^{2}\approx 1$ near the band edge, whereas for gauge disorder in Eq. (12) $\sin^{2}\theta_{\bm{k}}=4|\langle u^{c}_{\bm{k}^{\prime}}|\sigma_{+}|u^{c}_{\bm{k}}\rangle|^{2}\approx 0$, which makes a significant difference. Physical meaning for this difference is that close to the band edge, pseudospin of electrons is almost fixed, and the probability of a spin-flipping transition driven by gauge disorder is vanishingly small (see Fig. 1 (b)).

Once the coordinate shift $\delta\bm{r}_{\bm{k}^{\prime}\bm{k}}$ is derived, one can calculate its contribution, known as the side jump, to the valley Hall conductivity. There are two different types of side-jump effects: the direct side-jump contribution $\sigma_{xy}^{direct}$ and the anomalous distribution-induced contribution $\sigma_{xy}^{adist}$ Sinitsyn et al. (2006, 2007). We can first write down the scattering rate

[TABLE]

and the transport time

[TABLE]

For monolayer or bilayer graphene, point group symmetry requires that random gauge potential follows $d_{\pm}(\bm{q})=F(\bm{q})q^{2}e^{-2i\varphi_{\bm{q}}}$, which means $\langle d_{-}(\bm{q})d_{+}(\bm{q})\rangle_{\text{dis}}$ becomes a function of $\bm{q}=\bm{k}-\bm{k}^{\prime}$. Therefore the transport time $\tau^{tr}$ is isotropic for all $\bm{k}$ on the Fermi surface.

The coordinate shift $\delta\bm{r}^{c}_{\bm{k}^{\prime}\bm{k}}$ leads to an average side-jump velocity $\bm{v}^{sj}(\bm{k})$

[TABLE]

where $k_{F}$ is the Fermi wave vector. In the presence of an external electric field $E_{y}$, a nonequilibrium correction to the distribution function is given by

[TABLE]

where $n_{e}=1/[\exp((\epsilon_{c,\bm{k}}-\epsilon_{F})/k_{B}T)+1]$ is the Fermi distribution function, and $v_{y}^{c}(\bm{k})=v\sin\theta_{\bm{k}}\sin\phi_{\bm{k}}$ is the bare velocity of electrons along the electric field. As a result, at $T=0$ K, the valley Hall conductivity $\sigma_{xy}^{direct}$ (for each valley and spin) reads

[TABLE]

The physical process can be understood as follows: in the weak disorder limit, the scattering rate $\omega_{\bm{k}^{\prime}\bm{k}}$ is tiny, which gives rise to a small anomalous velocity $v^{sj}_{x}(\bm{k})$ on average during the scattering events. On the other hand, weak disorder means long lifetime, i.e., electrons can be accelerated by an electric field for more time until they are stopped by disorder scattering. This creates a large correction to the Fermi distribution $g_{\bm{k}}$, that is, more electrons contribute to the transverse transport. As a result, a product of the small anomalous velocity and the large number of electrons lead to a disorder-independent valley Hall conductivity $\sigma_{xy}^{direct}$ as the leading-order term of the disorder potential.

In addition, $\delta\bm{r}^{c}_{\bm{k}^{\prime}\bm{k}}$ can cause an anomalous distribution $g^{adist}_{\bm{k}}$ that also contributes to the Hall current, i.e., $\sigma_{xy}^{adist}$ term. To find its form, let us solve the equation

[TABLE]

to derive the nonequilibrium distribution function $g^{adist}_{\bm{k}}$. Take the ansatz $g^{adist}_{\bm{k}}=\gamma_{\bm{k}}k_{x}$, we find

[TABLE]

Then at $T=0$ K, $\sigma_{xy}^{adist}$ is given by

[TABLE]

Finally, the total side-jump valley Hall conductivity $\sigma_{xy}^{sj}$ is given by

[TABLE]

where the factor of 4 counts the valley and spin degeneracy. We notice that $\sigma_{H}^{\text{sj}}$ reaches its maximum value at the band edge, then decreases gradually as the Fermi energy moves away (see Fig. 1 (c)). By Eq. (4) and (23), the total valley Hall conductivity then becomes $\sigma_{H}^{v}=\sigma_{H}^{\text{int}}+\sigma_{H}^{\text{sj}}=2e^{2}/h$, which is quantized. We have also obtained the same valley Hall conductivity by adopting a different diagrammatic approach in Appendix A. Intriguingly, this gives us the same result as obtained for short-range gauge disorder using the diagrammatic approach Yang et al. (2011a, b). This indicates that valley Hall conductivity when a Dirac electron meets with gauge disorder is actually a universal quantity, which is in striking contrast to the case of scalar disorder. Note that at zero temperature our finding proposes a new type of geometric quantization of the Fermi surface, in contrast to the well-known topological quantization of the Fermi sea. This is one of the main results of our paper.

We may also apply our theory to bilayer graphene, which can be modeled by

[TABLE]

After some algebra, we find $\sigma_{H}^{\text{sj}}=(4e^{2}/h)\cos\theta_{F}$ and the total valley Hall conductivity $\sigma_{H}^{v}=4e^{2}/h$. The doubling of $\sigma_{H}^{v}$ can be traced back to the phase winding number of $2$ of the electrons around the Dirac point in bilayer graphene. Note that we have neglected the weaker fluctuation of next-nearest-neighbor interlayer coupling in bilayer graphene Son et al. (2011).

IV Temperature dependence

So far we have studied the valley Hall effect at zero temperature and its dependence on the Fermi energy. Next we will focus on the charge neutral point and consider the dependence of the valley Hall effect on the band gap and temperature. We will consider bilayer graphene systems to compare with the experiment Sui et al. (2015). If the system is perfectly uniform (chemical potential $\mu=0$), $\sigma_{H}^{v}$ due to thermally-activated carriers will be exponentially small. For example, take $2\Delta\sim 100$ meV in a bilayer graphene, then $\sigma_{H}^{v}\sim 2\times 10^{-3}e^{2}/h$ at $T=70$ K, several orders of magnitude smaller than experimental values. However, in graphene systems the charge density typically fluctuates due to the formation of electron-hole puddles or gate voltage fluctuation Dean et al. (2010). By assuming a small residue charge density $\delta n_{0}$, we can move the chemical potential into the conduction or valence bands according to

[TABLE]

where $n_{e/h}(\epsilon_{\bm{k}},\mu,T)=1/[\exp(\pm(\epsilon_{c/v,\bm{k}}-\mu)/k_{B}T)+1]$ is the Fermi distribution function for electrons (holes) and $\mu$ is the chemical potential. At finite temperatures, both electrons and holes contribute to the transport. Including all these effects, we find the extrinsic and intrinsic valley Hall conductivity

[TABLE]

where the two terms in the square brackets originate from electrons and holes, respectively. For monolayer graphene, we obtain a much enhanced $\sigma_{H}^{v}$:

[TABLE]

For bilayer graphene, we will acquire an extra factor 2, as compared to monolayer case.

Figure 2 shows the calculated $\sigma_{H}^{v}$ as a function of temperature for $\delta n_{0}\sim 1.0\times 10^{10}$ cm*-2*, a value well below the density resolution of the nonlocal peak in the experiment Gorbachev et al. (2014). The downwards trend at medium temperature range $T\sim 60-80$ K reproduces precisely the temperature dependence in bilayer graphene reported in one of the experiments (see Supplementary Figure 9b of reference  Sui et al. (2015)). Another trend found in the experiment Sui et al. (2015), i.e., smaller $\sigma_{H}^{v}$ with larger gap, is also reproduced. At higher temperature we find that $\sigma_{H}^{v}$ will increase again, although there is no available experimental data to compare with. As far as we know, there are no counterparts of such peculiar temperature behaviors reported in anomalous Hall physics Nagaosa et al. (2010); Xiao et al. . The origin of such non-monotonic temperature behavior is the competition between intrinsic and side-jump contributions under thermal activation. The intrinsic contribution, as a summation of Berry curvature over occupied states, favors high-temperature regime with more occupied states, whereas the side-jump contribution, suffering from a suppression of thermally activated carriers by $1/k_{B}T$, favors low-temperature regime, as shown in Eq. (LABEL:temperature2). Nevertheless, neither a gap nor temperature dependence of valley Hall conductivity has been discussed by previous theoretical papers Lensky et al. (2015); Kirczenow (2015); Li et al. (2011); Zhu et al. (2017); Song and Vignale ; Brown et al. (2018), and our work represents a first step towards understanding these peculiar behaviors. Moreover, the newly predicted non-monotonic temperature dependence can be used as a test for our theory.

V Numerical analysis

A crucial approximation we made in the semiclassical approach is assuming that $\langle d_{\pm}(\bm{q})d_{\pm}(-\bm{q})\rangle_{\text{dis}}\propto\exp(\mp 4i\varphi_{\bm{q}})$ is negligible after angular average. The accuracy of this approximation needs to be examined since it introduces a long-range anisotropic correlation function. The effect of such correlation function is unclear yet, and an analytical treatment seems impractical. To this end, in the following part we develop a new numerical approach which treats the correlation function rigorously in Appendix B. We consider a particular type of random strain, i.e., out-of-plane corrugations with $\bm{d}(\bm{q})=g_{1}|\bm{q}|^{-4}\mathcal{F}(\bm{q})(q_{x}^{2}-q_{y}^{2},-2q_{x}q_{y})$, where $\langle\mathcal{F}(\bm{q})\mathcal{F}(-\bm{q})\rangle_{\text{dis}}$ $\propto|\bm{q}|^{2}$, and $g_{1}$ quantifies the electron-phonon coupling strength in graphene Couto et al. (2014). We find that $\sigma_{H}^{sj}$ is very close to $\frac{2me^{2}}{h}\cos\theta_{F}$ ($m=1$ for monolayer; $m=2$ for bilayer), differing by $1.69\%$ for monolayer graphene and $0.91\%$ for bilayer graphene. Such small deviation from the quantized values justifies our previous treatment of dropping the fast-oscillating part $\langle d_{\pm}(\bm{q})d_{\pm}(-\bm{q})\rangle_{\text{dis}}$.

VI Effect of long-range scalar potential

Besides the gauge potential $\bm{d}(\bm{q})$, a long-range scalar potential may also arise due to random strain fluctuations. For out-of-plane corrugations Couto et al. (2014), the scalar potential is given by $V(\bm{q})=g_{2}|\bm{q}|^{-2}\mathcal{F}(\bm{q})$, where $g_{2}$ is an electron-phonon coupling parameter. Different from the gauge potential $\bm{d}(\bm{q})$, which represents a pseudomagnetic field, $V(\bm{q})$ is electrostatic in nature and thus it needs to be screened by the static dielectric function $\epsilon(\bm{q})=1+V_{c}(\bm{q})N_{F}$ within the Thomas-Fermi approximation Ando (2006); Das Sarma et al. (2011), where $N_{F}$ is the density of states, $V_{c}(\bm{q})=2\pi e^{2}/\epsilon_{0}|\bm{q}|$ is the two-dimensional Coulomb potential, and $\epsilon_{0}$ is the background (including the substrate) dielectric constant. Here we will focus on the strong screening regime (see Appendix A.2): $\epsilon_{0}\ll e^{2}N_{F}/k_{F}$, which characterizes the charge neutral point of gapped graphene systems.

Taking into account both the gauge and scalar potentials, we find that the total valley Hall conductivity of a gapped graphene at $T=0$ K is (see Appendix A.3):

[TABLE]

Here $w_{0}=\sqrt{(k_{F}W/\pi)^{2}-1}/\pi$ is a cutoff-related factor with the device width $W$. For $W=1$ $\mu$m and residue charge density $\delta n_{0}=1.0\times 10^{10}$ cm*-2*, we estimate $w_{0}\approx 2.2$. The ratio $\tau_{sc}/\tau_{a}=(2e^{2}g_{1}/\hbar\epsilon_{0}g_{2})^{2}w_{0}$ defines the relative strength between the scalar and gauge disorder. In the limit $\tau_{sc}/\tau_{a}\rightarrow\infty$, i.e., in the absence of scalar disorder, $\sigma_{H}^{v}$ reduces to $2e^{2}/h$. Figure 3 shows the calculated $\sigma_{H}^{v}$ as a function of $\tau_{sc}/\tau_{a}$ at two different Fermi energies. From this figure, one can immediately understand why nonlocal signal is only measured in high-quality graphene on hBN rather than on SiO2 substrate zyb , since in the former (latter) case gauge (scalar) disorder is the dominant source of disorder scattering, corresponding to the limit $\tau_{sc}/\tau_{a}\rightarrow\infty$ ($\tau_{sc}/\tau_{a}\rightarrow 0$).

VII Conclusion and discussion

We have provided an alternative scenario to understand the large valley Hall conductivity observed in experiments, based on scattering from random strain fluctuations. The origin is intimately related to an enhanced coordinate shift under gauge disorder scattering in Dirac systems. Temperature and gap dependence is qualitatively reproduced. Our work paves the way for studying the effect of classical strain modes, or phonon modes Zhang and Niu (2015), on the transverse transport of broad classes of 2D materials and van der Waals heterostructures.

A few remarks are in order. Note that our theory is only valid for the band transport regime; the valley Hall effect in the phonon-assisted variable-range-hopping regime Zou and Zhu (2010) still remains an open question. In addition, we would like to point out that a direct comparison of our result to the experiments requires a careful extraction of the valley Hall conductivity from the nonlocal measurement Gorbachev et al. (2014); Sui et al. (2015); Shimazaki et al. (2015). In particular, an accurate determination of the valley diffusion length $\ell_{v}$ is crucial Sui et al. (2015) since the nonlocal signals depend on $\ell_{v}$ exponentially Abanin et al. (2009). Then the “smoking gun” validation of our theory, i.e., a non-monotonic temperature behaviors of valley Hall conductivity, can be examined. As a supplement, weak-localization magneto-resistance measurement Couto et al. (2014) and Raman spectroscopy Neumann et al. (2015) can also be used to uncover the role of strain fluctuations. A complete understanding of the valley Hall effect thus requires further experimental and theoretical efforts.

A universal transverse Imbert-Fedorov (IF) shift of electrons was also discovered in Weyl semimetals Jiang et al. (2015); Yang et al. (2015). The universal coordinate shift proposed in our work differs from the IF shift in two ways. First, the IF shift appears in the three-dimensional massless Weyl fermions, and has no counterpart in 2D, while the coordinate shift is relevant to the 2D massive Dirac fermions. Second, the IF shift occurs at normal interface, while the coordinate shift becomes universal only when scattered by gauge disorder. Appealingly, the considerations in this work can be generalized to describe other 2D systems with topological properties, such as superconductors Yu et al. (2018), excitons Onga et al. (2017), plasmons Shi and Song (2018), polaritons Gutiérrez-Rubio et al. (2018), or under magnetic field Komatsu et al. (2018).

Recently, experimental observations of valley Hall transport have also been made in atomically thin MoS2 systems Wu et al. (2019); Hung et al. (2019). The fact that monolayer and trilayer MoS2 share qualitatively similar behaviors of valley Hall signals Wu et al. (2019) indicates that the details of Berry curvature distribution in the conduction bands may have little influence on the final result. This is actually consistent with our theoretical prediction of the article.

Acknowledgements

We are grateful to Qian Niu, Shengyuan A. Yang, Mengqiao Sui, and Yuanbo Zhang for stimulating discussions. This work is supported by DOE BES Pro-QM EFRC (DE-SC0019443). W. Y. S. also acknowledges the support of a startup grant from Guangzhou University.

Appendix A Diagrammatic approach to mixed gauge and scalar disorder

In this section we present a full quantum mechanical treatment using the diagrammatic approach Sinitsyn et al. (2007); Shan et al. (2013) to study the valley Hall effect in the presence of out-of-plane corrugations Couto et al. (2014). The purpose of this section is two fold. The diagrammatic approach provides an additional check of the semiclassical result. In addition, this approach is systematic, and can be applied to mixed gauge and scalar disorder. It is also convenient for numerical calculations when we consider the anisotropic long-range correlation functions. In the following we will focus on monolayer graphene.

A.1 Correlation function

For out-of-plane corrugation the induced vector potential $\bm{d}(\bm{q})$ and scalar potential $V(\bm{q})$ are given by Couto et al. (2014)

[TABLE]

where

[TABLE]

$h(\bm{q})$ is the Fourier transform of the height field $h(\bm{r})$, and $g_{1,2}$ quantify the electron-phonon coupling strength in graphene. We assume the height correlation is $\langle h(\bm{q})h(-\bm{q})\rangle_{\text{dis}}\propto|\bm{q}|^{-4}$, from which one finds $\langle\mathcal{F}(\bm{q})\mathcal{F}(-\bm{q})\rangle_{\text{dis}}=C|\bm{q}|^{2}$, where $C$ is a material-dependent parameter. The Born scattering amplitude is given by

[TABLE]

where $d_{\pm}(\bm{q})=g_{1}q^{2}_{\mp}\mathcal{F}(\bm{q})/|\bm{q}|^{4}$ and $V(\bm{q})=g_{2}\mathcal{F}(\bm{q})/|\bm{q}|^{2}$. The identity $\bm{q}=\bm{k}^{\prime}-\bm{k}$ leads to a useful relation $q_{\pm}=k(e^{\pm i\phi_{\bm{k}^{\prime}}}-e^{\pm i\phi_{\bm{k}}})$. The correlation functions are

[TABLE]

for gauge and scalar disorder, respectively. Similar to what we did in the main text, we will ignore the first two terms in the square bracket on the right-hand side of Eq. (32) and (33). These terms originate from $\langle d_{\pm}(\bm{q})d_{\pm}(\bm{q})\rangle_{\text{dis}}$, whose contributions are vanishingly small and can be neglected. The validity of this approximation will be demonstrated in the next section. Under this approximation, the correlation function for gauge disorder becomes

[TABLE]

A.2 Relaxation time and longitudinal conductivity

We can use the correlation function in Eq. (33) and (34) to derive the relaxation time. According to Fermi’s golden rule, the relaxation time for gauge disorder reads

[TABLE]

where $N_{F}=\epsilon_{F}/2\pi\hbar^{2}v^{2}$ is the density of states per spin and valley. The cutoff factor $w_{0}$, given by

[TABLE]

is introduced to remove the divergence at small momentum Vozmediano et al. (2010), whose physical origin is due to the finite-size effect of samples. Consider a sample with width $W$, then $|\bm{q}|\geq q_{0}=\frac{\pi}{W}$. This corresponds to a cutoff angle $\phi_{0}=2\arcsin(q_{0}/2k_{F})$, and hence $w_{0}=\sqrt{(k_{F}W/\pi)^{2}-1}/\pi$.

On the other hand, for scalar potential, we need to take into account the screening effect. Based on the Thomas-Fermi approximation Ando (2006); Das Sarma et al. (2011), we can replace $U^{V}(\bm{q})$ by $U^{V}(\bm{q})/\epsilon(\bm{q})$, where $\epsilon(\bm{q})$ is the dielectric function, $\epsilon(\bm{q})=1+g_{s}g_{v}V_{c}(\bm{q})N_{F}$, $V_{c}(\bm{q})=2\pi e^{2}/\epsilon_{0}|\bm{q}|$ is the two-dimensional Coulomb potential, $\epsilon_{0}$ is the background (including the substrate) dielectric constant, and $g_{s}=2$ $(g_{v}=2)$ refers to the spin (valley) degeneracy. Now we can write down the relaxation time

[TABLE]

In the last step of above derivation, we have taken the strong-screening limit: $\epsilon_{0}\ll e^{2}N_{F}/k_{F}$. The physics around the charge neutral point belongs to this limit, within which the correlation effectively becomes a short-range one.

Next we need to figure out the modified velocity $(\tilde{v}_{x}^{\bm{k}})_{cc}$ due to the intraband vertex correction. According to Fig 4, we have

[TABLE]

and

[TABLE]

for gauge and scalar disorder, respectively. Here the bare velocity is given by $(v_{x}^{\bm{k}})_{cc}\equiv\langle u^{c}_{\bm{k}}|\hat{v}_{x}|u^{c}_{\bm{k}}\rangle=v\sin\theta_{\bm{k}}\cos\varphi_{\bm{k}}$, and retarded (advanced) Green’s function is $G_{c,\bm{k}}^{r/a}=1/(\epsilon_{F}-\epsilon_{c,\bm{k}}\pm i\hbar/2\tau_{j})$, with $j=A,V$. To solve the equations, we take the ansatz $(\tilde{v}_{x}^{\bm{k}})_{cc}=\eta(v_{x}^{\bm{k}})_{cc}$ and substitute it into Eqs. (38) and (39). We find $\eta=w_{0}$ for gauge disorder and $\eta=\frac{4(1-\frac{1}{2}\sin^{2}\theta_{\bm{k}})}{1+3\cos^{2}\theta_{\bm{k}}}$ for scalar disorder, respectively. Based on this, we can obtain the transport time

[TABLE]

which agree with previous results by using the Boltzmann approach Couto et al. (2014).

A.3 Valley Hall conductivity

For scalar disorder in the strong-screening limit $\epsilon_{0}\ll e^{2}N_{F}/k_{F}$, the correlation function $\langle U^{V}(\bm{q})U^{V}(-\bm{q})\rangle_{\text{dis}}/\epsilon^{2}(\bm{q})$ exhibits a short-range behavior, whose extrinsic valley Hall conductivity has been obtained before Sinitsyn et al. (2007); Yang et al. (2011a):

[TABLE]

for each spin and valley. In this expression, higher-order side-jump or skew-scattering contributions have been ignored since we are interested in the low doping regime. Next we turn to the gauge disorder, which has not been studied before. Fig. 4 (b)-(e) correspond to the leading-order side-jump contributions to the valley Hall conductivity in the weak scattering limit, given by

[TABLE]

[TABLE]

By making use of $[U^{i}_{nn^{\prime}}(\bm{q})]^{*}=U^{i}_{n^{\prime}n}(-\bm{q})$, $i=A,V$, $n,n^{\prime}=c,v$, we find the following symmetry properties:

[TABLE]

After some algebra, we find that for gauge disorder

[TABLE]

and thus

[TABLE]

Moreover, there are equivalent contributions from diagrams by rotating Fig. 4 (b)-(e) by $180^{\circ}$, then exchanging the subscript $x$, $y$ Sinitsyn et al. (2007). Therefore the total extrinsic valley Hall conductivity (including spin degeneracy) reads $\frac{2e^{2}}{h}\cos\theta_{\bm{k}}$, which reproduces Eq. (23) by applying the semiclassical approach.

Furthermore we can study the situation with mixed gauge and scalar disorder. In this case, the total relaxation time is given by

[TABLE]

For convenience, we can introduce a Fermi-energy-independent relaxation time $\tau_{a}$, $\tau_{sc}$ for gauge and scalar disorder, respectively:

[TABLE]

For such mixed disorder, we can follow the same procedure as above and derive

[TABLE]

where the correction factor

[TABLE]

Finally, by adding the intrinsic term, we find the total contribution is

[TABLE]

which gives Eq. (28) in the main text.

Appendix B Numerical treatment of anisotropic correlation function

In our derivation, we have made use of the approximation that terms $\langle d_{\pm}(\bm{q})d_{\pm}(\bm{q})\rangle_{\text{dis}}$ in the correlation function vanish after angular average. In this section, we test this approximation by treating the correlation function exactly, i.e., keeping all the terms in Eq. (32) and (33). We focus on a generic type of gauge disorder satisfying $\langle\mathcal{F}(\bm{q})\mathcal{F}(-\bm{q})\rangle_{\text{dis}}=C|\bm{q}|^{2\epsilon+2}$, where the value of $\epsilon$ depends on the microscopic details. For example, $\epsilon=0$ and $1$ corresponds to thermally excited and substrate-induced ripples, respectively Vozmediano et al. (2010). We also consider a generic chiral model

[TABLE]

where $m=1$ and $2$ correspond to monolayer and bilayer graphene, respectively. This leads to the eigenvalue and eigenstates

[TABLE]

where $\cos\theta_{\bm{k}}=\Delta/\epsilon_{c,\bm{k}}$, $\sin\theta_{\bm{k}}=Ak^{m}/\epsilon_{c,\bm{k}}$. The density of states is given by $N_{F}=\epsilon_{F}/(2\pi mA^{2}k_{F}^{2m-2})$. We can write down a complete form of the correlation function

[TABLE]

For monolayer graphene with thermal ripples: $m=1$, $\epsilon=0$, it reduces to Eq. (32). Similar to Eq. (35), we can evaluate the relaxation time by

[TABLE]

where a cutoff factor is introduced

[TABLE]

For such anisotropic problem, it is convenient to establish a self-consistent equation for the mean free path $(L_{x}^{\bm{k}})_{cc}$ Tokura (1998)

[TABLE]

Since the bare velocity follows $(v_{x}^{\bm{k}})_{cc}=\frac{mA}{\hbar}k^{m-1}\sin\theta_{\bm{k}}\cos\phi_{\bm{k}}$, we can take the following ansatz,

[TABLE]

By substituting it into Eq. (55), we find that the coefficients $f_{n}$ satisfy

[TABLE]

where the matrix $T$ has a non-closed form

[TABLE]

Fortunately, numerics indicates that the result converges really fast. Base on this observation, we find that the total valley Hall conductivity become

[TABLE]

To gain some insight, we consider two special cases that may be relevant to experiments.

B.1 $\epsilon=0$

First we consider $\epsilon=0$, i.e., out-of-plane corrugations defined by Eq. (29) and (30). By definition (54), we have an iterative relation

[TABLE]

Since $w_{1}=w_{0}-1$, we obtain $w_{m}=w_{0}-m$. Then the matrix $T$ can be simplified, and the valley Hall conductivity reads

[TABLE]

for $m=1$ (monolayer graphene) and

[TABLE]

for $m=2$ (bilayer graphene). Note that these results are very close to $\frac{2me^{2}}{h}\cos\theta_{\bm{k}}$, differing by less than $2\%$, confirming the validity of the approximation used in Sec. A and the main text.

B.2 $\epsilon=0.821$

In this section, we consider $\epsilon=0.821$, which is a more realistic value for thermally excited ripples Le Doussal and Radzihovsky (1992), i.e., out-of-plane corrugations. By definition (54), we have an iterative relation

[TABLE]

and

[TABLE]

This leads to a solution

[TABLE]

where

[TABLE]

By numerics, we find that

[TABLE]

for $m=1$ and

[TABLE]

for $m=2$. Again the results are very close to $\frac{2me^{2}}{h}\cos\theta_{\bm{k}}$, differing by less than $0.1\%$.