# Dynamical properties of quasiparticles in a gapped graphene sheet

A. Qaiumzadeh  Institute for Advanced Studies in Basic Sciences (IASBS),
Zanjan, 45195-1159, Iran  School of Physics, Institute for research in
fundamental sciences, IPM 19395-5531 Tehran, Iran  F. K. Joibari  School of
Physics, Institute for research in fundamental sciences, IPM 19395-5531
Tehran, Iran  Reza Asgari  School of Physics, Institute for research in
fundamental sciences, IPM 19395-5531 Tehran, Iran

###### Abstract

We present numerical calculations of the impact of charge carriers-carriers
interactions on the dynamical properties of quasiparticles such as
renormalized velocity and quasiparticle inelastic scattering lifetime in a
gapped graphene sheet. Our formalism is based on the many-body
$G_{0}W$-approximation for the self-energy. We present results for the many-
body renormalized velocity suppression and the renormalization constant over a
broad range of energy gap values. We find that the renormalized velocity is
almost independence of the carrier densities at large density regime. We also
show that the quasiparticle inelastic scattering lifetime decreases by
increasing the gap value. Finally, we present results for the mean free path
of charge carriers suppression over the energy gap values.

###### pacs:

71.10.Ay, 81.05.Uw, 71.45.Gm

##  I Introduction

The latest rival to succeed silicon’s status is graphene, a single atomic
layer of graphite make a truly tiny transistor to decrease the size and to
improve the operational speed of the electronic devices. Silicon lost it’s
brilliant electronic properties in pieces smaller than about $10$nm and
practically the smallest silicon chips which has been used in silicon-based
electronics is $45$nm. Furthermore, silicon has some limitations in speed of
operations. These restrictions lead to serious challenges for the Moore’s law
which states that the number of transistors can be placed inexpensively on an
integrated circuit has increased exponentially, doubling approximately every
two years. This growth cannot be maintained forever and thus the search is on
to find and use new materials which may be able to produce higher performance
and better functionality.

The recent discovery novoselov  of graphene in $2004$, and its fabrication
into a field-effect transistor novoselov , has opened up a new field of
physics and offers exciting prospects for new electronic devices and
apparently possible to come over those aforementioned limitations. Graphene
has instructive and unique physics with special intriguing electronic
properties which has attracted remarkable attentions. Barth  First, the
electronic properties of graphene are improved in sizes less than $10$nm .
Second, the massless Dirac-like electrons move through graphene with almost
near-ballistic transport behavior with less resistance because back-scattering
is suppressed. Third, graphene is itself a good thermal conductor such that
graphene’s thermal conductivity is about $\sim 5.3\times 10^{3}$ W/mK at room
temperature which is greater than the thermal conductivity of carbon
nanotubes. thermal  Interestingly, the mobility of carriers in graphene is
quite high and it is about $10^{5}$ cm2/Vs at room temperature. morozov  It is
important to note that the highest electron mobility recorded on the
semiconductor junction H-Si(111)-vacuum FET is $8\times 10^{3}$cm2/Vs at $4.2$
K or the mobility of electrons in junction Si-SiO2(100) MOSFET systems is
$25\times 10^{3}$cm2/Vs at low temperature eng , make graphene promising for
different applications in devices.

Providing capability to control a type and density of charge carriers by gate
voltage or by the chemical doping dop  made graphene instructive for novel
nano-electronic devices. However, a gapped semiconducting behavior would be
more suitable for electronic applications. There have been some proposed in
literature for a gap generation in graphene due to breaking of the sublattice
symmetry by some substrates (such as SiC 3 , graphite 4  and boron nitride5 ),
to adsorbe some molecules (such as water, ammonia 6  and CrO3 cro ), spin-
orbit interaction spin  and finite size effect.size  In case, we are
interested to carry out the microscopic theory to calculate some physical
quantities of gapped graphene.

Theoretical calculations of quasiparticle properties of electron in
conventional two-dimensional electron liquid are performed within the
framework of Landau’s Fermi liquid theory landau  whose key ingredient is the
quasiparticle concept and its interactions. As applied to the electron liquid
model this entails the calculation of effective quasiparticle-quasiparticle
interactions which enter the many-body formalism allowing the calculation of
various physical properties. A number of calculations considered different
variants of the $G_{0}W$-approximation for the self-energy in two-dimensional
electron gasem ; yarlagadda_1994_2 ; em_bohm ; em_dassarma ; zhang ; asgari ;
asgari12  from which density, spin-polarization, and temperature dependence of
quasiparticle properties are obtained.

There is a mechanism for quasiparticle scattering against quasiparticles
because they interact through the Coulomb interaction. This is an inelastic
process and induced a finite lifetime of the quasiparticles. The carrier
lifetime in an epitaxial graphene layers grown on SiC wafers has been recently
measured. tau  Since experiments carried out their measurements on graphene
placed on $\rm SiC$, we expect that graphene was gapped. The experimental
measurements are relevant for understanding carrier intraband and interband
scattering mechanisms in graphene and their impact on electronic and optical
devises. williams ; gu

In this paper we focus on the effect of energy gap on the renormalized
velocity, the inelastic scattering lifetime of quasiparticles and the
inelastic mean free path in gapped graphene sheets over the broad range of
energy gap. Our formalism is based on the Landau-Fermi liquid theory
incorporating the $G_{0}W$-approximation for the self-energy. These quantities
are related to some important physical properties of both theoretical and
practical applications such as the band structure of ARPES spectra im , the
energy dissipation rate of injected carriers tau  and the width of the
quasiparticle spectral function. martin

The contents of the paper are described briefly as follows. In Section II we
discuss about our theoretical model which contains the effect of gap in the
renormalized velocity of quasiparticles and the inelastic scattering lifetime
$\tau_{in}$, of gapped graphene due to electron-electron interactions by using
$G_{0}W$-approximation. Our numerical results are given in Section III.
Finally, Section V contains the summery and conclusions.

##  II Theoretical model

Among the methods designed to deal with the intermediate correlation effects,
of particular interest for its physical appeal and elegance is Landau’s
phenomenological theory landau  dealing with low-lying excitations in a Fermi-
liquid. Landau called such single-particle excitations quasiparticles and
postulated a one-to-one correspondence between them and the excited states of
a non-interacting Fermi gas. He wrote the excitation energy of the Fermi-
liquid in terms of the energies of the quasiparticles and of their effective
interaction. The quasiparticle-quasiparticle interaction function can in turn
be used to obtain various physical properties of the system and can be
parameterized in terms of experimentally measurable data. In this paper, we
will compute the energy gap dependence of the renormalized velocity,
renormalization constant and the inelastic scattering lifetime of
quasiparticle in a gapped graphene sheet.

###  II.1 Quasiparticle renormalized velocity

The dynamics of quasiparticles in a gapless graphene are described by two-
dimensional (2D) massless Dirac Hamiltonian $\hat{H}=\hbar
v{\bf\sigma}\cdot{\bf k}$, with eigenvalues $\varepsilon_{s\bf k}=s\hbar vk$,
where $s=+(-)$ representing right- and left-handed helicity or chirality for
the electrons and holes, respectively. Note that chirality is the same as
helicity for the massless particles. $v=10^{6}$ m$/$s is the Fermi velocity.
As it has been shown before velocity , contrary to conventional 2D electron
systems, the interactions increase the velocity of quasiparticles in graphene
because of interband exchange interactions and the difference between positive
and negative energy branches due to the chirality.

The dynamics of quasiparticles in a gapped graphene are described by 2D
massive Dirac Hamiltonian given by $\hat{H}=\hbar v_{F}\sigma\cdot{\bf
k}+mv^{2}\sigma_{3}$ with eigenvalues $E_{s\bf k}=s\sqrt{(\hbar
vk)^{2}+\Delta^{2}}$ where $\Delta=mv^{2}$ is the gap energy. Due to massive
term in the Hamiltonian, the chirality differs from the helicity and also the
helicity is conserved but is frame dependence.

From the microscopic point of view, the quasiparticle energy can be calculated
by solving the Dyson equation,

$\displaystyle\delta\varepsilon_{s{\bf k}}^{QP}=\xi_{s{\bf k}}+\Re
e[\delta\Sigma_{s}^{ret}({\bf{\bf
k}},\omega)]|_{\omega=\delta\varepsilon_{s{\bf k}}^{QP}/\hbar},$ (1)

where $\xi_{s{\bf k}}=E_{s{\bf k}}-E_{\rm F}$ is the energy of a quasiparticle
relative to the Fermi energy. The Fermi wave vector in graphene is given by
$k_{\rm F}=(4\pi n/g_{s}g_{v})^{1/2}$ where $g_{s}=g_{v}=2$ are spin and
valley degeneracy, respectively. The Fermi energy of gapless graphene is
$\varepsilon_{\rm F}=\hbar vk_{\rm F}$. The retarded self-energy of gapped
graphene is $\Sigma_{s}^{ret}$ and we define $\delta\Sigma_{s}^{ret}({\bf
k},\omega)=\Sigma_{s}^{ret}({\bf k},\omega)-\Sigma_{s}^{ret}(k_{F},0)$. In the
on-shell approximation Giuliani , on the other hand, the above equation must
be solved by setting $\omega=\xi_{s{\bf k}}/\hbar$.

In the $G_{0}W$\- approximation Giuliani , the self-energy of gapped graphene
at finite temperature ($\beta=1/(k_{B}T)$) is given by

$\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 the dynamic screened effective interaction is $W({\bf
q},i\Omega_{m})=V_{q}/\epsilon(q,i\Omega_{m})$ and $\epsilon(q,i\Omega_{m})$
is the dynamical dielectric function and the bare Coulomb interaction is
$V_{q}=2\pi e^{2}/\kappa q$ where $\kappa$ is the averaged background
dielectric constant of graphene is placed on a substrate.
$G_{s}^{(0)}(q,i\Omega_{m})=1/(i\Omega_{m}-\xi_{s{\bf k}}/\hbar)$ is the
standard noninteracting Green’s function. The overlap function for gapped
graphene $F^{ss^{\prime}}({\bf k,k+q})$, is given by alireza

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

To evaluate of the zero temperature retarded self-energy, we decompose the
self-energy into the line which is purely a real function and residue
contributions, $\Sigma^{ret}_{s}({\bf k},\omega)=\Sigma^{line}_{s}({\bf
k},\omega)+\Sigma^{res}_{s}({\bf k},\omega)$, where $\Sigma^{line}$ is
obtained by performing the analytic continuation before summing over the
Matsubara frequencies, and $\Sigma^{res}$ is the correction which must be
taken into account in the total self-energy,

$\displaystyle\Sigma^{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^{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})$ (5) $\displaystyle\times$
$\displaystyle[\Theta(\omega-\xi_{s^{\prime}}({\bf
k+q})/\hbar)-\Theta(-\xi_{s^{\prime}}({\bf k+q})/\hbar)],$

where the dynamic dielectric function is given by $\epsilon({\bf
q},\omega)=1-V_{q}\chi^{(0)}(q,\omega)$ in the random phase approximation
(RPA) and $\chi^{(0)}(q,\omega)$ is the noninteracting polarization function
for gapped graphene. The noninteracting polarization function has been
recently calculated on both along the imaginary and real frequency axis
alireza ; pyatkovskiy . The noninteracting polarization function expressions
along the real frequency axis pyatkovskiy  are given in appendix A.

Note that there are two independent parameters in the self-energy. One of them
is the Fermi energy $E_{\rm F}$, and the other is the dimensionless coupling
constant $\alpha_{gr}=g_{s}g_{v}e^{2}/\kappa\hbar v$. The coupling constant in
graphene depends only on the substrate dielectric constant while in the
conventional 2D electron systems the coupling constant is density dependent.
For graphene placed on SiC or graphite substrates, the coupling constant is
about $\alpha_{gr}\simeq 1$.

