Stabilities and novel electronic structures of three carbon nitride bilayers.

We predict three novel phases of the carbon nitride (CN) bilayer, denoted α-C2N2, β-C2N2 and γ-C4N4, respectively. All of them consist of two CN sheets connected by C-C covalent bonds. The phonon dispersions reveal that all these phases are dynamically stable, because no imaginary frequency is present. The transition pathway between α-C2N2 and β-C2N2 is investigated, which involves bond-breaking and bond-reforming between C and N. This conversion is difficult, since the activation energy barrier is 1.90 eV per unit cell, high enough to prevent the transformation at room temperature. Electronic structure calculations show that all three phases are semiconductors with indirect band gaps of 3.76/5.22 eV, 4.23/5.75 eV and 2.06/3.53 eV, respectively, by PBE/HSE calculation. The β-C2N2 has the widest band gap among the three phases. All three bilayers can become metallic under tensile strain, and the indirect gap of γ-C4N4 can turn into a direct one. γ-C4N4 can become an anisotropic Dirac semimetal under uniaxial tensile strain. Anisotropic Dirac cones with high Fermi velocity of the order of 105 m/s appear under 12% strain. Our results suggest that the three two-dimensional materials have potential applications in electronics, semiconductors, optics and spintronics.

Introduction

Two-dimensional (2D) materials demonstrate many novel electronic and magnetic properties such as high mobility and optical characteristics[1,2]. The low-buckled honeycomb lattice that consists of silicon and germanium atoms, which are called silicene and germanene, have been predicted[3], and germanene exhibits the quantum-spin Hall effect[4] and the materials were synthetized soon after their prediction[5,6]. Plumbene is a normal insulator in its free state, and it can turn into topological insulator by electron doping[7].

The 2D honeycomb monolayer, consisting of nitrogen atoms denoted nitrogene, has been proposed and its electronic properties have been deeply investigated[8,9]. The electronic properties of nitrogene with vacancy and adsorbed adatoms were also analyzed[10]. Around the same time, another allotrope, octagon-nitrogene, was proposed, and its stability and electronic structure have been systematically studied[11,12].

The number of 2D materials has been restricted due to the limited number of possible geometric structures, and new compound materials that contain more than one elements are increasing gaining attention. With the rapid development of new low-dimensional materials and devices in experiments and applications, the exploration of new 2D materials is in drastic need. Binary compounds based on two types of atoms may possess new phenomenon compared to their element counterparts[13], for example, the fullerene-like carbon nitride has been investigated deeply by both experiments[14] and first-principles calculations[15,16], and the bulk structures formed by carbon and nitrogen atoms C3N4 have been predicted in recent years. Their electronic and optical properties have been studied by the first-principles method[17]. Recently, an interesting ferromagnetic ordering has been found in a prospective spintronic material, BCN monolayers[18], as well as the amusing Dirac cone have been found in the organic C4N3H monolayer theoretically[19]. Some other 2D carbon nitride materials can be used in ‘post-silicon electronics’[20]. However, carbon nitride bilayers with 1:1 stoichiometry have not been proposed or studied until now.

We are trying to explore new 2D materials composed of carbon and nitrogen atoms with novel properties based on density functional theory (DFT). As we know, there are four outermost electrons in the carbon atoms, and five in the nitrogen atoms. In order to satisfy the octect rule, the carbon atoms should form four covalent bonds, while the nitrogen atoms should form three covalent bonds. In this paper, we report a systematic study of the stabilities and electronic structures of the previously unknown phases of 2D binary compounds that have formula C2N2, or C4N4, and we name them α-C2N2 and β-C2N2, and γ-C4N4, respectively. Our results indicate that all of the three phases are dynamically stable, and all of them are indirect band gap insulators. Some other similar binary compounds with wide band gaps have been predicted in our other studies[21].

Results

Structure and Stabilities

