Current-induced birefringent absorption and non-reciprocal plasmons in graphene

The photonic, optoelectronic, and plasmonic properties of graphene have been attracting a great deal of attention for the past ten years [1–10]. This interest was sparked by a series of key experimental results that were demonstrated soon after the isolation of graphene [11]. Early on, it was shown [12, 13] that the optical conductivity of graphene is universal in the limit ${\hslash }\omega \gt 2{E}_{{\rm{F}}}$ and given by ${\sigma }_{{\rm{uni}}}={e}^{2}/(4{\hslash })$, where ω is the photon frequency and ${E}_{{\rm{F}}}$ is the Fermi energy. In the same frequency range, absorption is given by $\pi {\alpha }_{{\rm{QED}}}\simeq 2.3\%$, where ${\alpha }_{{\rm{QED}}}={e}^{2}/({\hslash }c)\simeq 1/137$ is the quantum electrodynamics fine structure constant. Importantly, absorption can be tuned by gating [14] and can be enhanced by using a variety of plasmonic protocols [15–20] or by embedding graphene in a optical cavity [21, 22]. At low photon energies (i.e. for ${\hslash }\omega \ll 2{E}_{{\rm{F}}}$), absorption occurs only due to the excitation of intra-band electron-hole pairs, which is allowed by momentum-non-conserving collisions.

More recently, Dirac plasmons, the self-sustained density oscillations of the two-dimensional (2D) electron system in a doped graphene sheet, have been studied experimentally by a variety of spectroscopic methods [4]. Fei et al.[23] and Chen et al.[24] carried out seminal scattering-type scanning near-field optical microscopy (s-SNOM) experiments in which Dirac plasmons were launched and imaged in real space. They showed that, in the mid infrared, the plasmon wavelength ${\lambda }_{{\rm{p}}}$ can be ∼40-60 times smaller than the free-space excitation wavelength ${\lambda }_{0}=2\pi c/\omega$, allowing an extreme concentration of electromagnetic energy, and that Dirac plasmon properties are easily gate tunable. These figures of merit have been tremendously improved by utilizing vertical stacks comprising graphene sheets encapsulated in hexagonal boron nitride crystals [25, 26], in the quest for a low-loss and tunable platform for infrared and THz plasmonics. Mid-infrared plasmons in encapsulated graphene sheets display [27] ultra-large field confinement (i.e. ${\lambda }_{{\rm{p}}}/{\lambda }_{0}\sim 1/150$), small group velocity, and a long lifetime exceeding $500\;{\rm{fs}}$.

In this article we revisit the optical and plasmonic properties of single-layer graphene which carries a dc current. Our interest is motivated by the simple observation that graphene sheets can support ultra-high current densities without being damaged. Indeed, early on Moser et al.[28] proposed a simple method to suppress contamination of graphene at low temperatures inside a cryostat. Applying a source-drain bias of a few volts across their samples, these authors induced a large current flow of up to $4\;{\rm{mA}}$, thereby removing contamination adsorbed on the surface. For a typical device width of $4\;\mu m$, such current translates into an extremely large current density $J\simeq 10\;{\rm{A}}/{\rm{cm}}$.

Here, we present a detailed theoretical study of the optical and plasmonic properties of current-carrying graphene sheets. Previous authors [29, 30] have considered the same problem, but have only reported approximate solutions that are valid for low current densities. Here, we describe a semi-analytical approach, which can be applied for any current density (i.e. for drift velocities $| {\boldsymbol{v}}|$ as large as the graphene Fermi velocity ${v}_{{\rm{F}}}$) and at an arbitrary electron temperature T. The combined functionalities offered by electrical (or chemical) doping and dc currents may further enrich the graphene optoelectronic and plasmonic portfolio, possibly leading to interesting discoveries and/or potential applications.

Our article is organized as follows. In section 2 we present analytical expressions for the imaginary part of the density-density response function of a current-carrying graphene sheet at zero temperature. The calculation of the real part is reduced to a one-dimensional quadrature. Results at finite T are then obtained by using an integral identity due to Maldague [31, 32]. In section 3 we use the results of section 2 to calculate the absorption as a function of photon frequency ω, drift velocity ${\boldsymbol{v}}$, and temperature T. We also calculate analytically the dependence of this quantity on the angle between the polarization of light and current direction. In section 4 we present our results for the plasmon modes of a current-carrying graphene sheet and quantify the degree of in-plane anisotropy induced by the current flow. In section 5 we present mathematical equations relating the chemical potential μ and drift velocity ${\boldsymbol{v}}$ to the carrier density ${n}_{{\rm{c}}}$ and current density ${\boldsymbol{J}}$. Finally, in section 6 we summarize our main findings and discuss potential obstacles against their experimental observation.

The optical properties of a current-carrying graphene sheet can be calculated by using linear response theory [32] with respect to the electric field of the incident light beam. We treat the effect of an external dc current in a non-perturbative fashion. To this end, we first find the distribution function of a homogeneous equilibrium state carrying a finite current, at a finite temperature T. The distribution function of such a state can be calculated by requiring that it nullifies the collision integral for electron–electron scattering processes, which conserve energy and momentum [33]. The end result of this calculation is the following Fermi–Dirac distribution function [33, 34]:

$$\begin{eqnarray}&&{n}_{{\rm{F}},\lambda }({\boldsymbol{k}},{\boldsymbol{v}};T)={\left[\mathrm{exp}\left(\displaystyle \frac{{\varepsilon }_{{\boldsymbol{k}},\lambda }-{\hslash }{\boldsymbol{v}}\cdot {\boldsymbol{k}}-\mu }{{k}_{{\rm{B}}}T}\right)+1\right]}^{-1},\end{eqnarray} \tag{ 1 }$$

where ${\varepsilon }_{{\boldsymbol{k}},\lambda }=\lambda {\hslash }{v}_{{\rm{F}}}| {\boldsymbol{k}}|$ is the Dirac band energy, ${v}_{{\rm{F}}}\sim {10}^{6}\;{\rm{m}}/{\rm{s}}$ is the Fermi velocity, and $\lambda =+$ ($\lambda =-$) denotes conduction (valance) band states. The function ${n}_{{\rm{F}},\lambda }({\boldsymbol{k}},{\boldsymbol{v}};T)$ is parametrized by the drift velocity ${\boldsymbol{v}}$ and chemical potential μ. The drift velocity will be related to the current density ${\boldsymbol{J}}=J\hat{{\boldsymbol{v}}}$ with $\hat{{\boldsymbol{v}}}={\boldsymbol{v}}/v$ and $v=| {\boldsymbol{v}}|$ below. Since we are dealing with an open system, i.e. a graphene sheet attached to leads, we do not follow the usual grand-canonical procedure [35] in which the chemical potential μ is forced to readjust at each value of T to make sure that the average number of particles is conserved. In this article, μ is fixed by the leads (and external gating) and does not change with temperature. (For the sake of definiteness, we consider an n-doped graphene sheet: results for p-doped systems are identical since the quantities we are interested in are particle-hole symmetric.)

As we will see below, the drift velocity satisfies the following inequality: $v\lt {v}_{{\rm{F}}}$. Our results, however, are not perturbative in the dimensionless parameter