The quasiparticle energy depends on the magnitude of ${\bf k}$ for the
isotropic systems. Expanding $\delta\varepsilon_{+k}^{QP}$ to first order in
$k-k_{\rm F}$, we obtain $\delta\varepsilon_{+k}^{QP}\simeq\hbar
v^{*}(k-k_{\rm F})$ which effectively defines the renormalized velocity as
$\hbar v^{*}=d\delta\varepsilon_{+k}^{QP}/dk|_{k=k_{\rm F}}$. The renormalized
velocity in the Dyson scheme is thus given by

$\displaystyle\frac{v^{*}}{v}=\frac{(1+\Delta^{2})^{-1/2}+v^{-1}\partial_{k}\Re
e[\delta\Sigma_{+}^{ret}({\bf k},\omega)]|_{\omega=0,k=k_{\rm
F}}}{1-\partial_{\omega}\Re e[\delta\Sigma_{+}^{ret}({\bf
k},\omega)]|_{\omega=0,k=k_{\rm F}}}~{}.$ (6)

In the on-shell approximation, on the other hand, the renormalized velocity is
given by $v^{*}/v=(1+\Delta^{2})^{-1/2}+v^{-1}\partial_{k}\Re
e[\delta\Sigma_{+}^{ret}({\bf k},\omega)]|_{\omega=0,k=k_{\rm
F}}+(1+\Delta^{2})^{-1/2}\partial_{\omega}\Re e[\delta\Sigma_{+}^{ret}({\bf
k},\omega)]|_{\omega=0,k=k_{\rm F}}$. The renormalized velocity in this
approximation demonstrates qualitatively the same behavior obtained by the
Dyson equation, Eq. (6) but its magnitude is larger than the one calculated
within the Dyson scheme. asgari  There is an ultraviolet divergence in the
wave vector integrals of the line contribution in a continuum model formulated
as discussed above. velocity  We introduce an ultraviolet cutoff for the wave
vector integrals, $k_{c}=\Lambda k_{F}$ which is the order of the inverse
lattice spacing and $\Lambda$ is dimensionless quantity. For definiteness we
take $\Lambda=k_{c}/k_{F}$ to be such that $\pi(\Lambda
k_{F})^{2}=(2\pi)^{2}/{\cal A}_{0}$, where ${\cal A}_{0}=3\sqrt{3}a^{2}_{0}/2$
is the area of the unit cell in the honeycomb lattice, with $a_{0}\simeq 1.42$
Å the carbon-carbon distance. With this choice,
$\Lambda\simeq{(gn^{-1}\sqrt{3}/9.09)^{1/2}}\times 10^{2}$, where $n$ is the
electron density in units of $10^{12}~{}{\rm cm}^{-2}$.

An important quantity in the Fermi-liquid theory is the renormalization
constant $Z$, defined as the square of the overlap between the state of the
system after adding (or removing) of an electron with the Fermi wave vector
and the ground-state of the system. The non-zero renormalization constant
value is always smaller than the one for the normal Fermi-liquid systems and
can be calculated explicitly as followGiuliani

$\displaystyle Z=\frac{1}{1-\partial_{\omega}\Re e[\delta\Sigma_{+}^{ret}({\bf
k},\omega)]|_{\omega=0,k=k_{\rm F}}}.$ (7)

We will show that $Z$ is a finite number for gapped graphene and it confirms
as well that the system is a Fermi-Liquid.

###  II.2 Inelastic scattering lifetime

In this subsection, we compute the inelastic scattering lifetime of
quasiparticles due to carriers-carriers interactions at zero temperature and
disorder-free for gapped graphene sheets. This is obtained through the
imaginary part of the self-energy inelastic  when the frequency evaluated at
the on-shell energy.

$\displaystyle\tau_{in}^{-1}({\bf k})=\Gamma_{in}({\bf k},{\xi_{+\bf
k}}/\hbar)=-\frac{2}{\hbar}\Im m\Sigma_{+}^{ret}({\bf k},\xi_{+\bf k}/\hbar),$
(8)

where $\Gamma_{in}({\bf k},{\xi_{s\bf k}}/\hbar)$ is the quantum level
broadening of the momentum with eigenstate $|s{\bf k}>$. It is worthwhile to
note that the expression of $\tau_{in}^{-1}$ is identical with a result
obtained by the Fermi’s golden rule summing the scattering rate of electron
and hole contributions at wave vector ${\bf k}$. Giuliani  Note again that the
total contribution of the imaginary part of the retarded self-energy comes
from the residue term both intra- and interband contributions, $\Im
m\Sigma^{ret}_{+}({\bf k},\omega)=\Im m\Sigma^{res}_{intra}({\bf
k},\omega)+\Im m\Sigma^{res}_{inter}({\bf k},\omega)$. However, the total
contribution of the imaginary part of the retarded self-energy evaluated at
the on-shell energy comes only from intraband term, $\Im
m\Sigma_{+}^{ret}({\bf k},\xi_{\bf k}/\hbar)=\Im m\Sigma^{res}_{intra}({\bf
k},\xi_{\bf k}/\hbar)$. We will discuss about that with more details in the
appendix B and C.

We turn our attention to investigate the imaginary part of the retarded self-
energy with more details. By starting from Eq. (5), we end up to an expression
for the imaginary part of self-energy which is given by,

$\displaystyle\Im m\Sigma^{ret}_{+}({\bf k},\omega)$ $\displaystyle=$
$\displaystyle\Im m\Sigma^{res}_{intra}({\bf k},\omega)+\Im
m\Sigma^{res}_{inter}({\bf k},\omega)$ (9) $\displaystyle=$
$\displaystyle\int\frac{d^{2}{\bf q}}{(2\pi)^{2}}V_{q}\Im m[\epsilon^{-1}({\bf
q},\omega-\xi_{+}({\bf k+q})/\hbar)]F^{++}({\bf k,k+q})$ $\displaystyle\times$
$\displaystyle[\Theta(\omega-\xi_{+}({\bf k+q})/\hbar)-\Theta(-\xi_{+}({\bf
k+q})/\hbar)]$ $\displaystyle+$ $\displaystyle\int\frac{d^{2}{\bf
q}}{(2\pi)^{2}}V_{q}\Im m[\epsilon^{-1}({\bf q},\omega-\xi_{-}({\bf
k+q})/\hbar)]F^{+-}({\bf k,k+q})$ $\displaystyle\times$
$\displaystyle[\Theta(\omega-\xi_{-}({\bf k+q})/\hbar)-\Theta(-\xi_{-}({\bf
k+q})/\hbar)].$

where the imaginary part of the inverse dielectric function in RPA level is
obtained by

$\Im m[\epsilon^{-1}({\bf q},\omega)]=\frac{V_{q}\Im m\chi^{(0)}({\bf
q},\omega)}{[1-V_{q}\Re e\chi^{(0)}({\bf q},\omega)]^{2}+[V_{q}\Im
m\chi^{(0)}({\bf q},\omega)]^{2}}.$ (10)

It is worth to note that the plasmon contributions in the imaginary part of
self-energy comes from the zero-solutions of denominator in Eq. (10).

##  III Numerical results

We turn to a presentation of our main numerical results. We present some
illustrative results for the quasiparticle dynamic properties such
renormalized velocity, renormalization constant and inelastic scattering
lifetime. All numerical data are calculated in the Dyson scheme at $\alpha=1$.

The Fermi liquid phenomenology of Dirac electrons in gapless graphene im ;
velocity  and conventional 2D electron liquid asgari  have the same structure,
since both systems are isotropic and have a single circular Fermi surface. The
strength of interaction effects in a conventional 2D electron liquid increases
with decreasing carrier density. At low densities, the quasiparticle
renormalization constant $Z$ is small, the renormalized velocity is suppressed
asgari , the charge compressibility changes sign from positive to negative,
and the spin-susceptibility is strongly enhanced asgari2 . These effects
emerge from an interplay between exchange interactions and quantum
fluctuations of charge and spin in the 2D electron liquid.

In the 2D massless electron graphene, on the other hand, it has been shown
yafis ; velocity ; im  that interaction effects also become noticeable with
decreasing density, although more slowly, the quasiparticle renormalization
constant, $Z$ tends to larger values, that the renormalized velocity is
enhanced rather than suppressed, and that the influence of interactions on the
compressibility and the spin-susceptibility changes sign. These qualitative
differences are due to exchange interactions between electrons near the Fermi
surface and electrons in the negative energy sea and to interband
contributions to Dirac electrons from charge and spin fluctuations.

In this paper we have shown the results for gapped graphene which are
determined values between the gapless graphene evaluated at $\triangle=0$ and
the conventional 2D electron liquid where $\triangle\rightarrow\infty$.

In Fig. 1, we have plotted the renormalized velocity as a function of carrier
density for the various energy gap. As a result, we see that the impact of
energy gap on quasiparticles velocity which is similar to the effect of
impurity to that on grapheneimpurity . The renormalized velocity is almost
density independent in gapped graphene at large carrier densities. The
renormalized velocity reduces dramatically by increasing the energy gap
especially in the low carrier densities. Importantly, the renormalized
velocity becomes less than the bare velocity at large energy gap and low
density values. It is physically accepted since the system tends to
conventional 2D electron liquid by increasing the energy gap values. Note that
in the conventional 2D electron systems, the renormalized velocity is
suppressed by increasing the coupling constant or reducing the density.

We have shown the renormalization constant $Z$, as a function of the energy
gap in Fig. 2. The renormalization constant enormously reduces by increasing
the energy gap in mild densities, however it decreases quite slowly in high
densities.

Fig. 3(a) is shown the absolute value of $\Im m\Sigma^{\rm\scriptstyle
ret}_{+}({\bf k},\omega)$ as from Eq. (9), evaluated at $\omega=\xi_{\bf
k}/\hbar$. By increasing the gap value, this function takes a finite jump at
the wave number of the plasmon dip and at large $\triangle$ values, a
discontinuity appears. The discontinuity is peculiar to $2D$ electron liquid.
giuliani_quinn  It is absent in gapless graphene and starts to arise from the
fact that the oscillator strength of the plasmon pole is non-zero at special
$k$ value for gapped graphene.

Fig. 3(b) is clearly shown the behavior of the energy gap dependence of the
inverse inelastic scattering lifetime. As it is argued in the Appendix B, the
imaginary part of self-energy evaluated at the on-shell energy start from
$\triangle-\sqrt{\varepsilon_{\rm F}^{2}+\triangle^{2}}$ and in case the
results are truncated below that. The quasiparticle lifetime decreases by
increasing the gap value and it is a clear difference between 2D massless
Dirac electron and gapped graphene. Consequently, the inelastic scattering
lifetime in graphene is always larger than the conventional 2D electron
liquid. In the case of gapless graphene, scattering rate is a smooth function
because of the absence of both plasmon emission and interband
processes,inelastic  nevertheless with generating a gap and increasing the
amount of it, plasmon emission causes to arise a discontinuity in the
scattering time, similar to conventional 2D electron liquid. asgari  We have
thus two mechanisms for scattering of the quasiparticles. The excitation of
electron-hole pairs which is dominant process at low wave vectors and the
excitation of plasmon appears in a specific wave vector. We also see in Fig.
3(b) that the scattering rate is quite sensitive to the gap energy and the
scattering rate increases by increasing the energy gap.