The fully relaxed geometric structures of three carbon nitride bilayers are shown in Fig. 1. The α-C2N2 and β-C2N2 bilayers have stable honeycomb structure, as shown in Fig. 1(a,b). Therefore, the Brillouin zone (BZ) is a regular hexagon, as shown in Fig. 2(g). In Fig. 1(a,b) the basis vectors are denoted by red arrows and the unit cells containing two carbon atoms and two nitrogen atoms, respectively, are denoted as green rhombs. Both of them have the same structure parameters with the in-plane lattice constant 2.35 Å, C-C bond length 1.62 Å, and C-N bond length 1.44 Å. However, they belong to the different point groups: D3h for α-C2N2 and D3d for β-C2N2, respectively. The thickness of bilayer structure is 2.60 Å. The bond angle between the two C-N bonds is 108° 57′, and the bond angle between the C-N bond and C-C bond is 109°58′. Compared with the bond angle of a regular tetrahedral carbon with sp3 hybridization, 109°28′, the α-C2N2 and β-C2N2 bilayers show a slight deviation from the standard tetrahedral carbon structure. For the γ-phase, we find that a1 is 2.35 Å, a2 is 3.99 Å, the N-N bond length is 1.47 Å, the intralayer C-C bond length is 1.54 Å, and the C-N bond length is 1.45 Å, and the thickness of the bilayer is 2.87 Å. The interlayer C-C bond length is 1.55 Å, which corresponds to the interlayer distance. The basis vectors and unit cell are denoted in Fig. 1(c). Its unit cell contains four carbon atoms and four nitrogen atoms. The BZ has a rectangular shape, as shown in Fig. 2(h). Unlike the α- and β-phases, the γ-C4N4 has D2h symmetry. All three bilayers can satisfy the octect rule.

Figure 1

Geometric structures of three carbon nitride bilayers. Top and side views of (a) α-C2N2, (b) β-C2N2 and (c) γ-C4N4. The grey and light blue spheres denote carbon and nitrogen atoms, respectively. Green rhombs are the unit cells, and the red arrows labeled a and a are the basis vectors. 

Figure 2

Phonon dispersions of (a) α-C2N2, (b) β-C2N2, and (c) γ-C4N4 bilayers along the HS lines in the BZ calculated by combining VASP with Phonopy, respectively. Phonon dispersions of (d) α-C2N2, (e) β-C2N2, and (f) γ-C4N4 bilayers along the HS lines in the BZ calculated by RESCU. (g) is the BZ and the HS points for α-C2N2 and β-C2N2, while (h) is for γ-C4N4. 

In order to explore the stabilities of three carbon nitride bilayers, their phonon dispersions have been calculated by both of VASP + Phonopy and RESCU, as shown in Fig. 2(a–c),(d–f), respectively. The absence of vibrational modes with imaginary frequency in the whole BZ, for all of the three phases, suggests that all of them are dynamically stable, and that no frequency gap is found in the dispersions. Moreover, the out-of-plane mode (ZA) around Γ point from RESCU is much stiffer than that from VASP + Phonopy. The wave velocities obtained are listed in Table 1, which shows that both the α-C2N2 and β-C2N2 bilayers are mostly isotropic. However, the γ-C4N4 bilayer exhibits great anisotropy in its mechanical properties, based on their different velocities from the Γ point to the other high symmetry (HS) points.

Table 1

Wave velocities of acoustic modes (m/s).

acoustic modes	α-C2N2		β-C2N2		γ-C4N4
	Γ→M	Γ→K	Γ→M	Γ→K	Γ→X	Γ→Y	Γ→M
ZA	1042(7214)	1208(8297)	1222(8222)	1196(9299)	607(3924)	2101(4729)	2052(4081)
TA	11554(8797)	11554(9782)	11492(9772)	11498(11782)	15706(5095)	11049(5040)	12631(6006)
LA	18605(15245)	18605(17306)	18411(14446)	18409(16066)	18789(8249)	17759(8689)	18022(7577)

The parenthesized numbers are results fitted from RESCU calculations, while the ones that are not parenthesized are calculated by VASP + Phonopy.

Due to the similarity between the α-C2N2 and β-C2N2, we have proposed a possible conversion path way between the two phases. From the CINEB calculation, the energy change during the transformation is shown in Fig. 3(a). Since the CINEB can optimize one of the images to the transition state (TS), we study its vibrational mode with imaginary frequency, which corresponds to the conversion pathway, as shown in Fig. 3(b). Obviously, the transition state during the conversion involves the breaking and reforming of C-N bonds, and the transition has a comparatively high energy of 1.90 eV per unit cell with respect to the α-C2N2 and β-C2N2 bilayers, as shown in Fig. 3(a), which is equivalent to an activation energy of 183 kJ/mol, a formidable energy barrier that can prevent reactions under normal conditions.

Figure 3

