Electronic and Vibrational Properties of Adsorbed and Embedded Graphene and Bigraphene with Defects

1. Introduction

The perennial interest in studying the physical properties of nanofilms has increased substantially over the last few years due to the development of nanotechnologies and the synthesis of new compounds – especially those based on carbon, which are extremely interesting for both fundamental research and potential applications.

An important feature of carbon nanofilms (including those with defects) is a close relation between the electronic and phonon properties, which is exhibited, for example, in the graphene-based systems with superconducting properties [1,2].

It is well known that graphene monolayers cannot exist as planar objects in the free state, because in flat 2D-crystals the mean-square amplitudes of the atoms in the direction normal to the layer plane diverge even at $$T=0$$ (see, e.g., [3]). So we can study and practically apply only such graphene, which is deposited on a certain substrate providing the stability of the plane carbon nanofilms (see, e.g., [4-6]). Only small flakes can be detached from the substrate and these flakes immediately acquire a corrugated shape [7]. When studying the electronic properties of graphene a dielectric substrate is often used. The presence of the substrate greatly increases the occurrence of various defects in graphene and carbon nanofilms. Our investigations make it possible to predict the general properties of phonon and electron spectra for graphene and bigraphene containing different defects.

This chapter consists of three sections: first section is devoted to the calculation of local discrete levels in the electron spectra of graphene with different defects. In the second section we describe the electronic properties of bilayered graphene and, finally, the third section deals with the influence of defects on electron spectra of bigraphene.

2. Impurity levels in the electron spectra of graphene

The exceptionality of graphene is manifested in the phonon and electron properties. Graphene is a semimetal whose valence and conduction bands touch at the points K and K’ of the Brillouin zone [8,9]. In the pure graphene unique electronic properties are manifested by the charge carriers behaving as massless relativistic particles - the dependence of energy on the momentum is linear rather than - as in ordinary solids - quadratic. Thus, the lower-dimensionality affects the formation of phonon localized states [10] and also the formation of localized states in the electronic spectrum. Absent gap between the valence and conduction bands is a consequence of the symmetry between two equivalent sublattices in graphene [11]. Presence of impurities lowers the symmetry. The influence of vacancies placed into one of the graphene sublattices was investigated in [12], where it was shown that the equivalence of the sublattices is broken.

In this section we describe the characteristics of localized and local states present in graphene due to the impurities of nitrogen and boron, respectively. The presence of the substrate greatly increases the possibility to introduce various defects into graphene. For example, in the graphene deposited on silicon, vacancies can occur [13, 14], whereas in graphite (a set of weakly interacting graphene monolayers) vacancies heal and form a stacking fault with local fivefold symmetry axis [15]. Impurity atoms embedded in graphene may lead to the appearance of impurity states outside the band of quasi-continuous spectrum. At low impurity concentrations (when impurity is considered as an isolated defect) these states appear in the form of local discrete levels (LDL).

Although such levels in various quasiparticle spectra have been known and studied over 60 years, an adequate description is still absent, even in the harmonic approximation for sufficiently realistic models of the crystal lattice. The dependence of the appearance conditions and characteristics of LDL on the parameters of a perfect lattice and defect was identified only in the most general terms. However, LDL may be used as an important source of information about the defect structure and force interactions in real crystals. To extract such information it is useful to have analytical expressions that relate main characteristics of LDL to the parameters of both the defect and the host lattice.

Here we present the results of our calculations and analyses of the characteristics of the electronic local discrete levels for substitutional impurities in graphene, especially for a boron substitutional impurity, using an analytical approximation based on the Jacobi matrices method [16, 17].

The fact that the charge carriers in graphene are formally described by the Dirac equation and not by the Schrödinger equation is due to the symmetry of the crystal lattice of graphene, which consists of two equivalent carbon sublattices. Electronic subbands formed by the symmetric and antisymmetric combinations of wave functions in the two sublattices intersect at the edge of the Brillouin zone, which leads to a cone-shaped energy spectrum near the K and K´ points of the first Brillouin zone. The electrons obey the linear dispersion law (in ordinary metals and semiconductors the dispersion law is parabolic).

The electronic spectrum of graphene can be described by a strong coupling approximation, and it is sufficient to consider the interaction between nearest neighbors only (see, e.g., [5,6,18-20]). The corresponding Hamiltonian is

$$\widehat{H}={\displaystyle \sum _i{\epsilon }_i|i〉}〈i|-{\displaystyle \sum _{i,j}{J}_{ij}}|i〉〈j|$$ E1

where i and j are the labels of the nodes of the two-dimensional lattice, $${\epsilon }_i$$ is the energy of electron at node i, and $$J_{ij}$$ is the so-called overlap integral.

Curve 1 in Figure 1 shows the density of electronic states of graphene as calculated using the method of Jacobi matrices [16, 17]. In a perfect graphene the local Green's function $$G\left(\epsilon ,i\right)=〈i|{\left(\epsilon \widehat{I}-\widehat{H}\right)}^{-1}|i〉$$ coincides with the total Green’s function $$G\left(\epsilon \right)=\underset{N\to \infty }{\lim }\frac{1}{N}{\displaystyle \sum _{i=1}^N〈i|{\left(\epsilon \widehat{I}-\widehat{H}\right)}^{-1}|i〉}$$ because of the physical equivalence of the atoms of both sublattices. Peculiarity of the density of states at $$\epsilon =\epsilon \left(K\right)$$ (the value $$\epsilon \left(K\right)$$ corresponds to the Fermi energy $${\epsilon }_F$$ in graphene) determines the behavior of the real part of the Green's function near $${\epsilon }_F$$ . For a wide class of perturbations caused by defects we can find, using the Lifshitz equation [21], quasilocalized states in the interval $$\left[-\epsilon \left(M\right),\epsilon \left(M\right)\right]$$ (in this model $$\epsilon \left(M\right)=J$$ ). This equation, which determines the energy of these states, can be written as (see, e.g., [3,17])

$$\mathrm{Re}G\left(\epsilon \right)=S\left(\epsilon ,{\Lambda }_{ik}\right)$$ E2

