Stacking-dependent optical absorption in multilayer graphene

We study the infrared optical absorption properties of ABA- and ABC-stacked graphene multilayers using the effective mass approximation. We calculate the optical absorption spectrum at various carrier densities, and find the characteristic features that identify the stacking types and the number of layers. We fully include the band parameters and discuss detailed features such as trigonal warping and electron–hole asymmetry.

In this work, we study the dynamical conductivity of AB A 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 AB A and ABC multilayers (from three-to five-layer) in sections 4 and 5, respectively.

Monolayer graphene
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][68][69][70][71][72][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 where p ± = p x ± i p y , p = −ih∇ and v ≈ 1 × 10 6 m s −1 is the band velocity [20,74]. The velocity v is expressed as v = ( √ 3/2)aγ 0 /h, where a = 0.246 nm is the lattice constant of the honeycomb lattice and γ 0 = 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 with electron momentum p = ( p x , p y ) as quantum numbers and p = p 2 x + p 2 y . The carrier density at a given Fermi energy ε F is given by 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 where S is the area of the system, v x = ∂H/∂ 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 constanth/τ ≡ 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] For graphene, σ (ω) is written at zero temperature as [59] σ (ω) = g v g s 16 where the first term in brackets represents the Drude (intra-band) part and the second the interband part. In the limit of → 0, its real part goes to Re σ (ω) = g v g s 16 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 whenhω > 2|ε F |, at the universal value, The transmission of light incident perpendicular to the graphene sheet is then showing that the absorption is given by πα ≈ 0.023 independent of frequency or wavelength, with the fine structure constant α ≡ e 2 /(ch) ≈ 1/137. This small absorption was experimentally observed [47][48][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 γ 1 ≈ 0.4 eV, and the unit of charge density to where γ 1 is the interlayer coupling in bulk graphite to be explained in the next section. The constant energy broadening is set to = 0.01γ 1 throughout the paper. We see the Drude peak at ω = 0 and the step increase of inter-band transition athω = 2|ε F |, as expected from equation (7). 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
Bilayer graphene [20][21][22][23][24] is a pair of graphene layers arranged in AB (Bernal) stacking as shown in figure 1(b). A unit cell includes A 1 and B 1 atoms on layer 1 and A 2 and B 2 on layer 2, and the layers are stacked with the interlayer spacing of d ≈ 0.334 nm, such that pairs of sites B 1 and A 2 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][26][27][28][29][30][31][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][12][13][14][15][16]. The parameter γ 0 describes nearest-neighbor coupling between A i and B i within each layer, as already introduced for monolayer graphene. γ 1 describes nearest-layer coupling between sites (B 1 and A 2 ) that lie directly above or below each other, γ 3 describes nearest-layer coupling between sites A 1 and B 2 and γ 4 is another nearest-layer coupling between A 1 and A 2 , and between B 1 and B 2 . The parameters γ 2 and γ 5 couple vertically located atoms on the next-nearest neighboring layers, and are thus absent in bilayer graphene. In AB A multilayer graphene illustrated in figure 1(c), γ 5 connects a pair of atoms such as B 1 and B 3 belonging to the vertical column connected by γ 1 , whereas γ 2 is for those which are not involved in the column, such as A 1 and A 3 . We define as the energy difference between the sites which are in the γ 1 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 = − γ 2 + γ 5 . For typical values of bulk graphite, we quote [16] γ 1 = 0.39 eV, γ 3 = 0.315 eV, γ 4 = 0.044 eV, γ 2 = −0.02 eV, γ 5 = 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 1 , |B 1 , |A 2 , |B 2 ) is given by [25] where v and p ± are the same as for monolayer graphene, and v 3 and v 4 are defined as v 3 = ( √ 3/2)(aγ 3 /h) and v 4 = ( √ 3/2)(aγ 4 /h). 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 gateinduced 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 3 , v 4 and are switched off in equation (11), the eigenenergies become [25,26] 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 ±γ 1 . In the following, we use the notation µ = 2, 1 instead of +, − and specify the four energy bands as (s, µ) = (±, 1), (±, 2). In ε γ 1 , the subbands (±, 1) touching at zero energy are approximately expressed as a quadratic form In this model (equation (12)), the inter-band part of the dynamical conductivity is explicitly written as [38,62] Re σ x x (ω) = g v g s 16 where we assumed that |ε F | < γ 1 . 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 γ 1 . When the parameters v 3 , v 4 and are included, the low-energy band equation (13) is modified to [44]   where tan ϕ = p y / p x and ξ = ± is the valley index (K and K ). The parameter v 3 produces trigonal warping where the energy band around each valley is stretched in 120 • symmetry [25]. The tiny parameter v 4 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 /γ 1 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 σ x x (ω) calculated for the bilayer graphene of equation (11) with the band parameters fully included. The curve for the nondoped case n s = 0 has essentially no prominent structure except for step-like increases near hω = γ 1 , 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 nearhω = γ 1 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 γ 3 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 athω ∼ 2ε F such as monolayer graphene's, which corresponds to the first term in equation (14).