$$\begin{eqnarray}&&\beta =\displaystyle \frac{v}{{v}_{{\rm{F}}}}.\end{eqnarray} \tag{ 2 }$$

To find the optical response of a current-carrying two-dimensional (2D) massless Dirac fermion (MDF) system [36], we calculate the dynamical density-density response function, which is defined by [32]

$$\begin{eqnarray}&&{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T)={N}_{{\rm{f}}}\displaystyle \sum _{\lambda ,{\lambda }^{\prime }}\displaystyle \int \displaystyle \frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}{F}_{\lambda {\lambda }^{\prime }}({\boldsymbol{k}},{\boldsymbol{k}}+{\boldsymbol{q}})\\ &&\quad \times \displaystyle \frac{{n}_{{\rm{F}},\lambda }({\boldsymbol{k}},{\boldsymbol{v}};T)-{n}_{{\rm{F}},{\lambda }^{\prime }}({\boldsymbol{k}}+{\boldsymbol{q}},{\boldsymbol{v}};T)}{{\hslash }\omega +{\varepsilon }_{{\boldsymbol{k}},\lambda }-{\varepsilon }_{{\boldsymbol{k}}+{\boldsymbol{q}},{\lambda }^{\prime }}+i\eta }.\end{eqnarray} \tag{ 3 }$$

Here, ${N}_{{\rm{f}}}=4$ is the number of fermion flavors in graphene,

$$\begin{eqnarray}&&{F}_{\lambda {\lambda }^{\prime }}({\boldsymbol{k}},{\boldsymbol{k}}+{\boldsymbol{q}})=\displaystyle \frac{1}{2}[1+\lambda {\lambda }^{\prime }\mathrm{cos}({\theta }_{{\boldsymbol{k}}}-{\theta }_{{\boldsymbol{k}}+{\boldsymbol{q}}})]\end{eqnarray} \tag{ 4 }$$

is the usual chirality factor [36], and $\eta ={0}^{+}$ is a positive infinitesimal.

2.1. The imaginary part of ${\chi }_{{nn}}^{(0)}$ at zero temperature

Following standard practice [32], we first calculate analytically the imaginary part of ${\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T)$ at T = 0.

It is convenient to introduce the following dimensionless wave number

$$\begin{eqnarray}&&\bar{q}=\displaystyle \frac{{\hslash }{v}_{{\rm{F}}}q}{\mu }\end{eqnarray} \tag{ 5 }$$

(with $q=| {\boldsymbol{q}}|$) and dimensionless frequency

$$\begin{eqnarray}&&\bar{\omega }=\displaystyle \frac{{\hslash }\omega }{\mu }.\end{eqnarray} \tag{ 6 }$$

After lengthy but straightforward calculations we find the following result:

$$\begin{eqnarray}&&\mathrm{Im}[{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;0)]={f}_{\mu }(\bar{q},\bar{\omega })\displaystyle \sum _{\gamma =\pm 1}U(\bar{q},\bar{\omega };\beta ,\gamma \psi ),\end{eqnarray} \tag{ 7 }$$

where ψ is the angle between ${\boldsymbol{q}}$ and $\hat{{\boldsymbol{v}}}$, i.e.

$$\begin{eqnarray}&&\mathrm{cos}(\psi )=\displaystyle \frac{{\boldsymbol{q}}\cdot \hat{{\boldsymbol{v}}}}{q},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{f}_{\mu }(x,y)=\displaystyle \frac{D(\mu )}{16}\displaystyle \frac{{x}^{2}}{\sqrt{| {x}^{2}-{y}^{2}| }},\end{eqnarray} \tag{ 9 }$$

$D(\mu )={N}_{{\rm{f}}}| \mu | /(2\pi {{\hslash }}^{2}{v}_{{\rm{F}}}^{2})$ is the 2D MDF density of states, and

$$\begin{eqnarray}&&U(\bar{q},\bar{\omega };\beta ,\psi )\\ &&\quad =-{\rm{\Theta }}(\bar{q}-\bar{\omega })\displaystyle \sum _{\alpha =\pm 1}\alpha {H}_{\mathrm{intra}}({A}_{\alpha },B,{C}_{\alpha })\\ &&\qquad +{\rm{\Theta }}(\bar{\omega }-\bar{q})[{H}_{\mathrm{inter}}({A}_{-1},B,{C}_{-1})-\pi ].\end{eqnarray} \tag{ 10 }$$

Analytical expressions for the real functions ${H}_{\mathrm{intra}}(x,y,z)$ and ${H}_{\mathrm{inter}}(x,y,z)$ in equation (10) are reported in tables 1 and 2, while the quantities ${A}_{\alpha }$, B and ${C}_{\alpha }$ are defined as following:

$$\begin{eqnarray}{A}_{\alpha } & = & 2+\alpha \bar{\omega }-\alpha \beta \mathrm{cos}(\psi )\bar{q},\\ B & = & \beta \bar{\omega }\mathrm{cos}(\psi )-\bar{q},\\ {C}_{\alpha } & = & \alpha \beta \mathrm{sin}(\psi )\sqrt{| {\bar{q}}^{2}-{\bar{\omega }}^{2}| }.\end{eqnarray} \tag{ 11 }$$

Table 1.  This table summarizes the expressions that the function ${H}_{{\rm{intra}}}(x,y,z)$ takes depending on the sign of the quantity ${A}_{\alpha }+B$. The left (right) column applies if ${A}_{\alpha }+B\gt 0$ (${A}_{\alpha }+B\lt 0$). The function ${G}_{\mathrm{intra}}(x)=x\sqrt{{x}^{2}-1}-\mathrm{log}(x+\sqrt{{x}^{2}-1})$ is evaluated at $x={x}_{i}^{\alpha }$, where ${x}_{i}^{\alpha }$ is the root of the function ${g}_{\mathrm{intra}}(x)={A}_{\alpha }+{Bx}-{C}_{\alpha }\sqrt{{x}^{2}-1}$. Here, the subscript $i=\pm$ of the root ${x}_{i}^{\alpha }$ is given by the sign of the derivative ${g}_{\mathrm{intra}}^{\prime }(x)$ evaluated at $x={x}_{i}^{\alpha }$, i.e. ${g}_{\mathrm{intra}}^{\prime }({x}_{+}^{\alpha })\gt 0$, while ${g}_{\mathrm{intra}}^{\prime }({x}_{-}^{\alpha })\lt 0$.

${A}_{\alpha }+B\gt 0$	${A}_{\alpha }+B\lt 0$
${G}_{\mathrm{intra}}({x}_{-}^{\alpha })$	${\rm{\Theta }}({A}_{\alpha }-\sqrt{{B}^{2}-{C}_{\alpha }^{2}})[{G}_{\mathrm{intra}}({x}_{-}^{\alpha })-{G}_{\mathrm{intra}}({x}_{+}^{\alpha })]$

2.2. The real part of ${\chi }_{{nn}}^{(0)}$ at zero temperature

The real part of ${\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T)$ at T = 0 can be calculated from the following integral

