\section{Introduction}

The potential of graphene~\cite{geim_naturemater_2007,castroneto_rmp_2009,peres_rmp_2010,dassarma_rmp_2011,kotov_rmp_2012,Katsnelsonbook, bonaccorso_matertoday_2012}---a two-dimensional (2D) crystal of Carbon atoms tightly packed in a honeycomb lattice---in optoelectronics, photonics, and plasmonics is attracting a truly considerable attention~\cite{bonaccorso_naturephoton_2010,koppens_nanolett_2011,grigorenko_naturephoton_2012}.

The particle-hole-symmetric spectrum of massless Dirac fermions~\cite{geim_naturemater_2007,castroneto_rmp_2009,peres_rmp_2010,dassarma_rmp_2011,kotov_rmp_2012,Katsnelsonbook} in graphene and its structureless optical conductivity~\cite{universal} are particularly suitable~\cite{bonaccorso_naturephoton_2010} to realize detectors of radiation in a wide range of photon energies from visible to terahertz (THz) frequencies.
As a consequence, the photoresponse of a graphene sheet has been the subject of truly intense investigations~\cite{xia_naturenano_2009, mueller_prb_2009,mueller_naturephoton_2010,xu_nanolett_2010,peters_apl_2010,echtermeyer_naturecommun_2011,lemme_nanolett_2011,song_nanolett_2011,gabor_science_2011,rao_acsnano_2011,konstantatos_naturenano_2012,vora_apl_2012,yan_naturenano_2012,vicarelli_naturemat_2012,freitag_naturephoton_2013,freitag_nanolett_2013}.
Three main photoresponse mechanisms have been identified in the literature to date.
(i) A finite dc response to an oscillating radiation field can simply stem from the ordinary {\it photovoltaic} effect.
Photons impinging on the sample excite electron-hole pairs and their basic constituents (the electron and the hole) are separated by the electric field across a p-n junction.
(ii) It has been recognized that hot-carrier-assisted transport plays an important role in graphene~\cite{xu_nanolett_2010}.
Due to the large optical phonon energy scale in this material, hot carriers created by the radiation field can remain at a temperature higher than that of the lattice for several tens of ${\rm ps}$.
Equilibration with the lattice occurs indeed mainly because of scattering between carriers and acoustic phonons~\cite{bistritzer_prl_2009,tse_prb_2009}.
These processes take place on a ${\rm ns}$ timescale, although they can experience a speed-up in the case of disorder-assisted collisions~\cite{song_prl_2012, graham_naturephys_2012, betz_naturephys_2012}.
Hot carriers can therefore significantly alter the photoresponse of a graphene sheet~\cite{xu_nanolett_2010,song_nanolett_2011,gabor_science_2011} by virtue of the {\it photo-thermoelectric} effect.
Carrier multiplication~\cite{brida_naturecommun_2013,tomadin_arxiv_2013} can greatly enhance the performance of graphene photodetectors operating on the basis of photovoltaic and photo-thermoelectric effects.
(iii) Finally, {\it bolometric} effects can play an important role in the photoresponse of a graphene sheet.
In this case, the radiation impinging on the sample produces heating, which, in turn, affects the magnitude of the resistance.

\begin{figure}
\includegraphics[width=\linewidth]{fig01}
\caption{\label{fig:one}
Schematics of a Dyakonov-Shur single-layer-graphene field-effect transistor for photodetection.
An oscillating radiation field coupled to a graphene field-effect transistor creates an oscillating potential difference $U(t) = U_0 + U_{\rm a}\cos(\Omega t)$ between the back gate and source.
The graphene flake is deposited on a substrate with dielectric constant $\epsilon_1$ and is separated by a distance $d$ from the back gate.
The dielectric constant of the material above the graphene sheet is $\epsilon_2$.
The separation between source and drain is $L$.
Because of the intrinsic nonlinear response of the two-dimensional electron fluid in the field-effect transistor channel, a finite dc potential difference $\Delta U$, which is proportional to the power of the incident radiation (i.e.~$\propto U^2_{\rm a}$), is measured between the source and drain at {\it zero} source-drain bias. }
\end{figure}

In a series of pioneering papers~\cite{dyakonov_prl_1993, dyakonov_prb_1995, dyakonov_ieee_1996a, dyakonov_ieee_1996b}, which appeared in the mid nineties, Dyakonov and Shur (DS) proposed a very elegant mechanism that yields a finite dc response to an oscillating radiation field.
The DS photodetection mechanism is based on the fact that a field-effect transistor (FET) hosting a 2D electron gas (2DEG) acts as a cavity for plasma waves.
When these are weakly damped, i.e.~when a plasma wave launched at the source can reach the drain in a time shorter than the momentum relaxation time $\tau$, the detection of radiation exploits constructive interference of the plasma waves in the cavity, which results in a resonantly enhanced response.
This is the so-called {\it resonant regime} of plasma-wave photodetection.
Plasma waves propagating in the FET channel cannot be simply identified with the well-known ``plasmons'' of a 2DEG~\cite{Giuliani_and_Vignale}.
Indeed, plasma waves are collective oscillations that occur in a {\it gated} 2DEG, whereby the long-range tail of the Coulomb interaction among electrons is screened by the presence of a metal gate.
Plasma waves with wave vector $k = |{\bm k}|\lesssim 1/d$ (where $d$ is the distance between the 2DEG and the gate) have a gapless dispersion relation $\omega = s k$ (at zero temperature and neglecting friction and viscosity~\cite{svintsov_jap_2012}) which resembles that of sound waves in ordinary gases and liquids.
On the contrary, plasmons in a 2DEG with long-range Coulomb interactions have a dispersion relation which is proportional to $\sqrt{k}$ in the long-wavelength limit~\cite{Giuliani_and_Vignale}.

DS showed~\cite{dyakonov_ieee_1996a} that the photovoltage response of the 2DEG in a FET, i.e.~the electric potential difference between drain and source, contains a dc component even if the incoming field is ac, and thus provides rectification of the signal.
In the resonant regime, the dc photoresponse is characterized by peaks at the odd multiples of the lowest plasma-wave frequency.
This rectification mechanism is of purely dynamical origin and is not related to other rectification mechanisms (occurring, for the example, at the contacts) which could also be present in a real device.
Note that rectification of the signal is necessary to detect incoming radiation that exceeds the typical cutoff frequencies of circuit elements.
The DS mechanism is therefore particularly useful to detect THz radiation.

The DS mechanism relies on two facts.
(i) The reflection symmetry corresponding to the exchange of source with drain in the FET channel (see Fig.~\ref{fig:one}) is broken by the DS boundary conditions.
These boundary conditions are unusual because DS fixed the value of the current at the drain and the value of the potential at the source (instead of operating the device by fixing the current or the potential both at the source and at the drain, as is more customary).
(ii) The fact that the photovoltage averaged over a cycle of the oscillating radiation field is finite ultimately stems from the nonlinearity of the continuity equation and can be viewed as a result of the gate modulating {\it both} the electron density and the drift velocity in the channel.

A substantial amount of experimental work has been carried out on DS photodetection in 2DEGs in ordinary semiconductor heterojunctions: for recent reviews see e.g.~Refs.~\onlinecite{knap_jimtw_2009,knap_nanotech_2013}.
Recent experimental work~\cite{vicarelli_naturemat_2012} has demonstrated that the DS photodetection mechanism is active also in the case of the 2D massless Dirac fermion (MDF) fluid in a graphene FET.
Vicarelli {\it et al.}~\cite{vicarelli_naturemat_2012} have demonstrated room-temperature THz detectors based on antenna-coupled graphene FETs, which exploit the DS mechanism but display also contributions of photovoltaic and photo-thermoelectric origin.
The plasma waves excited by THz radiation in Ref.~\onlinecite{vicarelli_naturemat_2012} are overdamped and the fabricated THz detectors, although enabling large area, fast imaging of macroscopic samples, do not {\it yet} operate in the aforementioned resonant regime.

In this Article we present a theory of DS plasma-wave photodetection in a graphene FET in the resonant regime.
We take into  account fundamental differences  between ordinary parabolic-band 2DEGs~\cite{dyakonov_prl_1993,dyakonov_prb_1995,dyakonov_ieee_1996a, dyakonov_ieee_1996b} and  2D MDF  fluids.
By employing the theory of relativistic hydrodynamics, we demonstrate in a fully analytical fashion that nonlinearities of purely relativistic origin substantially influence the response of graphene-based plasma-wave photodetectors.

This Article is organized as follows.
In Sec.~\ref{sect:model} we present a hydrodynamic theory of transport for MDFs in a gated graphene sheet.
In Sec.~\ref{sect:instability} we discuss the plasma-wave instability which arises when a graphene FET is subject to a small dc current bias.
In Sec.~\ref{sect:photoresponse} we use the hydrodynamic theory outlined in Sec.~\ref{sect:model} to calculate the photoresponse of a graphene FET to an oscillating electromagnetic field when DS boundary conditions are applied.
Finally, in Sec.~\ref{sect:conclusions} we summarize our main findings and draw our main conclusions.

\section{Hydrodynamic theory for massless Dirac fermions in a gated graphene sheet}
\label{sect:model}

The theory of hydrodynamics~\cite{Landau06} can be applied to describe transport in a 2D electron fluid when electron-electron (e-e) collisions take place on a time scale $\tau_{\rm ee}$ which is much shorter than the typical time scale of the evolution of macroscopic variables (density, current, and energy).
At the same time, collisions with impurities and phonons, which spoil the conservation of momentum and energy, must be assumed to occur much less frequently.

In the case of present interest, we require $\omega_{\rm P} \tau_{\rm ee} \ll 1$, where $\omega_{\rm P}$ is the lowest frequency of a plasma wave in the FET.
This requirement implies that plasma waves are collective modes which can be excited in the ``collisional'' regime, when e-e interactions dominate the dynamics.
On the contrary, plasmons are collective modes which can be excited in the ``collisionless'' regime, because they can be understood, in the first instance and at long wavelengths, as the sloshing mode of the center of mass of the electron liquid.

