**2.2. The equation for the envelope wavefunction**

To describe some phenomena in graphene-based heterostructures, an equation containing a mass term should be written for the envelope wavefunction. A bandgap opening in the energy spectrum of graphene results from the lack of symmetry between the two triangular sublattices of its hexagonal lattice. The corresponding tight-binding Hamiltonian taking into account nearest-neighbor hopping has the form [34]

$$\begin{split} \widehat{H} &= -t \sum\_{\mathbf{B},i,\sigma} \left[ a\_{\sigma}^{\dagger}(\mathbf{B} + \mathbf{d}\_{i}) b\_{\sigma}(\mathbf{B}) + b\_{\sigma}^{\dagger}(\mathbf{B}) a\_{\sigma}(\mathbf{B} + \mathbf{d}\_{i}) \right] \\ &+ \Delta \sum\_{\mathbf{B},\sigma} \left[ a\_{\sigma}^{\dagger}(\mathbf{B} + \mathbf{d}\_{1}) a\_{\sigma}(\mathbf{B} + \mathbf{d}\_{1}) - b\_{\sigma}^{\dagger}(\mathbf{B}) b\_{\sigma}(\mathbf{B}) \right], \end{split} \tag{2}$$

where *t* ≈ 2.8 eV is the nearest-neighbor hopping energy; the sum runs over the position vectors **<sup>B</sup>** of all *<sup>B</sup>* sublattice atoms; *<sup>σ</sup>* <sup>=</sup>↑, <sup>↓</sup> is the (pseudo)spin index; *<sup>a</sup><sup>σ</sup>* (*a*† *<sup>σ</sup>*) and *<sup>b</sup><sup>σ</sup>* (*b*† *σ*) are the annihilation (creation) operators of *A* and *B* sublattice electrons, respectively (see Fig. 1); and the parameter ∆ quantifies the on-site energy difference between the two sublattices (setting ∆ = 0 restores the symmetry between sublattices so that graphene becomes gapless, whereas nonzero ∆ equals the half-gap width in gapped graphene as shown below).

<sup>1</sup> Massless states can be characterized by two quantum numbers: helicity and sign of energy or helicity and eigenvalue of the operator *iγ*<sup>5</sup> [21]. Parity is analogous to the eigenvalue of *iγ*5.

**Figure 1.** Part of hexagonal lattice, with highlighted vectors **d***<sup>i</sup>* from a *B* sublattice atom to the three nearest-neighbor *A* sublattice atoms.

Performing a Fourier transform, we change to the momentum representation

Massless states with opposite helicities are decoupled [20]. In addition, charge carriers have chiral symmetry (helicity is conserved), and parity can be defined for both massless and massive carriers1. In other words, a higher symmetry of massless charge carriers implies the existence of an additional quantum number: helicity. Whereas parity distinguishes only between the valleys where carrier states belong (*λ* = +1 and −1 for states close to the *K* and *K* points, respectively), helicity differs between a particle (electron) and an antiparticle (hole). However, chiral symmetry is broken for massive charge carriers. So, helicity is not a good quantum number any longer. Carrier states in a planar heterostructure combining

Recall that the Dirac equation describing massless carriers in graphene in terms of 4×4 matrices is derived by assuming that they are spinless and have two valley degrees of freedom [10]. When analysis is restricted to charge carriers in one valley, the Dirac equation can be reduced to a 2×2 matrix representation by Weyl's equation for a massless fermion analogous to neutrino in two Euclidean dimensions [22]. The carrier energy spectrum with a pseudospin splitting in a planar heterostructure combining gapless and gapped graphene cannot be correctly analyzed in the 2×2 representation. For similar reasons, the representation of the Dirac algebra in terms of 2×2 matrices is not sufficient for describing

Using the two-dimensional 4×4 Dirac equation to describe charge carriers in a graphene-based nanostructure, we can study pseudospin effects following an approach to narrow-gap semiconductor heterostructures based on the Dirac model [24]. This makes methods developed for solving problems in the spintronics of narrow-gap semiconductor

