Graphene coherent states

TL;DR

This paper constructs and analyzes coherent states for a Dirac electron in graphene under a magnetic field, exploring their physical properties and quantum characteristics.

Contribution

It introduces a method to derive coherent states for graphene electrons in magnetic fields, including analysis of their physical and quantum properties.

Findings

Coherent states are explicitly constructed as eigenstates of the annihilation operator.

Heisenberg uncertainty product and probability density are analyzed.

Mean energy values of the states are calculated.

Abstract

In this paper, we will construct the coherent states for a Dirac electron in graphene placed in a constant homogeneous magnetic field which is orthogonal to the graphene surface. First of all, we will identify the appropriate annihilation and creation operators. Then, we will derive the coherent states as eigenstates of the annihilation operator, with complex eigenvalues. Several physical quantities, as the Heisenberg uncertainty product, probability density and mean energy value, will be as well explored.

Equations

$$-i\hslash v_F\vec{\sigma }\cdot \nabla \Psi (\mathbf{r})=E\Psi (\mathbf{r}).$$

$$v_F\vec{\sigma }\cdot \left(\vec{p}+\frac{e\vec{A}}{c}\right)\Psi (x,y)=E\Psi (x,y),$$

$$\Psi(x,y)=e^{iky}\left(\begin{array}[]{c}\psi^{+}(x)\\ i\psi^{-}(x)\\ \end{array}\right),$$

$$\left(\pm \frac{d}{dx}+\frac{e}{c\hslash }A(x)+k\right){\psi }^{\mp }(x)=\frac{E}{\hslash v_F}{\psi }^{\pm }(x),$$

$$H^{\pm }=-\frac{d^2}{d{x}^2}+V^{\pm },V^{\pm }={\left(\frac{eA(x)}{c\hslash }+k\right)}^2\pm \frac{e}{c\hslash }\frac{dA(x)}{dx},\mathcal{E}=\frac{E^2}{{\hslash }^2{v}_F^2}.$$

$$B_0=\frac{c\hslash }{2e}\omega \to \omega =\frac{2e{B}_0}{c\hslash },$$

$$V^{\pm }=\frac{{\omega }^2}{4}{\left(x+\frac{2k}{\omega }\right)}^2\pm \frac{1}{2}\omega .$$

$${\mathcal{E}}_0^{-}=0,{\mathcal{E}}_n^{-}={\mathcal{E}}_{n-1}^{+}=n\omega ,n=1,2,\dots$$

$${\psi }_n^{\pm }(x)=\sqrt{\frac{1}{2^{n}n!}{\left(\frac{\omega }{2\pi }\right)}^{1/2}}{H}_n\left[\sqrt{\frac{\omega }{2}}\left(x+\frac{2k}{\omega }\right)\right]e^{-\frac{\omega }{4}{\left(x+\frac{2k}{\omega }\right)}^2},$$

$$E_n^{\pm }=\pm \hslash v_F\sqrt{n\omega },n=0,1,\dots ,$$

$$\Psi_{0}(x,y)=e^{iky}\left(\begin{array}[]{c}0\\ i\,\psi^{-}_{0}(x)\\ \end{array}\right),\quad\Psi_{n}(x,y)=\frac{e^{iky}}{\sqrt{2}}\left(\begin{array}[]{c}\psi^{+}_{n-1}(x)\\ i\,\psi^{-}_{n}(x)\\ \end{array}\right),\quad n=1,2,\ldots$$

$$\hat{A}^{-}=\left(\begin{array}[]{cc}f_{1}(\hat{N})\hat{\vartheta}^{-}&0\\ 0&f(\hat{N}+\hat{1})\hat{\vartheta}^{-}\\ \end{array}\right),$$

$${\hat{\vartheta }}^{-}=\frac{1}{\sqrt{2}}\left(z+{\partial }_z\right),{\hat{\vartheta }}^{+}=\frac{1}{\sqrt{2}}\left(z-{\partial }_z\right),\hat{N}={\hat{\vartheta }}^{+}{\hat{\vartheta }}^{-},$$

$${\hat{A}}^{-}{\Psi }_n$$

$${\hat{A}}^{-}{\Psi }_n=c_n{\Psi }_{n-1}$$

$$\sqrt{n-1}{f}_1(n-2)=\sqrt{n}f(n),n=2,3,\dots$$

