An analytical approach for the energy spectrum and optical properties of gated bilayer graphene

Abstract

First, we present an analytical approach to access the exact energy spectrum and wave functions of gated Bernal bilayer graphene (BBLG), with all the tight-binding parameters included. To tackle the broken mirror symmetry caused by a gated voltage (V_g) and interlayer interactions, we create a unitary transformation to reduce the Hamiltonian matrix of BBLG to a simple form, which can offer the analytical energy spectrum. The formulae generate the gated tunable energy bands and reveal that V_g changes the subband spacing, produces the oscillating bands, and increases the band-edge states. Then, we employ the analytical model to revisit the optical dipole matrix element and optical absorption spectra. In the absence of V_g, the anisotropic dipole matrix element exhibits a maximum around the point M and zero value along the high symmetry line ΓK in the first Brillouin zone. V_g effectively induces the nonzero dipole matrix element along the high symmetry line ΓK, which makes a significant contribution in the absorption spectra. Moreover, the application of V_g opens an optical gap and gives rise to a profound low-energy peak in the absorption spectra. The dependence upon the gated bias V_g for the location and height of this peak clearly emerges through the analytical model. Our exact analytical model can be further used to study the many-body effect and exciton effect on the electronic and optical properties of BBLG.

---

1 Introduction

Graphene, a hit topic since 2004,¹ having attracted a surge of interest for fundamental and experimental studies,^2–9 is a one atom-thick layer made up of carbon atoms arranged in a honeycomb lattice. Stemming from the specular geometrical structure, graphene is an extraordinary two dimensional material with many unique properties, e.g., linear energy dispersions crossing at the Dirac point, electron–hole symmetry, Klein tunneling, high mobility at room temperature, and the novel quantum Hall effect.^7–13 Its unusual electric properties make graphene a possible candidate for fabricating future electronic devices. Graphene is a zero band gap semiconductor due to linear energy dispersions crossing at the Dirac point. The presence of valence and conduction bands in its energy spectrum limits the on-off current ratios achievable and prevents the use of graphene in making field transistors. It is significant for the opening of a band gap in graphene to remove the limitation in using graphene material for semiconductor applications. Many breakthrough methods are proposed to generate a band gap in monolayer graphene, e.g., using molecular doping,^14,15 the application of strain,^16,17 and patterning a graphene sheet into a nanoribbon.¹⁸ On the other hand, a tunable band gap is induced in Bernal bilayer graphene (BBLG) through the application of a perpendicular electric field. BBLG is made up of two monolayer graphenes, held by weak Van der Waals forces in the AB-stacking. Experimental works demonstrate that a tunable band gap of up to 0.25 eV is achieved for electrically gated bilayer graphene by a variable external electric field.^19–22

The energy dispersions of the gated BBLG are usually explored within the tight-binding model.^23–31 The minimal model, considering only the intralayer and main interlayer interactions, is widely employed to study the energy dispersions of BBLG.^30,31 Around the Dirac points, the analytical low-energy dispersions of gated BBLG are obtained through the continuum model. A 4 × 4 Hamiltonian matrix is solvable and gives out the analytical low-energy dispersions. A full tight-binding Hamiltonian model takes into account all hopping integrals, including the additional skew interaction (γ₃) and electron–hole asymmetry (γ₄) in the model.^23,27 γ₃ is the trigonal warping and γ₄ is the origin of electron–hole asymmetry. The trigonal warping can significantly modify the physical properties of BBLG.^32–34 Cserti has shown that trigonal warping has a great effect on the minimal conductivity of BBLG.³² Moreover, γ₃ gives rise to the electron–hole asymmetry and distorts the circular equi-energy contour to a trigonal symmetry.^33,34 A full tight-binding Hamiltonian matrix is a 4 × 4 matrix and its four energy bands are usually obtained by numerical diagonalization. Many important physical properties do not directly emerge through the numerical calculation because of a lack of the exact energy spectrum. To fully and exactly describe the electric and optical properties in the full energy region, an analytical description of the energy dispersions of the four-band model is expected.

Here, we present a model to access the analytical form of the energy spectrum and eigenstates of gated BBLG based on the tight-binding framework and all the tight-binding parameters, γ₀, γ₁, γ₃ and γ₄, are considered too (see the inset of Fig. 1(a)). Unlike AA-stacked grapheme,^35,36 the skew interlayer interactions (γ₃) and gated bias (V_g) destroy the inversion symmetry of gated BBLG. The broken symmetry increases the difficulty in solving the eigenvalue problem. Through a rotation operator, we generate a new set of tight-binding basis functions to deal with the broken symmetry resulting from the γ₃ and V_g. The renormalized V_g-dependent intralayer and interlayer interactions are derived based on the new set of basis functions. Most importantly, the new basis functions exactly reduce the Hamiltonian matrix (H) into a band-storage matrix. As a result, the energy dispersions and wave functions are analytically solvable. Furthermore, the present work aims to provide full-frequency energy bands, which are beyond the effective mass approximation. The analytical form of the energy spectrum and eigenstates are further applied to study the optical absorption spectra and explore the origin of optical transition channels.

