# Ab Initio Study of Graphene/hBN Van der Waals Heterostructures: Effect of Electric Field, Twist Angles and p-n Doping on the Electronic Properties

Simone Brozzesi, Claudio Attaccalite, Francesco Buonocore, Giacomo Giorgi, Maurizia Palummo, Olivia Pulci

<figure><img src="image_1.png" /></figure>

To cite this version:

Simone Brozzesi, Claudio Attaccalite, Francesco Buonocore, Giacomo Giorgi, Maurizia Palummo, et al. Ab Initio Study of Graphene/hBN Van der Waals Heterostructures: Effect of Electric Field, Twist Angles and p-n Doping on the Electronic Properties. Nanomaterials, 2022, 12 (12), pp.2118. 10.3390/nano12122118. [https://doi.org/10.3390/nano12122118](https://doi.org/10.3390/nano12122118)

HAL Id: hal-03700706

https://hal.science/hal-03700706

Submitted on 21 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Ab Initio Study of Graphene/hBN Van der Waals Heterostructures: Effect of Electric Field, Twist Angles and p-n Doping on the Electronic Properties

Simone Brozzesi <fnref n="1" />,* , Claudio Attaccalite <fnref n="2" /> , Francesco Buonocore <fnref n="3" /> , Giacomo Giorgi <fnref n="4" />,<fnref n="5" /> , Maurizia Palummo <fnref n="1" /> and Olivia Pulci <fnref n="1" />,*

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Brozzesi, S.; Attaccalite, C.; Buonocore, F.; Giorgi, G.; Palummo, M.; Pulci, O. Ab Initio Study of Graphene/hBN Van der Waals Heterostructures: Effect of Electric Field, Twist Angles and p-n Doping on the Electronic Properties.
*Nanomaterials* **2022**, *12*, 2118. https://doi.org/10.3390/nano12122118

Academic Editors: Béla Pécz, Carlos Miguel Costa and Guowei Yang

Received: 17 May 2022

Accepted: 16 June 2022

Published: 20 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2022 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** In this work, we study the structural and electronic properties of boron nitride bilayers sandwiched between graphene sheets. Different stacking, twist angles, doping, as well as an applied external gate voltage, are reported to induce important changes in the electronic band structure near the Fermi level. Small electronic lateral gaps of the order of few meV can appear near the Dirac points K. We further discuss how the bandstructures change applying a perpendicular external electric field, showing how its application lifts the degeneracy of the Dirac cones and, in the twisted case, moves their crossing points away from the Fermi energy. Then, we consider the possibility of co-doping, in an asymmetric way, the two external graphene layers. This is a situation that could be realized in heterostructures deposited on a substrate. We show that the co-doping acts as an effective external electric field, breaking the Dirac cones degeneracy. Finally, our work demonstrates how, by playing with field strength and p-n co-doping, it is possible to tune the small lateral gaps, pointing towards a possible application of C/BN sandwich structures as nano-optical terahertz devices.

**Keywords:** graphene; h-BN; Van der Waals heterostructures; density functional theory; THz nanodevices

## 1. Introduction

The exfoliation of graphene, a single layer of graphite sheet in 2004 and the successive Nobel Prize in Physics to Geim and Novoselov (2010), mark the birth of the 2D materials (2DM) age [1–4]. Since then, a multitude of one or few layers atomic thick materials, with different mechanical, electronic, magnetic, and optical properties is under focus, being interesting both at a fundamental and applicative point of view [5–13] Beyond the study of 2DM monolayers [14], a very exciting development in the field is the possibility to stack them, like the bricks of an atomic-scale Lego play, in artificial Van der Waals heterostructures (vdW HTs), with strong covalent bonds within the atomic layers and weak vdW interaction among them, whose resulting properties can be, in principle, controlled, tuned, and manipulated on demand [15,16]. In this context, vdW HTs composed of graphene and other 2DM are largely investigated. Indeed, despite graphene having unique chemical and physical properties, where surely the most peculiar one is the presence of Dirac-like linear

<footnote n="1">1 Dipartimento di Fisica and INFN, Università di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Rome, Italy; maurizia.palummo@roma2.infn.it</footnote>

<footnote n="2">2 Centre Interdisciplinaire de Nanoscience de Marseille UMR 7325 Campus de Luminy, CNRS/Aix-Marseille Université, CEDEX 9, 13288 Marseille, France; attaccalite@cinam.univ-mrs.fr</footnote>

<footnote n="3">3 Energy Technologies and Renewable Sources (TERIN) Department, Italian National Agency for New Technologies, Energy and Sustainable Economic Development (ENEA), Casaccia Research Centre, 00123 Rome, Italy; francesco.buonocore@enea.it</footnote>

<footnote n="4">4 Department of Civil & Environmental Engineering (DICA), Università degli Studi di Perugia, Via G. Duranti 93, 06125 Perugia, Italy; giacomo.giorgi@unipg.it</footnote>

<footnote n="5">5 CNR-SCITEC, 06123 Perugia, Italy</footnote>

<footnote>* Correspondence: simone.brozzesi@roma2.infn.it (S.B.); olivia.pulci@roma2.infn.it (O.P.); Tel.: +39-06-7259-4548 (O.P.)</footnote>
band dispersion, its semi-metallic nature prevents the use in several nanoelectronics and photonics' devices. The realization of HTs [17] is then one of the explored strategies to open the gap of graphene but to keep, as much as possible, high electronic carrier mobility values [18–20]. Due to similar atomic structure and almost perfect lattice match [21,22], hexagonal boron nitride (h-BN) has been used in combination with graphene to build up new HTs, from bilayer to multilayer forms [23–29]. With h-BN being an insulator with a wide band-gap of about 6 eV [30–32] and due to the large difference in the work functions of the two materials, a type-I band-alignment occurs. The resulting electronic properties are roughly the sum of those of the two separated systems, with low-energy graphene states falling always in the gap of h-BN. Nevertheless, as we show here and consistently with existing literature, the specific features near the Fermi level are strongly influenced by the specific stacking among graphene and h-BN sheets.

Motivated by interesting physical properties such as being a good platform for exciton condensation [33], tunable metal-insulator transition [34], observation of highly confined and topological plasmons modes [27,35] and possible applications in opto-electronic terahertz devices [36–38], we focus here on graphene/h-BN HTs, where two boron-nitride layers are inserted between two graphene ones. Our aim is to discuss how stacking order, twisting angle, external vertical gate voltage, and doping of graphene layers modify the electronic properties with respect to those of an isolated graphene sheet. The results presented here are obtained by means of first-principles Density Functional Theory (DFT) calculations. The analysis of the optical properties, as well as of the role of many-body effects, such as electronic self-energy and e-h interaction, will be the subject of future work.

## 2. Materials and Methods

All DFT calculations reported in this manuscript are obtained using the Quantum ESPRESSO suite [39–41] with plane wave expansion and scalar-relativistic norm-conserving pseudopotentials [42–44]. A kinetic energy cutoff of 60 Ry is adopted. PBE exchange-correlation [45] with vdW-DF2-b86r [46,47] exchange functional is exploited in order to take into account Van der Waals interaction. For the sake of comparison, we also perform several calculations of the structural properties using LDA exchange-correlation functional. A uniform Monkhorst-Pack k-point mesh of 60 × 60 × 1 is used for the non twisted systems, whereas a mesh of 12 × 12 × 1 is used for twisted HT. We introduce a vacuum of 30 Å in the non-periodic direction (z) to ensure periodic replicas' decoupling. In order to simulate *n* and *p* doping of the graphene layers, a set of modified pseudopotentials are used, where C (Z = 4) atoms are substituted by *C*<sub>*p*</sub> pseudo-atoms with Z = 3.98 (p doping) or *C*<sub>*n*</sub> atoms with Z = 4.02 (n doping). In this way, we introduce an infinitesimal electronic negative charge in the top graphene layer (+0.02 electrons for each *C*<sub>*n*</sub> atom) and an excess positive charge in the bottom one (−0.02 electrons for each *C*<sub>*p*</sub> atom), still preserving the total charge neutrality of the system. In the non-twisted case, all the C atoms in the unit cell have been modified. In the twisted case, only four of the twenty-eight C atoms in the unit cell have been substituted (two with *C*<sub>*p*</sub> atoms in the bottom graphene layer, and two with *C*<sub>*n*</sub> in the top graphene layer). This choice determines different doping levels in the two cases.

Structure relaxation is assumed at convergence when the maximum component of the residual forces on the ions is smaller than $10^{-4}$ Ry/Bohr. Thanks to the small lattice mismatch of ~1.5% between the two separated systems, we devise all the considered HTs using as a starting lattice parameter the average of the two lattice parameters and then fully relaxing the cell to find the equilibrium in-plane *a* lattice constant and inter-planar distances. Results of the structure optimization calculations confirm that the average of the lattice parameter was a valid guess.

## 3. Results

### 3.1. *Structural Properties and Stacking Energies*

In this section, we present the structural and energetic properties of three designed C-BN-BN-C quadrilayers, whereas the electronic properties are discussed in the next one.
The atomic structures are shown in Figure 1, where the selected nomenclatures, AB-AA'-AB, AB-AA'-AB', and T-AA'-T (T stands for twisted), are based on the stacking order of the central BN-BN and external C-BN interfaces. In the AA and AA' configuration, a hexagon in one plane is on top of a hexagon in the plane below, whereas, in the AB Bernal type configuration, the hexagons are shifted, with one atom in the upper plane aligned to the center of the hexagon in the lower plane. Since 2D-hBN has two different atoms in the unit cell, two AA patterns are feasible: simple AA stacking, with B on B and N on N, and AA', with B on N and vice versa. Preliminary stability analysis of the bilayers was performed in order to determine the lower energy configuration for HT interfaces. In particular, we calculate the stacking energies Δ*E*^st of bilayers and of quadrilayers, with Δ*E*^st = *E*^tot - Σᵢ*E*ᵢ, where *E*^tot and *E*ᵢ are the HT and *i*-th mono-layer total energy, respectively [48–52].

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Crystal structures of C-BN-BN-C quadrilayers: (a,d) top and side view of AB-AA'-AB stacking. (b,e) the same for AB-AA'-AB' stacking; (c,f) top and side view of the 21.8° twisted case T-AA'-T. In the side views, the stacked atoms are circled in black, whereas the red dots indicate the vertical alignment of one atom in the upper layer with the center of the hexagon in the bottom layer.</figcaption></figure>

Our study reveals that AB stacking for C/BN interface and AA' stacking for BN/BN
interface are the most stable configurations. As shown in Table 1, C-BN interface in AB
stacking, with C on top of B, is more stable than AB', with C on top of N, by ~24 meV/atom.
The same analysis for BN interface shows that the AA' is the more stable stacking. For
what concerns the quadrilayers, in the first HT (AB-AA'-AB, left panels of Figure 1), two
completely equivalent C/BN (AB) interfaces are present, half of the carbon atoms of the top
(bottom) graphene layer located above (below) boron atoms, and half at the center of h-BN
hexagons. In the second one (AB-AA'-AB', central panels of Figure 1), two non-equivalent
C/BN interfaces occur, one AB' (top) and the other AB (bottom): while, for the bottom
graphene layer, the situation is similar to the previous case, for the top layer, half carbon
atoms are located above the nitrogen atoms and half at the centers of hexagonal h-BN rings.
<table><caption>**Table 1.** Stacking energy $\Delta E^{st}$ calculated for quadrilayers and bilayers, obtained with PBE + vdW exchange correlation functional.</caption><thead><tr><th>System</th><th>Stacking</th><th>ΔE<sup>st</sup>[meV/atom]</th></tr></thead><tbody><tr><td rowspan="3">C-BN-BN-C</td><td>AB-AA'-AB</td><td>-57.8</td></tr><tr><td>AB-AA'-AB'</td><td>-34.7</td></tr><tr><td>T-AA'-T</td><td>-34.3</td></tr><tr><td>C-BN</td><td>AB</td><td>-42.6</td></tr><tr><td>C-BN</td><td>AB'</td><td>-18.7</td></tr><tr><td>BN-BN</td><td>AA'</td><td>-22.7</td></tr></tbody></table>

The twisted HT (T-AA'-T) is devised rotating clockwise by 21.8° the two external graphene layers with respect to the two internal BN layers and has a $\sqrt{7} \times \sqrt{7}$ unit cell. Each layer contains 14 atoms. In the twisted HT, as evident looking at the right panels of Figure 1, only a limited number of C atoms are below or above B or N, and mixed types of stacking are present. In Table 1, we compare the stacking energies $\Delta E^{st}$ of the three quadrilayers: all considered structures show a gain in energy with respect to their monolayer counterparts, with the AB-AA'-AB HT the most stable, and the twisted the least stable one. $\Delta E^{st}$ of bilayers are also reported for comparison.

The equilibrium structural data of the three quadrilayers, together with those of graphene, BN monolayer, BN-BN, and C-BN bilayers, are listed in Table 2. Structural optimizations confirm the initial guess for lattice parameters of C-BN and C-BN-BN-C as being the average of those of single monolayers. The interplanar distances are slightly affected by the numbers of layers in the HT, especially for internal BN-BN AA' layers. In fact, while graphene/boron-nitride distance is the same in C-BN and C-BN-BN-C, the distance between the two BN layers decreases from 3.30 Å (BN-BN freestanding) to 3.23(3.26) in twisted (AB-AA'-AB') due to the interaction between the p_z orbitals of B/N atoms and those of carbon atoms in the adjacent graphene layer. Our results are in agreement with those reported by Giovannetti et al. [53], where a C-BN-BN 3-layer structure has been considered. Test calculations using LDA exchange-correlation provided in-plane lattice parameters very similar to PBE+vdW but slightly shorter inter-planar distances.

<table><caption>**Table 2.** Lattice parameters and distances between the layers in the systems studied. The lattice constant of the twisted $\sqrt{7} \times \sqrt{7}$ HT (a = 6.55 Å) corresponds to an effective 1 × 1 lattice constant of about 2.47 Å. Calculations are done using PBE-vdW b86r functional.</caption><thead><tr><th>System</th><th>Stacking</th><th>a [Å]</th><th>d<sub>C-BN(bottom)</sub> [Å]</th><th>d<sub>BN-BN</sub> [Å]</th><th>d<sub>C-BN-C(top)</sub> [Å]</th></tr></thead><tbody><tr><td rowspan="3">4L C-BN-BN-C</td><td>AB-AA'-AB</td><td>2.49</td><td>3.23</td><td>3.24</td><td>3.23</td></tr><tr><td>AB-AA'-AB'</td><td>2.48</td><td>3.29</td><td>3.26</td><td>3.42</td></tr><tr><td>T-AA'-T</td><td>6.55</td><td>3.40</td><td>3.23</td><td>3.40</td></tr><tr><td rowspan="3">2L C-BN</td><td>AB</td><td>2.49</td><td>3.28</td><td></td><td></td></tr><tr><td>AB'</td><td>2.48</td><td>3.44</td><td></td><td></td></tr><tr><td>AA'</td><td>2.52</td><td></td><td>3.30</td><td></td></tr><tr><td rowspan="2">1L</td><td>Graphene</td><td>AB</td><td>2.46</td><td></td><td></td></tr><tr><td>hBN</td><td>AA'</td><td>2.50</td><td></td><td></td></tr></tbody></table>

## 3.2. Electronic Properties

The three quadrilayer C-BN-BN-C electronic band structures are reported in Figure 2. In order to highlight the important differences near the Fermi energy, the plots are here zoomed around the Dirac point K, and all the calculations are performed using the same $\sqrt{7} \times \sqrt{7}$ supercell, allowing an immediate comparison.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Band structure of the quadrilayers C-BN-BN-C in (a) non twisted AB-AA'-AB, (b) non twisted AB-AA'-AB', and (c) twisted T-AA'-T configuration. All bands are calculated in a $\sqrt{7} \times \sqrt{7}$ cell.</figcaption></figure>

<table><caption>**Table 3.** Fermi velocities and electron-hole effective masses for the C-BN-BN-C quadrilayers.</caption><thead><tr><th rowspan="2">System</th><th colspan="2"><b>v<sub>F</sub></b><br/>[10<sup>5</sup>m/s]</th><th><b>m<sup>*</sup></b><br/>[10<sup>-3</sup>m<sub>c</sub>]</th></tr><tr><th>M → K</th><th>K → Γ</th><th></th></tr></thead><tbody><tr><td>AB-AA'-AB</td><td>7.71</td><td>8.01</td><td>7.3</td></tr><tr><td>AB-AA'-AB'</td><td>7.82</td><td>8.40</td><td>0.9</td></tr><tr><td>T-AA'-T</td><td>8.20</td><td>8.14</td><td>0.05</td></tr><tr><td>Graphene</td><td>8.30</td><td>8.40</td><td>0</td></tr></tbody></table>

Notably, the behavior is significantly different in the three cases.

The AB-AA'-AB stacking (left panel) shows a clear parabolic dispersion with a band gap of 51 meV and a very small (not visible) degeneracy lift off, both in valence and conduction band-edges. This is due to the (weak) interaction between the graphene layers. The AB-AA'-AB' stacking (central panel) presents an interpenetration of the parabolic bands, with the opening of a gap at K of 0.126 eV and of two small lateral band gaps (≈7 meV) near K. The degeneracy of the bands is lifted. Finally, the bands of the twisted case (right panel) show a linear dispersion even for k-points very near to the Dirac point, with a central gap at K of 7.3 meV and two lateral gaps (at two k-points very close to K) of 0.4 meV. For the three quadrilayers, the Fermi velocity $v_F$ has been estimated, and the effective carrier mass $m^*$ calculated from the energy gap $E_g$ through Equation [14]:

$$m^* = \frac{E_g}{2v_F^2} \qquad (1)$$

The results are shown in Table 3.

In order to explain the different electronic properties of the tetralayers, it is useful to analyze the k-resolved Projected Density of States (k-PDOS) shown in Figure 3 calculated in a 1 × 1 unit cell along the M-K-Γ path, for the AB-AA'-AB HT.

Green and orange points in the figure are associated to $p_z$ and $p_{xy}$ orbitals of the carbon atoms, whereas red, blue, cyan, and grey ones are due to $p_z$ and $p_{xy}$ orbitals of B and N atoms.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** (a) k-resolved projected density of states (k-DOS) for the AB-AA'-AB HT. Colored dots indicate atomic orbitals; (b) zoom of the k-DOS near of the Dirac cone, after subtraction of the C states contribution. In both of the plots, the size of the dots is proportional to the weight of the contribution. In (b), the dot size is enlarged by a 200 factor in order to make them clearly visible.</figcaption></figure>

From Figure 3a, it is clear that HT's band structure can be roughly considered as the sum of graphene and BN-BN AA' bilayer bandstructures. (compare with Figure S1 in the Supplementary Material) and, as expected, the Dirac cone seems to be entirely due to C $p_z$ orbitals. However, zooming on the region near the Dirac point (see Figure 3b), and subtracting the C orbitals contribution, not only a small band gap opening and a parabolic dispersion become evident, but also a very weak contribution of B $p_z$ and, to a minor extent, of N $p_z$ orbitals, are unveiled. In other words, the band structure of the C-BN-BN-C near K is not just the sum of those of the starting monolayers, but there is a (weak) interaction between the layers, which results in a very small hybridization of C with B and N $p_z$ orbitals. This interaction makes the C atoms non-equivalent and causes the opening of the small gap of 51 meV at the Dirac point (visible in Figure 2a, and in Figure 3b). [53] In order to complete our analysis, we report in Figure 4 the band structure and the corresponding PDOS of a C-BN bilayer with AB stacking. As for the AB-AA'-AB HT case, a parabolic band dispersion is found with a very similar band gap opening of ~49 meV.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Electronic band structure and Projected Density of States (PDOS) of CBN bilayer in AB stacking in 1 × 1 unit cell. Inset: zoom at the Dirac point showing the opening of a small band gap.</figcaption></figure>

By comparison with Figure 3, we conclude that, in the AB-AA'-AB HT, the presence of one (C-BN) or two (C-BN-BN-C) BN layers mainly influences the electronic properties far from the Fermi energy, while the low energy band dispersion and gap opening can be already explained by the interaction of one graphene with one single BN layer.

A simple explanation of the low energy bands dispersion of all the quadrilayers, reported in Figure 2, can be extracted using first-neighbor 2 × 2 and 4 × 4 tight-binding (TB) models (see Supporting Information for a complete discussion). The 2 × 2 TB hamiltonian describes the C-BN bilayer, modeling it as a single layer of graphene with two non-equivalent carbon sites; the 4 × 4 TB matrix describes the quadrilayer as two layers of interacting graphene, both of them made up of two non-equivalent carbon sites. The quadrilayer TB solutions are schematically depicted in Figure 5. In the AB-AA'-AB HT case (Figure 5a), the parabolic band dispersion and gap opening are practically identical to those found in each isolated bilayer (in this case, both bilayers having AB stacking). The interaction between the layers has only a minor role and does not change the qualitative behavior. In the AB-AA'-AB' HT, where two non-equivalent bilayers (AB and AB') are present, the corresponding parabolic bands are shifted down (AB) and up (AB'). When the bilayers get closer, the parabola interpenetrate due to the build up of higher total dipole moment at the interface. Moreover, due to the small but not negligible interaction between the bilayers, two hybridization gaps open at the left and at the right of K, at the degenerate crossing points (see Figure 5b). Finally, in the twisted configuration, the situation is intermediate: while in AB (AB') C-BN bilayer half of C atoms has a B(N) as vertical first neighbor, the 21.8° rotation angle introduces several different vertical stackings, globally reducing the interaction with BN (see Figure 5c) and recovering a linear band dispersion behavior very close to the Dirac points.

In order to capture from ab-initio calculations, how the electronic bandstructure of AB-AA'-AB' HT quadrilayer evolves from very distant AB and AB' separated bilayers to the real composite, we perform several calculations decreasing the distance (from very large to the real equilibrium distance) between the two bilayers C-BN, using a supercell with very large size along z. Figure S5 of Supplementary Materials shows how, by decreasing the distance between the bilayers, the two parabolas start to compenminate, because of the intrinsic dipoles, present in both the separated bilayers, sum up. In other words, the net effect of decreasing the distance between the AB and AB' bilayers is to increase the
intrinsic electric field (which causes a vertical shift of the parabolas) and to switch on the inter-bilayer interaction (which causes the opening of the small lateral gaps).

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Schematic representation of the contribution of the C-BN bilayers energy states to the band structure of C-BN-BN-C quadrilayers. In the figure, we show both cases with and without the interlayer interaction for the different configurations: AB-AA'-AB (**a**), AB-AA'-AB' (**b**) and twisted (**c**). C(B) script highlights the fact that in the AB stacking a B atom is on top of a C atoms, whereas C(N) indicates that in AB' stacking a N atom is on top of a C atom.</figcaption></figure>

## Doping and Electric Field

To complete our study, we discuss, for the AB-AA'-AB and T-AA'-T HTs, how vertical gate voltage and doping of the external graphene layers can modulate their electronic properties.
Three different electric fields, with strengths of 0.003, 0.005, and 0.007 a.u. (corresponding to 0.15, 0.26, and 0.36 V/Å) perpendicular to the *xy* plane of graphene and BN layers, have been chosen and applied to AB-AA'-AB and T-AA'-T HTs, while the p-n doping has been implemented using ad-hoc modified carbon pseudopotentials, as discussed in Section 2. The band structures near *K* (black lines), compared with the corresponding ones obtained without an electric field (red lines), are shown in Figure 6. In the presence of an external field, the Dirac cones split (the degeneracy is lifted) and interpenetrate, displacing new small band gaps on either side close to the *K* point of 1BZ. The amount of the cones interpenetration is proportional to the magnitude of the electric field (see Table 4).

<table><caption>**Table 4.** Electronic gaps at K for the quadrilayers. In addition, the values of the small lateral gaps $E_{\text{gap}}^{\text{lat}}$ are reported. Electric field strengths are in a.u.</caption><thead><tr><th>C-BN-BN-C</th><th>E<sup>(K)</sup><sub>gap</sub></th><th>E<sup>lat</sup><sub>gap</sub></th></tr><tr><th>Stacking</th><th>[eV]</th><th>[meV]</th></tr></thead><tbody><tr><td>AB-AA'-AB</td><td></td><td></td></tr><tr><td>Pristine</td><td>0.051</td><td>—</td></tr><tr><td>E&#773; 0.003</td><td>0.430</td><td>8.0</td></tr><tr><td>E&#773; 0.005</td><td>0.639</td><td>3.8</td></tr><tr><td>E&#773; 0.007</td><td>0.819</td><td>5.8</td></tr><tr><td>doped</td><td>1.245</td><td>7.3</td></tr><tr><td>T-AA'-T</td><td></td><td></td></tr><tr><td>pristine</td><td>0.0073</td><td>0.4</td></tr><tr><td>E&#773; 0.003</td><td>0.406</td><td>0.1</td></tr><tr><td>E&#773; 0.005</td><td>0.599</td><td>2.0</td></tr><tr><td>E&#773; 0.007</td><td>0.681</td><td>1.8</td></tr><tr><td>doped</td><td>0.261</td><td>3.1</td></tr><tr><td>AB-AA'-AB'</td><td>0.126</td><td>6.7</td></tr></tbody></table>

HTs as the ones studied in the present manuscript are usually deposited on a substrate; therefore, it often happens that one of the two interfaces is doped with respect to the other one due to impurities present in the environment. For this reason, we investigate the effect of doping on one of the two carbon layers of the AB-AA'-AB HT. We found that doping induces a strong interpenetration of the Dirac cones. In fact, in this system, the doping is quite large (all C atoms in the top layer have been substituted with $C_n$ pseudoatoms, and all C atoms in the bottom layer have been substituted with $C_p$ atoms.) By exploiting the linear trend between the gap at K and the applied field, we can estimate that p-n-doping has, on the AB-AA'-AB HT, a net effect similar to an applied electric field of magnitude 0.01 a.u.

The case of the twisted HT is more complicated. Although the gap at K increases with the strength of the field, no clear linear trend is observed. Anyway, a rough linear fit gives an electric field of about 0.002 a.u. This field is smaller than the value obtained for the co-doped AB-AA'-AB HT, in agreement with the fact that, in the twisted quadrilayer, less than one third of the C atoms have been substituted with the pseudo-atoms $C_p$ and $C_n$. A further rough estimation, performed considering the four-layers system as a planar capacitor, confirms this order of magnitude for the electric field.

Finally, we notice that in the twisted case the Fermi level shifts downwards under the effect of the applied external electric field, and upon doping. This shift is due to the presence of an empty band (visible in the upper right corners in Figure 6b), which is very sensitive to perturbations, and becomes occupied for increasing electric fields. For what concerns the small lateral gaps on the left and on the right of the K point, we report in Table 4 their values. Their amplitudes do not show a linear trend with the applied electric field. They are in the range of few meV and below, hence it is interesting for possible generation and detection of THz radiation.
<figure><img src="image_10.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_11.png" /></figure>

**Figure 6.** Black lines: The effect of three different external electric fields in a.u. and p-n doping on the band structure of C-BN-BN-C in (a) AB-AA'-AB configuration and (b) twisted configuration. With red lines, we report results without electric field as a reference. The zero of the energy is set at the Fermi level. The scale on the x-axis in (a,b) is not the same, hence the dissimilar slope of the bands does not indicate variations in the Fermi velocity of the carriers (see Table 3).
## 4. Conclusions

We have performed a comprehensive study of the structural and electronic properties of C-BN-BN-C quadrilayers, for three different configurations obtained by different stacking and by twisting. The most stable structure results in being the AB-AA'-AB, but all possess a negative stacking energy, suggesting that all of them can be experimentally obtained. Important differences arise in the electronic band structures near the Fermi energy, which can be traced back to the different first neighbor interaction of the non-equivalent carbon atoms in each graphene plane. Moreover, we have investigated the effect of an external electric field on these HTs. Our results show that it is possible to finely tune the band gap of a graphene/BN multilayered heterostructure by applying an external electric field perpendicular to the atomic planes. Finally, we investigate the presence of an asymmetric p-n doping on single carbon layer, which is expected in experimental realization of these HTs due to the presence of impurities, induced by disorder or edges [54]. We found that asymmetric doping acts as an effective electric field on these HTs, and this means that its effect can be compensated by an opposite electric field to restore the original band structure of the HT.

Our work suggests a possible use of graphene/BN multilayers as a detector of THz radiation, by playing with stacking order, number of BN layers, twist angle, external gate voltage, and selective doping.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/nano12122118/s1, Figure S1: Band structure of graphene monolayer and 2D-BN bilayer AA' in 1 × 1 unit cell.; Figure S2: Band structure of C-BN-BN-C AB-AA'-AB in a $\sqrt{7} \times \sqrt{7}$ supercell; Figure S3: Plot of the homo-1, homo, lumo, lumo+1 level at K for CBNBNC in a) AB-AA'-AB and b) twisted case; Figure S4: Electronic band structure and projected density of states of CBN bilayer in AB stacking (lower panel) and in twisted configuration (upper panel); Figure S5: Band structure calculations performed for C-BN-BN-C in AB-AA'-AB' configuration varying the distance between the two C-BN AB and AB' bilayers; Figure S6: Plot of the two energy eigenvalues $E_{\pm}$ of the 2 × 2 matrix for the C-BN bilayers; Figure S7: Plot of the four energy eigenvalue of the 4 × 4 matrix for AB-AA'-AB stacking evaluated for different values of $\Delta\epsilon(a)$, $V(b)$ and $t_0(c)$ with $\Delta\epsilon^{(2)} > \Delta\epsilon^{(1)} > 0$, $V^{(2)} > V^{(1)} > 0$ and $t_0^{(2)} > t_0^{(1)} > 0$.

**Author Contributions:** Conceptualization, O.P.; methodology, M.P. and F.B.; validation, C.A. and O.P.; formal analysis, S.B.; investigation, S.B. and O.P.; resources, O.P.; data curation, O.P.; writing—original draft preparation, O.P., F.B. and G.G.; writing—review and editing, G.G., C.A. and O.P.; visualization, S.B.; supervision, O.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** O.P., S.B. and M.P. acknowledge financial funding from the EU MSCA-RISE project DiSeT-Com (HORIZON2020 GA 823728) and INFN project TIME2QUEST. M.P. acknowledges Tor Vergata University Project 2021 TESLA. F.B. acknowledges financial funding from the European Union's Horizon 2020 research and innovation programme under Grant No. 824158 (EoCoE-II). The research leading to these results has received funding from the European Union Seventh Framework Program under Grant No. 785219 Graphene Core2. This publication is based upon work from COST Action TUMIEE CA17126, supported by COST (European Cooperation in Science and Technology).

**Data Availability Statement:** Not applicable.

**Acknowledgments:** CPU time was granted by CINECA (ISCRA-B and ISCRA-C) and CRESCO ENEA HPC centers.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Geim, A.K.; Novoselov, K.S. The rise of graphene. *Nat. Mater.* **2007**, *6*, 183–191. [CrossRef] [PubMed]

2. Castro Neto, A.H.; Guinea, F.; Peres, N.M.R.; Novoselov, K.S.; Geim, A.K. The electronic properties of graphene. *Rev. Mod. Phys.* **2009**, *81*, 109–162. [CrossRef]

3. Novoselov, K.S.; Geim, A.K.; Morozov, S.V.; Jiang, D.; Zhang, Y.; Dubonos, S.V.; Grigorieva, I.V.; Firsov, A.A. Electric field effect in atomically thin carbon films. *Science* **2004**, *306*, 666–669. [CrossRef]
4. Novoselov, K.S.; Neto, A.H.C. Two-dimensional crystals-based heterostructures: materials with tailored properties. *Phys. Scr.* **2012**, *T146*, 014006. [CrossRef]

5. Bhimanapati, G.R.; Lin, Z.; Meunier, V.; Jung, Y.; Cha, J.; Das, S.; Xiao, D.; Son, Y.; Strano, M.S.; Cooper, V.R.; et al. Recent Advances in Two-Dimensional Materials beyond Graphene. *ACS Nano* **2015**, *9*, 11509–11539. [CrossRef] [PubMed]

6. Bernardi, M.; Ataca, C.; Palummo, M.; Grossman, J. Electronic Properties of Two-Dimensional Layered Materials. *Nanophotonics* **2016**, *6*, 479–493. [CrossRef]

7. Das, S.; Robinson, J.A.; Dubey, M.; Terrones, H.; Terrones, M. Beyond Graphene: Progress in Novel Two-Dimensional Materials and van der Waals Solids. *Annu. Rev. Mater. Res.* **2015**, *45*, 1–27. [CrossRef]

8. Chiodo, L. Two-dimensional innovative materials for photovoltaics. *Curr. Opin. Green Sustain. Chem.* **2019**, *17*, 49–56. [CrossRef]

9. Bartolomei, M.; Carmona-Novillo, E.; Hernández, M.I.; Campos-Martínez, J.; Pirani, F.; Giorgi, G. Graphdiyne pores: "Ad Hoc" openings for helium separation applications. *J. Phys. Chem.* **2014**, *118*, 29966–29972. [CrossRef]

10. Bartolomei, M.; Carmona-Novillo, E.; Hernández, M.I.; Campos-Martínez, J.; Pirani, F.; Giorgi, G.; Yamashita, K. Penetration barrier of water through graphynes' pores: first-principles predictions and force field optimization. *J. Phys. Chem. Lett.* **2014**, *5*, 751–755. [CrossRef]

11. Bartolomei, M.; Carmona-Novillo, E.; Giorgi, G. First principles investigation of hydrogen physical adsorption on graphynes' layers. *Carbon* **2015**, *95*, 1076–1081. [CrossRef]

12. Bartolomei, M.; Giorgi, G. A novel nanoporous graphite based on graphynes: first-principles structure and carbon dioxide preferential physisorption. *ACS Appl. Mater. Interfaces* **2016**, *8*, 27996–28003. [CrossRef] [PubMed]

13. Palummo, M.; Yamashita, K.; Giorgi, G. Advances in two-dimensional green materials for organic electronics applications. In *Sustainable Strategies in Organic Electronics*; Elsevier: Amsterdam, The Netherlands, 2022; pp. 391–422.

14. Bechstedt, F.; Gori, P.; Pulci, O. Beyond graphene: Clean, hydrogenated and halogenated silicene, germanene, stanene, and plumbene. *Prog. Surf. Sci.* **2021**, *96*, 100615. [CrossRef]

15. Geim, A.K.; Grigorieva, I.V. Van der Waals heterostructures. *Nature* **2013**, *499*, 419–425. [CrossRef] [PubMed]

16. Novoselov, K.S.; Jiang, D.; Schedin, F.; Booth, T.; Khotkevich, V.; Morozov, S.; Geim, A.K. Two-dimensional atomic crystals. *Proc. Natl. Acad. Sci. USA* **2005**, *102*, 10451–10453. [CrossRef] [PubMed]

17. Gao, G.; Gao, W.; Cannuccia, E.; Taha-Tijerina, J.; Balicas, L.; Mathkar, A.; Narayanan, T.; Liu, Z.; Gupta, B.K.; Peng, J.; et al. Artificially stacked atomic layers: toward new Van der Waals solids. *Nano Lett.* **2012**, *12*, 3518–3525. [CrossRef]

18. Morozov, S.; Novoselov, K.; Katsnelson, M.; Schedin, F.; Elias, D.; Jaszczak, J.A.; Geim, A. Giant intrinsic carrier mobilities in graphene and its bilayer. *Phys. Rev. Lett.* **2008**, *100*, 016602. [CrossRef]

19. Jian-Hao, C.; Jang, C.; Xiao, S.; Masa, I.; Fuhrer, M. Intrinsic and extrinsic performance limits of graphene devices on SiO₂. *Nat. Nanotechnol.* **2008**, *3*, 206–209.

20. Dean, C.R.; Young, A.F.; Meric, I.; Lee, C.; Wang, L.; Sorgenfrei, S.; Watanabe, K.; Taniguchi, T.; Kim, P.; Shepard, K.L.; et al. Boron nitride substrates for high-quality graphene electronics. *Nat. Nanotechnol.* **2010**, *5*, 722–726. [CrossRef]

21. Coraux, J.; Plasa, T.N.; Busse, C.; Michely, T.; N'Diaye, A.T. Structure of epitaxial graphene on Ir (111). *New J. Phys.* **2008**, *10*, 043033. [CrossRef]

22. Levinshtein, M.E.; Rumyantsev, S.L.; Shur, M.S. *Properties of Advanced Semiconductor Materials: GaN, AlN, InN, BN, SiC, SiGe*; John Wiley & Sons: Hoboken, NJ, USA, 2001.

23. Ju, L.; Velasco, J.; Huang, E.; Kahn, S.; Nosiglia, C.; Tsai, H.Z.; Yang, W.; Taniguchi, T.; Watanabe, K.; Zhang, Y.; et al. Photoinduced doping in heterostructures of graphene and boron nitride. *Nat. Nanotechnol.* **2014**, *9*, 348–352. [CrossRef] [PubMed]

24. Zomer, P.; Guimarães, M.; Brant, J.; Tombros, N.; Van Wees, B. Fast pick up technique for high quality heterostructures of bilayer graphene and hexagonal boron nitride. *Appl. Phys. Lett.* **2014**, *105*, 013101. [CrossRef]

25. Mishchenko, A.; Tu, J.; Cao, Y.; Gorbachev, R.V.; Wallbank, J.; Greenaway, M.; Morozov, V.; Morozov, S.; Zhu, M.; Wong, S.; et al. Twist-controlled resonant tunnelling in graphene/boron nitride/graphene heterostructures. *Nat. Nanotechnol.* **2014**, *9*, 808–813. [CrossRef] [PubMed]

26. Ni, G.; Wang, H.; Wu, J.; Fei, Z.; Goldflam, M.; Keilmann, F.; Özyilmaz, B.; Castro Neto, A.; Xie, X.; Fogler, M.; et al. Plasmons in graphene moiré superlattices. *Nat. Mater.* **2015**, *14*, 1217–1222. [CrossRef] [PubMed]

27. Woessner, A.; Lundeberg, M.B.; Gao, Y.; Principi, A.; Alonso-González, P.; Carrega, M.; Watanabe, K.; Taniguchi, T.; Vignale, G.; Polini, M.; et al. Highly confined low-loss plasmons in graphene–boron nitride heterostructures. *Nat. Mater.* **2015**, *14*, 421–425. [CrossRef]

28. Li, Q.; Liu, M.; Zhang, Y.; Liu, Z. Hexagonal boron nitride–graphene heterostructures: Synthesis and interfacial properties. *Small* **2016**, *12*, 32–50. [CrossRef]

29. Pan, W.; Xiao, J.; Zhu, J.; Yu, C.; Zhang, G.; Ni, Z.; Watanabe, K.; Taniguchi, T.; Shi, Y.; Wang, X. Biaxial compressive strain engineering in graphene/boron nitride heterostructures. *Sci. Rep.* **2012**, *2*, 893. [CrossRef]

30. Wickramaratne, D.; Weston, L.; Van de Walle, C.G. Monolayer to bulk properties of hexagonal boron nitride. *J. Phys. Chem.* **2018**, *122*, 25524–25529. [CrossRef]

31. Arnaud, B.; Lebègue, S.; Rabiller, P.; Alouani, M. Huge excitonic effects in layered hexagonal boron nitride. *Phys. Rev. Lett.* **2006**, *96*, 026402. [CrossRef]

32. Schué, L.; Sponza, L.; Plaud, A.; Bensalah, H.; Watanabe, K.; Taniguchi, T.; Ducastelle, F.M.C.; Loiseau, A.; Barjon, J. Bright Luminescence from Indirect and Strongly Bound Excitons in h-BN. *Phys. Rev. Lett.* **2019**, *122*, 067401. [CrossRef]
33. Pikalov, A.A.; Fil, D.V. Graphene bilayer structures with superfluid magnetoexcitons. *Nanoscale Res. Lett.* **2012**, *7*, 145. [CrossRef] [PubMed]

34. Ponomarenko, L.A.; Geim, A.K.; Zhukov, A.A.; Jalil, R.; Morozov, S.V.; Novoselov, K.S.; Grigorieva, I.V.; Hill, E.H.; Cheianov, V.V.; Fal'ko, V.I.; et al. Tunable metal—Insulator transition in double-layer graphene heterostructures. *Nat. Phys.* **2011**, *7*, 958–961. [CrossRef]

35. Profumo, R.E.V.; Asgari, R.; Polini, M.; MacDonald, A.H. Double-layer graphene and topological insulator thin-film plasmons. *Phys. Rev. B* **2012**, *85*, 085443. [CrossRef]

36. Farmani, A.; Yavarian, M.; Alighanbari, A.; Miri, M.; Sheikhi, M.H. Tunable graphene plasmonic Y-branch switch in the terahertz region using hexagonal boron nitride with electric and magnetic biasing. *Appl. Opt.* **2017**, *56*, 8931–8940. [CrossRef]

37. Yadav, D.; Tombet, S.B.; Watanabe, T.; Arnold, S.; Ryzhii, V.; Otsuji, T. Terahertz wave generation and detection in double-graphene layered Van der Waals heterostructures. *2D Mater.* **2016**, *3*, 045009. [CrossRef]

38. Viti, L.; Purdie, D.G.; Lombardo, A.; Ferrari, A.C.; Vitiello, M.S. HBN-encapsulated, graphene-based, room-temperature terahertz receivers, with high speed and low noise. *Nano Lett.* **2020**, *20*, 3169–3177. [CrossRef]

39. Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G.L.; Cococcioni, M.; Dabo, I.; et al. Quantum Espresso: A modular and open-source software project for quantum simulations of materials. *J. Phys. Condens. Matter* **2009**, *21*, 395502. [CrossRef]

40. Giannozzi, P.; Andreussi, O.; Brumme, T.; Bunau, O.; Nardelli, M.B.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M.; et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. *J. Phys. Condens. Matter* **2017**, *29*, 465901. [CrossRef]

41. Giannozzi, P.; Baseggio, O.; Bonfà, P.; Brunato, D.; Car, R.; Carnimeo, I.; Cavazzoni, C.; De Gironcoli, S.; Delugas, P.; Ferrari Ruffino, F.; et al. Quantum ESPRESSO toward the exascale. *J. Chem. Phys.* **2020**, *152*, 154105. [CrossRef]

42. Fuchs, M.; Scheffler, M. Ab initio pseudopotentials for electronic structure calculations of poly-atomic systems using density-functional theory. *Comput. Phys. Commun.* **1999**, *119*, 67–98. [CrossRef]

43. Hamann, D. Generalized norm-conserving pseudopotentials. *Phys. Rev. B* **1989**, *40*, 2980. [CrossRef] [PubMed]

44. Troullier, N.; Martins, J.L. Efficient pseudopotentials for plane-wave calculations. *Phys. Rev. B* **1991**, *43*, 1993. [CrossRef] [PubMed]

45. Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. *Phys. Rev. Lett.* **1996**, *77*, 3865. [CrossRef]

46. Hamada, I. Van der Waals density functional made accurate. *Phys. Rev. B* **2014**, *89*, 121103. [CrossRef]

47. Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D.C.; Lundqvist, B.I. Van der Waals density functional for general geometries. *Phys. Rev. Lett.* **2004**, *92*, 246401. [CrossRef] [PubMed]

48. Liao, J.; Sa, B.; Zhou, J.; Ahuja, R.; Sun, Z. Design of high-efficiency visible-light photocatalysts for water splitting: $MoS_2$/AlN (GaN) heterostructures. *J. Phys. Chem.* **2014**, *118*, 17594–17599. [CrossRef]

49. Zhang, Z.; Huang, B.; Qian, Q.; Gao, Z.; Tang, X.; Li, B. Strain-tunable III-nitride/ZnO heterostructures for photocatalytic water-splitting: A hybrid functional calculation. *APL Mater.* **2020**, *8*, 041114. [CrossRef]

50. Guo, H.; Zhang, Z.; Huang, B.; Wang, X.; Niu, H.; Guo, Y.; Li, B.; Zheng, R.; Wu, H. Theoretical study on the photocatalytic properties of 2D InX (X= S, Se)/transition metal disulfide ($MoS_2$ and $WS_2$) Van der Waals heterostructures. *Nanoscale* **2020**, *12*, 20025–20032. [CrossRef]

51. Guo, H.; Yin, Y.; Niu, H.; Robertson, J.; Zhang, Z.; Guo, Y. Tunable contacts and device performances in graphene/group-III monochalcogenides MX (M= In, Ga; X= S, Se) Van der Waals heterostructures. *J. Appl. Phys.* **2021**, *130*, 144303. [CrossRef]

52. Jiang, H.L.; Pan, J.; Zhou, W.; Li, H.M.; Liu, S. Fabrication and application of arrays related to two-dimensional materials. *Rare Met.* **2021**, *41*, 262–286. [CrossRef]

53. Giovannetti, G.; Khomyakov, P.A.; Brocks, G.; Kelly, P.J.; Van Den Brink, J. Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations. *Phys. Rev. B* **2007**, *76*, 073103. [CrossRef]

54. Siegel, D.A.; Zhou, S.Y.; El Gabaly, F.; Fedorov, A.V.; Schmid, A.K.; Lanzara, A. Self-doping effects in epitaxially grown graphene. *Appl. Phys. Lett.* **2008**, *93*, 243119. [CrossRef]