# Modeling satellite-based open water fraction via flexible Beta regression:
An application to wetlands in the north-western Pacific coast of Mexico

Inder Tecuapetla-Gómez, Julia Trinidad Reyes  I. Tecuapetla-Gómez, Dirección
de Cátedras, Consejo Nacional de Ciencia y Tecnología (CONACyT), Av.
Insurgentes Sur 1582, Crédito Constructor, Benito Juárez 03940, Ciudad de
México. I. Tecuapetla-Gómez, Dirección de Geomática, Comisión Nacional para el
Conocimiento y Uso de la Biodiversidad (CONABIO), Liga Periférico-Insurgentes
Sur 4903, Parques del Pedregal, Tlalpan 14010, Ciudad de México (e-mail:
itecuapetla@conabio.gob.mx). J. Trinidad Reyes, Licenciatura en Matemáticas
Aplicadas, Facultad de Ciencias, UNAM, Ciudad de México
(e-mail:july_treyes@ciencias.unam.mx).

###### Abstract

Carbon sequestration and water filtering are two examples of the several
ecosystem services provided by wetlands. Open water mapping is an effective
means to measure any wetland extension as these are comprised of many open
water bodies. An economical, though indirect, approach towards mapping open
water bodies is through applying geo-computational methods to satellite
images. In this work we propose the flexible Beta regression (FBR) model to
predict open water fraction from measurements of a water index. We focus on
observations derived from two MODIS images acquired during the dry season of
2008 in Marismas Nacionales, a wetland located in the north-western Pacific
coast of Mexico. A Bayesian estimation procedure is presented to estimate the
FBR model; in particular, we provide details of a nested Metropolis-Hastings
and Gibbs sampling algorithm to carry out parameter estimation. Our results
show that the FBR model produces valid predictors of water fraction unlike the
standard model. Our work is complemented by software developed in the R
language and available through a GitHub repository.

###### Index Terms:

Flexible Beta regression, open water fraction, Metropolis-Hastings, Gibbs
sampling, Satellite images, MODIS

##  I Introduction and Data

Wetlands are ecosystems which serve as feeding and refuge areas for fish and
crustaceans. Also, they function as biological water filters, act as natural
systems to control floods, provide a natural environment for the sequestration
and long-term storage of carbon dioxide from the atmosphere; they are the
largest natural source of methane to the atmosphere [1, 2].

In light of the above, modeling the dynamic of these ecosystems has gained a
lot of attention in recent years. An indirect and economical way to study the
dynamic of extensive wetlands is through modeling biophysical variables
derived from satellite images [3]. For instance, in order to determine the
seasonal dependent extension of a wetland, time series of moisture or water
indices can be generated and analyzed for the detection of open water.
Typically, these indices are obtained by combining different spectral bands of
satellite products such as the Moderate Resolution Imaging Spectroradiometer
(MODIS) or the Landsat mission.

Since the spatial resolution of a Landsat image (30 m) is higher than that of
many MODIS products (250, 500 or 1000 m), it is sensible to classify open
water in a Landsat image and then calculate the corresponding water fraction
in the MODIS image. For instance, if a MODIS pixel contains 25 Landsat pixels
and 12 of the latter are classified as open water, then the fraction of water
in the MODIS pixel is 12/25. The success of this seemingly straightforward
procedure depends on acquiring Landsat images with good quality data (e.g.,
cloudless and free of shadow clouds and artifacts). It is common that during
an entire dry season, only one Landsat image meets this quality requirements.
This precludes the generation of time series of water fraction of MODIS images
from previously classified Landsat products. Generating MODIS time series of
water indices, however, is a simple task. In [4] the authors found that among
14 moisture and water indices, the MODIS-based Modified Normalized Difference
Water Index $\mbox{B}_{6}$ (MNDWI6), see [5], has the best performance in
appropriately defining open water bodies based on moderate spatial resolution
products (such as Landsat).

Figure 1: Marismas Nacionales (in yellow). (Left) Pixels of water fraction
with good quality assessment (in blue). (Right) Pixels of MNDWI6 with good
quality assessment (in red).

It is of interest to assess whether the MNDWI6 is an appropriate predictor of
MODIS-based open water fraction which, in turn, was derived from a reference
Landsat image. Should the MNDWI6 be a suitable predictor of open water then we
shall use time series of this product for open water mapping thus
circumventing the acquisition of spotless Landsat images.

In this work we provide an approach for modeling water fraction as a function
of MNDWI6. Our application is tested on 2270 observations with good quality
assessment. These observations were taken from two MODIS images with spatial
resolution of 500 m, see Figure 1. These images were acquired on the 82 day of
2008, during the dry season of that year. The water fraction image was derived
from a Landsat one (acquired on the same date) which was used to classify
water presence; expert classification was applied to the Landsat image. These
images include Marismas Nacionales whose 133,854 ha of estuaries, small
patches of cedars, oaks, amapas and oil palm, coconut palm and white, red,
black and Chinese mangrove, among other vegetation types, makes it the most
important wetland in north-western Mexico; Marismas Nacionales is located at
$22^{\circ}07\mbox{'}N$ and $105^{\circ}31\mbox{'}W$ [6].

A first modeling approach consists of fitting a linear regression between the
water fraction and the MNDWI6. This model yields the line displayed in Figure
2-A; the intercept and slope estimates of this model are
$\hat{\beta}_{0}^{\mbox{\tiny{lrm}}}=0.6041$ and
$\hat{\beta}_{1}^{\mbox{\tiny{lrm}}}=1.1142$, respectively. A scatterplot of
the fitted values of this model against the water fraction variable, see
Figure 2-B, reveals that this approach is not appropriate for modeling
purposes; some of the values predicted by the linear model lie outside the
interval $(0,1)$.

Figure 2: A Linear regression of MNDWI6 and water fraction. B Scatterplot of
water fraction vs. fitted values.

In what follows we describe the novel flexible Beta regression (FBR) model
which, by construction, takes into account that the response variable (water
fraction) is bounded, hence, improving upon the linear model. Section II
presents the basics of the FBR model, as well as details of the Bayesian
estimation procedure employed to fit the FBR; a nested Metropolis-Hastings and
Gibbs sampling algorithm is discussed. The routines needed for this algorithm
were developed in the R language [7] and are available at the Github
repository https://github.com/inder-tg/FBR. Results are presented in Section
III and we provide an outlook of future work in the concluding Section IV.

##  II Methodology

Let $\boldsymbol{Y}^{\top}=(Y_{1},\ldots,Y_{n})$ denote $n$ independent
response variables, each one taking values on the interval $(0,1)$. Flexible
Beta regression was introduced by [8] as an alternative for modeling responses
like $\boldsymbol{Y}$, which are bounded, as a function of some covariates.
Briefly, we will discuss the main features of the flexible Beta regression
model.

###  II-A On flexible Beta random variables

The density function of a flexible Beta random variable (rv) is given by

$f_{FB}^{\ast}(y;\lambda_{1},\lambda_{2},\phi,p)=p\,f_{B}^{\ast}(y;\lambda_{1},\phi)+(1-p)\,f_{B}^{\ast}(y;\lambda_{2},\phi),$
(1)

where $y\in(0,1)$, $0<\lambda_{2}<\lambda_{1}<1$, $\phi>0$, $p\in(0,1)$ and
$f_{B}^{\ast}$ denotes the density of a mean-precision parametrized Beta rv:

$f_{B}^{\ast}(x;a,b)=\frac{\Gamma(b)}{\Gamma(ab)\,\Gamma((1-a)b)}\,x^{ab-1}(1-x)^{(1-a)b-1},$

where $x\in(0,1)$, $a\in(0,1)$, $b>0$.

That is, $f_{FB}^{\ast}$ is a mixture of two mean-precision Beta rvs with
different means ($\lambda_{1}$ and $\lambda_{2}$) but with the same precision
parameter ($\phi$). It is not difficult to see that the expected value of a rv
with density function given by Eq.(1) is equal to
$p\lambda_{1}+(1-p)\lambda_{2}$. The latter, gives rise to a reparametrization
of the flexible Beta: $\mu=p\lambda_{1}+(1-p)\lambda_{2}$, $\phi=\phi$,
$\tilde{\omega}=\lambda_{1}-\lambda_{2}$ and $p=p$. This reparametrization has
proven useful in enhancing parameter interpretability in the flexible Beta
regression model.

The flexible Beta density has a wide range of shapes, from unimodal, monotone
and U-shaped to heavy-tailed. This characteristic makes it attractive to
describe the main features of a broad group of variables with applications in
several scientific fields. In addition to this, under mild conditions, it is
known that two elements of the FB parametric family are equal if and only if
the corresponding parameters are the same (see Proposition 1 of [8]), this
type of strong identifiability property is typically found in non-mixture
models.

###  II-B The flexible Beta regression model

Let us suppose that the $i$-th response variable, $Y_{i}\in(0,1)$, is linked
to the vector of covariates $\boldsymbol{x}_{i}\in\mathbb{R}^{k}$. The
flexible Beta regression (FBR) model consists of assuming that each $Y_{i}$ is
independently distributed flexible Beta; in notation, $Y_{i}\sim
FB(\mu_{i},\phi,\omega,p)$ where for some _link_ function
$g:[0,1]\to\mathbb{R}$,

$g(\mu_{i})=\mathbb{E}[Y_{i}\mid\boldsymbol{x}_{i}]=\boldsymbol{x}_{i}^{\top}\,\boldsymbol{\beta},\quad\boldsymbol{\beta}\in\mathbb{R}^{k},$
(2)

and $\tilde{\omega}=\omega\,\min\\{\mu_{i}/p,(1-\mu_{i})/(1-p)\\}$. Observe
that this representation allows the parameters $\mu_{i}$, $\omega$ and $p$ to
vary freely on $(0,1)$ and $\phi>0$. The parameter $\boldsymbol{\beta}$ is
unknown.

Note also that the representation induced by the FBR model is equivalent to
requiring that

$\displaystyle\lambda_{1}$ $\displaystyle=\mu_{i}+(1-p)\,\tilde{\omega}$
$\displaystyle\lambda_{2}$ $\displaystyle=\mu_{i}-p\,\tilde{\omega},$ (3)

in Eq.(1). So that, in what follows we will use the notation
$f_{FB}^{\ast}(\cdot;\mu_{i},\phi,\omega,p)$ to refer to the density
associated to an observation of the FBR model. For our application, we will
consider the _logit_ as the link function in Eq.(2). We remark that the
flexible Beta regression is not a generalized linear model as the distribution
of mixture of Betas does not belong to the exponential family.

###  II-C Bayesian estimation

Let $\boldsymbol{y}^{\top}=(y_{1},\ldots,y_{n})$ denote $n$ independent
observations of the FBR model just introduced. The likelihood function of this
model is by definition

$\mathscr{L}(\boldsymbol{\eta}\mid\boldsymbol{y})=\prod_{i=1}^{n}\,f_{FB}^{\ast}(y_{i}\mid\mu_{i},\phi,\omega,p),$
(4)

where $\boldsymbol{\eta}=(\boldsymbol{\beta},\phi,\omega,p)$. According to
Proposition 2 of [8], this function is bounded from above almost surely, which
ensures the existence of a finite global maximum on the parameter space. This
result provides computational tractability to the FBR model, in particular for
Bayesian estimation algorithms.

The underlying assumption of the FBR is that there are two groups of
observations and the mean of one of these groups dominates the other, see
Eq.(II-B). Given the $i$-th observation, however, we do not know to which of
the two mixture components belongs. This makes the optimization of Eq.(4)
unfeasible. To ameliorate this, a Bayesian approach is chosen, in particular,
a data augmentation one is devised. More precisely, define the $n$-dimensional
random vector $\boldsymbol{v}$ with entries

$v_{i}=\begin{cases}1&\mbox{ if }y_{i}\mbox{ belongs to the first mixture
component,}\\\ 0&\mbox{ otherwise }\end{cases}.$

Although essentially these latent variables ($\boldsymbol{v}$) are missing
data, they can be included in the model and used for parameter estimation in a
2-step algorithm. In the first step, the parameter posterior distribution is
obtained conditional on $\boldsymbol{v}$, and in the second step, the
observations are classified ($\boldsymbol{v}$ is updated) conditional on
knowing the parameter.

A key element in this algorithm is the _complete-data_ likelihood function:

$\displaystyle\mathscr{L}_{CD}(\boldsymbol{\eta}\mid\boldsymbol{y},\boldsymbol{v})=\prod_{i=1}^{n}\,\left[pf_{B}^{\ast}(y_{i};\lambda_{1},\phi)\right]^{v_{i}}\left[(1-p)f_{B}^{\ast}(y_{i};\lambda_{2},\phi)\right]^{(1-v_{i})},$
(5)