where the $$S\left(\epsilon ,{\Lambda }_{ik}\right)$$ function is determined by the perturbation operator $$\widehat{\Lambda }$$ ( $${\Lambda }_{ik}$$ are matrix elements of this operator on defined basis).

The local spectral densities $$\rho \left(\epsilon ,i\right)\equiv \frac{1}{\pi }\underset{\gamma \downarrow 0}{\lim }\mathrm{Im}G\left(\epsilon +i\gamma ,i\right)$$ of impurity atoms are calculated in [6]. For an isolated substitutional impurity with the energy $${\epsilon }_0=\tilde{\epsilon }$$ of the impurity node $$i=0$$ and with the overlap integral $$J_{i0}=\left(1+\eta \right)J$$ , the function $$S\left(\epsilon ,\tilde{\epsilon },\eta \right)$$ has the form

Figure 1.

Electronic density of states of perfect graphene (curve 1) and the local density of states for an isolated boron substitutional impurity (curve 2).

$$S\left(\epsilon ,\tilde{\epsilon },\eta \right)=\frac{{\left(1+\eta \right)}^2}{\tilde{\epsilon }+\epsilon \eta \left(2+\eta \right)}$$ E3

For a nitrogen impurity $$\tilde{\epsilon }-\epsilon \left(K\right)\approx -0.525J$$ and $$\eta \approx -0.5$$ (according to [4]). As shown in [22], equation (2) has a solution for both interval $$\left[-\epsilon \left(M\right),\epsilon \left(K\right)\right]$$ and interval $$\left[\epsilon \left(K\right),\epsilon \left(M\right)\right]$$ .

The local density of states of the nitrogen substitutional impurity calculated in [11] has quasi-local maxima in both intervals. For an boron substitutional impurity ( $$\eta \approx 0.5$$ ) [6], quasi-localized states are absent in the $$\left[-\epsilon \left(M\right),\epsilon \left(M\right)\right]$$ interval [12]. Figure 2 shows the graphical solution of the Lifshitz equation (2) for a given impurity atom. In this case the Lifshitz equation has no solutions in interval $$\left[-\epsilon \left(M\right),\epsilon \left(M\right)\right]$$ (corresponding dependences $$S\left(\epsilon \right)$$ are shown as curves 3 in Figure 2). The local Green's function of the boron impurity (curve 2 in Figure 1) has two poles outside the band of quasi-continuous spectrum, which are called local discrete levels and which are also solutions of equation (2). As is clearly seen in Figure 1 the area under the curve 2 is smaller (by the sum of the residues at these poles) than the area under the curve 1.

Local discrete levels can be an important source of information about defective structure and force interactions in real crystals. To extract this useful information we should have analytical expressions that relate the main characteristics of LDL (primarily their energy) to the parameters of the defect and the host lattice.

Figure 2.

Graphical solution of equation (2) for boron substitutional impurity in graphene. Curve 1 is the electronic density of states of ideal graphene, curve 2 is the corresponding real part of the Green's function. Curves $$1^{\prime }$$ and $$2^{\prime }$$ are the “approximations of two moments” of these functions. Curve 3 represents the function $$S\left(\epsilon \right)$$ .

Such expressions were obtained in [23] for localized vibrations in the phonon spectrum of a three-dimensional crystal. Authors proposed an analytical approximation of the basic characteristics of local vibrations based on the rapid convergence of the real part of the Green's function outside the band of quasi-continuous spectrum using the method of Jacobi matrices [16,17].

Let us, briefly, to the extent necessary to understand the use of the classification of the eigenfunctions of Hamiltonian (1), to present the basics of the method of Jacobi matrices. This method allows, without finding the dispersion laws, to calculate directly the local partial Green's functions of the system, corresponding to the perturbation of one or more atoms. This perturbation is described by the so-called generating vector $${\overrightarrow{h}}_0\in H$$ , where H is the space of electronic excitations of atoms. Its dimension is $$qN$$ , where N is the number of atoms in the system, and q is the dimension of the displacement of a single atom (q = 1, 2, 3). Vectors of the space H are denoted by an arrow above the symbol, and “ordinary” q-dimensional vectors are in bold italics.

If, using the generating vector (p is the number of excited atoms) and the Hamiltonian (1), we construct the sequence $${\left\{{\widehat{H}}^n{\overrightarrow{h}}_0\right\}}_{n=0}^{\infty }$$ , then the linear envelope covering the vectors of this sequence forms, in the H space, a cyclic subspace invariant to the operator $$\widehat{H}$$ . This subspace contains, within itself, all the atomic displacements generated by the vector $${\overrightarrow{h}}_0$$ . The corresponding partial Green's function is determined as a matrix element $$G_{00}\left(\lambda \right)\equiv \left({\overrightarrow{h}}_0,{\left[\epsilon \widehat{I}-\widehat{H}\right]}^{-1}{\overrightarrow{h}}_0\right)$$ , where $$\lambda \equiv {\omega }^2$$ is the eigenvalue. Quantity $$\rho \left(\lambda \right)\equiv {\pi }^{-1}\mathrm{Im}{G}_{00}\left(\lambda \right)$$ is called the spectral density generated by the initial displacement $${\overrightarrow{h}}_0$$ . In the basis $${\left\{{\overrightarrow{h}}_n\right\}}_{n=0}^{\infty }$$ which is obtained by the orthonormalization of the sequence $${\left\{{\widehat{H}}^n{\overrightarrow{h}}_0\right\}}_{n=0}^{\infty }$$ , the operator (1) is represented in the form of a tridiagonal Jacobi matrix (or J-matrix). This matrix has a simple spectrum, what greatly simplifies finding the partial Green's functions and spectral densities. As can be seen, this method does not use explicitly the translational symmetry of the crystal, making it extremely effective for treating systems in which such symmetry is broken. The method of Jacobi matrices is particularly effective for treating systems with a simply connected quasi-continuous band of spectrum D. In this case, with increasing rank of the J-matrix ( $$n\to \infty$$ ), its diagonal elements $$a_n$$ converge to a corresponding to the middle of the bandwidth D, and nondiagonal elements b _n converge to b corresponding to the one-quarter of the bandwidth D.

For the local Green's function (LGF), corresponding to the excitations of one or more atoms, which are determined by the generating vector $${\overrightarrow{h}}_0$$ , we get following expression using the J-matrix method