To estimate the order of magnitude of $\omega_{\rm P} \tau_{\rm ee}$ for plasma waves, we use the following result~\cite{svintsov_jap_2012,principi_ssc_2011} for the plasma-wave group velocity
\begin{equation}
s  = v_{\rm F} \sqrt{N_{\rm f} \alpha_{\rm ee} d k_{\rm F}}~,
\end{equation}
which will be re-derived below---see Eq.~(\ref{eq:plasmaspeed}).
Here $v_{\rm F} \sim 1~{\rm nm}/{\rm fs}$ is the Fermi velocity of MDFs in graphene, $d$ is the distance between graphene and the gate (see Fig.~\ref{fig:one}), $k_{\rm F} = \sqrt{4\pi {\bar n}/N_{\rm f}}$ is the Fermi momentum corresponding to an average carrier (electron or hole) concentration ${\bar n}$, and $N_{\rm f} = 4$ is the number of fermion flavors in graphene (due to spin and valley degrees of freedom).
Finally, $\alpha_{\rm ee}$ is a dimensionless parameter that controls the strength of e-e interactions~\cite{kotov_rmp_2012}.
A microscopic definition of $\alpha_{\rm ee}$ will be given below after Eq.~(\ref{eq:plasmaspeed}).

We take $\alpha_{\rm ee} \sim 1$, $d \sim 100~{\rm nm}$, $k_{\rm F} \sim 0.17~{\rm nm}^{-1}$ (corresponding to a typical doping ${\bar n}\sim 10^{12}~{\rm cm}^{-2}$), and we find $s \sim 8.2~{\rm nm} / {\rm fs}$.
We notice that $s \gg v_{\rm F}$ with these parameters.
The largest wave vector $k$ of a plasma-wave resonance supported by the electron liquid in the FET channel is $\pi/(2 L)$, which corresponds to a standing wave with a node at the source and an antinode at the drain.
Taking $L \sim 5~\mu{\rm m}$ for a typical device, we have $\omega_{\rm P} = s k \sim 2.6 \times 10^{-3}~{\rm fs}^{-1}$ ($\hbar \omega_{\rm P} = 1.7~{\rm meV}$), corresponding to $\nu_{\rm P} = \omega_{\rm P} / (2\pi) \sim 410~{\rm GHz}$.
The transit time of a plasma wave in the device is $L / s \sim 0.61~{\rm ps}$.
A rough estimate of $\tau_{\rm ee}$ can be obtained by dividing the typical distance between two electrons ($\sim k^{-1}_{\rm F}$) by the Fermi velocity, i.e.~$\tau_{\rm ee} \sim (v_{\rm F} k_{\rm F})^{-1} \sim 5.9~{\rm fs}$.
This is in agreement with recent ultrafast two-color pump-probe measurements that have reported a time scale of the order of tens of femtoseconds to establish thermal equilibrium in single-layer graphene~\cite{brida_naturecommun_2013}.
For plasma waves in a graphene FET channel we therefore obtain $\omega_{\rm P} \tau_{\rm ee} \sim 1.5 \times 10^{-2}$, which shows that the electron system in a typical graphene FET is in the collisional regime.
On the contrary, to verify that, in the absence of a gate, plasmon excitations cannot be described within a collisional hydrodynamic theory, we consider the dispersion~\cite{grigorenko_naturephoton_2012} $\omega = v_{\rm F} \sqrt{2 \alpha_{\rm ee} k_{\rm F} k}$, evaluated at $k \sim 0.2~k_{\rm F}$.
Using the same parameters as above, we find $\omega \sim 0.11~{\rm fs}^{-1}$ ($\hbar \omega \sim 71~{\rm meV}$) and then $\omega \tau_{\rm ee} \sim 0.63$.

When the hydrodynamic assumption holds true, the system is in local equilibrium, i.e.~the electron distribution is a Fermi-Dirac distribution and does not evolve in time due to e-e collisions.
However, the parameters of the electron distribution change slowly in time and space and determine thermodynamic and hydrodynamic observables, which are not necessarily steady and evolve in time according to macroscopic conservation laws in the presence of suitable boundary conditions and slowly-varying external potentials.

Before concluding, we emphasize that for typical device lengths, the plasma-wave frequency $\nu_{\rm P}$ is in the THz regime.
For example, changing $L$ from $1~\mu{\rm m}$ to $10~\mu{\rm m}$, the plasma-wave frequency in a graphene FET changes from $\nu_{\rm P} \sim 2.1~{\rm THz}$ to $\nu_{\rm P} \sim 210~{\rm GHz}$.
For this reason, photodetectors based on the DS mechanism are naturally useful in the context of THz light detection.

\subsection{Continuity and Euler equations}

Hydrodynamic equations~\cite{Landau06} can be derived from the Boltzmann semiclassical equation following a standard procedure~\cite{KadanoffBaym, Landau05}, which relies on the conservation of the particle number, total momentum, and total energy in e-e scattering processes.
An outline of the derivation of hydrodynamic equations for 2D MDFs in graphene is reported in the Appendix.

The continuity equation takes the familiar form
\begin{equation}\label{eq:continuity}
\partial_{t}n({\bm r},t) + \nabla_{\bm r} \cdot \lbrack n({\bm r},t) {\bm v}({\bm r},t) \rbrack  = 0~,
\end{equation}
where $n({\bm r}) = n_{\rm e}({\bm r}) + n_{\rm h}({\bm r})$ is the total carrier density and $n_{\rm e,h}({\bm r},t)$ is the electron (hole) density.
Carriers drift with the local average velocity ${\bm v}({\bm r},t)$.
The continuity equation (\ref{eq:continuity}) is identical in form to the corresponding equation for an ordinary parabolic-band 2DEG.

The Euler and energy equations, on the contrary, are dramatically different from the 2DEG case and are reported in the Appendix.
With the aim of investigating the dynamics of plasma waves in gated graphene flakes, we restrict here to the limit in which the drift velocity ${\bm v}({\bm r},t)$ is much smaller than the Fermi velocity $v_{\rm F}$.
In this case the Euler equation reads
\begin{equation}\label{eq:euler}
\begin{split}
\frac{3 P({\bm r},t)}{v_{\rm F}^{2}} \lbrace  \partial_{t} {\bm v}({\bm r},t) & + \lbrack {\bm v}({\bm r},t) \cdot \nabla_{\bm r} \rbrack {\bm v}({\bm r},t) \rbrace  = \\
& - \lbrack n_{\rm e}({\bm r},t) - n_{\rm h}({\bm r},t) \rbrack  \nabla_{\bm r} U_{\rm eff}({\bm r},t) \\
& - \nabla_{\bm r} P({\bm r},t) - \frac{{\bm v}({\bm r},t)}{v_{\rm F}^{2}} \partial_{t} P({\bm r},t)~.
\end{split}\end{equation}
Here, $P({\bm r}, t)$ is the pressure and $U_{\rm eff}({\bm r},t)$ is the effective potential which acts on electrons, including the electrostatic potential generated by the plasma wave and the fields of nearby conductors.

Note that the ``free'' parts of the Euler and energy equations [Eqs.~(\ref{eq:eulerfullrel}) and~(\ref{eq:energyfullrel}) in the Appendix, respectively] can be written in a very compact form by making use of a covariant notation, commonly employed in relativistic problems~\cite{Landau06}, where the Fermi velocity $v_{\rm F}$ plays the role of the speed of light $c$.
In the covariant notation, the Euler and energy equations correspond to the space- and time-like component of the divergence of the relativistic energy-momentum tensor~\cite{Landau06}.
This formal analogy stems from the highly-nonlinear relation [see Eq.~(\ref{eq:momentumrudin}) in the Appendix] between the average drift momentum ${\bm p}({\bm r}, t)$ and the drift velocity ${\bm v}({\bm r}, t)$, which, in turn, stems from the linear MDF dispersion of carriers in graphene~\cite{geim_naturemater_2007,castroneto_rmp_2009,peres_rmp_2010,dassarma_rmp_2011,kotov_rmp_2012,Katsnelsonbook}.
The average drift velocity ${\bm v}({\bm r},t)$ is smaller than the Fermi velocity due to the large number of e-e scattering events that take place in the hydrodynamic limit in a period of a plasma wave.

Of course, the formal analogy with relativistic fluid dynamics does not mean that the Fermi velocity is the largest velocity at which signals can causally propagate in a 2D MDF fluid.
Indeed, as shown above, the plasma-wave group velocity $s$ is typically larger than $v_{\rm F}$.
Furthermore, the Coulomb interaction between electrons and all the external potentials propagate instantaneously and do not obey relativistic causality.
The hydrodynamic theory of MDFs in graphene is therefore fundamentally different from the hydrodynamic theory of truly relativistic fluids, which is routinely employed to describe e.g.~plasmas of astrophysical interest.

We now focus our analysis on the limit $n_{\rm e}({\bm r},t) \gg n_{\rm h}({\bm r},t)$ (limit of vanishingly small hole concentration) and approximate $n({\bm r},t) \simeq n_{\rm e}({\bm r},t)$.
The local density determines the local chemical potential $\mu({\bm r},t)$ and the local Fermi energy $\varepsilon_{\rm F}({\bm r},t)$,
\begin{equation}\label{eq:fermienergy}
\varepsilon_{\rm F}({\bm r},t) = \hbar v_{\rm F} \sqrt{4 \pi n({\bm r},t) / N_{\rm f}}~.
\end{equation}
We assume that the temperature $T$ is uniform and constant in the sample and that $k_{\rm B} T \ll \varepsilon_{\rm F}({\bm r},t)$.
The pressure [see Eq.~(\ref{eq:pressurefull}) in the Appendix] can then be approximated as
\begin{equation}\label{eq:pressure}
P({\bm r},t) \simeq \frac{1}{3} \varepsilon_{\rm F}({\bm r},t) n({\bm r},t)~.
\end{equation}

The Euler equation reduces to
\begin{equation}\label{eq:Euler_complete}
\begin{split}
\partial_{t} {\bm v}({\bm r},t) & + \lbrack {\bm v}({\bm r},t) \cdot \nabla_{\bm r} \rbrack {\bm v}({\bm r},t) = - \frac{v_{\rm F}^{2}}{\varepsilon_{\rm F}({\bm r},t)} \nabla_{\bm r} U_{\rm eff}({\bm r},t) \\
& - \frac{v_{\rm F}^{2}}{2 n({\bm r},t)} \nabla_{\bm r} n({\bm r},t)
- \frac{{\bm v}({\bm r},t)}{2 n({\bm r},t)} \partial_{t} n({\bm r},t)~.
\end{split}
\end{equation}

\subsection{Electrostatics of a graphene FET channel and the ``gradual channel'' approximation}
\label{ssec:geometry}

Let us now consider a graphene flake, which lies in the plane $z=0$ and is subject to the presence of a gate located at $z = -d$.
The gate will be modeled as a perfect conductor.
According to Fig.~\ref{fig:one}, the dielectric constant $\epsilon(z)$ is equal to $\epsilon_{1}$ ($\epsilon_{2}$) for $-d < z <0$ ($z > 0$).