In Fig. 4, we have depicted the inelastic mean free path $l_{in}({\bf
k})=v^{*}\tau_{in}({\bf k})$, as a function of the on-shell energy for various
gap energies. To this purpose we multiplied the results of $\tau_{in}({\bf
k})$ to a proper renormalized velocity. As a result the mean free path of a
gapped graphene is shorter than that obtained for gapless graphene.
Furthermore, the massless graphene has larger $l_{in}$ and it decreases by
increasing the energy gap values. Note that the typical value of energy gap
due to breaking sublattice symmetry is $\triangle=10-100$ meV corresponding
the inelastic mean free path is $l_{in}=20-50$ nm which implies that the
system remains in the semi-ballistic regime. 3 ; 4 ; 5 ; 6 .

##  IV Summery and concluding remarks

In summary, we have studied the problem of the microscopic calculation of the
quasiparticle self-energy and many-body renormalized velocity suppression over
the energy gap in a gapped graphene. We have carried out calculations of both
the real and the imaginary part of the quasiparticle self-energy within
$G_{0}W$-approximation. We have also presented results for the renormalized
velocity suppression and for the renormalization constant over a wide range of
energy gap. We have shown that the renormalized velocity for a gapped graphene
is almost independent of the carrier density at high density. We have finally
presented results for the quasiparticle inelastic scattering lifetime
suppression over the energy gap and show that the mean free path of the charge
carriers of a gapless graphene is larger than a gapped graphene one. In case
the mean free path of charge carriers decrease by increasing the energy gap.

A possible role of correlations including the charge-density fluctuations
beyond the Random Phase Approximation, remains to be examined.

###### Acknowledgements.

R. A. would like to thank the International Center for Theoretical Physics,
Trieste for its hospitality during the period when part of this work was
carried out. A. Q is supported by IPM grant.

##  Appendix A The dynamic polarization function for a gapped graphene

In this appendix we present the real and imaginary part of the noninteracting
polarization function for a gapped graphene, which is calculated recently by
Pyatkovskiy. pyatkovskiy  The dynamic polarization function for gapped
graphene in the imaginary frequency axis is also calculated by us in Ref.
[alireza, ]. Importantly, the noninteracting polarization function along the
imaginary frequency axis can be obtained by performing analytical continuation
from real axis and those results are the same. private

First, by introducing some following notations,

$\displaystyle f(k,\omega)$ $\displaystyle=$
$\displaystyle\frac{g_{s}g_{v}k^{2}}{16\pi\sqrt{|\hbar^{2}v^{2}k^{2}-\hbar^{2}\omega^{2}|}},$
$\displaystyle g_{\pm}$ $\displaystyle=$ $\displaystyle\frac{2E_{\rm
F}\pm\hbar\omega}{\hbar vk},$ $\displaystyle x_{0}$ $\displaystyle=$
$\displaystyle\sqrt{1+\frac{4\Delta^{2}}{\hbar^{2}v^{2}k^{2}-\hbar^{2}\omega^{2}}},$
$\displaystyle G_{<}(x)$ $\displaystyle=$ $\displaystyle
x\sqrt{x_{0}^{2}-x^{2}}-(2-x_{0}^{2})\cos^{-1}(x/x_{0}),$ $\displaystyle
G_{>}(x)$ $\displaystyle=$ $\displaystyle
x\sqrt{x^{2}-x_{0}^{2}}-(2-x_{0}^{2})\cosh^{-1}(x/x_{0}),$ $\displaystyle
G_{0}(x)$ $\displaystyle=$ $\displaystyle
x\sqrt{x^{2}-x_{0}^{2}}-(2-x_{0}^{2})\sinh^{-1}(x/\sqrt{-x_{0}^{2}}),$ (11)

the real part of noninteracting polarization function is given by,

$\displaystyle\Re e\chi^{(0)}(k,\omega)$ $\displaystyle=$
$\displaystyle-\frac{g_{s}g_{v}E_{\rm F}}{2\pi
v_{F}^{2}}+f(k,\omega)\times\left\\{\begin{array}[]{ll}0,&\hbox{1A}\\\
G_{<}(g_{-}),&\hbox{2A}\\\ G_{<}(g_{+})+G_{<}g_{-}),&\hbox{3A}\\\
G_{<}(g_{-})-G_{<}(g_{+}),&\hbox{4A}\\\
G_{>}(g_{+})-G_{>}(g_{-}),&\hbox{1B}\\\ G_{>}(g_{+}),&\hbox{2B}\\\
G_{>}(g_{+})-G_{>}(-g_{-}),&\hbox{3B}\\\
G_{>}(-g_{-})+G_{>}(g_{+}),&\hbox{4B}\\\
G_{0}(g_{+})-G_{0}(g_{-}),&\hbox{5B}\\\ \end{array}\right.$ (21)

and the imaginary part of noninteracting polarization function is given by,

$\displaystyle\Im m\chi^{(0)}(k,\omega)$ $\displaystyle=$ $\displaystyle
f(k,\omega)\times\left\\{\begin{array}[]{ll}G_{>}(g_{+})-G_{>}(g_{-}),&\hbox{1A}\\\
G_{>}(g_{+}),&\hbox{2A}\\\ 0,&\hbox{3A}\\\ 0,&\hbox{4A}\\\ 0,&\hbox{1B}\\\
-G_{<}(g_{-}),&\hbox{2B}\\\ \pi(2-x_{0}^{2}),&\hbox{3B}\\\
\pi(2-x_{0}^{2}),&\hbox{4B}\\\ 0,&\hbox{5B}\\\ \end{array}\right.$ (31)

with the followings regions in the $(k,\omega)$ space,

$\displaystyle\begin{array}[]{cc}1A&\hbar\omega<E_{\rm
F}-\sqrt{\hbar^{2}v^{2}(k-k_{\rm F})^{2}+\Delta^{2}},\\\ 2A&|E_{\rm
F}-\sqrt{\hbar^{2}v_{F}^{2}(k-k_{\rm F})^{2}+\Delta^{2}}|<\hbar\omega<-E_{\rm
F}+\sqrt{\hbar^{2}v^{2}(k+k_{\rm F})^{2}+\Delta^{2}},\\\
3A&\hbar\omega<-E_{\rm F}+\sqrt{\hbar^{2}v^{2}(k-k_{\rm
F})^{2}+\Delta^{2}},\\\ 4A&-E_{\rm F}+\sqrt{\hbar^{2}v^{2}(k+k_{\rm
F})^{2}+\Delta^{2}}<\hbar\omega<\hbar vk,\\\ 1B&k<2k_{\rm
F},~{}~{}~{}\sqrt{\hbar^{2}v^{2}k^{2}+4\Delta^{2}}<\hbar\omega<E_{\rm
F}+\sqrt{\hbar^{2}v^{2}(k-k_{\rm F})^{2}+\Delta^{2}},\\\ 2B&E_{\rm
F}+\sqrt{\hbar^{2}v^{2}(k-k_{\rm F})^{2}+\Delta^{2}}<\hbar\omega<E_{\rm
F}+\sqrt{\hbar^{2}v^{2}(k+k_{\rm F})^{2}+\Delta^{2}},\\\ 3B&\hbar\omega>E_{\rm
F}+\sqrt{\hbar^{2}v^{2}(k+k_{\rm F})^{2}+\Delta^{2}},\\\ 4B&k>2k_{\rm
F},~{}~{}~{}\sqrt{\hbar^{2}v^{2}k^{2}+4\Delta^{2}}<\hbar\omega<E_{\rm
F}+\sqrt{\hbar^{2}v^{2}(k-k_{\rm F})^{2}+\Delta^{2}},\\\ 5B&\hbar
vk<\hbar\omega<\sqrt{\hbar^{2}v^{2}k^{2}+4\Delta^{2}},\\\ \end{array}$ (41)

##  Appendix B The intraband contribution of self-energy

Since we are interested in quasiparticle properties, we therefore need only
$s=+$ contribution. Let us focus on the intraband contribution of the retarded
self-energy. The second argument of the dielectric function in Eq. (5) ( by
setting $\hbar=v=1$) is

$\displaystyle\omega-\xi_{+}({\bf k+q})=\omega+E_{\rm
F}-\sqrt{k^{2}+q^{2}+2kq\cos\phi+\Delta^{2}}~{}.$ (42)

In this case, we change the variable $\phi$ and integrate it over
$y=\sqrt{k^{2}+q^{2}+2kq\cos\phi+\Delta^{2}}$. Using the new variable, the
intraband contribution of self-energy changes to

$\displaystyle\Sigma_{intra}^{res}({\bf k},\omega)$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{2\pi\kappa\sqrt{k^{2}+\Delta^{2}}}\int_{0}^{+\infty}dq\int_{\sqrt{(k-q)^{2}+\Delta^{2}}}^{\sqrt{(k+q)^{2}+\Delta^{2}}}\frac{dy}{\sqrt{4k^{2}q^{2}-(y^{2}-k^{2}-q^{2}-\Delta^{2})^{2}}}\frac{(\sqrt{k^{2}+\Delta^{2}}+y)^{2}-q^{2}}{\epsilon({\bf
q},\omega+E_{\rm F}-y)}$ (43) $\displaystyle\times$
$\displaystyle[\Theta(\omega+E_{\rm F}-y)-\Theta(E_{\rm F}-y)].$

We can now simplify the $\Theta$-functions further in Eq. (43) by considering
the positive and negative regions of $\omega$ as follow

$\displaystyle 1)$ $\displaystyle\omega$ $\displaystyle+E_{\rm
F}-y>0~{}and~{}~{}E_{\rm F}-y<0:~{}It~{}implies~{}that~{}\omega>0~{},$
$\displaystyle 2)$ $\displaystyle\omega$ $\displaystyle+E_{\rm
F}-y<0~{}and~{}~{}E_{\rm F}-y>0:~{}It~{}implies~{}that~{}\omega<0~{}.$

To consider the first case where $\omega>0$, the difference between the two
$\Theta$-functions in Eq. (43) is equal to $+1$ if

$\displaystyle E_{\rm F}<y<\omega+E_{\rm
F}~{}~{}and~{}~{}\sqrt{(k-q)^{2}+\Delta^{2}}<y<\sqrt{(k+q)^{2}+\Delta^{2}}~{}.$
(44)

Now we do need to find the overlap between these two intervals. We simply end
up to inequivalent conditions which are $q>k-\sqrt{\omega^{2}+k_{\rm
F}^{2}+2\omega\sqrt{k_{\rm F}^{2}+\Delta^{2}}}$, $q<k+\sqrt{\omega^{2}+k_{\rm
F}^{2}+2\omega\sqrt{k_{\rm F}^{2}+\Delta^{2}}}$ and $q>k_{\rm F}-k$.
Collecting everything together and using the fact that $q\geq 0$, we finally
find

$\displaystyle\Sigma_{intra}^{res}({\bf k},\omega>0)$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{2\pi\kappa\sqrt{k^{2}+\Delta^{2}}}\int_{max(0,k_{\rm
F}-k,k-\sqrt{\omega^{2}+k_{\rm F}^{2}+2\omega\sqrt{k_{\rm
F}^{2}+\Delta^{2}}})}^{k+\sqrt{\omega^{2}+k_{\rm F}^{2}+2\omega\sqrt{k_{\rm
F}^{2}+\Delta^{2}}}}dq\int^{min(\omega+\sqrt{k_{\rm
F}^{2}+\Delta^{2}},\sqrt{(k+q)^{2}+\Delta^{2}})}_{max(\sqrt{k_{\rm
F}^{2}+\Delta^{2}},\sqrt{(k-q)^{2}+\Delta^{2}})}dy$ (45) $\displaystyle\times$
$\displaystyle\frac{(\sqrt{k^{2}+\Delta^{2}}+y)^{2}-q^{2}}{\epsilon({\bf
q},\omega+E_{\rm F}-y)\sqrt{4k^{2}q^{2}-(y^{2}-k^{2}-q^{2}-\Delta^{2})^{2}}}$

By considering of the second case where $\omega<0$, the difference between the
two $\Theta$-functions in Eq. (5) is equal to $-1$ if