where $\lambda_{1}$ and $\lambda_{2}$ are given by (II-B). Thus, given an
appropriate _prior_ distribution $\pi(\boldsymbol{\eta})$, the resulting
_posterior_ distribution is:

$\pi(\boldsymbol{\eta},\boldsymbol{v}\mid\boldsymbol{y})\propto\mathscr{L}_{CD}(\boldsymbol{\eta}\mid\boldsymbol{y},\boldsymbol{v})\,\pi(\boldsymbol{\eta}).$

Since there is no prior information available about the response variables
utilized in our application, and, as we argued above, the parameters of the
FBR model are strongly identifiable (also known as variation independent), it
is sensible to assume a factorized joint prior distribution:

$\pi(\boldsymbol{\eta})=\pi(\boldsymbol{\beta})\,\pi(\phi)\,\pi(\omega)\,\pi(p).$
(6)

Following [8] and [9], in our application, we used a $k$-dimensional Gaussian
prior
$\boldsymbol{\beta}\sim\mathscr{N}_{k}(\boldsymbol{0},\Sigma_{\boldsymbol{\beta}})$
for the regression parameters; we employed a gamma prior distribution
$\phi\sim\mathscr{G}(g,g)$ for the precision parameter, a rather usual choice,
see e.g. [10]; and we chose non-informative uniform priors for the parameters
$\omega\sim\mathscr{U}(0,1)$ and $p\sim\mathscr{U}(0,1)$. The $k\times k$
matrix $\Sigma_{\boldsymbol{\beta}}$ has _large_ values in its main diagonal
and zeroes otherwise. In the next section we provide details about the
algorithm that allows us to fit the FBR model based on the complete-data
approach just presented.

###  II-D Nested Metropolis-Hastings and Gibbs sampling

Some parameters of the FBR model have known conditional posterior
distributions, whereas the posterior distribution of others have to be
calculated numerically. For instance, it is not difficult to see that

$\pi(v_{i}\mid y_{i},\boldsymbol{\eta})=\mbox{Bernoulli}(\pi(v_{i}=1\mid
y_{i},\boldsymbol{\eta})),$ (7)

where

$\pi(v_{i}=1\mid
y_{i},\boldsymbol{\eta})=\frac{p\,f_{B}^{\ast}(y_{i}\mid\lambda_{1},\phi)}{p\,f_{B}^{\ast}(y_{i}\mid\lambda_{1},\phi)+(1-p)\,f_{B}^{\ast}(y_{i}\mid\lambda_{2},\phi)}.$

Note that by setting $n_{0}=\\#\\{t\mid v_{t}=0\\}$ and $n_{1}=\\#\\{s\mid
v_{s}=1\\}$,

$\pi(p\mid\boldsymbol{v})=\mbox{Beta}(n_{1}+1,n_{0}+1).$ (8)

Also, we get

$\displaystyle\pi(\tilde{\omega}\mid\boldsymbol{\beta},p)$
$\displaystyle=\mathscr{U}(0,r(\boldsymbol{\beta},p)),$ $\displaystyle
r(\boldsymbol{\beta},p)$
$\displaystyle=\min\left\\{\frac{\mu_{i}(\boldsymbol{\beta})}{p},\frac{1-\mu_{i}(\boldsymbol{\beta})}{1-p}\right\\},$
$\displaystyle\mu_{i}(\boldsymbol{\beta})$
$\displaystyle=\frac{\mbox{exp}(\boldsymbol{x}_{i}^{\top}\,\boldsymbol{\beta})}{1+\mbox{exp}(\boldsymbol{x}_{i}^{\top}\,\boldsymbol{\beta})}:=\mbox{logit}^{-1}(\boldsymbol{x}_{i}^{\top}\,\boldsymbol{\beta}),$
$\displaystyle\mbox{logit}(x)$ $\displaystyle=\log\left(x/(1-x)\right),\quad
x\in(0,1).$ (9)

Samples from the conditional posterior distribution of
$(\phi,\boldsymbol{\beta})$ given
$(\boldsymbol{y},\boldsymbol{v},p,\tilde{\omega})$ are obtained through
Metropolis-Hastings (M-H) algorithms, see Ch. 10 of [11]. In particular, for
the conditional posterior of $\phi$ we employed a random walk as the states
transfer (jumping) distribution. For the conditional posterior of
$\boldsymbol{\beta}$ the jumping distribution is a bivariate normal with mean
equal to the value of the chain in the previous step and covariance matrix
$\Sigma_{J}$, a diagonal matrix.

The following pseudocode summarizes the main steps needed to obtain samples
from the conditional posterior distribution
$\pi(\boldsymbol{\eta},\boldsymbol{v}\mid\boldsymbol{y})$.

0: $\boldsymbol{x}$, $\boldsymbol{y}$, nSamples, $p_{0}$, $\omega_{0}$,
$\phi_{0}$, $\boldsymbol{\beta}_{0}$, $\sigma_{\phi}$, $g$, $\Sigma_{J}$,
$\Sigma_{\boldsymbol{\beta}}$

0:
$\boldsymbol{\mu}^{(0)}=(\mbox{logit}^{-1}(\boldsymbol{x}_{1}^{\top}\,\boldsymbol{\beta}_{0}),\ldots,\mbox{logit}^{-1}(\boldsymbol{x}_{n}^{\top}\,\boldsymbol{\beta}_{0}))$,$\boldsymbol{\omega}^{(0)}=\omega_{0}\cdot\min\\{\boldsymbol{\mu}^{(0)}/p_{0},(1-\boldsymbol{\mu}^{(0)})/(1-p_{0})\\}$

for $j=1$ to nSamples do

$v_{\mbox{\tiny{posterior}}}\leftarrow\pi(\boldsymbol{v}\mid\boldsymbol{y},\boldsymbol{\omega}^{(j-1)},p_{j-1},\boldsymbol{\mu}^{(j-1)},\phi_{j-1})$

$p_{j}\leftarrow\pi(p\mid v_{\mbox{\tiny{posterior}}})$

$\boldsymbol{\mu}^{(j)}\leftarrow(\mbox{logit}^{-1}(\boldsymbol{x}_{1}^{\top}\,\boldsymbol{\beta}_{j-1}),\ldots,\mbox{logit}^{-1}(\boldsymbol{x}_{n}^{\top}\,\boldsymbol{\beta}_{j-1}))$

$\boldsymbol{\omega}^{(j)}\leftarrow\pi(\tilde{\omega}\mid\boldsymbol{\mu}^{(j)})$

$\phi_{i}\leftarrow\pi_{\mbox{\tiny{M-H}}}(\phi\mid\boldsymbol{y},v_{\mbox{\tiny{posterior}}},\boldsymbol{\mu}^{(j)},\boldsymbol{\omega}^{(j)},p_{j},\phi_{j-1},\sigma_{\phi},g)$

$\boldsymbol{\beta}_{j}\leftarrow\pi_{\mbox{\tiny{M-H}}}(\boldsymbol{\beta}\mid\boldsymbol{x},\boldsymbol{y},v_{\mbox{\tiny{posterior}}},p_{j},\omega_{i},\phi_{j},\boldsymbol{\mu}^{(j)},\boldsymbol{\beta}_{j-1},\Sigma_{J},\Sigma_{\boldsymbol{\beta}})$

end for

Details about burning periods, thinning, chain size and values of initial
parameters of this algorithm will be provided in the next section.

Having the sequence of vectors $\boldsymbol{\beta}_{j}$,
$j=1,\ldots,\mbox{nSamples}$, which are appropriate samples of the posterior
distribution of the regression parameter $\boldsymbol{\beta}$, we computed the
median of each column to estimate the coefficients of the FBR model (2):

$\displaystyle\hat{\beta}_{r}$
$\displaystyle=\mbox{median}(\boldsymbol{\beta}_{j}[,r]),\quad r=1,\ldots,k.$

Similarly, $\hat{p}=\mbox{median}(p_{1},\ldots,p_{\mbox{\tiny{nSamples}}})$,
$\hat{\boldsymbol{\mu}}=(\mbox{logit}^{-1}(\boldsymbol{x}_{1}^{\top}\,\hat{\boldsymbol{\beta}}),\ldots,\mbox{logit}^{-1}(\boldsymbol{x}_{n}^{\top}\,\hat{\boldsymbol{\beta}}))$
and
$\tilde{\omega}_{\mbox{\tiny{opt}}}=\min\\{\hat{\boldsymbol{\mu}}/\hat{p},(1-\hat{\boldsymbol{\mu}})/(1-\hat{p})\\}$.
Consequently, substituing these estimators in (II-B) we obtained the
regression functions $\hat{\lambda}_{1}$ and $\hat{\lambda}_{2}$.

##  III Results

For our application we considered $n=2270$ observations of the water fraction
product as the response variable ($y_{i}$) and an equal number of values of
the MNDWI6 index as their corresponding covariates ($z_{i}$); each $y_{i}$ and
$z_{i}$ were recorded at the same geographical position (pixel). We included
an intercept in our regression model, hence,
$\boldsymbol{x}_{i}^{\top}=(1,z_{i})$, that is, in the notation above, $k=2$.

For sampling from the conditional posterior distributions, we ran initially 20
chains each having size $\mbox{nSamples}=8000$. We used
$\boldsymbol{\beta}_{0}=(\hat{\beta}_{0}^{\mbox{\tiny{lrm}}},\hat{\beta}_{1}^{\mbox{\tiny{lrm}}})$.
Our results were stable across different values of $p_{0}$ and $\omega_{0}$;
we used small values for the diagonal of $\Sigma_{J}$ whereas large values for
the diagonal of $\Sigma_{\boldsymbol{\beta}}$. We tried different starting
values for the chain of the precision parameter $\phi$. We found stable chains
when $\phi_{0}$ was in a vicinity of 3, $\sigma_{\phi}=0.125$ and instead of
using a $\mbox{Gamma}(g,g)$ we employed a $\mbox{Gamma}(\kappa\,g,g)$ where
$\kappa=30,35$ and $g=0.1,0.2,0.3$.

As a burning strategy we discarded the first 4000 samples of each chain. The
remaining samples passed through the following thinning procedure to ensure no
serial dependence between the posterior samples. The empirical autocorrelation
function (ACF) applied to each chain revealed different dependence levels
(significant lags); we found 15 as the largest dependence level. Our thinning
strategy was to segment the sequence of samples in non overlapping intervals
such that the last member of a given interval and the first element of the
next one were separated at least 15 units. Then, at random we selected some
samples from these strategically separated intervals; we took 500 non
correlated samples on each chain. Thus, for each parameter of the FBR model we
ended up producing 10000 posterior samples. Figure 3 depicts posterior samples
of $p$ and $\phi$ along with their corresponding ACFs. Similarly, Figure 4
portraits the histogram of the posterior distribution of the flexible Beta
regression parameters $\beta_{0}$ and $\beta_{1}$ accompanied by their ACFs.

Figure 3: A-B Posterior samples of parameter $p$ with ACF. C-D Posterior
samples of precision parameter $\phi$ with ACF.

Figure 4: A-D Histograms and ACFs of the posterior distributions of the FBR
parameters $\beta_{0}$ and $\beta_{1}$.

We consider that our implementation of the nested Metropolis-Hastings and
Gibbs sampling algorithm and the subsequent application of our thinning
approach produce appropriate posterior samples of the FBR model. In order to
fit the FBR model we proceeded as explained at the end of Section II-D.

Figure 5: A Regression function $\hat{\lambda}_{2}$ (in orange) fitted via the
FBR model. B Scatterplot of water fraction vs. FBR fitted values.

Figure 5-A shows the $\hat{\lambda}_{2}$ regression function fitted by the FBR
model. As it is shown in Figure 5-B the predicted values of this model lie on
the interval (0,1), satisfying the minimal requirement of a predictor of water
fraction. The FBR model presented in this work is suitable for independent
observations, which is a strong assumption to hold true on satellite-derived
variables. In particular, the scatterplot of Figure 5-A shows levels of
heteroscedasticity in the measurements. As argued by [8] the FBR can be used
on heteroscedastic observations by allowing some parameters, for instance the
precision parameter $\phi$, to depend on covariates. This modeling approach
requires further study and it will be explored in the future.

##  IV Conclusions

Modeling water fraction measurements derived from satellite images can be
construed as the first step in measuring the dynamic of open water bodies. In
this work we presented the flexible Beta regression model as an alternative to
predict water fraction as a function of the water index MNDWI6.

We have provided an algorithm based on Bayesian estimation principles to fit
the FBR model. We consider that our implementation yields appropriate results
and is available to any interested user. We have shown that the predicted
values of the FBR model lie on the interval $(0,1)$. In this regard, this
model improves the prediction made by the linear model.