$$c_n$$

$$\displaystyle\hat{A}^{-}=\left(\begin{array}[]{cc}\frac{\sqrt{\hat{N}+\hat{2}}}{\sqrt{\hat{N}+\hat{1}}}f(\hat{N}+\hat{2})\hat{\vartheta}^{-}&0\\ 0&f(\hat{N}+\hat{1})\hat{\vartheta}^{-}\end{array}\right).$$

$${\hat{A}}^{-}{\Psi }_{\alpha }=\alpha {\Psi }_{\alpha },\alpha \in \mathbb{C}.$$

$$\Psi_{\alpha}(x,y)=\sum_{n=0}^{\infty}a_{n}\Psi_{n}(x,y)=\left[a_{0}\left(\begin{array}[]{c}0\\ i\psi^{-}_{0}(x)\end{array}\right)+\sum_{n=1}^{\infty}\frac{a_{n}}{\sqrt{2}}\left(\begin{array}[]{c}\psi^{+}_{n-1}(x)\\ i\psi^{-}_{n}(x)\\ \end{array}\right)\right]e^{iky}.$$

$$a_{n+1}=\frac{{\alpha }^n{a}_1}{\sqrt{(n+1)!}f(n+1)\cdots f(2)}=\frac{\sqrt{2}{\alpha }^{n+1}{a}_0}{\sqrt{(n+1)!}[f(n+1)]!},$$

