Optical conductivity of ABA stacked graphene trilayer: mid-IR resonance due to band nesting

Since the isolation of monolayer graphene almost a decade ago [1], there has been a high interest in the low energy transport and optical properties of not only the monolayer [2, 3] but also a few layer graphene systems [4–8]. These properties are determined by the electronic band structure near the $\mathbf K$ point. The undoped monolayer graphene (MLG) is characterized by the universal optical conductivity, σ₀ = e²/(4ℏ). This implies that the transmittance depends solely on the fine structure constant and originates the quantized visible opacity of suspended monolayer graphene [9, 10]. As far as doped graphene is concerned, there are several effects that arise in the optical properties, related to the restrictions introduced on the interband transitions by the state filling and also to the onset of intraband transitions [11]. The latter correspond to plasmons and give rise to the interesting and promising field of graphene plasmonics [3, 12].

Graphene multilayers offer a new ingredient to the interesting physics and potential applications. The relatively weak interlayer coupling, on the one hand, implies that they should inherit some properties of the parent material [4], on the other hand, it introduces a new energy scale, of the order of few tenths of the electron-volt, that should yield some new properties. It has been shown [13] that there is also a universal optical conductivity in a undoped N-layer graphene, equal to σ_0N = Nσ₀, that is reached in undoped graphene in the low frequency limit. At the same time, the optical response of doped bilayer graphene reveals intense strongly doping-dependent features in the mid-infrared (around 0.4 eV) [5, 6]. The origin of these experimentally observed features has been considered theoretically in these works and also, in more detail, in [14], where the band structure and the optical conductivity of bilayer graphene were calculated. Some novel plasmonic effects in Bernal-stacked bilayer graphene were predicted in the recent work [15]. As the number of layers increases beyond N = 2, the band structure and the optical conductivity become dependent not only on N but also upon the stacking arrangement. There are three distinct planar projections of the honeycomb lattice (usually denoted A, B and C) and, consequently, 2^N−2 distinct N−layer sequences [13]. In particular, the stacking of three layers in a graphene trilayer can be either ABA (also called Bernal-type) and ABC (rhombohedral) [16]. These two different stacking arrangements lead to strikingly different electronic band strucures [17]. For instance, it was found that undoped graphene ABC trilayer shows many-body correlations with an energy gap, while the Bernal-type stacking (taking place in graphite) does not lead to a gap [18]. Non-Bernal stacked multiple graphene layers have attracted a considerable attention related to the prospect of further enhancement of capabilities of graphene-based optoelectronic devices, in particular, THz and IR photodiodes [19]. The optical conductivity of ABC trilayers has been considered in a number of works [20–22], in particular, the effect of doping has been analysed [22]. Even though trilayer graphene contains, on average, regions of ABC and ABA stacking in an 15 : 85 ratio [23], apparently the optoelectronic properties of the latter attracted less attention and we are aware of only one work [24] devoted to this topic. This is notwithstanding the possibility of using far-infrared (FIR) spectroscopy, along with the common Raman scattering technique [25] in order to distinguish different trilayer graphene species [23, 24].

Therefore, the main purpose of this article is to analyse the spectral characteristics of the optical conductivity of intrinsic and doped ABA stacked trilayers. The mirror-symmetric Bernal stacking is the most common in graphene multilayers that can be exfoliated from natural graphite since it shares its crystalline structure [7]. Intuitively, one can expect that the electrons in ABA trilayers can have a monolayer- or bilayer-like character [26]. We shall present the analytical dispersion relation for the electrons near the Dirac point, from which it follows that the band structure of the ABA trilayer indeed looks like a superposition of those characteristic of a monolayer and a bilayer, although the latter corresponds to an effective interlayer hopping constant $\sqrt {2}$ times larger than the true one. Based on this band structure, we calculate both the interband and intraband (Drude) contributions to the optical conductivity. The most interesting result is that the optical absorption of a doped ABA trilayer is dominated by a narrow resonant peak at the frequency corresponding to this effective interlayer hopping energy (approximately 0.56 eV). It is caused by the fact that the dispersion curves corresponding to two bands are nearly parallel for a considerable range of wavevectors near the Dirac point. This effect sometimes is called 'band nesting' [27] in order to distinguish from van Hove singularities in the single-particle density of states. We will show that the intensity of this resonant absorption is approximately proportional to the Fermi energy and, therefore, can be controlled by adjusting the gate voltage applied to the graphene layer.

2.1. Band structure

The tight-binding Hamiltonian for non-interacting electrons in the ABA stacked trilayer involves three A-type and three B-type sites and includes the essential in-plane (t₀ ≈ 2.7 eV) and interlayer hoppings (t₁ ≈ 0.4 eV) as shown in figure 1. These two parameters, connecting atoms that are right on top of each other in adjacent layers, are sufficient to describe the main features of the band structure, such as the type of dispersion of the energy bands and their separation, as confirmed by recent DFT calculations [17]. Keeping only t₀ and t₁ hoppings permits to obtain simple formulae for band gaps, effective masses and the Fermi velocity.

