# Meteorological factors and non-pharmaceutical interventions explain local
differences in the spread of SARS-CoV-2 in Austria

Katharina Ledebur1,2 Michaela Kaleta1,2 Jiaying Chen1,2,3 Simon Lindner1,2
Caspar Matzhold1,2 Florian Weidle4 Christoph Wittmann4 Katharina Habimana5
Linda Kerschbaumer5 Sophie Stumpfl5 Georg Heiler2,6 Martin Bicher6,7 Nikolas
Popper6,7,8 Florian Bachner5 Peter Klimek1,2,∗

###### Abstract

The drivers behind regional differences of SARS-CoV-2 spread on finer spatio-
temporal scales are yet to be fully understood. Here we develop a data-driven
modelling approach based on an age-structured compartmental model that
compares 116 Austrian regions to a suitably chosen control set of regions to
explain variations in local transmission rates through a combination of
meteorological factors, non-pharmaceutical interventions and mobility. We find
that more than 60% of the observed regional variations can be explained by
these factors. Decreasing temperature and humidity, increasing cloudiness,
precipitation and the absence of mitigation measures for public events are the
strongest drivers for increased virus transmission, leading in combination to
a doubling of the transmission rates compared to regions with more favourable
weather. We conjecture that regions with little mitigation measures for large
events that experience shifts toward unfavourable weather conditions are
particularly predisposed as nucleation points for the next seasonal SARS-CoV-2
waves.

Medical University of Vienna, Section for Science of Complex Systems, CeMSIIS,
Spitalgasse 23, 1090 Vienna, Austria;

Complexity Science Hub Vienna, Josefstädter Straße 39, 1080 Vienna, Austria;

Division of Insurance Medicine, Department of Clinical Neuroscience,
Karolinska Institutet, SE-171 77, Stockholm, Sweden;

Zentralanstalt für Meteorologie und Geodynamik, Hohe Warte 38, 1190 Vienna,
Austria;

Austrian National Public Health Institute, Stubenring 6, A-1010 Vienna,
Austria;

Institute of Information Systems Engineering, TU Wien, Favoritenstraße 8-11,
A-1040 Vienna, Austria;

dwh simulation services, dwh GmbH, Neustiftgasse 57-59, A-1070 Vienna,
Austria;

Association for Decision Support Policy and Planning, DEXHELPP, Neustiftgasse
57-59, A-1070 Vienna, Austria. ∗Correspondence<EMAIL_ADDRESS>

##  1 Introduction

The SARS-CoV-2 pandemic impacted different world regions in a highly
heterogeneous way. A vast body of research has sought to explain these
differences by factors ranging from varying stringency of non-pharmaceutical
interventions (NPIs) [1, 2, 3, 4], socio-economic and demographic factors [5,
6, 7] to different virus variants. [8, 9] Relatively less attention has been
paid to more fine-scaled spatio-temporal variations in virus spread. [10, 11]
Next to regional differences in NPIs, meteorological [12, 13, 14, 15, 17, 16]
and behavioural factors have been proposed to account for such variations.
[18, 19, 20, 21] Infection waves typically start with localized outbreaks in
specific regions before case numbers start to soar on larger geographic
scales. [22] Such localized outbreaks are often initiated by singular events
in which transmission is dominated by a small number of individuals, [23, 24]
so-called superspreading events. [25] As these events are stochastic and
therefore near-impossible to predict, the question arises to which extent
fine-scaled spatio-temporal variations can actually be explained or whether
they are irreducibly random.

Here we seek to understand the factors that determine on finer geographic
scales the degree to which the infection dynamics in one region in Austria
deviated from the dynamics observed in neighboring regions. We consider
factors from three different domains to account for these variations, namely
meteorology (temperature, cloudiness, humidity, precipitation, wind), non-
pharmaceutical interventions (targeting schools, gastronomy, healthcare, or
events) as well as individual-level mobility (as inferred from
telecommunications data).

Previous literature reported inconsistent associations of weather with
transmission dynamics, [12, 13, 14, 15, 17, 16] suggesting that the assessment
of the impact of meteorological factors might depend on geography, epidemic
phase and the spatio-temporal scale on which the analysis is conducted. In
this work, we use meteorological forecast data derived from the mesoscale
numerical weather model AROME (Application of Research to Operations at
Mesoscale [26, 27]) to quantify regional weather related factors. AROME is
operated by ZAMG (Zentralanstalt für Meteorologie und Geodynamik) for a domain
covering the Alpine area. More details are given in section. 2.2

The literature is more consistent with regard to the transmission-reducing
effects of restricted mobility. [18, 19, 20, 21] However, such effects cannot
be fully disentangled from the effects of physical distancing policies. [28]
For Austria it was observed that the same regime of governmental interventions
induced quite heterogeneous mobility changes across regions. [29] Here we
measure mobility by means of the median Radius of Gyration of cell phone users
in a region, which is a measure for the daily area travelled by these users.

In autumn 2020, Austria adopted a tiered pandemic management plan in which
each region was assigned one of four alert levels (green, yellow, orange, red)
signalling its epidemiological risk. [30] The assessment was performed based
on a set of indicators including confirmed cases, available hospital
capacities and the functionality of contact tracing. Regional authorities
decided on appropriate response measures in reaction to the weekly risk
assessment, meaning that the same alert level could and did translate into
different NPIs across regions. We used an extensive dataset on all regional
response measures in Austria curated by the Austrian National Public Health
Institute. We clustered the response measures into four different categories,
based on whether they affect schools (bans on singing in class rooms, cloth
masks when not seated, …), gastronomy (limits for the number of persons
allowed to share a table, opening hours, …), healthcare settings (e.g.,
visitor bans), or events (size restrictions).

Our methodological approach is outlined in Fig. 1. There we show a map of
Austria divided into its nine federal states. The state of Tyrol (blue) is
further divided into its finer administrative regions of districts. In Fig. 1
we show the epidemiological curve observed in the district of Innsbruck
(orange) and compare it with two models. First, we compare it to results from
an age-structured compartmental model that assumes the infection rate observed
in Innsbruck (orange) to be the same as the rate observed in all other
districts taken together in the same federal state (blue, the control set for
the orange district), giving the red epidemiological curve. Second, we augment
this model by assuming that the age-dependent transmission rates are also a
function of the district-level meteorological, intervention and mobility time-
series with district-independent effect sizes learned from data by means of a
cross-validated nonlinear least squares solver. This approach provides us with
estimates for the effect of each of the ten input time series on regional
virus spread and for the degree to which the combined effects of weather,
interventions and individual-level mobility can be used to explain differences
in the spread between regions.

Figure 1: Visual representation of the methodological approach. For every one
of the 116 districts in Austria the epidemiological curve of the district is
compared to the epidemiological curve of the corresponding federal state. This
figure shows the example of the district Innsbruck (orange) in the federal
state Tyrol. The epidemiological curves are calculated employing an age-
structured compartmental model. The red curve is the epidemiological curve for
Innsbruck assuming transmission rates that equal those observed in Tyrol
(blue). By including a dependence of the transmission rate on weather,
interventions and mobility (green), the district-independent effect sizes of
these district-specific input variables can be calculated.

##  2 Data and Methods

###  2.1 Case data.

Our study period ranges from July 1 2020 to May 15 2021. The daily confirmed
Covid-19 cases of 116 districts in Austria were collected through the official
Austrian COVID-19 disease reporting system (EMS) [31]. A unique identification
number was assigned to each positive Covid-19 test. People can be tested
positive multiple times. We calculate daily numbers of positive tests per
district and age group. The population size per district was available from
the national statistics office and linked by the unique district ID.

###  2.2 Meteorological data.

We further used detailed information on the weather situation in Austria,
including gridded forecast data on five meteorological parameters: total
cloudiness, precipitation, 2m temperature, 2m relative humidity and 10m wind
speed. The gridded data sets were provided by the mesoscale weather model
AROME operated at ZAMG eight times per day producing forecasts up to 60 hours
ahead. AROME is a spectral limited area model currently operating on a 2.5km
horizontal grid covering the Alpine area using 90 vertical levels up to a
height of approximately 35hPa. The 3D atmospheric initial conditions for the
model runs are created using a variational data assimilation system (3DVAR,
[32]) which is fed by a large number of observations (surface station data,
radiosondes, satellite data, aircraft data, etc.). To derive the surface
initial conditions, an optimum interpolation method is implemented.

For this study, the meteorological model data is provided for the study period
on an hourly basis. The gridded AROME data was further aggregated to create
information on a district level, where the district value of an input time
series is given by the population-weighted average of the values observed in
the corresponding grid cells. We calculate the daily mean values of these
parameters in each district and subtract the mean value of the daily mean
values of all other districts in the same federal state (control set) as input
for the model.

###  2.3 Data on regional NPIs.

A dataset on regional Austrian NPIs was collected and curated by the Austrian
National Public Health Institute. This data contains precise information on
federal state, district, type and details of intervention, start and end date
of registered regional measures. The data set contains 11 individual types and
an additional 83 subtypes of interventions. We further grouped these types and
subtypes of NPIs into 4 larger groups of restrictions: 1) restrictions in
schools, 2) restrictions in gastronomy, 3) restrictions in healthcare, 4)
restrictions for larger events. Intervention measures dependent on the corona
traffic light were matched to the districts in appropriate traffic light
colors in the given time frame. For each day in the observed time period and
each district we assigned binary values of true/false for each intervention
group.

###  2.4 Mobility data.

The mobility data contains 1,683,270 entries, including information on the
date of travel, regions by postal code (which were mapped to districts), the
number of devices in this particular region, and the radius of gyration over
the study period. For each region only k-anonymized data is available. Any
region with less than k unique devices on a specific day is redacted. The
radius of gyration is defined as the square root of the time-weighted centroid
of the squared distances of the mobility events (locations) $\vec{x}_{i\mu}$
to the most prominent location (night location, center of mass)
$\overline{x}_{i}=\frac{\sum_{\mu}\vec{x}_{i\mu}t_{i\mu}}{\sum_{\mu}t_{i\mu}}$
of each day:

$R_{\mathrm{G},i}=\sqrt{\frac{\sum_{\tau}d(\overline{x}_{i},\vec{x}_{i\tau})^{2}}{\sum_{\tau}t_{i\tau}}}$
(1)

Finally, as input time series for the model we consider the logarithmic radius
of gyration and subtract the average logarithmic radius of gyration observed
in the control set.

###  2.5 Regional age-structured SIR null model without effects

To identify the impact of meteorological factors, mitigation strategies and
mobility on the infection dynamics across Austria, the epidemic curves of the
116 districts of Austria are compared. We employ a parsimonious age-structured
SIR model [33] to facilitate robust calibration.

Our approach makes use of the fact that when the confirmed cases per day are
known, the “effective” transmission rate $\alpha$ can be calculated from data
alone within the SIR model. For each district $b$, we compute a transmission
rate for day $n$ and age group $a$ in the null model (without effects),
$\alpha_{n,a}^{0}(b)$, as follows. Our null hypothesis is that the
epidemiological curve of a given district $b$ (shown as orange in Fig. 1) has
values of $\alpha_{n,a}^{0}(b)$ that are identical to the transmission rate
observed in all other districts of the same federal state (blue in Fig. 1). We
call these districts the control set of districts for $b$, $\mathcal{C}_{b}$.
For this control set we observe the time series of cumulative confirmed cases,
$C_{n,a}(\mathcal{C}_{b})$.

Let $S(\mathcal{C}_{b})$, $I(\mathcal{C}_{b})$, $R(\mathcal{C}_{b})$, and
$N(\mathcal{C}_{b})$ be the daily numbers of susceptible, infected, recovered
and total number of people in the control set, respectively. We have,

$\alpha_{n,a}(\mathcal{C}_{b})=\frac{\left(C_{n+1,a}(\mathcal{C}_{b})-C_{n,a}(\mathcal{C}_{b})\right)N_{a}(\mathcal{C}_{b})}{S_{n,a}(\mathcal{C}_{b})\sum_{a^{\prime}}c_{aa^{\prime}}I_{n,a^{\prime}}(\mathcal{C}_{b})}\quad,$
(2)

where $c_{aa^{\prime}}$ represents social mixing by age. We compute a social
mixing matrix for every district to define the number of infected individuals
a susceptible individual from one age group is exposed to. The social mixing
matrix by age for Austria is obtained from Prem et al. (2017) [34] with the
population being grouped into four age brackets (0-19, 20-39, 40-64 and 65+
year old individuals). We calculate a social mixing matrix for every district
and federal state. To calculate the social mixing matrix entry for one of the
four new age groups, we build the population-weighted sum over the
corresponding social mixing matrix entries provided by Prem et al. (2017). The
social mixing matrices for each federal state (used for the control sets) are
again given by the population-weighted sum over the individual districts.