Fig. 1 (a) The energy spectra, obtained using eqn (11), are illustrated by the red dashed curves. The black curves are the energy spectrum by numerical diagonalization. The geometric structure of the AB-stacked bilayer graphene and the intralayer and interlayer interactions are shown in the inset. (b) The same plots as in (a) but with the gated voltage V_g = 200 meV. The inset exhibits the first Brillouin zone. M is the saddle point and K and K′ are the Dirac points.

The rest of the paper is organized as follows: we develop a model to derive the analytical form of the energy spectrum of gated BBLG in Section 2. Subsequently, the electronic properties are revisited in Section 3 through the analytical energy spectrum. In Section 4, the energy spectrum and eigenstates are used to explore the dipole matrix element and optical absorption spectra of gated BBLG. Finally, the conclusions are drawn in Section 5.

2 Theory and model

BBLG is a stack of two identical graphene layers with a AB-stacking, as shown in the inset of Fig. 1(a). On each graphene plane, two basic atoms, A and B, are arranged in a honeycomb lattice and the bond length between atoms A and B is b = 1.42 Å. The distance between two graphene layers is assumed to be 3.35 Å.³⁷ The lower or upper graphene layers are denoted by the index l = 1 or 2. Atoms A₂ and A₁ have the same (x, y) projection and they are denoted as the dimer sites. The projections of the other half, the B atoms, lie in the center of the hexagons in the adjacent sheets. BBLG is a periodic system in the xy plane and has four atoms, A₁, B₁, A₂ and B₂, in its primitive cell. The first Brillouin zone of BBLG, as illustrated in the inset of Fig. 1(b), is the same as that of graphene, which is a hexagon with points Γ located at (0, 0), K at , and M at .

The electronic properties of BBLG subjected to a perpendicular electric field are studied within the tight-binding method framework. The Hamiltonian equation of BBLG is HΦ = EΦ, where the Bloch function is Φ = c₁|A₁〉 + c₂|A₂〉 + c₃|B₁〉 + c₄|B₂〉, which is the linear combination of the periodic wave functions |A₁〉, |A₂〉, |B₁〉, and |B₂〉. HΦ = EΦ is then transformed into the matrix equation H |c〉 = E |c〉, here |c〉 = |c₁, c₂, c₃, c₄〉^T. The representation of the Hamiltonian matrix reads

where H is a 4 × 4 matrix. The elements H_AA, H_AB, H_BA, and H_BB are 2 × 2 matrices. The hopping integrals (γ_i) usually used to describe Bernal bilayer graphene are illustrated in the inset of Fig. 1(a). γ₀ is the intralayer interaction. The hopping integral, γ₁, is the interaction between the dimer sites A₁ and A₂. The trigonal warping integral, γ₃, is the skew interlayer hopping between atoms B₁ and B₂. γ₄ represents γ_A1,B2 and γ_B2,A1. The matrices H_AA, H_AB and H_BA are further expressed aswhere the on-site energy (V_g) is the electric potential caused by the external perpendicular electric field. is the structure factor, where k = (k_x, k_y) is the wavevector and b_l represents the nearest neighbor on the same graphene plane. The three neighboring atoms are located at , b₂ = (0, b) and , and thus .

To find the analytical energy dispersions of the 4 × 4 matrix above, we now construct the symmetrized basis functions (|ψ₁〉, |ψ₂〉, |ψ₃〉, |ψ₄〉)^T, which are linear combinations of the periodic wave functions |A₁〉, |A₂〉, |B₁〉, and |B₂〉. The symmetrized basis functions are

where the coefficients s^±₁ and s^±₂ are the components of the eigenvector of the matrix H_AA (eqn (2a)). The eigen-equation of H_AA iswhere the eigenenergy is . The related eigen-vectors (s^±₁, s^±₂)^T areSimilarly, the eigenenergy related to H_BB is .

The Bloch function, spanned by the new basis functions (|ψ₁〉, |ψ₂〉, |ψ₃〉, |ψ₄〉)^T, is Φ = d₁|ψ₁〉 + d₂|ψ₂〉 + d₃|ψ₃〉 + d₄|ψ₄〉. According to the symmetrized basis functions |ψ_i〉, the matrix equation has the form

where the matrix elements are [script letter H]_ij = 〈ψ_i|H|ψ_j〉.The diagonal matrix elements areThe off-diagonal terms areNotably, in the absence of the gated bias, V_g = 0, [script letter H]₁₃, [script letter H]₃₁, [script letter H]₂₄ and [script letter H]₄₂ are all equal to zero.³⁸ Moreover, a straightforward calculation shows that [script letter H]₂₂² + [script letter H]₂₃[script letter H]₃₂ = ε_BB² = V_g² + γ₃²|f|² because the sub-matrix
is the transformation of H_BB.

2.1 In the limit of γ₄ = 0