$$\begin{eqnarray}&&{\rm{Re}}[{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;0)]={N}_{{\rm{f}}}\displaystyle \sum _{\lambda ,{\lambda }^{\prime }}\\ &&\quad \times { \mathcal P }\displaystyle \int \displaystyle \frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}{F}_{\lambda {\lambda }^{\prime }}({\boldsymbol{k}},{\boldsymbol{k}}+{\boldsymbol{q}})\\ &&\quad \times \displaystyle \frac{{n}_{{\rm{F}},\lambda }({\boldsymbol{k}},{\boldsymbol{v}};0)-{n}_{{\rm{F}},{\lambda }^{\prime }}({\boldsymbol{k}}+{\boldsymbol{q}},{\boldsymbol{v}};0)}{{\hslash }\omega +{\varepsilon }_{{\boldsymbol{k}},\lambda }-{\varepsilon }_{{\boldsymbol{k}}+{\boldsymbol{q}},{\lambda }^{\prime }}},\end{eqnarray} \tag{ 12 }$$

where ${ \mathcal P }$ denotes a principal-value integration.

Without loss of generality, we can take ${\boldsymbol{q}}=q\hat{{\boldsymbol{x}}}$ and denote by θ the polar angle of ${\boldsymbol{k}}$. We first compute analytically the integral over $k=| {\boldsymbol{k}}|$ in equation (12) and then treat the angular integration over θ numerically. By combining terms from the sum over λ and ${\lambda }^{\prime }$, the integral in equation (12) can be written as the sum of two contributions,

$$\begin{eqnarray}{\rm{Re}}[{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;0)] & = & {\rm{Re}}[\delta {\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;0)]\\ & & +{\rm{Re}}[{\chi }_{{\rm{u}}}({\boldsymbol{q}},\omega )],\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&\mathrm{Re}[{\chi }_{{\rm{u}}}({\boldsymbol{q}},\omega )]=-\displaystyle \frac{{N}_{{\rm{f}}}{q}^{2}{\rm{\Theta }}({v}_{{\rm{F}}}q-{\hslash }\omega )}{16\sqrt{{({\hslash }{v}_{{\rm{F}}}q)}^{2}-{({\hslash }\omega )}^{2}}},\end{eqnarray} \tag{ 14 }$$

is the undoped (or vacuum) contribution [37] and

$$\begin{eqnarray}&&{\rm{Re}}[\delta {\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;0)]\\ &&\quad ={N}_{{\rm{f}}}\displaystyle \sum _{\alpha ,{\alpha }^{\prime }=\pm 1}{ \mathcal P }\displaystyle \int \displaystyle \frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}{F}_{{\alpha }^{\prime }}({\boldsymbol{k}},{\boldsymbol{k}}+{\boldsymbol{q}})\\ &&\quad \times \displaystyle \frac{\alpha {n}_{{\rm{F}}}({\boldsymbol{k}},\alpha {\boldsymbol{v}};0)}{{\hslash }\omega +\alpha ({\varepsilon }_{{\boldsymbol{k}}}-{\alpha }^{\prime }{\varepsilon }_{{\boldsymbol{k}}+{\boldsymbol{q}}})}\end{eqnarray} \tag{ 15 }$$

is the contribution due to a finite doping. In equation (15) the distribution functions and band energies are evaluated for $\lambda =+$, i.e. ${\varepsilon }_{{\boldsymbol{k}}}={\varepsilon }_{{\boldsymbol{k}},+}$ and ${n}_{{\rm{F}}}({\boldsymbol{k}},{\boldsymbol{v}};0)={n}_{{\rm{F}},+}({\boldsymbol{k}},{\boldsymbol{v}};0)$.

To calculate the integral in equation (15) we observe that, in the absence of a drift velocity ($\beta =0$), the upper limit of the integral over k is the Fermi wave number ${k}_{{\rm{F}}}=\mu /({\hslash }{v}_{{\rm{F}}})$. Indeed, for $\beta =0$ the Fermi–Dirac distribution function in equation (15) has the form:

$$\begin{eqnarray}&&{n}_{{\rm{F}}}({\boldsymbol{k}},0;0)={\rm{\Theta }}(\mu -{\hslash }{v}_{{\rm{F}}}k)={\rm{\Theta }}({k}_{{\rm{F}}}-k),\end{eqnarray} \tag{ 16 }$$

where ${\rm{\Theta }}(x)$ is the Heaviside step function.

For $\beta \gt 0$ the distribution function changes into

$$\begin{eqnarray}{n}_{{\rm{F}}}({\boldsymbol{k}},\alpha {\boldsymbol{v}};0) & = & {\rm{\Theta }}(\mu -{\hslash }{v}_{{\rm{F}}}k(1-\alpha \beta \mathrm{cos}(\psi -\theta ))),\\ & = & {\rm{\Theta }}({k}_{{\rm{F}}}(\theta )-k),\end{eqnarray} \tag{ 17 }$$

where, in analogy with equation (16), we have introduced an angle-dependent Fermi wave number:

$$\begin{eqnarray}&&{k}_{{\rm{F}}}(\theta )=\displaystyle \frac{\mu }{{\hslash }{v}_{{\rm{F}}}[1-\alpha \beta \mathrm{cos}(\psi -\theta )]}.\end{eqnarray} \tag{ 18 }$$

This allows us to rewrite the real part of the response function as