$$G\left(\epsilon ,{\overrightarrow{h}}_0\right)=\left({\overrightarrow{h}}_0,{\left[\epsilon \widehat{I}-\widehat{H}\right]}^{-1}{\overrightarrow{h}}_0\right)=\underset{n\to \infty }{\lim }\frac{Q_n\left(\epsilon \right)-b_{n-1}{Q}_{n-1}\left(\epsilon \right)K_{\infty }\left(\epsilon \right)}{P_n\left(\epsilon \right)-b_{n-1}{P}_{n-1}\left(\epsilon \right)K_{\infty }\left(\epsilon \right)}$$ E4

where $$\widehat{I}$$ is the unit operator and polynomials $$P_n\left(\epsilon \right)$$ and $$Q_n\left(\epsilon \right)$$ are determined by the following recurrence relations

$$b_n{\left\{P,Q\right\}}_{n+1}\left(\epsilon \right)=\left(\epsilon -a_n\right){\left\{P,Q\right\}}_n\left(\epsilon \right)-b_{n-1}{\left\{P,Q\right\}}_{n-1}\left(\epsilon \right)$$ E5

The initial conditions are $$P_{-1}\left(\epsilon \right)=Q_0\left(\epsilon \right)\equiv 0$$ , $$P_0\left(\epsilon \right)\equiv 1$$ , $$Q_0\left(\epsilon \right)\equiv 0$$ , and $$Q_1\left(\epsilon \right)\equiv b_0^{-1}$$ , function $$K_{\infty }\left(\epsilon \right)$$ corresponds to the LGF operator, with all elements of its J-matrix being equal to their limit values a and b:

$$K_{\infty }\left(\epsilon \right)=2b^{-2}\left[\epsilon -a-Z\left(\epsilon \right)\cdot \sqrt{\left(\epsilon -a+2b\right)\left(\epsilon -a-2b\right)}\right]$$ E6

$$Z\left(\epsilon \right)\equiv \Theta \left(a-2b-\epsilon \right)+i\Theta \left(\epsilon -a+2b\right)\cdot \Theta \left(a+2b-\epsilon \right)-\Theta \left(\epsilon -a-2b\right)$$ E7

The method of Jacobi matrices can treat as a regular singular perturbation a much larger number of perturbations of the phonon spectrum due to the presence of various crystal defects than the traditional methods [19, 20]. In addition, perturbations do not change the bandwidth of the quasi-continuous spectrum, and consequently, the asymptotic values of the elements of the J-matrix can be regarded as an asymptotically degenerated regular perturbation [23]. This type of perturbations covers virtually all perturbations of the phonon spectrum caused by local defects. The calculation of vibration characteristics of such systems is performed, using the method of J-matrices, with the same accuracy as for the initial ideal system.

In practice, it is usually possible to calculate the Jacobi matrix of the Hamiltonian of a finite rank. The expression

$$G\left(\epsilon ,{\overrightarrow{h}}_0\right)\approx \frac{Q_n\left(\epsilon \right)-b_{n-1}{Q}_{n-1}\left(\epsilon \right)K_{\infty }\left(\epsilon \right)}{P_n\left(\epsilon \right)-b_{n-1}{P}_{n-1}\left(\epsilon \right)K_{\infty }\left(\epsilon \right)}$$ E8

is called analytical approximation of LGF. All dependences in Figure 1 and curves 1 and 2 in Figure 2 were calculated by the formula (8), using the Jacobi matrix of the Hamiltonian (1) with rank $$n=600$$ . If we count the energy from the Fermi energy level, then all diagonal elements of Jacobi matrices are zero ( $$a_n=a=0{;}^{}b={\epsilon }_0/2$$ , where $${\epsilon }_0=3J$$ is the half-width of the quasi-continuous spectrum). A good accuracy of the approximations shown in figures is confirmed also by the fact that they show the nonanalyticity effects corresponding to the densities of states of systems with the dimension larger than unity (so-called van Hove singularities). In the vicinity of these singularities the expression (8) slowly converges to the true values of the really and imaginary parts of LGF.

Curves $$1^{\prime }$$ and $$2^{\prime }$$ in Figure 2 show local density of states and their corresponding real parts of the LGF calculated by formula (8) with n=1. As can be seen in the band of the quasi-continuous spectrum these relationships have very little in common with curves 1 and 2. Thus, the curve $$1^{\prime }$$ does not even hint at the V-shaped “Dirac” singularity at $$\epsilon =\epsilon \left(K\right)={\epsilon }_F$$ , and on the curve $$2^{\prime }$$ in the interval $$\left[-\epsilon \left(M\right),\epsilon \left(M\right)\right]$$ both non-monotonous parts and the logarithmic singularities at the edges of the band of quasi-continuous spectrum, characteristic for the 2D systems are absent. However, outside the band of quasi-continuous spectrum (also in the area of intersection of the real part of LGF with curve (3)) curves 2 and $$2^{\prime }$$ practically coincide, and if we put in the Lifshitz equation (2) instead of LGF its approximation (8) for n=1, the obtained solutions give the energies of LDL with quite high accuracy. Moreover, these solutions can be easily found analytically. The LGF approximation by formula (8) for n=1 was named the approximation of two moments in [23]. Indeed, it follows from the orthonormality of the polynomials defined in [16,17] that

$$\begin{array}{l}{\displaystyle \underset{-{\epsilon }_0}{\overset{-{\epsilon }_0}{\int }}{P}_1\left(\epsilon \right)\rho \left(\epsilon \right)}d\epsilon =0\Rightarrow a_0={\displaystyle \underset{-{\epsilon }_0}{\overset{-{\epsilon }_0}{\int }}\epsilon \rho \left(\epsilon \right)}d\epsilon =M_1,\\ {\displaystyle \underset{-{\epsilon }_0}{\overset{-{\epsilon }_0}{\int }}{P}_1^2\left(\epsilon \right)\rho \left(\epsilon \right)}d\epsilon ={\displaystyle \underset{-{\epsilon }_0}{\overset{-{\epsilon }_0}{\int }}\frac{\left({\epsilon }^2-a_0^2\right)}{b_0^2}\rho \left(\epsilon \right)}d\epsilon =1\Rightarrow b_0=\sqrt{M_2-M_1^2}\end{array}$$