The exclusion of electron–hole asymmetry, γ₄ = 0, allows us to obtain a simple and exact analytical solution to the energy spectrum. The interactions [script letter H]₁₃, [script letter H]₃₁, [script letter H]₂₄ and [script letter H]₄₂ in H_red are all equal to zero as γ₄ = 0. The reduced Hamiltonian matrix (H_red) is a band-storage matrix and is expressed as follows:where [script letter H]₁₂ = [script letter H]₃₄ = γ₀|f_(k)|. The analytical energy dispersions, associated with H_red, readwhere
where [script letter H]₂₂² + [script letter H]₂₃[script letter H]₃₂ = γ₃²|f|² + V_g² is used. and the angle ϕ is . There are four branches of energy bands, which are in the sequence E(k)₊₊ > E(k)₊₋ > E(k)₋₋ > E(k)₋₊. The former tow, E(k)₊₊ and E(k)₊₋, are the energy dispersions associated with the conduction bands. The latter two, E(k)₋₋ and E(k)₋₊ are the energy spectra related to the valence bands.

At the point K, [scr B, script letter B] = γ₁²/2 + V_g² and [capital script C] = V_g²(γ₁² + V_g²). The state energies are exactly E(K) = ±V_g and . Around the point K and in the absence of the electric field, [scr B, script letter B] ≫ [capital script C] and [capital script C] ≈ (γ₀²|f_k|² − γ₁γ₃|f_k|)². The two branches of the low energy dispersions of BBLG are . They are the eigenenergies of the effective Hamiltonian matrix

2.2 The effect of γ₄ on the energy dispersions

A plainly analytical energy dispersion of the Hamiltonian matrix (eqn (6)), involving all hopping integrals, is in fact inaccessible. The energy dispersions are generally obtained by the numerical diagonalization method. To obtain a simply analytical solution to the energy spectrum, we treat γ₄ as a perturbation and neglect the matrix elements [script letter H]₁₃, [script letter H]₃₁, [script letter H]₂₄ and [script letter H]₄₂ because they make nonsignificant contributions to the energy spectrum. [script letter H]₁₃, [script letter H]₃₁, [script letter H]₂₄ and [script letter H]₄₂ are right dependent on the magnitude of the gated voltage (V_g) and interlayer interaction (γ₄). The magnitude of γ₄ is smaller than that of γ₀ and γ₃, so that [script letter H]₁₃, [script letter H]₃₁, [script letter H]₂₄ and [script letter H]₄₂ are negligible. γ₄ affects the matrix elements [script letter H]₁₂ and [script letter H]₃₄ in such a manner as and . The Hamiltonian equation with the band-storage Hamiltonian matrix is expressed as follows:The eigenenergies, determined by the criterion Det|[script letter H]| = 0, are the roots of the secular equation E⁴ − BE² − B′E + C = 0, where B = [script letter H]₁₁² + [script letter H]₂₂² + [script letter H]₂₃[script letter H]₃₂ + [script letter H]₁₂² + [script letter H]₃₄², B′ = ([script letter H]₁₁ + [script letter H]₂₂)([script letter H]₁₂² − [script letter H]₃₄²), and C = [script letter H]₁₂²[script letter H]₃₄² − [script letter H]₁₁[script letter H]₂₂([script letter H]₁₂² + [script letter H]₃₄²) + [script letter H]₁₁²([script letter H]₂₂² + [script letter H]₂₃[script letter H]₃₂). The analytical form of the eigenenergies looks so complicated and tedious just because the coefficient B′ ≠ 0. To acquire the simple analytical result, we further approximate B′ ≈ 0 because [script letter H]₁₂ ≈ [script letter H]₃₄. As a result, the eigenenergies, ε(k), are

To check out the accuracy of the afore-presented model, we first calculate energy dispersions ε(k)_±± at V_g = 0 in the high symmetry line K–Γ–M–K, shown by the dashed red curves in Fig. 1(a). There are four energy subbands ε(k)_±±. For comparison, energy dispersions, E_±±, obtained through the numerical diagonalization of the Hamiltonian matrix (eqn (1)), are also drawn by the black curves in Fig. 1(a). The tight binding parameters (γ_is) used to calculate the energy bands are: γ₀ = −3.12 eV, γ₁ = 0.38 eV, γ₃ = 0.280 eV, and γ₄ = 0.12 eV.³⁹ Obviously, the red curves are not identical to those in black in the higher energy region ε > 3.0 eV. The higher the energy (ε) is, the more discrepancy there is between the two curves. Plus, the energy spectrum ε in the presence of the gated voltage V_g = 200 meV are presented in Fig. 1(b). The applied gated voltage does not significantly affect the profile of energy bands in the high symmetry line K–Γ–M–K. The disagreement between the red and black curves reveals that the derived analytical formula can not exactly describe the energy spectra in the energy region E(eV) > 3.0 eV. Energy dispersions around the point K at different V_gs are shown in Fig. 2 for comparison. The calculation result illustrates that the presented analytical model can not replicate the exact energy dispersions around the point K. Notably, in the absence of the gate voltage, the presented analytical model ε_±± can not reproduce the analytical formula of BBLG at V_g = 0.³⁸

Fig. 2 The dashed curves are the low energy spectra at different V_gs acquired by eqn (11). For comparison, the energy spectra by numerical diagonalization are drawn in the black curves.