$$\begin{eqnarray}&&{\rm{Re}}[\delta {\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;0)]=D(\mu )\displaystyle \sum _{\alpha =\pm 1}\\ &&\qquad {\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{d\theta }{2\pi }\displaystyle \frac{I(\theta )}{\alpha \bar{\omega }-\bar{q}\mathrm{cos}\theta },\end{eqnarray} \tag{ 19 }$$

where

$$\begin{eqnarray}&&I(\theta )={ \mathcal P }{\displaystyle \int }_{0}^{{\bar{\kappa }}_{{\rm{F}}}(\theta )}{dkk}\displaystyle \frac{k+X(\theta )}{k+Y(\theta )}.\end{eqnarray} \tag{ 20 }$$

In equation (20) we have defined:

$$\begin{eqnarray}&&{\bar{\kappa }}_{{\rm{F}}}(\theta )=\displaystyle \frac{{\hslash }{v}_{{\rm{F}}}{k}_{{\rm{F}}}(\theta )}{\mu },\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&X(\theta )=\displaystyle \frac{\alpha \bar{\omega }+\bar{q}\mathrm{cos}\theta }{2}\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&Y(\theta )=\displaystyle \frac{{\bar{\omega }}^{2}-{\bar{q}}^{2}}{2(\alpha \bar{\omega }-\bar{q}\mathrm{cos}\theta )}.\end{eqnarray} \tag{ 23 }$$

The quadrature in equation (20) can be carried out analytically. The end result is

$$\begin{eqnarray}&&I=\displaystyle \frac{1}{2}{\bar{\kappa }}_{{\rm{F}}}(2X-2Y-{\bar{\kappa }}_{{\rm{F}}})\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&+Y(X-Y)\mathrm{log}\left|\displaystyle \frac{Y}{Y+{\bar{\kappa }}_{{\rm{F}}}}\right|.\end{eqnarray} \tag{ 25 }$$

Having found an expression for $I(\theta )$, the angular integration in equation (19) can be easily carried out numerically.

2.3. Extension to finite temperature

The temperature dependence of the response function in equation (3) is best evaluated by using the following elegant identity due to Maldague [31, 32]:

$$\begin{eqnarray}&&\displaystyle \frac{1}{\mathrm{exp}(x)+1}={\displaystyle \int }_{-\infty }^{\infty }{dy}\displaystyle \frac{{\rm{\Theta }}(y-x)}{4{\mathrm{cosh}}^{2}(y/2)}.\end{eqnarray} \tag{ 26 }$$

Because the Fermi–Dirac distribution at T = 0—see equation (17)—has the form of a Heaviside step function, we can use equation (26) to write the following integral representation:

$$\begin{eqnarray}&&\quad {n}_{{\rm{F}},\lambda }({\boldsymbol{k}},{\boldsymbol{v}};T)\\ &&={\displaystyle \int }_{-\infty }^{\infty }{dy}\displaystyle \frac{{[{n}_{{\rm{F}},\lambda }({\boldsymbol{k}},{\boldsymbol{v}};0)]}_{\mu \to y}}{4{k}_{{\rm{B}}}T{\mathrm{cosh}}^{2}[(y-\mu )/(2{k}_{{\rm{B}}}T)]}.\end{eqnarray} \tag{ 27 }$$

In equation (27) ${n}_{{\rm{F}},\lambda }({\boldsymbol{k}},{\boldsymbol{v}};T)$ is expressed as an integral over μ of the zero-temperature Fermi step, accompanied by an appropriate weighting function. We then observe that the difference between distribution functions in the numerator of equation (3) can be rewritten as:

$$\begin{eqnarray}&&{\displaystyle \int }_{-\infty }^{\infty }{dy}\displaystyle \frac{{[{n}_{{\rm{F}},\lambda }({\boldsymbol{k}},{\boldsymbol{v}};0)]}_{\mu \to y}-{[{n}_{{\rm{F}},{\lambda }^{\prime }}({\boldsymbol{k}}+{\boldsymbol{q}},{\boldsymbol{v}};0)]}_{\mu \to y}}{4{k}_{{\rm{B}}}T{\mathrm{cosh}}^{2}[(y-\mu )/(2{k}_{{\rm{B}}}T)]}.\end{eqnarray} \tag{ 28 }$$

Since μ appears only in the numerator of equation (3), we can integrate first over ${\boldsymbol{k}}$ and then over y obtaining:

$$\begin{eqnarray}&&{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T)\\ &&\quad ={\displaystyle \int }_{-\infty }^{\infty }{dy}\displaystyle \frac{{[{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;0)]}_{\mu \to y}}{4{k}_{{\rm{B}}}T{\mathrm{cosh}}^{2}[(y-\mu )/(2{k}_{{\rm{B}}}T)]}.\end{eqnarray} \tag{ 29 }$$

This formula expresses the finite temperature response function ${\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T)$ in terms of a weighted average of response functions ${\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;0){| }_{\mu \to y}$ at zero temperature and different chemical potentials.

In this section we use the results for the density–density response function reported in section 2 to calculate how the optical absorption spectrum of graphene is changed by a dc current.

Because the drift velocity ${\boldsymbol{v}}=({v}_{x},{v}_{y})$ defines a preferential direction in the graphene sheet, the conductivity ${\sigma }_{{ij}}(\omega )$ stemming from the density-density response function ${\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T)$ in equation (29) transforms like a rank-2 tensor with respect to the Cartesian indices $i,j=x,y$. Following common practice—see Section 3.4 in [32]—in the case of similar rank-2 tensors (like the current-current response function of a homogeneous and isotropic electron system), we can decompose the tensor ${\sigma }_{{ij}}(\omega )$ according to

$$\begin{eqnarray}&&{\sigma }_{{ij}}(\omega )={\sigma }_{{\rm{L}}}(\omega )\displaystyle \frac{{v}_{i}{v}_{j}}{{v}^{2}}+{\sigma }_{{\rm{T}}}(\omega )\left({\delta }_{{ij}}-\displaystyle \frac{{v}_{i}{v}_{j}}{{v}^{2}}\right),\end{eqnarray} \tag{ 30 }$$

where the longitudinal ${\sigma }_{{\rm{L}}}(\omega )$ and transverse ${\sigma }_{{\rm{T}}}(\omega )$ components are defined by:

$$\begin{eqnarray}&&{\sigma }_{{\rm{L}}}(\omega )\equiv {{ie}}^{2}\omega {\left.\underset{q\to 0}{\mathrm{lim}}\displaystyle \frac{{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T)}{{q}^{2}}\right|}_{\psi =0}\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&{\sigma }_{{\rm{T}}}(\omega )\equiv {{ie}}^{2}\omega {\left.\underset{q\to 0}{\mathrm{lim}}\displaystyle \frac{{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T)}{{q}^{2}}\right|}_{\psi =\pi /2}.\end{eqnarray} \tag{ 32 }$$

Of course, ${\sigma }_{{\rm{L}}}(\omega )$ and ${\sigma }_{{\rm{T}}}(\omega )$ depend on v.

The optical absorption ${A}_{{ij}}(\omega )$ is defined by

$$\begin{eqnarray}&&{A}_{{ij}}(\omega )\equiv \displaystyle \frac{4\pi }{c}\mathrm{Re}[{\sigma }_{{ij}}(\omega )].\end{eqnarray} \tag{ 33 }$$

Using equation (30)–(33) we find the following compact expression for the optical absorption:

$$\begin{eqnarray}&&{A}_{{ij}}(\omega )={A}_{{\rm{L}}}(\omega )\displaystyle \frac{{v}_{i}{v}_{j}}{{v}^{2}}+{A}_{{\rm{T}}}(\omega )\left({\delta }_{{ij}}-\displaystyle \frac{{v}_{i}{v}_{j}}{{v}^{2}}\right).\end{eqnarray} \tag{ 34 }$$

Equivalently, we can write equation (34) more explicitly in terms of the angle ψ in equation (8) between the polarization of the incident light and the direction of the current:

$$\begin{eqnarray}&&{A}_{{ij}}(\omega )={A}_{{\rm{L}}}(\omega ){\mathrm{cos}}^{2}(\psi )+{A}_{{\rm{T}}}(\omega ){\mathrm{sin}}^{2}(\psi ).\end{eqnarray} \tag{ 35 }$$

In the T = 0 limit we have

$$\begin{eqnarray}&&{A}_{{\rm{L}}({\rm{T}})}(\omega )\\ &&\quad =\pi {\alpha }_{{\rm{QED}}}{\rm{\Theta }}({\hslash }\omega -2\mu )+{\alpha }_{{\rm{QED}}}\\ &&\qquad {\rm{\Theta }}\left({\hslash }\omega -\displaystyle \frac{2\mu }{1+\beta }\right){\rm{\Theta }}\left(\displaystyle \frac{2\mu }{1-\beta }-{\hslash }\omega \right)\\ &&\qquad \times \left[\pm \displaystyle \frac{{\hslash }\omega -2\mu }{\beta {\hslash }\omega }\sqrt{1-{\left(\displaystyle \frac{{\hslash }\omega -2\mu }{\beta {\hslash }\omega }\right)}^{2}}\right.\\ &&\qquad \left.-\mathrm{sgn}\left(\displaystyle \frac{{\hslash }\omega -2\mu }{\beta {\hslash }\omega }\right)\mathrm{arccos}\left|\displaystyle \frac{{\hslash }\omega -2\mu }{\beta {\hslash }\omega }\right|\right].\end{eqnarray} \tag{ 36 }$$

Equations (35)–(36) are the most important results of this Section. In writing these equations we have restored physical units for the sake of transparency. In the limit $\beta =0$ equation (36) yields ${A}_{{\rm{L}}}(\omega )\;={A}_{{\rm{T}}}(\omega )=\pi {\alpha }_{{\rm{QED}}}{\rm{\Theta }}({\hslash }\omega -2\mu )$.

In figure 1 we show absorption as a function of photon energy ${\hslash }\omega$ (in units of μ) for $\psi =0$ (light polarization parallel to drift velocity) and $\pi /2$ (light polarization perpendicular to drift velocity) and different values of drift velocities (i.e. of the ratio $\beta =v/{v}_{{\rm{F}}}$), at zero and finite temperature. We clearly see that the net effect of the finite drift velocity is to smooth the sharp absorption step $\pi {\alpha }_{{\rm{QED}}}{\rm{\Theta }}({\hslash }\omega -2\mu )$ that occurs at $\beta =0$. The latter stems from Pauli blocking, namely the lack of available empty states in conduction band for vertical inter-band transitions assisted by photons with energy ${\hslash }\omega \lt 2\mu$. Absorption is allowed only for ${\hslash }\omega \gt 2\mu$. Now, for a finite drift velocity β, a window of photon energies opens up between the β-dependent boundaries over which the chemical potential is tilted by the applied dc current:

$$\begin{eqnarray}&&{\mu }^{\pm }\equiv \displaystyle \frac{\mu }{1\mp \beta },\end{eqnarray} \tag{ 37 }$$

with ${\mu }^{+}\gt \mu$ and ${\mu }^{-}\lt \mu$. As illustrated in figure 2, it is possible to excite electron-hole pairs if the photon energy surpasses the lower threshold $2{\mu }^{-}\lt 2\mu$. On the other hand, Pauli blocking extends its action above $2\mu$ in the window $2\mu \lt {\hslash }\omega \lt 2{\mu }^{+}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 3 we show that the way in which the sharp $\beta =0$ absorption step $\pi {\alpha }_{{\rm{QED}}}{\rm{\Theta }}({\hslash }\omega -2\mu )$ is smoothed depends on the angle ψ between in-plane light polarization and drift velocity $\hat{{\boldsymbol{v}}}$ in two major directions, $\psi =0$ and $\psi =\pi /2$. This results in a birefringent absorption in a current-carrying graphene sheet.

Zoom In Zoom Out Reset image size

In order to quantify the degree of absorption anisotropy and its dependence on temperature, we calculate the quantity

$$\begin{eqnarray}&&\delta A(\omega )=\displaystyle \frac{{A}_{{\rm{L}}}(\omega )-{A}_{{\rm{T}}}(\omega )}{{A}_{{\rm{L}}}(\omega )+{A}_{{\rm{T}}}(\omega )}.\end{eqnarray} \tag{ 38 }$$

figures 4(a) and (b) show the temperature and photon frequency dependence of $\delta A(\omega )$. Note that a finite temperature yields a degree of optical anisotropy that is spread over a wider spectral range, compared with the zero-temperature case. For example, figure 3 shows that at ${\hslash }\omega =2.7\;\mu$ the absorption is isotropic at T = 0, while figure 4(b) shows that the anisotropy $\delta A(\omega )$ is still significant at that photon frequency and finite temperature.

Zoom In Zoom Out Reset image size

We emphasize that figures 1 and 4 contain absorption data at finite temperature T, which have been obtained by using the Maldague approach described in section 2.3.

Having calculated the finite-temperature response function in equation (29), one finds the dielectric screening function in the random phase approximation [32]:

$$\begin{eqnarray}&&\varepsilon ({\boldsymbol{q}},\omega ;T)=1-{v}_{{\boldsymbol{q}},\omega }{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T).\end{eqnarray} \tag{ 39 }$$

Here ${v}_{{\boldsymbol{q}},\omega }=2\pi {e}^{2}/[\epsilon (\omega )q]$ is the 2D Fourier transform of the Coulomb potential, with $\epsilon (\omega )$ an effective frequency-dependent environmental dielectric constant and $-e\lt 0$ the electron charge. For the case of graphene on hexagonal boron nitride (hBN) [38, 39], for example, $\epsilon (\omega )=[1+\sqrt{{\epsilon }_{x}(\omega ){\epsilon }_{z}(\omega )}]/2$, while for the case of graphene encapsulated between two hBN slabs $\epsilon (\omega )=\sqrt{{\epsilon }_{x}(\omega ){\epsilon }_{z}(\omega )}$. In writing these equations we have neglected, for the sake of simplicity, finite-thickness effects [39–41] (i.e. Fabry-Pérot phonon-polariton modes).

Plasmons can be found as the roots of the real part of equation (39) on the real-frequency axis, when (Landau) damping is absent or small. Alternatively, they show up as strong absorption peaks in the loss function

$$\begin{eqnarray}&&L({\boldsymbol{q}},\omega ;T)\equiv -{\rm{Im}}\left[\displaystyle \frac{1}{\varepsilon ({\boldsymbol{q}},\omega ;T)}\right],\end{eqnarray} \tag{ 40 }$$

which is directly measured by electron-energy-loss spectroscopy [42]. The loss function displays sharp peaks at the plasmon poles and carries also precious information on inter-band transitions and Landau damping.

4.1. Non-reciprocity of plasmons

In the limit $\bar{q}\lt \bar{\omega }\ll \mu$, the real part of the response function in equation (12) can be expanded as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{Re}}[{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;0)]}{D(\mu )}\approx \displaystyle \frac{W(\beta )}{2\beta }\displaystyle \frac{{\bar{q}}^{2}}{{\bar{\omega }}^{2}}\\ &&\quad -\displaystyle \frac{{\bar{q}}^{2}}{8}+\gamma \displaystyle \frac{W{(\beta )}^{2}}{4\beta }\displaystyle \frac{{\bar{q}}^{3}}{{\bar{\omega }}^{3}}\\ &&\quad +\gamma \displaystyle \frac{\beta }{32}\displaystyle \frac{{\bar{q}}^{3}}{\bar{\omega }}+\displaystyle \frac{{\bar{q}}^{4}}{8\bar{\omega }}+\displaystyle \frac{3}{2}\displaystyle \frac{W{(\beta )}^{2}}{4{\beta }^{2}}\displaystyle \frac{{\bar{q}}^{4}}{{\bar{\omega }}^{4}}+...,\end{eqnarray} \tag{ 41 }$$