The number of susceptible, infected and removed individuals in one of the four
age groups $a$, can be calculated for the control set as,

$\displaystyle S_{n+1,a}(\mathcal{C}_{b})$ $\displaystyle=$ $\displaystyle
S_{n,a}(\mathcal{C}_{b})\left(1-\lambda_{n,a}(\mathcal{C}_{b})\right)$ (3)
$\displaystyle I_{n+1,a}(\mathcal{C}_{b})$ $\displaystyle=$ $\displaystyle
I_{n,a}(\mathcal{C}_{b})+\lambda_{n,a}(\mathcal{C}_{b})S_{n,a}(\mathcal{C}_{b})-\beta
I_{n,a}(\mathcal{C}_{b})$ (4) $\displaystyle R_{n+1,a}(\mathcal{C}_{b})$
$\displaystyle=$ $\displaystyle R_{n,a}(\mathcal{C}_{b})+\beta
I_{n,a}(\mathcal{C}_{b})\quad,$ (5)

where we defined
$\lambda=\frac{\alpha_{n,a}(\mathcal{C}_{b})}{N_{a}(\mathcal{C}_{b})}\sum_{a^{\prime}}c_{aa^{\prime}}I_{n,a^{\prime}}(\mathcal{C}_{b})$.
Initial conditions are given by $I_{0,a}$ as the day with the first case(s) in
age group $a$ and $S_{0,a}=1-I_{0,a}$.

The null model (without effects) for district $b$ can then be obtained by
assuming that

$\alpha_{n,a}^{0}(b)=\alpha_{n,a}(\mathcal{C}_{b})\ \forall\ n,a,b\quad,$ (6)

holds for an age-structured SIR model for district $b$.

###  2.6 Augmented regional model with effects

To quantify the impact of meteorological factors, mitigation strategies and
overall mobility on the epidemic curves of all districts in Austria, we
augment the null model in the following way.

In total we consider 10 input time series $X_{n}(b,i)$, $i=1,\dots,10$. The
NPI time series contain binary variables, $X_{n}(b,i)\in\\{0,1\\}$, indicating
whether a given measure was implemented in the region or not. For the
meteorological and mobility time series the input signal is given by the
deviation of $X_{n}(b,i)$ between the district and its control set, i.e., we
map $X_{n}(b,i)\to X_{n}(b,i)-\langle
X_{n}(b^{\prime},i)\rangle_{b^{\prime}\in\mathcal{C}_{b}}$ where
$\langle\cdot\rangle_{b^{\prime}\in\mathcal{C}_{b}}$ denotes the arithmetic
average over all districts in the control set. To each input time series we
assign an effect on the transmission rate by assuming that

$\alpha_{n,a}(b)=\alpha_{n,a}^{0}(b)\Pi_{i=1}^{10}(1+\alpha_{a}(i)X_{n}(b,i))\quad,$
(7)

such that $\alpha_{a}(i)$ quantifies by how many percent changes in the input
time series impact on the transmission rate. We further assume that
$\alpha_{a}(i)\equiv\alpha(i)$ (i.e., effects are not age-dependent) for all
input time series except NPIs related to schools, for which we assume an
effect size in the first age group ($<20y$) that is different from all other
age groups. Observe that $\alpha(i)$ are district-independent effects.

###  2.7 Hyperparameter search and cross-validation

We solve the models for $\alpha(i)$ by means of a Levenberg-Marquandt (LM)
algorithm optimizing the residual sum of squares (RSS) between age-specific
incidence time series in data and model; incidence is measured as fraction of
the total population. The model has one remaining hyperparameter, namely the
recovery rate $\beta$. To fix $\beta$, we perform a cross-validated
hyperparameter search over the range $1/\beta=4,\dots,35d$. We perform
randomized cross-validation by splitting the time series into blocks of 28d
and randomly picking 80% of these blocks for training (fitting $\alpha(i)$ via
LM) and evaluating the RSS on the remaining 20% test data. We then choose
$\beta$ such that the difference in RSS between null and augmented model is
maximized by means of the elbow method. To obtain the effect sizes we finally
evaluate the model on the full dataset.

##  3 Results

###  3.1 Model calibration.

Results for the cross-validated hyperparameter search to fix the recovery
time, $\beta$, are shown in Fig. 2. The larger the recovery time, the higher
the percentage of regional variations that the augmented model is able to
explain with respect to the null model. For recovery times of >20d, the
explained variation saturates at a bit more than 60%. Values for the training
data are only marginally higher than for the test data, suggesting that
overfitting is not much of an issue. In the following, we fix $\beta$ to 25d.

Figure 2: Results for the cross-validated hyperparameter search. For different
recovery times, $1/\beta$, we show the percent change in RSS, $\Delta RSS$,
between null and augmented model for training (blue) and test (yellow) data.
For recovery times of more than 20d, the augmented model explains more than
60% of the regional variations in both test and training data.

###  3.2 Effect sizes.

Results for the effect sizes of meteorological, intervention and mobility
variables are shown in Fig. 3. Confidence intervals are obtained from a
100-fold cross-validation. For the continuous variables (weather, mobility),
see Tab. 1, we report the percent change in transmission rates in units of one
standard deviation (SD) of the variable in the district with respect to its
control set of districts. All meteorological variables show a significant
effect. Strong effects can be observed for humidity, where one SD decreases
the transmission rate by $17.1$% ($95$%CI $-17.3$, $-16.9$), as well as for
cloudiness and precipitation, where one SD increases transmission rates by
$15.5$% ($95$%CI $-15.2$, $-15.8$) and $18.9$% ($95$%CI $-18.4$, $-19.5$)
respectively. The transmission rate is also inversely associated with
temperature, while wind slightly increases transmissions. Mobility, as
measured by the logarithmic radius of gyration, increases transmissions by
$7.7$% ($95$% CI $7.4$, $8.0$).

Considering the NPIs, restrictions targeting mass events show the strongest
effect on transmissions with reductions of about $40$%, see also Tab. 2.
Restrictions in healthcare settings and gastronomy reduce transmissions by
about $20\%$, respectively. For school measures we observe age-dependent
effects. While restrictions for schools reduced transmissions by about $8$% in
the population younger than 20y, we observe no significant impact of school
measures on older age groups.

variable | unit | SD | transmission rate change [%]
---|---|---|---
temperature | ${}^{\circ}C$ | $2.4$ |  $-6.9$ ($-7.1$, $-6.8$)
cloudiness | [$0$-$1$] | $0.11$ |  $15.5$ ($15.2$, $15.8$)
humidity | [$0$-$1$] | $0.063$ |  $-17.1$ ($-17.3$, $-16.9$)
precipitation | mm/h | $0.21$ |  $19.0$ ($18.4$, $19.5$)
wind | m/s | $0.95$ |  $3.0$ ($2.9$, $3.1$)
log. radius of gyration | $m$ | $7.9$ |  $7.7$ ($7.4$, $8.0$)
Table 1: Summary of effects of meteorological and mobility time series on the transmission rate. For each variable we give its unit, the standard deviation (SD) of the input time series and the percent change with its 95% confidence intervals of the transmission rate associated with a unit SD change in the input. category | $N$ |  examples |  transmission rate change [%]
---|---|---|---
schools | $144$ |  cloth masks when entering schools, no indoor singing, sports only outdoors, measures to avoid mixing of school classes, … |  $<20y$: $-7.9$ ($-10.6$, $-5.1$)
$\geq 20y$: $0.7$ ($-0.1$, $1.5$)
gastronomy | $123$ |  closing time at 10pm, limits for number of people seated at table, mandatory registration, … |  $-18$ ($-17$, $-19$)
healthcare | $161$ |  visitor ban or a maximum of one visitor per week, mandatory registration, FFP2 masks, … |  $-20.6$ ($-21.2$, $-20.1$)
events | $69$ |  ban or size limits of seated and unseated indoor and outdoor events |  $-37.5$ ($-38.6$, $-36.5$)
Table 2: Summary of effects of NPIs on the transmission rate. For each
category of NPIs we give the number of implementations observed in our data,
list typical examples of what the NPI consists of and the percent change of
the transmission rate associated with a unit SD change in the input. Figure 3:
Summary of effect sizes of NPIs. Impacts on the transmission rate are shown in
percent for weather variables (blue), NPIs (green) and mobility (magenta).
Results for weather and mobility timeseries refer to changes in $\alpha$ for a
unit change of one SD in the input. NPIs targeting large gatherings,
temperature and humidity show the strongest transmission rate reductions
whereas cloudiness leads to the strongest increase. Error bars denote the CI.

##  4 Discussion

In this work we aimed to understand the extent to which fine-scaled spatio-
temporal variations in the spread of SARS-CoV-2 can be explained by
differences in meteorological factors, NPIs, or mobility. We assume that
deviations in time- and age-dependent transmission rates of a district from
the transmission rate in other districts of the same federal state result from
different temperatures, cloudiness, humidity, precipitation, wind, mobility,
or measures targeting schools, gastronomy, healthcare or mass events. We found
that taken together these factors account for more than 60% of the variation
observed in regional transmission. Humidity, cloudiness and precipitation
turned out to be the dominating meteorological factors, and restrictions
targeting mass events showed the greatest reduction in local transmission.

As our approach aims to minimize the quadratic distance between data and model
incidences, the model is optimized to explain the data in pandemic phases of
high incidence. Fig. 1 shows this overall tendency in that the augmented model
fits the first peak (Nov 2020) substantially better than the second peak
(April 2021). Austria experienced its highest incidences within the study
period in most districts in November 2020, when the Austrian 7-day-incidence
per 100,000 reached values of above 500. Our effect estimates are therefore
particularly valid to describe the factors that might explain why certain
districts were more (e.g., the district Rohrbach with a maximal 7-day-
incidence of more than 1,500) or less (e.g., Gänserndorf where the incidence
peaked at about 240) affected in the growth phase of this wave, i.e. the early
growth behaviour of the 2020 seasonal SARS-CoV-2 wave. Note that here we do
not address factors that impact the spread of SARS-CoV-2 in all districts of a
federal state homogeneously, such as NPIs that where implemented on national
scales, seasonal influences (as opposed to fine-scaled weather influences), or
dominant virus variants.

Previous literature reported inconsistent associations of weather with
transmission dynamics. [12, 13, 14, 15, 17, 16] For instance, analysis of an
early outbreak in China found a positive association of the infection rate
with temperature, [12] as did a correlation study in Norway. [13] whereas an
analysis for Spain [17] and a meta-analysis of 202 locations in 8 countries
[15] observed no significant correlation. On the other hand, for the US [16]
and another early, multi-city study in China [35] negative associations were
reported. Also regarding the other meteorological factors the literature
reports contrary results. The previously mentioned analysis of an early
outbreak in China [12] reports that increase in humidity and precipitation
enhances the number of confirmed cases whereas the multi-city study in China
[35] reports that low humidity likely favors the transmission of Covid-19.
Additionally, the correlation study in Norway [13] stated that precipitation
is negatively related with Covid-19 cases. Wind may be a crucial factor in the
spread of infectious diseases. [36] Several studies had examined the relation
of wind speed and the transmission of coronavirus and showed heterogeneous
effects in different countries. [13, 37] Nevertheless, an observational study
using data from 190 countries suggested an inverse association of wind speed
with transmission. [38]

The reasons for this divergence in the literature regarding meteorological
influences are not entirely clear. The studies mentioned above differ greatly
in the included covariates, methodological approaches as well as geographic
and temporal scales on which the analysis was performed, which all hinders
comparability. A nonlinear dependence on temperature has also been proposed
for the transmission rate. [39] In our study that adjusts for larger-scale
regional trends, regional response measures and mobility we observe that
increases in temperature and humidity both and independently lead to decreased
transmission rates. These findings are in line with results from studies
showing that higher temperature and humidity both lead to faster inactivation
of SARS-CoV-2 on surfaces and in aerosols. [40, 41, 42] Airborne SARS-CoV-2 is
also known to be rapidly inactivated by sunlight [45, 43, 44], in line with
our finding that transmission rates increase with cloudiness.

Meteorological factors might impact SARS-CoV-2 transmission also through
behavioural effects. For example, it is feasible that people prefer to meet
inside rather than outside during rainy days which leads to an increase in
COVID-19 cases.