We here modify and improve the aforementioned analytical model to approach the exact energy dispersions at any V_g. In the absence of the gated voltage V_g = 0, the matrix elements [script letter H]₂₃ and [script letter H]₃₂ in eqn (10) play a minor role and can then be neglected.³⁸ E₊₊ and E₋₋ (E₊₋ and E₋₊) belong to the symmetrical states (the anti-symmetrical states). E₊₊ and E₋₋ depend only upon while E₊₋ and E₋₊ depend on .³⁸ We thus replace both [script letter H]₃₄ and [script letter H]₄₃ with [script letter H]₁₂ and [script letter H]₂₁ in eqn (10) during the calculation of E₊₊ and E₋₋ at any V_g. Then, we obtain the energy dispersions

where [scr B, script letter B] = [script letter H]₁₂² + ([script letter H]₁₁² + [script letter H]₂₂² + [script letter H]₂₃[script letter H]₃₂)/2 and [capital script C] = [script letter H]₁₂⁴ − 2[script letter H]₁₂²[script letter H]₁₁[script letter H]₂₂ + [script letter H]₁₁²([script letter H]₂₂² + [script letter H]₂₃[script letter H]₃₂). The other two energy bands E₊₋ and E₋₊ are acquired by utilizing to calculate [scr B, script letter B]^′ = [script letter H]₃₄² + ([script letter H]₁₁² + [script letter H]₂₂² + [script letter H]₂₃[script letter H]₃₂)/2 and [capital script C]^′ = [script letter H]₃₄⁴ − 2[script letter H]₃₄²[script letter H]₁₁[script letter H]₂₂ + [script letter H]₁₁([script letter H]₂₂² + [script letter H]₂₃[script letter H]₃₂). E₊₋ and E₋₊ are, respectively, expressed as

The analytical formulae (eqn (12a–d)) are now employed to calculate the energy dispersions of gated BBLG. We first present the analytical energy dispersions E(k)_±± at V_g = 0 in Fig. 3(a). E(k)_±± in the high-symmetry line K–Γ–M–K are drawn in red to be distinguished from those obtained through the numerical diagonalization of the Hamiltonian matrix in black. The red curves are exactly identical to the black ones except for the conduction bands in the high-symmetry line KM, where 3.0 eV < E < 3.5 eV (Fig. 3(b)). The energy spectrum in the presence of the gated voltage V_g = 200 meV is also exhibited in Fig. 3(c). The identity between the red and black curves reveals that the derived analytical formula can well describe the energy spectrum of the gated BBLG, except for the conduction bands in the energy region 2.6 eV < E < 3.4 eV (Fig. 3(d)). The results of the calculation show that the analytical model is able to provide us with the exact energy dispersions of the gated BBLG as E < 2.6 eV or E > 3.5 eV. The analytical energy dispersions and conjugated eigenstates are useful for exploring the electronic properties, such as energy dispersions, density of states (DOS), and absorption spectra.

Fig. 3 (a) The energy dispersions at V_g = 0 illustrated in the red dashed curves are obtained from the exact analytical model. The black curves are the energy spectra from numerical diagonalization. (b) The conduction bands around the point M. (c) The same plot as in (a) but for V_g = 200 meV. (d) The conduction bands around the point M.

3 Discussion: electronic properties

The characteristics of the low energy spectrum are explored. BBLG owns two pairs of parabolic conduction and valence bands, as shown in Fig. 4 in the dashed and solid curves, which are obtained, respectively, by the analytical formulae and the numerical diagonalization method. The dashed curves are in good agreement with the solid curves. With the comparison shown above, the exact solution to the low energy spectrum of BBLG emerges through our analytical model. While the high energy bands arouse some focus of our discussion, we are more interested in the lowest-energy electronic properties. At V_g = 0, E(k)₊₋ and E(k)₋₊, the first pair located near the chemical potential μ = 0, exhibit a tiny overlap. E(k)₊₋ and E(k)₋₋ are asymmetrical about μ = 0, i.e., electron–hole asymmetry, which is caused by the interlayer interaction γ₄. The minimum of E(k)₊₋ and the maximum of E(k)₋₊ are located at almost the same wavevector (∼the K point). E(k)₊₋ and E(k)₋₋ along the high symmetry lines Γ–K and M–K are strongly anisotropic, which originates in the vanishing of [script letter H]₂₃ and [script letter H]₃₂ in eqn (12a–d) along the line Γ–K.

Fig. 4 The dashed and solid curves are energy dispersions obtained from the exact analytical formulae and numerical diagonalization. For comparison, the exact energy dispersions are illustrated by the colored curves.

The gated voltage (V_g) has a great influence on the low energy spectrum. First, the low energy subbands are significantly modified by V_g and change from monotonically parabolic dispersions into oscillating ones, referred to as a Mexican hat. Then, the asymmetry between the conduction E(k)₊₋ and the valence E(k)₋₊ bands along the Γ–K direction about the chemical potential is also enhanced. In addition to the K point, two extra band-edge states emerge in all directions. There exists a band-edge state with the minimum (maximum) energy along the Γ–K direction. The energy difference between the lowest band-edge state of E₊₋(k) and the highest one of E₋₋(k) decides the band gap, E_g. The locations of the band-edge state are determined by and . By taking the main interactions, γ₀ and γ₁, into consideration, E₊₋ (E₋₋) has a minimum (maximum) at