where $\gamma =+1\;(-1)$ denotes 'upstream' ('downstream') plasmon propagation, i.e. ${\boldsymbol{q}}\cdot \hat{{\boldsymbol{v}}}/q=+1$ (${\boldsymbol{q}}\cdot \hat{{\boldsymbol{v}}}/q=-1$), and

$$\begin{eqnarray}&&W(\beta )=2\displaystyle \frac{1-\sqrt{1-{\beta }^{2}}}{\beta }.\end{eqnarray} \tag{ 42 }$$

Note that $W(\beta )=\beta$ for $\beta \to 0$. Equation (41) is correct up to ${ \mathcal O }({\bar{q}}^{4})$.

Using this result, it is possible to find an approximate expression for the plasmon dispersion in the presence of current flow and in the long-wavelength $\bar{q}\ll 1$ limit. For the case of a frequency-independent dielectric constant , we find

$$\begin{eqnarray}{\omega }_{{\rm{pl}}}(q) & = & \sqrt{\displaystyle \frac{2{{ \mathcal D }}_{0}W(\beta )}{\epsilon \beta }}\sqrt{q}\left[1+\gamma \displaystyle \frac{\sqrt{2\beta W(\beta )}}{4}\sqrt{\displaystyle \frac{q}{{k}_{{\rm{TF}}{\rm{}}}}}\right.\\ & & +\left.\displaystyle \frac{12-{N}_{{\rm{f}}}^{2}{\alpha }_{{\rm{ee}}}^{2}-3\beta W(\beta )}{16}\displaystyle \frac{q}{{k}_{{\rm{TF}}}}+\ldots \right],\end{eqnarray} \tag{ 43 }$$

