Stacking-dependent optical absorption in multilayer graphene

In multilayer graphene systems, the weak interlayer coupling allows a variety of stacking structures, giving different electronic properties depending on the atomic configuration [1–10]. The most stable stacking arrangement between adjacent layers is known to be AB-stacking, where A sublattice on one honeycomb layer comes right below (or above) B sublattice on the other as illustrated in figure 1(b). In a graphene stack of more than three layers, there are two possible forms called ABA (hexagonal or Bernal) [11–16] and ABC (rhombohedral) stacking [17–19] as shown in figures 1(c) and (d), respectively. The ABA phase is the most stable and common, while some portion in natural graphite takes the ABC phase partly. The electronic band structure of multilayer graphene sensitively depends on the stacking arrangement. In AB-stacked bilayer graphene [20–24], the low-energy band dispersion becomes quadratic unlike the linear Dirac cone in a single layer [25–32]. In ABA-stacked multilayer, the spectrum consists of quadratic bands analogous to bilayer graphene and a linear band-like monolayer [26, 28, 29, 31–39], whereas the ABC multilayer has a completely different spectrum with a pair of flat low-energy bands at zero energy [26, 33, 35, 36, 40–46].

Zoom In Zoom Out Reset image size

To investigate the electronic properties of layered graphene systems, optical measurement techniques are widely adopted. Raman spectroscopy is commonly used to determine the number of layers as well as the stacking structure [1, 4]. The angle-resolved photoemission spectroscopy is used to observe the electronic dispersion directly [21, 23]. Infrared spectroscopy is also widely applied to probe the low-energy band structure [47–58]. Theoretically, the optical absorption is related to the dynamical conductivity, and has been extensively studied for graphene and its multilayers [37, 38, 59–66]. The dependence of the infrared optical absorption on the stacking sequence has been observed in recent experiments [58] and also theoretically studied for graphene multilayers at neutral charge [66].

In this work, we study the dynamical conductivity of ABA and ABC graphene multilayers with various charge densities using the effective mass model. For each system, we calculate the spectral evolution as a function of charge density, and find several characteristic features that specify the stacking structure and the number of layers. We fully include possible band parameters and discuss the effect of additional features such as trigonal warping and electron–hole asymmetry. We study the graphene monolayer and AB-stacked bilayer in sections 2 and 3, respectively, and then consider the ABA and ABC multilayers (from three- to five-layer) in sections 4 and 5, respectively.

Graphene is composed of a honeycomb network of carbon atoms as shown in figure 1(a), where a unit cell contains a pair of sublattices denoted by A and B. Electronic states in the vicinity of the K and K' points in the Brillouin zone are well described by the effective mass approximation [12, 13, 67–73]. Let |A〉 and |B〉 be the Bloch functions at the K point, corresponding to the A and B sublattices, respectively. In a basis (|A〉,|B〉), the Hamiltonian around the K point becomes

$$\begin{equation} {\mathcal H}^K = \left(\begin{array}{@{}cc@{}} 0 & v p_- \\ v p_+ & 0 \end{array}\right), \end{equation} \tag{ 1 }$$

where p_± = p_x ± ip_y, ${\bf p} = -\mathrm{i}\hbar \boldsymbol{\nabla}$ and v ≈ 1 × 10⁶ m s⁻¹ is the band velocity [20, 74]. The velocity v is expressed as $v = (\sqrt {3}/2)a\gamma_0/\hbar$, where a = 0.246 nm is the lattice constant of the honeycomb lattice and γ₀ = 3.16 eV is a band parameter that describes the electron hopping energy between the nearest A and B sites. The Hamiltonian at the K' point is obtained by exchanging p_± in equation (1). The energy band is given by

$$\begin{eqnarray} &&\varepsilon_{\mathrm{s}}(p) = s v p\quad (s=\pm1) \end{eqnarray} \tag{ 2 }$$

with electron momentum p = (p_x,p_y) as quantum numbers and $p=\sqrt {p_x^2+p_y^2}$. The carrier density at a given Fermi energy ε_F is given by

$$\begin{equation} n_{\mathrm{s}}(\varepsilon_{\mathrm{F}}) = \frac{g_{\mathrm{v}} g_{\mathrm{s}}}{4\pi\hbar^2 v^2}\, \varepsilon_{\mathrm{F}}^2\,\, \mathrm{sgn}(\varepsilon_{\mathrm{F}}), \end{equation} \tag{ 3 }$$