Finding the characteristics of LDL is more convenient without using equation (2), looking for them as the poles of the LGF perturbed Hamiltonian $$\tilde{G}\left(\epsilon ,{\overrightarrow{h}}_0\right)=\left({\overrightarrow{h}}_0,{\left[\epsilon \widehat{I}-\widehat{H}-\widehat{\Lambda }\right]}^{-1}{\overrightarrow{h}}_0\right)$$ . In the approximation of two moments for the subspace generated by the excitation of an impurity atom, we get

$$\tilde{G}\left(\epsilon \right)=\frac{1}{\epsilon -a_0-b_0^2{K}_{\infty }\left(\epsilon \right)}$$ E9

For the case of an isolated substitutional impurity

$$a_0=\tilde{\epsilon };b_0=\sqrt{3}(1+\eta )J=\frac{1+\eta }{\sqrt{3}}{\epsilon }_0$$ E10

from where we get

$$\tilde{G}\left(\epsilon \right)=\frac{\left(1-\gamma \right)\epsilon -a_0+Z\left(\epsilon \right)\gamma \sqrt{\left|{\epsilon }^2-{\epsilon }_0^2\right|}}{R\left(\epsilon \right)}$$ E11

where $$R\left(\epsilon \right)=\left(1-2\gamma \right){\epsilon }^2-2\left(1-\gamma \right)a_0\epsilon +a_0^2+{\gamma }^2{\epsilon }_0^2$$ and

$$\gamma \equiv b_0^2/2b^2=2{\left(1+\eta \right)}^2/3$$ E12

Local discrete levels are poles (11), i.e. the roots of $$R\left(\epsilon \right)$$ are

$${\epsilon }_d^{(\pm )}=\frac{\left(\gamma -1\right)a_0\pm \gamma \sqrt{a_0^2+\left(2\gamma -1\right){\epsilon }_0^2}}{2\gamma -1}$$ E13

Residues at these poles $${\mu }_0^{(\pm )}=\underset{\epsilon ={\epsilon }_d^{(\pm )}}{r{}^{‵}es}\tilde{G}\left(\epsilon \right)$$ are called intensities of LDL and they determine the relative LDL “amplitude” on the impurity atom: $${\mu }_0^{(\pm )}=1-{\pi }^{-1}{\displaystyle \underset{-{\epsilon }_0}{\overset{{\epsilon }_0}{\int }}\mathrm{Im}\tilde{G}\left(\epsilon \right)}d\epsilon$$ . The condition that the intensity differs from zero defines the existence region of LDL. In this case

$$\begin{array}{c}{\mu }_0^{(+)}=\frac{\gamma a_0+\left(\gamma -1\right)\sqrt{a_0^2+\left(2\gamma -1\right){\epsilon }_0^2}}{\left(2\gamma -1\right)\sqrt{a_0^2+\left(2\gamma -1\right){\epsilon }_0^2}}\cdot \Theta \left(\frac{a_0}{{\epsilon }_0}-1+\gamma \right);\\ {\mu }_0^{(-)}=\frac{-\gamma a_0+\left(\gamma -1\right)\sqrt{a_0^2+\left(2\gamma -1\right){\epsilon }_0^2}}{\left(2\gamma -1\right)\sqrt{a_0^2+\left(2\gamma -1\right){\epsilon }_0^2}}\cdot \Theta \left(\gamma -1-\frac{a_0}{{\epsilon }_0}\right).\end{array}$$ E14

It was shown in [23] that

$$G_{mn}\left(\epsilon ,{\overrightarrow{h}}_0\right)=\left({\overrightarrow{h}}_m,{\left[\epsilon \widehat{I}-\widehat{H}\right]}^{-1}{\overrightarrow{h}}_n\right)=-P_m\left(\epsilon \right)Q_n\left(\epsilon \right)+P_m\left(\epsilon \right)P_m\left(\epsilon \right)G\left(\epsilon ,{\overrightarrow{h}}_0\right)\left(m<n\right).$$

This implies that the damping of LDL, i.e. the decay of its intensity with the increasing distance from the impurity atom (i.e. with the increase of n) follows the equation $${\mu }_n^{(\pm )}=P_n^2\left({\epsilon }_d^{(\pm )}\right)\cdot {\mu }_0^{(\pm )}$$ . Using the method of mathematical induction we can prove that

$$P_n\left({\epsilon }_d^{(\pm )}\right)=\sqrt{2\gamma }\cdot {\left[\pm \frac{\sqrt{a_0^2+\left(2\gamma -1\right){\epsilon }_0^2}\mp a_0}{{\epsilon }_0\left(2\gamma -1\right)}\right]}^n$$ E15

The intensities $${\mu }_n^{(\pm )}$$ decay with increasing n according to $${\mu }_{n0}^{(\pm )}=2\gamma \cdot q^n\cdot {\mu }_0^{(\pm )}$$ , that is, starting from n=1 they form an infinitely decreasing geometric progression whose denominator

$$q^{(\pm )}={\left[\frac{\sqrt{a_0^2+\left(2\gamma -1\right){\epsilon }_0^2}\mp a_0}{{\epsilon }_0\left(2\gamma -1\right)}\right]}^2$$ E16

Summing these progressions we see that $${\displaystyle \sum _{n=0}^{\infty }{\mu }_n^{(+)}=}{\displaystyle \sum _{n=0}^{\infty }{\mu }_n^{(-)}=}1$$ , that is the formation of each LDL is the formation of one quasi-particle outside the band of the quasi-continuous electron spectrum.

Formulas (10, 12-14, 16) give simple analytical expressions of the local conditions of the existence of discrete levels due to the presence of a substitutional impurity in graphene.

Regions of the LDL existence for $${\epsilon }_d^{(+)}$$ (in this case $$\tilde{\epsilon }{\epsilon }_0\left(1-\gamma \right)={\epsilon }_0\left[1-2{\left(1+\eta \right)}^2/3\right]$$ ) and for $${\epsilon }_d^{(-)}$$ (in this case $$\tilde{\epsilon }{\epsilon }_0\left(\gamma -1\right)={\epsilon }_0\left[2{\left(1+\eta \right)}^2/3-1\right]$$ ) lie above and under curves in Figure 3, as indicated by arrows.