(a) Change of total energy of per unit cell during the transformation from α-C2N2 to β-C2N2. The energy of α-C2N2 is set to zero. Blue dots correspond to the inserted images during the CINEB calculation, the point A corresponds to the α-C2N2 bilayer, the point B corresponds to the β-C2N2 bilayer, and the point C corresponds to the transition state. The red curve is the spline interpolation. (b) The geometric structure of the transition states. The vibration mode with imaginary frequency is presented by the double-ended arrow, corresponding to the conversion path along the transition.

Electronic Structures

In this section, the electronic structures of the three bilayers will be discussed. In order to interpret the orbital properties of the phases, the projected density of states (PDOS) have been calculated. The PDOS of orbitals of α-C2N2 are shown in Fig. 4(a), with the upper and lower panels being the PDOS of carbon and nitrogen, respectively. As shown in Fig. 4(a), most of the low-energy states come from the s orbits of the nitrogen atoms, while states from valance bands near the Fermi energy come from the pz orbit of the nitrogen atoms. In the conducting bands, the px+y orbit of both elements make a greater contribution than the s orbits and the pz orbits, and the carbon plays a major role in forming the conducting band. The PDOS of the β-C2N2 and γ-C4N4 are similar.

Figure 4

Electronic structures of the α-C2N2 bilayer. (a) The density of states: the upper and lower part of the graph represents the C and N components, respectively, and the red, blue and cyan lines represent the s, px+y, pz orbital projected densities of states, respectively. (b) The band structures: blue solid/red dotted lines represent the PBE/HSE calculation, the red shaded area indicates −2~0 eV with respect to the Fermi level of the PBE calculation. (c) Top and (d) side view of the charge density of the states in the shaded area of (b). Brown and silvery joints correspond to the carbon and nitrogen atoms, respectively. The Fermi energy is set to zero.

The electronic bands of the α-C2N2 calculated with PBE and HSE functionals are plotted in Fig. 4(b). The results indicate that the α-C2N2 is an indirect band gap semiconductor, with a 3.76 eV band gap by PBE calculation and a 5.22 eV band gap by HSE calculation. The PBE calculation underestimates the band gap by 1.46 eV compared with HSE calculation. The PBE and HSE eigenvalues are nearly ‘degenerate’ near the Fermi level. Figure 5(b) shows that β-C2N2 is also a semiconductor with an indirect band gap of 4.23 eV by PBE calculation and 5.75 eV by HSE calculation. The PBE calculation underestimates the band gap by 1.52 eV compared with HSE calculation. The γ-C4N4 has an indirect band gap of 2.06 eV by PBE calculation, and the HSE calculation gives a wider one with 3.53 eV, which is 1.47 eV wider than the PBE band from Fig. 6(b). The band gap of α-C2N2 is close to the band gap of nitrogene[8,9], and the β-C2N2 has the widest band gap, while the γ-C4N4 has the narrowest one among them. Furthermore, while the conducting band minimum (CBM) of the three carbon nitride bilayers are all located at the M point, and the valence band maximum (VBM) of the α-C2N2 and β-C2N2 bilayers are located along the K-Γ line, while the VBM of the γ-C4N4 is located at the Γ point.

Figure 5

Same as Figure 4 for  β-C2N2 bilayer.

Figure 6

Same as Figure 4 for γ-C4N4 bilayer.

In order to clarify the bonding mechanism of the structures, the partial charge densities (PCD) are calculated. The charge densities from states from −2~0 eV with respect to the Fermi level of α-C2N2, as indicated in the red shaded area in Fig. 4(b), are shown in Fig. 4(c,d). The charge density distributed around the nitrogen atoms exhibits the character of a pz oribital that points out of the surface. By contrast, the charge density around the C-C bond shows the character of the σ bonds. Most states are located around the nitrogen atoms, which are mostly formed by the pz orbit, because the nitrogen atom possesses one more valance electron than the carbon atom. The other two phases exhibit similar phenomenon, see Figs 5(c,d) and 6(c,d). The nitrogen atoms and carbon atoms in the three structures prefer the sp3 hybridization by combining the geometric structures with orbital characteristics.

Effect of Tensile Strain