where g_s = 2 and g_v = 2 represent the degrees of freedom associated with spin and valley, respectively, and n_s > 0 and n_s < 0 represent electron concentration and hole concentration, respectively.

The optical absorption is related to the dynamical conductivity of the system. Within the linear response, the dynamical conductivity is generally given by

$$\begin{eqnarray} &&\sigma(\omega) = \frac{e^2\hbar}{\mathrm{i}S} \sum_{\alpha,\beta} \frac{f(\varepsilon_\alpha)-f(\varepsilon_\beta)}{\varepsilon_\alpha- \varepsilon_\beta} \frac{| \langle\alpha | v_x | \beta\rangle |^2} {\varepsilon_\alpha- \varepsilon_\beta+\hbar\omega+\mathrm{i}\delta}, \end{eqnarray} \tag{ 4 }$$

where S is the area of the system, $v_x = \partial {\cal H}/\partial p_x$ is the velocity operator, δ is the positive infinitesimal, f(ε) is the Fermi distribution function and |α〉 and ε_α describe the eigenstate and the eigenenergy of the system. In the following, we include the disorder effect by replacing δ with the phenomenological constant ℏ/τ ≡ 2Γ throughout the paper. The optical absorption intensity is related to the real part of σ(ω). The transmission of the light incident perpendicular to the two-dimensional system is given by [75]

$$\begin{equation} T = \Big| 1+ \frac{2\pi}{c}\sigma(\omega) \Big|^{-2} \approx 1- \frac{4\pi}{c}\mathrm{Re}\,\sigma(\omega). \end{equation} \tag{ 5 }$$

For graphene, σ(ω) is written at zero temperature as [59]

$$\begin{eqnarray} &&\sigma(\omega) = \frac{g_{\mathrm{v}}g_{\mathrm{s}}}{16}\frac{e^2}{\hbar} \left[\frac{4}{\pi}\frac{\mathrm{i}\varepsilon_{\mathrm{F}}}{\hbar\omega+2{\rm i}\Gamma}+\left(1+\frac{\mathrm{i}}{{\pi}}\ln\frac{\hbar\omega-2\varepsilon_{\mathrm{F}}+2\mathrm{i}\Gamma}{\hbar\omega+2\varepsilon_{\mathrm{F}}+2\mathrm{i}\Gamma}\right)\right], \end{eqnarray} \tag{ 6 }$$

where the first term in brackets represents the Drude (intra-band) part and the second the inter-band part. In the limit of Γ → 0, its real part goes to

$$\begin{equation} \mathrm{Re\,} \sigma(\omega) = \frac{g_{\mathrm{v}}g_{\mathrm{s}}}{16}\frac{e^2}{\hbar}\left[4\varepsilon_{\mathrm{F}} \delta(\hbar\omega)+\theta(\hbar\omega-2|\varepsilon_{\mathrm{F}}|)\right], \end{equation} \tag{ 7 }$$

where the Drude part becomes a delta function centered at ω = 0, and the inter-band part becomes a step function which is non-zero only when ℏω > 2|ε_F|, at the universal value,

$$\begin{equation} \sigma_{\mathrm{mono}} = \frac{g_{\mathrm{v}}g_{\mathrm{s}}}{16}\frac{e^2}{\hbar}. \end{equation} \tag{ 8 }$$

The transmission of light incident perpendicular to the graphene sheet is then

$$\begin{equation} T \approx 1- \frac{4\pi}{c}\sigma_{\mathrm{mono}} = 1 - \frac{g_{\mathrm{v}}g_{\mathrm{s}}}{4}\frac{\pi e^2}{c\hbar} \end{equation} \tag{ 9 }$$

showing that the absorption is given by πα ≈ 0.023 independent of frequency or wavelength, with the fine structure constant α ≡ e²/(cℏ) ≈ 1/137. This small absorption was experimentally observed [47–49], and is often used to identify the number of layers of graphene flakes created by mechanical exfoliation.

The real part of σ(ω) for several different charge densities n_s is plotted in figure 2(b). The Hamiltonian of monolayer graphene is essentially scaleless, but to compare the multilayer's result in later sections, we set the unit of energy to γ₁ ≈ 0.4 eV, and the unit of charge density to