It is seen that such levels exist in a very wide range of variables $$\tilde{\epsilon }$$ and $$\eta$$ . The absence of LDL is possible only in a narrow range of values $$\tilde{\epsilon }$$ at $$\eta \sqrt{3/2}-1$$ . If $$\eta \sqrt{3/2}-1$$ at least one local discrete level exists. In fact, the lines delineating the area of existence of LDL in graphene in the plane $$\left\{\tilde{\epsilon },\eta \right\}$$ must pass through the origin of coordinates, since $$\mathrm{Re}G\left(\epsilon \right)\to \infty$$ for $$\epsilon \to \pm {\epsilon }_0$$ . However, since this divergence is logarithmic, for any appreciable splitting of LDL from the boundary of quasi-continuous spectrum there is a certain threshold. Curves 2 and $$2^{\prime }$$ in Figure 2 merge at $$\left|\epsilon \right|>{\epsilon }_0$$ .

Figure 3.

Regions of the existence of discrete levels for substitutional impurities in graphene.

Figure 4 shows, for $$\tilde{\epsilon }\approx 0.525J$$ (for boron, from [6]), the dependences of energies, the LDL intensities at the boron impurity and the damping parameters of the value $$\eta$$ that characterizes the change in the overlap integral of the boron impurity atom (10). Solid lines show the characteristics of LDL, calculated using the approximation of the Green's function (9), i.e. according to the analytical formulas (13), (14) and (16). Open circles show the results of numerical calculations of dependences $${\epsilon }_d^{(\pm )}\left(\eta \right)$$ and $${\mu }_0^{(\pm )}\left(\eta \right)$$ , using the Green's function in the form of (8), calculated by the Jacobi matrices of the $$n\ge 100$$ rank.

Figure 4.

The basic characteristics of LDL in the presence of a boron substitutional impurity in graphene.

It is seen that at the threshold values $${\epsilon }_d^{(\pm )}=\pm {\epsilon }_0$$ of the LDL formation, the intensities of LDL equal zero and the parameters of damping are equal to unit. Further increase of $$\left|{\epsilon }_d^{(\pm )}\right|$$ is accompanied by the increase of $${\mu }_0^{(\pm )}$$ and by the strengthening of the containment level.

So, a good agreement between the results of numerical calculation of the LDL characteristics using Jacobi matrices of high rank, and their analytical description by the Green's function (9), which relates these characteristics to the parameters of the defect (13) (14), makes it relatively easy to extract the defect parameters from the known characteristics of LDL. Experimental measurement (e.g. by scanning tunneling microscopy) of values $${\epsilon }_d^{(\pm )}$$ should lead, using (10), (12) and (13), to the determination of the parameters $$\tilde{\epsilon }$$ and $$\eta$$ and this might represent a significant advance in creating nanomaterials with predetermined spectral characteristics. As can be seen from Figure 4, with increasing $$\eta$$ the intensity of LDL $${\mu }_0^{(+)}\to {\mu }_0^{(-)}\to 1/2$$ . That is, the impurity levels can not be completely localized on the impurities, but they also appear in the spectra of surrounding carbon atoms. This greatly increases the probability of experimental detection of such levels, even at low concentrations of impurities.

3. The electronic spectrum of bilayer graphene

Bilayer graphene is a carbon film consisting of two graphene monolayers, linked together by (as in bulk graphite) van der Waals forces. Since the distance between the layers (film thickness) $$h~3.5$$ Å the bilayer graphene can be considered not a nanofilm but a subnanofilm. The constants of the interatomic interaction of bilayer graphene were determined and its phonon density of states and partial contributions to this quantity from the atomic displacements along different crystallographic directions were calculated [7,24]. On the basis of the analyses of the mean-square amplitudes of atomic displacements calculated using data from the spectral densities, we have shown that the flat shape of a free bilayer graphene remains stable up to the temperatures much higher than the room temperature, which makes this compound promising for nanoelectronics. In this section we calculate and analyze the electronic spectrum of a defect-free bilayer graphene. Naturally of greatest interest is its behavior in the energy range close to $${\epsilon }_F$$ where there are characteristic Dirac points on the spectrum of graphene monolayer (whose plane shape is unstable).

The unit cell of graphene contains two physically equivalent atoms and therefore local Green's function and the local density of states (LDOS) of the atoms of different sublattices are identical. On the other hand, bilayer graphene unit cell contains four atoms, and atoms of different sublattices of a single graphene layer interact differently with the atoms of the other layer and their physical equivalence is disrupted (Figure 5a).

The electronic spectrum of bilayer graphene, as well as the electronic spectrum of graphene can be described in the strong coupling approximation. Corresponding Hamiltonian has form (1). For graphene and bigraphene we assume that the electron hopping within the layer is possible only between nearest neighbors $$J_{ij}=J\approx 2.8$$ eV (see for example [25]). Electron hopping between layers is also assumed to be possible only between nearest neighbors from different layers, that is, between those which lie at a distance h from each other. Denote the corresponding hopping integral $$J^{\prime }$$ . Note that only half of the bilayer graphene atoms have such neighbors (sublattice AI and AII, see Figure 1). In sublattices BI and BII no such neighbors exist, since nearest neighbors from different layers are at a distance $$\sqrt{h^2+a^2}$$ ( $$a\approx 1.415$$ Å is the distance between nearest neighbors in the layer plane). Since this distance is only by less than 10% greater than $$h$$ , we can neglect the interaction with the atoms of the sublattices BI and BII and this does not lead to qualitative changes in the behavior of the spectra near $${\epsilon }_F$$ (see, e.g. [6]). Then the dispersion relation of each of the four branches of the electronic spectrum of bilayer graphene can be written as