where the Fermi velocity (v) is ħv = 3bγ₀/2. Therefore, the energy gap is
The size of the band gap, E_g = E₋₊ (k_c) − E₋₋ (k_v), at first grows rapidly but then gradually declines with the increase in V_g. Another band-edge state along the KM direction belongs to the saddle point. Such states are the critical points in the energy-wavevector space and thus have a high DOS.

The DOS is useful for understanding the essential physical properties of BBLG. DOS is defined as

where η₁, η₂ = ±, ± represents the subband index. The broadening energy width (Γ) is set as 3.0 meV. The full density of state, D(ω), is the summation of D_η1,η2(ω), resulting from each subband, E_η1,η2. The D(ω) and D_η1,η2(ω) of BBLG at V_g = 0 are illustrated in Fig. 5(a). The DOS exhibits logarithmic peaks, originating in the saddle point M around ω = ±3 eV, which directly reflects the main features of the energy dispersions. The detailed structure of the DOS at low energy is presented in Fig. 5(b). D₊₋(ω) (D₋₊(ω)) makes contributions to D(ω) as 0 < ω < γ₁ (−γ₁ < ω < 0). BBLG at V_g = 0 is a semi-metal due to the non-vanished DOS at ω = 0. The gated voltage (V_g) can open a band gap and induce two van Hove singularities in the energy region −γ₁ < ω < γ₁, as shown in Fig. 5(c–e), where the band gap is labeled by vertical red bars and van Hove singularities are indicated by the blue arrows. One can modulate the band gap in BBLG using an external electric field. The increase in V_g not only enlarges the size of the band gap, but also enhances the height of van Hove singularities. The predicted size of the band gap and locations of van Hove singularities can be verified by experimental measurements.

Fig. 5 (a) DOS D_η1η2, associated with each subband E_η1η2, of BBLG at V_g = 0. (b)-(e) are the low-energy DOS at V_g = 0, 100, 200, and 300 meV. The vertical red bars indicate the energy band gap. The blue arrows label the van Hove singularities in the DOS.

4 Absorption spectra

Utilizing the analytical model, we revisit the low-energy absorption spectra of gated BBLG, which has been studied through numerical methods.^40–42 The absorption function of bilayer graphene at zero temperature (T = 0) is given by³⁸where Γ is the broadening parameter owing to various nonradiative processes, f (i) denotes the final (initial) state and F_f (F_i) is the Fermi-Dirac distribution function. E_ex = E_f(k) − E_i(k) is the excitation energy. , the dipole matrix element, is the velocity operator between the initial and final wave functions Φ_i(k) and Φ_f(k). The absorption spectrum A_x(ω) originates in electronic transitions that correspond to excitations from the occupied valence bands to the unoccupied conduction bands, excited by the electromagnetic field. In this work, the direction of the electromagnetic field is assumed to be parallel to the x-axis. Within the gradient approximation,^43,44 the dipole matrix element , where is derived from the gradient of the Hamiltonian operator, . Because of the zero momentum of photons, the optical selection rule is Δk = 0. The optical transition channels are determined by the thermal factor F_F − F_i and the dipole matrix element.

The dipole matrix element, by inserting the tight-binding wave functions Φ = c₁|A₁〉 + c₂|A₂〉 + c₃|B₁〉 + c₄|B₂〉, is expressed as

where and H_lm is the element of the Hamiltonian matrix (eqn (1)). After the action of the rotation operator (eqn (3)), the dipole matrix element is then changed intowhere the analytical form of the matrix element [scr V, script letter V]_lm is given explicitly in the Appendix. |d^i,f₁, d^i,f₂, d^i,f₃, d^i,f₄〉 are the initial (final) eigenstates associated with the Bloch function, eqn (6). The eigenstate corresponding to each energy dispersion (E_±±) can be exactly specified asThe analytical expression for the dipole matrix element can be used to efficiently evaluate the magnitude of M_fi and the optical absorption spectra A(ω). M_fi depends on energy spectrum, velocity operator, and eigenstates.

The square of the absolute value of dipole matrix element |M_fi(k)|² in the high-symmetry line Γ–K–M for the transition E₋₋ → E₊₋ at different V_gs is drawn in Fig. 6, where the solid (dashed) curves are calculated with the energy spectrum and the eigenstates are obtained through the numerical method (analytical model). The dashed curves are exactly in agreement with the solid ones along the high-symmetry line ΓK, while they exhibit some discrepancy around the point M. In the absence of the gated voltage (V_g = 0), |M_fi(k)|² exhibits a highly anisotropic characteristic with a maximal value around the point M in the high-symmetry line Γ–K–M. |M_fi| is equal to zero in the line Γ–K, along which the structure factor f(k) is a real number and the matrix elements [script letter H]₂₃ and [script letter H]₃₂ in eqn (6) are equal to zero. The 4 × 4 Hamiltonian matrix is decomposed into two 2 × 2 blocks. The eigen-vector corresponding to the final (initial) state is |d^f〉 = |d₁, d₂, 0, 0〉 (|dⁱ〉 = |0, 0, d₃, d₄〉). The six matrix elements [scr V, script letter V]₁₃, [scr V, script letter V]₂₃, [scr V, script letter V]₂₄, [scr V, script letter V]₃₁, [scr V, script letter V]₃₂, and [scr V, script letter V]₄₂ are all equal to zero as Im[f(k)] = 0 and V_g = 0. The straightforward calculation exhibits that the magnitude of the dipole matrix element vanishes, i.e., |M_fi| = 0.