$$\begin{equation} n_0 = \frac{g_{\mathrm{v}}g_{\mathrm{s}}}{2\pi}\left(\frac{\gamma_1}{\hbar v}\right)^2, \end{equation}\noindent \tag{ 10 }$$

where γ₁ is the interlayer coupling in bulk graphite to be explained in the next section. The constant energy broadening is set to Γ = 0.01γ₁ throughout the paper. We see the Drude peak at ω = 0 and the step increase of inter-band transition at ℏω = 2|ε_F|, as expected from equation (7).

Zoom In Zoom Out Reset image size

In a real system, the energy broadening Γ is not actually a constant but generally depends on the Fermi energy ε_F, and this becomes crucial when considering the static conductivity σ(ω = 0) at the Dirac point ε_F = 0. In the Boltzmann approximation, for example, [70] Γ is shown to be proportional to ε_F, so that the Drude part (the first term in equation (6)) at ω = 0 approaches a non-zero value in the limit of ε_F → 0, while it just vanishes at non-zero ω. Therefore, σ(ω;ε_F → 0) as a function of ω has a singularity at ω = 0 [59]. This anomaly remains even in more realistic calculations properly including the level broadening effects such as self-consistent Born approximation [70].

Bilayer graphene [20–24] is a pair of graphene layers arranged in AB (Bernal) stacking as shown in figure 1(b). A unit cell includes A₁ and B₁ atoms on layer 1 and A₂ and B₂ on layer 2, and the layers are stacked with the interlayer spacing of d ≈ 0.334 nm, such that pairs of sites B₁ and A₂ lie directly above or below each other. There the interlayer coupling drastically changes the linear band structure of the monolayer, leaving a pair of quadratic energy bands touching at zero energy [25–32].

The effective mass Hamiltonian for bilayer graphene can be built using the Slonczewski–Weiss–McClure parameterization of tight-binding couplings of bulk graphite [11–16]. The parameter γ₀ describes nearest-neighbor coupling between A_i and B_i within each layer, as already introduced for monolayer graphene. γ₁ describes nearest-layer coupling between sites (B₁ and A₂) that lie directly above or below each other, γ₃ describes nearest-layer coupling between sites A₁ and B₂ and γ₄ is another nearest-layer coupling between A₁ and A₂, and between B₁ and B₂. The parameters γ₂ and γ₅ couple vertically located atoms on the next-nearest neighboring layers, and are thus absent in bilayer graphene. In ABA multilayer graphene illustrated in figure 1(c), γ₅ connects a pair of atoms such as B₁ and B₃ belonging to the vertical column connected by γ₁, whereas γ₂ is for those which are not involved in the column, such as A₁ and A₃. We define Δ' as the energy difference between the sites which are in the γ₁ column and the sites which are not. In the usual parameterization of the Slonczewski–Weiss–McClure model, we have Δ instead of Δ', and the two are related by Δ' = Δ − γ₂ + γ₅. For typical values of bulk graphite, we quote [16] γ₁ = 0.39 eV, γ₃ = 0.315 eV, γ₄ = 0.044 eV, γ₂ = −0.02 eV, γ₅ = 0.04 eV and Δ' = 0.047 eV.

Using the above parameters, the effective Hamiltonian of bilayer graphene at the K point for the basis (|A₁〉,|B₁〉, |A₂〉,|B₂〉) is given by [25]

$$\begin{eqnarray} {\cal H} = \left(\!\begin{array}{@{}cccc@{}} 0 & v p_- & -v_4 p_- & v_3 p_+ \\ v p_+ & \Delta' & \gamma_1 & -v_4 p_- \\ -v_4 p_+ & \gamma_1 & \Delta' & v p_- \\ v_3 p_- & -v_4 p_+ & v p_+ & 0 \end{array}\!\right), \end{eqnarray} \noindent \tag{ 11 }$$

where v and p_± are the same as for monolayer graphene, and v₃ and v₄ are defined as $v_3 = (\sqrt {3}/2)(a \gamma _3/\hbar )$ and $v_4 = (\sqrt {3}/2)(a \gamma _4/\hbar )$. The effective Hamiltonian for K' can be obtained again by exchanging p₊ and p₋, giving the equivalent spectrum.

