Graphene, a single layer of sp²- bonded carbon atoms having a 2D honeycomb lattice, is the primary element for all other graphitic forms like fullerenes and carbon nanotubes (CNTs) etc.^[1]. Si-based technology is approaching a physical limitation as a result of extensive continuous scaling due to dominant short channel effects (SCEs) and hot carrier effects (HCEs)^[2]. Thus graphene based devices are expected to be the substitute owing to high mobility, excellent thermal conductivity, prodigious mechanical attributes and stability^[3]. Despite its astounding properties, the lack of bandgap in its pristine form restrains its actual application in semiconductor electronics. This has led to extensive research for methods to open and tune the bandgap in graphene. Trimming graphene into nanoribbons^[4], deposition of graphene on top of a substrate^[5], doping with other elements^[6], hydrogenation of graphene^[7] and applying strain^[8], are examples of such methods.

The application of external perpendicular electric field is the most efficacious method to open and tune bandgap because it tears inversion symmetry without carrier mobility being affected notably in bi-layer graphene (BLG)^[9–11]. Because of this, the application of an external perpendicular electric field to BLG and graphene heterostructures, is the method used in this work to opening the bandgap.

Graphitic carbon nitride (g-C₃N₄), appraised to be the most stable allotropes of C₃N₄ under ambient conditions, has enticed substantial considerations because of its propitious application in electronic devices^[12]. The bandgap opening and tuning of graphene/g-C₃N₄ hetero bi-layer (HBL), through the application of an external electric field, has been done previously^[13]. g-C₃N₄/graphene/g-C₃N₄ sandwich heterostructure under applied external electric field has been also studied^[14]. The band structures under an external electric field of the two heterostructures: (i) graphitic carbon nitride-graphene-hexagonal Boron Nitride (g-C₃N₄/G/h-BN) in hexagonal stack (AAA)^[15], and (ii) graphene-graphitic carbon nitride-hexagonal boron nitride (G/g-C₃N₄/h-BN) in Bernal stack (ABA), as shown by Fig. 1, has been proposed and analyzed. The interlayer distance, effective mass and binding energy has also been calculated for the proposed heterostructures and comparison has been made with BLG, graphene-hexagonal boron nitride (G/h-BN) hetero bi-layer, graphene-graphitic carbon Nitride (G/g-C₃N₄) hetero bi-layer and graphitic carbon nitride-graphene- graphitic carbon nitride (g-C₃N₄/G/g-C₃N₄) heterostructures in Bernal and hexagonal stack.

All the calculations have been performed on Quantumwise Atomistix Toolkit (ATK) simulation package^[16]. The structures have been optimized using density functional theory (DFT) with the generalized gradient approximation exchange correlation in the parameterization of Perdew-Burke-Ernzerhof (GGA-PBE). The electronic structures are calculated using an LCAO calculator which uses DFT and norm-conserving pseudopotentials by expansion of single-particle wave functions in a basis of numerical atomic orbitals with compact support^[17–19]. When the external perpendicular electric field is applied the charge density is reorganized which is also dependent on number of layers. This field causes a charge density difference which results in bandgap opening. In DFT, the electron density is given by filled Eigen states of Kohn-Sham Hamiltonian: , where represents the occupancy of energy level α given by the Fermi-Dirac distribution for finite temperature:

$ {f_\alpha } = {\left[ {1 + \exp \left( {\frac{{{\varepsilon _\alpha } - {\varepsilon _ {\rm{F}}}}}{{kT}}} \right)} \right]^{ - 1}} . $  (1)

The electron difference density, a comparison between the electron density of a many-body system to the superposition of individual atom-based electron density is given by: , where is the position of atom μ in the many-body system; The effective potential is: , where is the Hartree potential due to mean field interaction between electrons, is the exchange correlation potential arising from quantum mechanical nature of electrons, and potential is due to other field in the system comprised of electrostatic potential of ions (due to norm-conserving potentials) and external electrostatic fields (due to external sources). The total energy is: , where is kinetic energy of non-interacting electron gas having density n, denotes the exchange correlation energy, denotes Hartree potential energy and denotes the interaction energy of electrons.

Band structure and bandgap calculation of optimized structures have been done using ATK-SE (semi-empirical) Extended Huckel method^[20]. In ATK-SE, total energy is given by: , where is one-electron energy of the non-self-consistent Hamiltonian, is electrostatic difference energy, is electrostatic interaction between electrons and external field, is the spin polarization energy and is repulsive energy from a pair-potential between each atom pair.

To calculate the distance between two graphene layers, van der Waals interactions are included. The long-range van der Waals interaction is included through the semi-empirical Grimme correction (DFT-D2) by adding a term to the DFT total energy and given by: . An attractive semi-empirical pair potential gives dispersion correction :