where ${{ \mathcal D }}_{0}=4{E}_{{\rm{F}}}{\sigma }_{{\rm{uni}}}/{\hslash }$ is the Drude weight of non-interacting 2D MDFs [4], ${k}_{{\rm{TF}}}={N}_{{\rm{f}}}{\alpha }_{{\rm{ee}}}{k}_{{\rm{F}}}$ is the T = 0 Thomas-Fermi screening wave vector [36], and ${\alpha }_{{\rm{ee}}}={e}^{2}/(\epsilon {\hslash }{v}_{{\rm{F}}})$ is graphene's fine structure constant [4, 36].

Equation (43), which is the most important result of this section, shows that the dc current flow affects the Dirac plasmon dispersion in two ways. On the one hand, it renormalizes the Drude weight, i.e. ${{ \mathcal D }}_{0}\to {{ \mathcal D }}_{0}W(\beta )/\beta$. On the other hand, it adds a linear-in-q subleading term to the Dirac plasmon dispersion, whose sign depends on $\gamma =\pm 1$. Figure 5 shows a comparison between the analytical result (43) (dashed lines) and the plasmon dispersion (solid lines) calculated numerically from the root of equation (39).

Zoom In Zoom Out Reset image size

4.2. Plasmon collimation

The particle-hole excitation spectrum of a current-carrying 2D MDF differs from that in the absence of a current, i.e. for $\beta =0$. A geometrical analysis shows that the region where single-particle excitations cannot exist, i.e. where ${\rm{Im}}[{\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;0)]=0$, is formed by the intersection of two cones in the $({\boldsymbol{q}},\omega \gt 0)$ space:

$$\begin{eqnarray}&&{\hslash }\omega ={\hslash }{v}_{{\rm{F}}}| {\boldsymbol{q}}| ,\end{eqnarray} \tag{ 44 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{2\mu }{1-{\beta }^{2}}-{\hslash }\omega =\left|{\hslash }{v}_{{\rm{F}}}{\boldsymbol{q}}-\displaystyle \frac{2\mu }{1-{\beta }^{2}}{\boldsymbol{\beta }}\right|,\end{eqnarray} \tag{ 45 }$$

where ${\boldsymbol{\beta }}=\beta \hat{{\boldsymbol{v}}}$. Equation (44) defines the usual boundary between the inter- and intra-band electron-hole continua, which is a cone with its apex pointing down. Equation (45) represents another cone with its apex pointing up and shifted from the origin of the $({\boldsymbol{q}},\omega )$ plane by $2\mu {\boldsymbol{\beta }}/(1-{\beta }^{2})$ along the $\hat{{\boldsymbol{v}}}$ direction and $2\mu /(1-{\beta }^{2})$ along the ω axis.

In figure 6 we illustrate the relation between the $\beta =0$ and $\beta \ne 0$ inter- and intra-band continua. In particular, in figures 6(a) and (b) we show the inter- and intra-band continua in the plane $({q}_{\parallel },\omega \gt 0)$ with ${q}_{\parallel }=q\mathrm{cos}(\psi )$, for $\psi =0$ and $\psi =\pi$. Figures 6(c) and (d) show the same quantities for a fixed value of $\omega \gt 0$ and arbitrary $({q}_{\parallel },{q}_{\perp })$, where ${q}_{\perp }=q\mathrm{sin}(\psi )$. The intersection between the two cones (44) and (45) yields the circular arcs in figure 6(d).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 6 we also show plasmon modes obtained from the solution of the equation ${\rm{Re}}[\varepsilon ({\boldsymbol{q}},\omega ;T)]=0$ at T = 0. For the sake of simplicity (to avoid the complications due to the hybridization between plasmons and the hBN optical phonon modes) we have set $\epsilon (\omega )=1$. We clearly see that a finite value of β implies plasmon collimation in a window of angles $-{\psi }_{{\rm{c}}}\leqslant \psi \leqslant {\psi }_{{\rm{c}}}$ around the direction $\hat{{\boldsymbol{v}}}$ of current flow. This is because the lower edge of the inter-band continuum strongly depends on ψ: above a certain frequency threshold—see figure 8—only plasmon modes inside the window $-{\psi }_{{\rm{c}}}\leqslant \psi \leqslant {\psi }_{{\rm{c}}}$ cannot decay by emitting single electron-hole pairs.

Zoom In Zoom Out Reset image size

In figures 7(a) and (c) we show the loss function $L({\boldsymbol{q}},\omega ;T)$ calculated for graphene on hBN at ${\hslash }\omega =140\;{\rm{meV}}$, for a chemical potential $\mu \;=115\;{\rm{meV}}$, and in the absence of a current. The two panels report results for different values of the electron temperature T. The role of an applied current is shown in figures 7(b) and (d), where we present results for the same parameters but for a finite current $J=3\;{\rm{A}}/{\rm{cm}}$. Comparing figure 7(a) with figure 7(b), we clearly see that plasmons that are fully Landau damped at $T=30\;{\rm{K}}$ in the absence of an applied current become sharply defined modes in a relatively narrow range of angles $-{\psi }_{{\rm{c}}}\leqslant \psi \leqslant {\psi }_{{\rm{c}}}$ in the presence of a current. A similar conclusion is also reached at $T=300\;{\rm{K}}$.

The collimation angle ${\psi }_{{\rm{c}}}$, which is defined at T = 0 by the intersection of the solution of ${\rm{Re}}[\varepsilon ({\boldsymbol{q}},\omega ;T)]=0$ (at a fixed energy ${\hslash }\omega$) with the lower threshold of the inter-band electron-hole continuum, is plotted in figure 8 as a function of ${\hslash }\omega /\mu$. The results shown in this figure have been obtained for the simplified case of a frequency-independent effective dielectric constant . In the case of an applied current, plasmons propagating in the range of angles $-{\psi }_{{\rm{c}}}\leqslant \psi \leqslant {\psi }_{{\rm{c}}}$ around the direction of the applied current suffer no inter-band Landau damping. We therefore conclude that an applied dc current allows to collimate the propagation of Dirac plasmons.

Throughout this article we have used the chemical potential μ and the dimensionless drift velocity β as fundamental parameters. For practical purposes it may be useful to have results in terms of the carrier density ${n}_{{\rm{c}}}(\mu ,T,\beta )$ and the current density $J(\mu ,T,\beta )$. The latter quantities are obtained from μ, T, and β as

$$\begin{eqnarray}&&{n}_{{\rm{c}}}(\mu ,T,\beta )=\displaystyle \frac{{n}_{{\rm{c}}}(\mu ,T)}{{(1-{\beta }^{2})}^{3/2}}\end{eqnarray} \tag{ 46 }$$

and

$$\begin{eqnarray}&&{\boldsymbol{J}}(\mu ,T,\beta )=-{{ev}}_{{\rm{F}}}{n}_{{\rm{c}}}(\mu ,T,\beta ){\boldsymbol{\beta }},\end{eqnarray} \tag{ 47 }$$

where ${n}_{{\rm{c}}}(\mu ,T)$ is the carrier density in a graphene sheet at temperature T and chemical potential μ, in the absence of an applied dc current:

$$\begin{eqnarray}&&{n}_{{\rm{c}}}(\mu ,T)={n}_{{\rm{e}}}(\mu ,T)+{n}_{{\rm{h}}}(\mu ,T).\end{eqnarray} \tag{ 48 }$$

Here, ${n}_{{\rm{e}}}(\mu ,T)$ is the density of electrons in conduction band and ${n}_{{\rm{h}}}(\mu ,T)$ is the density of holes in valence band. In this article we consider $\mu \gt 0$ so holes are present only due to thermal excitation. The electron and hole densities in equation (48) are obtained from the $\beta =0$ Fermi–Dirac distribution function (1):

$$\begin{eqnarray}&&{n}_{{\rm{e}}}(\mu ,T)={S}^{-1}\displaystyle \sum _{{\boldsymbol{k}}}{n}_{{\rm{F}},+}({\boldsymbol{k}},0;T)\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{n}_{{\rm{h}}}(\mu ,T)={S}^{-1}\displaystyle \sum _{{\boldsymbol{k}}}[1-{n}_{{\rm{F}},-}({\boldsymbol{k}},0;T)],\end{eqnarray} \tag{ 50 }$$

where S is the sample area. The integrals in equations (49)–(50) can be written as

$$\begin{eqnarray}&&{n}_{\alpha }(\mu ,T)={N}_{{\rm{f}}}\int \displaystyle \frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}{M}_{\alpha }(k,\mu ,T),\end{eqnarray} \tag{ 51 }$$