A recent meta-analysis found that school closing was the most effective NPI in
reducing the spread of SARS-CoV-2, followed by workplace closing, business
closing and public event bans. [46] The school measures evaluated here were
substantially less disruptive than full closures and included “soft”
restrictions such as bans in indoor singing and sport activities, the
requirement to wear a cloth mask when entering or leaving the school building
(but not when seated in the class room), or measures to mitigate contacts
between pupils from different classes. We find that such less disruptive
measures nevertheless coincided with significantly reduced transmission rates
in the age group $<20y$ by about $-7.9$% ($-10.6$, $-5.1$). Our finding of a
non-significant effect of school measures on transmission rates in age groups
above $20y$ should not be interpreted to mean transmissions in school settings
are decoupled from the population-level spread. Rather, transmission rates in
older population groups are indirectly impacted by school measures in our
model through the age-structured social mixing matrices.

Regarding gastronomy, we again note that our analysis focuses on less
disruptive measures that did not consist of full closures, but rather of
restrictions such as mandatory registration of visitors, limits for the
opening hours or for the number of people seated at a table. While the
effectiveness of fully closing such venues has been repeatedly established in
the literature, [46] it is maybe surprising to note that less disruptive
interventions also coincide with noticeably reduced transmission rates by
about $-18$% ($-17$, $-19$).

Measures targeting the healthcare sector such as visitor bans and mandatory
FFP2 masks showed a stronger effect than the less disruptive measures
targeting schools and gastronomy, with reductions of$-20.6$% ($-21.2$,
$-20.1$). This result emphasizes once more the necessity to protect hospitals,
[47] long-term care facilities [48] and other health and social institutions
[49] as one of the first lines of defense against SARS-CoV-2.

The NPI with the strongest effect size concerns public events with a reduction
of $-37.5$ ($-38.6$, $-36.5$). In our dataset, this NPI also includes event
bans, particularly on large unseated indoor events. While public event bans
have consistently been identified as one of most effective NPIs, [46] our
effect size exceeds previous estimates which are maximally in the range of
$-25$%. [4] This difference in effect sizes can be explained by the fact that
our effect estimates are susceptible to factors influencing the rapid rise in
case numbers observed at the onset of the seasonal wave in autumn 2020. It is
feasible that the regional onset and early growth behaviour of this seasonal
wave was heavily driven by super-spreading events that have a disproportionate
impact on regional transmission rates while case numbers are still relatively
low. Regions with event bans already in place were not susceptible to such
large increases in transmission rates in the early growth phase, which might
explain our rather large effect estimate.

We find a comparably small but significant transmission-increasing effect of
mobility by about $7.7$% (95% CI $7.4$, $8.0$), even controlling for regional
measures. While it is well-established that mobility serves as a surrogate
measure to quantify the effectiveness of the corresponding NPI regime, [50]
our results indicate that telecommunications data derived mobility estimates
might capture additional behavioural differences.

Our work is subject to several limitations. As our parsimonious
epidemiological modelling approach was designed to be easy and robust to
calibrate, it does not account for incubation periods, distinguish between
symptomatic and asymptomatic infections or contain undetected or quarantined
cases. More concrete, our estimates for the recovery time are much larger than
the reported duration of infectiousness. However, our recovery time should be
understood as giving a time span within which changes in NPIs or weather
events might influence transmission rates. A recent analysis reported that
social distancing measures show a delayed impact of up to 18d to take full
effect. [51] Hence our estimates for $\beta$ effectively combine the
incubation time and duration of infectiousness with a delayed onset of NPI
effectiveness. We also experimented with models that introduce an additional
variable for the delay between input time series and their effect on
transmission rates. This additional parameter did not result in substantial
changes in model quality or effect sizes so we removed it per Occam’s razor.

In terms of limitations we further note that our approach cannot detect
spatio-temporal variations on finer scales than districts or days (e.g.,
whether rain is spread out over multiple hours or in a short and intense
burst) and that we do not address nonlinear dependencies or interaction
effects between individual input time series. There is also a number of NPIs
that was implemented much less than 50 times (e.g., regional quarantine
measures) for which we could not estimate effect sizes in a statistically
robust way.

In conclusion, we find that regional differences in SARS-CoV-2 spread can in
large parts be explained by a combination of meteorological factors and
regional NPIs. Our approach focuses particularly on factors influencing the
peak of the seasonal autumn 2020 wave in Austria, which appears to be mostly
driven by differences in temperature, cloudiness, humidity, and policies
targeting large public events. These findings have implications for what to
expect for upcoming seasonal SARS-CoV-2 waves. In particular, based on our
results we would expect the next wave to commence in regions where low
mitigation measures targeting large public events combine with shifts toward
unfavourable environmental conditions. If no mitigation measures for public
events are in place, and precipitation, cloudiness and humidity move one SD in
the direction of winter conditions, our model expects transmission rates to be
more than twice as high compared to a region with control measures for public
events and more favourable weather. While it has been previously noted that
epidemic forecasting is like weather forecasting due to uncertainties stemming
from nonlinear dynamics calibrated to noisy data streams, [52] our results
suggest that epidemic forecasting to some extent is weather forecasting, with
the added difficulty of needing to control for human behaviour.

##  5 Acknowledgments

We thank Wolfang Knecht for help with the figures. PK acknowledges financial
support from the Vienna Science and Technology Fund WWTF under MA16-045, from
the Medizinisch-Wissenschaftlichen Fonds des Bürgermeisters der
Bundeshauptstadt Wien, no. CoVid004, and the Austrian Science Promotion Agency
FFG under 857136.

##  6 Author Contributions

KL and PK developed the analytical framework and analyzed the data. MK, JC, SL
and CM contributed to developing the analytical framework and analyzing the
data. FW and CW curated the weather data. KH, LK, SS and FB curated the NPI
data. GH curated the mobility data. MB and NP contributed to analyzing the
weather data. KL, MK and PK wrote the first draft of the manuscript. All
authors reviewed and edited the manuscript. PK conceived and designed the
study.

##  7 References

## References

  * [1] Haug, N., Geyrhofer, L., Londei, A. et al. Ranking the effectiveness of worldwide COVID-19 government interventions. Nat Hum Behav 4, 1303–1312 (2020).
  * [2] Islam N, Sharp SJ, Chowell G, et al. Physical distancing interventions and incidence of coronavirus disease 2019: natural experiment in 149 countries. BMJ 370 (2020).
  * [3] Brauner, JM, Mindermann S, Sharma M, et al. Inferring the effectiveness of government interventions against COVID-19. Science 371, 6531 (2021).
  * [4] Li Y, Campbelt H, Kulkarni D, et al. The temporal association of introducing and lifting non-pharmaceutical interventions with the time-varying reproduction number (R) of SARS-CoV-2: a modelling study across 131 countries. The Lance Infectious Diseases 21(2), 193 (2021).
  * [5] Plümper T, Neumayer E. The pandemic predominantly hits poor neighbourhoods? SARS-CoV-2 infections and COVID-19 fatalities in German districts. European Journal of Public Health 30(6) 1176 (2020).
  * [6] Buja A, Paganini M, Cocchio S, et al. Demographic and socio-economic factors, and healthcare resource indicators associated with the rapid spread of COVID-19 in Northern Italy: An ecological study. PLoS One, 15(12), e0244535 (2020).
  * [7] Ehlert, A. The socio-economic determinants of COVID-19: a spatial analysis of German county level data. Socio-Economic Planning Sciences, 101083 (2021).
  * [8] Campbell, F, Archer, B, Laurenson-Schafer, H, et al. Increased transmissibility and global spread of SARS-CoV-2 variants of concern as at June 2021. Eurosurveillance, 26(24), 2100509 (2021).
  * [9] Gómez, CE, Perdiguero, B, Esteban, M. Emerging SARS-CoV-2 variants and impact in global vaccination programs against SARS-CoV-2/COVID-19. Vaccines, 9(3), 243 (2021).
  * [10] Moreno, GK, Braun, KM, Riemersma, KK, et al. Revealing fine-scale spatiotemporal differences in SARS-CoV-2 introduction and spread. Nature Communications, 11(1), 1-13 (2020).
  * [11] Gross, B, Zheng, Z, Liu, S, et al. Spatio-temporal propagation of COVID-19 pandemics. EPL (Europhysics Letters), 131(5), 58003 (2020).
  * [12] Chen Y, Li Q, Karimian H, Chen X, Li X. Spatio-temporal distribution characteristics and influencing factors of COVID-19 in China. Sci Rep. 2021 Feb 12;11(1):3717.
  * [13] Menebo MM. Temperature and precipitation associate with Covid-19 new daily cases: A correlation study between weather and Covid-19 pandemic in Oslo, Norway. Sci Total Environ. 2020 Oct 1;737:139659.
  * [14] Ujiie M, Tsuzuki S, Ohmagari N. Effect of temperature on the infectivity of COVID-19. Int J Infect Dis. 2020 Jun;95:301-303.
  * [15] Pan J, Yao Y, Liu Z, Meng X, Ji JS, Qiu Y, Wang W, Zhang L, Wang W, Kan H. Warmer weather unlikely to reduce the COVID-19 transmission: An ecological study in 202 locations in 8 countries. Sci Total Environ. 2021 Jan 20;753:142272.
  * [16] Smith TP, Flaxman S, Gallinat AS, Kinosian SP, Stemkovski M, et al. Temperature and population density influence SARS-CoV-2 transmission in the absence of nonpharmaceutical interventions. Proceedings of the National Academy of Sciences USA. 2021 Jun 22;118:e2019284118.
  * [17] Briz-Redón Á, Serrano-Aroca Á. A spatio-temporal analysis for exploring the effect of temperature on COVID-19 early evolution in Spain. Sci Total Environ. 2020 Aug 1;728:138811.
  * [18] Badr HS, Du H, Marshall M, Dong E, Squire MM, Gardner LM. Association between mobility patterns and COVID-19 transmission in the USA: a mathematical modelling study. The Lancet Infectious Diseases. 2020 Nov 1;20:1247-1254.
  * [19] Chang S, Pierson E, Koh PW, Gerardin J, Redbird B, Grusky D, Leskovec J. Mobility network models of COVID-19 explain inequities and inform reopening. Nature, 2020 Nov 10;589:82-87.
  * [20] Nouvellet P, Bhatia S, Cori A, Ainslie KEC, Baguelin M, et al. Reduction in mobility and COVID-19 transmission. Nature Communications. 2021 Feb 17;12:1090.
  * [21] Yechezkel M, Weiss A, Rejwan I, Shahmoon E, Ben-Gal S, Yamin D. Human mobility and poverty as key drivers of COVID-19 transmission and control. BMC Public Healh. 2021 Mar 25;21:596.
  * [22] Althouse, BM, Wenger, EA, Miller, JC, et al. Superspreading events in the transmission dynamics of SARS-CoV-2: Opportunities for interventions and control. PLoS biology, 18(11), e3000897 (2020).
  * [23] Popa, A, Genger, JW, Nicholson, MD, et al. Genomic epidemiology of superspreading events in Austria reveals mutational dynamics and transmission properties of SARS-CoV-2. Science translational medicine, 12(573) (2020).
  * [24] Lemieux, JE, Siddle, KJ, Shaw, BM, et al. Phylogenetic analysis of SARS-CoV-2 in Boston highlights the impact of superspreading events. Science, 371(6529) (2021).
  * [25] Liu, Y, Eggo, RM, Kucharski, AJ. Secondary attack rate and superspreading events for SARS-CoV-2. The Lancet, 395(10227), e47 (2020).
  * [26] Seity, Y., Brousseau, P., Malardel, S., Hello, G., Bénard, P., Bouttier, F., Lac, C., Masson, V. The AROME-France Convective-Scale Operational Model, Monthly Weather Review, 139(3), 976-991 (2011).
  * [27] Termonia, P., Fischer, C., Bazile, E., Bouyssel, F., Brožková, R., Bénard, P., Bochenek, B., Degrauwe, D., Derková, M., El Khatib, R., Hamdi, R., Mašek, J., Pottier, P., Pristov, N., Seity, Y., Smolíková, P., Španiel, O., Tudor, M., Wang, Y., Wittmann, C., and Joly, A.: The ALADIN System and its canonical model configurations AROME CY41T1 and ALARO CY40T1, Geosci. Model Dev., 11, 257–281, https://doi.org/10.5194/gmd-11-257-2018 (2018).
  * [28] Petherick, A, Goldszmidt, RG, Andrade, EB, et al. A worldwide assessment of COVID-19 pandemic-policy fatigue. Nature Human Behaviour (2021).
  * [29] Heiler G, Reisch T, Hurt J, Forgahani J, Omani A, et al. Country-wide mobility changes observed using mobile phone data during COVID-19 pandemic. In 2020 IEEE International Conference on Big Data (Big Data) (pp. 3123-3132). IEEE.
  * [30] GOEG/AGES, Manual Corona-Kommission, https://corona-ampel.gv.at/sites/corona-ampel.gv.at/files/Manual%20Ampelsystem_3.0_bf.pdf, accessed: 2021-08-04.
  * [31] https://www.ages.at/en/wissen-aktuell/publikationen/epidemiologische-parameter-des-covid19-ausbruchs-oesterreich-2020/, accessed: 2021-08-04.
  * [32] Brousseau, P., Berre, L., Bouttier, F. and Desroziers, G. Background-error covariances for a convective-scale data-assimilation system: AROME–France 3D-Var. Q.J.R. Meteorol. Soc., 137: 409-422. https://doi.org/10.1002/qj.750 (2011)
  * [33] Anderson, RM, May, RM. Infectious diseases of humans: dynamics and control. (Oxford University Press, England, 1992).
  * [34] Prem, K., Cook, A., and Jit, M.. Projecting social contact matrices in 152 countries using contact surveys and demographic data. PLOS Computational Biology, 13(9):e1005697, 10.1371/journal.pcbi.1005697 (2017)
  * [35] Liu, J., Zhou, J., Yao, J., Zhang, X., Li, L., Xu, X., He, X., Wang, B., Fu, S., Niu. T., Yan, J., Shi, Y., Ren, X., Niu, J., Zhu, W., Li, S., Luo, B., Zhang, K. Impact of meteorological factors on the COVID-19 transmission: A multi-city study in China. Science of The Total Environment (2020), 726, 138513, https://doi.org/10.1016/j.scitotenv.2020.138513
  * [36] Ellwanger JH, Chies JAB. Wind: a neglected factor in the spread of infectious diseases. Lancet Planet Health. 2018 Nov;2(11).
  * [37] Pani SK, Lin NH, RavindraBabu S. Association of COVID-19 pandemic with meteorological parameters over Singapore. Sci Total Environ. 2020;740:140112.
  * [38] Guo C, Bo Y, Lin C, et al. Meteorological factors and COVID-19 incidence in 190 countries: An observational study. Sci Total Environ. 2021;757:143783.
  * [39] Gao, M, Zhou, Q, Yang, X, et al. Nonlinear modulation of COVID-19 transmission by climate conditions. Meteorological Applications, 28(2), e1985 (2021).
  * [40] Biryukov, J, Boydston, JA, Dunning, RA, et al. Increasing temperature and relative humidity accelerates inactivation of SARS-CoV-2 on surfaces. mSphere, 5(4), e00441-20 (2020).
  * [41] Dabisch, P, Schuit, M, Herzog, A, et al. The influence of temperature, humidity, and simulated sunlight on the infectivity of SARS-CoV-2 in aerosols. Aerosol Science and Technology, 55(2), 142-153 (2021).
  * [42] Kwon, T, Gaudreault, NN, Richt, JA. Environmental stability of SARS-CoV-2 on different types of surfaces under indoor and seasonal climate conditions. Pathogens, 10(2), 227 (2021).
  * [43] Ratnesar-Shumate, S, Williams, G, Green, B, et al. Simulated sunlight rapidly inactivates SARS-CoV-2 on surfaces. The Journal of infectious diseases, 222(2), 214-222 (2020).
  * [44] Schuit, M, Ratnesar-Shumate, S, Yolitz, J, et al. Airborne SARS-CoV-2 is rapidly inactivated by simulated sunlight. The Journal of infectious diseases, 222(4), 564-571 (2020).
  * [45] Gunthe, S.S., Swain, B., Patra, S.S. et al. On the global trends and spread of the COVID-19 outbreak: preliminary assessment of the potential relation between location-specific temperature and UV index. J Public Health (Berl.) (2020). https://doi.org/10.1007/s10389-020-01279-y
  * [46] Mendez-Brito, A, El Bcheraoui, C, Pozo-Martin, F. Systematic review of empirical studies comparing the effectiveness of non-pharmaceutical interventions against COVID-19. Journal of Infection, https://doi.org/10.1016/j.jinf.2021.06.018 (2021).
  * [47] Xiang, B, Li, P, Yang, X, et al. The impact of novel coronavirus SARS-CoV-2 among healthcare workers in hospitals: An aerial overview. American journal of infection control, 48(8), 915-917 (2020).
  * [48] Krutikov, M, Hayward, A, Shallcross, L. Spread of a variant SARS-CoV-2 in long-term care facilities in England. New England Journal of Medicine, 384(17), 1671-1673 (2021).
  * [49] Lee, EC, Wada, NI, Grabowski, MK, et al. The engines of SARS-CoV-2 spread. Science, 370(6515), 406-407 (2020).
  * [50] García-Cremades, S, Morales-García, J, Hernández-Sanjaime, R, et al. Improving prediction of COVID-19 evolution by fusing epidemiological and mobility data. Scientific Reports, 11(1), 1-16 (2021).
  * [51] Nader, IW, Zeilinger, EL, Jomar, D, Zauchner, C. Onset of effects of non-pharmaceutical interventions on COVID-19 infection rates in 176 countries. BMC Public Health, 21(1), 1-7 (2021).
  * [52] Moran, KR, Fairchild, G, Generous, N, et al. Epidemic forecasting is messier than weather forecasting: the role of human behavior and internet data streams in epidemic forecast. The Journal of infectious diseases, 214, 404-408 (2016).