To describe some phenomena in graphene-based heterostructures, an equation containing a mass term should be written for the envelope wavefunction. A bandgap opening in the energy spectrum of graphene results from the lack of symmetry between the two triangular sublattices of its hexagonal lattice. The corresponding tight-binding Hamiltonian taking into

*<sup>σ</sup>*(**<sup>B</sup>** + **<sup>d</sup>***i*)*bσ*(**B**) + *<sup>b</sup>*†

*<sup>σ</sup>*(**<sup>B</sup>** <sup>+</sup> **<sup>d</sup>**1)*aσ*(**<sup>B</sup>** <sup>+</sup> **<sup>d</sup>**1) <sup>−</sup> *<sup>b</sup>*†

where *t* ≈ 2.8 eV is the nearest-neighbor hopping energy; the sum runs over the position

are the annihilation (creation) operators of *A* and *B* sublattice electrons, respectively (see Fig. 1); and the parameter ∆ quantifies the on-site energy difference between the two sublattices (setting ∆ = 0 restores the symmetry between sublattices so that graphene becomes gapless,

<sup>1</sup> Massless states can be characterized by two quantum numbers: helicity and sign of energy or helicity and eigenvalue

vectors **<sup>B</sup>** of all *<sup>B</sup>* sublattice atoms; *<sup>σ</sup>* <sup>=</sup>↑, <sup>↓</sup> is the (pseudo)spin index; *<sup>a</sup><sup>σ</sup>* (*a*†

whereas nonzero ∆ equals the half-gap width in gapped graphene as shown below).

*<sup>σ</sup>*(**B**)*aσ*(**B** + **d***i*)

*<sup>σ</sup>*(**B**)*bσ*(**B**)

 , (2)

*σ*)