$\displaystyle E_{\rm F}+\omega<y<E_{\rm
F}~{}~{}and~{}~{}\sqrt{(k-q)^{2}+\Delta^{2}}<y<\sqrt{(k+q)^{2}+\Delta^{2}}$
(46)

As what we did before, we calculate overlap between intervals and thus we find
$q>k-k_{\rm F}$ and $q<k+k_{\rm F}$, $q>\sqrt{\omega^{2}+k_{\rm
F}^{2}+2\omega\sqrt{k_{\rm F}^{2}+\Delta^{2}}}-k$. Putting everything together
and using the fact that $q\geq 0$ we finally find

$\displaystyle\Sigma_{intra}^{res}({\bf k},\Delta-E_{\rm F}<\omega<0)$
$\displaystyle=$
$\displaystyle-\frac{e^{2}}{2\pi\kappa\sqrt{k^{2}+\Delta^{2}}}\int_{max(0,k-k_{\rm
F},\sqrt{\omega^{2}+k_{\rm F}^{2}+2\omega\sqrt{k_{\rm
F}^{2}+\Delta^{2}}}-k)}^{k+k_{\rm F}}dq$ (47) $\displaystyle\times$
$\displaystyle\int^{min(\sqrt{k_{\rm
F}^{2}+\Delta^{2}},\sqrt{(k+q)^{2}+\Delta^{2}})}_{max(0,\omega+\sqrt{k_{\rm
F}^{2}+\Delta^{2}},\sqrt{(k-q)^{2}+\Delta^{2}})}dy$ $\displaystyle\times$
$\displaystyle\frac{(\sqrt{k^{2}+\Delta^{2}}+y)^{2}-q^{2}}{\epsilon({\bf
q},\omega+E_{\rm F}-y)\sqrt{4k^{2}q^{2}-(y^{2}-k^{2}-q^{2}-\Delta^{2})^{2}}}$
$\displaystyle\Sigma_{intra}^{res}({\bf k},\omega<-\Delta-E_{\rm F})$
$\displaystyle=$
$\displaystyle-\frac{e^{2}}{2\pi\kappa\sqrt{k^{2}+\Delta^{2}}}\int_{max(0,k-k_{\rm
F})}^{k+k_{\rm F}}dq\int^{min(\sqrt{k_{\rm
F}^{2}+\Delta^{2}},\sqrt{(k+q)^{2}+\Delta^{2}})}_{max(0,\omega+\sqrt{k_{\rm
F}^{2}+\Delta^{2}},\sqrt{(k-q)^{2}+\Delta^{2}})}dy$ (48) $\displaystyle\times$
$\displaystyle\frac{(\sqrt{k^{2}+\Delta^{2}}+y)^{2}-q^{2}}{\epsilon({\bf
q},\omega+E_{\rm F}-y)\sqrt{4k^{2}q^{2}-(y^{2}-k^{2}-q^{2}-\Delta^{2})^{2}}}$

The real and imaginary part of intraband contributions can be computed.

##  Appendix C The interband contribution of self-energy

Now, we focus on the interband contribution of the retarded self-energy. The
second argument of the dielectric function in Eq. (5) is

$\displaystyle\omega-\xi_{-}({\bf k+q})=\omega+E_{\rm
F}+\sqrt{k^{2}+q^{2}+2kq\cos\phi+\Delta^{2}}~{}.$ (49)

We change variable $y=\sqrt{k^{2}+q^{2}+2kq\cos\phi+\Delta^{2}}$, then we find

$\displaystyle\Sigma_{inter}^{res}({\bf k},\omega)$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{2\pi\kappa\sqrt{k^{2}+\Delta^{2}}}\int_{0}^{+\infty}dq\int_{\sqrt{(k-q)^{2}+\Delta^{2}}}^{\sqrt{(k+q)^{2}+\Delta^{2}}}\frac{dy}{\sqrt{4k^{2}q^{2}-(y^{2}-k^{2}-q^{2}-\Delta^{2})^{2}}}\frac{q^{2}-(y-\sqrt{k^{2}+\Delta^{2}})^{2}}{\epsilon({\bf
q},\omega+E_{\rm F}+y)}$ (50) $\displaystyle\times$
$\displaystyle[\Theta(\omega+E_{\rm F}+y)-1].$

Note that $\Sigma_{inter}^{res}$ can be non-zero if $\omega+E_{\rm
F}+y<0~{}~{}and~{}~{}y>0$. It means that $\omega<-E_{\rm F}$. In this case the
difference between the two $\Theta$-functions in Eq. (5) becomes -1 if

$\displaystyle 0<y<-(\omega+E_{\rm
F})~{}~{}and~{}~{}\sqrt{(k-q)^{2}+\Delta^{2}}<y<\sqrt{(k+q)^{2}+\Delta^{2}}~{}.$
(51)

Now we do need to find the overlap between these two intervals. We end up to
inequivalent conditions that $q>k-\sqrt{\omega^{2}+k_{\rm
F}^{2}+2\omega\sqrt{k_{\rm F}^{2}+\Delta^{2}}}$ and
$q<k+\sqrt{\omega^{2}+k_{\rm F}^{2}+2\omega\sqrt{k_{\rm F}^{2}+\Delta^{2}}}$.
Putting everything together and using the fact that $q\geq 0$ we finally get

$\displaystyle\Sigma_{inter}^{res}({\bf k},\omega<-E_{\rm F})$
$\displaystyle=$
$\displaystyle-\frac{e^{2}}{2\pi\kappa\sqrt{k^{2}+\Delta^{2}}}\int_{max(0,k-\sqrt{\omega^{2}+k_{\rm
F}^{2}+2\omega\sqrt{k_{\rm F}^{2}+\Delta^{2}}})}^{k+\sqrt{\omega^{2}+k_{\rm
F}^{2}+2\omega\sqrt{k_{\rm
F}^{2}+\Delta^{2}}}}dq\int_{\sqrt{(k-q)^{2}+\Delta^{2}}}^{min(\sqrt{(k+q)^{2}+\Delta^{2}},-(\omega+\sqrt{k_{\rm
F}^{2}+\Delta^{2}}))}dy$ (52) $\displaystyle\times$
$\displaystyle\frac{q^{2}-(y-\sqrt{k^{2}+\Delta^{2}})^{2}}{\epsilon({\bf
q},\omega+\varepsilon_{\rm
F}+y)\sqrt{4k^{2}q^{2}-(y^{2}-k^{2}-q^{2}-\Delta^{2})^{2}}}.$

If we want to calculate $\Im m\Sigma^{res}_{+}({\bf k},\xi_{+}(\bf k))$ needed
for computing the quasiparticle lifetime, we will only need the intraband
contribution of the self-energy since the interband contribution is zero.

## References

  * (1) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, A. A. Firsov, Science 306, 666 (2004); K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos and A. A. Firsov, Nature 438, 197 (2005); K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V. Khotkevich, S. V. Morozov, and A. K. Geim, Proc. Nat. Acad. Sci. 102, 10451 (2005) .
  * (2) Andreas Barth, Werner Marx, arXiv:0808.3320 (2008) .
  * (3) A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao, C. N. Lau, Nano Lett., 8, 902 (2008) .
  * (4) S.V. Morozov, K.S. Novoselov, M.I. Katsnelson, F. Schedin, D.C. Elias, J.A. Jaszczak, A.K. Geim, Phys. Rev. Lett. 100, 016602 (2008); K. I. Bolotin, K. J. Sikes, J. Hone, H. L. Stormer, and P. Kim, Phys, Rev, Lett. 101, 096802 (2008); Xu Du, Ivan Skachko, Anthoy Barker and Eva Y. Andrei, Nature Nanotech. 3, 491 (2008); K.I. Bolotin, K.J. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, H.L. Stormer, Solid State Commun. 146, 351 (2008) .
  * (5) K. Eng, R. N. McFarland, and B. E. Kane, Appl. Phys. Lett. 87, 052106 (2005); E. H. Hwang and S. Das Sarma, Phys. Rev. B 75, 073301 (2007) .
  * (6) T. O. Wehling, K. S. Novoselov, S. V. Morozov, E. E. Vdovin, M. I. Katsnelson, A. K. Geim, A. I. Lichtenstein, Nano Lett. 8, 173-177 (2008); S.Y. Zhou, D.A. Siegel, A.V. Fedorov, A. Lanzara, Phys. Rev. Lett. 101, 086402 (2008); Isabella Gierz, Christian Riedl, Ulrich Starke, Christian R. Ast, Klaus Kern, arXiv:0808.0621 .
  * (7) 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, Nat. Mater. 6, 770 (2007) .
  * (8) G. Li, A. Luican, and E. Y. Andrei, arXiv:0803.4016 (2008) .
  * (9) Gianluca Giovannetti, Pet A. Khomyakov, Geert Brocks, Paul J. Kelly, and Jeoroen Van den Brink, Phys. Rev. B 76, 073103 (2007) .
  * (10) R. M. Ribeiro, N. M. R. Peres, J. Coutinho, and P. R. Briddon, Phys. Rev. B 78, 075442 (2008) .
  * (11) I. Zanella, S. Guerini, S. B. Fagan, J. Mendes Filho, and A. G. Souza Filho, Phys. Rev. B 77, 073404 (2008) .
  * (12) Yugui Yao, F. Ye, Xiao-Liaug Qi, Shou-Cheng Zhang and Zhang Fang, Phys. Rev. B 75, 041401(R) (2007); C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005); Hongki Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, Leonard Kleinman, and A. H. MacDonald, Phys. Rev. B 74, 165310 (2006) .
  * (13) Y. W. Son, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 97, 216803 (2006); M. Y. Han, B. Ozyilmaz, Y. Zhang, and P. Kim, ibid 98, 206805 (2007) .
  * (14) L.D. Landau, Sov. Phys. JEPT, 3, 920 (1957) .
  * (15) I.K. Marmorkos and S. Das Sarma, Phys. Rev. B 44, R3451 (1991); J.D. Lee and B.I. Min, Phys. Rev. B 53, 10 988 (1996); H.-J. Schulze, P. Schuck, and N. Van Giai, Phys. Rev. B 61, 8026 (2000) .
  * (16) S. Yarlagadda and G.F. Giuliani, Phys. Rev. B 49, 7887 (1994); 61, 12556 (2000); C.S. Ting, T.K. Lee, and J.J. Quinn, Phys. Rev. Lett. 34, 870 (1975).
  * (17) H. M. Böhm and K. Schörkhuber, J. Phys.: Condens. Matter 12, 2007 (2000) .
  * (18) Y. Zhang and S. Das Sarma, Phys. Rev. B 71, 045322 (2005); S. Das Sarma, V. M. Galitski, and Y. Zhang, Phys. Rev. B 69, 125334 (2004) .
  * (19) Y. Zhang, V.M. Yakovenko, and S. Das Sarma, Phys. Rev. B 71, 115105 (2005) .
  * (20) R. Asgari, B. Davoudi, M. Polini, Gabriele F. Giuliani, M. P. Tosi, and G. Vignale, Phys. rev. B 71, 045323 (2005); R. Jalabert and S. Das Sarma, Phys. Rev. B 40, 9723 (1989) .
  * (21) R. Asgari, B. Davoudi and B. Tanatar, Solid State Commun. 130, 13 (2004).
  * (22) Jahan M. Dawlaty, Shriram Shivaraman, Mvs Chandrasekhar, Farhan Rana, Micheal G. Spencer, App. Phys. Lett 92, 042116 (2008) .
  * (23) J. R. Williams, L. DiCarlo and C. M. Marcus Since 317, 638 (2007) .
  * (24) G. Gu, S. Nie, R. M. Feenstra, R. P. Devaty, W. J. Choyke, W. K. Chan and M. G. Kane, Appl. Phys. Lett. 90, 253507 (2007) .
  * (25) Marco Polini, Reza Asgari, Giovanni Borghi, Yafis Barlas, T. Pereg-Barnea, and A.H. MacDonald, Phys. Rev. B 77, 081411(R) (2008); E. H. Hwang and S. Das Sarma, Phys. Rev. B 77, 081412(R) (2008) .
  * (26) Z. Q. Li, E. A. Henriksen, Z. Jiang, Z. Hao, M. C. Martin, P. Kim, H. L. Stromer, D. N. Basov, Nature Phys. 4, 532 (2008) .
  * (27) M. Polini, R. Asgari, Y. Barlas, T. Pereg-Barnea, A. H. MacDonald, Solid State Commun. 143, 58 (2007) .
  * (28) G. F. Giuliani and G. Vignale, Quantum Theory of The Electron Liquid (Cambridge University Press, Cambridge, England, 2005).
  * (29) A. Qaiumzadeh and R. Asgari, arXiv:0807.3183 (2008) .
  * (30) P. K. Pyatkovskiy, arXiv:0808.0931 (2008) .
  * (31) E. H. Hwang, BenYu-Kaung Hu, and S. Das Sarma, Phys. Rev. B 76, 115434 (2007) .
  * (32) R. Asgari and B. Tanatar, Phys. Rev. B 65 (2002) 085311 .
  * (33) Y. Barlas, T. Pereg-Barnea, M. Polini, R. Asgari and A. H. MacDonald, Phys. Rev. Lett. 98 (2007) 236601 .
  * (34) A. Qaiumzadeh, N. Arabchi, R. Asgari, Solid State Commun. 147, 172 (2008) .
  * (35) G.F. Giuliani and J. J. Quinn, Phys. Rev. B 26, 4421 (1982) .
  * (36) P. K. Pyatkovskiy, Private communication .