$$\begin{array}{l}{\epsilon }_{\mathrm{1,2}}\left(k\right)=\pm \sqrt{{\epsilon }_0^2\left(k\right)+\frac{{J^{\prime }}^2}{2}-J^{\prime }\sqrt{{\epsilon }_0^2\left(k\right)+\frac{{J^{\prime }}^2}{4}}};\\ {\epsilon }_{\mathrm{3,4}}\left(k\right)=\pm \sqrt{{\epsilon }_0^2\left(k\right)+\frac{{J^{\prime }}^2}{2}+J^{\prime }\sqrt{{\epsilon }_0^2\left(k\right)+\frac{{J^{\prime }}^2}{4}}},\end{array}$$ E17

where $${\epsilon }_0\left(k\right)$$ is the electronic spectrum of graphene, calculated in the strong coupling approximation:

$${\epsilon }_0\left(k\right)=\pm J\sqrt{1+4\cos \left(k\cdot \frac{a_1-a_2}{2}\right)\left[\cos \left(k\cdot \frac{a_1+a_2}{2}\right)+\cos \left(k\cdot \frac{a_1-a_2}{2}\right)\right]}$$ E18

where $$a_1=\left(\frac{a}{2},\frac{a\sqrt{3}}{2}\mathrm{,0}\right)$$ and $$a_2=\left(\frac{a}{2},-\frac{a\sqrt{3}}{2}\mathrm{,0}\right)$$ are two-dimensional Bravais lattice vectors (see Figure 5b).

Figure 5.

The structure of bilayer graphene (a); Bravais lattice and the first Brillouin zone (b).

Dispersion curves along highly symmetric directions ΓK, ΓM and KM for J’ = 0.1J are shown in Figure 6. Region near the K-point $${\rm K}=\left(\mathrm{0,}\pm 4\pi /3a\sqrt{3}\mathrm{,0}\right)$$ and $${\rm K}=\left(\pm 2\pi /3a,\pm 2\pi /3a\sqrt{3}\mathrm{,0}\right)$$ is shown in the inset. The same inset shows also the dispersion curves for graphene (18). We clearly see the quasi-relativistic nature of the electronic spectrum of graphene as well as ordinary quadratic dispersion curves $${\epsilon }_{\mathrm{1,2}}\left(k\right)$$ for bigraphene. Spectral branches $${\epsilon }_{\mathrm{3,4}}\left(k\right)\notin \left(-J^{\prime },J^{\prime }\right)$$ are determined in the $$\left(-J^{\prime },J^{\prime }\right)$$ interval. Indeed, if k takes value along ΓK, then $${\epsilon }_0\left(k\right)=\pm \left(1+2\cos \frac{ak\sqrt{3}}{2}\right)$$ and putting $$k={\rm K}+\kappa$$ ( $$\kappa 1$$ ) we find $${\epsilon }_0\left({\rm K}+\kappa \right)\approx \mp \frac{3a\kappa }{2}$$ , i.e. a linear (relativistic) dispersion relation. For electronic modes of bilayer graphene $${\epsilon }_{\mathrm{1,2}}\left(k\right)$$ near the $${\rm K}$$ point, we can write

$$\begin{array}{l}{\epsilon }_{\mathrm{1,2}}^2\left({\rm K}+\kappa \right)\approx \frac{J_1^2}{2}+{\epsilon }_0^2\left({\rm K}+\kappa \right)-\sqrt{{\left(\frac{J_1^2}{2}+{\epsilon }_0^2\left({\rm K}+\kappa \right)\right)}^2-{\epsilon }_0^4\left({\rm K}+\kappa \right)}\approx \frac{J_1^2}{2}+{\epsilon }_0^2\left({\rm K}+\kappa \right)-\\ -\left[\frac{J_1^2}{2}+{\epsilon }_0^2\left({\rm K}+\kappa \right)\right]\cdot \left\{1-\frac{{\epsilon }_0^4\left({\rm K}+\kappa \right)}{2{\left[J_1^2/2+{\epsilon }_0^2\left({\rm K}+\kappa \right)\right]}^2}\right\}=\frac{{\epsilon }_0^4\left({\rm K}+\kappa \right)}{J_1^2+2{\epsilon }_0^2\left({\rm K}+\kappa \right)}\approx \frac{81J^4{a}^4{\kappa }^4}{16J_1^2},\end{array}$$ E19

that is, we get an ordinary quadratic dispersion law

$${\epsilon }_{\mathrm{1,2}}^{}\left({\rm K}+\kappa \right)\approx \pm \frac{9J^2{a}^2{\kappa }^2}{4J_1^{}}$$ E20

Figure 6.

Dispersion curves along high-symmetry directions of bilayer graphene.

The electron effective mass in considered branches, determined from the relation $$\epsilon ={\left(\hslash \kappa \right)}^2/2m^{*}$$ , is equal to

$$m^{*}=\frac{2{\hslash }^2{J}^{\prime }}{9J^2{a}^2}$$ E21

In the case considered above $$J^{\prime }=0.1J$$ the effective mass is $$m^{*}\approx 2.75\cdot {10}^{-32}$$ kg, and if $$J^{\prime }\to J$$ it tends to a value close to the free electron mass. Since $$m^{*}~J^{\prime }$$ , then by changing the interlayer hopping integral we can change the effective mass of charge carriers.

It should be noted that for $$J^{\prime }\to 0$$ , that is for the transition from bilayer graphene to two non-interacting graphene monolayers the effective mass, $$m^{*}\to 0$$ , and formulas (19 - 21) cannot correctly describe this transition, since they were obtained under the assumption $${\epsilon }_0\left(k\right)<<J^{\prime }$$ .

Electron density of states for values of energy near ε_F is determined by branches of $${\epsilon }_1$$ and $${\epsilon }_2$$ only, and it follows from (3) that $$g_1\left(\epsilon \right)=g_2\left(-\epsilon \right)$$ . Then

$$g\left(\epsilon \right)=\frac{{\Sigma }_0}{{\left(2\pi \right)}^2}{\displaystyle \underset{\epsilon \left(\kappa \right)=\epsilon }{\oint }\frac{d{l}_{\mathrm{1,2}}}{\left|\partial {\epsilon }_{\mathrm{1,2}}/\partial \kappa \right|}}$$ E22