Zoom In Zoom Out Reset image size

Then the Hamiltonian is given by [18]:

$$\begin{eqnarray} &\fl \hat{H}\,{=}\,{-}t_{0}\sum_{n,\delta_{j}}\vert A_{1},\vec R_{n}+\vec \delta_{j}\rangle \langle B_{1}, \vec R_{n} \vert \nonumber\\ &\quad -t_{0}\sum_{n,\vec \delta_{j}} \vert A_{3},\vec R_{n}+\vec \delta_{j}\rangle \langle B_{3},\vec R_{n}\vert\nonumber -t_{0}\sum_{n,\vec \delta_{j}}\vert A_{2},\vec R_{n}\rangle \nonumber\\ &\quad \times\,\langle B_{2},\vec R_{n}-\vec \delta_{j}\vert +\,t_1\sum_{n}\vert A_{2},\vec R_{n}\rangle \langle B_{1},\vec R_{n}\vert \nonumber\\ &\quad +\,t_1\sum_{n}\vert A_{2},\vec R_{n}\rangle \langle B_{3},\vec R_{n}\vert \:+ \: \mathbf {H. c.} {Ham} \end{eqnarray} \tag{ 1 }$$

Here n = (n₁, n₂) is the composite index, which determines the atomic positions in the lattice, $\vec R_{n}=n_1\vec g_1+n_2\vec g_2$, with $\vec g_1=(\frac{\sqrt {3}}{2},\frac{3}{2})a_0$, $\vec g_2=(-\frac{\sqrt {3}}{2},\frac{3}{2})a_0$ being the lattice vectors. The three vectors that connect the B atom to its three nearest neighbors are $\vec{\delta}_{1}=(\frac{\sqrt {3}}{2},-\frac{1}{2})a_0$, $\vec{\delta}_{2}=(-\frac{\sqrt {3}}{2},-\frac{1}{2})a_0$, and $\vec{\delta}_{3}=(0,1)a_0$, where a₀ is the C–C interatomic distance. The positions of the A atoms relative to the B atoms in each of the three layers are shown in figure 1.

The energy spectrum of the Hamiltonian (A.4) is composed of six energy bands given by (see appendix A for details):

$$\begin{eqnarray} &&E_{\pm1}\left(\vec k\right)=\pm\sqrt {t_1^2+\vert \phi\left(\vec k\right)\vert^2- t_1\sqrt {t_1^2+2\vert \phi\left(\vec k\right)\vert^2}}\:; \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&E_{\pm2}\left(\vec k\right)=\pm\sqrt {\vert \phi\left(\vec k\right)\vert^2}\:; \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&E_{\pm3}\left(\vec k\right)=\pm\sqrt {t_1^2+\vert \phi\left(\vec k\right)\vert^2+ t_1\sqrt {t_1^2+2\vert \phi\left(\vec k\right)\vert^2}}\:, \end{eqnarray} \tag{ 4 }$$

where

$$\begin{equation} \phi(\vec{k})=-t_{0}\left[\!\exp\left({\rm i}k_{y}a_{0}\right)+2\exp\left(\frac{-{\rm i}k_{y}a_{0}}{2}\right)\cos \left(\!\frac{k_{x}\sqrt 3}{2}a_{0}\!\right)\!\right]. {eq:phi} \end{equation} \tag{ 5 }$$

The band structure is depicted in figure 2. Note that the gap between the bands ±3 and the Dirac point is ${2\it{\Delta}=2\sqrt 2}t_1$, while its counterpart in graphene bilayer is just equal to 2t₁ [14].

Zoom In Zoom Out Reset image size

In the vicinity of $\mathbf K_{\pm}$ points [Dirac points, $\vec k\,{=}\,{\pm}(4\pi/$ $(3\sqrt {3}a_0),0) ]$ the band energies can be approximated in the following way:

$$\begin{eqnarray} &&E_{\pm1}=\pm\frac{\hbar^2K^2}{2m_{1}}; \ E_{\pm2}=\pm\hbar \acute{v_{f}}K; \ E_{\pm3}=\pm\left(\frac{\hbar^2K^2}{2m_{3}}+\it{\Delta}\right), \nonumber\\ \end{eqnarray} \tag{ 6 }$$

where $K=\sqrt {k_y^2+k_y^2}$. Note that the effective masses at the bottom of the bands 1 and 3 are equal within the present model, $m_1=m_3=t_1/(\sqrt 2 v_f^2)$, and the Fermi velocity of the Dirac-type band coincides with that of monolayer graphene, $\acute{v_{f}}=3t_{0}a_{0}/(2\hbar)=v_f$.

2.2. Optical conductivity: interband part

The finite frequency (optical) conductivity is calculated through a standard procedure using the Kubo formula [29],