$$[f(n)]!\equiv \{\begin{array}{ll}1& \text{for}n=0\\ f(1)\cdots f(n)& \text{for}n=1,2,\dots \end{array}$$

$${\Psi }_{\alpha }(x,y)={\left[1+\sum _{n=1}^{\infty }\frac{2|\alpha {|}^{2n}}{n!([f(n)]!{)}^2}\right]}^{-1/2}\left[{\Psi }_0(x,y)+\sum _{n=1}^{\infty }\frac{\sqrt{2}{\alpha }^n}{\sqrt{n!}[f(n)]!}{\Psi }_n(x,y)\right].$$

$$a_{n+1}\sqrt{n+1}f(n+1)=\alpha a_n,n=1,2,\dots$$

$$a_{n+1}=\frac{{\alpha }^n}{\sqrt{(n+1)!}[g(n)]!}{a}_1.$$

$${\Psi }_{\alpha }(x,y)={\left[\sum _{n=0}^{\infty }\frac{|\alpha {|}^{2n}}{(n+1)!([g(n)]!{)}^2}\right]}^{-1/2}\sum _{n=0}^{\infty }\frac{{\alpha }^n}{\sqrt{(n+1)!}[g(n)]!}{\Psi }_{n+1}(x,y).$$

$${\Psi }_{\alpha }(x,y)={\left[\sum _{n=0}^{\infty }\frac{|\alpha {|}^{2n}}{(n+2)!([h(n)]!{)}^2}\right]}^{-1/2}\sum _{n=0}^{\infty }\frac{{\alpha }^n}{\sqrt{(n+2)!}[h(n)]!}{\Psi }_{n+2}(x,y),$$

$$\hat{z}=\frac{1}{\sqrt{2}}({\hat{\vartheta }}^{+}+{\hat{\vartheta }}^{-}),\hat{p}=\frac{i}{\sqrt{2}}({\hat{\vartheta }}^{+}-{\hat{\vartheta }}^{-}).$$

$${\sigma }_z^2{\sigma }_p^2\ge \frac{1}{4},$$

$${\Psi }_{\alpha }(x,y)=\frac{1}{\sqrt{2e^{r^2}-1}}\left[{\Psi }_0(x,y)+\sum _{n=1}^{\infty }\frac{\sqrt{2}{\alpha }^n}{\sqrt{n!}}{\Psi }_n(x,y)\right],$$

Full text

1 Introduction

Graphene is a single layer of carbon atoms arranged in a hexagonal honeycomb lattice, which is the basic structural element of other allotropes including graphite, charcoal, carbon nanotubes and fullerenes.

Graphene is a zero-gap semiconductor, because its conduction and valence bands meet at the Dirac points which are six locations in momentum space, on the edge of the Brillouin zone, divided into two non-equivalent sets of three points, typically labeled as $K$ and $K^{\prime}$. By contrast, for traditional semiconductors the point of primary interest is denoted as $\Gamma$, where the momentum is zero [1, 2, 3, 4, 5, 6] (see Figure 1).

Thus, even neglecting their spin, at low energies the electrons can be described by an equation that is formally equivalent to the massless Dirac equation:

[TABLE]

Here $v_{F}\sim 10^{6}$ m/s (.003 $c$) is the Fermi velocity in graphene, which replaces the velocity of light in Dirac theory, $\vec{\sigma}$ is the vector of Pauli matrices, $\Psi(\mathbf{r})$ is the two-component wave function of the electrons and $E$ is its energy [7].

As a consequence, the electrons and holes are called Dirac fermions; their appearance was predicted in the silicene, germanene or dichalcogenides [8, 9, 10, 11, 12, 13, 14] and emerge naturally from the tight-binding model for a generic hexagonal lattice in the low-energy regime [15], as shown in Figure 1. This pseudo-relativistic description is restricted to the chiral limit, i.e., to vanishing rest mass leading to interesting additional features: when Dirac fermions are compared with ordinary electrons placed in magnetic fields, their behavior leads to new physical phenomena such as the anomalous integer quantum Hall effect, the Zitterbewegung and the Klein paradox [8, 16].

It is important to stress that graphene belongs to a class of systems in condensed matter for which the low-energy quasi-particles behave like massless or massive Dirac fermions. These systems are known as Dirac materials in the literature [17].

As we shall see below, under particular physical conditions, a problem similar to that considered in [18] arises naturally. Due to this, it seems obvious the need to build up the coherent states for the graphene, and then to analyze their properties.

In order to do that, this paper is organized as follows. In section 2 the Dirac-Weyl equation will be introduced and the physical problem to be considered will be briefly discussed. In section 3 the annihilation operator associated to the system will be defined, and the corresponding coherent states will be constructed as eigenstates of that operator. We will analyze as well several physical quantities for these states. Our conclusions will be presented in section 4.

2 Dirac-Weyl equation

Let us suppose now that the graphene is placed in a static magnetic field which is orthogonal to the material surface (the $x-y$ plane) [7, 19, 20]. The interaction of a Dirac electron with such a field close to a Dirac point $K$ in the Brillouin zone is described by the Dirac-Weyl equation, which is obtained by replacing in Eq. (1) the momentum operator $\vec{p}=-i\hbar\,\nabla$ by $\vec{p}+e\vec{A}/c$, leading to:

[TABLE]

where $-e$ is the charge of the electron. Landau gauge is conveniently chosen, with the vector potential given by $\vec{A}=A(x)\hat{e}_{y}$ and $\vec{B}=\nabla\times\vec{A}=B(x)\hat{e}_{z}$, $B(x)=A^{\prime}(x)$. Taking into account the translational invariance along the $y$ direction, the two-component spinor $\Psi(x,y)$ is expressed as:

[TABLE]

being $k$ the wave number in the $y$ direction and $\psi^{\pm}(x)$ describing the electron amplitude on two adjacent sites in the unit cell of graphene. Thus, the Dirac-Weyl equation (2) yields two coupled first-order linear differential equations

[TABLE]

which can be easily decoupled into two Schrödinger equations $H^{\pm}\psi^{\pm}(x)=\mathcal{E}\psi^{\pm}(x)$, where

[TABLE]

For a constant magnetic field, orthogonal to the graphene surface and pointing in the positive $z$ direction ($\vec{B}=B_{0}\hat{e}_{z}$ with $B_{0}>0$), the vector potential is selected as $\vec{A}=B_{0}x\hat{e}_{y}$. Introducing now the constant $\omega$ as:

[TABLE]

whose dimensions are $($lenght$)^{-2}$, the potentials in (5) become two shifted oscillators of the form

[TABLE]

Thus, the eigenvalues $\mathcal{E}^{\pm}_{n}$ for the Hamiltonians $H^{\pm}$ are related as follows:

[TABLE]

and the associated eigenfunctions are those of the standard harmonic oscillator:

[TABLE]

where $H_{n}$ denotes the Hermite polynomial of degree $n\in\mathbb{N}$.

We conclude that the complete solution of the corresponding Dirac-Weyl equation in a constant magnetic field consists of the eigenvalues

[TABLE]

where the plus (minus) sign refers to the enegy electrones (holes), and the normalized eigenvectors

[TABLE]

Note that the solutions of the Dirac-Weyl equation have an internal degree of freedom that mimics the spin, called pseudospin. It admits the following interpretation: each component is the projection of the particle wavefunction onto the sublattice A (spin up) or B (spin down).

3 Annihilation operator

Since the eigenstates of the previous Dirac-Weyl equation are expressed in terms of the eigenfunctions of the standard harmonic oscillator, it seems natural to look for an annihilation operator for the Hamiltonian in Eq. (2). In fact, let $\hat{A}^{-}$ be the operator defined by

[TABLE]

where $\hat{\vartheta}^{\pm}$, $\hat{N}$ are given by

[TABLE]

with $z=\sqrt{\omega/2}\,(x+2k/\omega)$, and $f,\ f_{1}$ are two real adjustable functions which will be used to guarantee that $\hat{A}^{-}\Psi_{n}=c_{n}\Psi_{n-1}$. Then:

[TABLE]

In order to ensure that

[TABLE]

it must happen that

[TABLE]

In such a case it is obtained that:

[TABLE]

and the explicit expression for the annihilation operator $\hat{A}^{-}$ turns out to be:

[TABLE]

Equation (15) indicates that the explicit form of the function $f(n)$ is required to determine the complete action of the annihilation operator $\hat{A}^{-}$ onto the eigenstates $\Psi_{n}$. It is also quite important for the properties of the graphene coherent states, as it will be immediately seen.

3.1 Coherent states as eigenvectors of $\hat{A}^{-}$

Let us define the coherent states $\Psi_{\alpha}$ as eigenstates of the annihilation operator $\hat{A}^{-}$ with complex eigenvalue $\alpha$:

[TABLE]

Expressing $\Psi_{\alpha}$ as a linear combination of the states $\Psi_{n}$, we have

[TABLE]

Using equation (19) we find a recurrence relation for the coefficients $a_{n}$ which depends on the value taken by $f(1)$. We can identify two different cases.

3.1.1 Case with $f(1)\neq 0$

First of all suppose that $f(n)\neq 0\ \forall\ n=1,2,\dots$ Thus we get $a_{1}={\sqrt{2}\alpha a_{0}}/{f(1)}$ and

[TABLE]

where

[TABLE]

The free constant $a_{0}$ is used to normalize $\Psi_{\alpha}$; we obtain:

[TABLE]

3.1.2 Case with $f(1)=0$

If $f(1)=0$ we obtain that $a_{0}=0$ and the following recurrence relationship:

[TABLE]

Now, depending on the value of $f(2)$, two possibilities appear once again.

A. Case with $f(2)\neq 0$.

If we suppose that $f(n)\neq 0$ $\forall$ $n=2,3,\ldots$ and define $g(n)\equiv f(n+1)$, equation (23) leads to:

[TABLE]

Substituting this expression in equation (20) and then normalizing we obtain:

[TABLE]

B. Case with $f(2)=0$.

On the other hand, if $f(2)=0$ and $f(n)\neq 0$ $\forall\ n=3,4,\dots$ the normalized coherent states turn out to be now:

[TABLE]

where $h(n)\equiv f(n+2)$.

Let us notice that the graphene coherent states of Eqs. (22, 25, 26) look similar to the so-called vector coherent states. For information concerning the last states, the reader can seek e.g. [21, 22, 23] and references therein.

3.2 Mean values and Heisenberg uncertainty relation

Let the dimensionless position and momentum operators be given by:

[TABLE]

In units of $\hbar$, the Heisenberg uncertainty relation is expressed by

[TABLE]

where $\sigma_{S}^{2}\equiv\langle\hat{S}^{2}\rangle-\langle\hat{S}\rangle^{2}$ for an arbitrary observable $\hat{S}$.

We will calculate next these quantities for some examples of the coherent states, which will stress the important role played by the function $f(n)$ in our treatment.

3.2.1 The case with $f(1)\neq 0$.

Let us consider first the particular choice $f(\hat{N})=\hat{1}$. Thus, expression (22) leads to:

[TABLE]

where $r=|\alpha|$.

Using these coherent states, the mean values for the operators $\hat{z},\hat{p}$ of Eq. (27) as well as their squares become:

[TABLE]

Through them it is straightforward to calculate $(\sigma_{z})^{2}_{\alpha}(\sigma_{p})^{2}_{\alpha}$ (see Figure 2). Note that in the limit $\alpha\rightarrow 0$ we have $(\sigma_{z})^{2}_{\alpha}(\sigma_{p})^{2}_{\alpha}\rightarrow 1/4$.

3.2.2 The case with $f(1)=0$.

As we saw in section 3.1.2, when $f(1)=0$ two options appear, which depend on the value taken by $f(2)$.

A. The case with $f(2)\neq 0$.

Let us choose now $f(\hat{N}+\hat{1})=g(\hat{N})=\frac{\sqrt{\hat{N}}}{\sqrt{\hat{N}+\hat{1}}}$, so that $f(n)\neq 0$ $\forall$ $n=2,3,\ldots$. From Eq. (25), the explicit form for the normalized coherent states becomes:

[TABLE]

The mean values for the operators of Eq. (27) and their squares become:

[TABLE]

In the limit $\alpha\rightarrow 0$ it turns out that $(\sigma_{z})^{2}_{\alpha}(\sigma_{p})^{2}_{\alpha}\rightarrow 1$ (see Figure 3).

B. The case with $f(2)=0$.

Let us consider finally that $f(\hat{N}+\hat{2})=h(\hat{N})=\frac{\hat{N}\sqrt{\hat{N}+\hat{1}}}{\sqrt{\hat{N}+\hat{2}}}$. The explicit expression for the normalized coherent states arises from equation (26):

[TABLE]

where ${}_{p}F_{q}$ is a generalized hypergeometric function defined by

[TABLE]

The mean values for the operators in Eq. (27) and their squares are now:

[TABLE]

In the limit $\alpha\rightarrow 0$ we get $(\sigma_{z})^{2}_{\alpha}(\sigma_{p})^{2}_{\alpha}\rightarrow 4$ (see Figure 4).

As we can see, the Heisenberg uncertainty relation depends strongly on the coherent states under consideration. Thus, for the states in Eq. (29) it takes its minimum at $\alpha=0$, while for those in Eqs. (31) and (33) their maxima are reached at the same point; this is so since the lowest energy eigenstate involved in the corresponding linear combination is different if different families of coherent states are taken into account (see also [24, 25, 26]).

3.3 Magnetic field and probability density

The probability density $\rho=\Psi_{\alpha}^{\dagger}\Psi_{\alpha}$ will be used to analyze the properties of the graphene coherent states. It will depend on the following matrix elements [19]:

[TABLE]

Note also that, according to Eq. (8), it depends on the magnetic field intensity $B_{0}$ through the parameter $\omega$.

3.3.1 Probability density for $f(1)\neq 0$.

A straightforward calculation using Eq. (29) leads to:

[TABLE]

where $\theta=$Arg$[\alpha]$, $r=|\alpha|$. Plots of this probability density for two magnetic field intensities and different $\theta$’s and $r$’s are shown in Figure 5.

3.3.2 Probability density for $f(1)=0$

A. Case with $f(2)\neq 0$.

On the other hand, for the states of Eq. (31) we get (see Figure 6):

[TABLE]

B. Case with $f(2)=0$.

Finally, by employing the states of Eq. (33) we arrive at (see Figure 7)

[TABLE]

3.3.3 Discussion

For the graphene coherent states under study, the probability density reaches a maximum which displaces along the $x$ direction as the parameter $\theta\in[0,\pi/2]$ grows. Also, when $r\rightarrow\infty$ we have $\rho(r)\rightarrow 0$ $\forall\,\theta$.

The previous behavior can be interpreted as follows. According to the mean values of the position and momentum operators, they can be expressed in terms of the complex eigenvalue $\alpha$ in the way:

[TABLE]

where $F_{1,2}(|\alpha|)$ are certain functions that depend on the graphene coherent states under study. In particular, when $\alpha$ is real (for $\theta=0$) we have $\langle\hat{p}\rangle=0$ which means that, on average, the electron moves as many times to the right as to the left, canceling out at the end the positive momentum contributions with the negative ones. Meanwhile, when $\alpha$ is purely imaginary (for $\theta=\pi/2$) we have that $\langle\hat{z}\rangle=0$. This can be interpreted as if the system would perform symmetric oscillations around the equilibrium position $z_{0}$ (or potential center), which is determined by the magnetic field intensity.

On the other hand, when $B_{0}$ increases the maximum of the probability density also grows up while their width decreases (due to probability conservation). This means that the electron is to be found in a more bounded region as $B_{0}$ grows. An opposite interpretation can be formulated when $B_{0}$ decreases.

3.4 Mean energy value

The mean energy value, $\langle\hat{H}\rangle_{\alpha}$, is another quantity useful to characterize the graphene coherent states.

According to the expansion in Eq. (20), where $\Psi_{n}$ are the eigenfunctions of the Dirac-Weyl Hamiltonian (see Eq. (2)), for the graphene coherent states we have that

[TABLE]

The mean energy value is calculated for each coherent state using this expression, which leads to the following results.

3.4.1 $\langle\hat{H}\rangle_{\alpha}$ for $f(1)\neq 0$.

For the coherent states of Eq. (29) we obtain (see Figure 8):

[TABLE]

3.4.2 $\langle\hat{H}\rangle_{\alpha}$ for $f(1)=0$

A. Case with $f(2)\neq 0$.

For the coherent states of Eq. (31) we get (see Figure 9):

[TABLE]

B. Case with $f(2)=0$.

Finally, for the coherent states of Eq. (33) we arrive at (see Figure 10):

[TABLE]

While Figures 8 and 9 have a similar qualitative behavior for $\langle\hat{H}\rangle_{\alpha}$ as function of $\alpha$, Figure 10 shows that the mean energy value for the states of Eq. (33) grows more slowly than the previous ones. These differences depend once again on the structure of the coherent states taken into account. In addition, according to Eq. (39) the mean energy value depends as well on the magnetic field intensity as $\sqrt{B_{0}}$.

4 Conclusions

Dirac electrons in graphene placed in homogeneous magnetic fields which are orthogonal to the material surface are ideal systems to start implementing the coherent states treatment in solid state physics. In particular, for constant magnetic fields the problem has been addressed for the first time quite recently [27]. In fact, in [27] the same physical configuration of this paper was considered, with the assumption that the magnetic field strenght is strong, in order that the Dirac electron stays always at the $n=0$ Landau energy level. On the other hand, in this article we are supposing that the magnetic field strenght is not strong, so that the state of the electron can be a coherent linear combination of all the eigenstates for the Landau energy levels. That is the reason why in this paper we required first to identify the appropriate annihilation and creation operators, in order to build then the coherent states as eigenstates of the former operator. Due to its non-uniqueness, however, it was possible to build different sets of coherent states. Although some of them could look similar to the standard coherent states for the harmonic oscillator, our graphene coherent states in general involve generalized hypergeometric functions. This dependence is more apparent when calculating the Heisenberg uncertainty relation for each set of this paper. This uncertainty achieves a minimum, equal to 1/4, for the coherent states of Eq. (29), since the ground state $\Psi_{0}$ is involved in this linear combination, while it reaches a maximum for the coherent states of Eqs. (31) and (33), depending on the minimum excited state energy involved in the corresponding linear combination (see Figures 2-4).

It is important to remark that, in a sense, the graphene coherent states remind the multiphoton coherent states [28, 29, 30, 31, 32, 33], which appear from realizations of the Polynomial Heisenberg Algebras (PHA) for the harmonic oscillator [24, 25, 26, 34, 35, 36]. In that formalism, the Hilbert space decomposes as a direct sum of $m$ orthogonal subspaces, on each of which it is possible to construct the corresponding coherent states as superpositions of standard coherent states, while in the case of this paper the minimum energy states can be isolated from the remaining Hilbert subspace, depending on the values taken by $f(n)$.

On the other hand, the analysis of the probability density allows to characterize some physical properties of the graphene coherent states. This function indicates that the description for these states remains simple for finite $r$, whatever the value of the parameter $\theta$ is. However, the probability density reaches a maximum whose position along the $x$ axis actually depends on $\theta$ (see Figures 5-7). Meanwhile, the behavior of the mean energy value suggests the possibility of using the graphene coherent states in semi-classical treatments.

Finally, it is important to stress that the non-uniqueness of the annihilation operator leaves open the possibility of exploring more complicated expressions for this operator. As a consequence, plenty of new sets of coherent states can be generated; some of them could be more useful than others for describing interesting physical phenomena in graphene and other carbon allotropes (see e.g. [37, 38, 39]).