$ {E_ {\rm{disp}}} = {S_6}\sum\limits_{\mu < \;\mu '} {{V^{\rm PP}}} \left( {{Z_\mu },{Z_{\mu '}},{R_{\mu ,\mu '}}} \right) , $  (2)

where the value of is 0.75 for GGA-PBE functional.

The density mesh cutoff of 75 Hartree and tolerance of 10⁻⁵ Hartree has been used with the Pulay mixer algorithm^[21] with a maximum of 200 steps and double zeta polarized basis set. The damped van der Waals (vdW) correction, proposed by Grimme (PBE-D2), has been adopted as weak interactions are not described in standard PBE function well, whereas weak van der Waals interactions are considered to be critical. Scale factor, damping factor and cutoff distance (Å) are set to 0.75, 20 and 30, respectively, in the parameters used for Grimme DFT-D2. Force tolerance of 0.02 eV/Å has been used for geometry optimization. Band structure and bandgap calculation of optimized structures have been done using ATK SE (semi-empirical) Extended Huckel method with density mesh cutoff 10 Hartree, unpolarized spin and Wolfsberg weighting scheme. Cerda Carbon (graphite), Cerda Boron (BN hexagonal) and Cerda Nitrogen has been used as basis type under Huckel basis sets. Hamiltonian variable is used as mixing variable. The sampling of 9 × 9 × 1 K-points is used for structure optimization and 11 × 11 × 1 K-points is used for band structure calculation. A multi-grid Poisson solver has been used with boundary condition set to Dirichlet on both sides in the C direction. The above parameters are used for optimization as well as inter-layer distance and binding energy calculation.

The hexagonal boron nitride h-BN, sometimes referred to as ‘white graphene’, has lattice structure similar to grapheme. The lattice mismatch is about of 1.5% only^[22]. The g-C₃N₄ lattice parameter is almost three times the lattice parameter of grapheme. Hence commensurability is imposed for G/g-C₃N₄ heterostructure. Similarly commensurability is imposed for hetero tri-layers to match the lattice parameters. The strain introduced in the process does not affect the conclusions^{[13, 14]}.

There is no bandgap in the BLG in Bernal (AB) stack at 0 V/nm as shown by Fig. 2(a). It can be seen from the figure that there are two almost parallel conduction bands beyond the two almost parallel valence bands near the Fermi level of the Bernal stacked BLG without any gating i.e. absence of an external electric field. The conduction band minimum (CBM) and the valence band maximum (VBM) touch each other around the K-point, if no electric field is applied, resulting in zero bandgap. An external perpendicular electric field is applied by means of two metallic plates (one at the top E_t and one at the bottom E_b to induce the bandgap in BLG. As the field is applied, two effects are produced: (i) i.e. difference of the two fields results in a net carrier doping and i.e. the average of the two fields breaks the inversion symmetry^[9]. The CBM shifts upwards while VBM shifts downwards causing bandgap in the BLG as shown by Fig. 2(b). This holds true for all the structures i.e. bi-layer as well as tri-layer graphene structures.

The interlayer distance of various optimized graphene structures are shown in Tables 1 and 2. The vectors d₁ and d₂ represent distance between neighboring in-plane carbon atoms. The ‘–’ in tables represent not-applicable cases. It has been found that interlayer distance of various graphene structures decrease with increasing electric field in comparison to zero electric field. It implies that due to the application of an electric field, the structures get slightly distorted. Bandgap comparison for structures in Bernal stack and hexagonal stack are shown in Figs. 3 and 4.

It is observed from these figures that bandgap is dependent on interlayer spacing as well as stacking pattern. It can be seen from Figs. 3(b) and 4(b) that as the interlayer spacing approaches 4.0 Å, the bandgap also approaches zero for all the structures. This is because of such large interlayer spacing’s; there is hardly any interaction between the layers. As smaller interlayer spacing is approached (towards 2 Å), the bandgap is found to increase with decrease in interlayer spacing.

The layers having these charges interact with each other and build a field themselves, depending on the interlayer spacing. The π-orbitals overlap of neighboring layers increases with high charge density. This overlap is greater in hexagonal stack structures because of its geometry and less in Bernal stack structures. Hence the bandgap is greater in hexagonal stack structures as shown in Figs. 3(b) and 4(b).

The structural stability of the heterostructures, is assessed by the binding energy (E_b). Binding energies for different structures have been calculated using formulae given in Table 3. ‘n’ is the number of C atoms in the graphene, E_C is the total energy of carbon atoms. E_{heterostructure} is the energy of the heterostructure, E_G is the energy of graphene, E_BN is the energy of hexagonal boron nitride and E_C3N4 is the energy of graphitic Carbon Nitride with lattice parameter being same. The negative values of binding energy indicate the stable heterostructures. The heterostructures are more favorable if the binding energy is more negative i.e. higher |E_b| value. The C₃N₄/G/C₃N₄-AAA has the highest stability among the various structures studied as represented by Fig. 5. g-C₃N₄/G/h-BN-AAA is the best choice, since the effective mass is lowest of all the structures (Table 4) implying that it exhibits high carrier mobility, while the band gap is significantly opened. The results of the present work indicate that the properties of graphene heterostructures are affected by stacking pattern and the applied electric fields.

The effective mass is closely related to the carrier mobility. The effective mass (of electrons and holes) is calculated as where k denotes wave vector and E(k) is the corresponding energy which is obtained by fitting the parabolic functions to the valence band maximum and conduction band minimum^[13]. It has been found that the electron and hole effective masses for C₃N₄/G/BN heterostructure under applied external fields are lower than other structures, from Table 4.

Sandwich heterostructures C₃N₄/G/BN-AAA and G/C₃N₄/BN-ABA have been proposed and investigated and their properties have been compared using DFT with other reported structures. C₃N₄/G/BN-AAA turns out to be the best choice among all the structures studied because it offers lowest effective mass i.e. higher carrier mobility with significant bandgap opening. The results suggest that effectively controlling the stacking patterns and applying a perpendicular electric field significantly affect the properties of the heterostructures and open up exciting opportunities for the development of electronic and optoelectronic devices.