$$\begin{eqnarray} &\fl \sigma (\omega)=\frac{-2{\rm i}e^2}{\omega S}\sum_{\vec k,l^\prime\neq l}\left|\left\langle \vec k,l\right|\hat{V}_x\left|\vec k,l^\prime\right\rangle\right|^2 \nonumber\\ &\quad \times \frac{n_{f}\left[E_{l}(\vec k)\right]-n_{f}\left[E_{l^\prime}(\vec k)\right]}{\hbar\omega-E_{l}(\vec k)+E_{l^\prime}(\vec k)+{\rm i}\Gamma}, \end{eqnarray} \tag{ 7 }$$

where e is the electron charge, $S=\frac{3\sqrt {3}a_{0}^{2}}{2}$ is the unit cell area, n_f is the Fermi function, $\left|\vec k,l\right\rangle$ is the eignestate corresponding to energy $E_{l}(\vec k)$ (l, l' are the band indices), Γ is a damping parameter and $\hat{V}_x$ denotes the x-component of in-plane velocity operator, defined by

$$\begin{eqnarray} \nonumber &\fl \hat{\vec{V}}=\frac{{\rm i}t_{0}}{\hbar}\sum_{\delta_{j}}\vec{\}delta_{j}\left[\left|B_{1},\vec R_{n}\right\rangle\left\langle A_{1},\vec R_{n}+\vec \delta_{j}\right|\right.\nonumber\\ &-\left|A_{1},\vec R_{n}+\vec \delta_{j}\right\rangle\left\langle B_{1},\vec R_{n}\right|+\left|B_{2},\vec R_{n}-\vec{\}delta_{j}\right\rangle\left\langle A_{2},\vec R_{n}\right|\nonumber\\ &-\left|A_{2},\vec R_{n}\right\rangle\left\langle B_{2},\vec R_{n}-\vec \delta_{j}\right|+\left|B_{3},\vec R_{n}\right\rangle\left\langle A_{3},\vec R_{n}+\vec \delta_{j}\right| \nonumber\\ &\left.-\left|A_{3},\vec R_{n}+\vec \delta_{j}\right\rangle\left\langle B_{3},\vec R_{n}\right|\right]. \end{eqnarray} \tag{ 8 }$$

The summation over $\vec k$ in equation (7) involves wavevectors in the first Brilluion zone, which is the usual hexagon with side $\frac{4\pi}{(3\sqrt 3 a_{0})}$. In practice, this summation is replaced by a 2D integration over the triangle formed by the points $\mathbf \Gamma$ (0; 0), $\mathbf K$ $\left(\frac{4\pi}{(3\sqrt 3 a_{0})};\:0\right)$, and $\mathbf M$ $\left(\frac{\pi}{(3 a_{0})};\:\frac{\pi}{(\sqrt 3 a_{0})}\right)$, and the result is multiplied by twelve.

2.3. Optical conductivity: Drude part

For doped graphene (μ ≠ 0) it is necessary to include in the optical conductivity also the term related to intraband transitions, often referred to as Drude term. Although it can be calculated through the Kubo formula, it is easier to derive this term using the Boltzmann transport equation [12, 28]. Thus, the electric current is written as

$$\begin{equation} \vec{J}=\frac{4e}{(2\pi)^{2}}\sum_{l=1}^3\int {\rm d}\vec{k}\,\,\delta n_{\vec{k}}^{(l)}\,\vec{v}_l(\vec{k})\,, \end{equation} \tag{ 9 }$$

where $\delta n_{\vec{k}}^{(l)}$ is the deviation of the carriers distribution from the equilibrium Fermi–Dirac function, $n_{f}[E_l(\vec k)]$. The former is readily obtained from the Boltzmann equation,

$$\begin{equation} \delta n_{\vec{k}}^{(l)}=\frac{e\vec{\mathcal E}\cdot \vec{v}_l(\vec{k})}{{\tau(\vec{k})}^{-1}-{\rm i}\omega}\left (-\frac {\partial n_{f}[E_l(\vec k)]}{\partial E_l(\vec k)} \right )\:, \end{equation} \tag{ 10 }$$

where $\vec{v}_l(\vec{k})$ is the group velocity of the charge carriers in the l-th band and $\tau(\vec{k})$ denotes the carrier relaxation time. The factor of 4 in equation (9) is due to the spin and valley degeneracy.

Substitution of equation (10) into equation (9) yields the Drude conductivity and for zero temperature we have:

$$\begin{eqnarray} &&\sigma _D =\frac{4\sigma_0}{\pi}\left [\frac{3\mu}{\hbar(\gamma-{\rm i}\omega)} +\frac{2(\mu-\it{\Delta})}{\hbar(\Gamma-{\rm i}\omega)}\theta(\mu-\it{\Delta})\right ]\:, {Drude} \end{eqnarray} \tag{ 11 }$$

where the damping parameters γ and Γ are defined as the inverse of the corresponding relaxation time at the Fermi level (we make no distinction between bands 1 and 2). This Drude term has to be added to the optical conductivity (7). Note that the first term is precisely three times the Drude conductivity of monolayer graphene.

We shall now concentrate on the frequency dependence of the real part of the derived optical conductivity, σ'(ω), that determines the absorption, for different values of the chemical potential, μ. The spectra of real, σ'(ω), and imaginary, σ^''(ω), parts of the conductivity for three different values of μ (conveniently expressed in units of t₁) are presented in figures 3 and 4, respectively. Note that the first value corresponds to undoped graphene (μ = 0), the second one is μ = 0.2t₁ < Δ and the last μ > Δ. For intrinsic trilayer graphene, either ABA or ABC stacked, σ'(ω) tends to 3σ₀ ([22, 24]). We clearly see the effects of doping present as the Fermi step at ℏω = 2μ (for μ > 0), known in semiconductors as Burstein–Moss effect [31] and familiar in monolayer graphene [28]. The feature characteristic of trilayer graphene, located at ℏω ⩾ Δ, splits into two for μ > 0 (compare red and blue curves to the black one in figure 3). But the most impressive effect of the doping is the onset of the resonant peak at ℏω = Δ, whose intensity grows strongly with μ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to understand these spectral changes caused by doping, we analyzed all possible optical transitions listed in table 1 (note that none of the velocity matrix elements vanishes, they are all allowed!). As it can be seen from figure 2, there are two Dirac-type and four approximately parabolic bands that arise from six atoms in the unit cell of the ABA trilayer. Let us consider the allowed transitions that correspond to the features of the optical conductivity. For the undoped graphene (μ = 0), we have nine possible transitions, with energy conservation restrictions imposed on some of them (see table 1). There are four transitions allowed without any restriction for all frequencies, namely, l' = −1 to l = 1 (denoted as −1 → 1, −1 → 2, −2 → 1 and −2 → 2. The onset of transitons involving bands ±3 is at $\hbar\omega>\sqrt 2t_1$ and $\hbar\omega>2\sqrt 2t_1$. That is why there are four possible transitions including −1 → 3, −2 → 3, −3 → 1 and −3 → 2, with a threshold frequency $\hbar\omega=\sqrt 2{t_1}$. Still another possible transition is from −3 to 3, with a threshold at $\hbar\omega=2\sqrt 2t_1$, although it is less pronounced as clearly seen in figure 3 (bold black curve).³ With doping the system so as μ = 0.2t₁, the behavior of the optical conductivity is altered, caused by the changes in the allowed transitions with applying more state–filling restrictions as well as the onset of new transitions.

Table 1. Optical transitions and threshold frequencies in units of t₁.

Transition	μ/t1 = 0	μ/t1 = 0.2	μ/t1 = 1.6	Equation
−1 → 1	0	2μ/t1	2μ/t1	(12)
−1 → 2	0	0.22	2.64	(12)
−1 → 3	$\sqrt 2$	$\sqrt 2$	$2\mu/t_1-\sqrt 2=1.8$	(13)
−2 → 1	0	0.77	3.8	(14)
−2 → 2	0	2μ/t1 = 0.4	2μ/t1 = 3.2	(15)
−2 → 3	$\sqrt 2$	$\sqrt 2$	2.15	(15)
−3 → 1	$\sqrt 2$	1.8	4.8	(16)
−3 → 2	$\sqrt 2$	1.64	4.05	(16)
−3 → 3	2 $\sqrt {2}$	2 $\sqrt {2}$	2μ/t1 = 3.2	(17)
1 → 2	Not possible	0.17	0.56	(18)
1 → 3	Not possible	$\sqrt 2$	$\sqrt 2$	(18)
2 → 3	Not possible	$\sqrt 2$	0.86	(19)

The general expression for the threshold energy for the l' → l transition (for the case l' < 0, l > 0) is given by:

$$\begin{eqnarray*}&&\nonumber \it{\Delta}_{l^\prime \rightarrow l}=E_{l}\left(K_{f}^{(l)}\right)-E_{l^\prime}\left(K_{f}^{(l)}\right), \end{eqnarray*}$$

where $K_{f}^{(l)}$ is the Fermi wavevector in the l-th band (the root of equation $E_{l}\left(K_{f}^{(l)}\right)=\mu$, or $K_{f}^{(l)}=0$, if the equation does not have a real root). At the same time, for the case where l', l > 0 the expression for the threshold energy is given by

$$\begin{eqnarray*}&\fl\it{\Delta}_{l^\prime \rightarrow l}={\rm min}\left[E_{l}\left(K_{f}^{(l)}\right)-E_{l^\prime}\left(K_{f}^{(l)}\right),\:E_{l}\left(K_{f}^{(l^\prime)}\right)\right. \nonumber \\ &\left.-E_{l^\prime}\left(K_{f}^{(l^\prime)}\right)\right]. \nonumber \end{eqnarray*}$$

There is also an upper cut-off frequency in this case, given by

$$\begin{eqnarray*}&\fl\Omega_{l^\prime \rightarrow l}={\rm max}\left[E_{l}\left(K_{f}^{(l)}\right)-E_{l^\prime}\left(K_{f}^{(l)}\right),\:E_{l}\left(K_{f}^{(l^\prime)}\right)\right. \nonumber \\ &\left.-E_{l^\prime}\left(K_{f}^{(l^\prime)}\right)\right], \nonumber \end{eqnarray*}$$

so that the allowed transition frequencies lie inside the domain Δ_l' → l ⩽ ℏω ⩽ Ω_l' → l. It is possible to obtain simple expressions for the threshold energies using the approximation (6). Thus, we obtain the following values:

$$\begin{eqnarray} &&\it{\Delta}_{-1 \rightarrow 1}=2\mu,\ \ \it{\Delta}_{-1 \rightarrow 2}=\mu+\frac{\mu^2}{\sqrt 2t_1}, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&\it{\Delta}_{-1 \rightarrow 3}=2\left(\mu-\sqrt 2t_1\right)\theta(\mu-\sqrt 2t_1)+\sqrt 2t_1, \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&\it{\Delta}_{-2 \rightarrow 1}=\mu+\sqrt {\sqrt 2t_1\mu}, \ \ \it{\Delta}_{-2 \rightarrow 2}=2\mu, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &\fl \it{\Delta}_{-2 \rightarrow 3}=\sqrt 2t_1+\theta(\mu-\sqrt 2t_1) \nonumber\\ &\qquad \times \left(\mu-\sqrt 2t_1+\sqrt {\sqrt 2t_1\mu-2t_1^2}\right), \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&\it{\Delta}_{-3 \rightarrow 1}=2\mu+\sqrt 2{t_1},\ \ \it{\Delta}_{-3 \rightarrow 2}=\mu+\sqrt 2{t_1}+\frac{\mu^2}{\sqrt 2{t_1}}, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&\it{\Delta}_{-3 \rightarrow 3}=2\left(\mu-\sqrt 2{t_1}\right)\theta(\mu-\sqrt 2t_1)+2\sqrt 2{t_1}. \end{eqnarray} \tag{ 17 }$$

Several transitions can occur only for nonzero μ:

$$\begin{eqnarray} &&\it{\Delta}_{1 \rightarrow 2}=\mu-\frac{\mu^2}{\sqrt 2t_1},\ \ \it{\Delta}_{1 \rightarrow 3}=\sqrt 2, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} &&\it{\Delta}_{2 \rightarrow 3}=\sqrt 2t_1+\frac{\mu^2}{\sqrt 2{t_1}}-\mu. \end{eqnarray} \tag{ 19 }$$