where $${\Sigma }_0=3a^2\sqrt{3}/2$$ is the area of the two-dimensional Bravais cell, and integration is done along a closed isoenergetic line $${\epsilon }_{\mathrm{1,2}}\left(k\right)=\epsilon$$ . At $$\epsilon =0$$ (Fermi level) the line contracts into a point and near ε_F the contour of integration is a circle. Taking into account (20) we may write

$$\begin{array}{c}{g}_{\mathrm{1,2}}\left({\epsilon }_F\right)\approx \frac{S_0}{{\left(2\pi \right)}^2}{\displaystyle \underset{\epsilon \left(\delta \right)=\epsilon }{\oint }\frac{\delta d\phi }{\left|\partial {\epsilon }_{\mathrm{1,2}}/\partial \kappa \right|}}=\frac{3a^2\sqrt{3}}{2{\left(2\pi \right)}^2}\cdot 2\pi \frac{\delta }{\frac{18J^2{a}^2\kappa }{4J^{\prime }}}=\frac{J^{\prime }\sqrt{3}}{6\pi J^2}=const\\ g_{\mathrm{3,4}}\left({\epsilon }_F\right)=0;\Rightarrow g\left({\epsilon }_F\right)\approx \frac{J^{\prime }}{J^2\sqrt{3}}.\end{array}$$ E23

This means that the electron density of states (DOS) is constant and different from zero. As follows from (20), DOS is an analytical function and has minimum at $$\epsilon =0$$ and near $${\epsilon }_F$$ the function $$g\left(\epsilon \right)~{\epsilon }^2$$ .

Total electron DOS can be represented as a mean-arithmetic function of the two LDOS corresponding to atoms of sublattices A and B: $${\rho }_{AI}\left(\epsilon \right)={\rho }_{AII}\left(\epsilon \right)\equiv {\rho }_A\left(\epsilon \right);$$ $${\rho }_{BI}\left(\epsilon \right)={\rho }_{BII}\left(\epsilon \right)\equiv {\rho }_B\left(\epsilon \right)$$ and $$g\left(\epsilon \right)=\left[{\rho }_A\left(\epsilon \right)+{\rho }_B\left(\epsilon \right)\right]/2.$$ For each sublattice with perfect structure the LDOS may be written as

$${\rho }_s\left({\epsilon }_F\right)\approx \frac{S_0}{{\left(2\pi \right)}^2}{\displaystyle \sum _{\alpha =1}^4{\displaystyle \underset{\epsilon \left(\delta \right)=\epsilon }{\oint }\frac{\delta d\phi {\left|{\psi }_s\left(\alpha ,\kappa \right)\right|}^2}{\left|\partial {\epsilon }_{\mathrm{1,2}}/\partial \kappa \right|}}}$$ E24

where index s is the designation of sublattice, index α is the designation of branch, and $${\psi }_s\left(\alpha ,\kappa \right)$$ are the eigenwave functions corresponding to atoms from sublattices. LDOS are calculated by the method of Jacobi matrix [16, 17] and are shown in Figure 7.

In this figure we clearly see two-dimensional van-Hove peculiarities at energy values corresponding to points $$\Gamma$$ and $${\rm M}$$ of the first Brillouin zone ( see Figure 5b). These eigenvalues are given in the top inset in Figure 7. The bottom inset shows LDOS, in enlarged scale, near $${\epsilon }_F$$ . It is seen that near the Fermi level the local density of states as well the total density of states are analytical and their dependences on energy are essentially nonlinear (for comparison the DOS of graphene is also shown in the same inset). Besides, $${\rho }_A\left(\epsilon \right)$$ differs from both $${\rho }_B\left(\epsilon \right)$$ and the total density of states and it approaches to zero for $$\epsilon \to 0$$ .

Figure 7.

LDOS for atoms of different sublattices of bigraphene.

Indeed, putting zero eigenvalue in the equation of eigenfunctions of Hamiltonian (1) we get the values $${\psi }_{AI}~{\psi }_{AII}=O\left({\kappa }^2\right)$$ , $${\psi }_{BI}~{\psi }_{BII}~1$$ . Therefore, near the Fermi level $${\rho }_A\left(\epsilon \right)~{\epsilon }^2{\rho }_B\left(\epsilon \right)$$ .

Peculiarities at both ends, i.e. at $$\epsilon \approx \pm 0.1J=\pm J^{\prime }$$ , originate from the contribution of modes ε_3,4 which are not represented on interval $$\epsilon \in \left[-J^{\prime },J^{\prime }\right]$$ . Peculiarities at $$\epsilon \approx \pm 0.05J$$ are due to the fact that beginning from these energies the anisotropy of isofrequency lines becomes essential.

So we can conclude that, in contrary to graphene, the bigraphene has an ordinary non-relativistic form of the electronic dispersion law. The effective mass of electron in the bigraphene strongly depends on the value of integral describing the hopping between two layers, and this value may be changed by external conditions (for example by pressure). Near the Fermi level the LDOS of atoms of different sublattices qualitively differ from each other. If the LDOS for the atoms of sublattice A at the Fermi level equals to zero and it slowly increases near this level (“a quasi-gap” appears), then the LDOS of sublattice B for the same energy values differs from zero and increases very quickly.

4. The influence of defects on the electron spectrum of bigraphene

Some peculiarities in the behavior of the bigraphene electron spectrum near $$\epsilon =\epsilon \left({\rm K}\right)$$ , that is near the Fermi level, indicate the possibility of a strong influence of various defects [26-28].

Let us first consider point defects [26]. Figure 8 shows the real parts of the local Green's functions $$\mathrm{Re}{G}_A\left(\epsilon \right)$$ and $$\mathrm{Re}{G}_B\left(\epsilon \right)$$ of atoms of the sublattices A and B.

Figure 8.

Real parts of the local Green’s functions for the atoms of the two sublattices of bigraphene. The dashed lines are $$S\left(\epsilon \right)$$ functions correspond to substitutional impurities of Boron (B-curves) and Nitrogen (N-curves).

This figure also shows the dependences S(ε) appearing in the Lifshitz equation (2) and corresponding correspond to the presence of substitutional impurities of boron and nitrogen in bigraphene.