where $\alpha =\pm 1$ corresponds to electrons (holes) and

$$\begin{eqnarray}&&{M}_{\alpha }(k,\mu ,T)=\displaystyle \frac{1}{\mathrm{exp}\left(\displaystyle \frac{{\hslash }{v}_{{\rm{F}}}k-\alpha \mu }{{k}_{{\rm{B}}}T}\right)+1}.\end{eqnarray} \tag{ 52 }$$

The angular integration in equation (51) yields $2\pi$, while the radial integration gives:

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{dk}}{2\pi }{{kM}}_{\alpha }(k,\mu ,T)\\ &&\quad =-\displaystyle \frac{{({k}_{{\rm{B}}}T)}^{2}}{2\pi {({\hslash }{v}_{{\rm{F}}})}^{2}}{{\rm{Li}}}_{2}[-{e}^{\alpha \mu /({k}_{{\rm{B}}}T)}].\end{eqnarray} \tag{ 53 }$$

Here, ${{\rm{Li}}}_{2}[x]$ is the dilogarithm function. This provides us with the well known expression for the carrier density in a graphene sheet

$$\begin{eqnarray}{n}_{{\rm{c}}}(\mu ,T) & = & \displaystyle \frac{{N}_{{\rm{f}}}{({k}_{{\rm{B}}}T)}^{2}}{2\pi {({\hslash }{v}_{{\rm{F}}})}^{2}}\{-{{\rm{Li}}}_{2}[-{e}^{\mu /({k}_{{\rm{B}}}T)}]\\ & & -{{\rm{Li}}}_{2}[-{e}^{-\mu /({k}_{{\rm{B}}}T)}]\}.\end{eqnarray} \tag{ 54 }$$

To obtain equation (46) we need to calculate equations (49)–(50) with the distribution function (1). We find

$$\begin{eqnarray}&&{n}_{{\rm{c}}}(\mu ,T,\beta )={N}_{{\rm{f}}}\displaystyle \sum _{\alpha =\pm 1}\int \displaystyle \frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}{N}_{\alpha }({\boldsymbol{k}},\mu ,T,\beta ),\end{eqnarray} \tag{ 55 }$$

where

$$\begin{eqnarray}&&{N}_{\alpha }({\boldsymbol{k}},\mu ,T,\beta )\\ &&\quad =\displaystyle \frac{1}{\mathrm{exp}\left\{\displaystyle \frac{{\hslash }{v}_{{\rm{F}}}k[1-\alpha \beta \mathrm{cos}(\theta )]-\alpha \mu }{{k}_{{\rm{B}}}T}\right\}+1}.\end{eqnarray} \tag{ 56 }$$

In equation (56), θ is the angle between ${\boldsymbol{k}}$ and ${\boldsymbol{v}}$. Without loss of generality, the integral in equation (55) can be calculated by choosing ${\boldsymbol{v}}$ to be oriented along the $\hat{{\boldsymbol{x}}}$ direction. To calculate the integral one can replace the radial integration variable k by introducing ${k}^{\prime }=k[1-\alpha \beta \mathrm{cos}(\theta )]$ such that the integral becomes

$$\begin{eqnarray}&&\displaystyle \int \displaystyle \frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}{N}_{\alpha }({\boldsymbol{k}},\mu ,T,\beta )\\ &&={\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{d\theta }{2\pi }\displaystyle \frac{1}{{[1-\alpha \beta \mathrm{cos}(\theta )]}^{2}}\\ &&\qquad \times {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{{dk}}^{\prime }}{2\pi }{k}^{\prime }{M}_{\alpha }({k}^{\prime },\mu ,T).\end{eqnarray} \tag{ 57 }$$

The radial integral in equation (57) is the same as in equation (53), while the angular integration can be carried out by using the residue theorem: for $0\lt \beta \lt 1$ we find

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{d\theta }{2\pi }\displaystyle \frac{1}{{[1-\alpha \beta \mathrm{cos}(\theta )]}^{2}}=\displaystyle \frac{1}{{(1-{\alpha }^{2}{\beta }^{2})}^{3/2}}.\end{eqnarray} \tag{ 58 }$$

The desired result (46) follows from equations (57)–(58).

To obtain equation (47) one proceeds in a similar way, by carrying out an integral over ${\boldsymbol{k}}$ of the expression for the current density in wave-vector space, i.e.

$$\begin{eqnarray}&&{\boldsymbol{J}}(\mu ,T,\beta )=-{{ev}}_{{\rm{F}}}{N}_{{\rm{f}}}\displaystyle \sum _{\alpha =\pm 1}\\ &&\quad \times \displaystyle \int \displaystyle \frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}\displaystyle \frac{{\boldsymbol{k}}}{k}\alpha {N}_{\alpha }({\boldsymbol{k}},\mu ,T,\beta ),\end{eqnarray} \tag{ 59 }$$