Fig. 6 (a)-(d) The solid curves are the square of the absolute values of the optical matrix |M_fi|, obtained by the numerical diagonalization method. The same plots obtained using analytical formulae are given in the dashed curves.

In the presence of the gated voltage, |M_fi|² shows different aspects, as illustrated in Fig. 6(b–d). V_g alters the Hamiltonian matrix elements, energy spectrum and eigenstates |dⁱ〉 and |d^f〉. Moreover, the V_g-dependent [scr V, script letter V]_lms are also modified by the gated voltage. As a result, at V_g = 100 meV, |M_fi|² (blue solid curve in Fig. 6(b)) around the point M exhibits a different aspect and magnitude. The maximum is located in the vicinity of the point M. The increase in V_g does not alter the profile of |M_fi|² around the point M. Most importantly, the application of V_g enhances |M_fi|² in the high-symmetry line KΓ. V_g introduces the non-zero matrix elements [script letter H]₂₃ and [script letter H]₃₂ in eqn (6), which mixes the symmetrical and anti-symmetrical wave functions, i.e., the alternation of the eigenstates. Moreover, V_g also produces the six matrix elements [scr V, script letter V]₁₃, [scr V, script letter V]₂₃, [scr V, script letter V]₂₄, [scr V, script letter V]₃₁,[scr V, script letter V]₃₂, and [scr V, script letter V]₄₂, being unequal to zero. As a result, V_g intensifies |M_fi|² in the symmetry line KΓ. Such a characteristic strongly modifies the low energy absorption spectra in the presence of V_g.

The low-energy absorption spectra A(ω) of AB BLG reflect the characteristics of the energy dispersions. There are six possible interband excitations resulting from four branches of energy dispersion. The optical excitation from the subband E₋₊ to E₋₋ and transition between E₊₋ and E₊₊ are forbidden by the thermal factor F_f − F_i. Four allowed excitations, A₁(ω), A₂(ω), A₃(ω) and A₄(ω), are shown by the vertical lines in the inset on the left panel of Fig. 7(a). A₁(ω) stems from the transition between the highest valence subband (E₋₋) and the top conduction subband (E₊₊), and A₂(ω) originates in the excitation E₋₊ → E₊₋. Moreover, the transition from the subband E₋₋ to E₊₋ generates the spectra A₃(ω) and the excitation between the subbands E₋₊ and E₊₊ leads to A₄(ω). A₁(ω) [A₂(ω)] shows a discontinuous structure at ω ≈ 380 meV, as illustrated in Fig. 7(a). Such an optical gap originates in the interband transition between the E₋₋ and E₊₊ subbands [E₋₊ and E₋₊ subbands]. A₄(ω), resulting from the transition from the E₋₊ to E₊₊ subband, exhibits an absorption edge at ω ≈ 800 meV. There is no optical gap found in the spectra A₃(ω) (red curve) because of the overlap between subbands E₋₋ and E₊₋. As a result of the superposition of the sub-spectra, the total spectra exhibit two discontinuities at ω ≈ 380 and ω ≈ 800 meV, and no optical gap is found in A(ω).

Fig. 7 (a) A₁(ω), A₂(ω), A₃(ω), and A₄(ω), the sub-optical spectra, are shown by the blue, green, red and brown curves. Four allowed optical transitions are indicated by the vertical lines in the inset in the left panel. (b) The low-energy absorption spectra (A_i(ω)) at V_g = 0, 100, 200, and 300 meV are drawn in red, green, brown and blue. (c) The low energy absorption spectra with different tight-binding parameters.

The application of V_g has a great effect on the low-energy absorption spectra A₃(ω), as depicted in Fig. 7(b). V_g not only opens a band gap, but also gives rise to a sharp peak in A₃(ω). In the presence of V_g, the steep rise of the spectra, due to the transition between the band-edge states of E₋₋ and that of the E₊₋ subband in the ΓK direction, can be used to determine the size of the energy gap, E_g. The energy gaps are equal to ω = 0.17, 0.25, and 0.28 eV at V_g = 100, 200, and 300 meV. The first absorption peak is a compound one, consisting of a shoulder on the low energy side and a main peak. The former (shoulders) result from the transition between the E₋₊ and E₊₋ subbands, denoted by the vertical red arrow in Fig. 4. The excitation indicated by the vertical blue arrow in Fig. 4 brings about the main peak. At V_g = 100, 200, and 300 meV, the main peaks are located at ω = 0.18, 0.30, and 0.35 eV. The main peak is in the logarithmic form, stemming from the saddle point in energy dispersions in the high-symmetry line MK. The frequency of the main peak depends on the magnitude of V_g. These peaks make a blue-shift with the increase in V_g. The height of the first peak also enhances with the increase in V_g because the DOS around the band-edge state of E₋₊ (E₊₋) and |M_fi|² in the high-symmetry line ΓK are enhanced by V_g, as shown in Fig. 5 and 6.