The effective potential is given by
\begin{equation}\label{eq:effpot}
U_{\rm eff}({\bm r},t) = -e \varphi({\bm r}',t)|_{z = 0}~,
\end{equation}
where $\varphi({\bm r}',t)$ is the three-dimensional (3D) electric potential, ${\bm r}' = ({\bm r},z)$ being a coordinate in 3D space.

At this stage it is convenient to introduce the so-called ``gate-to-channel swing''~\cite{dyakonov_prl_1993, dyakonov_prb_1995, dyakonov_ieee_1996a, dyakonov_ieee_1996b},
\begin{equation}
U({\bm r},t) = U_{0} - \varphi({\bm r}',t)|_{z = 0}~,
\end{equation}
where $U_{0} = \varphi({\bm r}',t)\vert_{z = -d}$ is the uniform electric potential at the gate.
From this definition, $\nabla_{\bm r} U_{\rm eff}({\bm r},t) = e \nabla_{\bm r} U({\bm r},t)$ and ${\bm E}({\bm r},t) = \nabla_{\bm r} U({\bm r},t)$ is the electric field in the graphene sheet.

The relation between the electric potential and the electron density is found by solving the Poisson equation $\nabla_{\bm r'} \cdot \lbrack \epsilon(z) {\bm E}({\bm r},t) \rbrack = -4 \pi e n({\bm r},t) \delta(z)$.
The average density is $\bar{n} = C U_{0} / e$, where $C = \epsilon_{1} / (4 \pi d)$ is the geometrical capacitance per unit area.
Electric potential and density fluctuations are related by
\begin{equation}\label{eq:channelpotential}
\varphi({\bm k},t)\vert_{z = 0} = - \frac{2 \pi e}{k}\frac{1 - \exp{(-2 d k)}}{\displaystyle \frac{\epsilon_{1} + \epsilon_{2}}{2} - \frac{\epsilon_{2} - \epsilon_{1}}{2} \exp(-2 d k)}n({\bm k},t)~,
\end{equation}
where ${\bm k}$ is a 2D wave vector in the graphene sheet and $\varphi({\bm k},t)\vert_{z = 0}$ and $n({\bm k},t)$ are the 2D Fourier transforms of the electric potential and of the electron density, respectively.
In the limit $d k \to 0$, in which the distance between the graphene sheet and gate is much smaller than the typical wavelength of density oscillations, the relation (\ref{eq:channelpotential}) between the density and gate-to-channel swing assumes the local form
\begin{equation}\label{eq:GCA}
n({\bm r},t) = \frac{C}{e} U({\bm r},t)~.
\end{equation}
The latter expression is known as the ``gradual channel'' approximation~\cite{dyakonov_prl_1993, dyakonov_prb_1995, dyakonov_ieee_1996a, dyakonov_ieee_1996b}.
In this regime, the long-range tail of the Coulomb interaction is completely screened by the charges on the gate.
Using the gradual channel approximation, the Euler and continuity equations read
\begin{equation}\label{eq:eulerGCA}
\begin{split}
\partial_{t} {\bm v}({\bm r},t) & + \lbrack {\bm v}({\bm r},t) \cdot \nabla_{\bm r} \rbrack {\bm v}({\bm r},t) = - \frac{v_{\rm F}^{2} e}{\varepsilon_{\rm F}({\bm r},t)} {\bm E}({\bm r},t) \\
& - \frac{v_{\rm F}^{2}}{2 U({\bm r},t)} \nabla_{\bm r} U({\bm r},t)
- \frac{{\bm v}({\bm r},t)}{2 U({\bm r},t)} \partial_{t} U({\bm r},t)
\end{split}
\end{equation}
and
\begin{equation}\label{eq:continuityGCA}
\partial_{t}U({\bm r},t) + \nabla_{\bm r} \cdot \lbrack U({\bm r},t) {\bm v}({\bm r},t) \rbrack  = 0~.
\end{equation}

\subsection{Making contact with the hydrodynamics of 2DEGs}
\label{ssec:2DEGlike}

In this Section we make a crude approximation on the Euler equation (\ref{eq:eulerGCA}) for a 2D MDF fluid to recover earlier results pertaining to the 2DEG literature.

When the two terms on the second line of Eq.~(\ref{eq:eulerGCA}) are neglected and the Fermi energy in the denominator of the first term on the right-hand side of Eq.~(\ref{eq:eulerGCA}) is assumed to be constant in space, $\varepsilon_{\rm F}({\bm r},t) \to \hbar v_{\rm F} k_{\rm F}$, where $k_{\rm F} = \sqrt{4 \pi \bar{n}/N_{\rm f}}$ is the Fermi wave number corresponding to the average electron density $\bar{n}$, Eq.~(\ref{eq:Euler_complete}) reduces to
\begin{equation}\label{eq:euler2DEG}
\partial_{t} {\bm v}({\bm r},t) + \lbrack {\bm v}({\bm r},t) \cdot \nabla_{\bm r} \rbrack {\bm v}({\bm r},t) = - \frac{e}{m_{\rm c}} \nabla_{\bm r} U({\bm r},t)~,
\end{equation}
which is identical to the Euler equation used by DS~\cite{dyakonov_prl_1993, dyakonov_prb_1995, dyakonov_ieee_1996a, dyakonov_ieee_1996b} with the bare electron mass $m$ replaced by the density-dependent cyclotron mass~\cite{geim_naturemater_2007,castroneto_rmp_2009} $m_{\rm c} = \hbar k_{\rm F}/v_{\rm F}$.
The validity of this 2DEG-like approximation to describe the plasma-wave photoresponse of a graphene FET will be critically examined below.

In the following, we assume that our system is translationally invariant in one direction (the one perpendicular to the source-drain direction) and therefore specify Eqs.~(\ref{eq:eulerGCA})-(\ref{eq:continuityGCA}) to a one-dimensional (1D) geometry in the $x$ direction, with $0\leq x \leq L$.

Linearizing Eqs.~(\ref{eq:continuityGCA})-(\ref{eq:euler2DEG}) we obtain
\begin{equation}
\partial_{t}U(x,t) + U_{0} \partial_{x}v(x,t)  = 0
\end{equation}
and
\begin{equation}
\partial_{t}v(x,t) + \frac{e}{m_{\rm c}} \partial_{x} U(x,t)  = 0~.
\end{equation}
It is well known that these two coupled equations admit solutions in the form of traveling waves
\begin{equation}
v(x,t), \, U(x,t) \propto e^{-i \omega t} e^{i k x} + \mbox{c.c.}~,
\end{equation}
with $\omega = s k$.
The plasma-wave velocity is given by~\cite{principi_ssc_2011,svintsov_jap_2012}
\begin{equation}\label{eq:plasmaspeed}
s \equiv \sqrt{\frac{e U_{\rm 0}}{m_{\rm c}}} = v_{\rm F} \sqrt{N_{\rm f} \alpha_{\rm ee} d k_{\rm F}}~,
\end{equation}
where
\begin{equation}\label{eq:alphaeegate}
\alpha_{\rm ee} \equiv \frac{e^{2}}{\hbar v_{\rm F} \epsilon_{1}}~
\end{equation}
is a dimensionless coupling constant that controls the strength of e-e interactions~\cite{footnote1}.
As shown above, for typical carrier densities and device sizes, $s \gg v_{\rm F}$.
Our derivation, however, does not formally rely on this inequality and applies also to the regime in which the Fermi velocity and the plasma wave  velocity are comparable, a situation which is achieved when $d k_{\rm F} \sim 1$.

\subsection{Expansion in powers of density fluctuations}

The Euler equation (\ref{eq:eulerGCA}) is not amenable to analytic treatment.
We therefore introduce an expansion of this equation in powers of the deviation $\delta n(x,t) \equiv n(x,t) -{\bar n}$ of the local density $n(x,t)$ from its average value $\bar{n}$.
Because of the gradual channel approximation (\ref{eq:GCA}), a similar expansion can be carried out in powers of $\delta U(x,t) \equiv U(x,t) - U_0$.

At this stage, it is also convenient to introduce {\it dimensionless} quantities by scaling lengths, electrical potentials, velocities, and time with $L$, $U_{0}$, $s$, and $L/s$, respectively.
For the sake of notational simplicity, dimensionless variables will be denoted with the same symbols used for dimensionful variables.

Expanding the Euler equation (\ref{eq:eulerGCA}) up to {\it second order} in $\delta U(x,t)$ we obtain
\begin{equation}\label{eq:eulerMDF}
\begin{split}
\partial_{t} v(x,t) & + v(x,t) \partial_{x} v(x,t) = - \partial_{x} \delta U(x,t)  \\
& + \frac{1}{2} \delta U(x,t)  \partial_{x} \delta U(x,t)\\
& - \frac{1}{2} v_{\rm F}^{2} \left \lbrack  \partial_{x} + \frac{v(x,t)}{v_{\rm F}^{2}} \partial_{t} \right  \rbrack \delta U(x,t)  \\
& + \frac{1}{2} v_{\rm F}^{2} \delta U(x,t)  \left \lbrack  \partial_{x} + \frac{v(x,t)}{v_{\rm F}^{2}} \partial_{t} \right  \rbrack  \delta U(x,t) ~.
\end{split}
\end{equation}
Once again, all the quantities appearing in the previous equation ($x$, $t$, $\delta U$, etc.) are dimensionless.
The first two terms in the right-hand side of Eq.~(\ref{eq:eulerMDF}) stem from the expansion in powers of $\delta U(x,t)$ of the term in Eq.~(\ref{eq:eulerGCA}) which is proportional to ${\bm E}({\bm r},t)/\varepsilon_{\rm F}({\bm r},t) = \nabla_{\bm r}U({\bm r},t)/\varepsilon_{\rm F}({\bm r},t)$.
The third and the fourth term originate from $\nabla_{\bm r} P({\bm r},t)$ and $\partial_t P({\bm r},t)$.
The first line in Eq.~(\ref{eq:eulerMDF}) reproduces the 1D version of Eq.~(\ref{eq:euler2DEG}).
The terms contained in the last three lines of Eq.~(\ref{eq:eulerMDF}) are qualitatively new terms pertaining to the hydrodynamic theory of the 2D MDF fluid.
They do not have an equivalent in a 2DEG.
The implications of these nonlinear terms on the plasma-wave photoresponse of a gated graphene sheet will be the subject of the reminder of this Article.

We now proceed to estimate the relative magnitude of each term in the right-hand side of Eq.~(\ref{eq:eulerMDF}).
We assume a profile $\delta U(x,t)  = U(k,\omega)\cos(k x - \omega t)$ with $\omega = k$.
In this case, the differential operators simplify to $\partial_{x} \to k$ and $\partial_{t} \to \omega$, and $k$ can be collected as a common prefactor.
We remind the reader that $v(x,t) \ll v_{\rm F}$.
As discussed above, typically $v_{\rm F} \ll 1$ (the Fermi velocity is much smaller than the plasma-wave group velocity when $d k_{\rm F}\gg 1$).
The first nonlinear term on the right-hand side of Eq.~(\ref{eq:eulerMDF}) can be estimated as
\begin{equation}\label{eq:fluctuationsmass}
\frac{1}{2} \delta U(x,t)  \partial_{x} \delta U(x,t) \approx \frac{k}{2}[U(k,\omega)]^2~.
\end{equation}
Similarly,
\begin{equation}
\begin{split}
\frac{1}{2} v_{\rm F}^{2} \left\lbrack \partial_{x} \vphantom{\frac{v(x,t)}{v_{\rm F}^{2}}} \right. & + \left.\frac{v(x,t)}{v_{\rm F}^{2}} \partial_{t} \right\rbrack \delta U(x,t) \\
& \approx \frac{k}{2}\left[v^2_{\rm F} + v(x,t)\right]U(k,\omega)~.
\end{split}
\end{equation}
We therefore see that the term (\ref{eq:fluctuationsmass}) does not contain velocity factors.
Keeping only this term in Eq.~(\ref{eq:eulerMDF}) can be seen as a ``zeroth-order expansion'' in powers of both $v(x,t) \ll 1$ and $v_{\rm F} \ll 1$.

\section{Generation of plasma waves}
\label{sect:instability}

\subsection{Plasma waves in a 2D MDF fluid moving with a constant speed}
\label{ssec:linearization}

In this Section we study the dispersion of plasma waves in a 2D MDF fluid that moves with an average constant speed $v_0$, measured in units of the plasma-wave speed $s$ defined in Eq.~(\ref{eq:plasmaspeed}).
We repeat the analysis of Sect.~\ref{ssec:2DEGlike}, allowing, however, for $v_0 \neq 0$.

To this aim, we discard all the terms in Eq.~(\ref{eq:eulerMDF}) which are ${\cal O}(\delta U^2)$ obtaining
\begin{equation}\label{eq:eulerLin}
\begin{split}
\partial_{t} v(x,t) & + v(x,t) \partial_{x} v(x,t) = - \partial_{x} \delta U(x,t) \\
& - \frac{1}{2} v_{\rm F}^{2} \left \lbrack  \partial_{x} + \frac{v(x,t)}{v_{\rm F}^{2}} \partial_{t} \right  \rbrack  \delta U(x,t)~.
\end{split}
\end{equation}
As usual, the previous equation must be supplemented by the continuity equation.
We now seek solutions of the form $\delta U(x,t) = U_1(x,t)$ and $v(x,t) = v_0 + v_1(x,t)$, with
\begin{equation}\label{eq:travelingwaves}
v_1(x,t), \, U_1(x,t) \propto e^{-i \omega t} e^{i k x} + \mbox{c.c.}~.
\end{equation}
We find that for each value of $\omega$, two values of $k$ are possible, corresponding to waves traveling in opposite directions.
In the limit $v_0 \ll v_{\rm F}$, the dispersion reads
\begin{equation}\label{eq:relvelsum}
\omega \simeq \left \lbrack \frac{3}{4}v_{0} \pm s(v_{\rm F}) \right \rbrack k~,
\end{equation}
where
\begin{equation}\label{eq:relsoundspeed}
s(v_{\rm F}) \equiv \sqrt{1 + v_{\rm F}^{2}/2}~.
\end{equation}
The previous result differs from Eq.~(\ref{eq:plasmaspeed}) because we have taken into account the contribution of the pressure to the restoring force, i.e.~the term $-v^2_{\rm F}\partial_x \delta U(x,t)/2$ in Eq.~(\ref{eq:eulerLin}).
The correction $v_{\rm F}^{2} / 2$ in Eq.~(\ref{eq:relsoundspeed}) can therefore be understood from the following general argument.
From Eqs.~(\ref{eq:fermienergy}) and (\ref{eq:pressure}) we see that the 2D MDF fluid fulfills a {\it polytropic} equation of state, $P = C \rho^{\gamma}$, $\rho = m_{\rm c} {\bar n}$ being the mass density and $\gamma = 3/2$.
Then, the usual thermodynamic expression for the sound speed gives $s_{\rm P} = \gamma P /\rho = v_{\rm F}/\sqrt{2}$.

Before concluding, we would like to comment on the factor $3/4$ which appears in front of $v_0$ in Eq.~(\ref{eq:relvelsum}).
From Eq.~(\ref{eq:relvelsum}) with $v_{0} = 0$, we see that $s(v_{\rm F})$ is the plasma-wave speed in a reference system which moves with velocity $v_{0}$ with respect to the laboratory rest frame.
Contrary to the 2DEG case, however, the plasma-wave speed in the laboratory rest frame cannot simply be obtained by performing a Galilean transformation $v \to v_{0} + v$.
In the case of present interest, we find that there is a correction of magnitude $v_0/4$ which reduces (increases) the plasma-wave speed, with respect to the result of the Galilean transformation $v_{0} \pm s(v_{\rm F})$.
This fact is well-known in the context of relativistic hydrodynamics~\cite{font_aa_1994}.

For future purposes, we introduce the parameter
\begin{equation}\label{eq:xi}
\xi \equiv \frac{v_{\rm F}}{s(v_{\rm F})} = \frac{1}{\sqrt{N_{\rm f}\alpha_{\rm ee} k_{\rm F} d +1/2}}~.
\end{equation}

\subsection{Plasma-wave instability}

The dispersion relation (\ref{eq:relvelsum}) characterizes the collective modes of a 2D MDF fluid in the presence of a metal gate and in the thermodynamic limit.

DS showed~\cite{dyakonov_prl_1993} that appropriate boundary conditions in a device of length $L$ define a resonator which supports discrete modes with frequencies $\lbrace \omega_{n} \rbrace_{n \in {\mathbb N}}$.
Moreover, when the average drift velocity $v_{0}$ varies in a certain interval, the resonating modes become unstable and their amplitude increases exponentially with a rate $\Gamma = \Im m \lbrack \omega_{n} \rbrack$, which is independent of $n$.

The DS boundary conditions for the instability are
\begin{equation}\label{eq:bcinst}
U(x = 0,t) = 1, \quad j(x=1,t) = v_{0}~,
\end{equation}
where $j(x,t) = U(x,t)v(x,t)$ is the current density.
The boundary conditions correspond to a constant density at the source ($x = 0$) and a constant current at the drain ($x = 1$).
We study the DS instability with both space and time derivatives of the pressure taken into account.
However, we neglect momentum relaxation, which takes place on a time scale $\tau$.
This approximation is consistent if we assume that the exponential growth of the amplitude of the modes is quenched by nonlinear terms in the equations of motion on a time scale shorter than $\tau$.

The boundary conditions for the fluctuations, which follow from Eq.~(\ref{eq:bcinst}), read
\begin{equation}\label{eq:bclininst}
U_{1}(x = 0,t) = 0, \quad v_{1}(x=1,t) + v_{0} U_{1}(x=1,t) = 0~.
\end{equation}
To find the allowed values of $\omega$, we insert the ansatz~(\ref{eq:travelingwaves}) into Eq.~(\ref{eq:bclininst}) and solve for $\omega$ to linear order in $v_{0}$.
We find that the real part of the frequency is of the form
 \begin{equation}\label{spagnasuca}
\Re e \lbrack \omega_{n} \rbrack = s(v_{\rm F}) \frac{\pi}{2} n~,
\end{equation}
where $n$ is a positive integer, while the imaginary part is
$\Gamma = \Im m \lbrack \omega_{n} \rbrack = 3 v_{0} / 4$.
The fundamental $n=1$ plasma resonance defines the Fermi-velocity-dependent plasma frequency
\begin{equation}\label{eq:omegaPvF}
\omega_{\rm P}(v_{\rm F}) \equiv \frac{\pi}{2} \sqrt{1 + \frac{v_{\rm F}^{2}}{2}}~.
\end{equation}
With respect to the ordinary 2DEG~\cite{dyakonov_prl_1993}, the instability rate $\Gamma$ is reduced by a factor $3/4$.
As noted above, this effect is of purely relativistic origin and stems from the non-Galilean nature of electron dynamics in graphene.

\section{Resonant plasma-wave photoresponse}
\label{sect:photoresponse}

\begin{figure}
\begin{overpic}[width=\linewidth]{fig02}\put(2,60){}\end{overpic}
\caption{\label{fig:photovoltage}
(Color online) The dc photovoltage $\Delta U$ (measured in units of the incident power) as a function of the frequency $\Omega$ of the incoming radiation [measured in units of the plasma-wave frequency (\ref{eq:omegaPvF})]. The solid line represents the dc photovoltage calculated within the 2DEG-like approximation discussed in Sect.~\ref{ssec:2DEGlike} [$C_{1}$,$C_{2}\to 0$ in Eq.~(\ref{eq:photovoltage})]. The dashed line is the photovoltage obtained by including the contribution of the local fluctuations of the cyclotron mass [$C_{1}\to 0$ and $C_2 \neq 0$ in Eq.~(\ref{eq:photovoltage})]. Finally, the dash-dotted line is the full expression in Eq.~(\ref{eq:photovoltage}) ($C_1, C_2 \neq 0$).
In the last two cases, $v_{\rm F}$ is obtained from Eq.~(\ref{eq:plasmaspeed}) with $\alpha_{\rm ee} \times d k_{\rm F} = 0.18$.
We remind the reader that the ratio reported on the horizontal axis depends on the magnitude of $v_{\rm F}$.
In all three cases, $\tau = 1.0~L/s$. }
\end{figure}

We now consider the response of the system to a periodic modulation of the gate-to-channel swing at the source, with vanishing current at the drain.
The periodic modulation at the source can be induced by an impinging radiation (see Fig.~\ref{fig:one}), collected, for instance, by means of an appropriate antenna~\cite{vicarelli_naturemat_2012}.

We want to solve Eq.~(\ref{eq:eulerMDF}) together with the continuity equation
\begin{equation}\label{eq:continuity1D}
\partial_{t}U(x,t) + \partial_x\lbrack U(x,t)v(x,t)\rbrack  = 0~,
\end{equation}
with the DS boundary conditions for detection:
\begin{equation}\label{eq:bcphotoresp}
U(x = 0,t) = 1 + U_{\rm a}\cos(\Omega t), \quad j(x=1,t) = 0~.
\end{equation}
Here $\Omega$ is the frequency of the incoming radiation.
Our aim is to demonstrate that a dc potential difference $\Delta U$ is generated between the drain and source in response to the oscillating radiation and {\it in the absence of a source-drain bias}.
Most importantly, we will elucidate the contributions to $\Delta U$ arising from the peculiar relativistic corrections to the Euler equation of motion pertaining to graphene.

Following DS~\cite{dyakonov_ieee_1996a}, we seek solutions of the Euler (\ref{eq:eulerMDF}) and continuity (\ref{eq:continuity1D}) equations in the form of power expansions:
\begin{equation}\label{eq:expansionv}
v(x,t) = \epsilon v_{1}(x,t) + \epsilon^{2}\lbrack \delta v(x) + v_{2}(x,t) \rbrack + \dots
\end{equation}
and
\begin{equation}\label{eq:expansionU}
U(x,t) = 1 + \epsilon U_{1}(x,t) + \epsilon^{2}\lbrack \delta U(x) + U_{2}(x,t) \rbrack + \dots~.
\end{equation}
Here, $\epsilon$ is a dimensionless parameter that helps with the bookkeeping while solving Eqs.~(\ref{eq:eulerMDF}) and~(\ref{eq:continuity1D}) perturbatively.
The functions $v_1(x,t)$, $U_1(x,t)$ and $v_2(x,t)$, $U_2(x,t)$ are assumed to oscillate in time with frequency $\Omega$ and $2\Omega$, respectively, while $\delta v(x)$ and $\delta U(x)$ are constant in time.
With the ansatz (\ref{eq:expansionU}) the dc photovoltage is $\Delta U \equiv \delta U(1) - \delta U(0)$.
This perturbative approach is fully justified when the external perturbation is sufficiently weak, i.e.~when $U_{\rm a} = {\cal O}(\epsilon)$.
We find that $\Delta U  = {\cal O}(\epsilon^2)$, i.e.~the dc photovoltage is proportional to the power of the incoming radiation.

To linear order, the equations of motion read:
\begin{equation}\label{eq:eulerlin}
\partial_{t} v_{1}(x,t) = - (1 + v_{\rm F}^{2}/2) \partial_{x} U_{1}(x,t) - v_1(x,t) / \tau
\end{equation}
and
\begin{equation} \label{eq:contlin}
\partial_{t} U_{1}(x,t) + \partial_{x} v_{1}(x,t) = 0~,
\end{equation}
with boundary conditions
\begin{equation}\label{eq:bclinphotoresp}
U_{1}(x=0,t) = U_{\rm a} \cos{(\Omega t)},\quad v_{1}(x=1,t) = 0~.
\end{equation}

Note that in Eq.~(\ref{eq:eulerlin}) we have included a linear phenomenological friction term, proportional to the inverse momentum relaxation time $\tau^{-1}$.

The solution of the linear system of differential equations (\ref{eq:eulerlin})-(\ref{eq:contlin}) is given by
\begin{equation}\label{eq:linearSolutionDS}
v_{1}(x,t) = \frac{U_{\rm a}}{2} \frac{\Omega}{K} \left \lbrack  \frac{e^{i K x}}{1+e^{2 i K}} - \frac{e^{-i K x}}{1+e^{-2 i K}} \right  \rbrack  e^{-i \Omega t} + \mbox{c.c.}
\end{equation}
and
\begin{equation}
U_{1}(x,t) = \frac{U_{\rm a}}{2} \left \lbrack  \frac{e^{i K x}}{1+e^{2 i K}} + \frac{e^{-i K x}}{1+e^{-2 i K}} \right  \rbrack  e^{-i \Omega t} + \mbox{c.c.}~,
\end{equation}
with
\begin{equation}\label{eq:dispersionK}
K = K(\Omega) = \frac{\Omega}{s(v_{\rm F})} \sqrt{1 + \frac{i}{\Omega \tau}}~.
\end{equation}
Note that in the limit $\Omega \tau \gg 1$ Eq.~(\ref{eq:dispersionK}) coincides with Eq.~(\ref{eq:relvelsum}) with $v_{0} = 0$.

\begin{figure}
\begin{overpic}[width=\linewidth]{fig03a}\put(2,60){(a)}\end{overpic}
\begin{overpic}[width=\linewidth]{fig03b}\put(2,60){(b)}\end{overpic}
\caption{\label{fig:visibility}
(Color online) (a) The visibility ${\cal V}$ of the photoresponse as a function of $\xi$ as defined in Eq.~(\ref{eq:xi}).
Different curves correspond to different values of the momentum relaxation time $\tau$ measured in units of the channel length $L$ divided by the plasma wave speed $s$: from bottom to top the value of $\tau$ increases.
The top horizontal axis reports the density ${\bar n}$ corresponding to the value of
$\xi$ reported in the bottom horizontal axis: in this case we have fixed $d = 10~{\rm nm}$ and $\alpha_{\rm ee} = 0.9$.
(b) The average carrier density ${\bar n}$ (in units of $10^{12}~{\rm cm}^{-2}$) as a function of $\xi$.
Different curves correspond to different values of $d$ (in units of ${\rm nm}$): from top to bottom $d$ increases. }
\end{figure}

We then write the Euler equation to second order in $\epsilon$ and average it over one period of the incoming radiation.
We find
\begin{equation}\label{eq:euleraverage}\begin{split}
\langle v_{1}(x,t) \partial_{x}v_{1}(x,t) \rangle_{t} & = - (1 + v_{\rm F}^{2}/2) \partial_{x} \delta U(x) - \frac{1}{\tau} \delta v(x) \\
& + \frac{1}{2}(1 + v_{\rm F}^{2}) \langle U_{1}(x,t) \partial_{x} U_{1}(x,t) \rangle_{t} \\
& - \frac{1}{2} \langle v_{1}(x,t) \partial_{t} U_{1}(x,t) \rangle_{t}~,
\end{split}\end{equation}
where $\langle f(t) \rangle_{t} \equiv T^{-1} \int_{0}^{T} dt f(t)$ denotes the time average over one period $T = 2\pi/\Omega$ of the external radiation.

Eq.~(\ref{eq:euleraverage}) depends on $\delta v(x)$---see Eq.~(\ref{eq:expansionv}).
This quantity can be easily obtained by averaging over time the continuity equation, written up to second order in $\epsilon$.
We find
$\delta v(x) = - \langle v_{1}(x,t) U_{1}(x,t) \rangle_{t}$.
The final expression for the photovoltage is
\begin{equation}\begin{split}
\Delta U  = & \frac{1}{1 + v_{\rm F}^{2} / 2} \left \lbrace \frac{1}{2} \langle v_{1}(0,t)^{2} \rangle_{t} \right . \\
& + \frac{1}{4}(1 + v_{\rm F}^{2}) \lbrack \langle U_{1}(1,t)^{2} \rangle_{t} - \langle U_{1}(0,t)^{2} \rangle_{t} \rbrack \\
& + \frac{1}{\tau} \int_{0}^{1} dx \langle v_{1}(x,t) U_{1}(x,t) \rangle_{t} \\
& \left . - \frac{1}{2} \int_{0}^{1} dx \langle v_{1}(x,t) \partial_{t} U_{1}(x,t) \rangle_{t} \vphantom{\frac{1}{2}}  \right \rbrace~.
\end{split}\end{equation}
Remarkably, we have found an analytical expression for $\Delta U$ which can be given in the following rather compact form:
\begin{equation}\label{eq:photovoltage}
\begin{split}
\Delta U & = \left ( \frac{U_{\rm a}}{2} \right )^{2} \left \lbrace 1 - \frac{C_{2}}{2} - \frac{2 - C_{2}}{\cos{2 K_{1}} + \cosh{2 K_{2}}} \right . \\
& + \left . \beta (1 - C_{1}) \frac{\cosh{2 K_{2}} - \cos{2 K_{1}}}{\cos{2 K_{1}} + \cosh{2 K_{2}}} \right \rbrace~,
\end{split}\end{equation}
where $C_{1} = 1/4$, $C_{2} = (1 + v_{\rm F}^{2}) / (1 + v_{\rm F}^{2} / 2)$, $\beta = 2 \Omega \tau / \sqrt{1 + (\Omega \tau)^{2}}$, and $K_{1}$ ($K_{2}$) is the real (imaginary) part of $K$, which reads
\begin{equation}
K_{1,2} = \frac{\Omega}{s(v_{\rm F})} \frac{1}{\sqrt{2}} \sqrt{ \sqrt{1 + \frac{1}{(\Omega \tau)^{2}}} \pm 1}~.
\end{equation}
Eq.~(\ref{eq:photovoltage}) reduces to the DS result for an ordinary 2DEG~\cite{dyakonov_ieee_1996a} (with the cyclotron mass $m_{\rm c}$ playing the role of the electron mass $m$) in the limit $C_1 = C_2 =0$.
The terms controlled by $C_2$ stem from the local fluctuations of the cyclotron mass.
This coefficient tends to unity in the limit $v_{\rm F} \ll 1$.
The terms controlled by $C_1$ originate from the derivatives of the pressure and suppress the photoresponse with respect to the case $C_1 =0$.

Illustrative plots of $\Delta U$ as a function of the incoming radiation frequency $\Omega$ are presented in Fig.~\ref{fig:photovoltage} for a given value of the momentum relaxation time $\tau$.
The photoresponse features maxima at odd multiples of the plasma frequency (\ref{eq:omegaPvF}).
Nervous readers can easily convert dimensionless values of the momentum relaxation time $\tau$ into values of the mobility $\mu$.
The mobility $\mu$ is obtained from Eq.~(\ref{eq:eulerlin}) in the steady-state regime (i.e.~when the time derivative on the left-hand side is neglected) and reads $\mu = (1 + v_{\rm F}^{2}/2) \tau$.
Restoring, for a moment, physical units, the previous relation reads:
\begin{equation}
\mu =  \frac{e \tau}{m_{\rm c}} \left(1 + \frac{v_{\rm F}^{2}}{2s^2} \right) = \frac{e}{\hbar} \frac{v_{\rm F} \tau}{\pi n d} \left (d\sqrt{\pi n} + \frac{1}{8 \alpha_{\rm ee}} \right )~.
\end{equation}
More conveniently,
\begin{equation}
\mu = 5.0 \times 10^{5} \frac{\tau}{nd} \left ( 0.18 d \sqrt{n} + \frac{1}{8 \alpha_{\rm ee}} \right )~{\rm cm}^2/({\rm V} {\rm s})~,
\end{equation}
where $\tau$, $n$, and $d$ must be expressed in ${\rm ps}$, $10^{12}~{\rm cm}^{-2}$, and ${\rm nm}$, respectively.

The photovoltage (\ref{eq:photovoltage}) depends sensibly on the ratio $\xi$ between the Fermi velocity $v_{\rm F}$ and the effective plasma-wave speed $s(v_{\rm F})$.
From Eq.~(\ref{eq:xi}) we see that this ratio can be tuned in the range $[0,\sqrt{2}]$ by changing the distance $d$ between graphene and the gate, the dielectric constant $\epsilon_1$, and the average carrier density ${\bar n}$.
Fig.~\ref{fig:photovoltage} shows that the photoresponse of a system that is described by relativistic hydrodynamics can be tuned to regimes where its resonant maxima have larger amplitude and its minima are shallower.
We are therefore naturally led to introduce the visibility ${\cal V}$ of the resonant photovoltage profile according to
\begin{equation}\label{eq:visibility}
{\cal V} \equiv \frac{\Delta U_{\rm max} - \Delta U_{\rm min}}{\Delta U_{\rm max} + \Delta U_{\rm min}}~,
\end{equation}
where $\Delta U_{\rm max}$ ($\Delta U_{\rm min}$) is the maximum (minimum) of the photoresponse evaluated at $\Omega \approx \Re e[\omega_1]$ ($\Omega \approx \Re e[\omega_2]$)---see Eq.~(\ref{spagnasuca}).
Illustrative plots of ${\cal V}$ are shown in Fig.~\ref{fig:visibility}.
We see that, for all values of the momentum relaxation time $\tau$, the visibility drastically increases as $\xi$ increases.
In particular, in the limit $\xi \to \sqrt{2}$ the visibility reaches its maximum value. Reaching values of $\xi$ of the order of $\sqrt{2}$ requires to minimize the product $dk_{\rm F}$. Values of this product as small as $0.1$ can be reached by placing graphene on high-quality substrates like h-BN, where carrier densities as low as $\sim 10^{10}~{\rm cm}^{-2}$ can be reached without entering the regime where disorder
(electron-hole puddles) dominates. In these samples, gold metal gates can be placed as close as $10~{\rm nm}$ to the graphene sheet, as recently shown in a quantum capacitance measurement carried out in single-layer graphene~\cite{manchester}.

\section{Conclusions and perspectives}
\label{sect:conclusions}

In this Article we have presented a theory of a photoresponse mechanism which is based on the excitation of plasma waves in the channel of a graphene field-effect transistor (FET).
The FET is coupled to a source of radiation which periodically modulates the voltage difference between the gate and source.
Our main result is a fully analytical expression---Eq.~(\ref{eq:photovoltage})---for the dc voltage difference $\Delta U$ between the drain and source, in response to the oscillating radiation field.

Plasma waves in the channel are described within a hydrodynamic model which takes into account the linear dispersion of massless Dirac fermions in graphene.
Formally, the hydrodynamic theory of two-dimensional massless Dirac fermions in graphene presents strong analogies with the hydrodynamic theory of a relativistic fluid~\cite{Landau06}, with the Fermi velocity $v_{\rm F}$ playing the role of the speed of light $c$.
In particular, we have shown that this yields a breakdown of Galilean invariance whereby the speed of plasma waves in the frame of the moving fluid is not connected by a Galilean transformation to that in the laboratory frame.
Moreover, the photoresponse substantially depends on the peculiar nonlinearities of the equations of motion which arise due to the formal relativistic nature of the energy-momentum dispersion of carriers in graphene.
We have shown that it is possible to leverage the effect of these nonlinearities to increase the photoresponse of a gated graphene sheet with respect to that occurring in an ordinary two-dimensional electron gas.
Measuring the photoresponse of a graphene sheet in a FET geometry offers the opportunity to demonstrate the importance of relativistic corrections to the ordinary hydrodynamic theory.

Finally, we emphasize that the present theory assumes that the gate of the FET is as long as the channel.
In certain experimental situations~\cite{vicarelli_naturemat_2012} it may be more convenient to exploit gates which are much shorter than the FET channel.
In this case the theory of this Article needs to be changed.
The dependence of the DS-like dc photovoltage on the position of the gate with respect to the source, say, may present non-trivial features.
In the case of a short gate, however, the dc photovoltage will suffer~\cite{vicarelli_naturemat_2012} from contributions of thermoelectric and photoconductive origin.
Last but not least, we would like to emphasize that Vicarelli {\it et al.}~\cite{vicarelli_naturemat_2012} have also shown that {\it bilayer graphene} FETs can be used to fabricate plasma-wave THz photodetectors with low noise-equivalent power.
In the case of bilayer graphene, we note that modulating the gate-to-source voltage will not only modulate the density in the FET channel but also the band gap~\cite{fogler_prb_2010}.
These issues are well beyond the scope of the present Article and will be investigated in forthcoming publications.

\begin{acknowledgments}
It is a pleasure to thank Vittorio Pellegrini and Alessandro Tredicucci for useful discussions.
This work was supported by the Italian Ministry of Education, University, and Research (MIUR) through the program ``FIRB - Futuro in Ricerca 2010'' Grant No.~RBFR10M5BT (``PLASMOGRAPH: plasmons and terahertz devices in graphene'').
We have made use of free software (www.gnu.org, www.python.org).
\end{acknowledgments}


\section{Introduction}

\noindent Characterizing how the brain organizes its activity to carry out complex cognitive tasks is a highly non trivial pursuit. While early neuroimaging studies typically aimed at identifying patches of task-specific activation or local time-varying patterns of activity, it is now recognized that both at rest and during behaviour, that brain activity is characterized by the formation and dissolution of functionally meaningful integrated activity \cite{bre2002}.

Over the last ten years, functional activity has started been investigated using a statistical physics understanding of a very old branch of pure mathematics: graph theory. This approach has extensively been applied to complex physical and biological
systems \cite{new2003,boc2006,alb2007}, and represents a shift from a parallel distributed, computer-like to a complex system vision of the brain. In this approach, functional activity is thought of as a network. A network consists of a number of units (nodes) interacting by links. The specific pattern of connections between the nodes defines the network’s topology.  Network theory provides suitable descriptors at all scales from local microscopic, to mesoscopic, to global macroscopic ones. Models to reproduce global topological properties of natural systems e.g. the “small-world” \cite{wat1998} and “scale-free” \cite{bar1999} properties, and measures to properly quantify these and other observed features have been proposed \cite{boc2006}.
The existence of central regions in brain functional networks and the way of quantifying their importance in the network structure has been the object of numerous studies \cite{spo2007,ach2006,buc2009,zam2010,joy2010,loh2010,hwa2012,kuh2012}. The centrality (i.e. importance) of a node $i$ can {\it prima facie} be measured by the number of connections $K(i)$ (degree) in the functional network, indicating how many regions are coordinated with it during a certain task. Once the amount of synchrony is quantified, a weighted network of connections can be constructed. In this case, the node strength $S(i)$, i.e. the sum of the weights of its links, becomes a more accurate indicator of the importance of a node in the network. Nevertheless, both the degree and the strength are local measures, which do not take into account the global topology of the functional network.
To overcome this issue, global measures of centrality have been proposed. Node closeness $C(i)$, takes into account the number of steps that we have to do to go from one to another: the lower number of steps, the higher closeness a node has. The eigenvector centrality $E(i)$, is another indicator of the global relevance of a node and it is measured by computing the eigenvector $v_1$ associated to the first eigenvalue of the connectivity matrix (i.e., the matrix of interactions between nodes). A more intuitive measure is represented
by node betweenness $B(i)$, which quantifies the number of shortest paths (minimum number of steps between two nodes) that cross a certain node, thus reflecting its importance in the transmission of information in the whole network. Closeness, eigenvector centrality and betweenness, make use of information about the structure of the whole network, nevertheless node centrality could be restricted solely to a certain region or community. To measure the community importance of nodes, Guimer\`a {\em et al.} \cite{gui2005} introduced two indicators, the within-module-degree $z_i$ and the participation coefficient $p_i$. The combination of both parameters not only measures the importance of a node inside its community but allows classifying the role played by the central nodes (hubs) in the overall community structure.

The application of these measures to functional brain networks may help in quantifying how relevant nodes are impaired by the emergence of different brain diseases. Here, we analyze the alteration of functional networks' hubs caused by mild cognitive impairment (MCI), a brain syndrome involving cognitive impairments beyond those expected based on the age and education of the individual often representing a transitional stage between normal aging and dementia.
The progressive accumulation of the beta amyloid protein and the loss of cells and synapses along the Alzheimer's disease (AD) continuum impair cognitive function as a reflection of network malfunctioning. Thus, although hubs organization has been studied for fully-fledged AD \cite{buc2009,sta2009,deh2012}, whether hub organization impairment at the early stages of the disease such as in MCI.
The MEG signal of seventeen patients suffering from MCI and seventeen control patients, was recorded during the execution of a memory task (see \cite{bul2011} for details). Next, we measured the synchronization likelihood ($SL$) between all pairs of nodes, and obtained a weighted correlation matrix that was analyzed by using Complex Networks Theory methods.
In Buld\'u {\em et al.}, the attention was devoted to the global properties of the functional networks. Nevertheless, how hub organization is affected by the disease was not considered. In the present work we focus on how the role of the hubs is distorted by MCI.
Our hub characterization relies on various parameters, but we show how the eigenvector centrality $E(i)$ of the nodes, a global measure which quantifies the importance of a node in the whole structure of the network, is the more adequate indicator, in contrast with local measures such as degree or strength centrality.
Therefore, we use this measure of centrality to identify (and quantify) network hubs.
Eigenvector centrality allows detecting the dominant regions in the functional connectivity network and, next, how these regions are affected by the disease. Our results show that the leading role played by network hubs is attenuated by MCI, which transforms the network into a more homogeneous one. These results are consistent with previously reported analyses showing that MCI increases the randomness of the global structure of the functional network \cite{bul2011}, thus, leading to a more homogeneous connectivity network.

The manuscript is organized as follows: in Section \ref{functional} we explain the methodology followed to obtain functional networks from magnetoencephalography measurements. Section \ref{methods} includes all
information related with the experimental procedure. Section \ref{centrality} contains a detailed explanation of the different centrality measures used in this work together with the main results. Finally, in Section
\ref{discussion}, we discuss the implications of the results obtained and the relation with other studies
of hub characterization in functional brain networks.

\section{Obtaining functional networks}
\label{functional}

Functional brain networks rely on the coordinated activity between brain regions, which requires to quantify the
synchronization of each pair of regions, i.e. nodes, within the network. There exist a number of different
measures of synchronization \cite{boc2002,per2005} that can be applied to the recorded time series of each node (see Section \ref{methods} for details on how the dynamics of the nodes are obtained). Among them, we have
chosen Synchronization Likelihood ($SL$) \cite{sta2002}, a nonlinear measure of the synchronized activity that has been proven to be a suitable quantifier for datasets obtained from
magnetoencephalographic recordings \cite{sta2006,bul2011}. This index, which is closely related to the concept of generalized mutual information \cite{buz1994}, relies on the detection of simultaneously occurring patterns, which can be complex and widely different for
every pair of signals (see the example of Fig. \ref{fig:fig01}). Let \textbf{X} denote the matrix containing all the $M$ signals (one per each channel that measures a certain cortical region) of $n$ time steps, and let $X_n=[x_n, y_n,...,z_n]$, where $X_{1,n}=x_n$, $X_{2,n}=y_n$ and so on. In the time series of a given channel, for each time step $n$, we define the probability that embedded vectors are closer to each other than a distance $\varepsilon$:

\begin{equation}
P_{n}^{\varepsilon}\left( X\right) =\dfrac {1} {2\left( w_{2}-w_{1}\right)} \sum_{\substack{m=1\\w_1<\vert n-m\vert<w_2}}^N\theta \left( \varepsilon -\left| x_{n}-x_{m}\right| \right)
\end{equation}

where $\vert.\vert$ is the Euclidean distance and $\theta$ is the Heaviside step function, being $\theta(x)=0$ if $x\leq0$ and $\theta(x)=1$ for $x>0$. Here $w_1$ and $w_2$ are two windows: $w_1$ is the Theiler correction for autocorrelation effects and should be at least of the order of the autocorrelation time; $w_2$ is a window that sharpens the time resolution of the synchronization measure and is chosen such that $w_1\ll w_2\ll N$.

\begin{figure}[htbp]
\begin{center}
\vskip1.2cm
\psfig{file=fig01.eps,width=0.7\textwidth}
\vskip0.5cm
\end{center}
\caption{Synchronization likelihood ($SL$) between two brain regions. Time series correspond to the recorded electromagnetic field
of two brain sites. The $SL$ algorithm detects windows of repeated patterns within the time series of channel A and, next, checks whether the
channel B also shows a repeated pattern at the same time windows, no matter if it is the same or different to that observed in channel A.
Values of $SL$ close to one (zero), indicate a high (low) correlation.}
\label{fig:fig01}
\end{figure}

Now, for each signal and each time $n$, the critical distance $\varepsilon_n$ is determined for which $P_{n}^{\varepsilon}(x)= p_{ref}$, where $p_{ref}\ll1$. We can determine for each discrete time pair $(n,m)$ within our considered window $(w_1 < \vert n-m \vert <w_2)$ the number of channels $H_{n,m}$ where the embedded vectors $x_{k,n}$ and $x_{k,m}$ will be closer together than this critical distance $\varepsilon_{k,i}$:

\begin{equation}
H_{n,m}=\sum _{k=1}^{M}\theta \left( \varepsilon _{k,n}-\left| x_{k,n}-x_{k,m}\right| \right)
\end{equation}

This number lies in a range between $0$ and $M$, and reflects how many of the embedded signals ``resemble'' each other.

We can now define a \textit{synchronization likelihood} $SL_{n,m}(X)$ for each channel $k$ and each discrete time pair $(n,m)$ as:

\begin{equation}
\begin{aligned}
& \mathrm{if} \left| x_{n}-x_{m}\right| < \varepsilon_{n}:SL_{n,m}\left( X\right) =\dfrac {H_{n,m}-1} {M-1}\\
& \mathrm{if} \left| x_{n}-x_{m}\right| \geq \varepsilon_{n}:SL_{n,m}\left( X\right) =0
\end{aligned}
\end{equation}

By averaging over all $m$, we finally obtain the synchronization likelihood, $SL_n(X)$:

\begin{equation}
SL_{n}\left(X\right) =\dfrac {1} {2\left( w_{2}-w_{1}\right)} \sum_{\substack{m=1\\w_1<\vert n-m\vert<w_2}}^N SL_{n,m} \left( X \right)
\end{equation}

Synchronization likelihood $SL_n(X)$ describes how strongly channel $x$ at time $n$ is synchronized to all the other $ M-1 $ channels.
The range of values of $SL$ is $0\leq SL\leq 1$, being ($p_{ref}$) when all $M$ time series are uncorrelated, and 1 for maximal synchronization of all $M$ time series. The value of $p_{ref}$ can be set at an arbitrarily low level, and does not depend on the properties of the time series, nor it is influenced by the embedding parameters \cite{sta2002}.

The $SL$ yields a symmetric and weighted correlation matrix $w_{ij}$, which can be analyzed using Complex Networks metrics \cite{boc2006}.
Before computing their network parameters, we follow the normalization technique proposed in \cite{bul2011} in order to avoid intrinsic differences from different individuals. The off-diagonal weights in the correlation matrices are rescaled to the interval $[0,1]$ by means of

\begin{equation}
SL_{ij}=\dfrac{w_{ij}-min(w_{ij})}{ max(w_{ij})-min(w_{ij})}
\end{equation}

\noindent while the diagonal is set to zero after normalization.
Figure \ref{fig:fig02} shows an example of a functional brain network obtained from a control group of seventeen healthy individuals
(see Section \ref{methods} for details).
The normalized weighted correlation matrix $W_{ij}$ leads to a fully connected network (since all pair of nodes have a $SL$ $\geq 0$), so only
the $5\%$ of the links with higher weights (i.e., $SL$) have been plotted to ease visualization of the network. The size of a node $i$ is proportional to its strength $S(i)$, defined as the sum of the weights of all its connections.
We can observe how the density of connections is higher at the peripheral regions and, specially, at the occipital lobe, while the central
cortical region is sparsely connected (note that we are only considering the links with higher correlations). Interestingly, nodes with higher strengths,
i.e., the network hubs, are mainly localized in the occipital lobe.

\begin{figure}[htbp]
\begin{center}
\psfig{file=fig02.eps,width=0.4\textwidth}
\end{center}
\caption{Averaged functional network of the control group. Only $5\%$ of the links with higher weights have been plotted. Colors indicate the
lobe that a node belongs to: frontal-left (blue), frontal-right (light green), central (dark green), temporal left (red), temporal-right
(light orange) and occipital (dark grey). Node size is proportional to the node strength S(i). Note the higher density of connections
at the occipital and frontal lobes. Node number is indicated for those nodes with higher strength.}
\label{fig:fig02}
\end{figure}

\section{Materials and Methods}
\label{methods}

The magnetoencephalographic (MEG) signal was recorded with a 256 $Hz$ sampling frequency and a band pass of [$0.5,50$] $Hz$, using a $148$-channel whole-head magnetometer (MAGNES\copyright $2500$ WH, $4-D$ Neuroimaging) confined in a magnetically shielded room. Seventeen MCI patients and seventeen healthy volunteers were recorded during a Sternberg's letter-probe task \cite{det1991,mae2001}.

MCI diagnosis was established according to the criteria proposed by Petersen {\it et al.} \cite{gru2004,pet2004}. Thus, MCI patients fulfilled the following criteria: (1) cognitive complaint corroborated by an informant (a person who stays with the patient at least for half a day at least 4 days a week); (2) objective cognitive impairment, documented by delayed recall in the logical memory II sub-test of the revised Wechsler Memory Scale (score $\leq16/50$ for patients with more than 15 years of education; score $\leq8/50$ for patients with $8-15$ years of education); (3) normal general cognitive function, as assessed by a clinician during a structured interview with the patient and an informant and, additionally, a mini mental state examination (MMSE) score greater than 24; (4) relatively preserved daily living activities as measured by the Lawton scale; (5) not sufficiently impaired, cognitively and functionally to meet criteria for dementia. Age and years of education were matched to the SMC group. According to their clinical and neuropsychological profile, all patients in this group were considered multi-domain MCI patients (see \cite{pet2004}). As for the geriatric depression scale, none of the MCI showed depression (score lower than 9) \cite{yes1991}.

In the Sternberg's letter-probe task \cite{det1991,mae2001} a set of five letters was presented and participants were asked to keep the letters in mind. A series of single letters ($1 s$ in duration with a random ISI between $2$ and $3$ $s$) was then introduced one at a time, and participants were asked to press a button with their right hand when a member of the previous set was detected.

Thereafter, single trial epochs were visually inspected by an experienced investigator, and epochs containing visible blinks, eye movements or muscular artifacts were excluded from further analysis. Thirty-five epochs corresponding to each subject were used in order to calculate the functional connectivity values (i.e., their synchronization likelihood). This lower bound was determined by the participant with least epochs. To have an equal number of epochs across participants, thirty-five epochs were randomly chosen from each of the other participants.
The effect of plasticity in the evaluation of network synchronization is not dealt with in this work, although it could be a subject of analysis in future ones.

Before the MEG recordings, all participants or legal representatives gave written consent to participate in the study, which was approved by the local ethics committee of the Hospital Cl\'inico San Carlos (Madrid, Spain).

\section{Measuring node centrality}
\label{centrality}

The emergence and evolution of MCI has been studied extensively during the last years \cite{pet2004,dic2005,bab2006,sch2007,sch2009,baj2010} since it is known to be related to prodromal Alzheimer's disease (AD).
Neuropathological studies indicate that MCI patients have clear pathophysiological characteristics,
such as the presence of neurofibrillary tangles, loss of dendritic
spines and the accumulation of beta-amyloid protein in the
associative cortex \cite{mar2010}. From the point of view of brain connectivity, MEG recordings of patients
suffering from MCI revealed
an enhancement of the synchronization between cortical regions when memory tasks were
performed \cite{baj2010}. The increase of the synchronized behaviour was accompanied with a change
of the topological structure of the associated functional network, which turned into more random configurations
\cite{bul2011}. A reduction of the modular behaviour of the network together with an increase
of the long-range functional connections has also been associated with the appearance of MCI \cite{bul2011}.

Nevertheless, how functional hubs are affected by the disease is still unclear. It is known
that network hubs are strongly affected by brain diseases like schizophrenia \cite{bas2008} or
Alzheimer's disease \cite{buc2009,sta2009,deh2012}. This is bad news, since
the targeted attack to leading nodes leads to a fast damage of the whole network properties \cite{alb2000} and, eventually, to cascading failures \cite{mot2002}.

We used four different centrality measures from complex weighted-network theory
to detect network hubs and to evaluate how their topological roles are affected by the disease. These measures were computed over the matrices
obtained from the patients and control individuals.

The most straightforward centrality measure is the strength of nodes, $S(i)$, defined as the sum of the weights
of the links to all neighbors of a node $i$. If a node has strong connections with its neighbors, it will have higher influence on the functional network, i.e. it will be more central. $S(i)$ can be easily obtained as:
\begin{equation}
S(i)=\sum_{j\in N(i)}W_{ij}
\end{equation}
where $N(i)$ is the set of neighbors of $i$, in our case, all nodes of the network, since ${\bf W}$ matrices are fully connected. Figure \ref{fig:fig03}A shows the strength $S(i)$ of both control and MCI groups, averaged
over the whole groups. Network hubs can be easily identified as those with higher network strength, while the influence of the disease in the hubs is directly the difference of $S(i)$ between both
groups. We can observe that strength of the network hubs (highest peaks of the $S(i)$ distribution) is not specially altered, while from node 1 to 40 (localized at the central lobe) there is a significant increase of
strength, although they are not hubs of the network. Therefore, node strength seems to indicate that MCI is not specially severe with the
hubs.

\begin{figure}[tp]
\begin{center}
\vskip1.0cm
\psfig{file=fig03.eps,width=0.7\textwidth}
\vskip0.5cm
\end{center}
\caption{Averaged centrality measures for the control (black) and MCI (red) groups. Specifically, we calculate the node strength
$S(i)$ (a), node closeness $C(i)$ (b), node betweenness $B(i)$ (c) and node eigenvector centrality $E(i)$ (d). All measures have been calculated and averaged over each epoch and for each subject of both groups. Note that $E(i)$ is the measure that better captures the differences between groups in those nodes with higher centrality. See Fig. \ref{fig:fig02} to locate the position of the hubs.
The periodicity reported in the four figures is only a consequence of the node numbering. See Fig. 2 for details on where nodes are placed in the Euclidean space.
}
\label{fig:fig03}
\end{figure}

Nevertheless, network strength $S(i)$ is a local measure and
does not accounts for correlations of link strength or specific structural distributions of the weights in the network.
To overcome this issue, two measures of centrality have been applied, with the aim of including
the influence of the network structure from a point of view of information
transfer: the node closeness $C(i)$ and the node betweenness
$B(i)$. Both measures are based on the geodesic distance from node to node $d_{ij}$, which is obtained
as the shortest distance to arrive from one node $i$ to any other ($j$) by jumping through nodes of the network. At each jump, the traveled distance is increased by the inverse of the weight of the link that is used to pass from one node to the other (e.g., when jumping from node $k$ to node $l$ the travelled distance increases $\Delta d (kl) = 1/w_{kl}$). Finally, the mean geodesic distance of a node $i$, is defined as $l_i=\frac{1}{n}\sum_{j}d_{ij}$. If the mean of geodesic paths of vertex $i$ is low, it means that this vertex can interchange information faster and farther within the network. Closeness centrality $C(i)$ is then defined as the inverse of the mean geodesic path,
\begin{equation}
C(i)=\frac{1}{l_i}=\dfrac{n}{\sum_{j}d_{ij}}
\end{equation}

A closely related, but not equal, measure of node importance is the betweenness centrality $B(i)$, which measures the ability of a node to  act as transmitter of information through the network. $B(i)$ is defined as
\begin{equation}
B(i)=\sum_{j\ne k\ne i}\frac{n_{jk}^{i}}{g_{jk}}
\end{equation}

where $n_{jk}^{i}$ accounts for the number of shortest paths between every pair of nodes $j$ and $k$ that pass through node $i$, properly normalized by the total number of geodesic paths $g_{jk}$ from $j$ to $k$.

In Fig. \ref{fig:fig03}B-C we plot the average values of $C(i)$ and $B(i)$ for control and MCI groups. While both measures have their proper characteristics, a common behavior appears in the high correlation present among the centralities of control and MCI for all nodes. The position of the network hubs perfectly matches when
comparing the local measure of centrality $S(i)$ with $C(i)$, while $B(i)$ seems to suggest the existence
of hubs but with a more noisy distribution. Interestingly, the peaks reported in the three distributions correspond to nodes placed at the occipital lobe (see node number in Fig. \ref{fig:fig02}). Again, there are no
significant differences in the centrality of the network hubs, and the most clear signature of MCI is an
increase of centrality in the nodes belonging to the central lobe (from 1 to 40). This alteration is not
captured by $B(i)$, thus indicating that it is not the most suitable measure to evaluate changes in network centrality.

\begin{figure}[tp]
\begin{center}
\vskip1.0cm
\psfig{file=fig04.eps,width=0.7\textwidth}
\vskip0.5cm
\end{center}
\caption{Correlation between the eigenvector centrality $E(i)$ and its variation $\Delta E(i)=E(i)_{MCI}-E(i)_{control}$. Note the
negative correlation indicating that nodes with higher eigenvector centrality are those that are affected the most by the disease.
}
\label{fig:fig04}
\end{figure}

It is worth noting that although $C(i)$ and $B(i)$ contain information about how weights are distributed within the network, they only refer to shortest paths distribution and disregard other structural properties of the network.
To overcome this issue, we compute the eigenvector centrality $E(i)$ of the nodes. $E(i)$ is a measure
obtained from the spectral analysis of the connectivity matrix ${\bf W}$. Specifically, it is calculated as the eigenvector $v_1$ of the first eigenvalue of the matrix ${\bf W}$. Note that any modification of the weight of the links, will result in a change of the matrix ${\bf W}$ and, therefore, will be reflected in the value
of the eigenvector centrality $E(i)$. On the contrary, if an alteration of a link does not imply any modification in the distribution of shortest paths, it will not be captured by $C(i)$ or $B(i)$ centralities.
The mathematical definition of the eigenvector centrality
is related to an iterative process were the centrality of a node $i$ is calculated as the sum of the centralities of its neighbors:
\begin{equation}
 E(i)= \lambda^{-1} \sum_j W_{ij} E(j),
\label{eq:centrality}
\end{equation}
where $\lambda$ is a constant, $E(j)$ is the eigenvector centrality of node $j$, being $j$ a neighbour of $i$ and $W_{ij}$ are the components of the connectivity matrix ${\bf W}$. In matrix notation, Eq. \ref{eq:centrality} reads
{\bf $\lambda E= W E$} so that {\bf $E$} can be expressed as a linear combination of the eigenvectors {\bf $v_i$} of the adjacency
matrix. Since the final value of $E(i)$ is obtained when $t \rightarrow \infty$, the value of $E(\infty)$ is equal to the eigenvector $v_1$ associated to the dominant eigenvalue $\lambda_1$. Therefore, the eigenvector centrality $E(i)$ is equal to the eigenvector $v_1(i)$ of the connectivity matrix $W$ \cite{new2004}.

Figure \ref{fig:fig03}D shows the eigenvector centrality $E(i)$ for both groups. It detects the position of the network hubs which, as in $S(i)$ and $C(i)$, are mainly placed at the occipital lobe. Interestingly, this measure captures a clear decay in the centrality of the hubs, as indicated by the decrease of the peaks height. The consequence of this decrease is a network were the hubs play a less relevant role. In Fig. \ref{fig:fig04} we
plot how the variation of the eigenvector centrality $\Delta E (i)$ due to MCI is related with the initial centrality of the nodes. We can observe a negative correlation, indicating that those nodes with higher centrality are, in turn, the nodes that decrease its importance the most. It is also interesting to see what is the position of the nodes in the whole network. Figure \ref{fig:fig05} shows the increase (green) and decrease
(red) of eigenvector centrality with the node size proportional to the value of $|\Delta E(i)|$. It is the occipital lobe the one that is mostly affected by the decrease of centrality, while the central lobe takes
advantage of this reduction.

\begin{figure}[tbp]
\begin{center}
\psfig{file=fig05.eps,width=0.4\textwidth}
\end{center}
\caption{Distribution of the variation of eigenvector centrality $E(i)$. The size of the nodes is proportional
to $|\Delta E(i)|=|E(i)_{MCI}-E(i)_{control}|$, while node color indicates an increase (green) or decrease (red) in the node
centrality. Node number is indicated for those nodes with high $|\Delta E(i)|$.}
\label{fig:fig05}
\end{figure}

\section{Conclusions}
\label{discussion}

The analysis of functional and anatomical brain networks using complex networks analysis have revealed
interesting information about how these networks are organized \cite{bul2009}. All studies agree on the fact that
brain networks are highly heterogeneous, giving rise to the existence of hubs, i.e. leading nodes
from the point of view of the network structure. In this manuscript we have shown a detailed analysis of how
the importance of network hubs in a functional network associated with the execution of a memory task
is affected by the
emergence of MCI. We report how the disease particularly affects network hubs, reducing their importance
in the network.
This reduction is captured by all different measures of centrality, independently of whether they rely on the local or global properties of the network.
Among all centrality measures, eigenvector centrality $E(i)$ is the one that better captures
the effect of the disease on the network hubs.
The occipital region, containing the majority of the hubs in the case of healthy individuals,
is the region where centrality decreased the most, while nodes belonging to the central lobe benefit from the hub deterioration. Interestingly, anatomical networks of patients suffering from schizophrenia also show
a reduction of the hub importance, together with a dispersion in the location of network hubs \cite{bas2008}.
This reduction of hub centrality can be used as a signature of the existence of pathology since it is known
that the hub structure is quite stable in healthy adult individuals \cite{hwa2012,zuo2012}.
In addition, it reinforces the hypothesis that MCI increases the randomness and
homogeneity of the functional network \cite{bul2011}, since the reduction of hub importance leads to more homogeneous networks.
The implications of hub deterioration are still unclear. For example, in Alzheimer's disease, it has been shown that amyloid-beta deposition in the locations of cortical hubs could act critically in the severity of the disease
\cite{buc2009}. In the case of anatomical networks in schizophrenia, the reduction of centrality of the
frontal hubs has been related to the disorganization of the anatomical network.
The reduction in hub centrality in posterior scalp regions could be related to the high levels of beta amyloid accumulation in those regions found in MCI and Alzheimer's disease.
Conversely, the increased hub centrality found for anterior sensors could reflect enhanced engagement of frontal regions compensating for decreased capacity to tackle the demands of the memory task (see \cite{baj2010} for a similar interpretation).
We believe that our work will shed some light on the
role played by hubs in brain networks, which may have strong influence on the network
robustness \cite{alb2000}, cascading processes \cite{mot2002} and network controllability \cite{tan2012}.
Similar studies could be carried out in other neurodegenerative diseases where the existence of functional hubs has been reported.

\nonumsection{Acknowledgments}
This work was supported by the Spanish Ministry of S\&T [FIS2009-07072], the Community of Madrid
under the R\&D Program of activities MODELICO-CM [S2009ESP-1691], the Spanish Ministry of economy and competitiveness
[PSI2012-38375-C03-01], as well to a grant from the CAM [S2010/BMD-2460].

\newpage