Figure 1: ( Color online) (a): The renormalized velocity as a function of
density for the various energy gaps at $\alpha_{gr}=1$. (b): The renormalized
velocity as a function of the energy gap for the various densities.  Figure 2:
( Color online) The renormalization constant as a function of the energy gap
for the various densities.

Figure 3: ( Color online) (a): The absolute value of the imaginary part of the
retarded self-energy on the energy shell as a function of the wavevector for
the various energy gaps; (b): The inelastic quasiparticle
lifetime($\tau_{in}$) in graphene as a function of the on-shell energy for the
various energy gaps at $n=5\times 10^{12}$cm-2.  Figure 4: ( Color online) The
quasiparticle mean free path as a function of the on-shell energy for the
various energy gaps at $n=5\times 10^{12}$cm-2.


# Ultra-low latency recurrent neural network inference on FPGAs for physics
applications with hls4ml

Elham E Khoda1 , Dylan Rankin2 , Rafael Teixeira de Lima3
Philip Harris2, Scott Hauck1, Shih-Chieh Hsu1, Michael Kagan3, Vladimir
Loncar4,
Chaitanya Paikara1, Richa Rao1, Sioni Summers4, Caterina Vernieri3, Aaron
Wang1

1University of Washington, Seattle, WA 98195, USA
2Massachusetts Institute of Technology Cambridge, MA 02139, USA
3SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
4European Organization for Nuclear Research (CERN), Geneva, Switzerland
ekhoda@uw.edudrankin<EMAIL_ADDRESS>
###### Abstract

Recurrent neural networks have been shown to be effective architectures for
many tasks in high energy physics, and thus have been widely adopted. Their
use in low-latency environments has, however, been limited as a result of the
difficulties of implementing recurrent architectures on field-programmable
gate arrays (FPGAs). In this paper we present an implementation of two types
of recurrent neural network layers—long short-term memory and gated recurrent
unit—within the hls4ml framework. We demonstrate that our implementation is
capable of producing effective designs for both small and large models, and
can be customized to meet specific design requirements for inference latencies
and FPGA resources. We show the performance and synthesized designs for
multiple neural networks, many of which are trained specifically for jet
identification tasks at the CERN Large Hadron Collider.

##  1 Introduction

Machine learning (ML) has seen a huge expansion in its range of uses over the
last decade. It is difficult to find a field of industry or science that has
not at least explored ML in some capacity. One particular field where ML usage
has seen widespread interest is in high energy physics, which benefits from
complex multidimensional problems, large datasets of accurate simulation, and
substantial existing computing infrastructure. These all contribute to a field
which has adopted ML algorithms for many aspects of research. While most ML
algorithms in high energy physics are run using CPUs and graphics processing
units (GPUs) which provide inference latencies in the milliseconds, field-
programmable gate arrays (FPGAs) and application-specific integrated circuits
(ASICs) have begun to be used for those applications that demand low latencies
[1, 2, 3, 4]. Up to now recurrent neural networks (RNNs) have received
relatively little attention for these low latency applications, despite their
success in many physics tasks and prevalence in the field at large.

RNNs are neural network architectures that treat their inputs as a sequence
with a well-defined order. RNNs act successively on each entry of the input
sequence, utilizing the same set of weights in each step, allowing for RNNs to
act on sequences of variable length. Most modern RNN applications utilize one
of two recurrent layer implementations: long-short term memory layer (LSTMs)
[5] or gated recurrent unit layers (GRUs) [6]. Both LSTMs and GRUs contain
forget gates, which avoid vanishing gradient problems [7] and allows for long-
distance correlations to be learned by the algorithm. These RNN models, often
employed to treat time-series signal processing problems, have been
successfully employed by physicists to handle different types of data. For
example, RNN architectures have been utilized to represent hadronic showers
(jets) in collider events, aiding on tasks of identifying different jet types
by using their constituents to form sequences [8, 9, 10]; have been applied to
finding interaction vertices in lepton colliders [11]; and have been explored
as monitoring tools for the CERN Large Hadron Collider (LHC) superconducting
magnets [12]. In general, RNNs have critical applications also outside the
realm of collider physics, having been applied, among others, to gravitational
waves detection experiments [13], to neutrino detectors from nuclear reactors
[14], and to the reconstruction of quantum dynamics of superconducting qubits
[15].

In this paper, we focus on particle physics datasets produced by the LHC at
CERN [16]. Collisions within the LHC occur at 40 MHz, making it impossible to
readout and store the entire collisions record. To handle this rate, LHC
experiments filter out only interesting events through an online selection
system called the trigger. These systems are usually set in two stages, where
the first stage (Level-1 trigger or L1T) needs to operate at 40 MHz with a
latency of O(1 $\mu$s). The selections performed at these stages need to
ensure that interesting events are kept, while discarding common, non-
interesting events, which occur several orders of magnitude more frequently
than the former. Therefore, utilizing complex algorithms such as RNNs is of
utmost importance.

The severe constraints of the trigger prevent the usage of CPUs and GPUs.
Instead, custom low-latency hardware such as FPGAs and ASICs must be deployed
to meet the latency requirements and offer the flexibility to adapt to
changing conditions. These devices are also able to take advantage of high
parallelism making their designs both efficient and fast. ML inference in this
regime has not seen much support due to its specialized nature, but some tools
specifically designed for ultra-low latency inference have emerged [17, 18].
Support in this area has focused primarily on dense and convolutional layers,
owing to their versatility and popularity. In this work we present a
generalized and flexible implementation of RNNs written in High Level
Synthesis (HLS) for the hls4ml package [19]. The implementation supports a
wide range of RNN sizes, design requirements, and is capable of translating
both GRUs and LSTMs trained in the Keras framework [20]. Using three different
benchmark neural network models of varying size, we show that ultra-low
latency inference can be achieved within the resource limitations of modern
FPGAs. This integration into hls4ml opens the door for much wider usage of
RNNs for low latency applications. While the focus in this work is on
applications in physics, we note that there is also a demand for low-power
efficient RNNs in industry as well.

##  2 Related Work

Previous work in the realm of fast RNN inference on FPGAs has focused largely
on millisecond-latency inferences [21, 22, 23, 24]. However, for the
applications discussed above we are interested primarily in latencies in the
microsecond range, or faster. Some previous work has explored this design
space in the context of low-power sparse LSTMs [25], small LSTMs for real-time
energy reconstruction [26], RNNs for gravitiation-wave experiments [27], and
highly quantized RNNs [28]. In contrast, the work in this paper is focused on
general support for both large and small LSTMs and GRUs for problems with a
range of latency and device constraints. The examples we use are largely
chosen from the high energy physics domain, but the applicability is by no
means limited to this field. The work is built on top of high level synthesis
(HLS) and the hls4ml framework.

HLS tools are designed to simplify the use of FPGAs by automatically
transforming algorithms written in C into the register-transfer level (RTL)
[29]. There are multiple advantages of these tools. One substantial advantage
is that they allow users without a knowledge of highly technical Verilog/VHDL
langauges to generate effective RTL [30]. Additionally, they can greatly
simplify the effort required in prototyping designs, especially those that are
complex. FPGA manufacturers like Xilinx and Altera have their own HLS
compilers for their devices. There also exist open-source compilers, such as
Catapult HLS [31]. In this work we use the Vivado/Vitis compiler from Xilinx
[32, 33].

The hls4ml framework is built on top of HLS compilers, and is capable of
converting neural network models into fully-ready HLS projects [17]. The
details of the HLS design, in particular the resources and latency, can be
controlled through multiple tunable parameters in hls4ml. These are important
to allow a flexible design flow that performs well for a wide range of network
sizes, architectures, and FPGAs. They also allow the design to be optimized
for the target use case. hls4ml already has support for MLPs, CNNs, graph NNs,
and several other architectures. Building on hls4ml allows for models with
these architectures to be interfaced with the models in this paper.
Furthermore, extensive work has been done to hls4ml to ensure that matrix
multiplications and other core ML components are optimized. Our work uses the
hls4ml design flow in order to leverage these existing framework capabilities.

##  3 Implementation Details

Recurrent neural networks at their core are comprised of many standard
operations in machine learning. Our implementation relies on this fact to
avoid re-implementing any unnecessary operations in the hls4ml framework.
Taking the LSTM, we see that each state update requires 4 distinct matrix
multiplications, given by

$\displaystyle i_{t}$ $\displaystyle=\sigma(W_{i}x_{t}+U_{i}h_{t-1}+b_{i})$
$\displaystyle f_{t}$ $\displaystyle=\sigma(W_{f}x_{t}+U_{f}h_{t-1}+b_{f})$
$\displaystyle o_{t}$ $\displaystyle=\sigma(W_{o}x_{t}+U_{o}h_{t-1}+b_{o})$
$\displaystyle c_{t}$ $\displaystyle=\tanh(W_{c}x_{t}+U_{c}h_{t-1}+b_{c})$ (1)

where $W$ and $U$ are the weight matrices (denoted as the kernel and recurrent
kernel, respectively), $b$ are the biases for the input gate, forget gate,
output gate, and cell state, respectively, $x_{t}$ is the input at time-step
$t$, and $h_{t-1}$ is the computed recurrence value from the previous state.
Each of the operations involving $W$ and $U$ in Eq. 1 are standard matrix-
vector multiplications and the activation functions can be taken directly from
the hls4ml library and integrated into the LSTM-specific layer
implementations. The GRU is composed of two gates (update and reset) and a
hidden state, but the weights of the kernel and recurrent kernel for the gates
are again packaged together and can thus be handled together with one dense
layer call each. The remaining operations to complete a state update for both
an LSTM and a GRU are Hadamard products, which is not part of the existing
hls4ml library of operations. In this paper, we implemented an HLS-optimized
Hadamard product. Because many of the operations in our implementation are
taken from the existing hls4ml framework, we are able to trivially support the
standard tuning knobs for reuse and precision.