The numerical values of the threshold energies (in units of t₁) are given in table 1 where the numerical values correspond to the exact band structure (see appendix B) because the parabolic/linear approximation becomes insufficient for higher values of μ. It becomes clear from table 1 that almost only the 1 → 3 and 2 → 3 transitions (not possible for μ = 0) are responsible for the onset of the resonant peak at ℏω = Δ. We calculated these contributions to the optical conductivity and found that the intensity of the 1 → 3 transition exceeds by far that of the other one. It can be understood by the fact that the dispersion curves E₁(K) and E₃(K) are nearly parallel for a broad range of K values (band nesting) and therefore the joint density of states for this transition is large, as can be seen clearly from figure 5. This figure takes into account only the 1 → 3 transition and shows that it is responsible for the intensity of the resonant absorption band. The increasing intensity of the 1 → 3 transition as μ grows is related to the growing number of occupied states that are depopulated by absorbing electromagnetic radiation. This is the principal spectral feature observed in [24] (where it was denoted 'band C'), which can be considered as characteristic of ABA graphene [23, 25] and we return to it below. For heavily doped p−type layers studied in [24], two weak absorption bands (denoted there as A and B bands) were observed at lower frequencies and attributed to the −3 → −2 and −2 → −1 transitions. In our consideration of n−type layers these transitions correspond to 1 → 2 and 2 → 3, respectively, both involving the Dirac-type band. Indeed, they are allowed (even though only for doped samples) but, according to our results, their intensity is rather low (for instance, we cannot see any feature at ω/t₁ ≈ 0.86, the threshold frequency of 2 → 3 transition for μ/t₁ = 1.6 (see table 1). Possibly some further effects can enhance transitions to and from the Dirac-type electronic bands.

Zoom In Zoom Out Reset image size

The approximately triangular shaped feature in the spectrum for μ = 0 (ℏω ⩾ Δ) (that could be anticipated intuitively [26]) splits into two sub-bands corresponding to the −1 → 3 and −3 → 1 transitions. At μ = 0, they have the same energy (= Δ) but for μ > 0 the former is shifted to higher energy (see table 1). Note that, by chance, in figures 3 and 4 the features related to the −3 → 1 transition for μ/t₁ = 0.2 and to the −1 → 3 one for μ/t₁ = 1.6 appear nearly at the same frequency, ω/t₁ ≈ 1.8. Such a feature was not observed in [24] (whose experimental situation qualitatively corresponds to our case of μ/t₁ = 0.2), possibly because it was hindered by Fabri-Perot interference in the substrate.

The intensity of the 1 → 3 band as function of μ can be evaluated analytically if we assume that the matrix element does not depend on K:

$$\begin{eqnarray*}\fl I(\omega)=&\frac{\rm const}{\omega^2S}\sum_{\vec{k}}\delta\left\{E_3(\vec k)-E_1(\vec k)-\hbar\omega\right\}\nonumber \\ &\times\left\{n_f \left[E_1(\vec k)\right]-n_f \left[E_3(\vec k)\right]\right\}. \nonumber \end{eqnarray*}$$

Changing from sum to integration and using the parabolic approximation (6) near the Dirac point, in the limit of low energies for $E_1(\vec k)$ and $E_3(\vec k)$ we obtain:

$$\begin{eqnarray} &\fl\;\; I(\omega) ={\rm const}\cdot \frac{m_1}{2\pi(\hbar\omega)^2}{\mathcal L} (\hbar \omega-\it{\Delta})\nonumber\\ & \times\!\!\int_0^\mu\!\!{\rm d}E_1\left\{\theta\left[\mu-E_1(\vec k)\right]- \theta\left[\mu-E_1(\vec k)- \hbar\omega\right]\right\}\!. \end{eqnarray} \tag{ 20 }$$

Here ${\mathcal L} (\hbar \omega-\it{\Delta})$ denotes a Lorentzian replacing the δ− function broadened because of natural reasons. The integral in equation (20) depends on the value of μ, that is, I(ω) ∼ μ for μ < Δ and I(ω) = const for μ > Δ. This analytical result is compared to the numerical data in figure 6.

Zoom In Zoom Out Reset image size

Higher order hopping parameters, in particular, those connecting more distant atoms in adjacent layers and usually denoted γ₃ and γ₄ [5, 17], are known to cause a distortion (trigonal warping) of the low-energy bands and a small electron-hole asymmetry [17]. Therefore they can slightly affect the intensity and shape of the resonant absorption band but its position, to a good approximation, is determined by the t₁ hopping [17]. This is confirmed by the results of our ab-initio DFT calculations for intrinsic graphene, shown in the inset of figure 3.

In summary, we calculated the spectral dependence of the real and imaginary parts of the optical conductivity of Bernal-stacked trilayer graphene. Even though the energy bands of this material look like a superposition of those of a monolayer and a bilayer (with a larger gap between the Dirac point and higher parabolic band), the optical spectra are rich and interesting, especially in the case of gated (doped) graphene because all kinds of interband transitions are allowed. In particular, there is a strong and narrow resonant band cause by transitions between two 'nested' parabolic bands. The intensity of this band is controlled by the Fermi level position and attains a maximum for μ ≈ 1.4 t₁ (figure 6), providing a strong light-matter coupling in the atomically thin layer [32]. Although the physical origin of the resonant absorption is clear and it could be predicted by just inspecting the band structure of the ABA stacked graphene, its dependence on the Fermi level is not evident without calculations. It would be interesting to show experimentally that the intensity of the resonant absorption attains its maximum for the Fermi energy of μ ≈ 0.6 eV. This can be achieved by combining usual doping and application of gate voltage [33]. The application of a large gate voltage can slightly modify the band structure but it can be taken into account by adding appropriate (unequal) constant potentials in the Hamiltonian (1) for atoms belonging to the different monolayers [24, 34]. The modulation of the peak absorption can be interesting for optoelectronic devices such as resonant photodetectors with adjustable sensitivity and it is robust with respect to the temperature. We performed calculations also for T = 300 K and the spectra are very similar to those presented in the figures, except for the obvious broadening of the Fermi steps seen in figure 3, i.e. the resonant band remains narrow. Compared to bilayer graphene, the resonant band occurs at a (1.4 times) higher frequency, which broadens the range of applications of Bernal-stacked multilayer graphene materials that can be directly exfoliated from natural graphite. Also, we would like to point out that the Drude conductivity is three times higher than for monolayer graphene with same doping level, accordingly, the surface plasmon frequency would be also higher (by a factor of $\sqrt 3$) [12], extending the spectral range of possible applications of graphene plasmonics.

Financial support from the COMPETE Programme (FEDER) and the Portuguese Foundation for Science and Technology (FCT) through Projects PEst-C/FIS/UI0607/2013 and PTDC/FIS/113199/2009, as well as MATEPRO Project (ON2 Program) is gratefully acknowledged. ZR thanks the hospitality of the Physics Center of Minho University during her stay in Portugal. YVB, NMRP, RR and MIV acknowledge support by the EC under Graphene Flagship (contract no. CNECT-ICT-604391).

After transforming the Hamiltonian (1) to the momentum space (by representing $\vert A_{m},\vec R_{n}\rangle=\frac{1}{\sqrt {N}}\sum_{k}\psi_{A_{m},k}\exp\left(-{\rm i}\vec{k}\cdot\vec R_{n}\right)$, $\vert B_{m},\vec R_{n}\rangle=\frac{1}{\sqrt {N}}\sum_{k}\psi_{B_{m},k}\exp\left(-{\rm i}\vec{k}\cdot\vec R_{n}\right)$, m = 1, 2, 3 and N is the number of unit cells), its matrix form is:

$$\begin{eqnarray*}\hat{H}= \left( \begin{array}{@{}cccccc@{}} 0 & \phi& 0 & 0 & 0 & 0 \\ \phi^{*}& 0 & t_1 & 0 & 0 & 0 \\ 0 & t_1 &0 & \phi & 0 & t_1 \\ 0& 0& \phi^{*}& 0 & 0 & 0 \\ 0& 0& 0& 0 & 0 & \phi \\ 0& 0& t_1& 0 & \phi^{*} & 0 \end{array} \right)\!, \end{eqnarray*}$$

where

$$\begin{equation} \phi(\vec{k})=-t_{0}\sum_{\delta_{j}}\exp({\rm i}\vec{k}\cdot \vec{\delta_{j}}), \nonumber \end{equation} \tag{ A.1 }$$

which is given explicitly by equation (5). Note that the relation of this φ to the frequently defined [2, 28] quantity $f(\vec{k})$ is the following: $\vert \phi(\vec{k})\vert^2=3+f(\vec{k})$. The Hamiltonian (A.1) is written in the basis of the atomic orbital eigenfunctions:

$$\begin{eqnarray} \fl \psi_{\vec k}= \left( \begin{array}{@{}cccccc@{}} \psi_{A_1,k} &\psi_{B_1,k} & \psi_{A_2,k} & \psi_{B_2,k} & \psi_{A_3,k} & \psi_{B_3,k}\\ \end{array} \right)^{T}. \end{eqnarray} \tag{ A.2 }$$

The ABA trilayer is mirror symmetric with respect to the middle layer, therefore one can transform the Hamiltonian into a tridiagonal form using a unitary transformation [7, 8]

$$\begin{eqnarray} \hat{U}= \left( \begin{array}{@{}cccccc@{}} \frac{1}{\sqrt 2} & 0 & 0 & 0 & -\frac{1}{\sqrt 2} & 0 \\ 0 & \frac{1}{\sqrt 2} & 0 & 0 & 0 & -\frac{1}{\sqrt 2} \\ 0 & 0 &1 & 0 & 0 & 0 \\ 0& 0& 0 & 1 & 0 & 0 \\ \frac{1}{\sqrt 2} & 0 & 0 & 0 & \frac{1}{\sqrt 2} & 0 \\ 0 & \frac{1}{\sqrt 2} & 0 & 0 & 0 & \frac{1}{\sqrt 2} \end{array} \right)\!. \end{eqnarray} \tag{ A.3 }$$

In the new basis the Hamiltonian reads

$$\begin{eqnarray} \fl \hspace{-12pt}\hat{H}^\prime=\hat{U}\hat{H}\hat{U}^{-1}= \left( \begin{array}{@{}cccccc@{}} 0 & \phi& 0 & 0 & 0 & 0 \\ \phi^{*}& 0 & 0 & 0 & 0 & 0 \\ 0 & 0 &0 & \phi & 0 & \sqrt 2t_1 \\ 0& 0& \phi^{*}& 0 & 0 & 0 \\ 0& 0& 0& 0& 0 & \phi \\ 0& 0& \sqrt 2t_1& 0 & \phi^{*} & 0 \end{array} \right)\!. \end{eqnarray} \tag{ A.4 }$$

This matrix is composed of two blocks containing only intralayer hoppings and an effective interlayer one, ${\sqrt 2 t_1}$. It can easily be diagonalised yielding equations (2)–(4).

Here we derive the expressions for the characteristic energies presented in table 1, beyond the parabolic (linear) approximation for the band spectra. First we consider the case where μ = 0.2t₁, i.e. lies below the bottom of the band $E_{+3}(\vec{k})$. Among the first three possible transitions from the $E_{-1}(\vec{k})$ to the upper bands, the −1 → 3 transition threshold remains the same (see equation (13)) and similar to the case of μ = 0. Transition from −1 → 1 become possible if ℏω > 2μ, then for μ = 0.2t₁ the allowed energies are ℏω > 0.4t₁.

The domain of allowed energies for the the transition −1 → 2 can be obtained as $\hbar\omega>E_{+2}-E_{-1}=\sqrt {\vert \phi\vert^2}+\sqrt {t_{1}^2+\vert \phi\vert^2-t_{1}\sqrt {t_{1}^2+2\vert \phi\vert^2}}$. Since $E_{+2}=\sqrt {\vert \phi\vert^2}=\mu$, we obtain

$$\begin{eqnarray} &&\hbar\omega>\mu+\sqrt {t_{1}^2+\mu^2-t_1\sqrt {t_{1}^2+2\mu^2}}=0.22t_{1}, \end{eqnarray} \tag{ B.1 }$$

the number corresponds to μ = 0.2t₁.

For the second group of transitions from the −2 band, the −2 → 3 threshold remains as given by equation (15). The next one, −2 → 2, can occur for energies higher than 2μ, namely ℏω > 0.4t₁. Transition −2 → 1 is allowed if $\hbar\omega>E_{+1}-E_{-2}=\sqrt {t_{1}^2+\vert \phi\vert^2-t_{1}\sqrt {t_{1}^2+2\vert \phi\vert^2}}+\sqrt {\vert \phi\vert^2}$. Since $\mu=E_{+1}=\sqrt {t_{1}^2+\vert \phi\vert^2-t_{1}\sqrt {t_{1}^2+2\vert \phi\vert^2}}$, by solving this equation for $\vert \phi\vert^2=\mu^2+\sqrt 2t_1\mu$ we obtain the relation:

$$\begin{eqnarray} &&\hbar\omega>\mu+\sqrt {\mu^2+\sqrt 2t_{1}\mu}=0.77t_1. \end{eqnarray} \tag{ B.2 }$$

Finally, for the last group of three possible transitions from E₋₃ to the three upper bands, for the −3 → 3 one equation (17)] holds beyond the parabolic approximation, −3 → 1 is allowed if $\hbar\omega >\sqrt 2t_1+2{\mu}$ (then ℏω > 1.8t₁), and the −3 → 2 transition is allowed if