*<sup>σ</sup>*) and *<sup>b</sup><sup>σ</sup>* (*b*†

gapless and gapped graphene should be characterized by parity.

the chiral symmetry breaking in (2+1)QED [23].

182 Graphene - New Trends and Developments

heterostructures applicable to graphene-based ones [25–33].

**2.2. The equation for the envelope wavefunction**

account nearest-neighbor hopping has the form [34]

*<sup>H</sup>* = − *<sup>t</sup>* ∑

**B**,*i*,*σ*

+ ∆ ∑ **B**,*σ a*†

of the operator *iγ*<sup>5</sup> [21]. Parity is analogous to the eigenvalue of *iγ*5.

 *a*†

$$a\_{\sigma}(\mathbf{A}) = \int\_{\Omega\_{\mathcal{B}}} \frac{d^2 k}{(2\pi)^2} a\_{\sigma}(\mathbf{k}) \varepsilon^{i\mathbf{k}\cdot\mathbf{A}} \, \, b\_{\sigma}(\mathbf{B}) = \int\_{\Omega\_{\mathcal{B}}} \frac{d^2 k}{(2\pi)^2} b\_{\sigma}(\mathbf{k}) \varepsilon^{i\mathbf{k}\cdot\mathbf{B}} \, \,\_2$$

where Ω*<sup>B</sup>* means integration over the first Brillouin zone (it is shown on Fig. 2). Hamiltonian (2) is rewritten as

$$
\hat{H} = \sum\_{\sigma} \int\_{\Omega\_{\mathcal{B}}} \frac{d^2 k}{(2\pi)^2} \left( a\_{\sigma}^\dagger(\mathbf{k}) \, b\_{\sigma}^\dagger(\mathbf{k}) \right) \begin{pmatrix} \Delta & -t \sum\_{\stackrel{i}{\hat{i}}} e^{-i\mathbf{k} \cdot \mathbf{d}\_{\hat{i}}} \\ -t \sum\_{\stackrel{i}{\hat{i}}} e^{i\mathbf{k} \cdot \mathbf{d}\_{\hat{i}}} & -\Delta \end{pmatrix} \begin{pmatrix} a\_{\sigma}(\mathbf{k}) \\ b\_{\sigma}(\mathbf{k}) \end{pmatrix}. \tag{3}
$$

Conduction and valence band extrema lie at the corners of the first Brillouin zone. In the case of gapless graphene, they touch each other and there are the cone-like energy surfaces at the *K* and *K* points (see Fig. 3). We use Hamiltonian (3) expanded around the *K* point with quasimomentum **q**<sup>1</sup> = 2*<sup>π</sup>* <sup>3</sup>*<sup>a</sup>* , <sup>2</sup>*<sup>π</sup>* 3 <sup>√</sup>3*<sup>a</sup>* or around the *K* point with **q**<sup>2</sup> = 2*<sup>π</sup>* <sup>3</sup>*<sup>a</sup>* , <sup>−</sup> <sup>2</sup>*<sup>π</sup>* 3 <sup>√</sup>3*<sup>a</sup>* 

$$
\hat{H} = \sum\_{\sigma} \int \frac{d^2k}{(2\pi)^2} \hat{\Psi}\_{\sigma}^{\dagger}(\mathbf{k}) \hat{\mathcal{H}} \hat{\Psi}\_{\sigma}(\mathbf{k}).
$$

where integration is performed over small neighbor-hoods of the *K* and *K* points. Near the corners, the Hamiltonian reduces to

$$
\hat{\mathcal{H}} = \begin{pmatrix}
v\_F \boldsymbol{\sigma} \cdot \mathbf{k} + \Delta \sigma\_z & 0 \\
0 & v\_F \boldsymbol{\sigma}^\* \cdot \mathbf{k} + \Delta \sigma\_z
\end{pmatrix} \tag{4}
$$

**Figure 2.** The first Brillouin zone of graphene, with linear energy spectrum at the corners (Dirac points). The reciprocal lattice vectors **b**<sup>1</sup> = (2*π*/3*a*, 2*π*/ <sup>√</sup>3*a*) and **<sup>b</sup>**<sup>2</sup> = (2*π*/3*a*, <sup>−</sup>2*π*/ <sup>√</sup>3*a*), where *<sup>a</sup>*=1.42 Å is the lattice spacing, combined with the dashed lines equivalently represent the first Brillouin zone as a rhombus.

where *vF* = <sup>3</sup> <sup>2</sup> *at* is the carrier Fermi velocity, σ = (*σx*, *σy*) and σ<sup>∗</sup> = (*σx*, −*σy*) are Pauli matrices in the sublattice space, and Ψ *<sup>σ</sup>*(**k**) is the bispinor defined as

$$
\hat{\Psi}\_{\sigma}(\mathbf{k}) = \begin{pmatrix}
\hat{\Psi}\_{\sigma}^{(1)}(\mathbf{k}) \\
\hat{\Psi}\_{\sigma}^{(2)}(\mathbf{k})
\end{pmatrix},
$$

in terms of

$$
\hat{\Psi}\_{\sigma}^{(1,2)}(\mathbf{k}) = \exp\left(\frac{5\pi i}{12}\sigma\_z\right)\sigma\_z \begin{pmatrix} a\_{\sigma}(\mathbf{q}\_{1,2} + \mathbf{k}) \\ b\_{\sigma}(\mathbf{q}\_{1,2} + \mathbf{k}) \end{pmatrix}.
$$

We write an equation for the envelope wavefunction in a planar heterostructure:

$$
\sigma\_{\mathbf{F}\mathbf{j}} \left( \tau\_{\mathbf{0}} \otimes \sigma\_{\mathbf{x}} \hat{\mathbf{p}}\_{\mathbf{x}} + \tau\_{\mathbf{z}} \otimes \sigma\_{\mathbf{y}} \hat{\mathbf{p}}\_{\mathbf{y}} \right) + \tau\_{\mathbf{0}} \otimes \sigma\_{\mathbf{z}} \Delta\_{\mathbf{j}} + \tau\_{\mathbf{0}} \otimes \sigma\_{\mathbf{0}} \left( V\_{\mathbf{j}} - E \right) \Big| \, \mathbf{\varPsi}(\mathbf{x}, \mathbf{y}) = \mathbf{0}. \tag{5}
$$

Here, ∆*<sup>j</sup>* = *Egj*/2 (*j* = 1, 2, . . .) denotes half-width of bandgap; the respective work functions *Vj* depend on the mid-gap energies relative to the Dirac points for the corresponding materials; the 2×2 unit matrix *σ*<sup>0</sup> acts in the sublattice space; the 2×2 unit matrix *τ*<sup>0</sup> and the matrix *τ<sup>z</sup>* defined similar to the Pauli matrix *σ<sup>z</sup>* act in the valley space; ⊗ is the Kronecker product symbol; and *<sup>p</sup><sup>x</sup>* <sup>=</sup> <sup>−</sup>*i∂<sup>x</sup>* and *<sup>p</sup><sup>y</sup>* <sup>=</sup> <sup>−</sup>*ih*¯ *<sup>∂</sup><sup>y</sup>* are the momentum operator components. Assuming that the carrier Fermi velocities may differ between the regions, we denote those for gapped *j*th regions by *vFj* for both gapless graphene and its gapped modifications.

Charge carriers move freely along the *y* axis:

$$\Psi(x,y) = \Psi(x)e^{ik\_y y}.$$

The wavefunction Ψ(*x*) is a bispinor:

**Figure 2.** The first Brillouin zone of graphene, with linear energy spectrum at the corners (Dirac points). The reciprocal

 <sup>Ψ</sup>(1) *<sup>σ</sup>* (**k**) <sup>Ψ</sup>(2) *<sup>σ</sup>* (**k**)

+ *τ*<sup>0</sup> ⊗ *σz*∆*<sup>j</sup>* + *τ*<sup>0</sup> ⊗ *σ*<sup>0</sup>

Here, ∆*<sup>j</sup>* = *Egj*/2 (*j* = 1, 2, . . .) denotes half-width of bandgap; the respective work functions *Vj* depend on the mid-gap energies relative to the Dirac points for the corresponding materials; the 2×2 unit matrix *σ*<sup>0</sup> acts in the sublattice space; the 2×2 unit matrix *τ*<sup>0</sup> and the matrix *τ<sup>z</sup>* defined similar to the Pauli matrix *σ<sup>z</sup>* act in the valley space; ⊗ is the Kronecker product symbol; and *<sup>p</sup><sup>x</sup>* <sup>=</sup> <sup>−</sup>*i∂<sup>x</sup>* and *<sup>p</sup><sup>y</sup>* <sup>=</sup> <sup>−</sup>*ih*¯ *<sup>∂</sup><sup>y</sup>* are the momentum operator components. Assuming that the carrier Fermi velocities may differ between the regions, we denote those for gapped *j*th regions by *vFj* for both gapless graphene and its gapped modifications.

Ψ(*x*, *y*) = Ψ(*x*)*e*

*iky <sup>y</sup>*.

<sup>2</sup> *at* is the carrier Fermi velocity, σ = (*σx*, *σy*) and σ<sup>∗</sup> = (*σx*, −*σy*) are Pauli

 ,

*aσ*(**q**1,2 + **k**) *bσ*(**q**1,2 + **k**)

> *Vj* − *E*

 .

Ψ(*x*, *y*) = 0. (5)

<sup>√</sup>3*a*), where *<sup>a</sup>*=1.42 Å is the lattice spacing, combined