Arguably, measurements of water fraction and water index are both spatially
related. It is sensible to expect that this relation is still present when the
measurements are vectorized (which is the type of data that we employed in
this work). Our approach lacks of accounting for this apparent dependence
between observations. It has been noticed by [8] that by allowing the
precision parameter to depend on covariates it is plausible that the FBR model
may be able to account for heteroscedasticity in the observations. We envision
to explore this idea in our future research.

Subsequent steps for modeling open water fraction include the acquisition of
biophysical and geographical variables such as precipitation, temperature,
elevation and hydrological basin models and include them into our model. We
believe that with this extra information along with the improvements discussed
above our model will become a useful tool to assess the dynamics of the open
water bodies of Mexico.

##  V Acknowledgment

We are grateful to Berenice Vázquez with CONABIO for kindly providing the
images analyzed in this paper. Julia Trinidad Reyes was partly funded by
Academia Mexicana de Ciencias through a summer grant during “Verano de la
Investigación Científica 2020”.

## References

  * [1] W. Mitsch, B. Bernal, A. Nahlik, Ü. Mander, L. Zhang, C. Anderson, S. Jørgensen, and H. Brix, “Wetlands, carbon, and climate change,” Landscape Ecology, vol. 28, pp. 583–597, 2013.
  * [2] S. Kirschke, P. and Bousquet, P. Ciais, M. Saunois, J. Canadell, E. Dlugokencky, P. Bergamaschi, D. Bergmann, D. Blake, L. Bruhwiler, and others, “Three decades of global methane sources and sinks, Nature Geoscience, vol. 6, pp. 813–823, 2013.
  * [3] C. Prigent, E. Matthews, F. Aires, W. Rossow, “Remote sensing of global wetland dynamics with multiple satellite data sets,” Geophysical Research Letters, vol. 28, pp. 4631–4634, 2001.
  * [4] R. Colditz, C. Troche-Souza, B. Vazquez, A. Wickel, and R. Ressl, “Analysis of optimal thresholds for identification of open water using MODIS-derived spectral indices for two coastal wetland systems in Mexico,” Int J of Appl Earth Obs Geoinf, vol. 70, pp. 13–24, 2018.
  * [5] Hanqiu Xu, “Modification of Normalised Difference Water Index (NDWI) to enhance open water features in remotely sensed imagery,” International Journal of Remote Sensing, vol. 27, pp. 3025–3033, 2006.
  * [6] National Commission for the Knowledge and Use of Biodiversity, “Mapa de uso del suelo y vegetación de la zona costera asociada a los manglares de México en 2015,” 1st. Ed. Mangrove Monitoring System of Mexico (SMMM), Mexico City, Mexico, 2016.
  * [7] R Core Team, “R: A Language and Environment for Statistical Computing,” R Foundation for Statistical Computing, https://www.r-project.org/, 2016
  * [8] S. Migliorati, A. Di Brisco, and A. Ongaro, “A new regression model for bounded responses,” Bayesian Analysis, vol. 13, pp. 845–872, 2018.
  * [9] N. Bouguila, D. Ziou, and E. Monga, “Practical Bayesian estimation of a finite beta mixture through Gibbs sampling and its applications,” Statistics and Computing, vol. 16, pp. 215–225, 2006.
  * [10] A. Branscum, W. Johnson, and M. Thurmond, “Bayesian beta regression: applications to household expenditure data and genetic distance between foot-and-mouth disease viruses,” Australian & New Zealand Journal of Statistics, vol. 49, pp. 287–301, 2007.
  * [11] A. Gelman, J. Carlin, H. Stern, D. Dunson, A. Vehtari, and D. Rubin, Bayesian Data Analysis, CRC press, 2013.


# Influence of interface-induced valley-Zeeman and spin-orbit couplings
on transport in graphene-on-WSe2 heterostructures

M. Zubair<EMAIL_ADDRESS>Department of Physics, Concordia
University, 7141 Sherbrooke Ouest, Montreal, Quebec H4B 1R6, Canada  P.
Vasilopoulos<EMAIL_ADDRESS>Department of Physics, Concordia
University, 7141 Sherbrooke Ouest, Montreal, Quebec H4B 1R6, Canada  M. Tahir
<EMAIL_ADDRESS>Department of Physics, Colorado State University, Fort
Collins, CO 80523, USA

###### Abstract

We investigate the electronic dispersion and transport properties of
graphene/WSe2 heterostructures in the presence of a proximity induced spin-
orbit coupling (SOC) using a low-energy Hamiltonian, with different types of
symmetry breaking terms, obtained from a four-band, first and second nearest-
neighbour tight-binding (TB) one. The competition between different
perturbation terms leads to inverted SOC bands. Further, we study the effect
of symmetry breaking terms on ac and dc transport by evaluating the
corresponding conductivities within linear response theory. The scattering-
independent part of the valley-Hall conductivity, as a function of the Fermi
energy $E_{F}$, is mostly negative in the ranges $-\lambda_{R}\leqslant E_{F}$
and $E_{F}\geqslant\lambda_{R}$ when the strength $\lambda_{R}$ of the Rashba
SOC increases except for a very narrow region around $E_{F}=0$ in which it
peaks sharply upward. The scattering-dependent diffusive conductivity
increases linearly with electron density, is directly proportional to
$\lambda_{R}$ in the low- and high-density regimes, but weakens for
$\lambda_{R}=0$. We investigate the optical response in the presence of a SOC-
tunable band gap for variable $E_{F}$. An interesting feature of this SOC
tuning is that it can be used to switch on and off the Drude-type intraband
response. Furthermore, the ac conductivity exhibits interband responses due to
the Rashba SOC. We also show that the valley-Hall conductivity changes sign
when $E_{F}$ is comparable to $\lambda_{R}$ and vanishes at higher values of
$E_{F}$. It also exhibits a strong dependence on temperature and a
considerable structure as a function of the frequency.

##  I Introduction

Two-dimensional (2D) materials have become a hot topic in solid state physics,
especially since the discovery of graphene, both theoretically and
experimentally because of their prominent mechanical, optical, electrical and
magnetic properties [1]. Recently graphene has attracted a lot of attention in
the field of spintronics due to its large electronic mobility, low spin-orbit
coupling (SOC), negligible hyperfine interaction and gate tunability [2]. For
a clear example, it has been proven that graphene exhibits a very long spin
relaxation length even at room temperature [3, 4]. Due to the weak SOC though,
it is not a suitable candidate for the observation of important spin-dependent
phenomena including the spin-Hall effect [5] and anomalous Hall effect [6].

To render graphene useful in spintronics, several experimental groups used
different techniques to tailor the SOC strength in it through coupling with
foreign atoms or materials [7, 8, 9, 10, 11, 12, 13], such as graphene
hydrogenation [14, 15] or fluorination [16] as well as heavy adatom decoration
[17, 18]. However, these approaches not only reduce the transport quality, but
also make it difficult to reproduce [14, 15] and detect [16, 17, 18] the
induced SOC. To overcome these difficulties, graphene is recently grown on
different novel 2D materials, which are ideal candidates to induce SOC via
proximity effects [19, 20, 21, 22, 23, 24, 25]. Hexagonal boron nitride (BN)
has a weak SOC, and therefore, is not a suitable substrate for the proximity
effect [26]. The family of 2D transition metal dichalcogenides (TMDCs) are the
next best candidates, which have large direct band gaps and giant intrinsic
SOC [27, 28]. In this respect, graphene on TMDCs has been investigated for
transport [29, 30, 31] as well as intriguing technological applications,
including field-effect tunnelling transistors (FETTs), radio-frequency
oscillators, and efficient phototransistors [32, 33, 34, 35, 36, 37, 38].
Also, the proximity-induced SOC in graphene/TMDCs heterostructures has
recently been shown to depend [39, 40] on the twist angle between the lattice
of graphene and that of the TMDC.

In addition, it has been found in room-temperature experimental studies of the
spin-Hall effect that few-layer WS2 induces a large SOC in graphene, about
$17$ meV [41] as compared to the very weak one in pristine graphene [42].
Also, it has been unambiguously demonstrated experimentally that a room-
temperature spin-Hall effect in graphene is induced by MoS2 proximity [43].
Moreover, when graphene is placed on a multilayer WS2 substrate, an additional
valley-Zeeman SOC, due to the broken sublattice symmetry, along with the
Rashba SOC have been predicted theoretically and observed experimentally [25,
44, 45, 46]. This SOC induces a spin splitting of degenerate bands, with out-
of-plane spin polarization at the $K$ and $K^{\prime}$ points, and an opposite
spin splitting in different valleys. Analogous to the Zeeman splitting, the
SOC is termed valley-Zeeman because the effective Zeeman fields are valley-
dependent. It is the dominant SOC in TMDCs and is also predicted to be induced
in graphene on TMDCs [25, 44, 45, 46]. To our knowledge though, apart from
some spin-transport studies [47] and two experimental magneto-transport
studies [48], neither ac and dc scattering-dependent charge transport nor the
simultaneous effect of valley-Zeeman and Rashba SOCs have been theoretically
studied in graphene on TMDCs.

In this work we study in detail the effect of the valley-Zeeman and Rashba-
type SOCs on ac and dc transport in graphene/WSe2 heterostructures. There
results a mexican hat dispersion [49] contrary to other family memebers of
TMDCs , e.g., MoS2, WS2 etc. [50]. Such a dispersion leads to more features in
the optical conductivity when the Fermi level moves between the minimum and
maximum of the mexican hat. Also, we compare our results with those for
pristine graphene.

In Sec. II we specify the Hamiltonian and obtain the eigenvalues and
eigenfunctions in the presence of symmetry breaking terms. In Sec. III we
present general expressions for the conductivities and provide numerical
results. Conclusions and a summary follow in Sec. IV.

Figure 1: (a) Real-space graphene with $\vec{a}_{1}$ and $\vec{a}_{2}$ the
primitive lattice vectors. (b) Graphene’s first Brillouin zone and high
symmetry points $\Gamma$, $K$, and $M$ in reciprocal space. Its primitive
lattice vectors are $\vec{b}_{1}$ and $\vec{b}_{2}$. (c) Schematic
representation of graphene on a WSe2 substrate.

##  II Formulation

Graphene is a 2D, one-atom thick planar sheet of bonded carbon atoms densely
packed in a honeycomb structure as shown in Fig. 1 $(a)$. The lattice
structure can be viewed as a triangular lattice with two sites $A$ (red filled
spheres) and $B$ (blue filled spheres) per unit cell. The arrows indicate the
primitive lattice vectors $\vec{a}_{1}=a(1,0)$ and
$\vec{a}_{2}=a(1/2,\sqrt{3}/2)$, with $a$ the triangular lattice constant of
the structure, and span the graphene lattice. Further, $\vec{a}_{1}$ and
$\vec{a}_{2}$ generate the reciprocal lattice vectors of the Brillouin zone,
cf. Fig. 1 $(b)$, given by $\vec{b}_{1}=4\pi/\sqrt{3}a(\sqrt{3}/2,-1/2)$ and
$\vec{b}_{2}=4\pi/\sqrt{3}a(0,1)$. From the explicit expressions of
$\vec{b}_{1}$ and $\vec{b}_{2}$ we find the two inequivalent Dirac points
(valleys) given by $\vec{K}=(4\pi/3a)(1,0)$ and
$\vec{K}^{\prime}=(4\pi/3a)(1/2,\sqrt{3}/2)$.

Figure 2:  Energy dispersion in a graphene/WSe2 heterostructure using the TB
model (1) along the path $-M\rightarrow-K\rightarrow\Gamma\rightarrow
K\rightarrow M$ for (a) $\lambda_{c_{i}},\lambda_{R}=0$, (b)
$\lambda_{c_{i}}\neq 0,\lambda_{R}=0$, (c) $\lambda_{c_{i}}=0,\lambda_{R}\neq
0$, and (d) $\lambda_{c_{i}},\lambda_{R}\neq 0$ .

The monolayer graphene system is described by the four-band, second nearest-
neighbour tight-binding (TB) Hamiltonian [23, 49, 51]

$\displaystyle H$ $\displaystyle=$ $\displaystyle\sum_{\langle
i,j\rangle,\alpha}tc_{i\alpha}^{\dagger}c_{j\alpha}+\sum_{i\alpha}\Delta\eta_{c_{i}}c_{i\alpha}^{\dagger}c_{i\alpha}+\sum_{\langle\langle
i,j\rangle\rangle}\Delta_{ij}c_{i\alpha}^{\dagger}c_{j\alpha^{\prime}}$ (1)
$\displaystyle+\dfrac{2i}{3}\sum_{\langle
i,j\rangle}\sum_{\alpha\alpha^{\prime}}c_{i\alpha}^{\dagger}c_{j\alpha^{\prime}}[\lambda_{R}(\bm{s}\times\mathbf{\hat{d}}_{ij})_{z}]_{\alpha\alpha^{\prime}}.$