How the tight-binding parameters (γ_is) influence the low energy absorption spectra (A(ω)) is deliberated. The dashed brown curve in Fig. 7(c) represents A(ω) at V_g = 200 meV, through the analytical model, including all the tight-binding parameters, γ₀, γ₁, γ₃, and γ₄. The close of the interlayer interaction γ₄ has no effect on A(ω), as shown by the cyan curve in Fig. 7(c). According to eqn (7b), the magnitude of gradually declines with an increase in V_g. As a result, the absorption spectra at V_g = 200 meV are not affected by the interlayer interaction γ₄. A(ω), illustrated by the red curve, exhibits a larger band gap (E_g = 270 meV) and a sharp peak located at ω = E_g as the main intralayer interaction (γ₀) and interlayer interaction (γ₁) are taken into consideration. Moreover, A(ω) presents a simple feature and weaker intensity of absorption relative to those in the brown curve. The chief cause is that the close of the interlayer interaction γ₃ not only changes the energy dispersions and wave functions, but also turns off some optical transition channels.

A comparison of the study results with those obtained through ab initio calculations is made. The analytical formulae can exactly describe four energy dispersions and generate the energy bands with “Mexican hat” structures. Our research demonstrates that the gated voltage, V_g, alters the subband spacing, changes the energy band gap, produces the oscillating bands, and induces more band-edge states. The aforementioned electronic properties, e.g., band feature, tunable band gap, and electric-field-modified oscillating bands, are also produced through the ab initio calculations.⁴¹ That is, the electronic properties derived from the analytical model are in qualitative agreement with those given by the ab initio calculations. To be more specific, we quantitatively compared the sizes of the gated tunable band gaps acquired by the two methods. The calculated tunable band gaps, shown in Fig. 8 in the cross symbols, are taken from the ab initio calculations.⁴⁵ The band gap first increases monotonically as the strength of the electric field increases. Then, the energy gaps saturate to a value E_g ≈ 270 meV when the field strength |F| is greater than 0.2 V Å⁻¹. We now evaluate E_g by employing the analytical model for γ₃ = γ₄ = 0 (black curve) and γ₃ ≠ 0 and γ₄ ≠ 0 (red curve). The red curve exhibits that E_g first increases sublinearly when the strength of the electric field (|F| = 2V_g/c and c is the layer distance) is less than 0.05 V Å⁻¹. Meanwhile, the energy gaps approach a value of E_g ≈ 300 meV when |F| > 0.2 V Å⁻¹. The behavior of E_g predicted by the analytical method (red curve) is similar to that of the ab initio method. The comparison between the black and red curves illustrates that when the field strength is weaker than 0.05 V Å⁻¹, the tunable band gaps are independent of the hopping integral γ₃. After that, the hopping integral γ₃ reduces the size of E_g.

Fig. 8 The cross symbols are band gaps acquired through the ab initio calculations (data are taken from APL 98, 2011, 263107). E_g is obtained from the analytical model for γ₃ = γ₄ = 0 (black curve) and γ₃ ≠ 0 and γ₄ ≠ 0 (red curve).

5 Conclusions

We develop a method to access the analytical form of the energy spectrum and eigenstates of gated BBLG within the tight-binding model with all the tight-binding parameters included. The trigonal warping (γ₃) and applied gated bias (V_g) destroy the mirror symmetry of BBLG. To deal with the broken symmetry caused by γ₃ and V_g, a new set of tight-binding basis functions is constructed through a rotation operator. The renormalized intralayer and interlayer interactions exhibit a gated bias dependence. Most importantly, the Hamiltonian matrix (H) is transformed into a band-storage matrix in the subspace spanned by the new basis functions. The formulae of the energy dispersions and wave functions are then obtained analytically. γ₄, the electron–hole asymmetry, is treated as a perturbation and included in the analytical formula. The electronic properties, density of states, optical dipole matrix element and optical absorption spectra are explored by employing this analytical formula. The study shows that V_g can open an optical gap and produce a profound low-energy peak in the absorption spectra. The dipole matrix element (M) and absorption spectra are strongly modified by V_g. Finally, the characteristics of the absorption spectra, e.g., the location and height of the low-energy peak, are significantly presented and explored through the analytical model.

Appendix

The element [scr V, script letter V]_lm is the operator sandwiched by the symmetrized basis functions 〈ϕ_l| and |ϕ_m〉, i.e. , where H is the Hamiltonian representation (eqn (1)). The analytical form of each [scr V, script letter V]_lm is listed as follows.

[scr V, script letter V]₁₁ = 0,

[scr V, script letter V]₁₄ = 0,

[scr V, script letter V]₂₁ = [scr V, script letter V]^*₁₂,

[scr V, script letter V]₃₁ = [scr V, script letter V]^*₁₃,

[scr V, script letter V]₃₂ = [scr V, script letter V]^*₂₃,

[scr V, script letter V]₃₃ = −[scr V, script letter V]₂₂,

[scr V, script letter V]₄₂ = [scr V, script letter V]^*₂₄,

[scr V, script letter V]₃₄ = [scr V, script letter V]^*₄₃,

[scr V, script letter V]₄₄ = 0.

Acknowledgements

The authors gratefully acknowledge the support of the Taiwan National Science Council under the Contract no. NSC 101-2112-M-165-001-MY3.