# The effect of sublattice symmetry breaking on the electronic properties of a
doped graphene

Alireza Qaiumzadeh  School of Physics, Institute for Research in Fundamental
Sciences (IPM), Tehran 19395-5531, Iran  Institute for Advanced Studies in
Basic Sciences (IASBS), Zanjan, 45195-1159, Iran  Reza Asgari 111Corresponding
author: Tel: +98 21 22280692; fax: +98 21 22280415.
E-mail address<EMAIL_ADDRESS> School of Physics, Institute for
Research in Fundamental Sciences (IPM), Tehran 19395-5531, Iran

###### Abstract

Motivated by a number of recent experimental studies, we have carried out the
microscopic calculation of the quasiparticle self-energy and spectral function
in a doped graphene when a symmetry breaking of the sublattices is occurred.
Our systematic study is based on the many-body G0W approach that is
established on the random phase approximation and on graphene’s massive Dirac
equation continuum model. We report extensive calculations of both the real
and imaginary parts of the quasiparticle self-energy in the presence of a gap
opening. We also present results for spectral function, renormalized Fermi
velocity and band gap renormalization of massive Dirac Fermions over a broad
range of electron densities. We further show that the mass generating in
graphene washes out the plasmaron peak in spectral weight.

###### pacs:

71.10.Ay, 73.63.-b, 72.10.-d, 71.55.-i

##  I INTRODUCTION

Graphene is a single atomic layer of crystalline carbon on the honeycomb
lattice consists of two interpenetrating triangular sublattices A and B, has
opened up a new field for fundamental studies and applications. novoselov1 ;
novoselov2 ; novoselov3 ; novoselov4  Peculiar electronic properties of
graphene give rise possibility to come over silicon-based electronics
limitations. carbon  The single-particle energy spectrum in graphene contains
two zero-energy at $K^{+}$ and $K^{-}$ points of the Brillouin zone which are
called as valleys or Dirac points. Due to the presence of the two carbon atoms
per unit cell, the quasiparticle (QP) need to be described by a two component
wave function.

The charge carriers in a pristine graphene show linear and isotropic energy
dispersion relation and massless chiral behavior for the energy scales up to 1
eV. Recently, graphene has revealed a variety of unusual transport phenomena
characteristics of two-dimensional (2D) Dirac Fermions such as an anomalous
integer quantum Hall effect at room temperature, a minimum quantum
conductivity, Klein tunneling paradox, weak and anti-localization, an absence
of Wigner crystallization phase and Shubnikov-deHaas oscillations that exhibit
a phase shift of $\pi$ due to Berry’s phase. graphene_rev1 ; graphene_rev2 ;
graphene_rev3 ; graphene_rev4 ; graphene_rev5 ; graphene_rev6 ; tomadin  One
important difference between conventional electron gas and Dirac Fermion
particle is that the contribution of exchange and correlation to the chemical
potential is an increasing rather than a decreasing function of carrier-
density. This property implies that exchange and correlation increase the
effectiveness of screening, in contrast to the usual case in which exchange
and correlation weakens screening. This unusual property follows from the
difference in sublattice pseudospin chirality between the Dirac model’s
negative energy valence band states and its conduction band states.

The massless Dirac-like carriers in graphene have almost semi-ballistic
transport behavior with small resistance due to the suppression of back-
scattering process, and moreover graphene is a good thermal conductor. thermal
The mobility of carriers in graphene is quite high morozov1 ; morozov2 ;
morozov3 ; morozov4  which is much higher than the electron mobility revealed
on the semiconductor hetrostructures. eng1 ; eng2  On the other hand, by
measuring the stiffness of materials it is shown that graphene is the
strongest material in two-dimension structures. strong1 ; strong2  These
properties as well as capability to control of the type and density of charge
carriers by gate voltage or the chemical doping dop1 ; dop2 ; dop_lanzara1 ;
dop_lanzara2  make graphene an ideal candidate for superior nano-electronic
devices operating at high frequencies.

Most electronic applications are based on the presence of a gap between the
valence and conduction bands in the conventional semiconductors. The band gap
is a measure of the threshold voltage and on-off ratio of the field effect
transistors (FETs). FET1 ; FET2  Therefore, for integrating graphene into
semiconductor technology, it is crucial to induce a band gap in Dirac points
to control the transport of carriers. Consequently, band gap engineering in
graphene is a hot topic with fundamental and applied significance. gap  In the
literature several routes have being proposed and applied to induce and
control a gap in graphene. One of them is using quantum confined geometries
such as quantum dots and nanoribbons. gap_ribbon1 ; gap_ribbon2 ; gap_ribbon3
; gap_ribbon4 ; gap_ribbon5  It is shown that the gap values increases by
decreasing of nanoribbon width. Another alternative way is spin-orbit coupling
whose origin is due to both intrinsic spin-orbit interactions and the Rashba
interaction. gap_spin1 ; gap_spin2 ; gap_spin3 ; gap_spin4  Another method to
generate a gap in graphene sheets is an inversion symmetry breaking of the
sublattices when the number of electrons on A and B atoms are different
gapsub1 ; gapsub2 ; gapsub3 ; gapsub4  or Kekulé kekule1  distortion, e.g.
graphene on proper substrates dop_lanzara1 ; dop_lanzara2 ; lanzara1 ;
lanzara2 ; eva ; gruneis1 ; gruneis2 ; giovannetti  or adsorb of some
molecules such as water, ammonia ribeiro1 ; ribeiro2  and CrO3 zanella  or an
alkali-metal sub-monolayer on graphene sheets.

Recently angle resolved photoemission spectroscopy ( ARPES) experiments on
graphene epitaxially grown on SiC and ab initio simulations reported a gap
opening in the band structure of graphene placed on proper substrates, and
suggested that interactions between the graphene sheet and the substrate leads
to symmetry breaking of the A and B sublattices and it consequences to induce
a gap in the band structure. Experimenters dop_lanzara1 ; dop_lanzara2 ;
lanzara1 ; lanzara2 ; kruczynski  observed a gap of 260 meV in band structure
of the epitaxial graphene on SiC substrate due to interaction with substrate.
In addition, Zhou et al. dop_lanzara1  found a reversible metal-insulator
transition and a fine tuning of the carriers from electron to hole by
molecular doping in gapped graphene. A Density Functional Theory (DFT)
calculation confirmed a substrate induced symmetry breaking. kim . Their
results showed a gap in the band spectra of graphene about 200 meV which is in
agreement with recent experimental observation. Their calculation determined
that there is a 140 meV on-site energy difference between two sublattices. In
addition, a band gap is observed in spectra of graphene on Ni(111) substrate
gruneis1 ; gruneis2  as well as a gap about of 10 meV in suspended graphene
above a graphite substrate eva  due to sublattice symmetry breaking mechanism.
Moreover, based on the ab initio calculations, it is suggested that boron
nitride substrate induced a gap of 53 meV. giovannetti  Note that the gap
value calculated within DFT is in general underestimating the true band gap
value.