Here $\Delta_{ij}=i\lambda_{c_{i}}\nu_{ij}s_{z}/3\sqrt{3}$,
$c_{i\alpha}^{\dagger}$ creates an electron with spin polarization $\alpha$ at
site $i$ that belongs to sublattice $A$ or $B$, and $\langle i,j\rangle$
$(\langle\langle i,j\rangle\rangle)$ runs over the nearest (second nearest)
neighbouring sites. The second term is a staggered on-site potential, which
takes into account the effective energy difference experienced by atoms at the
lattice sites $A$ $(\eta_{c_{i}}=+1)$ and $B$ $(\eta_{c_{i}}=-1)$,
respectively. The third and fourth terms represent the proximity-induced
enhancement of the SOC due to a weak hybridization with the heavy atoms in
WSe2. The third term is the valley-Zeeman SOC where $\nu_{ij}=+1$, if the
second nearest hopping is anticlockwise with respect to the positive $z$ axis,
and $\nu_{ij}=-1$ if it is clockwise. The last term is the Rashba SOC
parametrized by $\lambda_{R}$. It arises because the inversion symmetry is
broken when the graphene sheet is placed on top of WSe2 as shown in Fig.1
$(c)$. Also, $\mathbf{\hat{d}}_{ij}=\mathbf{d}_{ij}/|\mathbf{d}_{ij}|$, where
$\bm{s}=(s_{x},s_{y},s_{z})$ is the Pauli spin matrix and $\mathbf{d}_{ij}$
the vector connecting the sites $i$ and $j$ in the same sublattice.

In Fig. 2 we plot the numerically evaluated energy dispersion of Eq. (1) to
better understand the characteristics of the induced intrinsic SOCs. Near the
$K$ point, for $\lambda_{c_{i}}=\lambda_{R}=0$, the band structure has linear
band crossings near $k=0$ as can be seen from Fig. 2 (a). For
$\lambda_{c_{i}}\neq 0$ and $\lambda_{R}=0$ the spectrum is gapless and the
spin degeneracy is broken away from $k=0$, see Fig. 2 (b). Further, if only
$\lambda_{R}$ is present, the spectrum is also gapless, cf. Fig. 2 (c).
However, a gap is created when both $\lambda_{c_{i}}$ and $\lambda_{R}$ are
finite, cf. Fig. 2 (d).

We analyze the physics of electrons near the Fermi energy using a low-energy
effective Hamiltonian derived from Eq. (1) and a Dirac theory around the $K$
and $K^{\prime}$ valleys [23, 44, 46]. It reads

$H^{s}_{\eta}=v_{F}(\eta\sigma_{x}p_{x}+\sigma_{y}p_{y})+\Delta\sigma_{z}+\lambda\sigma_{0}s\eta+\lambda_{R}(\eta
s_{y}\sigma_{x}-s_{x}\sigma_{y}).$ (2)

Here $\eta=\pm 1$ denotes the valleys $K$ and $K^{\prime}$, $\Delta$ is the
mass term that breaks the inversion symmetry, $\lambda=\lambda_{c_{i}}$ is the
valley-Zeeman SOC strength, $\lambda_{R}$ the Rashba type SOC strength,
($\sigma_{x}$, $\sigma_{y}$, $\sigma_{z}$) the Pauli matrix that corresponds
to the pseudospin (i.e., $A-B$ sublattice), $\sigma_{0}$ is the unit matrix in
the sublattice space, and $v_{F}$ ($8.2\times 10^{5}$ m/s) denotes the Fermi
velocity of Dirac fermions. For simplicity, we neglect the intrinsic SOC
$\lambda_{i}$ and consider only the $\lambda_{R}>\lambda_{i}$ case. Also, we
expect that small but finite values of $\lambda_{i}$ do not qualitatively
affect our results as long as $\lambda>>\lambda_{i}$. Further, we will also
neglect the $\Delta$ term in our numerical treatment because
$\lambda>>\Delta$.

Figure 3: Low-energy dispersion in a graphene/WSe2 heterostructure for
$\Delta=0$ and different combinations of $\lambda$ and $\lambda_{R}$.

Upon diagonalizing Eq. (2) we obtain the dispersion

$\displaystyle E_{\xi}(k)$
$\displaystyle=l\big{[}\Delta^{2}+\lambda^{2}+\hslash^{2}v_{F}^{2}k^{2}+2\lambda_{R}^{2}+2s\sqrt{\Upsilon}\,\big{]}^{1/2},$
(3)

where
$\Upsilon=\lambda_{R}^{2}\left(\lambda_{R}^{2}-2\lambda\Delta\right)+\hslash^{2}v_{F}^{2}k^{2}\left(\lambda_{R}^{2}+\lambda^{2}\right)+\lambda^{2}\Delta^{2}$
and $\xi=\\{l,s\\}$. Further, $l=+1(-1)$ denotes the conduction (valence) band
and $s=+1(-1)$ represents the spin-up (spin-down) branches. Notice that Eq.
(3) has a valley degeneracy despite the valley-Zeeman term. The normalized
eigenfunctions for both valleys are

$\psi_{\xi}^{+}(k)=\dfrac{N_{\xi}^{+}}{\sqrt{S_{0}}}\begin{pmatrix}1\\\
A_{\xi}^{\eta}e^{i\phi}\\\ -iB_{\xi}^{\eta}e^{i\phi}\\\
-iC_{\xi}^{\eta}e^{2i\phi}\end{pmatrix}e^{i{\bf k}\cdot{\bf r}},$ (4)
$\psi_{\xi}^{-}(k)=\dfrac{N_{\xi}^{-}}{\sqrt{S_{0}}}\begin{pmatrix}-A_{\xi}^{\eta}e^{i\phi}\\\
1\\\ iC_{\xi}^{\eta}e^{2i\phi}\\\
-iB_{\xi}^{\eta}e^{i\phi}\end{pmatrix}e^{i{\bf k}\cdot{\bf r}},$ (5)

respectively, with

$\displaystyle N_{\xi}^{\eta}$
$\displaystyle=l\big{[}1+(A_{\xi}^{\eta})^{2}+(B_{\xi}^{\eta})^{2}+(C_{\xi}^{\eta})^{2}\big{]}^{-1/2},$
(6)

$S_{0}=L_{x}L_{y}$ the area of the sample, and $\phi=\tan^{-1}(k_{y}/k_{x})$.
Further, $A_{\xi}^{\eta}=(E_{\xi}^{\eta}-\eta\Delta-\eta\lambda)/\hslash
v_{F}k$,
$B_{\xi}^{\eta}=2\lambda_{R}\big{[}(E_{\xi}^{\eta})^{2}-(\Delta+\lambda)^{2}\big{]}\big{/}\hslash
v_{F}k\big{[}(E_{\xi}^{\eta}+\eta\lambda)^{2}-\Delta^{2}-\hslash^{2}v_{F}^{2}k^{2}\big{]}$,
and
$C_{\xi}^{\eta}=2\lambda_{R}(E_{\xi}^{\eta}-\eta\Delta-\eta\lambda)\big{/}[(E_{\xi}^{\eta}+\eta\lambda)^{2}-\Delta^{2}-\hslash^{2}v_{F}^{2}k^{2}]$.

Figure 4: Density of states for $(\lambda,\lambda_{R})=(0,0)$,
$(\lambda,\lambda_{R})=(3,0)$ meV, $(\lambda,\lambda_{R})=(0,6)$ meV and
$(\lambda,\lambda_{R})=(3,6)$ meV. All cases are for $\Delta=0$.

We plot Eq. (3) in Fig. 3 for different combinations of the $\lambda$ and
$\lambda_{R}$ terms whose realistic values fall in the ranges $5-6$ meV and
$10-15$ meV, respectively, as determined experimentally in Ref. 52. Here, the
larger values of SOCs are used just to see well-resolved bands splitting. For
$\lambda=\lambda_{R}=0$, the band structure has linear bands crossing near
$k=0$ for both valleys as can be seen from panel (a). For $\lambda\neq 0$ and
$\lambda_{R}=0$, the energy dispersion is spin non-degenerate and valley
degenerate with a gapless behaviour as shown in panel (b). Further, the energy
dispersion shows the gapless behaviour for $\lambda=0$ and $\lambda_{R}\neq 0$
whereas it is spin-split as seen from panel (c). However, for $\lambda$ and
$\lambda_{R}$ finite, the Rashba coupling not only creates a gap between the
conduction and valence band, by mixing the spin-up and spin-down states, but
also produces an avoided crossing, see Fig. 1 (d). The analytical form of the
momentum $k_{1}$, at which an avoided crossing occurs, and of the gap
$E_{g}=\Delta_{1}$ are

$k_{1}={1\over\hslash
v_{F}}\Big{[}\dfrac{(\lambda^{2}+\lambda\Delta)(\lambda^{2}+2\lambda_{R}^{2}-\lambda\Delta)}{\lambda^{2}+\lambda_{R}^{2}}\Big{]}^{1/2},$
(7)
$\Delta_{1}=2\lambda_{R}\Big{[}\dfrac{\lambda^{2}+\Delta(2\lambda+\Delta)}{\lambda^{2}+\lambda_{R}^{2}}\Big{]}^{1/2}.$
(8)

The density of states (DOS) per unit area corresponding to Eq. (3) is given by
$D(E)=\sum_{\zeta}\delta(E-E_{\zeta})$ with
$\left|\zeta\right\rangle=\left|\xi,\eta,k\right\rangle$. For $\lambda_{R}=0$
it takes the simple form

$D(E)=\dfrac{1}{2\pi\hslash^{2}v_{F}^{2}}\sum_{\xi}\left|\dfrac{E}{l}-s\lambda\right|\Theta(\dfrac{E}{l}-s\lambda-\Delta),$
(9)

and for $\Delta=\lambda=0$ the form

$D(E)=\dfrac{1}{2\pi\hslash^{2}v_{F}^{2}}\sum_{\xi}\left|\dfrac{E}{l}-s\lambda_{R}\right|\Theta(\dfrac{E}{l}-(s+1)\lambda_{R}).$
(10)

The DOS is shown in Fig. 4 for several values of $\lambda$ and $\lambda_{R}$.
The black curve is for monolayer graphene, with $\lambda=\lambda_{R}=0$, and
is included for comparison. The $E_{+-}$ and $E_{++}$ dispersions give rise to
a square root singularity at
$E=\lambda\lambda_{R}/\sqrt{\lambda^{2}+\lambda_{R}^{2}}$ and a step at
$E=\sqrt{\lambda^{2}+4\lambda_{R}^{2}}$, respectively, as shown by the black
dot-dashed curve of Fig. 4. The origin of the singularity is the mexican-hat
energy dispersion, cf. Fig. 3. In addition, the step emerges from the bottom
of the $E_{++}$ band and is a van Hove singularity associated with the
dispersion flattening at this point. The square root singularity is calculated
near the mexican-hat minimum
$E=\lambda\lambda_{R}/\sqrt{\lambda^{2}+\lambda_{R}^{2}}$ at which $D(E)$
reads

$D(E)=\dfrac{k_{1}}{4\pi\hslash}\sqrt{\dfrac{2m^{\ast}}{E-\Delta_{1}}},$ (11)

with
$m^{\ast}=\lambda_{R}(\lambda^{2}+\lambda_{R}^{2})^{3/2}/2v_{F}^{2}\lambda(\lambda^{2}+2\lambda_{R}^{2})$
the effective mass and
$E_{+,-}=\Delta_{1}+(\hslash^{2}/2m^{\ast})(k-k_{1})^{2}$ the energy. This
singularity is similar to that of the one-dimensional density of states. In
the limit $\lambda_{R}=0$ and $\lambda\neq 0$, the DOS has a finite value
$\lambda/2\pi\hslash^{2}v_{F}^{2}$ at $E=0$ (see blue dashed curve). For
$E\geqslant\lambda$, it increases linearly with $E$. Also, for $\lambda=0$ and
$\lambda_{R}\neq 0$, it is finite at $E=0$ but has a step at $E=2\lambda_{R}$,
see the red dotted curve.