where the factor α describes the fact that holes contribute to the current density with the opposite sign with respect to electrons. Similarly to equation (57), we can write

$$\begin{eqnarray}&&\displaystyle \int \displaystyle \frac{{d}^{2}{\boldsymbol{k}}}{{(2\pi )}^{2}}\displaystyle \frac{{\boldsymbol{k}}}{k}{N}_{\alpha }({\boldsymbol{k}},\mu ,T,\beta )\\ &&=\hat{{\boldsymbol{v}}}{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{d\theta }{2\pi }\displaystyle \frac{\alpha \mathrm{cos}(\theta )}{{[1-\alpha \beta \mathrm{cos}(\theta )]}^{2}}\\ &&\qquad \times {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{{dk}}^{\prime }}{2\pi }{k}^{\prime }{M}_{\alpha }({k}^{\prime },\mu ,T).\end{eqnarray} \tag{ 60 }$$

We then notice that

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{d\theta }{2\pi }\displaystyle \frac{\alpha \mathrm{cos}(\theta )}{{[1-\alpha \beta \mathrm{cos}(\theta )]}^{2}}=\displaystyle \frac{\partial {g}_{\alpha }(\beta )}{\partial \beta }\end{eqnarray} \tag{ 61 }$$

with

$$\begin{eqnarray}&&{g}_{\alpha }(\beta )={\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{d\theta }{2\pi }\displaystyle \frac{1}{1-\alpha \beta \mathrm{cos}(\theta )}=\displaystyle \frac{1}{\sqrt{1-{\alpha }^{2}{\beta }^{2}}}.\end{eqnarray} \tag{ 62 }$$

Combining equations (60)–(62) with equation (59) we finally find equation (47).

In figure 9(a) we show the dimensionless drift velocity β as a function of the current density J, for different values of the chemical potential. We clearly see that for low current densities the relation is linear as expected from equation (47). With increasing current density the function $\beta =\beta (J)$ flattens, and $\beta \to 1$ (i.e. the induced drift velocity reaches the graphene Fermi velocity) asymptotically, as expected from equation (47). Finite-temperature effects, shown in figure 9(b), tend to suppress the induced drift velocity.

Zoom In Zoom Out Reset image size

In this article we have presented extensive calculations of the optical and plasmonic properties of a graphene sheet carrying a dc current.

We have calculated in a fully analytical fashion the imaginary part of the density–density linear response function ${\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T)$ of a current-carrying graphene sheet at zero temperature: our main results are reported in section 2.1, equations (7)–(11) and tables 1- 2. The corresponding real part has been calculated by performing the one-dimensional quadrature in equation (19). Results at finite temperature have been obtained by utilizing the Maldague identity in equation (29). We emphasize that the latter results are not of mere academic interest, since the electron temperature can be significantly raised above the lattice temperature under the application of a large dc current.

Table 2.  Same as in table 1, but for the function ${H}_{{\rm{inter}}}(x,y,z)$. Here, ${G}_{\mathrm{inter}}(x)=x\sqrt{1-{x}^{2}}+\mathrm{arcsin}(x)$ and y_i with $i=\pm 1$ denote the roots of the function ${g}_{\mathrm{inter}}(y)={A}_{-}+{By}-{C}_{-}\sqrt{1-{y}^{2}}$. As in table 1, the subscript i is such that ${g}_{\mathrm{inter}}^{\prime }({y}_{+})\gt 0$, while ${g}_{\mathrm{inter}}^{\prime }({y}_{-})\lt 0$.

	${A}_{-}+B\gt 0$	${A}_{-}+B\lt 0$
${A}_{-}-B\gt 0$	$\begin{array}{c}{\rm{\Theta }}[{A}_{-}-\mathrm{sgn}({C}_{-})\sqrt{{B}^{2}+{C}_{-}^{2}}][{G}_{\mathrm{inter}}({y}_{-})-{G}_{\mathrm{inter}}({y}_{+})+\pi ]\\ +{\rm{\Theta }}[\mathrm{sgn}[{C}_{-}]\sqrt{{B}^{2}+{C}_{-}^{2}}-{A}_{-}]\end{array}$	${G}_{\mathrm{inter}}({y}_{-})+\pi /2$
${A}_{-}-B\lt 0$	$\pi /2-{G}_{\mathrm{inter}}({y}_{+})$	${\rm{\Theta }}[{A}_{-}-\mathrm{sgn}({C}_{-})\sqrt{{B}^{2}+{C}_{-}^{2}}][{G}_{\mathrm{inter}}({y}_{-})-{G}_{\mathrm{inter}}({y}_{+})]$

This accurate knowledge of the density-density response function of a current-carrying graphene sheet has allowed us to calculate the optical absorption spectrum ${A}_{{ij}}(\omega )$, which displays a birefringent character with respect to the angle ψ between the in-plane light polarization and current flow: see equation (36) and figures 1- 3.

Finally, we have used the non-interacting density-density response function ${\chi }_{{nn}}^{(0)}({\boldsymbol{q}},\omega ;T)$ to calculate the dielectric screening function in the random phase approximation, the plasmon dispersion, and the loss function of a current-carrying graphene sheet. We have showed that graphene plasmons acquire a degree of non-reciprocity under the influence of a dc current, see equation (43), and that their propagation can be collimated in a window of angles around the current direction, see figures 6- 8.

Before concluding, we would like to briefly discuss the possibility to test our predictions experimentally. To observe the discussed effects one needs to pass a large current through a graphene flake and measure polarized optical reflection or transmission through the sample. However, the large current passing through graphene would imply large Joule losses. Hence, an efficient thermal sink is needed in order to keep the temperature of graphene in equilibrium with the thermal bath and, therefore, constant. This implies that we cannot use the same set-up that was used previously to show that graphene membranes can host large current densities. We have to immobilize the graphene flake on an inert and thermally conductive substrate that acts as the thermal sink and need to perform the measurements in an inert gas atmosphere or vacuum environment. In addition, because low temperatures are beneficial for the observation of the effect, the usage of a suitable optical cryostat is needed. Currently, experiments to detect the current induced birefringence in graphene are under way.

This work was supported by the EC under the Graphene Flagship program (contract no. CNECT-ICT-604391) and MIUR through the program 'Progetti Premiali 2012' - Project 'ABNANOTECH'. B.V.D. wishes to thank the Scuola Normale Superiore (Pisa, Italy) for the kind hospitality while this work was carried out and Research Foundation Flanders (FWO-Vl) for a PhD Fellowship.