In this paper we consider the sublattice symmetry breaking mechanism for a gap
opening in a pristine doped graphene sheet and study the impact of gap upon
some electronic properties of QPs. To investigate the influence of gap in the
many-body properties of QP in graphene we use the random phase approximation
(RPA) and the G0W approximation. It should be noted that a detailed analysis
provided a framework for the microscopic evaluation of the QP-QP interaction
in the gapless graphene by means of the RPA was carried out by us in Ref.
[im1, ] At the beginning, we review briefly the results of the ground state
thermodynamic properties that we have already presented elsewhere. alireza
Our new results are based on the QP self-energy properties in the presence of
a gap opening in the electronic spectrum. From the self-energy we then obtain
the QP energies, renormalized Fermi velocity, spectral function which can be
compared with ARPES spectra and finally the band gap renormalization of
massive Dirac Fermions in doped graphene. We have shown that mass generating
in graphene washes out a satellite band in the spectral function in agreement
with recent experimental observations. lanzara1

This paper is organized as followed. In Section II we introduce our model
Hamiltonian and then review some ground state properties of gapped graphene.
In Section III we focus on the properties of imaginary and real parts of self-
energy for gapped graphene and then calculate QP spectral function,
renormalized Fermi velocity and band gap renormalization. Finally we conclude
in Section IV.

##  II GROUND STATE THERMODYNAMIC PROPERTIES

We consider the sublattice symmetry breaking mechanism in which the densities
of particles associated to on-site energy $\mu_{a(b)}$, for A(B) sublattice
are different. The electronic structure of graphene can be reasonably good
described using a rather simple tight-binding Hamiltonian, leading to
analytical solutions for their energy dispersion and related eigenstates. The
noninteracting tight binding Hamiltonian for $\pi$ band electrons is
determined by gapsub1 ; gapsub2 ; gapsub3 ; gapsub4

$\displaystyle\hat{H_{0}}$ $\displaystyle=$ $\displaystyle
t\sum_{i}(a_{i}^{\dagger}b_{i}+{\rm
c.c.})+\mu_{a}\sum_{i}a^{\dagger}_{i}a_{i}+\mu_{b}\sum_{i}b^{\dagger}_{i}b_{i}$
(1) $\displaystyle=$ $\displaystyle t\sum_{i}(a_{i}^{\dagger}b_{i}+{\rm
c.c.})+\frac{\mu_{a}-\mu_{b}}{2}\sum_{i}(a^{\dagger}_{i}a_{i}-b^{\dagger}_{i}b_{i})+\frac{\mu_{a}+\mu_{b}}{2}\sum_{i}(a^{\dagger}_{i}a_{i}+b^{\dagger}_{i}b_{i})$

where the sums run over unit cells, $t\simeq 2.7$ eV denotes the nearest
neighbor hopping parameter and $a_{i}(b_{i})$ is Fermi annihilation operator
acts on A(B) sublattice. The second term in the noninteracting Hamiltonian
breaks the inversion symmetry and causes to a band gap with value of
$2\Delta=|\mu_{a}-\mu_{b}|$ at the Dirac points. The last term is a constant
and we left it out. The effective Hamiltonian at low excited energies lead to
a 2D massive Dirac Hamiltonian, $\mathcal{\hat{H}}_{0}=\hbar v_{\rm
F}{\vec{\sigma}}\cdot{\bf k}+\Delta\sigma_{3}$, where $\vec{\sigma}$ are Pauli
matrices and $v_{\rm F}=3ta/2\hbar\simeq 10^{6}$ m/s is the Fermi velocity
where $a\simeq 1.42$ Å is the carbon-carbon distance in honeycomb lattice. The
two eigenvalues of noninteracting Hamiltonian are given by $E_{\bf
k}=\pm\sqrt{(\hbar v_{\rm F}k)^{2}+\Delta^{2}}$ for conduction band (+) and
valance band (-) which is a fully occupied. In addition, the model Hamiltonian
can be used as an approximated model for describing a graphene antidot lattice
in the vicinity of a band gap with a small effective mass value antidot , or
moreover used as an effective Hamiltonian for the intrinsic spin-orbit
interaction in graphene where $\Delta=\Delta_{\rm SO}$ is the strength of the
spin-orbit interaction. gap_spin1 ; gap_spin2 ; gap_spin3 ; gap_spin4  If
$\mu_{a}=\mu_{b}$ the Hamiltonian reduces to massless Dirac Hamiltonian with
two chiral eigenstates having the conical band structures
$\varepsilon_{k}=\pm\hbar v_{\rm F}k$.

We consider the long-range Coulomb electron-electron interaction. We left out
the intervalley scattering and use the two component Dirac Fermion model.
Accordingly, the total interacting Hamiltonian in a continuum model at $K^{+}$
point is expressed as yafis ; Giuliani

$\hat{H}=\sum_{{\bf k},\sigma}\Psi^{\dagger}_{{\bf
k},\sigma}\mathcal{\hat{H}}_{0}\Psi_{{\bf k},\sigma}+\frac{1}{2S}\sum_{{\bf
q}\neq 0}V_{q}({\hat{n}}_{\bf q}{\hat{n}}_{-{\bf q}}-{\hat{N}}),$ (2)

where $\Psi^{\dagger}_{{\bf k},\sigma}=(\psi^{a}_{+,\sigma}({\bf
k}),\psi^{b}_{+,\sigma}({\bf k}))$ is two component pseudospinors of the
noninteracting Hamiltonian, $S$ is the sample area, ${\hat{N}}$ is the total
number operator and $V_{q}=2\pi e^{2}/\epsilon q$ is the bare Coulomb
interaction where $\epsilon$ is an average dielectric constant of the
surrounding medium. The coupling constant in graphene is
$\alpha_{gr}=g_{s}g_{v}e^{2}/\epsilon\hbar v_{\rm F}$ where $g_{s}=g_{v}=2$
being the spin and valley degeneracy, respectively. The coupling constant in
graphene depends only on the substrate dielectric constant while in the
conventional 2D electron systems is density dependent. The typical value of
dimensionless coupling constant is 1 or 2 for graphene supported on a
substrate such a SiC or SiO2.

A central quantity in the many-body techniques is the noninteracting dynamical
polarizability function $\chi^{(0)}({\bf q},i\omega,\mu)$ where $\mu$ is the
chemical potential. The problem of linear density response is set up by
considering a fluid described by the Hamiltonian, $\hat{H}$, which is subject
to an external potential. The external potential must be sufficiently weak for
low-order perturbation theory to suffice. The induced density change has a
linear relation to the external potential through the noninteracting dynamical
polarizability function. This function is recently calculated along the
imaginary frequency axis and it is given by alireza

$\displaystyle\chi^{(0)}({\bf q},i\omega,\mu\geq\Delta)$ $\displaystyle=$
$\displaystyle-\frac{g_{s}g_{v}}{2\pi\hbar^{2}v_{\rm
F}^{2}}\\{\mu-\Delta+\frac{\varepsilon_{q}^{2}}{2}(\frac{\Delta}{{y^{2}}}+\frac{x_{-}^{2}}{2y}\tan^{-1}(\frac{y}{2\Delta}))$
(3) $\displaystyle-$ $\displaystyle\frac{\varepsilon_{q}^{2}}{4y}\Re
e\left[x_{-}^{2}(\sin^{-1}\frac{z(\mu)}{x_{+}}-\sin^{-1}\frac{z(\Delta)}{x_{+}})\right]$
$\displaystyle+$ $\displaystyle\frac{\varepsilon_{q}^{2}}{4y}\Re
e\left[z(\mu)\sqrt{x_{+}^{2}-z^{2}(\mu)}-z(\Delta)\sqrt{x_{+}^{2}-z^{2}(\Delta)}\right]\\},$

where $x_{\pm}=\sqrt{1\pm
4\Delta^{2}/(\varepsilon_{q}^{2}+\hbar^{2}\omega^{2})}$,
$y=\sqrt{\varepsilon_{q}^{2}+\hbar^{2}\omega^{2}}$ and
$z(x)=(2x+i\hbar\omega)/\varepsilon_{q}$. The Fermi energy of a 2D massive
Dirac Fermion system is given by $E_{\rm F}=\mu=\sqrt{(\hbar v_{\rm F}k_{\rm
F})^{2}+\Delta^{2}}$ and the Fermi wavevector is $k_{\rm F}=\sqrt{4\pi
n/g_{s}g_{v}}$ where $n$ is the density of carriers. The noninteracting
density of states (DOS) is determined by
$D(E)=g_{s}g_{v}|E|/2\pi\hbar^{2}v_{\rm F}^{2}\Theta(E^{2}-\Delta^{2})$ which
is density dependent at the Fermi surface. It should be noticed that $D(E_{\rm
F})$ equals to $m/2\pi\hbar^{2}$ in the conventional 2D electron gas system.
Here, $\Theta(x)$ is Heaviside step function.

We now turn to present our first numerical results which are based on the
noninteracting polarization function. The static polarization function as a
function of wavevector for various gap values is shown in Fig. 1(a). The
static polarization function in gapless case is a smooth function whereas a
kink at $q=2k_{\rm F}$ occurs for gapped graphene and thus the derivatives of
$\chi^{(0)}({\bf q},0,\Delta\neq 0)$ has a singular feature. The singular
behavior is the source of several phenomena such as the Friedel oscillations
and moreover the Ruderman-Kittel-Kasuya-Yoshida (RKKY) interaction which the
later is absent in gapless graphene. In Figs. 1(b) and (c) we have plotted the
dynamic polarization function as a function of frequency for wave vectors
smaller and larger than $q=2k_{\rm F}$, respectively. $\chi^{(0)}({\bf
q},i\omega)$ tends to zero like $\omega^{-1}$ at large frequency region.

The polarization function along the real $\omega$ axis can be obtained by
performing analytical continuation of Eq. (3). alireza ; pyatkovskiy1 ;
pyatkovskiy2  In Fig. 2 we have presented the real and imaginary parts of the
noninteracting polarization function as a function of frequency. Sharp cutoffs
in the imaginary part of $\chi^{(0)}({\bf q},\omega)$ are related to the rapid
swing in the real part of $\chi^{(0)}({\bf q},\omega)$. These behaviors are in
result of the fact that the real and imaginary parts of the polarization
function are related through the Kramers-Krönig relations. Importantly, the
sign change of the real part from negative to positive shows a sweep across
the electron-hole continuum. At very large gap values, the polarization
function of massive Dirac Fermions can be reduced to the polarization function
(the Lindhard’s function) of conventional two dimensional electron gas
systems, as they are determined in Figs. 1 and 2. Consequently, we settle
under situation that we can describe a range of band structures from the
Dirac’s cone (gapless graphene) to the parabolic (conventional semiconductors)
band structure behavior by tuning the gap values from zero to a large value,
respectively. We limited our calculations to the intermediate values of of
$\Delta$ and we thus expect wide range of the particular properties related to
unique behavior of the polarization function.

Figure 1: (Color online) (a): The static noninteracting polarization function
as a function of $q$ for various $\Delta$. The dimensionless noninteracting
dynamic polarization at (b): $q=0.5k_{\rm F}$ and (c): $q=2k_{\rm F}$ as a
function of $\omega$ for various $\Delta$.

Figure 2: (Color online) The gap dependence of the real and the imaginary part
of the noninteracting polarization function as a function of $\omega$ for
wavevectors (a), (c) $q=0.5k_{\rm F}$ and (b), (d) $q=2.5k_{\rm F}$

We can calculate the total ground state energy of gapped graphene within RPA.
alireza ; yafis  The ground-state energies can be calculated using the
coupling constant integration technique, which has the contributions
$\varepsilon_{tot}=\varepsilon_{kin}+\varepsilon_{\rm xc}$. The kinetic energy
per particle is given by $\varepsilon_{kin}=2(E_{\rm
F}^{3}-\Delta^{3})/3\varepsilon_{\rm F}^{2}$.

As discussed previously alireza ; yafis  we might subtract the vacuum energy
contribution from the total energy,
$\delta\varepsilon_{tot}=\varepsilon_{tot}(k_{\rm F})-\varepsilon_{tot}(k_{\rm
F}=0)$. Due to the number of states in the Brillouin zone must be conserved,
we do need a ultraviolet cut-off $k_{c}$, which is approximated by $\pi
k_{c}^{2}\simeq(2\pi)^{2}/A_{0}$, where $A_{0}$ is the area of the unit cell.
The dimensionless parameter $\Lambda$ is defined as $k_{c}/k_{\rm
F}\simeq(g_{s}g_{v}n^{-1}\sqrt{3}/9.09)^{1/2}\times 10^{2}$.

Figure 3: (Color online) (a) Exchange-correlation energy and (b)
compressibility as a function of $n^{-1/2}$ ( in units of 10-6 cm) for various
$\Delta$ value.