Figure 5: Valley-Hall conductivity versus Fermi energy at $T=0.5$ K. For further clarity, the range -$0.35$ meV $\leq E_{F}\leq 0.35$ meV is shown in the inset without the $\lambda_{R}=0$ curve. Table 1: Band gap energies involved in optical transitions, cf. Fig. 7, for $\lambda=8$ meV, $\lambda_{R}=6$ meV, and two values of $E_{F}$. Transition energies | Formula |  $E_{F}=6.6$ meV $E_{F}=9.6$ meV |
---|---|---|---
$\Delta_{1}$ | $2\lambda\lambda_{R}/\sqrt{\lambda^{2}+\lambda_{R}^{2}}$ | 9.6 | 9.6
$\Delta_{2}$ | $2\sqrt{(4\lambda^{4}+4\lambda_{R}^{4}+9\lambda^{2}\lambda_{R}^{2})/(\lambda^{2}+\lambda_{R}^{2})}$ | 41.2 | 41.2
$\Delta_{01}$ | $2\lambda$ | 16 | 16
$\Delta_{02}$ | $2\sqrt{\lambda^{2}+4\lambda_{R}^{2}}$ | 28.8 | 28.8
$\Delta_{a}$ | $2\sqrt{2\lambda^{2}+2\lambda_{R}^{2}+E_{F}^{2}-2M-2L}$ | 32.2 |
$\Delta_{b}$ | $2\sqrt{2\lambda^{2}+2\lambda_{R}^{2}+E_{F}^{2}+2M+2L}$ | 50 | 57.4
|  $M=\sqrt{(\lambda^{2}+\lambda_{R}^{2})E_{F}^{2}-\lambda^{2}\lambda_{R}^{2}}$, $L=\sqrt{\lambda_{R}^{4}+(\lambda^{2}+\lambda_{R}^{2})(E_{F}^{2}+\lambda^{2}\pm 2M)}$ |  |

##  III  Conductivities

We consider a many-body system described by the Hamiltonian
$H=H_{0}+H_{I}-\mathbf{R\cdot F}(t)$, where $H_{0}$ is the unperturbed part,
$H_{I}$ is a binary-type interaction (e.g., between electrons and impurities
or phonons), and $\mathbf{-R\cdot F}(t)$ is the interaction of the system with
the external field F(t) [53]. For conductivity problems we have
$\mathbf{F}(t)=e\mathbf{E}(t)$, where $\mathbf{E}(t)$ is the electric field,
$e$ the electron charge, $\mathbf{R}=\sum_{i}{\bf r}_{i}$ , and
$\mathbf{r}_{i}$ the position operator of electron $i$. In the representation
in which $H_{0}$ is diagonal the many-body density operator
$\rho=\rho^{d}+\rho^{nd}$ has a diagonal part $\rho^{d}$ and a nondiagonal
part $\rho^{nd}$. For weak electric fields and weak scattering potentials, for
which the first Born approximation applies, the conductivity tensor has a
diagonal part $\sigma_{\mu\nu}^{d}$ and a nondiagonal part
$\sigma_{\mu\nu}^{nd}$; the total conductivity is
$\sigma_{\mu\nu}^{T}=\sigma_{\mu\nu}^{d}+\sigma_{\mu\nu}^{nd},\mu,\nu=x,y$.

In general we have two kinds of currents, diffusive and hopping, with
$\sigma_{\mu\nu}^{d}=\sigma_{\mu\nu}^{dif}+\sigma_{\mu\nu}^{col}$, but usually
only one of them is present. If no magnetic field is present, the hop-
ping term $\sigma_{\mu\nu}^{col}$ vanishes identically and only the term
$\sigma_{\mu\nu}^{dif}$ survives. For elastic scattering it is given by [53]

$\sigma_{\mu\nu}^{d}(\omega)=\dfrac{\beta
e^{2}}{S_{0}}\sum_{\zeta}f_{\zeta}(1-f_{\zeta})\dfrac{v_{\nu\zeta}\,v_{\mu\zeta}\,\tau_{\zeta}}{1+i\omega\tau_{\zeta}},$
(12)

with $\tau_{\zeta}$ the momentum relaxation time, $\omega$ the frequency, and
$v_{\mu\zeta}$ the diagonal matrix elements of the velocity operator. Further,
$f_{\zeta}=\big{[}1+\exp[\beta(E_{\zeta}-E_{F})]\big{]}^{-1}$ is the Fermi-
Dirac distribution function, $\beta=1/k_{B}T$ and $T$ the temperature.

Figure 6: Same as in Fig. 5 but for different values of $T$.

Regarding the contribution $\sigma_{\mu\nu}^{nd}$ one can use the identity
$f_{\zeta}(1-f_{\zeta^{\prime}})\big{[}1-\exp[\beta(E_{\zeta}-E_{\zeta^{\prime}})]\big{]}=f_{\zeta}-f_{\zeta^{\prime}}$
and cast the original form in the more familiar one [53]

$\sigma_{\mu\nu}^{nd}(\omega)=\dfrac{i\hslash
e^{2}}{S_{0}}\sum_{\zeta\neq\zeta^{\prime}}\dfrac{(f_{\zeta}-f_{\zeta^{\prime}})\,v_{\nu\zeta\zeta^{\prime}}\,v_{\mu\zeta\zeta^{\prime}}}{(E_{\zeta}-E_{\zeta^{\prime}})(E_{\zeta}-E_{\zeta^{\prime}}+\hslash\omega-i\Gamma)},$
(13)

where the sum runs over all quantum numbers $\zeta$ and $\zeta^{\prime}$ with
$\zeta\neq\zeta^{\prime}$. The infinitesimal quantity $\epsilon$ in the
original form has been replaced by $\Gamma_{\zeta}$ to account for the
broadening of the energy levels. In Eq. (13) $v_{\nu\zeta\zeta^{\prime}}$ and
$v_{\mu\zeta\zeta^{\prime}}$ are the off-diagonal matrix elements of the
velocity operator. The relevant velocity operators are given by
$v_{x}=\partial H/\hslash\partial k_{x}$ and $v_{y}=\partial H/\hslash\partial
k_{y}$. With $\zeta=\\{l,s.k,\eta\\}=\\{\xi,k,\eta\\}$ for brevity, they read

$\left\langle\zeta\right|v_{x}\left|\zeta^{\prime}\right\rangle=v_{F}N_{\xi}^{\eta}N_{\xi^{\prime}}^{\eta}(D_{\xi,\xi^{\prime}}^{\eta}e^{i\phi}+F_{\xi,\xi^{\prime}}^{\eta}e^{-i\phi})\delta_{k,k^{\prime}},$
(14)
$\left\langle\zeta^{\prime}\right|v_{y}\left|\zeta\right\rangle=iv_{F}N_{\xi}^{\eta}N_{\xi^{\prime}}^{\eta}(D_{\xi,\xi^{\prime}}^{\eta}e^{-i\phi}-F_{\xi,\xi^{\prime}}^{\eta}e^{i\phi})\delta_{k,k^{\prime}},$
(15)

where
$D_{\xi,\xi^{\prime}}^{\eta}=A_{\xi^{\prime}}^{\eta}+B_{\xi}^{\eta}C_{\xi^{\prime}}^{\eta}$
and
$F_{\xi,\xi^{\prime}}^{\eta}=A_{\xi}^{\eta}+B_{\xi^{\prime}}^{\eta}C_{\xi}^{\eta}$.

We now calculate the conductivity $\sigma_{yx}^{nd}(i\omega)$ given by Eq.
(13). Further, the velocity matrix elements (14) and (15) are diagonal in $k$,
therefore $k$ will be suppressed in order to simplify the notation. The
summation in Eq. (13) runs over all quantum numbers $\xi$,$\xi^{\prime}$,
$\eta$, $\eta^{\prime}$, and $k$. The parameter
$\Gamma_{\eta\eta^{\prime}}^{\xi\xi^{\prime}}$, that takes into account the
level broadening, is assumed to be independent of the band and valley indices,
i.e., $\Gamma_{\eta\eta^{\prime}}^{\xi\xi^{\prime}}=\Gamma$. Using Eqs. (14)
and (15) we can express Eq. (13) as

$\displaystyle\sigma_{yx}^{nd}(i\omega)$ $\displaystyle=$
$\displaystyle\dfrac{e^{2}\hslash^{2}v_{F}^{2}}{h}\sum_{\xi\xi^{\prime}}\int
dkk\,\dfrac{(N_{\xi}^{\eta}N_{\xi^{\prime}}^{\eta})^{2}(f_{\xi
k}^{\eta}-f_{\xi^{\prime}k}^{\eta})}{\Delta_{\xi\xi^{\prime}}^{\eta}\big{[}(\Delta_{\xi\xi^{\prime}}^{\eta}+\hslash\omega)^{2}+\Gamma^{2}\big{]}}$
(16)
$\displaystyle\times\big{[}\Delta_{\xi\xi^{\prime}}^{\eta}+\hslash\omega-i\Gamma)\big{]}\big{[}D_{\xi,\xi^{\prime}}^{\eta})^{2}-(F_{\xi,\xi^{\prime}}^{\eta})^{2}\big{]}$

where $\Delta_{\xi\xi^{\prime}}^{\eta}=E_{\xi
k}^{\eta}-E_{\xi^{\prime}k}^{\eta}$. Further, in the limit $\Gamma=\omega=0$,
Eq. (16) reduces to

$\displaystyle\sigma_{yx}^{nd}$ $\displaystyle=$
$\displaystyle\dfrac{e^{2}\hslash^{2}v_{F}^{2}}{h}\sum_{\xi\xi^{\prime}}\int
dkk\,\dfrac{(N_{\xi}^{\eta}N_{\xi^{\prime}}^{\eta})^{2}(f_{\xi
k}^{\eta}-f_{\xi^{\prime}k}^{\eta})}{(\Delta_{\xi\xi^{\prime}}^{\eta})^{2}}$
(17)
$\displaystyle\times\big{[}(D_{\xi,\xi^{\prime}}^{\eta})^{2}-(F_{\xi,\xi^{\prime}}^{\eta})^{2}\big{]}$

In the valley-Hall effect electrons from regions near the inequivalent $K$ and
$K^{\prime}$ valleys flow to opposite transverse edges of the system, in the
presence of SOCs when a longitudinal electric field is applied [54, 55]. The
valley-Hall conductivity corresponding to Eq. (16) is defined by

$\sigma_{yx}^{v}=\sum_{ss^{\prime}}\sigma_{yx}^{nd}(\eta=+,s,s^{\prime})-\sigma_{yx}^{nd}(\eta=-,s,s^{\prime}).$
(18)

The spin-Hall conductivity $\sigma_{yx}^{s}$, corresponding to Eq. (16), is
finite only when both the Kane-Mele and valley- Zeeman SOCs are present.
Hence, even in the presence of Rashba SOC, $\sigma_{yx}^{s}$ vanishes [43].
Since a spin current is defined by
$\bf{J}_{s}=(\hslash/2e)(\bf{J}_{\uparrow}-\bf{J}_{\downarrow})$, we have to
multiply $\sigma_{yx}^{v}$ by $1/2e$ [5, 56]. Further, we find that charge
Hall conductivity always vanishes

$\sigma_{yx}^{c}=\sum_{\eta ss^{\prime}}\sigma_{yx}^{nd}(\eta,s,s^{\prime})=0$
(19)

Figure 7: Band structure near the Fermi energy $E_{F}$ in the presence of SOC
terms for $\lambda=8$ meV and $\lambda_{R}=6$ meV. The black dashed and dotted
lines show $E_{F}=6.6$ meV and $E_{F}=9.6$ meV. The various gap energies,
indicated by $\Delta_{1},\Delta_{2}$, etc. are displayed in table 1. Notice
that for $E_{F}=9.6$ meV the energy $\Delta_{a}$ does not contribute to any
transitions.

The component $\sigma_{xx}^{nd}(i\omega)$ is also obtained from Eq. (13):

$\displaystyle\sigma_{xx}^{nd}(i\omega)$ $\displaystyle=$
$\displaystyle\dfrac{ie^{2}\hslash^{2}v_{F}^{2}}{h}\sum_{\eta\xi\xi^{\prime}}\int
dkk\,\dfrac{(N_{\xi}^{\eta}N_{\xi^{\prime}}^{\eta})^{2}(f_{\xi
k}^{\eta}-f_{\xi^{\prime}k}^{\eta})}{\Delta_{\xi\xi^{\prime}}^{\eta}\big{[}(\Delta_{\xi\xi^{\prime}}^{\eta}+\hslash\omega)^{2}+\Gamma^{2}\big{]}}$
(20)
$\displaystyle\times\big{[}\Delta_{\xi\xi^{\prime}}^{\eta}+\hslash\omega-i\Gamma)\big{]}\big{[}D_{\xi,\xi^{\prime}}^{\eta})^{2}+(F_{\xi,\xi^{\prime}}^{\eta})^{2}\big{]}.$