For an isolated substitutional impurity that differs from the atom of the basic lattice by the values of energies at the impurity site and also by the overlap integral, function $$S\left(\epsilon \right)$$ has for each sublattice form (3). Also, as in the case of graphene, for the nitrogen impurity in the interval $$\left[{\epsilon }_2\left({\rm M}\right),{\epsilon }_1\left({\rm M}\right)\right]$$ equation (2) has a solution and this impurity forms quasilocalized states in this interval. For the boron impurity, equation (2) has two solutions outside the band of quasi-continuous spectrum $$\left[{\epsilon }_4\left({\rm M}\right),{\epsilon }_3\left({\rm M}\right)\right]$$ , corresponding to local discrete levels.

LDOS of isolated impurity atom of nitrogen in sublattice A or B are shown in Figure 9 (we remind that atoms AI and AII as well as BI and BII are physically equivalent).

As the inset shows, the nitrogen impurity does not forms a quasilocal maximum on LDOS and substantially changes it near the Fermi level. Figure 10 shows a LDOS of the boron impurity in sublattices A and B. As in the case of boron impurity in graphene (see Figure 1), the area under this curve is less than unity. Outside the band of the quasi-continuous spectrum local discrete energy levels are formed.

Figure 9.

LDOS of nitrogen impurity in the sublattices of bigraphene (red curve is for A sublattice, purple is for B sublattice). For comparison the figure shows the DOS of graphene (black dashed line) and LDOS (thin solid gray and dashed gray, respectively).

Figure 10.

LDOS of the boron impurity in different sublattices of bigraphene (designations are the same as in Figure 9).

The energies of these levels for the considered ratio between the overlap integrals J and J’ are slightly different. As can be seen from Figure 8, the energies of the local discrete levels for the substituted boron atom in sublattice A or sublattice B can be calculated using the two-moment approximation (13). Near the Fermi level the LDOS of the boron atom, substituting an atom of sublattice B, is considerably lower than the LDOS of the carbon atom of this sublattice. That is the boron impurity lowers the conductivity of bigraphene.

Because in bigraphene atoms of sublattices A and B are physically inequivalent, the influence on their electron spectra by various defects is different. In the first part of this section we have described the influence of substitutional impurities on the electron density of states. However, the influence of vacancies in sublattices on the electron DOS is even more profound. Figures 11 and 12 show the LDOS of neighbors of vacancies in sublattices A (Figure 11) and B (Figure 12).

Neighboring atoms are in the same layer as a vacancy, either in sublattice B (top) or in sublattice A (bottom). Insets show the arrangement of atoms. Atoms are shown in the same color as the corresponding LDOS.

Figure 11.

LDOS of bigraphene atoms which are neighbors of a vacancy in the sublattice A. Neighboring atoms are in the same layer as a vacancy, either in sublattice B (top) or in sublattice A (bottom). Insets show the arrangement of atoms. Atoms are shown in the same color as the corresponding LDOS.

Figure 12.

LDOS of bigrafene atoms which are neighbors of a vacancy in the sublattice B.

The reason for the specific evolution of the LDOS with increasing distance from the vacancy has been explained in [22]. In bigraphene there is an analogous situation. Electronic spectra of vacancy neighbors belonging to other sublattice also have sharp resonance peaks at ε = ε_F (see upper parts Figure 11 and 12). So, if the vacancy is in the sublattice A, then the maximum of the LDOS of the sublattice B atoms is sharp, and the maximum on the spectrum the sublattice A atoms is blurred. Therefore we can conclude that the vacancy in bigraphene should have a more pronounced effect.

5. Conclusion

Unique properties of both graphene and bigraphene are caused, above all, by an unusual symmetry leading to the absence of a gap between the valence and conduction bands. The quantum states of quasiparticle in these bands are described by the same wave function, i.e. quasiparticle in graphene and bigraphene have the so-called chiral symmetry. In graphene this leads to the fact that dispersion relation of electron spectrum is linear and is described by the Dirac equation, characteristic for massless ultrarelativistic quasiparticle in quantum electrodynamics. On other hand, in bigraphene the presence of chiral symmetry leads at low energies to an ordinary parabolic dispersion relation, i.e. quiasiparticle of a new type appear – massive chiral fermions having no analogy in quantum electrodynamics. At symmetry breaking a gap between the valence and conduction band appears, allowing to tune the conducting properties of these materials.

In this chapter, using the method of Jacobi matrices, we analyzed how different impurities affect the energy gap width and the local density of states. The method of Jacobi matrices enables to investigate heterogeneous systems and to calculate the densities of states for each atom in different sublattices. Such analysis is necessary for correct comparison with experimental results. As defects we considered the vacancies and the impurity of nitrogen or boron. The presence of nitrogen leads to the formation of sharp resonance peaks (quasilocal states) inside the continuous spectrum; on the other hand, the boron impurity leads to states outside the continuous spectrum (local states). Both quasilocal and local states can be investigated experimentally by, for example, a scanning tunneling microscope. In the presence of vacancies we have analyzed how the density of states in each of the sublattices A or B changes. Different situations were analyzed. For example, if the vacancy in this sublattice vanishes at Fermi level, whereas in sublattice B the density has maximum at this point. We also investigated the conditions for opening the gap and changing its width. Main attention was paid to the analysis of electronic properties of considered systems. Moreover, computational method we used has also been successfully applied to the analysis of vibrational states. For example, increasing the overlap integral between the boron impurity and carbon atom leads to the strengthening of force interaction contacts between them. In addition, boron atom is 16% lighter than carbon atom, i.e. all necessary conditions are sent for the appearance of local vibrations in the phonon spectrum of both graphene and bigraphene with boron impurity. We hope that predicted peculiarities in electronic and vibrational spectra of perfect graphene and bigraphene as well as of graphene and bigraphene with defects will be detected experimentally.

Acknowledgments

This work is the result of the project implementation: Research and Education at UPJŠ – Heading towards Excellent European Universities, ITMS project code: 26110230056, supported by the Operational Program Education funded by the European Social Fund (ESF). This work was also supported by the grant of the Ukrainian Academy of Sciences under the contract No 4/10-H and by the grant of the Scientific Grant Agency of the Ministry of Education, Science, Research and Sport of the Slovak Republic and Slovak Academy of Sciences under No. 1/0159/09.