In addition to the pre-existing hls4ml methods for adjusting resources and
latency, RNNs introduce an additional possibility which we refer to as
“static” and “non-static”, shown in Fig. 1. In static mode, a single RNN block
is created, which processes every input for every sequence. This block stores
the necessary state vectors internally, and outputs the final result at the
end of the sequence. Since there is only one RNN block in static mode, the
resources are kept to a minimum. However, the initiation interval (II) of the
design increases linearly with the length of the sequence since a new RNN
inference cannot begin until the previous inference is complete; in other
words, the II is equal to the latency. In contrast to static mode, non-static
mode creates RNN blocks for each input in the sequence, and the necessary
state vectors are passed from one RNN block to another. This results in a
resource utilization that is a factor of the sequence length larger than the
resource utilization in static mode. For large RNNs, or inference with long
sequences, this means that static mode can be the only viable option. However,
for those RNNs for which it is possible to use non-static mode, the II can be
dramatically reduced with respect to static mode since non-static mode allows
a new inference to begin once the first RNN block has finished processing the
first input from the previous inference. This reduces the II by a factor of
the length of the sequence, and thus increases the overall throughput of the
RNN layer by the same factor. Further increases in throughput are possible
when the II of a single RNN block can be made small since in this case each
individual block may be used simultaneously by distinct inferences. While not
implemented in this paper, we note that multiple inferences can be cached
during static mode when the initiation interval of a single RNN block is less
than its latency, thus allowing for higher throughput.

Figure 1: Schematics of the two RNN modes: static (left) and non-static
(right). A sketch of the latency and pipelining is shown at the bottom for
each mode, and illustrates the latency advantage of non-static mode at the
cost of resources.

##  4 Benchmark Studies

In order to measure the performance of the hls4ml translation of the RNN-based
architectures, we use three different tasks as benchmarks. The sizes of the
networks are chosen such that they span a range of use cases, input
dimensions, and numbers of weights. The first benchmark is a binary classifier
with approximately 4000 parameters, trained to classify jets of particles
coming from top-quark decays. The second benchmark is a multi-class
classifier, with approximately 50,000 parameters, trained for heavy-flavor jet
identification using the reconstructed trajectories of charged particles
(tracks) within a jet. The final benchmark is also a multi-class classifier
with approximately 130,000 parameters, trained to classify sequences of
strokes into five different image classes. Each benchmark is studied with two
models using LSTM and GRU recurrent layers, respectively. All the models are
trained using Keras and TensorFlow. A summary of each of these models is given
in Table 1, and details are given in the following sections.

Table 1: Network hyperparameters and total number of trainable parameters for different benchmark models. Benchmark | Sequence | Input | Hidden | Dense | output | Trainable parameters
---|---|---|---|---|---|---
length | vector | vector | layer | vector | Non-RNN | LSTM | GRU
| size | size | sizes | size | layers |  |
Top tagging | 20 | 6 | 20 | 64 | 1 | 1,409 | 2,160 | 1,680
Flavor tagging | 15 | 6 | 120 | 50/10 | 3 | 6,593 | 60,960 | 46,080
QuickDraw | 100 | 3 | 128 | 256/128 | 5 | 66,565 | 67,584 | 51,072

###  4.1 Top Quark Tagging

The top quark tagging algorithm is trained to classify top quarks from light-
flavor quarks using simulated events generated at $\sqrt{s}=13\,\text{TeV}$
for comparison to LHC performance. Algorithms designed for this task could be
utilized in the Level-1 trigger systems of LHC experiments to help increase
the acceptance for these types of interesting decays. Their use would require
algorithm latencies of less than approximately a few microseconds in order to
fit within system constraints.

The data used for training and testing consists of parton-level scattering
processes with top quark-antiquark pair ($\mathrm{t}\overline{\mathrm{t}}$)
and light-flavor quark-antiquark pair ($\mathrm{q}\overline{\mathrm{q}}$)
final states that are generated at leading-order using MadGraph [34] with the
NNPDF23LO1 parton distribution functions (PDFs) [35]. The transverse momenta
($p_{\mathrm{T}}$) of the partons are generated in a window with energy spread
given by $\delta p_{\mathrm{T}}/p_{\mathrm{T}}=0.01$, centered at 1 TeV. These
parton-level events are then decayed and showered using Pythia8 [36] (version
8.212) with the Monash 2013 tune [37], including the contribution from the
underlying event. Jets are clustered using the anti-$k_{\mathrm{T}}$
algorithm, with a distance parameter of $0.8$. Only low level features are
used in this study. Particles inside a jet are ordered according to their
$p_{\mathrm{T}}$ and up to 20 particles are used in this study. For each
particle six features are use: $p_{\mathrm{T}}$, pseudorapidity ($\eta$),
azimuthal angle ($\phi$), energy, relative angular distance from the jet axis,
particle ID given by the generator. The generated events are split into
training (95%) and testing (5%), and during training 20% of the training data
is used for validation.

Two networks are trained for identifying the jets coming from top-quark decay.
The padded sequence of particles, with maximum length of 20, is fed into a
recurrent layer with an output size of 20. The output from the final recurrent
layer is passed through a dense layer of size 64 before sending it to the
output layer. The recurrent layer used sigmoid and hyperbolic tangent
activation functions. The activation function for the hidden layers is ReLU
[38] while the output layer activation function is a sigmoid function. The
binary cross-entropy loss function is minimized with L1 ($10^{-5}$) and L2
($10^{-4}$) regularization of the weights using the Adam algorithm [39] with a
learning rate of 2$\times 10^{-4}$ and a batch size of 246. The two models use
different recurrent layers; one uses GRU and the other uses LSTM. There are
total 3,089 and 3,569 trainable parameters for the LSTM and GRU models,
respectively.

###  4.2 Jet Flavor Tagging

The jet flavor tagging algorithm was trained on CMS Open Data samples
containing top quark pairs decaying hadronically with center-of-mass energy of
7 TeV [40]. These events are rich in bottom quark jets (b jets), charm quark
jets (c jets) and jets from light quarks and gluons (light jets), and so are
optimal for training this class of algorithms. Jets are labeled b jets if they
contain bottom quarks, c jets if they do not contain bottom quarks but contain
charm quarks, and light jets if they do not contain bottom or charm quarks.
The main feature that separates b jets (and c jets) from light jets is the
presence of the displaced vertex corresponding to the decay of the hadron
containing the b (or c) quark. These hadrons can travel significant distances
before decaying due to their mass, depending on their momenta. The algorithm
proposed here aims to identify the presence of tracks that are consistent with
these displaced vertices with the usage of an RNN architecture. This strategy
is inspired by the RNNIP algorithm used by the ATLAS experiment [8].

In this study, we consider jets reconstructed with the anti-kt algorithm with
a distance parameter of $0.5$, with transverse momenta larger than 30 GeV and
absolute pseudorapidities less than 2.0. Tracks with transverse momenta larger
than 1 GeV are associated to the nearest jet according to the angular distance
$\Delta R$; a maximum $\Delta R$ of 0.5 is required for association. Tracks
within a jet are ordered by the significance of their transverse impact
parameter ($\mathcal{S}(d_{0})$), and only the first 15 tracks are used by the
algorithm. Each track is represented by a vector of the following features:
relative transverse momentum
($p_{\mathrm{T}}(\mathrm{track})/p_{\mathrm{T}}(\mathrm{jet})$); $\Delta
R(\mathrm{track,jet})$; transverse and longitudinal impact parameters
($d_{0}$, $d_{z}$) and their significances ($\mathcal{S}(d_{0})$,
$\mathcal{S}(d_{z})$).

Flavor tagging models were constructed using Keras/TensorFlow, using either
GRU or LSTM layers. The padded sequence of tracks, with maximum length of 15,
is fed into either one recurrent layer (GRU or LSTM) with 120 hidden units.
The recurrent layer outputs a representation of the padded sequence directly
into two dense layers with ReLU activation function, with 50 and 10 hidden
layers. The following layer outputs the probabilities of a jet to be
classified as either a b jet, c jet or light jet; it contains three output
nodes with softmax activation function. The training is performed with a
categorical cross-entropy loss, with 30% of the training data retained as the
validation dataset and used for early stopping based on the accuracy metric.
The GRU (LSTM) architecture contains 52,673 (67,553) trainable parameters, of
which 46,080 (60,960) are in the recurrent layer.

###  4.3 Quickdraw Dataset

The QuickDraw dataset [41] is a collection of 50 million drawings in 345
categories created by Google and contributed by players of the game Quick!
Draw. In this game, users are asked to draw a specified drawing in under 15
seconds. The drawings are recorded as a time-stamped sequence of the strokes
from which the drawing is created. For each 15 second stroke the $x$ and $y$
coordinates of the pen are recorded 100 times. The coordinates along with the
timestamp make up the network inputs. While 345 different drawing categories
exist, we use only 5 for our tests; these are ants, butterflies, bees,
mosquitos and snails. Contrary to other popular representations of images, the
QuickDraw dataset is completely stroke-based. We train two recurrent neural
networks to classify images into each respective category. Similar algorithms
could be used to identify tracks based on the sequence of their hits, or
incident particles based on their showers in finely-segmented calorimeters.
While not appropriate for the Level-1 trigger environment, running these
algorithms at later stages in the trigger for a whole event could still
require low latencies under a millisecond depending on the exact application.

The two networks process the sequence of 100 stroke inputs from the QuickDraw
dataset with a recurrent layer whose output size is 128. The final recurrent
layer output is passed through two dense layers of sizes 256 and 128,
respectively, before being sent to the final output layer. Dropout layers are
placed before the two dense layers to regularize during training. The
recurrent layer uses a hyperbolic tangent activation function, the dense
layers use ReLU activations, and the output is a five-class softmax layer. The
only difference in the two networks is that in one the recurrent layer is a
GRU and in the other it is an LSTM. These networks have total sizes of 117,637
and 134,149 trainable parameters, respectively. The two networks perform well
and show top-1 AUCs that are nearly identical, approximately 99% for each of
the five classes.

##  5 Performance, Resource and Latency Estimation

The models described in Sec. 4 are translated into HLS using the hls4ml
framework. Vivado HLS 2019.2 is used for HLS synthesis with the synthesis
clock frequency set to 200 MHz. For the top tagging and flavor tagging models
we use a Xilinx Kintex UltraScale FPGA (part number xcku115-flvb2104-2-i) as
the target device and for the Quickdraw models we use a Xilinx Alveo U250
(part number xcu250-figd2104-2-e). For each of the three benchmark models a
range of different settings are considered simultaneously to accurately
profile the full design space. These settings modify the quantization of the
model by adjusting the fixed-point data type and modify the degree of
parallelism of the design by using the hls4ml reuse parameter. Finally, we
also consider the static and non-static implementations discussed in Sec 3.

###  5.1 Quantization

The weights and biases in trained models are typically stored with 32-bit
floating-point precision. However, 32-bit floating-point calculations are
often not required for optimal network inference, and are costly to implement
on FPGAs. Other quantization techniques can offer more efficient ways of
compressing neural networks by reducing the number of bits used to represent
the weights and biases, ideally with no or minimal loss in performance. In
hls4ml, all the inputs, weights, biases, sums, and outputs of each layer are
represented as fixed-point numbers. In this scheme the amount of bits used to
store the integer and decimal components of the number are configured, such
that, for example, an unsigned fixed point number with 4 integer bits and 3
fractional points is capable of storing values between 0 and 15.875 with a
granularity of 0.125. The total numbers of bits is also referred to as the
precision of the fixed-point number. Hls4ml allows a different precision to be
chosen for the computations and internal values of each individual layer; for
the sake of consistency we fix the precision to be the same for all layers in
the scans below. We do find that it is necessary to increase the precision and
size of the LUT used for the softmax computation at the end of the flavor-
tagging and QuickDraw models, but this has a minimal impact on the overall
resource usage.