In the experimental situations, the top and the bottom layers may have different electrostatic potentials. This is particularly important when the charge carriers are doped by the gate electrode, which induces a perpendicular electric field at the same time. The interlayer potential difference generally modifies the band structure of the multilayer graphene system [25, 26, 28–30, 46, 76, 77] and affects the optical absorption spectrum [56, 78]. The effect of the gate-induced band distortion on the optical spectrum was investigated for the bilayer graphene [78]. Here, for simplicity, we neglect the interlayer potential difference even in the doped case, while such a situation can be actually realized in an experimental setup using both the top and bottom gate electrodes [79].

When the parameters v₃, v₄ and Δ' are switched off in equation (11), the eigenenergies become [25, 26]

$$\begin{equation} \varepsilon_{\mathrm{s},\mu}({\bf p}) = s\left[\frac{\mu}{2}\gamma_1+ \sqrt{\frac{1}{4} \gamma_1^2 + v^2p^2} \right]. \end{equation} \tag{ 12 }$$

The branch μ = − gives a pair of conduction (s = +) and valence (s = −) bands touching at zero energy. The other branch μ = + is another pair repelled away by ±γ₁. In the following, we use the notation μ = 2,1 instead of +,− and specify the four energy bands as (s,μ) = (±,1),(±,2). In ε ≪ γ₁, the subbands (±,1) touching at zero energy are approximately expressed as a quadratic form

$$\begin{eqnarray} &&\epsilon_{\pm,1}({\bf p}) \approx \pm \frac{v^2p^2}{\gamma_1}. \end{eqnarray} \tag{ 13 }$$

In this model (equation (12)), the inter-band part of the dynamical conductivity is explicitly written as [38, 62]

$$\begin{eqnarray} &&\fl \mathrm{Re} \, \sigma_{xx} (\omega) = \frac{g_{\mathrm{v}}g_{\mathrm{s}}}{16}\frac{e^2}{\hbar} \left\{ \frac{\hbar\omega+2\gamma_1}{\hbar\omega+\gamma_1} \theta(\hbar\omega - 2|\varepsilon_{\mathrm{F}}|) + \left(\frac{\gamma_1}{\hbar\omega}\right)^2 \left[ \theta(\hbar\omega -\gamma_1) + \theta(\hbar\omega -\gamma_1- 2|\varepsilon_{\mathrm{F}}|) \right] \right.\nonumber\\ &&\left. + \frac{\hbar\omega-2\gamma_1}{\hbar\omega-\gamma_1} \theta(\hbar\omega -2\gamma_1) + \gamma_1 \log\left[\frac{2|\varepsilon_{\mathrm{F}}|+\gamma_1}{\gamma_1}\right] \delta(\hbar\omega-\gamma_1) \vphantom{\frac{\hbar\omega+2\gamma_1}{\hbar\omega+\gamma_1} \theta(\hbar\omega - 2|\varepsilon_{\mathrm{F}}|) + \left(\frac{\gamma_1}{\hbar\omega}\right)^2 \left[ \theta(\hbar\omega -\gamma_1) + \theta(\hbar\omega -\gamma_1- 2|\varepsilon_{\mathrm{F}}|) \right]}\right\}, \end{eqnarray} \noindent \tag{ 14 }$$

where we assumed that |ε_F| < γ₁. The first three terms in brackets are the absorptions from the valence band to the conduction band, where the first term is from the band (−,1) to (+,1), the second from (∓,1) to (±,2) and the third from (−,2) to (+,2). The fourth term is present only in electron- or hole-doped cases, and is caused by a transition between the two conduction bands (+,1) and (+,2) or between the valence bands (−,1) and (−,2). It is expressed by a delta function in ω because in the present approximation (equation (12)), the bands (±,1) and (±,2) have identical dispersion except for an energy shift by γ₁.

When the parameters v₃, v₄ and Δ' are included, the low-energy band equation (13) is modified to [44]

$$\begin{eqnarray} &&\epsilon_{\pm,1}({\bf p}) \approx \pm \sqrt{\frac{v^4p^4}{\gamma_1^2} - 2 \xi \frac{v_3v^2p^3}{\gamma_1}\cos 3\varphi + v_3^2 p^2 } + \frac{2 v v_4 p^2}{\gamma_1}, \end{eqnarray} \tag{ 15 }$$