<sup>√</sup>3*a*) and **<sup>b</sup>**<sup>2</sup> = (2*π*/3*a*, <sup>−</sup>2*π*/

Ψ *<sup>σ</sup>*(**k**) =

5*πi* <sup>12</sup> *<sup>σ</sup><sup>z</sup> σz* 

We write an equation for the envelope wavefunction in a planar heterostructure:

with the dashed lines equivalently represent the first Brillouin zone as a rhombus.

<sup>Ψ</sup>(1,2) *<sup>σ</sup>* (**k**) = exp

*<sup>τ</sup>*<sup>0</sup> <sup>⊗</sup> *<sup>σ</sup><sup>x</sup> <sup>p</sup><sup>x</sup>* <sup>+</sup> *<sup>τ</sup><sup>z</sup>* <sup>⊗</sup> *<sup>σ</sup><sup>y</sup> <sup>p</sup><sup>y</sup>*

Charge carriers move freely along the *y* axis:

matrices in the sublattice space, and Ψ *<sup>σ</sup>*(**k**) is the bispinor defined as

lattice vectors **b**<sup>1</sup> = (2*π*/3*a*, 2*π*/

184 Graphene - New Trends and Developments

where *vF* = <sup>3</sup>

in terms of

 *vFj* 

$$\Psi(x) = \begin{pmatrix} \psi\_K(x) \\ \psi\_{K'}(x) \end{pmatrix}.$$

where the spinors *ψK*(*x*) and *ψ<sup>K</sup>*(*x*) represent charge carriers in the *K* and *K* valleys, respectively:

$$\psi\_K(\mathfrak{x}) = \begin{pmatrix} \psi\_{KA}(\mathfrak{x}) \\ \psi\_{KB}(\mathfrak{x}) \end{pmatrix}, \quad \psi\_{K'}(\mathfrak{x}) = \begin{pmatrix} \psi\_{K'A}(\mathfrak{x}) \\ \psi\_{K'B}(\mathfrak{x}) \end{pmatrix}.$$

In the present context, the parity operator is expressed as follows:

$$
\widehat{P} = \mathfrak{r}\_{\mathbb{Z}} \otimes \sigma\_0. \tag{6}
$$

Equation (5) is solved here in the parity basis. The eigenfunctions Ψ*λ*(*x*) of parity operator (6) are defined as follows:

$$\begin{aligned} P\Psi\_{\lambda}(\mathbf{x}) &= \lambda \Psi\_{\lambda}(\mathbf{x}), \\\\ \Psi\_{+1}(\mathbf{x}) &= \begin{pmatrix} \Psi\_{+1,K}(\mathbf{x}) \\ 0 \end{pmatrix}, \\ \Psi\_{-1}(\mathbf{x}) &= \begin{pmatrix} 0 \\ \Psi\_{-1,K'}(\mathbf{x}) \end{pmatrix}. \end{aligned} \tag{7}$$

Rewriting Eq. (5) as the 2×2 matrix equations

$$\left(-i\upsilon\_{Fj}\sigma\_{\mathbf{x}}\partial\_{\mathbf{x}} + \upsilon\_{Fj}k\_{\mathbf{y}}\sigma\_{\mathbf{y}} + \lambda\Delta\_{\mathbf{j}}\sigma\_{\mathbf{z}} + V\_{\mathbf{j}}\right)\psi\_{\lambda K}(\mathbf{x}) = E\_{\lambda}\psi\_{\lambda K}(\mathbf{x}),\tag{8}$$

$$\left(-i\upsilon\_{Fj}\sigma\_{\mathbf{x}}\partial\_{\mathbf{x}} - \upsilon\_{Fj}k\_{\mathbf{y}}\sigma\_{\mathbf{y}} - \lambda\Delta\_{\mathbf{j}}\sigma\_{\mathbf{z}} + V\_{\mathbf{j}}\right)\psi\_{\lambda K'}(\mathbf{x}) = E\_{\lambda}\psi\_{\lambda K'}(\mathbf{x}).\tag{9}$$

We see that setting ∆*<sup>j</sup>* = 0 and *Vj* = 0 brings us back to the spinor wavefunctions describing chiral states near the *K* or *K* point, where the operator *h* can be defined. However, chiral symmetry is broken when ∆ = 0 (see previous section). Defining parity *λ* as the eigenvalue of operator (6), we find that it indicates the valley where charge carriers belong: by virtue of (7), *λ* = +1 for the states near the *K* point described by Eq. (8) and *λ* = −1 for the states near the *K* point described by Eq. (9).

In both gapped and gapless graphene, the valleys transform into each other under time reversal. This is indicated by the opposite signs of the terms proportional to *ky* in Eqs. (8) and (9), since *ky* → −*ky* under time reversal. It can be shown directly by using the time reversal operator T in explicit form that *λ* → −*λ* under T . Indeed, if *ky* is parallel to the line *K* − *M* − *K* (see Fig. 2) and its origin is set at *M*, then *K* → *K* and *K* → *K* under T .

**Figure 3.** The energy surface for gapless graphene.

Equations (8) and (9) are equivalently rewritten as the 2×2 matrix equation

$$\left(-i\upsilon\_{\rm Fj}\sigma\_{\rm x}\partial\_{\rm x} + \lambda\upsilon\_{\rm Fj}k\_{\rm y}\sigma\_{\rm y} + \Delta\_{\rm j}\sigma\_{\rm z} + V\_{\rm j}\right)\psi\_{\lambda}(\mathbf{x}) = E\_{\lambda}\psi\_{\lambda}(\mathbf{x}).\tag{10}$$

Hereinafter, valley indices *K* and *K* are omitted as unnecessary since *λ* specifies the valley where charge carriers belong.

#### **2.3. The boundary conditions**

Now, let us discuss the boundary conditions at the interfaces between different graphene materials. At the outset, we note that they are easier to formulate than those used at the graphene–free-space interface in models of edge states [6, 35]. To derive boundary conditions in the present model, we must find a relation between *ψλ*(*l*) and *ψλ*(−*l*) as *l* → 0 in the neighbor-hood of x = 0, where *l* goes down to an atomic scale (condition at x = d is derived similarly). Multiplying Eq. (10) by *<sup>ψ</sup>*†(*x*) on the left, we integrate it over [−*l*, *<sup>l</sup>*]. Since a is small, we neglect all terms except those containing a derivative with respect to *x* to obtain<sup>2</sup>

$$
\psi\_{
\lambda}^{(-)\dagger}(-l)\upsilon\_F^{(-)}\psi\_{
\lambda}^{(-)}(-l) = \psi\_{
\lambda}^{(+)\dagger}(l)\upsilon\_F^{(+)}\psi\_{
\lambda}^{(+)}(l).
$$

where *<sup>ψ</sup>*(−) *<sup>λ</sup>* and *<sup>ψ</sup>*(+) *<sup>λ</sup>* are defined on the left- and right-hand sides of the boundary (at *x <* 0 and *x >* 0, respectively). Representing these functions as

$$\psi\_{\lambda}^{(\pm)} = \left| \psi\_{\lambda}^{(\pm)} \right| \exp\left(i\varphi^{(\pm)}\right) \cdot \lambda$$

we rewrite the equality above as

$$\sqrt{\upsilon\_F^{(-)}} \left| \psi\_\lambda^{(-)\dagger}(-l) \right| = \sqrt{\upsilon\_F^{(+)}} \left| \psi\_\lambda^{(+)\dagger}(l) \right|.$$

<sup>2</sup> From the given equality, we have the continuity of the current component normal to the interface in the heterostructure plane. It is necessary condition.

To formulate the boundary condition in final form, we assume that the difference between phases of *<sup>ψ</sup>*(−) *<sup>λ</sup>* and *<sup>ψ</sup>*(+) *<sup>λ</sup>* near the interface is a multiple of 2*π*:

$$
\varphi^{(+)} = \varphi^{(-)} + 2\pi n, \quad n \in \mathbb{Z}.
$$

As *l* goes to zero, we obtain the following wavefunction-matching condition [26, 27]

$$
\sqrt{v\_F^{(-)}} \psi\_\lambda^{(-)} = \sqrt{v\_F^{(+)}} \psi\_\lambda^{(+)}\, \, \, \, \, \tag{11}
$$

where minus and plus signs refer to the materials on the left- and right-hand sides of the boundary, respectively.