In Fig. 3, we have shown the exchange- correlation energy in units of
$\varepsilon_{\rm F}=\hbar v_{\rm F}k_{\rm F}$, as a function of $n^{-1/2}$ in
units of 10-6 cm for various $\Delta$ value. The exchange energy arises
entirely from the antisymmetry of the many-body wave function under exchange
of two electrons is positive while the correlation energy, the difference
between the ground state energy and the sum of the kinetic energy and the
exchange energy is negative. This has important implications on the
thermodynamic properties can be calculated from the derivative of the ground
state energy with respect to the density. The compressibility can be
calculated from its definition,
$\kappa^{-1}=n^{2}\partial^{2}(n\delta\varepsilon_{tot})/\partial n^{2}$. Fig.
3(b) shows the ratio between the noninteractiong value,
$\kappa_{0}=2/n\varepsilon_{\rm F}$ and the interaction value of
compressibility as a function of $n^{-1/2}$. The exchange tends to reduce the
compressibility while correlations tends to enhance it. At large $\Delta$, a
minimum structure occurs at the inverse of compressibility behavior and we
expect that at very large $\Delta$, it starts at $\kappa_{0}$ and reduces by
increasing $n^{-1/2}$ behaves like the compressibility of the conversional 2D
electron gas.

##  III THE QP SELF-ENERGY AND THE SPECTRAL FUNCTION

The generation of QPs in an electron liquid leads to two effects. First it
induces a decay of a particle losing momentum via inelastic scattering which
is determined by the imaginary part of self-energy and second is the
renormalization of the dispersion relation of the carriers which is described
by the real part of self-energy. $\Re e\Sigma^{\rm ret}({\bf k},\omega)$ is
defined as the difference between the measured carrier energy $\hbar\omega$,
and the energy of free particle, $\xi_{s{\bf k}}=sE_{{\bf k}}-E_{\rm F}$. To
satisfy causality, the real and imaginary parts of self-energy are related by
a Hilbert transformation. In this section, we first derive the imaginary and
the real part of QP self-energies and then calculate some important quantities
such as a renormalized Fermi velocity, a spectral function and a band gap
renormalization in the presence of a band gap value. These quantities are
related to some important physical properties of both theoretical and
practical applications like the band structure of ARPES, the energy
dissipation rate of injected carriers and the width of the QP spectral
function. kaminski1 ; kaminski2

In the G0W approximation, the self-energy of gapped graphene is given by
($\beta=1/(k_{B}T)$) fateme :

$\displaystyle\Sigma_{s}({\bf{k}},i\omega_{n})$ $\displaystyle=$
$\displaystyle-\frac{1}{\beta}\sum_{s^{\prime}}\int\frac{d^{2}{\bf
q}}{(2\pi)^{2}}F^{ss^{\prime}}({\bf k,k+q})$ $\displaystyle\times$
$\displaystyle\sum_{m=-\infty}^{+\infty}W({\bf
q},i\Omega_{m})G^{(0)}_{s^{\prime}}({\bf k+q},i\omega_{n}+i\Omega_{m}),$

where $W({\bf q},i\Omega_{m})=V_{q}/\epsilon(q,i\Omega_{m})$ is the dynamical
screened effective interaction and
$\epsilon(q,i\Omega_{m})=1-V_{q}\chi^{(0)}(q,i\Omega_{m})$ is dynamical
dielectric function in RPA. The overlap function for gapped graphene
$F^{ss^{\prime}}({\bf k,k+q})$ arises from the graphene band structure is
given by alireza

$\displaystyle F^{ss^{\prime}}({\bf
k,k+q})=\frac{1}{2}(1+ss^{\prime}\frac{\hbar^{2}v_{\rm F}^{2}{\bf k}\cdot({\bf
k}+{\bf q})+\Delta^{2}}{E_{\bf k}E_{\bf k+q}}).$ (5)

It should be noted that $F^{s=-s^{\prime}}({\bf q}=0)=0$. However, in gapless
graphene, intraband backward scattering should not be allowed, namely
$F^{s=s^{\prime}}({\bf q}=-2\textbf{k},\Delta=0)=0$, as well as
$F^{s=-s^{\prime}}({\bf q}=0,\Delta=0)=0$. In Eq. (III), $G_{s}^{(0)}({\bf
k},i\omega)=1/(i\omega-\xi_{s{\bf k}}/\hbar)$ is the noninteracting Green’s
function. Notice that in typical density of carriers in graphene namely
$n>10^{12}$cm-2, the Fermi temperature is about $T_{\rm F}=\varepsilon_{\rm
F}/k_{B}>10^{3}$K, and we therefore can eliminate temperature parameter in our
calculations. To evaluate the zero-temperature retarded self-energy we perform
the line-residue decompositions, $\Sigma^{\rm ret}_{s}({\bf
k},\omega)=\Sigma^{\rm line}_{s}({\bf k},\omega)+\Sigma^{\rm res}_{s}({\bf
k},\omega)$, where $\Sigma^{\rm line}$ is obtained by performing the analytic
continuation before summing over the Matsubara frequencies, and $\Sigma^{\rm
res}$ is the correction which must be taken into account in the total self-
energy. Giuliani  At zero temperature we have

$\displaystyle\Sigma^{\rm line}_{s}({\bf k},\omega)$ $\displaystyle=$
$\displaystyle-\sum_{s^{\prime}}\int\frac{d^{2}{\bf
q}}{(2\pi)^{2}}V_{q}F^{ss^{\prime}}({\bf k,k+q})$ $\displaystyle\times$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\Omega}{2\pi}\frac{1}{\epsilon({\bf
q},i\Omega)}\frac{1}{\omega+i\Omega-\xi_{s^{\prime}}({\bf k+q})/\hbar},$

and

$\displaystyle\Sigma^{\rm res}_{s}({\bf k},\omega)$ $\displaystyle=$
$\displaystyle\sum_{s^{\prime}}\int\frac{d^{2}{\bf
q}}{(2\pi)^{2}}\frac{V_{q}}{\epsilon({\bf q},\omega-\xi_{s^{\prime}}({\bf
k+q})/\hbar)}F^{ss^{\prime}}({\bf k,k+q})$ (7) $\displaystyle\times$
$\displaystyle[\Theta(\omega-\xi_{s^{\prime}}({\bf
k+q})/\hbar)-\Theta(-\xi_{s^{\prime}}({\bf k+q})/\hbar)].$

The line contribution of the self-energy is purely real. The imaginary part of
the self-energy has two contributions where $\Im m\Sigma_{+}^{\rm
ret}({\bf{k}},\omega)=\Im m\Sigma_{+,\rm intra}^{\rm res}({\bf{k}},\omega)+\Im
m\Sigma_{+,\rm inter}^{\rm res}({\bf{k}},\omega)$, and real part of the self
energy can be decomposed as $\Re e\Sigma_{+}^{\rm
ret}({\bf{k}},\omega)=\Sigma_{+}^{\rm line}({\bf{k}},\omega)+\Re
e\Sigma_{+,\rm inter}^{\rm res}({\bf{k}},\omega)+\Re e\Sigma_{+,\rm
intra}^{\rm res}({\bf{k}},\omega)$.

For $\omega>0$ and fixed ${\bf q}$, the RPA decay process represents
scattering of an electron from momentum ${\bm{k}}$ and energy $\omega$ to
${\bm{k}}+{\bm{q}}$ and $\xi_{s^{\prime}}({\bm{k}}+{\bm{q}})$, with all
energies in Eq. (7) measured from the Fermi energy of doped graphene. Since
the Pauli exclusion principle requires that the final state is unoccupied, it
must lie in the conduction band, i.e. $s^{\prime}=+1$. Furthermore since the
Fermi sea is initially in its ground state, the QP must lower its energy, i.e.
$\xi_{s^{\prime}}<\omega$, electrons decay by going down in energy. For
$\omega<0$, the self-energy expresses the decay of holes inside the Fermi sea,
which scatter to a final state, by exciting the Fermi sea. In this case the
final state must be occupied so both band indices are allowed for
$s^{\prime}$, and energy conservation requires that holes decay by moving up
in energy. Since photoemission measures the properties of holes produced in
the Fermi sea by photo ejection, only $\omega<0$ is relevant for this
experimental probe.

In what follows, we calculate the intraband and interband contributions of
self-energy. We have found the intraband term of residue part of self-energy
as following for various values of the frequencies,

$\displaystyle\Sigma_{+,\rm intra}^{\rm res}({\bf k},\omega>0)=C\int_{\rm
max(0,~{}k_{\rm F}-k,~{}k-\beta)}^{k+\beta}dq\int^{\rm min(\hbar\omega+E_{\rm
F},~{}\alpha_{+})}_{\rm max(E_{\rm F},~{}\alpha_{-})}dyf_{+}(y,q)$
$\displaystyle\Sigma_{+,~{}\rm intra}^{\rm res}({\bf k},\Delta-E_{\rm
F}<\hbar\omega<0)=-C\int_{\rm max(0,~{}k-k_{\rm F},~{}\beta-k)}^{k+k_{\rm
F}}dq\int^{\rm min(E_{\rm F},~{}\alpha_{+})}_{\rm max(0,~{}\hbar\omega+E_{\rm
F},~{}\alpha_{-})}dyf_{+}(y,q)$ $\displaystyle\Sigma_{+,\rm intra}^{\rm
res}({\bf k},\hbar\omega<-(\Delta+E_{\rm F}))=-C\int_{\rm max(0,~{}k-k_{\rm
F})}^{k+k_{\rm F}}dq\int^{\rm min(E_{\rm F},~{}\alpha_{+})}_{\rm
max(0,~{}\hbar\omega+E_{\rm F},~{}\alpha_{-})}dyf_{+}(y,q)~{},$ (8)

where

$f_{\pm}(y,q)=\frac{\pm(y\pm E_{k})^{2}\mp q^{2}}{\epsilon({\bf
q},\omega+E_{\rm F}\mp y)\sqrt{4k^{2}q^{2}-(y^{2}-E_{k}^{2}-q^{2})^{2}}}~{},$

$C=e^{2}/2\pi\epsilon E_{k}$, $\alpha_{\pm}=\sqrt{\hbar^{2}v^{2}_{\rm F}(k\pm
q)^{2}+\Delta^{2}}$ and
$\beta=\sqrt{\hbar^{2}\omega^{2}+\hbar^{2}v^{2}k^{2}_{\rm F}+2\hbar\omega
E_{\rm F}}$. On the other hand, the interband contribution of residue part of
the self energy is determined by

$\Sigma_{+,\rm inter}^{\rm res}({\bf k},\hbar\omega<-(\Delta+E_{\rm
F}))=-C\int_{\rm max(0,~{}k-\beta)}^{k+\beta}dq\int_{\alpha_{-}}^{\rm
min(\alpha_{+},-\hbar\omega-E_{\rm F})}dyf_{-}(y,q)~{},$ (9)

and eventually for the line contribution of self-energy we have

$\displaystyle\Sigma_{+,\rm intra}^{\rm line}({\bf{k}},\omega)$
$\displaystyle=$
$\displaystyle-\frac{e^{2}}{4\pi^{2}\epsilon}\int_{0}^{k_{c}}dq\int_{0}^{2\pi}d\phi
F^{++}({\bf
q,q+k},\Delta)\int_{-\infty}^{+\infty}d\Omega\frac{g_{+}(\phi,\Omega,q)}{\epsilon({\bf
q},i\Omega)}$ $\displaystyle\Sigma_{+,\rm inter}^{\rm line}({\bf{k}},\omega)$
$\displaystyle=$
$\displaystyle-\frac{e^{2}}{4\pi^{2}\epsilon}\int_{0}^{k_{c}}dq\int_{0}^{2\pi}d\phi
F^{+-}({\bf
q,q+k},\Delta)\int_{-\infty}^{+\infty}d\Omega\frac{g_{-}(\phi,\Omega,q)}{\epsilon({\bf
q},i\Omega)}$ (10)

where

$g_{\pm}(\phi,\Omega,q)=\frac{\hbar\omega+E_{\rm F}\mp
E_{\bf{k+q}}}{(\hbar\omega+E_{\rm F}\mp E_{{\bf
k+q}})^{2}+\hbar^{2}\Omega^{2}}$ (11)

and $\phi$ denotes an angle between ${\bf k}$ and ${\bf q}$. Note that the
real part of self-energy is $k_{c}$ dependent.