where tan φ = p_y/p_x and ξ = ± is the valley index (K and K'). The parameter v₃ produces trigonal warping where the energy band around each valley is stretched in 120° symmetry [25]. The tiny parameter v₄ gives the common quadratic term to both the conduction and valence bands, and thus breaks the electron–hole asymmetry. The parameter Δ' appears in the second order of Δ'/γ₁ in equation (15) and is neglected there, while it gives a positive energy shift Δ' to both (±,2). As a result, the energy distance between (+,1) and (+,2) bands becomes larger than between (−,1) and (−,2).

Figures 3(a) and (b) are the band structure and Re σ_xx(ω) calculated for the bilayer graphene of equation (11) with the band parameters fully included. The curve for the non-doped case n_s = 0 has essentially no prominent structure except for step-like increases near ℏω = γ₁, corresponding to transitions from (−,1) to (+,2) and from (−,2) to (+,1). The energy difference between the two transitions is due to the electron–hole asymmetry introduced by Δ'. When ε_F is shifted from the charge neutral point, a strong peak appears near ℏω = γ₁ given by the transition from (±,1) to (±,2) with ± for electron and hole doping regions, respectively. The peak corresponds to the delta function in equation (14), but it is now broadened because the trigonal warping caused by γ₃ distorts the lower energy bands (±,1) more seriously than (±,2), so that the transition energy depends on k. We also observe a step-like feature at ℏω ∼ 2ε_F such as monolayer graphene's, which corresponds to the first term in equation (14).

Zoom In Zoom Out Reset image size

The calculation of dynamical conductivity can be extended to graphene multilayers with more than two layers in a straightforward manner. Here we consider an ABA-stacked graphene multilayer composed of N layers illustrated in figure 1(c), which contains A_j and B_j atoms on the layer j = 1,...,N. In ABA stacking, the sites B₁,A₂,B₃,A₄,... are arranged along vertical columns normal to the layer plane, while the rest of the sites A₁,B₂,A₃,B₄,... are above or below the center of hexagons in the neighboring layers. If the bases are sorted as |A₁〉,|B₁〉; |A₂〉,|B₂〉; ...; |A_N〉,|B_N〉, the Hamiltonian around the K point becomes [26, 28, 29, 31, 32, 34, 38]

$$\begin{eqnarray} {\cal H} = \left(\begin{array}{@{}ccccccc@{}} H_1 & V & W & & \\ V^{\dagger} & H_2 & V^{\dagger} & W' & \\ W & V & H_3 & V & W & \\ & W' & V^{\dagger} & H_4 & V^{\dagger} & W' & \\ & & \ddots & \ddots & \ddots & \ddots & \ddots \end{array}\right) \end{eqnarray} \tag{ 16 }$$

with

$$\begin{eqnarray} H_j = \left\{ \begin{array}{@{}l@{}} \left(\begin{array}{@{}cc@{}} 0 & v p_- \\ v p_+ & \Delta' \end{array}\right) \quad (j=\mathrm{odd}), \\ \left(\begin{array}{@{}cc@{}} \Delta' & v p_- \\ v p_+ & 0 \end{array}\right) \quad (j=\mathrm{even}), \end{array} \right. \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} \fl V = \left(\begin{array}{@{}cc@{}} -v_4 p_- & v_3 p_+ \\ \gamma_1 & -v_4 p_- \end{array}\right),\quad W = \left(\begin{array}{@{}cc@{}} \gamma_2/2 & 0 \\ 0 & \gamma_5/2 \end{array}\right),\quad W' = \left(\begin{array}{@{}cc@{}} \gamma_5/2 & 0 \\ 0 & \gamma_2/2 \end{array}\right) \end{eqnarray} \noindent \tag{ 18 }$$

with the band parameters already introduced in the previous section. The diagonal blocks H_j describe intralayer coupling in layer j, V nearest-neighbor interlayer coupling and W and W' next-nearest-neighbor interlayer coupling. The effective Hamiltonian for K' is obtained by exchanging p₊ and p₋.

The Hamiltonian can be approximately decomposed into subsystems equivalent to monolayer or bilayer graphene [34, 37, 38, 80]. The dynamical conductivity of the system can then be written as a summation over each sub-Hamiltonian, so that the absorption spectrum becomes essentially equivalent to the collection of the effective monolayer and bilayer's spectra [37, 38]. The decomposed subsystems are labeled by index m, which is defined as

$$\begin{eqnarray} &&m = \left\{ \begin{array}{@{}l@{}} 1,3,5, \ldots , N-1,\quad N = \mathrm{even} \, , \\ 0,2,4, \ldots , N-1,\quad N = \mathrm{odd} \, . \end{array} \right. \end{eqnarray} \tag{ 19 }$$

The subsystem m is composed of four bases u_m = {|ϕ_m^(A,odd)〉,|ϕ_m^(B,odd)〉,|ϕ_m^(A,even)〉,|ϕ_m^(B,even)〉} with

$$\begin{eqnarray} &&|\phi_m^{(X, \mathrm{odd})}\rangle = \sum_{j=1}^N f_m(j) |X_j\rangle, \quad |\phi_m^{(X, \mathrm{even})}\rangle = \sum_{j=1}^N g_m(j) |X_j\rangle \, , \end{eqnarray} \tag{ 20 }$$

where X is A or B, and

$$\begin{equation} f_m(j) = c_m \sqrt{\frac{2}{N+1}}[1 - (-1)^j] \sin \kappa_m j, \end{equation} \tag{ 21 }$$

$$\begin{equation} g_m(j) = c_m \sqrt{\frac{2}{N+1}}[1 + (-1)^j] \sin \kappa_m j \end{equation} \tag{ 22 }$$

and

$$\begin{equation} \kappa_m = \frac{\pi}{2}-\frac{m\pi}{2(N+1)}, \end{equation} \tag{ 23 }$$

$$\begin{equation} c_m = \left\{ \begin{array}{@{}ll@{}} 1/2 & (m=0) \\ 1/\sqrt{2} & (m\neq 0). \end{array} \right. \end{equation} \tag{ 24 }$$

Obviously, f_m(j) is zero on even j layers, while g_m(j) is zero on odd j layers. A superscript such as (A, odd) indicates that the wave function has a non-zero amplitude only on |A_j〉 sites with odd j.

The Hamiltonian matrix between u_m and u_m' becomes [38, 80]

$$\begin{eqnarray} &&{\cal H}_{m'm} \equiv {\bf u}^\dagger_{m'} {\cal H} {\bf u}_m = {\cal U}(\lambda_m)\delta_{m'm} + {\cal W}(\alpha_{m'm},\beta_{m'm}), \end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray} {\cal U}(\lambda) = \left(\begin{array}{@{}cccc@{}} 0 & vp_- & -\lambda v_4p_- & \lambda v_3p_+ \\ vp_+ & \Delta' & \lambda \gamma_1 & -\lambda v_4p_- \\ -\lambda v_4p_+ & \lambda \gamma_1 & \Delta' & vp_- \\ \lambda v_3p_- & -\lambda v_4p_+ & vp_+ & 0 \end{array}\!\right), \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} {\cal W}(\alpha,\beta) = \left(\begin{array}{@{}cccc@{}} \alpha \gamma_2 & 0 & 0 & 0 \\ 0 & \alpha \gamma_5 & 0 & 0 \\ 0 & 0 & \beta \gamma_5 & 0 \\ 0 & 0 & 0 & \beta \gamma_2 \end{array}\right) \end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray} &&\fl \lambda_m = 2 \cos \kappa_m,\\ &&\fl \alpha_{m'm} = 2c_mc_{m'}\Biggl\{ \delta_{mm'}(1+\delta_{m0}) \cos 2\kappa_m + \frac{\sin\kappa_m \sin\kappa_{m'}}{N+1} \left[ 2+(-1)^{\frac{m-m'}{2}}(1-(-1)^N) \right] \Biggr\},\nonumber \end{eqnarray} \noindent \tag{ 28 }$$

$$\begin{eqnarray} &&\fl \beta_{m'm} = 2c_mc_{m'}\Biggl\{ \delta_{mm'}(1-\delta_{m0}) \cos 2\kappa_m + \frac{\sin\kappa_m \sin\kappa_{m'}}{N+1} (-1)^{\frac{m-m'}{2}}(1+(-1)^N) \Biggr\}. \end{eqnarray} \noindent \tag{ 29 }$$

${\cal U}(\lambda _m)$, appearing in diagonal blocks (m = m'), is equivalent to the Hamiltonian of bilayer graphene with nearest-layer coupling parameters multiplied by λ_m. ${\cal W}$ with m = m' adds on-site asymmetric potential to this effective bilayer system ${\cal U}(\lambda _m)$. For a pair of low-energy bands near zero energy, it causes an overall energy shift of (α + β)γ₂/2 and an energy gap of ∼|(α − β)γ₂| at the band touching point, as in asymmetric bilayer graphene.

The case of m = 0 is special in that g_m(j) is identically zero, so that only two basis states {|ϕ^(A,odd)₀〉, |ϕ^(B,odd)₀〉} survive in equation (20). The matrix for the m = 0 block for the two remaining bases is

$$\begin{eqnarray} {\cal H}_{0} = \left(\begin{array}{@{}cc@{}} 0 & vp_- \\ vp_+ & \Delta' \end{array}\right) - \frac{N-1}{N+1} \left(\begin{array}{@{}cc@{}} \gamma_2 & 0 \\ 0 & \gamma_5 \end{array}\right), \end{eqnarray} \tag{ 30 }$$

which, barring the diagonal terms, is equivalent to the Hamiltonian of monolayer graphene. Therefore, odd-layered graphene with N = 2M + 1 is composed of one monolayer-type and M bilayer-type subsystems, while even-layered graphene with N = 2M is composed of M bilayers but no monolayer. The subsystems with different m's are not exactly independent but mixed by off-diagonal matrix elements W with m ≠ m', which contains only the next-nearest interlayer couplings, γ₂ and γ₅. This is responsible for small band repulsion at crossing points between different m's, while the overall band structure is well described by only retaining the diagonal blocks [38].

The trilayer graphene N = 3 has a monolayer block (m = 0) and a bilayer block (m = 2) with $\lambda _2 = \sqrt {2}$. The band structure is shown in figure 4(a), where the symbols M and B represent the monolayer and the bilayer block. The dynamical conductivity for several n_s's is plotted in figure 4(b). This actually looks like a summation of monolayer in figure 2 and bilayer in figure 3, while the prominent peak of the bilayer block shifts from ℏω ∼ γ₁ to λ₂γ₁. Also, the peak is broader than the real bilayer graphene in figure 3, because the trigonal warping effect is more pronounced for larger λ_m [34].

Zoom In Zoom Out Reset image size

Figures 5 and 6 show similar plots for N = 4 and 5, respectively. The band structure of four-layer graphene is composed of a bilayer-like band with a wide splitting ($\lambda _3=(1{+}\sqrt {5})/2$, indicated as 'B') and another bilayer-like band with a narrow splitting ($\lambda _1=(-1+\sqrt{5})/2$, 'b'). As a result, the spectrum is composed of a broader peak at ℏω ∼ λ₃γ₁, and a narrower peak at λ₁γ₁. The same is true for five-layer graphene which contains two bilayer blocks $\lambda _4=\sqrt {3}$ ('B') and λ₂ = 1 ('b') and one monolayer block ('M').

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The electronic structure of the ABC-stacked multilayer is quite different from ABA and so is the optical absorption spectrum. In the ABC lattice structure, B_i and A_i+1 are always vertically located as illustrated in figure 1(d). In the bases |A₁〉,|B₁〉; |A₂〉,|B₂〉; ... ; |A_N〉,|B_N〉, the Hamiltonian around the K point becomes [19, 26, 40, 46]

$$\begin{eqnarray} {\cal H} = \left(\begin{array}{ccccccc} H_1 & V & W & & \\ V^{\dagger} & H_2 & V& W& \\ W^\dagger & V^\dagger & H_3 & V & W & \\ & \ddots & \ddots & \ddots & \ddots & \ddots \end{array}\right) \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} H_j = \left\{ \begin{array}{@{}l@{}} \left(\begin{array}{@{}cc@{}} 0 & v p_- \\ v p_+ & \Delta' \end{array}\right) \quad (j=1),\\ \left(\begin{array}{@{}cc@{}} \Delta' & v p_- \\ v p_+ & \Delta' \end{array}\right) \quad (j=2,\ldots,N-1),\\ \left(\begin{array}{@{}cc@{}} \Delta' & v p_- \\ v p_+ & 0 \end{array}\right) \quad (j=N), \end{array} \right. \end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray} V = \left(\begin{array}{@{}cc@{}} -v_4 p_- & v_3 p_+ \\ \gamma_1 & -v_4 p_- \end{array}\right), \quad W = \left(\begin{array}{@{}cc@{}} 0 & \gamma_2/2 \\ 0 & 0 \end{array}\right). \end{eqnarray} \tag{ 33 }$$

The effective Hamiltonian for K' is again obtained by exchanging p₊ and p₋.

The overall band structure can be understood by considering the simplest model only including v and γ₁. The system is then regarded as a one-dimensional chain depicted in the side panel of figure 1(d), in which A_j and B_j are coupled by vp_± and B_j and A_j+1 by γ₁ [26]. At p = 0, the connections between A_j and B_j are lost so that the eigenstate set is composed of (N − 1)-fold degenerate dimer states between B_j and A_j+1 at each of ε = ± γ₁, and two-fold degenerate states localized on A₁ and B_N at zero energy. When infinitesimal p is switched on, the dimer states are coupled by a term linear to p giving a set of (N − 1) linear bands at each of ε = ± γ₁, whereas zero-energy states at A₁ and B_N are coupled by Nth-order perturbation in p, leading to weak dispersion (vp)^N/γ^N−1₁ [26, 36, 42, 44]. We number the mth conduction band and valence band from the center as (+,m) and (−,m) (m = 1,2,...,N), respectively. The bands with m = 1 correspond to the nearly zero-energy bands and m = 2,...,N to the dimer-state bands at ±γ₁. The band structures of ABC multilayers N = 3,4,5 with all the band parameters included are presented in the left panels of figures 7–9, respectively. The effects of the additional parameters are similar to those in bilayer graphene. γ₃ and γ₂ cause trigonal warping of the band dispersion [44] and γ₅ and Δ' cause a shift in the dimer-state bands at ±γ₁ toward the upper energy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The absorption spectra of ABC multilayers shown in the right panels of figures 7–9 are obviously different from their ABA counterparts. At N = 3, the spectrum for the doped system is characterized by a strong peak due to the transition between (±,1) and (±,2) similar to bilayer graphene, whereas the peak center is significantly smaller than ℏω = γ₁. This is because the band energy of (+,2) decreases when p shifts from zero, so that the distance between (+,1) and (+,2) tends to be smaller than the value at p = 0, which is about γ₁. We also note that the peak of the transition between (∓,1) and (±,2) is much more pronounced here than in bilayer, owing to the large density of states of the flat surface state bands near p = 0. The transitions related to the highest band (±,3) are almost invisible, because the band dispersion is conical and contributes to only a small density of states. In ABC multilayers with N = 4 and 5, the peak of the transition between (±,1) and (±,2) moves even lower ω as N increases, reflecting deeper band bottom of (±,2) with respect to ±γ₁. We also see that the transitions involving (±,3) become gradually observable, as the band bottom of (±,3) becomes flatter with increasing N.

We calculated the optical absorption spectra of graphene multilayers with ABA and ABC stackings with several different thicknesses. We traced the spectral evolution with varying charge density, and observed the characteristic features that may be useful in identifying the stacking sequences as well as the number of layers. We found that the absorption spectrum of ABA multilayer is characterized by strong peaks corresponding to the transitions within the decomposed bilayer subband, and they can appear in a wide frequency region, 0 < ℏω < 2γ₁. In ABC multilayers, on the other hand, the main peaks are given by the transition from the zero-energy surface state band to the excited bands near ±γ₁, and the main peak energies are always around γ₁ or lower. Here we neglect the interlayer potential asymmetry given by the gate electric field, while it is expected to cause broadening and shift of the absorption peaks via the band deformation, as observed in the bilayer graphene [56, 78]. Detailed investigation of those effects for ABA and ABC multilayers is left for future works.

This work was funded by the JST–EPSRC Japan–UK Cooperative Programme under grant no. EP/H025804/1 and by Grants-in-Aid for Scientific Research No. 24740193 from JSPS.