For $\lambda=0$ and $\lambda_{R}\neq 0$, Eq. (16) vanishes because the factor
$(D_{\xi,\xi^{\prime}}^{\eta})^{2}-(F_{\xi,\xi^{\prime}}^{\eta})^{2}$ becomes
zero, whereas Eq. (20) survives. Moreover, in the limit
$\lambda=\lambda_{R}=0$, Eq. (20) reduces to the optical conductivity of
pristine graphene, which is independent of $\hslash\omega$ and given by
$e^{2}/2h$ [57].

We now consider the diagonal component $\sigma_{xx}^{d}$ given by Eq. (12).
Using Eq. (14), with $\xi=\xi^{\prime}$, we obtain

$\displaystyle\sigma_{xx}^{d}(i\omega)$ $\displaystyle=$
$\displaystyle\dfrac{e^{2}v_{F}^{2}\beta}{\pi}\sum_{\eta\xi}\int
dkk\,(N_{\xi}^{\eta})^{4}f_{\xi k}^{\eta}(1-f_{\xi k}^{\eta})$ (21)
$\displaystyle\times\dfrac{(A_{\xi}^{\eta}+B_{\xi}^{\eta}C_{\xi}^{\eta})^{2}\,\tau_{\xi
k}^{\eta}}{1+i\omega\tau_{\xi k}^{\eta}}$

At very low temperatures we can make the approximation $\beta f_{\xi
k}^{\eta}(1-f_{\xi k}^{\eta})\approx\delta(E_{\xi}-E_{F})$ and $\tau_{\xi
k}^{\eta}=\tau_{\xi k_{F}}^{\eta}$ because all states untill the Fermi level
are occupied.

Figure 8: Band structure near the Fermi energy $E_{F}$ in the presence of SOC
terms for $\lambda=8$ meV and $\lambda_{R}=6$ meV. Black and red arrows
represent possible interband transitions. Red arrows indicate the Drude type
intraband transitions. (a) $\lambda\neq 0$, $\lambda_{R},E_{F}=0$. (b)
$\lambda\neq 0$, $\lambda_{R}=0$, $E_{F}=1$ meV. (c) $\lambda_{R}\neq 0$,
$\lambda,E_{F}=0$. (d) $\lambda_{R}\neq 0$, $\lambda_{R}=0$, $E_{F}=2.8$ meV.

In Fig. 5 we plot Eq. (16) in the dc limit ($\omega=0$) as a function of
$E_{F}$ for $\Gamma=0.2$ meV, $\lambda=3$ meV and for different values of
$\lambda_{R}$. When $E_{F}$ is in the gap, i.e., in the range
$-\lambda\lambda_{R}/\sqrt{\lambda^{2}+\lambda_{R}^{2}}$ $\leqslant
E_{F}\leqslant\lambda\lambda_{R}/\sqrt{\lambda^{2}+\lambda_{R}^{2}}$, the
valley-Hall conductivity is quantized in units of $2e/2h$ similar to the case
of gapped graphene and topological insulators [54, 58]. The reason is that the
factor
$\sum_{\eta\xi\xi^{\prime}}(N_{\xi}^{\eta}N_{\xi^{\prime}}^{\eta})^{2}[(D_{\xi,\xi^{\prime}}^{\eta})^{2}-(F_{\xi,\xi^{\prime}}^{\eta})^{2}]/(\Delta_{\xi\xi^{\prime}}^{\eta})^{2}$,
called Berry curvature $\Omega(k)$, of Eq. (16) in the limit $\omega=0$ has a
peak, which is well covered by occupied states for
$E_{F}>\lambda\lambda_{R}/\sqrt{\lambda^{2}+\lambda_{R}^{2}}$. As a
consequence, the valley-Hall conductivity approaches the quantized value. For
$\lambda\lambda_{R}/\sqrt{\lambda^{2}+\lambda_{R}^{2}}$ $\leqslant
E_{F}\leqslant\lambda_{R}$, $\sigma_{yx}^{v}$ decreases with $E_{F}$. Further,
as can be seen, when $E_{F}$ becomes comparable to $\lambda_{R}$, a sign
change occurs in the conductivity which later vanishes at higher values of
$E_{F}$, $E_{F}>>\sqrt{\lambda^{2}+4\lambda_{R}^{2}}$. The change in sign is
due to the Rashba coupling between the spin-up and spin-down bands.
Furthermore, this off-diagonal term in spin space permits transitions between
two conduction spin subbands (see Eq. (3)), that could be interpreted as spin-
flip transitions near the band touching. In addition, the coupling strength
between opposite spin bands becomes weaker as $\lambda_{R}$ increases. As a
result, the negative part of the conductivity due to the spin-up band
diminishes and $\sigma_{yx}^{v}$ shows the usual behaviour of gapped graphene
and topological insulators [54, 58]. Further, as can be seen in the inset, the
band gap increases with $\lambda_{R}$. Also, the value of the conductivity at
$E_{F}=0$ is due to the finite one of $\Gamma$ (= 0.2 meV); if we take
$\Gamma=0$, the conductivity diverges at $E_{F}=0$ but its overall qualitative
behavior remains as shown.

Figure 9: Real part of longitudinal conductivity versus photon energy at
$T=0.5$ K. The upper panel is for $\lambda_{R}=0$ and the lower one is for
$\lambda_{R}\neq 0$.

We now take into account the effect of temperature $T$ on the valley-Hall
conductivity contained in the Fermi function, which is independent of
electron-phonon interaction in the first Born approximation [53]. The valley-
Hall conductivity is evaluated numerically with the help of Eq. (16) and
plotted in Fig. 6 for four values of $T$. We find a strong $T$ dependence,
particularly when the Fermi level is in the gap. The quantization of the
valley-Hall conductivity is destroyed at high values of $T$. This occurs when
the thermal broadening $k_{B}T$ becomes comparable to the energy gap. Notice
that the effect of temperature on $\sigma_{yx}^{v}$ is similar to that on the
spin-Hall conductivity in a graphene/MoS2 heterostructure by considering
valley-Zeeman and Kane-Mele SOCs in the absence of the Rashba SOC.

Various transition energies, which play an important role in the optical
conductivity, are shown in Fig. 7 for $\lambda,\lambda_{R}\neq 0$. Their
analytical expressions are displayed in table 1. Notice that for $E_{F}=6.6$
meV, the energies $\Delta_{a}$ and $\Delta_{b}$, indicated with black arrows,
become also important in optical transitions, since $E_{F}$ crosses the curve
$E_{+-}$ at two values of the momentum. However, for $E_{F}=9.6$ meV, only
$\Delta_{b}$ contributes to optical transitions because $E_{F}$ cuts $E_{+-}$
curve only at one value of the momentum. In Fig. 8, we show possible allowed
interband and intraband transitions by contrasting the case $\lambda\neq
0,\lambda_{R}=0$ in the upper panels and the case $\lambda=0,\lambda_{R}\neq
0$ in the lower panels. The blue arrows represent the interband transitions
$E_{+-}\rightarrow E_{++}$ for $0<E_{F}<\lambda$ and $0<E_{F}<\lambda_{R}$ as
can be seen in Fig. 8 (b) and (d). The black arrows represent the allowed
interband transitions $E_{-+}\rightarrow E_{+-}(E_{++})$ and
$E_{--}\rightarrow E_{+-}(E_{++})$ for $E_{F}=0$ and $E_{F}\neq 0$,
repectively, while the red arrows indicate intraband transitions that occur
near $E_{F}$.

Figure 10: Same as in Fig. 9 but for $\lambda,\lambda_{R}\neq 0$.

Now we present results for the real part of Eqs. (20) and (21)
$(\textrm{Re}\sigma_{xx}=\textrm{Re}\sigma_{xx}^{d}+\textrm{Re}\sigma_{xx}^{nd})$,
evaluated numerically, versus $\hslash\omega$ using a Lorentzian form of Dirac
delta function and taking $\Gamma=0.2$ meV for $T\neq 0$. We start from the
upper panel of Fig. 9 by considering the case $\lambda\neq 0$ and
$\lambda_{R}=0$. The transitions are vertical for photon’s momentum $q\sim 0$
and connect the filled valence band to empty conduction band, see Fig. 8 (a).
For the case of $E_{F}=0$, intraband response appears due to the transition
$E_{+-}\rightarrow E_{+-}$ and has a $\delta$ function form, centred around
$\hslash\omega=0$, which broadens the peak when any kind of scattering is
taken into account. Further, intraband responses occur when the Fermi level is
located away from the Dirac point. For $\hslash\omega=2\lambda$ we obtain
another Dirac delta peak due to the transition from $E_{--}\rightarrow
E_{+-}$, which is also broadened through $i\pi\delta(x)=\lim_{\Gamma\to\
0}(1/x-i\Gamma)$, cf. Eq. (20). For $0<E_{F}<\lambda$, the new absorption
peaks appear at $\hslash\omega=2E_{F}$ and $\hslash\omega=2(E_{F}+\lambda)$
due to the possible transitions $E_{--}\rightarrow E_{+-}$ and
$E_{-+}\rightarrow E_{+-}$. For $E_{F}>\lambda$, the absorption peaks
disappear below $\hslash\omega<2\lambda$ because the transition
$E_{--}\rightarrow E_{+-}$ is no longer possible due to the filling of states
below the Fermi level that are Pauli blocked. Further, the Drude peak persists
at low $\hslash\omega$, but now two other pieces of interband transitions
emerge with onsets at $\Delta_{a}+E_{F}$ and $\Delta_{b}+E_{F}$.

Figure 11: Valley-Hall conductivity versus photon energy for $\lambda_{R}=0$
at $T=0.5$ K.

In the lower panel of Fig. 9 we show the results for real part of the
longitudinal conductivity for $\lambda=0$, $\lambda_{R}\neq 0$ for different
values of $E_{F}$. For $E_{F}=0$, we can see that there is a peak at
$2\lambda_{R}$ which is the separation between $E_{--}$ and $E_{+-}$ bands. In
addition, there is a kink at $4\lambda_{R}$ due to the transition
$E_{-+}\rightarrow E_{++}$. As we increase the Fermi level, say,
$0<E_{F}<\lambda_{R}$ and $E_{F}>\lambda_{R}$, the peak becomes sharper and we
see a onset of a Drude contribution at low $\hslash\omega$ due to intraband
transitions $E_{+-}\rightarrow E_{+-}$ and $E_{++}\rightarrow E_{++}$ in
contrast to $E_{F}=0$ case (black dot-dashed curve). Further, for finite
values of $E_{F}$, we see the steps at $2E_{F}$ similar to monolayer graphene
($\lambda=\lambda_{R}=0$) as well as features at $\Delta_{a}+E_{F}$,
$\Delta_{b}-E_{F}$, and $\Delta_{b}+E_{F}$ above which we attain the flat
absorption like pristine graphene [57]. Note that our results are similar to
bilayer graphene [59, 60]. But here, the Rashba SOC, which allows the
interband transitions between opposite spin bands, gives rise to the
absorption peaks, while these peaks in bilayer graphene are due to interlayer
hopping between two graphene sheets.

Figure 12: Same as in Fig. 11 but with the upper panel for
$\lambda_{R}>\lambda$ and the lower one for $\lambda_{R}<\lambda$.

The real part of the longitudinal conductivity as a function of the photon
energy, for $\lambda,\lambda_{R}\neq 0$, is show in Fig. 10 for several values
of $E_{F}$: (i) just below the maximum of the mexican hat i.e.
$\lambda\lambda_{R}/(\lambda^{2}+\lambda_{R}^{2})^{1/2}<E_{F}<\lambda$ (ii)
just above the mexican hat, i.e., for
$\lambda<E_{F}<(\lambda^{2}+2\lambda_{R}^{2})^{1/2}$. For $E_{F}=0$ we find a
large absorption peak at approximately $2\lambda_{R}$, which corresponds to
transitions between the two square-root singularities of the DOS, see Fig. 4,
or transitions between the two minima of the mexican hat structures of the
$E_{--}$ and $E_{-+}$ bands. As $E_{F}$ moves into the mexican hat, this
feature disappears because states below $E_{F}$ are occupied and, therefore,
Pauli blocked. Further, the major peaks are due to the transitions
$E_{+-}\rightarrow E_{++}$, $E_{--}\rightarrow E_{-+}$, $E_{--}\rightarrow
E_{++}$ and $E_{-+}\rightarrow E_{++}$, respectively. The gap energies which
contribute to the onset of these transition peaks are indicated in Fig. 7 and
given analytically in table 1. Also, the conductivity retains the flat
absorption at sufficiently higher values of $\hslash\omega$ similar to
pristine graphene [57].