The optimal precision for each model depends on the training details, the
specific task, and the inputs. All the models are quantized only after
training, a method referred to as “post-training quantization” (PTQ). For each
model we profile the performance of the synthesized design from hls4ml as a
function of the bit precision of the weights and activation functions. Figure
2 shows the ratio of the AUC from the quantized model to the AUC from the
floating-point model as a function of fractional bits while keeping the
precision of the integer part fixed to 6, 8, 10, or 12 bits.

(a) Top quark tagging

(b) Jet flavor tagging

(c) QuickDraw

Figure 2: Ratios of the fixed-point and floating-point AUCs as function of
fractional bits for the (a) top quark tagging, (b) jet flavor tagging, and (c)
QuickDraw models. The precision of the integer part is kept fixed to 6 (blue),
8 (orange), 10 (green), and 12 (red) bits. All four lines overlap for the top-
tagging GRU (dashed) and LSTM (solid) models.

The best performance for each model, measured by the AUC ratio, is generally
achieved with at least 10 fractional bits irrespective of the values of the
integer bit. For the top quark and flavor tagging models, 6 integer bits are
sufficient, while the QuickDraw models require at least 10 integer bits. For
further results in this paper, we fix the integer bits for each model to these
values. We note that there is a small performance degradation in the GRU
models after quantization for all three benchmark cases. The difference is
particularly visible for the Top Tagging model, but it is less than 5%. It is
possible that quantization-aware training or sequence masking could
potentially enable models with lower precision to perform as well as or better
than the ones we present in this paper.

###  5.2 Parallelization

The other main tuning knob besides the precision is the amount of parallelism
employed during weight matrix multiplication. This is controlled in hls4ml
through a parameter called “reuse”. Specifically, reuse is the number of
multiplication operations each digital signal processing (DSP) block must do
for a given matrix multiply. Setting reuse to 1, i.e. the fully parallel case,
means that each multiplication is done by its own DSP and can happen
simultaneously. Increasing the reuse factor reduces the number of DSPs that
are required, but increases the latency and initiation interval of the layer
computation in proportion to the reuse. All three benchmark models are
synthesized with different values of the reuse factor ($R$) and fractional bit
precision. The results are expressed for different FPGA resource categories:
onboard FPGA memory (BRAM), DSPs, and registers and programmable logic like
flip-flops (FFs) and lookup tables (LUTs). In hls4ml, a model can either be
synthesized to minimize the latency (latency strategy) or the resource
utilization (resource strategy). For large models with 40k or more trainable
parameters it becomes difficult to synthesize the models with the latency
strategy, and so resource strategy must be used. With resource strategy the
design is optimized for low resource utilization by reusing existing hardware
to complete operations in multiple stages. Out of the three benchmark models
only the Top Tagging model is small enough to be synthesized with both latency
and resource strategies, whereas only resource strategy is used for the other
two models. The minimum and maximum latencies for each model are shown in
Table 2, 3, and 4. The amounts of DSPs, FFs, and LUTs for each model are shown
for different $R$ values in Fig. 3, 4, and 5, respectively. In these results
the reuse factor values are written in the form $R=(X,Y)$, where $X$ and $Y$
correspond to the reuse factors for the kernel and recurrent kernel matrix
multiplications discussed in Eq. 1. The the numbers shown in the square
brackets correspond to reuse factor for the LSTM layer in the case when the
reuse factor differs between the LSTM and GRU implementations.

We observe that all resources generally increase with smaller values of $R$
and increased precision. In the case of FFs and LUTs this increase is roughly
linear, while for DSPs the utilization remains flat until the precision
exceeds the DSP input width. The latency, on the other hand, follows a scaling
inverse to that of the FFs and LUTs with respect to the reuse. Thus, as with
other architectures supported under hls4ml, reuse can be used to reduce FFs
and LUTs at the expense of latency. This simple scaling is critical for
allowing users to tune the resource usage and latency such that the
synthesized designs to meet desired requirements. The latency strategy adds
another finer option to this tuning space for latency-limited tasks, but comes
at the cost of larger resource usage. As expected we find that the GRU models
use approximately $1/4$ less resources when compared to the LSTM models. This
is a result of the 3:4 ratio between the number of matrix multiplications in
GRU and LSTM models. Finally, it is important to note that the results shown
in this paper are from HLS synthesis. When running Vivado synthesis we observe
a reduction in LUT usage between 20% and 65% and in FF usage between 10% and
20%. This is particularly important to note in the case of the larger flavor
tagging and QuickDraw models where the estimated LUT usages from HLS synthesis
are quite large.

For the top tagging models we observe that designs with maximal quantized
performance can be implemented on one SLR of a Xilinx Virtex Ultrascale+ VU9P,
the planned future device for an upgrade to the CMS Level-1 trigger [42]. We
observe slightly larger resource usage for the flavor tagging models, as
expected, but still within the resource constraints of a single SLR of a VU9P.
In both cases the latencies for the designs are also within the task
requirements. For the QuickDraw models, the estimates from Vivado synthesis
suggest that maximal quantized performance could be implemented on a Xilinx
Alveo U250, a popular device for the types of coprocessor applications we
envision for these models. Extrapolating from the II, we find the average
throughput of the QuickDraw LSTM model is between 4300 to 9700 events/sec.
Tests of the batch 1 inference for the same model using an Nvidia Tesla V100
GPU yield a throughput of 660 events/sec. Increasing the batch size to 10
increases the throughput to 7700 events/sec, comparable with the FPGA
throughput. The throughput increases further by a factor of five to
approximately 30000 if the batch size is increased to 100. While it is
unsurprising that GPU inference at large batch sizes is able to outperform an
FPGA, many physics tasks are inherently low-batch problems. This is because
each event must be processed separately and latency is extremely important,
therefore the maximal batch size is dictated by the amount of inferences
necessary only for a single event. For example, algorithms to classify
particle-induced showers in a detector need only to be run once every event if
the algorithm can use the full detector information in one pass. Thus, a
factor of ten improvement in the FPGA inference for batch-1 inference is
highly relevant for future trigger applications.

(a) Top quark tagging

(b) Jet flavor tagging

(c) QuickDraw

Figure 3: DSP utilization as a function of total width for the (a) top quark
tagging, (b) jet flavor tagging (c) QuickDraw model. Performance of both GRU
(dashed) and LSTM (solid) models are shown. DSP utilization with latency
strategy is shown only for the top quark tagging models and all other lines
correspond to different reuse factors. The DSPs available in the target FPGA
for each model are shown in the black dashed horizontal line.

(a) Top quark tagging

(b) Jet flavor tagging

(c) QuickDraw

Figure 4: FF utilization as a function of total width for the (a) top quark
tagging, (b) jet flavor tagging (c) QuickDraw model. Performance of both GRU
(dashed) and LSTM (solid) models are shown. FF utilization with latency
strategy is shown only for the Top Tagging models and all other lines
correspond to different reuse factors. The FFs available in the target FPGA
for each model are shown in the black dashed horizontal line.

(a) Top quark tagging

(b) Jet flavor tagging

(c) QuickDraw

Figure 5: LUT utilization as a function of total width for the (a) top quark tagging, (b) jet flavor tagging (c) QuickDraw model. Performance of both GRU (dashed) and LSTM (solid) models are shown. LUT utilization with latency strategy is shown only for the Top Tagging models and all other lines correspond to different reuse factors. The LUTs available in the target FPGA for each model are shown in the black dashed horizontal line. Table 2:  Minimum and maximum latencies for the top quark tagging model. Model | Latency [$\mu$s] |  $R=(6,5)$ [$\mu$s] |  $R=(12,10)$ [$\mu$s] |  $R=(30,20)$ [$\mu$s] |  $R=(60,60~{}[40])$ [$\mu$s]
---|---|---|---|---|---
GRU | 1.7–1.7 | 2.4–6.5 | 3.2–7.3 | 5.0–9.1 | 8.0–12.1
LSTM | 1.7–1.7 | 2.7–6.8 | 3.5–7.6 | 5.3–9.4 | 8.3–12.4
Table 3:  Minimum and maximum latencies for the jet flavor tagging model. Model |  $R=(48,40)$ [$\mu$s] |  $R=(90,60)$ [$\mu$s] |  $R=(120,120)$ [$\mu$s] |  $R=(240,240)$ [$\mu$s]
---|---|---|---|---
GRU | 6.7–24.8 | 9.8–27.9 | 11.5–29.6 | 20.5–38.6
LSTM | 6.9–25.0 | 10.1–28.2 | 11.7–29.8 | 20.7–38.8
Table 4:  Minimum and maximum latencies for the QuickDraw model. Model |  $R=(48,32)$ [$\mu$s] |  $R=(96,64)$ [$\mu$s] |  $R=(192,128)$ [$\mu$s] |  $R=(384,384~{}[256])$ [$\mu$s]
---|---|---|---|---
GRU | 35.4–164.0 | 59.4-–188.0 | 107.0–235.0 | 203.0–331.0
LSTM | 35.9–164.0 | 59.9–188.0 | 107.0–236.0 | 203.0–332.0

###  5.3 Static and Non-static Comparison

In order to study the impacts of the static and non-static modes discussed in
Sec 3 we limit our consideration to the top quark tagging models. As shown in
Fig. 6, resource usage for non-static mode increases dramatically compared to
static mode. For even moderate-sized models, non-static mode requires too many
resources to be feasible. In the case of the top quark tagging model we see
that non-static mode is able to fit within the available resources of the chip
only for very small bitwidths. However, Table 5 confirms that although non-
static mode offers similar overall latency to static mode, the II in non-
static mode is reduced from 315 (314) to 1 for the GRU (LSTM) models. This
results in a increased throughput for non-static mode by a factor of more than
300. The increased throughput of non-static mode would be vital for L1 trigger
applications that run inferences at rates of up to 40 MHz. While this
particular top quark tagging model suffers in performance using a total
bitwidth of 10 (6 integer and 4 fractional bits in this case), there are
multiple options, such as per-layer quantization or quantization-aware
training, that were not considered for this study but could potentially allow
a performant version of this model to be synthesized in non-static mode.

(a) Top quark tagging: DSP

(b) Top quark tagging: FF

(c) Top quark tagging: LUT

Figure 6: Resource usage of DSPs (a), FFs (b), and LUTs (c) for the top quark tagging model in static and non-static mode with solid and dashed lines, respectively. The resources available in the target FPGA are shown in the black dashed horizontal line. Table 5:  Minimum and maximum latencies and initiation intervals for the top quark tagging model in both static and non-static mode. Model | Static latency [$\mu$s] | Non-static latency [$\mu$s] | Static II | Non-static II
---|---|---|---|---
GRU | 1.7–1.7 | 1.6–1.6 | 315 | 1
LSTM | 1.6–1.6 | 1.5–1.5 | 314 | 1

##  6 Summary and Outlook

RNNs have shown substantial success for many tasks in particle physics. They
are particularly well-suited to those problems involving sequences of particle
or detector signals, outperforming densely connected deep neural networks
(DNNs) [43] and convolution neural networks (CNNs) [44] on certain jet
classification tasks. In spite of this success, RNNs have not seen the
widespread adoption in ultra-low latency environments in physics when compared
to DNNs and CNNs. This difference is owed in part to tools such as hls4ml that
simplify the adaptation of the latter models from Keras to HLS. The support
for GRUs and LSTMs in hls4ml that we present in this work represents the
removal of a major barrier to the use of RNNs in ultra-low latency
environments. This has ramifications not only for high energy physics but also
other research areas where RNNs have become popular. While we have focused on
the usage of hls4ml with FPGAs, it is important to note that hls4ml can also
be used to create ASIC designs [45], and thus this work also allows for the
possibility of RNN usage on ASICs as well.