References

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science, 2004, 306, 666.
2. Y. B. Zhang, Y. W. Tan, H. L. Stormer and P. Kim, Nature, 2005, 438, 201.
3. V. P. Gusynin and S. G. Sharapov, Phys. Rev. Lett., 2005, 95, 146801.
4. A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov and A. K. Geim, Rev. Mod. Phys., 2009, 81, 109.
5. N. M. R. Peres, Rev. Mod. Phys., 2010, 82, 2673.
6. M. J. Allen, V. C. Tung and R. B. Kaner, Chem. Rev., 2010, 110, 132.
7. D. S. L. Abergela, V. Apalkovb, J. Berashevicha, K. Zieglerc and T. Chakrabortya, Adv. Phys., 2010, 59, 261.
8. K. S. Novoselov, V. I. Fal'ko, L. Colombo, P. R. Gellert, M. G. Schwab and K. Kim, Nature, 2012, 490, 192.
9. N. O. Weiss, H. Zhou, L. Liao, Y. Liu, S. Jiang, Y. Huang and X. Duan, Adv. Mater., 2012, 24, 5782.
10. K. S. Novoselov, E. McCann, S. V. Morozov, V. I. Fal'ko, M. I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin and A. K. Geim, Nat. Phys., 2006, 2, 177.
11. F. Guinea, M. I. Katsnelson and A. K. Geim, Nat. Phys., 2010, 6, 30.
12. S. Das Sarma, S. Adam, E. H. Hwang and E. Rossi, Rev. Mod. Phys., 2011, 83, 407.
13. D. A. Abanin, B. E. Feldman, A. Yacoby and B. I. Halperin, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 115407.
14. M. Kim, N. S. Safron, C. Huang, M. S. Arnold and P. Gopalan, Nano Lett., 2012, 12, 182.
15. B. M. Wong, S. H. Ye and G. O'Bryan, Nanoscale, 2012, 4, 1321.
16. G. Cocco, E. Cadelano and L. Colombo, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 241412.
17. Z. H. Ni, T. Yu, Y. H. Lu, Y. Y. Wang, Y. P. Feng and Z. X. Shen, ACS Nano, 2008, 2, 2301.
18. M. Y. Han, B. Ozyilmaz, Y. Zhang and P. Kim, Phys. Rev. Lett., 2007, 98, 206805.
19. Y. Zhang, T. T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen and F. Wang, Nature, 2009, 459, 820.
20. A. B. Kuzmenko, I. Crassee, D. van der Marel, P. Blake and K. S. Novoselov, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 165406.
21. T. Thiti and J. H. Pablo, Phys. Rev. Lett., 2010, 105, 166601.
22. M. T. Allen, J. Martin and A. Yacoby, Nat. Commun., 2012, 3, 934.
23. C. L. Lu, C. P. Chang, Y. C. Huang, R. B. Chen and M. L. Lin, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 144427.
24. M. Orlita, C. Faugeras, J. M. Schneider, G. Martinez, D. K. Maude and M. Potemski, Phys. Rev. Lett., 2009, 102, 166401.
25. E. V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R. Peres, J. M. B. Lopes dos Santos, J. Nilsson, F. Guinea, A. K. Geim and A. H. Castro Neto, Phys. Rev. Lett., 2007, 99, 216802.
26. E. V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R. Peres, J. M. B. Lopes dos Santos, J. Nilsson, F. Guinea, A. K. Geim and A. H. Castro Neto, J. Phys.: Condens. Matter, 2010, 22, 175503.
27. Y. H. Ho, Y. H. Chiu, D. H. Lin, C. P. Chang and M. F. Lin, ACS Nano, 2010, 4, 1465.
28. E. McCann and M. Koshino, Rep. Prog. Phys., 2013, 76, 056503.
29. R. Côte and M. Barrette, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 245445.
30. H. Cheraghchi and F. Adinehvand, J. Phys.: Condens. Matter, 2014, 26, 015302.
31. F. Kadi and E. Malic, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 045419.
32. J. Cserti, A. Csordas and G. David, Phys. Rev. Lett., 2007, 99, 066802.
33. M. Koshino, New J. Phys., 2009, 11, 095010.
34. L. M. Zhang, M. M. Fogler and D. P. Arovas, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 075451.
35. C. P. Chang, J. Appl. Phys., 2011, 110, 013725.
36. C. P. Chang, Carbon, 2013, 61, 209.
37. J. C. Charlier, J. P. Michenaud and X. Gonze, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 4531.
38. C. P. Chang, J. Appl. Phys., 2012, 111, 103711.
39. K. Nakao, J. Physiol. Soc. Jpn., 1976, 40, 761.
40. C. L. Lu, C. P. Chang, Y. C. Huang, R. B. Chen and M. L. Lin, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 144427.
41. H. Min, B. Sahu, S. K. Banerjee and A. H. MacDonald, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 155115.
42. L. Yang, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 155445.
43. J. G. Johnson and G. Dresselhaus, Phys. Rev. B: Solid State, 1973, 7, 2275.
44. F. L. Shyu, C. P. Chang, R. B. Chen, C. W. Chiu and M. F. Lin, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 67, 045405.
45. B. R. Wu, Appl. Phys. Lett., 2011, 98, 263107.

---