As previous studies have shown, the electronic structures of nitrogene[9,12] and phosphorene[22-24] are sensitive to tensile strain, which means, the band structure can be tuned by strain. In this section, the band structures of the three bilayers have been engineered by tensile strain. All of them undergo an insulator-mental transition under tensile strain. For α-C2N2, as shown by the red dots connected by blue line in Fig. 7(a), the band gap decreases monotonically as the biaxial tensile strain increases, and the band gap closes as the strain reaches 27%. Under the 27% strain, the material turns to two slightly buckled layers with buckling distance 0.007 Å, and the distance between these two layers is 4.446 Å, and the one of nitrogen atoms is 4.453 Å. The bilayer then becomes metallic as the strain continues to increase. For the β-C2N2, the band gap increases until 2%, and decreases monotonically as the biaxial tensile strain increases, and the band gap closes at 22%, as shown by the green dash line in Fig. 7(a). Similar to α-C2N2, the β-C2N2 suddenly becomes two slightly buckled layers under the critical strain of 22%, and the distance between the two layers is 3.67 Å, which is shorter than that for α-C2N2. Because of its different point group and its obvious structural anisotropy, we consider the effect of tensile strain along and , respectively. Strains along both directions can turn the bilayer into a direct semiconductor and ultimately drive an insulator-metal transition. The band gap decreases almost linearly as the tensile strain increases from 4% to 11% along the direction, then closes as the strain reaches 11%. Interestingly, the tensile strain along the direction can turn the bilayer to a direct gap semiconductor with a critical strain of 5.5%, as shown by the red dots connected by the blue line in Fig. 7(b). For the direction strain, the band gap decreases monotonically as the strain increases and closes at 11.5%. There are two near linear relations from 0% to 7% and from 7% to 11.5%, respectively, and the critical strain of the indirect-direct transition is 7.5%, as shown the green squares connected by the dark green line in Fig. 7(b). The γ-C4N4 bilayer becomes a puckered anisotropic Dirac semimetal as the strains in both directions reach the critical strain.

Figure 7

Dependence of band gap on the strains. (a) The red dots connected by the blue line and green squares connected by the dark green line represent the dependence of the gap on the biaxial strain of α-C2N2 and β-C2N2, respectively. (b) The red dots connected by the blue line and green blocks connected by the dark green line represent the dependence of the gap of γ-C4N4 under strain along and directions, respectively.

To understand the process of the band gap closing of the γ-C4N4 bilayer under the strains in detail, the band structures around the critical strain have been plotted in Fig. 8. The red line and blue dashed lines represent the band without and with considered SOC, respectively. The Dirac points are denoted as D0, D1, and D2, respectively. Under the strain along the a, as shown in Fig. 8(a–c), the bilayer becomes a gapless anisotropic Dirac semimetal at the 11% strain, and the CBM and the VBM are degenerate at the Γ point. As the strain along the a direction exceeds the critical value, the anisotropic Dirac point shifts away from the Γ point along the Γ-X line. As an example, under 12% strain, a linear dispersion relation appears. The Fermi velocity can be obtained by formula , where is the reduced Planck constant. The largest Fermi velocity is 2.69 × 105 m/s at D0 point as shown Fig. 8(c), while, for the strain along a, the band gap closes at 11.3% and the CBM and VBM are degenerate at the Γ point. As the strain along a increases further, two anisotropic Dirac points appear and move toward to the X and Y points in the BZ, respectively, as shown in Fig. 8(e,f). This is obviously different from the one obtained by stretching along a. Under the 12% strain along a, the largest Fermi velocity at the D1 point is 1.01 × 105 m/s and the largest velocity at the D2 point is 4.73 × 105 m/s. The energy bands with SOC and without SOC are degenerate in all of the BZ, which indicates the SOC can be neglected in this material. The linear dispersion relation always appears when the strains in the two directions reach some critical values or, in the other words, when the anisotropic Dirac points are robust under the strains.

Figure 8

Band structures of the γ-C4N4 bilayer under different strains. Bands of bilayer under strain with (a) 10.5%, (b) 11%, (c) 12% along the a direction. Bands of bilayer under strain with (d) 11%, (e) 11.3%, (f) 12% along direction. The red line and blue dashed line represent the bands without and with considered SOC, respectively. The Fermi energy is set to zero. The Dirac points are denoted as D0, D1, and D2, respectively.

Discussion