The implementation we present in hls4ml in this work maintains the main tuning
features of hls4ml, namely the reuse factor and per-layer bit precision. This
is necessary to allow the customization of the synthesized design to meet the
needs of a given task. The benchmark models chosen cover a range of sizes,
latencies, and problems, and showcase the quality of the hls4ml support for a
variety of realistic scenarios. We also add an RNN-specific tuning parameter
to hls4ml called the RNN mode, with static and non-static settings capable of
further adjusting the behavior of the synthesized design. While we show that
this work is capable of producing results with high accuracy, there are
multiple possibilities for future development. In particular, we observe that
even small RNN models can require a substantial amount of resources to
implement. While the post-training quantization scheme we have used here is
able to minimize resources to a certain extent, other methods, such as
quantization-aware training, have shown that even more resource reduction can
be possible with little to no cost to performance. This is perhaps even more
true for RNNs than dense neural networks due to the repeated use of the
recurrent layer weights. Other techniques such as masking are also a possible
method for reducing both resource usage and dependence on small weight values
(high bit precision).

The recurrent or repeating nature of many modern algorithms, such as RNNs,
transformers and graph neural networks, make them very difficult to be run,
particularly at low latency, on FPGAs. In this work, we present the successful
deployment of RNNs in models with number of trainable parameters ranging from
$\mathcal{O}$(1 k) to $\mathcal{O}$(100 k) achieving latencies of
$\mathcal{O}$(1 $\mu$s) to $\mathcal{O}$(100 $\mu$s). This represents an
important step in enabling support in hls4ml for more complex architectures
with recursive computations.

##  7 Acknowledgment

We acknowledge the Fast Machine Learning collective as an open community of
multi-domain experts and collaborators. This community was important for the
development of this project. We thank Javier Duarte, Maurizio Pierini,
Zhiqiang Que and Nhan Tran for their valuable comments and suggestions. M.
Kagan and R. Teixeira de Lima are supported by the US Department of Energy
(DOE) under grant DE-AC02-76SF00515. P. Harris and D. Rankin acknowledge DOE
grant DE-SC0021943 and National Science Foundation (NSF) grants No. 1934700
and 1931469. S.-C. Hsu and E. Khoda are supported by NSF grants No.2117997 and
1934360.

##  8 Code Availability Statement

The hls4ml library is available at
https://github.com/fastmachinelearning/hls4ml and RNN support is available as
of commit 59ed8249f4bbdb4b23ff0c6f0bfc976b44d3ac7e. For examples on how to use
hls4ml, the notebooks in https://github.com/fastmachinelearning/hls4ml-
tutorial serve as a general introduction.

## References

  * [1] J. Duarte et al., “Fast inference of deep neural networks in FPGAs for particle physics”,  JINST 13 (2018) P07027, doi:10.1088/1748-0221/13/07/P07027, arXiv:1804.06913.
  * [2] J. Ngadiuba et al., “Compressing deep neural networks on FPGAs to binary and ternary precision with hls4ml”,  Mach. Learn.: Sci. Technol. 2 (2020) 015001, doi:10.1088/2632-2153/aba042, arXiv:2003.06308.
  * [3] T. Åarrestad et al., “Fast convolutional neural networks on FPGAs with hls4ml”,  Mach. Learn.: Sci. Technol. 2 (7, 2021) 045015, doi:10.1088/2632-2153/ac0ea1, arXiv:2101.05108.
  * [4] D. S. Rankin et al., “FPGAs-as-a-Service Toolkit (FaaST)”, in  2020 IEEE/ACM International Workshop on Heterogeneous High-performance Reconfigurable Computing (H2RC).  10, 2020.  arXiv:2010.08556.  doi:10.1109/H2RC51942.2020.00010.
  * [5] S. Hochreiter and J. Schmidhuber, “Long Short-Term Memory”,  Neural Computation 9 (11, 1997) 1735–1780, doi:10.1162/neco.1997.9.8.1735, arXiv:https://direct.mit.edu/neco/article-pdf/9/8/1735/813796/neco.1997.9.8.1735.pdf.
  * [6] K. Cho et al., “Learning phrase representations using RNN encoder-decoder for statistical machine translation”, CoRR abs/1406.1078 (2014) arXiv:1406.1078.
  * [7] R. Pascanu, T. Mikolov, and Y. Bengio, “On the difficulty of training recurrent neural networks”, 2012.  doi:10.48550/ARXIV.1211.5063, https://arxiv.org/abs/1211.5063.
  * [8] ATLAS Collaboration, “Identification of Jets Containing $b$-Hadrons with Recurrent Neural Networks at the ATLAS Experiment”,.
  * [9] ATLAS Collaboration, “Identification of hadronic tau lepton decays using neural networks in the ATLAS experiment”,.
  * [10] R. Teixeira de Lima, “Sequence-based Machine Learning Models in Jet Physics”, arXiv:2102.06128.
  * [11] K. Goto et al., “Development of a Vertex Finding Algorithm using Recurrent Neural Network”, arXiv:2101.11906.
  * [12] M. Wielgosz, A. Skoczeń, and M. Mertik, “Using lstm recurrent neural networks for monitoring the lhc superconducting magnets”,  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 867 (2017) 40–50, doi:https://doi.org/10.1016/j.nima.2017.06.020.
  * [13] A. Schmitt, K. Fu, S. Fan, and Y. Luo, “Investigating deep neural networks for gravitational wave detection in advanced ligo data”, in  Proceedings of the 2nd International Conference on Computer Science and Software Engineering, p. 73–78.  Association for Computing Machinery, New York, NY, USA, 2019.  doi:10.1145/3339363.3339377.
  * [14] A. Li et al., “KamNet: An Integrated Spatiotemporal Deep Neural Network for Rare Event Search in KamLAND-Zen”, arXiv:2203.01870.
  * [15] E. Flurin, L. S. Martin, S. Hacohen-Gourgy, and I. Siddiqi, “Using a recurrent neural network to reconstruct quantum dynamics of a superconducting qubit from physical observations”,  Phys. Rev. X 10 (Jan, 2020) 011006, doi:10.1103/PhysRevX.10.011006.
  * [16] “CERN Accelerating science. (n.d.)”. https://home.cern/science/accelerators/large-hadron-collider, 2016.
  * [17] J. Duarte et al., “Fast inference of deep neural networks in FPGAs for particle physics”,  JINST 13 (2018), no. 07, P07027, doi:10.1088/1748-0221/13/07/P07027, arXiv:1804.06913.
  * [18] Y. Umuroglu et al., “FINN”, in  Proceedings of the 2017 ACM/SIGDA International Symposium on Field-Programmable Gate Arrays.  ACM, feb, 2017.  doi:10.1145/3020078.3021744.
  * [19] “hls4ml: keras-lstm”, June, 2020.  https://github.com/richarao/hls4ml/tree/keras-lstm.
  * [20] F. Chollet et al., “Keras”, 2015\.  https://keras.io.
  * [21] C. Heelan, A. V. Nurmikko, and W. A. Truccolo, “Fpga implementation of deep-learning recurrent neural networks with sub-millisecond real-time latency for bci-decoding of large-scale neural sensors (104 nodes)”,  2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) (2018) 1070–1073.
  * [22] A. X. M. Chang, B. Martini, and E. Culurciello, “Recurrent neural networks hardware implementation on FPGA”, CoRR abs/1511.05552 (2015) arXiv:1511.05552.
  * [23] M. Lee et al., “Fpga-based low-power speech recognition with recurrent neural networks”, CoRR abs/1610.00552 (2016) arXiv:1610.00552.
  * [24] J. Fowers et al., “A configurable cloud-scale dnn processor for real-time ai”, in  Proceedings of the 45th Annual International Symposium on Computer Architecture, ISCA ’18, p. 1–14.  IEEE Press, 2018.  doi:10.1109/ISCA.2018.00012.
  * [25] S. Han et al., “ESE: efficient speech recognition engine with compressed LSTM on FPGA”, CoRR abs/1612.00694 (2016) arXiv:1612.00694.
  * [26] G. Aad et al., “Artificial Neural Networks on FPGAs for Real-Time Energy Reconstruction of the ATLAS LAr Calorimeters”,  Comput. Softw. Big Sci. 5 (2021), no. 1, 19, doi:10.1007/s41781-021-00066-y.
  * [27] Z. Que et al., “Accelerating Recurrent Neural Networks for Gravitational Wave Experiments”, in  32nd IEEE International Conference on Application-specific Systems, Architectures and Processors.  6, 2021.  arXiv:2106.14089.  doi:10.1109/ASAP52443.2021.00025.
  * [28] V. Rybalkin et al., “Finn-l: Library extensions and design trade-off analysis for variable precision lstm networks on fpgas”, 2018.  doi:10.48550/ARXIV.1807.04093, https://arxiv.org/abs/1807.04093.
  * [29] P. Coussy and A. Morawiec, “High-level synthesis”.  Springer Netherlands, Dordrecht, 1st ed. 2008. edition, 2008.
  * [30] R. Nane et al., “A survey and evaluation of FPGA high-level synthesis tools”,  IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst. 35 (2016), no. 10, 1591, doi:10.1109/TCAD.2015.2513673.
  * [31] Mentor/Siemens, “Catapult High-Level Synthesis”. https://www.mentor.com/hls-lp/catapult-high-level-synthesis, 2020.
  * [32] AMD/Xilinx, “Vivado”. https://www.xilinx.com/support/university/vivado.html, 2022.
  * [33] T. Feist, “Vivado design suite”,  White Paper 5 (2012) 30.
  * [34] J. Alwall et al., “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations”,  Journal of High Energy Physics 2014 (2014), no. 7, 79, doi:10.1007/JHEP07(2014)079.
  * [35] R. D. Ball et al., “Parton distributions with lhc data”,  Nuclear Physics B 867 (2013), no. 2, 244–289, doi:https://doi.org/10.1016/j.nuclphysb.2012.10.003.
  * [36] T. Sjöstrand et al., “An introduction to pythia 8.2”,  Computer Physics Communications 191 (2015) 159–177, doi:https://doi.org/10.1016/j.cpc.2015.01.024.
  * [37] P. Skands, S. Carrazza, and J. Rojo, “Tuning pythia 8.1: the monash 2013 tune”,  The European Physical Journal C 74 (2014), no. 8, 3024, doi:10.1140/epjc/s10052-014-3024-y.
  * [38] V. Nair and G. E. Hinton, “Rectified linear units improve restricted boltzmann machines”, in  ICML 2010, pp. 807–814.  2010\.
  * [39] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization”, in  3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, Y. Bengio and Y. LeCun, eds.  2015\.
  * [40] “CERN Open Data Portal”. http://opendata.cern.ch.
  * [41] “QuickDraw dataset”. https://github.com/googlecreativelab/quickdraw-dataset.
  * [42] CMS Collaboration, “The Phase-2 upgrade of the CMS Level-1 trigger”, CMS Technical Design Report CERN-LHCC-2020-004. CMS-TDR-021, CERN, 2020.
  * [43] S. Egan et al., “Long Short-Term Memory (LSTM) networks with jet constituents for boosted top tagging at the LHC”, arXiv:1711.09059.
  * [44] K. Fraser and M. D. Schwartz, “Jet Charge and Machine Learning”,  JHEP 10 (2018) 093, doi:10.1007/JHEP10(2018)093, arXiv:1803.08066.
  * [45] G. Di Guglielmo et al., “A Reconfigurable Neural Network ASIC for Detector Front-End Data Compression at the HL-LHC”,  IEEE Trans. Nucl. Sci. 68 (2021), no. 8, 2179–2186, doi:10.1109/TNS.2021.3087100, arXiv:2105.01683.