Plots of the real part of $\sigma_{yx}^{v}$ for $E_{F}=0$ ( black dotdashed
curve) and $E_{F}\neq 0$ (red and blue dashed curves) in the absence of Rashba
SOC ($\lambda_{R}=0$) are shown in Fig. 11. In the dc limit, the expected
value of the valley-Hall conductivity is obtained as can be seen in Fig. 5
(black curve). If the system is illuminated by photons of frequency $\omega$,
the amplitude of the absorption peaks is suppressed for $E_{F}=0$, while an
increase in it is observed for $E_{F}\neq 0$. For $\hslash\omega=2|\lambda|$ a
strong valley-Hall response is observed for $E_{F}\neq 0$. Therefore, it can
be expected that a stronger valley-Hall response may be accessible when the
photon energy is tuned to the valley-Zeeman SOC. For
$\hslash\omega>2|\lambda|$, $\sigma_{yx}^{v}$ decreases rapidly and approaches
zero at sufficiently higher values of $\hslash\omega$.

Figure 13:  Longitudinal conductivity $\sigma^{d}_{xx}$ in units of $e^{2}/h$
versus electron concentration $n_{e}$ for different values of $\lambda$ and
$\lambda_{R}$. For further clarity, the range $0-5$ meV is shown in inset.

The real part of the valley-Hall conductivity is shown in Fig. 12 for several
values of $E_{F}$. In the dc limit $(\omega=0)$, we obtain the quantized value
of the valley-Hall conductivity (Re$\sigma_{yx}^{v}=e/h$) for $E_{F}=0$ (black
curve in the upper panel). If the system is subjected to photon of frequency
$\omega$, an increase in the magnitude of the valley-Hall response is
observed. The absorption peaks occur at the same onset energies as indicated
in Fig. 7. For example, the first peak appeared when
$\hslash\omega=2\Delta_{1}$ or transition between the minima of the $E_{--}$
and $E_{+-}$ bands. Further, the change in sign of the conductivity is due to
the Rashba SOC, which is responsible for the coupling between spin-up and
spin-down bands e.g., the transition from the maximum of mexican hat of
$E_{--}$ band to the minimum of $E_{++}$ band around $k=0$. Furthermore, for
finite values of $E_{F}$ we obtain new features in the optical spectrum due to
the emergence of new transitions such as $E_{+-}\rightarrow E_{++}$, e.g.,
some features are completely removed due to Pauli blocking. Also, the valley-
Hall response is diminished at sufficiently high frequencies. However, in the
case of $\lambda_{R}<\lambda$ (lower panel), the difference among the optical
transition energies is significantly enhanced due to larger values of
$\lambda$ and new features emerge at the momenta at which $E_{F}$ crosses the
$E_{+-}$ band (see Fig. 7). Moreover, some of the optical transitions are no
longer possible, e.g., $E_{--}\rightarrow E_{+-}$ when $E_{F}$ is just above
the mexican hat because the states below it are occupied and, therefore, Pauli
blocked (blue curve).

In Fig. 13 we plot $\sigma_{xx}^{d}$, from Eq. (21), by evaluating it
numerically versus electron concentration ($n_{e}$) and using the expression
of $\tau$ given in Appendix A but evaluated at the Fermi level, $k=k_{F}$. The
conductivity increases with $E_{F}$ and therefore with the carrier density
$n_{e}$. The diffusive conductivity increases linearly with $n_{e}$ but cusp-
like features appear when $E_{++}$ band begin to occupied at specific values
of $n_{e}$ in contrast to pristine graphene [61, 62]. This behaviour makes
graphene/WS2 a suitable candidate for charge switches contrary to pristine
graphene. The screening effect becomes significantly weaker when only the
$\lambda$ term is present. Moreover, the conductivity shown in Fig. 14
increases in the low-density regime for $\lambda=0$ and $\lambda_{R}\neq 0$ as
compared to the $\lambda\neq 0$, $\lambda_{R}=0$ and $\lambda,\lambda_{R}\neq
0$ case. In the limit $\lambda=\lambda_{R}=0$ we obtain the result similar to
pristine graphene [61, 62].

Figure 14:  Same as in Fig. 13 but for different values of $\lambda_{R}$.

##  IV summary and conclusion

We studied the energy dispersion of graphene/WSe2 heterostructures by using a
TB model in the presence of valley-Zeeman and Rashba SOCs. We found that the
effective Hamiltonian (2) derived from the TB one (1) nicely captures the low-
energy physics near the $K$ and $K^{\prime}$ valleys. We demonstrated that the
density of states has a finite value around $E=0$ in both cases $\lambda\neq
0,\lambda_{R}=0$ and $\lambda=0,\lambda_{R}\neq 0$. In addition, it has a
square root singularity when both $\lambda$ and $\lambda_{R}$ are present.
This singularity is similar to that in biased bilayer graphene; however, here
it is due to the Rashba SOC whereas in biased bilayer graphene it is due to
interlayer hopping. We also found that the ac and dc valley-Hall
conductivities change sign in the presence of the $\lambda_{R}$ term, which
leads to interband transitions. Also, the band gap is enhanced by increasing
the strength $\lambda_{R}$. Further, for $\lambda_{R}>>\lambda$ the valley-
Hall conductivity exhibits a behaviour similar to that in gapped graphene and
topological insulators [54, 58]. The screening effect in the diffusive
conductivity is dominant only when the Rashba SOC is present, whereas it is
significantly suppressed for $\lambda\neq 0,\lambda_{R}=0$. Also, the
conductivity increases with $\lambda_{R}$ in the low- and high-density
regimes, see Fig. 14.

The dc valley-Hall conductivity changes sign when $E_{F}$ is comparable to
$\lambda_{R}$ and vanishes at higher values of $E_{F}$, cf. Fig. 5. It also
exhibits a strong temperature dependence when the Fermi level in the gap, cf.
Fig. 6.

The intraband response of the ac longitudinal conductivity for $\lambda_{R}=0$
(see upper panel of Fig. 9) shifts towards lower photon energies when $E_{F}$
increases compared to $\lambda_{R}\neq 0$ (see lower panel of Fig. 9 and Fig.
10). We also noted the switching on and off of the Drude response when the
Fermi energy is varied (see Fig. 10), which may be of interest in
technological applications. In addition, for $\lambda,\lambda_{R}\neq 0$ new
onsets in the optical conductivity appear due to the shifting of the Fermi
level through the mexican hat structure (see Figs. 10 and 12), which may be a
promising feature in optical experiments. Our findings may be pertinent to
developing future spintronics and valleytronics devices such as field-effect
tunnelling transistors, memory devices, phototransistors, etc.

###### Acknowledgements.

M. Z. and P. V. acknowledge the support of the Concordia University Grant No.
VB0038 and a Concordia University Graduate Fellowship. The work of M. T. was
supported by Colorado State University. *

##  Appendix A Relaxation time

The relaxation time is generally a function of the incoming electron’s wave
vector and at low temperatures only states near the Fermi level will
contribute to transport and single-particle properties. Below we provide
expressions for the relaxation time at the Fermi energy in the limiting cases
$\Delta,\lambda\neq 0,\lambda_{R}=0$ and $\Delta,\lambda=0,\lambda_{R}\neq 0$,
because in these cases the summation over final states can be performed
analytically. Within the first Born approximation the standard formula for the
momentum relaxation time has the form

$\displaystyle\dfrac{1}{\tau_{\zeta}}=\dfrac{1}{\tau_{\xi
k}^{\eta}}=\dfrac{2\pi
n_{i}}{\hslash}\sum_{\xi^{\prime},\eta^{\prime},k^{\prime}}|\left\langle\xi,\eta,k\right|U(\mathbf{r})\left|\xi^{\prime},\eta^{\prime},k^{\prime}\right\rangle|^{2}\delta(E_{\xi
k}-E_{\xi^{\prime}k^{\prime}})(1-\cos\theta),$ (22)

where $U(\mathbf{r})$ is the impurity potential, $n_{i}$ the impurity density,
and $\theta$ the angle between the initial $k$ and final $k^{\prime}$ wave
vectors. Equation (22) holds only for elastic scattering
$(\xi=\xi^{\prime},\eta=\eta^{\prime},k=k^{\prime})$ and for central
potentials $U(\mathbf{r})$ i.e. $U(\mathbf{r})=U(r)$. The results for two
types of impurity potentials are as follows.

Short-range impurities. We have
$U(\mathbf{r})=U_{0}\delta(\mathbf{r}-\mathbf{r_{i}})$ where $\mathbf{r}$ and
$\mathbf{r_{i}}$ are the position vectors of the electron and impurity,
respectively, and $U_{0}$ is the strength of potential. In this case
$U(\mathbf{q})=U_{0}$ is the Fourier transform of
$U(\mathbf{r})=(1/\sqrt{L_{x}L_{y}})\sum_{q}U(q)e^{i\mathbf{q.r}}$ with
$|\mathbf{q}|=2k\sin(\theta/2)$. The results are:

i) $\lambda_{R}=0$

$\displaystyle\dfrac{1}{\tau_{sk_{F}}^{\eta}}=\dfrac{V_{0}^{2}n_{i}(N_{s}^{\eta})^{4}}{\hslash}\dfrac{\sqrt{\Delta^{2}+(\hslash
v_{F}k)^{2}}}{(\hslash
v_{F})^{2}}\,\big{[}(A_{s}^{\eta})^{4}-(A_{s}^{\eta})^{2}+1\big{]}.$ (23)

In the limit $\Delta,\lambda=0$, the above result reduces to graphene’s
scattering time Eq. (24) of Ref. [61]

$\dfrac{1}{\tau_{k_{F}}}=\dfrac{V_{0}^{2}n_{i}k}{4\hslash^{2}v_{F}}.$ (24)

Also, for $\lambda=0$ Eq. (23) agrees with the result for topological
insulators [58].

ii) $\Delta,\lambda=0$

$\displaystyle\dfrac{1}{\tau_{sk_{F}}^{\eta}}$ $\displaystyle=$
$\displaystyle\dfrac{V_{0}^{2}n_{i}(N_{s}^{\eta})^{4}\sqrt{\lambda_{R}^{2}+(\hslash
v_{F}k)^{2}}}{\hslash^{3}v_{F}^{2}}\,\Big{[}[(A_{s}^{\eta})^{2}+(B_{s}^{\eta}]^{2}]^{2}+(C_{s}^{\eta})^{4}-[1+(C_{s}^{\eta})^{2}][(A_{s}^{\eta})^{2}+(B_{s}^{\eta})^{2}]$
$\displaystyle\hskip 128.0374pt-2(C_{s}^{\eta})^{2}+1\Big{]}.$

Long-range impurities. We assume
$U(\mathbf{r})=eQe^{-k_{s}r}/4\pi\epsilon_{0}\epsilon r$, where $k_{s}$ is the
screening wave vector, $Q$ is the charge of the impurity, and $\epsilon$ the
dielectric constant. In this case $U(q)=2\pi U_{0}/\sqrt{k_{s}^{2}+q^{2}}$
with $U_{0}=eQ/4\pi\epsilon_{0}\epsilon$. The results are:

i) $\lambda_{R}=0$

$\displaystyle\dfrac{1}{\tau_{sk_{F}}^{\eta}}$ $\displaystyle=$
$\displaystyle\dfrac{V_{0}^{2}n_{i}(N_{s}^{\eta})^{4}\sqrt{\Delta^{2}+(\hslash
v_{F}k)^{2}}}{2\hslash^{3}v_{F}^{2}k^{2}}\,\Big{[}[1+(A_{s}^{\eta})^{4}]\big{[}1-\dfrac{a_{s}}{\sqrt{a_{s}^{2}+1}}\big{]}+2(A_{s}^{\eta})^{2}\big{[}2a_{s}^{2}-\dfrac{a_{s}(2a_{s}^{2}+1)}{\sqrt{a_{s}^{2}+1}}\big{]}\Big{]}.$
(26)

In the limit $\Delta=\lambda=0$ we set $a_{s}=k_{s}/2k$ and obtain the
relaxation time in pristine graphene [63]

$\dfrac{1}{\tau_{k_{F}}}=\dfrac{V_{0}^{2}n_{i}\,(a_{s}-\sqrt{a_{s}^{2}+1})^{2}}{4\hslash^{2}v_{F}k}.$
(27)

Moreover, for $\lambda=0$ Eq. (26) gives the relaxation time for topological
insulators [58].

ii) $\Delta,\lambda=0$