Now we are in a situation that can calculate some important physical
quantities. The QP lifetime or the single-particle relaxation time $\tau$, is
obtained by setting the frequency to the on-shell energy in imaginary part of
the self-energy, $\tau_{s}^{-1}=\Gamma_{s}({\bf k},{\xi_{s\bf
k}}/\hbar)=\frac{2}{\hbar}|\Im m\Sigma_{s}^{\rm ret}({\bf k},\xi_{s\bf
k}/\hbar)|$ where $\Gamma_{s}({\bf k},{\xi_{s\bf k}}/\hbar)$ is the quantum
level broadening of the momentum eigenstate $|s{\bf k}>$. This quantity is
identical with the Fermi’s golden rule expression for the sum of the
scattering rate of a QP and quasihole at wavevector ${\bf k}$. Giuliani  From
Eqs. 8 and 9, one can conclude that total contribution of the imaginary part
of the retarded self-energy on the energy shell comes from the intraband term,
$\Im m\Sigma_{+}^{\rm ret}({\bf k},\xi_{\bf k}/\hbar)=\Im m\Sigma^{\rm
res}_{\rm intra}({\bf k},\xi_{\bf k}/\hbar)$. fateme  In the case of gapless
graphene, scattering rate is a smooth function because of the absence of both
plasmon emission and interband processes. inelastic1 ; inelastic2  However,
with generating a gap and increasing the amount of it, plasmon emission cause
discontinuities in the scattering time, similar to conventional 2D electron
gas. asgari1 ; asgari2  We have thus two mechanisms for scattering of the QPs.
The excitation of electron-hole pairs which is dominant process at long
wavelength regions and the excitation of plasmon appears in a specific wave
vector. As discussed previously fateme , in clean graphene sheets the
inelastic mean free path reduces by increasing the gap whereas the mean free
path is large enough in the range of the typical gap values 10-130 meV, and
thus transport remains in the semi-ballistic regime.

The many-body interactions in graphene as a function of doping can be observed
by ARPES which plays as a central role to investigate QP properties such as
group velocity and lifetime of carriers on the Fermi surface. ARPES is a
useful complementary tool which capable of measuring the constant energy
surfaces for all partially occupied states and the fully occupied band
structure. The information of band dispersion and the Fermi surface can be
elicited from those data measured in ARPES experiments. The relation of the
Green’s function to the single-particle excitation spectrum in the interacting
fluid is expressed by its spectral function. The spectral function is related
to the retarded self-energy by the following expression Giuliani

$\displaystyle A_{s}({\bf k},\omega)=\frac{\hbar}{\pi}\frac{|\Im
m\Sigma_{s}^{\rm ret}({\bf k},\omega)|}{[\hbar\omega-\xi_{s}({\bf k})-\Re
e\delta\Sigma_{s}^{\rm ret}({\bf k,\omega})]^{2}+[\Im m\Sigma_{s}^{\rm
ret}({\bf k},\omega)]^{2}}$ (12)

where $\delta\Sigma_{s}^{\rm ret}({\bf k},\omega)=\Sigma_{s}^{\rm ret}({\bf
k},\omega)-\Sigma_{s}^{\rm ret}(k_{\rm F},0)$, and then ARPES intensity can be
described by $I({\bf k},\omega)=A({\bf k},\omega)n(\omega)$, where $n(\omega)$
is the Fermi-Dirac distribution. The spectral function is the Lorentzian
function where $\Re e\Sigma$ specifying the location of the peak of the
distribution, and $|\Im m\Sigma|$ is the linewidth. The amplitude of the the
Lorentzian function is proportional to $1/|\Im m\Sigma|$. This quantity is the
distribution of energies $\hbar\omega$, in the system when a QP with momentum
${\bf k}$, is added or removed from that. For the noninteracting system we get
$A^{(0)}({\bf k},\omega)=\delta(\omega-\xi({\bf k})/\hbar)$. The Fermi liquid
theory applies only when the spectral function at the Fermi momentum
$A^{(0)}(k=k_{\rm F},\omega)$, behaves as a delta function, and has a
broadened peak indicating damped QPs at $k\neq k_{\rm F}$.

To progress of the interband single particle excitation and plasmon effects on
the $\Im m\Sigma_{s}^{\rm ret}$, we must study the retarded self-energy on the
off-shell frequency which is $\omega\neq\xi_{s{\bf k}}/\hbar$. im1 ; im2  This
quantity gives the scattering rate of a QP with momentum ${\bf k}$ and kinetic
energy $\hbar\omega+E_{\rm F}$. The scattering rate or the linewidth raising
from electron-electron interactions is anisotropic and varies significantly
via wavevector at a constant energy. The imaginary part of self energy shows
the width of the QP spectral function.

In Fig. 4 we have shown the absolute value of the imaginary part of the self
energy in unit of $\varepsilon_{\rm F}$ for various gap values. It would be
noticed that there is an area of frequency is associated to the gap value,
$2\Delta$ in which no QP could enter in. In this case, there is a gap in the
$\Im m\Sigma$ between $\xi_{-,k=0}$ and $\xi_{+,k=0}$. We see that $\Im
m\Sigma_{+}$ vanishes as $\omega^{2}$ for $\omega$ tends to zero, a universal
properties of normal Fermi liquid. Moreover, at large frequency, $\Im
m\Sigma_{+}$ tend towards to $\omega$ linearly. Except from the Dirac point,
the conduction band $\Im m\Sigma_{+}$ peaks broaden because of the dependence
on scattering angle of $\xi({\bf k}+{\bf q})$. For low energy, only intraband
single particle excitation contributes to $\Im m\Sigma$ up to $E_{\rm F}$ and
then the interband single particle excitation contribution increases sharply
about $E_{\rm F}$. The interband contribution increases with increasing the
gap values.

To evaluate the scattering rate in interband channel, we have shown $\Im
m\Sigma_{inter(intra)}$ as a function of frequency in Fig. 5. The intraband
contribution of the imaginary part of self energy associated to scattering
rate of QP in the intraband contribution increases with increasing the gap
values while the interband contribution reduces, as we physically expected.
Moreover, by increasing of the electrons in the conduction band the interband
scattering rate reduces whereas the intraband scattering contribution
increases. The gap value suppresses the scattering rate at
$\omega=-\varepsilon_{\rm F}$.

In Fig. 6 we have plotted the real part of self-energy in unit of
$\varepsilon_{\rm F}$ as a function of the energy for various gap values.
Notice again that the real part of residue self-energy has a gap which is
associated the feature calculated in the imaginary part of self-energy. The
line part of self energy is a continues curve and then we have a jump near to
the boundary of gap values in the $\Re e\Sigma$ for gapped graphene. A kink
around $E_{\rm F}$ is associated to the interband plasmon contribution and it
is broaden due to the gap value. This feature affects noticeably in the
interacting electron density of states.

Figure 4: (Color online) The absolute value of the imaginary part of retarded
self-energy ( $+$ channel) as a function of (a) $\rm\omega$ and (b) $\rm k$,
for the various energy gaps.

Figure 5: (Color online) (a) The intraband and (b) the interband contributions
of the imaginary part of self-energy ( $+$ channel) as a function of $\omega$
for the various energy gaps at $\rm k=0.25k_{\rm F}$.

Figure 6: (Color online) The real part of retarded self-energy ( $+$ channel)
as a function of (a) $\omega$ and (b) $\rm k$, for the various energy gaps at
$\Lambda=100$. The $\Re e\Sigma_{+}$ are measured from the interaction
contribution of the chemical potential $\Re e\Sigma_{+}(k_{\rm F},\omega=0)$.

Figure 7: (Color online) The QP spectral function for $+$ channel as a
function of $\omega$ (a) and $\rm k$ (b), for the various energy gaps at
$\Lambda=100$.

As discussed before im1 ; im2  in a zero temperature and disorder free gapless
graphene, the peaks of the spectral function correspond to the nearly
solutions of Dyson’s equation in which the quasiparticle excitation energies
are obtained by $E=\xi_{+}+\Re e\delta\Sigma_{+}^{\rm ret}$. The intersection
of $\Re e\Sigma$ and the lines $E-\xi_{+}$ indicates a satellite long
wavelength plasmaron peak related to the electron-plasmon excitation due to
the long-range electron-electron Coulomb interaction and the Dyson equation
with $\Im m\Sigma=0$ corresponds to a QP peak related to the single particle
excitation. Importantly, in the presence of gap values, the plasmaron peak
suppressed. In Figs. 7(a) and (b) we have shown the energy distribution curves
(EDC) and momentum distribution curves (MDC), respectively. In the presence of
gap values, as shown in Fig. 7(a) there is only the single QP peak.

Figure 8: (Color online) The imaginary (a) and real part (b) of the $(-)$
channel of retarded self-energy as a function of $\omega$ for the various
energy gaps at $\Lambda=100$. The self-energy are measured from the
interaction contribution of the chemical potential $\Re e\Sigma_{-}(k_{\rm
F},\omega=0)$. The QP spectral function as a function of $\omega$ (c), for the
various energy gaps at $\Lambda=100$.

The valance band self-energy contributions are shown in Fig. 8. There is an
area of frequencies is associated to the gap value in which no QP could exist
in $\Im m\Sigma_{-}$ exactly the same as the conduction band. The $s=+1$ and
$-1$ peaks in $\Im m\Sigma_{s}$ in Figs. 4 and 8 separate at finite $k$
because of chirality factors which emphasize ${\bf k}$ and ${\bf q}$ in nearly
parallel directions for conduction band and ${\bf k}$ and ${\bf q}$ in nearly
opposite directions for valance band states. Consequently, at finite $\Delta$,
the QP peak of $A_{-}(k,\omega)$ which is broaden shifts toward the left in
the opposite behavior of $A_{+}(k,\omega)$. These feature have significant
effects in the interacting electron density of states.

Figure 9: (Color online) The renormalized velocity (+ channel) as a function
of the energy gap for various densities at $\alpha_{gr}=1$.

It is essential to note that the satellite band which is theoretically
predicted im1  for gapless garphene has not been seen in experiments. There
are several reasons that could wash out this feature. For example, the plasmon
damping, disorder effects, electron interactions with the buffer layer and
importantly the effect of gap at Dirac point.

One of the important information which can be extracted from ARPES spectra is
the renormalized Fermi velocity $v^{*}$. A consequence of the interaction is a
Fermi velocity renormalization from the backflow of the fluid around a moving
particle. The density of states at the Fermi energy is also changed. The QP
energy measured from the chemical potential of interacting system
$\delta\varepsilon_{sk}^{QP}$, can be calculated by solving self consistently
the Dyson equation $\delta\varepsilon_{s{\bf k}}^{QP}=\xi_{s{\bf k}}+\Re
e[\delta\Sigma_{s}^{\rm ret}({\bf{\bf
k}},\omega)]|_{\omega=\delta\varepsilon_{s{\bf k}}^{QP}/\hbar}$.Giuliani  In
the isotropic systems the QP energy, depends on the magnitude of ${\bf k}$.
Expanding $\delta\varepsilon_{sk}^{QP}$ to first order in $k-k_{\rm F}$ we can
write $\delta\varepsilon_{sk}^{QP}\simeq\hbar v^{*}_{s}(k-k_{\rm F})$ which
effectively defines the renormalized velocity as $\hbar
v^{*}_{s}={d\delta\varepsilon_{sk}^{QP}}/{dk}|_{k=k_{\rm F}}$. From the Dyson
equation we can calculated the renormalized Fermi velocity as fateme ;
velocity1 ; velocity2

$\displaystyle\frac{v^{*}_{s}}{v_{s}}=\frac{\varepsilon_{\rm F}E_{\rm
F}^{-1}+(\hbar v_{s})^{-1}\partial_{k}\Re e[\delta\Sigma_{s}^{\rm ret}({\bf
k},\omega)]|_{\omega=0,k=k_{\rm F}}}{1-\hbar^{-1}\partial_{\omega}\Re
e[\delta\Sigma_{s}^{\rm ret}({\bf k},\omega)]|_{\omega=0,k=k_{\rm F}}},$ (13)

where $v_{s}=sv_{\rm F}$. It is found before fateme ; velocity1 ; velocity2 ;
gw1 ; gw2  that electron-electron interaction increases the renormalized Fermi
velocity in gapless graphene sheets which this behavior is in contrast to
conventional 2DES. asgari1 ; asgari

Fig. 9 shows the renormalized Fermi velocity in unit of the bare Fermi
velocity as a function of band gap for various carrier densities. The
renormalized Fermi velocity decreasing with increasing the gap value. $v^{*}$
is density independent after $\Delta=0.8\varepsilon_{\rm F}$ which is in good
agreement with recent experiment observation. dop_lanzara1 ; dop_lanzara2