AB A-stacked multilayer graphene
The calculation of dynamical conductivity can be extended to graphene multilayers with more than two layers in a straightforward manner. Here we consider an AB A-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 AB A stacking, the sites B 1 , A 2 , B 3 , A 4 , . . . are arranged along vertical columns normal to the layer plane, while the rest of the sites A 1 , B 2 , A 3 , B 4 , . . . are above or below the center of hexagons in the neighboring layers. If the bases are sorted as |A 1 , |B 1 ; |A 2 , |B 2 ; . . . ; |A N , |B N , the Hamiltonian around the K point becomes [26,28,29,31,32,34,38]   with 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 where X is A or B, and and 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. where and U(λ m ), appearing in diagonal blocks (m = m ), is equivalent to the Hamiltonian of bilayer graphene with nearest-layer coupling parameters multiplied by λ m . W with m = m adds on-site asymmetric potential to this effective bilayer system U(λ m ). For a pair of low-energy bands near zero energy, it causes an overall energy shift of (α + β)γ 2 /2 and an energy gap of ∼|(α − β)γ 2 | 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) 0 , |φ (B,odd) 0 } survive in equation (20). The matrix for the m = 0 block for the two remaining bases is 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, γ 2 and γ 5 . 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 λ 2 = √ 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 fromhω ∼ γ 1 to λ 2 γ 1 . 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]. 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 (λ 3 = (1+ √ 5)/2, indicated as 'B') and another bilayer-like band with a narrow splitting (λ 1 = (−1 + √ 5)/2, 'b'). As a result, the spectrum is composed of a broader peak athω ∼ λ 3 γ 1 , and a narrower peak at λ 1 γ 1 . The same is true for five-layer graphene which contains two bilayer blocks λ 4 = √ 3 ('B') and λ 2 = 1 ('b') and one monolayer block ('M').

ABC-stacked multilayer graphene
The electronic structure of the ABC-stacked multilayer is quite different from AB A 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 1 , |B 1 ; |A 2 , |B 2 ; . . . ; |A N , |B N , the Hamiltonian around the K point becomes [19,26,40,46] 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 γ 1 . 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 γ 1 [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 ε = ±γ 1 , and two-fold degenerate states localized on A 1 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 ε = ±γ 1 , whereas zero-energy states at A 1 and B N are coupled by N th-order perturbation in p, leading to weak dispersion (vp) N /γ N −1 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 ±γ 1 . 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. γ 3 and γ 2 cause trigonal warping of the band dispersion [44] and γ 5 and cause a shift in the dimer-state bands at ±γ 1 toward the upper energy.
The absorption spectra of ABC multilayers shown in the right panels of figures 7-9 are obviously different from their AB A 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 thanhω = γ 1 . 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 γ 1 . 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 ±γ 1 . We also see that the transitions involving (±, 3) become gradually observable, as the band bottom of (±, 3) becomes flatter with increasing N .

Conclusion
We calculated the optical absorption spectra of graphene multilayers with AB A 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 AB A 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 <hω < 2γ 1 .
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 ±γ 1 , and the main peak energies are always around γ 1 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 AB A and ABC multilayers is left for future works.