$\displaystyle\dfrac{1}{\tau_{sk_{F}}^{\eta}}$ $\displaystyle=$
$\displaystyle\dfrac{V_{0}^{2}n_{i}(N_{s}^{\eta})^{4}\sqrt{\lambda_{R}^{2}+(\hslash
v_{F}k)^{2}}}{2\hslash^{3}v_{F}^{2}k^{2}}\,\Big{[}1+[(A_{s}^{\eta})^{2}+(B_{s}^{\eta})^{2}]^{2}+(C_{s}^{\eta})^{4}[1-\dfrac{a_{s}}{\sqrt{a_{s}^{2}+1}}]$
$\displaystyle+2[1+(C_{s}^{\eta})^{2}][A_{s}^{\eta})^{2}+(B_{s}^{\eta})^{2}]\big{[}2a_{s}^{2}-\dfrac{a_{s}(2a_{s}^{2}+1)}{\sqrt{a_{s}^{2}+1}}\big{]}$
$\displaystyle-2(C_{s}^{\eta})^{2}\big{[}1+\dfrac{a_{s}}{\sqrt{a_{s}^{2}+1}}+8a_{s}^{3}\sqrt{a_{s}^{2}+1}-8a_{s}^{2}(2a_{s}^{2}+1)\big{]}\Big{]}.$

## References

  * [1] A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * [2] N. Tombros, C. Jozsa, M. Popinciuc, H. T. Jonkman, and B. J. van Wees, Nature 448, 571 (2007).
  * [3] J. Ingla-Aynés, M. H. D. Guimarães, R. J. Meijerink, P. J. Zomer, and B. J. van Wees, Phys. Rev. B 92, 201410(R) (2015).
  * [4] M. Drögeler, C. Franzen, F. Volmer, T. Pohlmann, L. Banszerus, M. Wolter, K. Watanabe, T. Taniguchi, C. Stampfer, and B. Beschoten, Nano Lett. 16, 3533 (2016).
  * [5] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * [6] Z. Qiao, S. A. Yang, W. Feng, W.-K. Tse, J. Ding, Y. Yao, Y. Wang and Q. Niu, Phys. Rev. B 82, 161414 (2010); C.-X. Liu, S.-C. Zhang and X.-L. Qi, Annu. Rev. Condens. Matter Phys. 7, 301 (2016); Y. Ren, Z. Qiao and Q. Niu, Rep. Prog. Phys. 79, 066501 (2016); H. Weng, R. Yu, X. Hu, X. Dai and Z. Fang, Adv. Phys. 64, 227 (2015).
  * [7] A. H. C Neto, and F. Guinea, Phys. Rev. Lett. 103, 026804 (2009).
  * [8] C. Weeks, J. Hu, J. Alicea, M. Franz and R. Wu, Phys. Rev. X 1, 021001 (2011).
  * [9] J. Ding, Z. Qiao, W. Feng, Y. Yao and Q. Niu, Phys. Rev. B 84, 195444 (2011).
  * [10] J. Hu, J. Alicea, R. Wu and M. Franz, Phys.Rev. Lett. 109, 266801 (2012).
  * [11] D. Ma, Z. Li and Z. Yang, Carbon 50, 297 (2012).
  * [12] K.-H. Jin and S.-H. Jhi, Phys. Rev. B 87, 075442 (2013).
  * [13] A. Ferreira, T. G. Rappoport, M. A. Cazalilla and A. H. C. Neto, Phys. Rev. Lett. 112, 066601 (2014).
  * [14] A. A. Kaverzin and B. J. van Wees, Phys. Rev. B 91, 165412 (2015).
  * [15] J. Balakrishnan, G. K. W. Koon, M. Jaiswal, A. H. C. Neto and B. C. Ozyilmaz, Nat. Phys. 9, 284 (2013).
  * [16] X. Hong, S.-H. Cheng, C. Herding, and J. Zhu, Phys. Rev. B 83, 085410 (2011).
  * [17] Z. Jia, B. Yan, J. Niu, Q. Han, R. Zhu, D. Yu and X. Wu, Phys. Rev. B 91, 085411 (2015).
  * [18] U. Chandni, E. A. Henriksen and J. P. Eisenstein, Phys. Rev. B 91, 245402 (2015).
  * [19] Y.-C. Lin, N. Lu, N. Perea-Lopez, J. Li, Z. Lin, X. Peng, C. H. Lee, C. Sun, L. Calderin, P. N. Browning, M. S. Bresnehan, M. J. Kim, T. S. Mayer, M. Terrones, and J. A. Robinson, ACS Nano 8, 3715 (2014).
  * [20] M.-Y. Lin, C.-E. Chang, C.-H. Wang, C.-F. Su, C. Chen, S.-C. Lee, and S.-Y. Lin, Appl. Phys. Lett. 105, 073501 (2014).
  * [21] A. Azizi, S. Eichfeld, G. Geschwind, K. Zhang, B. Jiang, D. Mukherjee, L. Hossain,A. F. Piasecki, B. Kabius, J. A. Robinson, and N. Alem, ACS Nano 9, 4882 (2015).
  * [22] Y. Kim, D. Choi, W. J. Woo, J. B. Lee, G. H. Ryu, J. H. Lim, S. Lee, Z. Lee, S. Im, J.-H. Ahn, W.-H. Kim, J. Park, and H. Kim, Appl. Surf. Sci. 494, 591 (2019).
  * [23] A. M. Alsharari, M. M. Asmar, and S. E. Ulloa, Phys. Rev. B 98, 195129 (2018).
  * [24] L. A. Benitez, J. F. Sierra, W. S. Torres, A. Arrighi, F. Bonnell, M. V. Costache, and S. O. Valenzuera, Nat. Phys. 14, 303 (2018).
  * [25] A. W. Cummings, J. H. Garcia, J. Fabian, and S. Roche, Phys. Rev. Lett. 119, 206601 (2017).
  * [26] W. Han, R. K. Kawakami, M. Gmitra, and J. Fabian, Nat. Nanotechnol. 9, 794 (2014).
  * [27] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Phys. Rev. Lett. 105, 136805 (2010).
  * [28] A. Kormányos, G. Burkard, M. Gmitra, J. Fabian, V. Zólyomi, N. D. Drummond, and V. Falḱo, 2D Mater. 2, 022001 (2015).
  * [29] C.-P. Lu, G. Li, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Phys. Rev. Lett. 113, 156804 (2014).
  * [30] S. Larentis, J. R. Tolsma, B. Fallahazad, D. C. Dillen, K. Kim, A. H. MacDonald, and E. Tutuc, Nano Lett. 14, 2039 (2014).
  * [31] L. Banszerus,T. Sohier, A. Epping,F. Winkler, F. Libisch, F. Haupt, K. Watanabe, T. Taniguchi, K. Muller-Caspary, N. Marzari, F. Mauri, B. Beschoten, and C. Stampfer, arXiv: 1909.09523.
  * [32] S. Bertolazzi, D. Krasnozhon, and A. Kis, ACS Nano 7, 3246 (2013).
  * [33] K. Roy, M. Padmanabhan, S. Goswami, T. P. Sai, G. Ramalingam, S. Raghavan, and A. Ghosh, Nat. Nanotechnol. 8, 826 (2013).
  * [34] W. Zhang, C.-P. Chuu, J.-K. Huang, C.-H. Chen, M.-L. Tsai, Y.-H. Chang, C.-T. Liang, Y.-Z. Chen, Y.-L. Chueh, J.-H. He, M.-Y. Chou, and L.-J. Lib, Sci. Rep. 4, 3826 (2014).
  * [35] N. A. Kumar, M. A. Dar, R. Gul, and J. Baek, Mater. Today 18, 286 (2015).
  * [36] L. Britnel, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Schedin, A. Mishchenko, T. Georgiou, M. I. Katsnelson, L. Eaves, S. V. Morozov, N. M. R. Peres, J. Leist, A. K. Geim1, K. S. Novoselov and L. A. Ponomarenko, Science 335, 947 (2012).
  * [37] A. Mishchenko, J. S. Tu, Y. Cao, R. V. Gorbachev, J. R. Wallbank, M. T. Greenaway, V. E. Morozov, S. V. Morozov, M. J. Zhu, S. L. Wong, F. Withers, C. R. Woods, Y-J. Kim, K. Watanabe, T. Taniguchi, E. E. Vdovin, O. Makarovsky, T. M. Fromhold, V. I. Fal’ko, A. K. Geim, L. Eaves and K. S. Novoselov, Nat. Nanotechnol. 9, 808 (2014).
  * [38] K. Roy, M. Padmanabhan, S. Goswami, T. P. Sai, G. Ramalingam, S. Raghavan and A. Ghosh, Nat. Nanotechnol. 8, 826 (2013).
  * [39] A. David, P. Rakyta, A. Kormányos, and Guido Burkard, Phys. Rev. B 100, 085412 (2019).
  * [40] Y. Li and M. Koshino, Phys. Rev. B 99, 075438 (2019).
  * [41] A. Avsar, J. Y. Tan, T. Taychatanapat, J. Balakrishnan, G. K. W. Koon, Y. Yeo, J. Lahiri, A. Carvalho, A. S. Rodin, E. C. T. O´Farrell, G. Eda, A. H. C. Neto and B. Özyilmaz, Nat. Commun. 5, 4875 (2014).
  * [42] M. Gmitra, S. Konschuh, C. Ertler, C. Ambrosch-Draxl, and J. Fabian, Phys. Rev. B 80, 235431 (2009).
  * [43] C. K. Safeer, J. Ingla-Aynés, F. Herling, J. H. Garcia, M. Vila, N. Ontoso, M. Reyes Calvo, S. Roche, L. E. Hueso, and F. Casanova, Nano Lett. 19, 1074 (2019).
  * [44] B. Yang, M.-F. Tu, J. Kim, Y. Wu, H. Wang, J. Alicea, R. Wu, M. Bockrath and J. Shi1, 2D Mater. 3, 031012 (2016).
  * [45] S. Zihlmann, A. W. Cummings, J. H. Garcia, M. Kedves, K. Watanabe, T. Taniguchi, C. Schonenberger, and P. Makk, Phys. Rev. B 97, 075434(R) (2018).
  * [46] Z. Wang, D. K. Ki, H. Chen, H. Berger, A. H. MacDonald, and A. F. Morpurgo, Nat. Commun. 6, 8339 (2015).
  * [47] Jose H. Garcia, Marc Vila, Aron W. Cummings, and Stephan Roche, Chem. Soc. Rev. 47, 3359 (2018); A. Mreńca-Kolasińska, B. Rzeszotarski, and B. Szafran, Phys. Rev. B 98, 045406 (2018).
  * [48] T. Völkl, T. Rockinger, M. Drienovsky, K. Watanabe, T. Taniguchi, D. Weiss, and J. Eroms, Phys. Rev. B 96, 125405 (2017); B. Yang, E. Molina, J. Kim, D. Barroso, M. Lohmann, Y. Liu, Y. Xu, R. Wu, L. Bartels, K. Watanabe, T. Taniguchi, and Jing Shi, Nano Lett. 18, 3580 (2018).
  * [49] M. Gmitra, D. Kochan, P. Högl, and J. Fabian, Phys. Rev. B 93, 155104 (2016).
  * [50] M. Gmitra and J. Fabian, Phys. Rev. B 92, 155403 (2015).
  * [51] M. Gmitra, D. Kochan, and J. Fabian, Phys. Rev. Lett. 110, 246602 (2013).
  * [52] Z. Wang, D.-K. Ki, J. Y. Khoo, D. Mauro, H. Berger, L. S. Levitov, and A. F. Morpurgo1, Phys. Rev. X 6, 041020 (2016).
  * [53] M. Charbonneau, K. M. Van Vliet, and P. Vasilopoulos, J. Math. Phys. 23, 318 (1982).
  * [54] D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. 99, 236809 (2007).
  * [55] A. Rycerz, J. Tworzydlo, and C. W. J. Beenakker, Nat. Phys. 3, 172 (2007).
  * [56] Z. Li and J. P. Carbotte, Phys. Rev. B 86, 205425 (2012).
  * [57] V. P. Gusynin, S. G. Sharapov, and J. P. Carbotte, Phys. Rev. Lett. 96, 256802 (2006).
  * [58] V. Vargiamidis and P. Vasilopoulos, J. Appl. Phys. 116, 063713 (2014).
  * [59] D. S. L. Abergel and V. I. Fal´ko, Phys. Rev. B 75, 155430 (2007).
  * [60] E. J. Nicol and J. P. Carbotte, Phys. Rev. B 77, 155409 (2008).
  * [61] T. Stauber, N. M. R. Peres, and F. Guinea, Phys. Rev. B 76, 205423 (2007).
  * [62] K. Nomura and A. H. MacDonald, Phys. Rev. Lett. 98, 076602 (2007).
  * [63] A. A. Patel and S. Mukerjee, Phys. Rev. B 86, 075411 (2012).