$$\begin{eqnarray} &&\hbar\omega>\mu+\sqrt {t_1^2+\mu^2+t_1\sqrt {t_1^2+2\mu^2}}=1.64t_1. \end{eqnarray} \tag{ B.3 }$$

The non-zero chemical potential gives rise to new transitions that are not possible in the intrinsic material. The 1 → 3 transition is allowed for the frequencies in the vicinity of $\omega\equiv \sqrt 2 t_{1}/\hbar$. The 1 → 2 one occurs in the domain

$$\begin{eqnarray} &&\mu-\sqrt {t_{1}^2+\mu^2-t_1\sqrt {t_{1}^2+2\mu^2}}<\hbar\omega<\sqrt {\mu^2+\sqrt {2}\mu t_1}-\mu, \nonumber\\ \end{eqnarray} \tag{ B.4 }$$

For μ = 0.2t₁ this corresponds to 0.17t₁ < ℏω < 0.37t₁.

The situation is similar for the 2 → 3 transition because, in addition to a threshold energy, there is also a cut-off:

$$\begin{eqnarray} &&\sqrt {\mu^2+t_1^2+t_1\sqrt {t_1^2+2\mu^2}}-\mu<\hbar\omega<\sqrt 2 t_1\!, \end{eqnarray} \tag{ B.5 }$$