Three novel phases of carbon nitride bilayers, denoted α-C2N2, β-C2N2 and γ-C4N4, respectively, which satisfy the octect rule, have been predicted by the density functional theory. Their phonon dispersions, transition states and electronic structures have been calculated. Our results show that all of them are stable, while the conversion between the α-C2N2 and β-C2N2 bilayers is found to be hard due to the high energy barrier. The acoustic modes of the phonon dispersions, especially the out-of-plane mode (ZA), appear to be harder than in low-buckled silicene[3], and black and blue phosphorus[22]. The phonon modes of the three bilayers calculated by the new method RESCU are consistent with the VASP + Phonopy, but harder in the lowest acoustic modes. In terms of energy, both α-C2N2 and β-C2N2 are nearly equally stable. However, the β phase is 0.403 meV per atom higher than the α phase, and the γ phase is 259.59 meV per atom higher than the α phase. The formation energies of the proposed phases with respect to graphite and nitrogen gas can be obtained by the formulaHere Egraphite and are the energies of graphite and nitrogen gas, n and n are the number of C atoms and N atoms in the unit cell, respectively. For α-C2N2 and β-C2N2, n = n = 2, while for γ-C2N2, n = n = 4. The formation energies of α-C2N2, β-C2N2, and γ-C4N4 are 409.552 meV/atom, 409.955 meV/atom, and 669.142 meV/atom, respectively. Furthermore, the energy barrier between α-C2N2 and β-C2N2 is 5 meV/atom higher than that from black to blue phosphorus[23]. The electronic bands calculated both by PBE and HSE methods indicate that the three bilayer structures are all indirect semiconductors. Among them, the β-C2N2 bilayer has the widest band gap. Both the phonon dispersions and electronic bands show strong anisotropic properties due to the novel lattice structure of the γ-C4N4 bilayer. The van Hove singularity points from the DOS are similar to graphene[25], strong electron-electron interactions can be established by applying a strong perpendicular magnetic field. In this case, superconductivity may appear[26,27]. The three bilayers undergo the insulator-metal transitions under strain and the γ-C4N4 bilayer can become a distorted Dirac semimetal with high Fermi velocity under the axial tensile strain. The optimized structures under the strains reveal that the three bilayers are stable as well as flexible, which implies that they may be grown on the suitable substrates. For these bilayers, the critical strains under which the insulator-metal transition take places are greater than those of phosphorene[22,24]. The Fermi velocities under 12% strain in both directions have the same order of magnitudes as the octagons and pentagons graphene-Z[28], as well as the Dirac monolayer TiB2 structure[29]. Obviously, the γ phase under the strains is a distorted Dirac semimetal with high Fermi velocity[30], which possesses more extensive applications than the symmetrical one due to their anisotropic electronic properties at the Dirac point[19]. The α-C2N2 and β-C2N2 bilayers predicted in this work could perhaps be synthesized by the self-condensation scheme from g-C3N4 materials[31], or via the wet-chemical reaction as reported for C2N-h2D[32]. These novel bilayer materials may have promising applications in fields such as semiconductors, spintronics[18], batteries, supercapacitors[33] and energy transducers[34].

Methods

The optimization and electronic structure calculations have been performed by the Vienna Ab initio Simulation Package (VASP) code[35] within projector augmented wave with cut-off energy of 650 eV, and the exchange-correlation is treated by Perdew-Burke-Ernzerhof (PBE) functional[36]. In addition to PBE calculations, the Heyd-Scuseria-Ernzerhof (HSE06) hybrid functional with a screened Coulomb potential[37,38] is also used for electronic structure calculation. The energy convergence criterion for electronic iteration in VASP is set to be 10−8 eV. The vacuum between two different bilayers is no less than 20 Å. Structures are relaxed until the net force on each ion is less than 10−4 eV/Å. The Brillouin zone was sampled with a Γ-centered grid of 20 × 20 × 1 k points. Phonon dispersions were calculated using both VASP interface of Phonopy[39] and a Real space Electronic Structure CalcUlator (RESCU)[40], respectively, for comparison. In the phonon calculation with VASP + Phonopy, force constants were obtained from a 4 × 4 × 1 supercell for α-C2N2 and β-C2N2, and a 4 × 6 × 1 supercell for γ-C4N4. In the phonon calculation with RESCU, the dimensions of the supercell become 2 × 2 × 1 for α-C2N2, 4 × 4 × 1 for β-C2N2 and γ-C4N4. In RESCU calculation, which is performed in real space grid, double-zeta numerical atomic orbital basis is used to generate Hilbert subspace for the Kohn-Sham DFT solver, with Local Density Approximation (LDA) as the exchange-correlation functional. The grid spacing is no greater than 0.4 Bohr. Self-consistence iteration is performed until the global charge variation is less than 10−5 electrons for each valence electron, and the total energy variation per valence electron is less than 10−6 Hartree. The conversion path between α-C2N2 and β-C2N2 was calculated by the Climbing Image Nudged Elastic Band method (CINEB)[41-43], in which five images were inserted in-between.