Finally we calculated a band gap renormalization (BGR). bgr1 ; bgr2 ; bgr3 ;
bgr4  The BGR for conductance band is given by the QP self-energy at the band
edge, namely $\rm BGR=\Re e\Sigma_{+}^{\rm ret}({\bf k}=0,\omega=(\Delta-
E_{\rm F})/\hbar)$. Fig. 10 shows the BGR for the various gap values as a
function of the electron density. The BGR decreases by increasing of the
electron density and in the small energy gap values, it is less density
dependent respect to large energy gap values. In gapless case, we have
obtained a induced band gap or kink due to many-body electron-electron
interactions and it tends to a constant with increasing the electron density.
im1 ; im2 ; gw1 ; gw2  This feature is in agreement with the results obtained
within ab intio DFT calculation. gw1

Figure 10: (Color online) The band gap renormalization as a function of
electron density for various energy gap at $\alpha_{gr}=1$.

##  IV SUMMERY AND CONCLUSION

We have revisited the problem of the microscopic calculation of the QP self-
energy and many-body effective velocity suppression in a gapped graphene when
the conduction band is partially occupied. We have performed a systematic
study is based on the many-body G0W approach that is established upon the
random-phase-approximation and on graphene’s massive Dirac equation continuum
model. We have carried out extensive calculations of both the real and the
imaginary part of the QP self-energy and discussed about the interband and
intraband contributions in the scattering process in the presence of gap
value. We have also presented results for the effective velocity and for the
band gap renormalization over a wide range of coupling strength. Accordingly,
we have critically examined the merits of the gap values in dynamical QP
properties.

Most feature of mass generating in graphene is the washing out of the
plasmaron peak in the spectral weight. Increasing of the gap value makes
density independent behavior of the renormalized Fermi velocity. We have shown
that the band gap renormalization in gapped graphene decreases by increasing
the carrier density at large $\Delta$. This is in contrast with the gapless
case in which many body electron-plasmon interactions induce a very small gap
in band structure. These distinct features of the massive Dirac’s Fermions are
related to mixing of the chiralities and reduce of the interband transitions
in graphene sheets.

Acknowledgment

R. A. would like to thank the Scuola Normale Superiore, Pisa, Italy for its
hospitality during the period when the final stage of this work was carried
out. A. Q. supported by IPM grant.

## References

  * (1) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, A. A. Firsov 2004 Science 306 666 .
  * (2) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos and A. A. Firsov 2005 Nature 438 197 .
  * (3) K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V. Khotkevich, S. V. Morozov, and A. K. Geim 2005 Proc. Nat. Acad. Sci. 102 10451 .
  * (4) Y. Zhang, Y. Tan, H. L. Stormer, P. Kim 2005 Nature 438 201 .
  * (5) P. Avouris, Z. Chen, and V. Perebeinos, 2007 Nature Nanotech. 2 605 .
  * (6) A. H. Castro Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim 2009 Rev. Mod. Phys. 81 109 .
  * (7) C. W. Beenakker 2008 Rev. Mod. Phys. 80 1337 .
  * (8) A. K. Geim and P. Kim, 2008 Sci. Am. 298 90 .
  * (9) A. K. Geim and A. H. MacDonald 2007 Physics Today 60 35 .
  * (10) A. K. Geim and K. S. Novoselov 2007 Nature Mater. 6 183 .
  * (11) M. I. Katsnelson, 2007 Materials Today 10 20 .
  * (12) M. Polini, A. Tomadin, R. Asgari and A. H. MacDonald 2008 Phys. Rev. B 78 115426 .
  * (13) A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao, C. N. Lau 2008 Nano Lett. 8 902 .
  * (14) S.V. Morozov, K.S. Novoselov, M.I. Katsnelson, F. Schedin, D.C. Elias, J.A. Jaszczak, A.K. Geim 2008 Phys. Rev. Lett. 100 016602 .
  * (15) K. I. Bolotin, K. J. Sikes, J. Hone, H. L. Stormer, and P. Kim 2008 Phys, Rev, Lett. 101 096802 .
  * (16) Xu Du, Ivan Skachko, Anthoy Barker and Eva Y. Andrei 2008 Nature Nanotech. 3 491 .
  * (17) K.I. Bolotin, K.J. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, H.L. Stormer 2008 Solid State Commun. 146 351 .
  * (18) K. Eng, R. N. McFarland, and B. E. Kane 2005 Appl. Phys. Lett. 87 052106 .
  * (19) E. H. Hwang and S. Das Sarma 2007 Phys. Rev. B 75 073301 .
  * (20) C. Lee, X. Wei, J. W. Kysar, J. Hone 2008 Science 321 385 .
  * (21) M. Neek-Amal and R. Asgari 2009 arXiv:0903.5035
  * (22) T. O. Wehling, K. S. Novoselov, S. V. Morozov, E. E. Vdovin, M. I. Katsnelson, A. K. Geim, A. I. Lichtenstein 2008 Nano Lett. 8 173 .
  * (23) I. Gierz, C. Riedl, U. Starke, C. R. Ast, K. Kern 2008 Nano Lett. 8 4603 .
  * (24) S. Y. Zhou, D. A. Siegel, A. V. Fedorov, and A. Lanzara 2008 Phys. Rev. Lett. 101 086402 .
  * (25) D. A. Siegel, S. Y. Zhou, F. El Gabaly, A. V. Fedorov, A. K. Schmid, and A. Lanzara, 2008 Appl. Phys. Lett. 93 243119 .
  * (26) Y. Lin, K. A. Jenkins, A. Valdes-Garcia, J. P. Small, D. B. Farmer, and P. Avouris 2009 Nano Lett. 9 422 .
  * (27) J. Kedzierski, P. Hsu, P. Healey, P. W. Wyatt, C. L. Keast, M. Sprinkle, C. Berger, and W. A. de Heer 2008 IEEE Trans. Electron Devices 55 2078 .
  * (28) K. Novoselov 2007 Nature Mater. 6 720 .
  * (29) Y. W. Son, M. L. Cohen and S. G. Louie 2006 Phys. Rev. Lett. 97 216803 .
  * (30) M. Y. Han, B. Ozyilmaz, Y. Zhang and P. Kim 2007 Phys. Rev. Lett. 98 206805 .
  * (31) Li Yang, Cheol-Hwan Park, Young-Woo Son, M. L. Cohen, and S. G. Louie 2007 Phys. Rev. Lett. 99 186801 .
  * (32) D. Finkenstadt, G. Pennington, and M. J. Mehl 2007 Phys. Rev. B 76 121405(R) .
  * (33) Y.-W. Son, M. L. Cohen and S. G. Louie 2006 Nature 444 347 .
  * (34) Xue-Feng Wang and T. Chakraborty 2007 Phys. Rev. B 75 033408 .
  * (35) Y. Yao, F. Ye, X. L. Qi, S. C. Zhang, and Z. Fang 2007 Phys. Rev. B 75 041401(R) .
  * (36) C. L. Kane and E. J. Mele 2005 Phys. Rev. Lett. 95 226801 .
  * (37) H. Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, L. Kleinman, and A. H. MacDonald 2006 Phys. Rev. B 74 165310 .
  * (38) G. W. Semenoff 1984 Phys. Rev. Lett. 53 2449 .
  * (39) K. Ziegler, 1996 Phys. Rev. B 53 9653 .
  * (40) V. P. Gusynin, S. G. Sharapov, J. P. Carbotte 2007 Int. J. Mod. Phys. B 21 4611 .
  * (41) A. Bostowick, T. Ohta, J. L. McCesney, K. V. Emtsev, T. Seyller, K. Horn and E. Rotenberg 2007 New J. Phys. 9 385 .
  * (42) C.-Y. Hou, C. Chamon, and C. Mudry 2007 Phys. Rev. Lett. 98 186809 .
  * (43) S. Y. Zhou, G. H. Gweon, A. V. Federov, P. N. First, W. A. de Heer, D. H. Lee, F. Guinea, A. H. Castro Neto, and A. Lanzara 2007 Nature Mater. 6 770 .
  * (44) S.Y. Zhou, D.A. Siegel, A.V. Fedorov, and A. Lanzara 2008 Physica E 40, 2642 .
  * (45) A. Grüneis and D. V. Vyalikh 2008 Phys. Rev. B 77 193401
  * (46) A. Grüneis, K. Kummer and D. V. Vyalikh 2009 New J. Phys. 11 073050 .
  * (47) G. Li, A. Luican, and E. Y. Andrei 2009 Phys. Rev. Lett. 102 176804 .
  * (48) G. Giovannetti, P. A. Khomyako, G. Brocks, P. J. Kelly and J. Van den Brink 2007 Phys. Rev. B 76 073103 .
  * (49) M. Polini, R. Asgari, G. Borghi, Y. Barlas, T. Pereg-Barnea, and A. H. MacDonald 2008 Phys. Rev. B 77 081411(R) .
  * (50) A. Qaiumzadeh and R. Asgari 2009 Phys. Rev. B 79 075414 .
  * (51) R. M. Ribeiro, N. M. R. Peres, J. Coutinho and P. R. Briddon, 2008 Phys. Rev. B 78 075442 .
  * (52) Eduardo V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R. Peres, J. M. B. Lopes dos Santos, Johan Nilsson, F. Guinea, A. K. Geim and A. H. Castro Neto 2007 Phys. Rev. Lett. 99 216802 .
  * (53) I. Zanella, S. Guerini, S. B. Fagan, J. Mendes Filho and A. G. Souza Filho 2008 Phy. Rev. B 77 073404 .
  * (54) M. Mucha-Kruczyński, O. Tsyplyatyev, A. Grishin, E. McCann, Vladimir I. Fal’ko, Aaron Bostwick and Eli Rotenberg 2008 Phys. Rev. B 77 195403 .
  * (55) S. Kim, J. Ihm, H. J. Choi, and Y. Son 2008 Phys. Rev. Lett. 100 176802 .
  * (56) T. G. Pedersen, A. Jauho, and K. Pedersen 2009 Phy. Rev. B 79 113406 .
  * (57) G. F. Giuliani and G. Vignale 2005 Quantum Theory of The Electron Liquid (Cambridge University Press, Cambridge, England) .
  * (58) Y. Barlas, T. Pereg-Barnea, M. Polini, R. Asgari and A. H. MacDonald 2007 Phys. Rev. Lett. 98 236601 .
  * (59) P. K. Pyatkovskiy, 2009 J. Phys.: Condens. Matter  21 025506 .
  * (60) B. Wunsch, T. Stauber, F. Sols and F. Guinea 2006 New J. Phys. 8 318
  * (61) A. Kaminski and H. M. Fretwell 2005 New J. Phys. 7 98 .
  * (62) A. Damascelli, Z. Hussain, and Z.-X. Shen 2003 Rev. Mod. Phys. 75 473 .
  * (63) A. Qaiumzadeh, F. K. Joibari, and R. Asgari 2008 arXv: 0810.4681 .
  * (64) E. H. Hwang, BenYu-Kaung Hu, and S. Das Sarma 2007 Phys. Rev. B 76115434 .
  * (65) W. Tse, E. H. Hwang, and S. Das Sarma 2008 Appl. Phys. Lett. 93 023128 .
  * (66) R. Asgari, B. Davoudi, M. Polini, G. F. Giuliani, M. P. Tosi, and G. Vignale 2005 Phys. Rev. B 71 045323 .
  * (67) G. F. Giuliani and J. J. Quinn 1982 Phys. Rev. B 26 4421 .
  * (68) E. H. Hwang and S. Das Sarma 2008 Phys. Rev. B 77 081412(R) .
  * (69) M. Polini, R. Asgari, Y. Barlas, T. Pereg-Barnea, A. H. MacDonald 2007 Solid State Commun. 143 58 .
  * (70) A. Qaiumzadeh, N. Arabchi, R. Asgari 2008 Solid State Commun. 147 172
  * (71) Paolo E. Trevisanutto, Christine Giorgetti, Lucia Reining, Massimo Ladisa and Valerio Olevano 2008 Phys. Rev. Lett. 101 226405 .
  * (72) C. Park, F. Giustino, C. D. Spataru, M. L. Cohen, and S. G. Louie 2009 Phys. Rev. Lett. 102 076803 .
  * (73) R. Asgari and B. Tanatar 2006 Phys. Rev. B 74 075301 .
  * (74) Y. Zhang and S. Das Sarma 2005 Phys. Rev. B 72 125303 .
  * (75) S. Das Sarma, R. Jalabert, and S. R. Eric Yang 1990 Phys. Rev. B 41 8288 .
  * (76) K. F. Berggren and B. E. Sernelius 1984 Phys. Rev. B 29 5575 .
  * (77) K. F. Berggren and B. E. Sernelius 1981 Phys. Rev. B 24 1971 .