where the expression in the left hand side is equal to 1.24t₁ for μ = 0.2t₁.

The case of large chemical potential, (μ = 1.6t₁), above the bottom of the band E₊₃ corresponds to the situation where the third band is partially filled. For ℏω > 2μ, the −1 → 1 transition becomes possible. The frequency domain for transition −1 → 2 is similar to the previous case (equation (B.1)) but for μ = 1.6t₁ we have ℏω > 2.64t₁.

Transition from −1 to 3 is allowed for $\hbar\omega>E_{+3}-E_{-1}=\sqrt {t_{1}^2+\vert \phi\vert^2+t_{1}\sqrt {t_{1}^2+2\vert \phi\vert^2}}+\sqrt {t_{1}^2+\vert \phi\vert^2-t_{1}\sqrt {t_{1}^2+2\vert \phi\vert^2}}$. From the dispersion relation μ = E₊₃ we obtain $\vert \phi\vert^2=\mu^2-\sqrt 2t_1\mu$, and $\hbar\omega>2\mu-\sqrt 2t_{1}$. Notice that this equation is exactly the same as equation (13), obtained within the parabolic approximation.

Among the next three possible transitions, the energy domain for −2 → 1 is given by equation (B.2). For μ = 1.6t₁ we have ℏω > 3.8t₁. The −2 → 2 transition is possible when ℏω > 2μ. For −2 → 3, from the dispersion relation μ = E₊₃ we have $\vert \phi\vert^2=\mu^2-\sqrt 2t_1\mu$. Hence, $\hbar\omega>E_{+3}-E_{-2}=\mu+\sqrt {\mu^2-\sqrt {2}t_1\mu}=2.15t_1$.

The transitions −3 → 1 and −3 → 3 are allowed if $\hbar\omega>2\mu+\sqrt 2t_1=4.8{t_{1}}$ (similar to previous case) and ℏω > 2μ, respectively. Transition from −3 to 2 is described by equation (B.3), then ℏω > 4.05t₁.

Finally, 1 → 2 is described by equation (B.4) (then 0.56t₁ ⩽ ℏω ⩽ 0.6t₁); 1 → 3 becomes possible for $\hbar\omega\equiv \sqrt 2 t_1$ and 2 → 3 for the finite frequency domain defined by $\sqrt {t_{1}^2+\mu^2+t_1\sqrt {t_{1}^2+2\mu^2}}-\mu<\hbar\omega<\mu-\sqrt {\mu^2-\sqrt 2t_{1}{\mu}}$. For μ = 1.6t₁ it means 0.86t₁ < ℏω < 1.06t₁.