Ab initio study of the effect of 2D layer rippling on the electronic properties of 2D/H-terminated diamond (100) heterostructures

We report an ab initio study of the effect of rippling on the structural and electronic properties of the hexagonal Boron Nitride (hBN) and graphene two-dimensional (2D) layers and heterostructures created by placing these layers on the Hydrogen-terminated (H-) diamond (100) surface. Surprisingly, in graphene, rippling does not open a band gap at the Dirac point but does cause the Dirac cone to be shifted and distorted. For the 2D/H-diamond (100) heterostructures, a combined sampling and a clustering approach were used to find the most favorable alignment of the 2D layers. Heterostructures with rippled layers were found to be the most stable. A larger charge transfer was observed in the heterostructures with rippled hBN (graphene) than their planner counterparts. Band offset analysis indicates a Type-II band alignment for both the wavy and planar heterostructures, with the corrugated structure having stronger hole confinement due to the larger valence band offset between the hBN layer and the H-diamond (100) surface.


Introduction
Diamond has immense potential as a wide bandgap semiconductor for use in next-generation high-frequency and highpower electronics, owing to its high hole mobility along with its high thermal conductivity and breakdown field [1]. Unfortunately, however, it is very difficult to engineer high carrier concentrations in diamond using traditional impurity doping due to the low solid-solubility of most dopants. One promising strategy to overcome this challenge is to surface transfer dope diamond. In this technique, hydrogen-terminated diamond (H-diamond) (100) surface is allowed to interface with the adsorbates such as NO 2 , NO, O 2 , etc. Since the lowest unoccupied molecular orbitals (LUMO) of these adsorbates lie below the valance band edge of the H-diamond (100) surface, the charge equilibrium at the interface leads to the charge transfer from diamond to the adsorbates. This charge transfer process results in the formation of a thin, highly conductive hole layer in the diamond surface. This highly conductive hole layer acts as a channel layer in the field effect transistor (FET) [2,3]. Furthermore, the C-H dipoles at the surface of H-diamond induce a strong upward bending of the valence band at the surface, leading to easier charge transfer to the acceptor layer. Doping densities achieved through this approach have been reported in the range of 10 -14 to 10 -12 cm −3 , making surface transfer doping a viable solution for creating diamond-based FET devices [4].
Recently, transition metal oxide (MoO 3 , WO 3 , and V 2 O 5 ) acceptor layers garnered interest from the academic community for their stability and ability to sustain harsh experimental conditions while preserving hole channel in the diamond surface [5][6][7][8][9]. Though rich in chemistry and compatible to the existing device fabrication technology, these oxide-based acceptor layers suffer from the higher degree of surface inhomogeneity and dangling bonds, leading to lower channel mobility and charge density. In an earlier work [10], we reported that a sheet of hBN or graphene interfaced with the H-diamond (100) surface-the technologically relevant surface for diamond-based electronic Invited Paper devices-could act as a protective layer that preserves the p-type bulk-like conductivity in the H-diamond (100) surface. During the analysis of the structural alignment of the 2D layers on the H-diamond (100) surface, we also noticed a strain-mediated 2D layer corrugation or rippling. In this present manuscript, we extend our work to examine the effects of the 2D layers' rippling on the electronic structure of the isolated 2D layers and 2D/Hdiamond (100) heterostructure.
The absence of dangling bonds at the surface of 2D materials enables the stacking and construction of various van der Waals (vdW) heterostructures, even in cases with substantial lattice mismatch [11]. This allows combining different 2D layers with each other or three-dimensional (3D) materials to mitigate limitations in engineering high-power electronic devices [12]. In many cases, such as the presented hBN(graphene)/H-diamond (100) system, there is a lattice mismatch between the 2D layer and the 3D material [13,14]. In these cases, even though the weak vdW interfacial bond means that the energy penalty for the mismatch is low, the intrinsically low bending stiffness of 2D materials means that the energy penalty can sometimes be further reduced by buckling the 2D layer to accommodate any compressive strain needed to bring the 2D layer and substrate into registry. Moreover, rippling is a common occurrence in computer simulations of 2D heterostructures as an artifact of the periodic boundary conditions and the need to enforce lattice matching in a computationally tractable system size. In either case, elucidating the effects of rippling on electronic structure is essential for understanding the resulting heterostructure's electrical transport properties. If the ripples are real, we need to understand them, and if rippling is an artifact of computation, we need to be able to discount the ripple's effect from any predictions made or inferences drawn.
In experiment, ripples are observed both in free-standing graphene due to instability of the 2D lattice [15] or lattice perturbation [16], and in 2D-based heterostructures, due to the interaction with the underlying substrate such as SiC (0001) [17] and Ge [18]. Rippling in graphene interfaced with Ge substrate was found to result in a slight n-doping of graphene in the rippled regions [18]. Since the ripples are developed in response to an external perturbation, they could potentially affect the electronic properties and subsequently the device performance. Controlling the degree of rippling in graphene is suggested as a way to engineer devices with tailored bandgaps and functionality. A previous study reported generating controlled ripples in suspended graphene through thermally applied strains [16].
In computation, the generation of large strains is an unavoidable part of modeling heterostructures with lattice mismatch. Understanding how ripples affect a system's properties helps drawing a better comparison between the computational predictions and real-world applications. Previous theoretical studies have focused on elucidating the thermodynamic and mechanical properties of rippled hBN [19,20]. Others have studied graphene under an imposed modulated potential with a similar wavelength to the ripples in our work, finding that the imposed potential opens the Dirac cones at k-points that are not parallel to the wavevector of the modulation [21]. Going beyond these works, in this paper, in addition to studying the effect of the ripples on the 2D layer, we elucidate their impact on the computed interaction of the 2D layer with a substrate. One merit of this is to determine whether rippling works in favor of electronic devices. In this case, to obtain better charge transfer properties, we can experimentally devise device geometries with the layers forced to form ripples under some geometrically imposed strain.

Models and methods
We are specifically interested in the (100) surface of diamond, as this is the dominant surface in the experimental growth of diamond films via chemical vapor deposition [22,23] and exhibits smooth bulk-like to (2 × 1) reconstruction [24]. Besides, the high stability of the 2 × 1 reconstruction of the (100) surface enables the growth of defect-free facets during fabrication [25]. This surface undergoes a 2 × 1 surface reconstruction that affects the pattern of H termination. Computing the electronic structure of this substrate forming a heterostructure with a hexagonal 2D layer requires finding a supercell that is small enough for efficient computation, but that retains the 2 × 1 pattern of the diamond while accommodating the 2D layer with minimal strain. For this, lattice parameters of graphene and hBN are not well-matched with diamond. The best compromise structure is a 2 √ 2 × 1 √ 2 unit of 2 × 1 reconstructed H-diamond (100) surface, which was used to set the simulation cell size, and then to rotate and strain a 5 × 1 unit of the hexagonal layer to fit into this box. The stages for bringing the hBN layer into registry with the H-diamond (100) surface are illustrated in Fig. 1a-c. While the overall area dilation in the hBN is small, the hBN undergoes a significant shear strain after this process. If the atoms are displaced slightly to break symmetry before relaxation (from step 2 to 3a in Fig. 2), the sheet buckles as it relaxes to reduce the compressive strain, generating an undulation along the long axis over the compute cell.
The generated ripple has a wavelength equal to the length of the compute cell, 10.089 Å. Previous studies have suggested that since ripples are generated as a result of an interplay between the strain and the periodicity of the substrate, the wavelength will change with increasing the cell size until reaching a critical value where the structure is no more stable [26]. Under a 5% compressive strain, this critical rippling wavelength has been predicted as 1-2 nm [27], comparable to the ripples in our system. Although an analysis of layers with various degrees of rippling is out of the scope of this study, we believe that the ripples in our calculations are an artifact of the imposed compute geometry and not yet at the critical wavelength. Our goal is to understand the effect of rippling on the heterostructures' structural and electronic properties, and also provide a comparison against the generated models with a strained-planar hBN (graphene) layer. To obtain these, the structures were relaxed, with the atoms in the 2D layer permitted to move in plane but constrained from moving out of plane. To differentiate the former and latter heterostructures, hereon, we refer to them as wavy and planar hBN(graphene)/Hdiamond (100), respectively. The C-C bond length in graphene is 1.42 Å, which is similar to the 1.45 Å B-N bond length in  hBN, and so the same supercell structure was used for graphene. The slightly shorter bond length in the graphene means that it has a larger tensile strain but smaller compressive than the hBN, and as a result produces ripples with a smaller amplitude. Further details of the relaxation and structural transformation are provided in our earlier article [10].
The final step in constructing the strain-compensated heterostructures was to find the lowest energy alignment of the 2D layer on the H-diamond (100) surface. To achieve this goal efficiently, we performed a coarse sampling of the potential energy landscape by rastering the position of the 2D layer over the diamond substrate to generate 38 different starting structures. For each resulting structure, the initial 30 steps of structural relaxation were performed, after which the energies were mapped and compared (Fig. 1d). The aim of this was not to find the minimum energy position, but rather to find the watersheds between different potential energy wells. To gain a better understanding of the structure-energy distribution, we performed the k-means clustering method [28], a commonly used algorithm for cluster analysis that groups similar items into k clusters. The k-means clustering error converged with k = 4 , i.e., four clusters yield the most accurate partition of the structures, shown in Fig. 1e. The group of initial structures with the lowest final energies, denoted by "The low energy group", are located on the y = 0.33 line. In other words, we arrive at the lowest energy structures starting from structure for any x shift in the 2D sheet providing the y shift is 0.33 of the box width. This implies that there is little to prevent the sheet from sliding in the x direction during minimization. However, in the 2D/H-diamond (100) heterostructures, the hydrogens are positioned on the outer edges of the H-diamond (100) surface, which is expected to inhibit the atoms in the layer from shifting in y direction during the structure relaxation. Based on this, we selected the structure from this group with the lowest energy and continued its structural relaxation to minimize its energy. As the next step, we found the optimal vdW spacing between the constituent layer and the H-diamond surface. The procedure for finding the optimum vdW spacing is fully discussed in our earlier publication [10]. Accordingly, structures with a vdW spacing of 3.0 Å were found to yield convergence. The H-diamond (100) surfaces were separated by a vacuum spacing of 20 Å to avoid interactions between the adjacent blocks. In all structure relaxation calculations, the surface hydrogens and the first four carbon layers of the H-diamond (100) surface were allowed to relax while the remaining layers of diamond were fixed.

Results
The results and discussions are presented in two sections: (a) atomic structure and (b) electronic structure of the isolated layers and the resulting heterostructures.

Atomic structure
To understand the effect of the 2D layers' rippling on the heterostructures, we need to evaluate each contributing factor to the layers' transformation and the heterostructures' construction. Hence, we computed the energy at each step of the process for mating the 2D layer with the diamond substrate. Figure 2 shows the stages for transforming the hBN (graphene) layer and aligning the transformed layer on the H-diamond (100) surface. As illustrated, the contributing factors include imposing strain on the layer, full relaxation of the layer through which the layer develops a rippled structure, adhesion of the layer on the H-diamond surface, and the vdW interaction between the constituent layer and the substrate.
The strain energy is defined as the energy gained when the layer is strained (Fig. 2, subsets 1 and 2). The energy released when the strained layer develops ripples is considered as the rippling energy, equivalent to the energy difference between stages 2 and 3a in Fig. 2. The adhesion energies ( E ad ) and vdW energies ( E vdW ) were computed using the Here, E S/D denotes the total energy of the heterostructure. To compute the adhesion energy values, the fully relaxed atomic structure for the layer and H-diamond (100) surface prior to integration, illustrated in Figs. 2, 3a and 3b, were used as E S and E D , respectively. Whereas, to evaluate the vdW interaction strength, we used the frozen atomic structure of the constituent layer and the substrate isolated from the fully relaxed 2D/H-diamond (100) heterostructures. For consistent evaluation, the energy values were normalized with respect to the number of atoms in the layer.
To understand the effect of the vdW correction scheme on the structural properties, we evaluated the adhesion energies of simulations with DFT-D2, zero damping DFT-D3, and DFT-D3 with the Becke-Jonson (BJ) damping ( Table 1). The results indicate that the choice of vdW correction, at least within the tested approaches, does not significantly affect the adhesion energies. Quantitively, the change in adhesion energy was 0.008 eV. In addition, we did not observe a significant difference between the relaxed structures using different vdW corrections. Since all the correction schemes resulted in similar adhesion energy and vdW gap, we have selected the D2 correction scheme for our work.
The energy values for all contributing factors are presented in Table 2. The results show a higher energy cost to enforce epitaxy on the hBN layer compared with graphene due to a larger imposed strain on the hBN layer. Rippling of the strained layers lowers the energy and makes the structure more stable. The process of rippling in the hBN and graphene releases 0.112 and 0.047 eV/atom of energy, respectively. The larger energy loss during rippling in hBN is caused by the higher applied strain during modeling. Although the energy gain through imposing strain on the layers is not fully compensated by the rippling, the negative values of adhesion energies indicate that the adhesion is exothermic and hence energetically favorable for all the studied heterostructures.
Here, in addition to the adhesion energy, which incorporates all the contributing factors, including the lattice deformation (buckling) and the vdW interactions, we also computed the energies for the bare contribution of the vdW interactions. Although the adhesion energy values are slightly stronger for the heterostructures constituting of planar layers due to strain, the average vdW interactions between the constituent layer and the H-diamond surface are stronger for the heterostructures composed of wavy layer. In addition, for both planar and wavy models, the adhesion and vdW interactions are stronger in the heterostructures composing of hBN as compared with the ones with a graphene layer. The higher interaction between hBN and H-diamond surface can be related to the higher degree of rippling, as well as the polar nature of hBN. These results are consistent with an earlier study which predicted enhancement of binding strength for adhesion of Hydrogen to rippled graphene compared with a planar layer [26]. The bending of the bonds at the crest and trough of each ripple causes sp 2 to sp 3 orbitals rehybridization, leaving available (free) orbitals to participate in bonding with the substrate. Such improvement in reactivity due to ripples in graphene has been found beneficial for the hydrogenation and dehydrogenation process [26]. A previous study reported a quantified correlation between the graphene rippling degree and hydrogen bonding energy [29]. In addition to evaluating the energy path, to compare the overall stability of the planar    and wavy heterostructures, we evaluated their total energies. We also measured the average C-H bond lengths at the surface of H-diamond in the heterostructures, as illustrated in Fig. 3, the results of which are tabulated in Table 2. These show that in the wavy layer/H-diamond heterostructures, the surface C-H bonds are slightly contracted, whereas these bonds in both planar layer/H-diamond systems are slightly stretched compared with the C-H bonds at the bare H-diamond (100) surface with a bond length of 1.10 Å.

Electronic structure
In the following two sections, we present our results on the effect of rippling on the electronic properties of: 1. 2D layers, and 2. 2D/Hdiamond (100) heterostructures. We used PBE functional for computing the electronic properties. The PBE functional is well-known for underestimating the band gap. However, since our goal here is to draw a clear comparison rather than providing exact values, the PBE functional was found to provide sufficient accuracy while being computationally more affordable. To prove the validity of this claim, we compare the results computed using PBE [30] and HSE06 [31] functionals for the 2D layers, which require less computational resources as compared with the heterostructures.

2D layer
To elucidate the effect of strain and rippling on the electronic properties of the 2D layers, we calculated and compared the electronic band structures for the planar and wavy hBN (graphene) layers. In pristine (planner) graphene, the Dirac cone is located at the high symmetry K point at the corner of the first Brillouin zone (BZ). However, straining the graphene layer distorts the perfect honeycomb lattice and BZ of the planner graphene, shifting the K point to a new position, as shown in Fig. 4a. On computing the band structure at this position in the BZ, one finds a small apparent band gap opening in both the planar and rippled graphene. Though this gives an impression of strain-induced band gap opening, this is mainly due to a subtle shift in the location of the Dirac cone to new locations that we call K S in the strained planar unit cell and K R in the 5 × 1 supercell of rippled graphene. To find the true location of the Dirac cone, we computed the band energies on a fine mesh of k-points surrounding the corner of the strained BZ, and then fit the equation of a strained and biased cone to the conduction and valence band energies at these k-points. The fitting function had the form: where the vector s describes the bias or tilt to the cone, and matrix L = ε 11 ε 12 ε 12 ε 22 describes the cone distortion. Terms k o and E o give the location and energy of the cone's apex. After fitting, the location of the cone center at k o was taken as the center for generating a second refined k-point mesh and the process was repeated as shown in Fig. 4. By the third mesh refinement, we identify a true gapless Dirac cone. Independent fits to the valence and conduction bands (Fig. 4c) produce cones with E o and k o that are identical for both the strained-planar structure (K S ) and the rippled supercell (K R ), that is, they align perfectly in k and their apices meet. The elliptical distortion of the Dirac cone can be seen clearly in Fig. 4c. Although it is enough to see an apparent band gap opening, K S is only minutely shifted from the BZ corner at K, and it is likely that this difference is due to rounding error in the specification of the lattice vectors. However, for the rippled 5 × 1 supercell of graphene the K point is folded into the interior of the supercell's BZ, and so is no longer anchored to a high symmetry point of the BZ. In this case the location of the Dirac cone is no longer constrained by symmetry, and we find that K R is significantly displaced from the corner of the strained unit cell's BZ as can be seen in Fig. 4b. Figure 4d and e show the band structures of the strained and rippled graphene in isolation. Figure 4d shows a conventional tour through the high symmetry points in the strained BZ, while plot 4e shows the band tour that is directed through the true Dirac cones at K S , K ′ S , K R and K ′ R . It can be seen how these shifts in the Dirac cone can mislead one to believe that a band gap has opened when in fact it has not. It highlights that, when computing the band structure of graphene in other situations, one must take care to search around the K point before one can confidently identify a true band gap opening. This process may have been overlooked in some earlier works, and may require reevaluation [32][33][34]. It is also interesting here that, in contrast to the effect in a graphene superlattice formed with a modulated potential [21], the Dirac cones at all six of the high symmetry points in the BZ remain intact. The band structures in Fig. 4 also show that strain causes a significant distortion in the band structure away from the Dirac cone, inducing a large difference in the size of the vertical gap at the M a * and M b * points. With these large changes brought about by strain, it is all more remarkable that the Dirac cone is robust to both straining and rippling.
We evaluated the effect of including spin-orbit coupling (SOC) in calculating graphene's band structure, as illustrated in Fig. 4f-i. For each of these calculations, we repeated out iterative grid search and cone fitting procedure to pin down the location of the Dirac cone to better than nine digits of precision in the wave vector, as is required to resolve the high magnification of the cone apex. Inclusion of SOC might induce an additional gap opening at the Dirac cone and cause a band splitting [35]. Previous ab initio studies report an SOC induced gap in range of 25 to 50 μeV [36]. Including a magnetic field, addition of substrate, lattice distortion, and other changes that cause rehybridization of sp 3 orbitals affects the magnitude of SOC split [35,36]. We found that including SOC induces a gap and splitting of bands in unstrained, strained-planar, and rippled graphene ( Table 3). The magnitude of the gap and splits depend on the choice of parameters in DFT calculations. Though the reported values in this work are different from previous studies, we note that the evaluated gaps and splits were negligible or within the resolution of the numerical accuracy for all the cases.
The band structure of the strained-planar and wavy hBN layers are plotted in Fig. 5b and c, respectively. In both cases, the plots were generated from calculations of a 5 × 1 hBN supercell and show the effective (unfolded) band structure along a tour of the 1st BZ of the strained unitcell through the points indicated in Fig. 5a. The unfolded band structure of the strainedplaner 5 × 1 supercell exactly matched that computed directly in a single strained unitcell. In contrast, the unfolded band structure of the wavy hBN, while containing the same features as the planar hBN, also contains several extra bands that arise The band structure and Dirac cone of rippled and strained but planar graphene in isolation. Plot (a) shows the Brillouin zone (BZ) of the graphene unit cell rotated and distorted so that the 5 × 1 supercell of graphene is commensurate with the diamond substrate. Plot (b) shows a zoom in of the region around the K point that was covered by the first Dirac cone search mesh. The search points after the 1st and 2nd mesh refinement are plotted in green and red. Also shown at point K R is the location of the Dirac cone in the rippled graphene. Plot (c) shows a 3D image of the valence and conduction bands in the most refined k-point mesh. The spheres show the DFT computed band energies, and the surface plot is the fit of the distorted cone (Eq. 1) to these points. Plot (d) show the band structure along a path to the BZ corners for the strained unit cell (orange) and the effective (unfolded) band structure for the 5 × 1 rippled graphene (blue) along the same k-path (for clarity, the band weights for unfolded case are not shown as they are in Fig. 5). Plot (e) shows the same data as in plot (d) except that the path runs to the locations of Dirac cone at K S and K R of each material rather than to the BZ corner. From these latter band tours, it becomes clear that there is no band gap. As the two paths in (d) are slightly different, it can be seen that the way points along the path length do not line up exactly. Plot (f ) shows the 3D image of the valence and conduction bands with spin-orbit coupling. Plots (g-i) show the zoomed band structure in the vicinity of the Dirac point K K ′ on the Ŵ-K-M path for the unstained, strainedplanar, and rippled graphene (unfolded), respectively. In plots (g-i), orange lines indicate calculations including spin-orbit coupling, whereas blue lines represent calculations without spin-orbit coupling.  [37,38]. The PBE measured bandgap of the strained-planar, and wavy hBN layer are 3.96 eV and 4.18 eV, respectively. Using the HSE functional, the energy gaps of the strained-planar and wavy layers increased to 5.34 eV and 5.52 eV, respectively. As compared to the PBE functional, the HSE calculations show a 20-35% increase in the magnitude of the gap, which is consistent with our previous study [10] and Refs. [39,40]. Though the magnitude of EA depends on the choice of functional, using both PBE and HSE functionals yield an increase of EA as the hBN layer is strained and rippled, that is EA plan < EA str−plan < EA rip . One might argue that the observed decrease in EA is in contradiction with the increase of CBM energies when using HSE functionals. The non-monotonous decrements in the EA predicted by the HSE functional is a consequence of asymmetric scaling of the direct and indirect gaps due to the correction to the self-interaction error in PBE functional. Furthermore, the change in surface dipoles as a result of the applied strain and rippling could also alter the relative alignment of CBM band and vacuum level [41].

2D/H-diamond (100) heterostructures
To understand the effect of rippling on surface-layer interactions in the 2D/H-diamond (100) heterostructures, we analyzed the charge transfer between the H-diamond (100) surface and the constituent 2D layer. The interface charge was evaluated by computing the charge density difference using the equation �ρ = ρ (S/D) − ρ D − ρ S where ρ (S/D) , ρ D , and ρ S represent the electron density of the heterostructure, H-diamond (100) surface, and hBN (graphene) layer, respectively. The �ρ for the planar (top-panel) and wavy (bottom-panel) hBN/H-diamond (100) and graphene/H-diamond (100) heterostructures are shows the unfolded effective band structure for a 5 × 1 supercell of strained-planar hBN computed using vaspkit [55], and plot (c) shows that effective band structure for the same system when the sheet is free to form a ripple.
illustrated in Fig. 6a and d, respectively. The areas enclosed by red surfaces are the charge accumulation regions, while the blue regions represent areas with charge depletion. A higher degree of charge transfer is observed for the rippled systems, due to sp 2 to sp 3 rehybridization of orbitals, which leads to the appearance of additional free orbitals ready to participate in surface reactions, resulting in stronger vdW interactions between the rippled layers and the H-diamond (100) surface as reported earlier in the Atomic Structure section. As a result, in the rippled heterostructures, the charge transfer across the interface looks inhomogeneous, with more charge accumulation around the valley, indicating a stronger vdW interaction due to rippling. As compared to the graphene/H-diamond (100) interface, a higher degree of charge transfer at the valley of hBN is observed, which can be attributed to the larger degree of rippling in the hBN layer.
Since the 2D layer rippling influences the interlayer electronic interactions, controlled by vdW forces, we extract the planar charge profiles across the 2D/H-diamond (100) interface by integrating charges along the normal direction [14,42,43], as shown in Fig. 6b and e. As illustrated, for both the heterostructures, the overall feature of the out-of-plane ( z-direction) charge profiles are similar. However, variations in the charge profile along the z-direction are observed in the proximity of the 2D layer. Quantitatively, the percentage charge differences of 51% and 35% around the vicinity of the rippled and strainedplanar layers for the hBN and graphene layers, respectively, confirms a stronger interaction between the rippled 2D layers and H-diamond (100) surface. In addition, a noticeably larger charge difference above the layer is observed for the rippled hBN layer as compared with graphene. Although a noticeable difference in charge transfer is observed between the heterostructures containing planar and rippled layers, our quantitative charge transfer evaluation using the Bader analysis [44] shows that the amount of charge transferred to the substituent layers in all the heterostructures are negligible, and the charge carriers are mainly preserved across the vdW region. Therefore, for all the rippled 2D/H-diamond (100) heterostructures, rippling does not significantly affect charge extraction by the 2D layer. The rippling, if engineered precisely, could have a much larger effect on the transport properties of systems with smaller vdW spacing, particularly where a covalent bond is formed between the layer and the diamond substrate, as suggested in an earlier work [45].
In H-diamond (100), surface dipoles are generated around the C-H bonds due to the electronegativity difference between the C and H atoms. These surface dipoles are also responsible for the modification of the system's work-function ( ψ ), which is defined as ψ = V vac − E f , as illustrated in Fig. 6c and f. Since we are mainly focused on the relative shift in the ψ due to rippling, we utilize PBE functional to extract potential profiles and the Fermi levels for all the simulated systems. The work function of the H-diamond (100) computed using the PBE functional is 4.21 eV. For the planar and wavy hBN/Hdiamond (100) heterostructures, the ψ values were evaluated as 3.93 eV and 4.07 eV, respectively. The ψ of planar and wavy graphene/H-diamond (100) were found to be slightly higher with values of 3.98 eV and 4.13 eV, respectively. These results indicate that integrating the 2D layers and the H-diamond (100) surface reduces the ψ of the H-diamond (100) surface, regardless of the type and atomic composition of the 2D layer, mainly due to the suppression of the surface dipoles. This reduction of ψ due to the adsorption of 2D layers on the H-diamond (100) surface is consistent with the adsorption of organic molecule on Si (100) surface [46] and on diamond (100) surface [47]. In addition, for both types of 2D layers, the heterostructures comprising of a rippled layer have lower ψ. Quantitatively, a 3.6% and 3.8% reduction in the ψ of the hBN/H-diamond (100) and graphene/H-diamond (100) is observed when the constituent layer is rippled. This is mainly due to the rippled induced inhomogeneity in the dipoles along the x-y plane leading to a reduction in the vacuum level. The change in work function, �ψ , according to the Helmholtz equation, is defined as �ψ = �σ · (A · ǫ 0 ) −1 , here �σ , A, and ǫ 0 represent change in surface dipole moment, surface area, and vacuum permittivity. In the case of the simulated heterostructures, the bulk of the contributions to the �ψ (by corollary to the potential drop ( V )) should come from �σ , a parameter which includes all the changes in surface dipoles, vdW interaction and charge transfer between these systems. Therefore, as illustrated in Fig. 6c and f, a larger V is observed for the rippled systems as compared with the planar ones, resulting from the inhomogeneity of dipoles, higher degree of charge transfer and stronger vdW interactions in the rippled heterostructures. Further analysis of the potential profiles reveals some of the device related features, such as potential drop and depletion width at the interface, are larger/wider in the hBN/H-diamond (100) than the graphene/H-diamond (100) system. This is mainly attributed to the relatively larger amount of charge transfer aided by the larger degree of rippling in hBN, which further confirms the roles, albeit small, played by the structural features such as strain and rippling in modifying device related properties in 2D/3D heterostructures.
To obtain a better understanding of the effect of rippling on the electronic properties of the heterostructures, we analyzed and compared the site-resolved local density of states (LDOS) of the wavy and planar hBN(graphene)/H-diamond (100) systems, as illustrated in Fig. 7. The hBN/H-diamond (100) heterostructures with strained-planar and wavy constituent layers, shown in Fig. 7a and c, exhibit semiconducting characteristics with wide energy gaps. The PBE computed bandgap values for H-diamond (100), strained-planar, and wavy hBN/H-diamond (100) systems are 2.62 eV, 2.54 eV and 2.41 eV, respectively. Quantitatively, the energy gap of strained-planar and rippled hBN/Hdiamond (100) is smaller than the gap of H-diamond (100) by 3.1% and 8%, respectively. A larger reduction is observed in the energy gap of the rippled hBN/H-diamond (100) heterostructure as compared with the planar system. This could be due to the stronger vdW interactions and the larger charge transfer, in the rippled structure as compared with the planar system, causing a more significant change in the electronic structure of the heterostructure. In both planar and wavy heterostructures, the H-diamond (100) states mostly contribute to the valence band edge (VBE) while the hBN states mainly contribute to the conduction band edge (CBE). The VBE in the wavy hBN is lower than that of the strained-planar layer by 0.3 eV. This 23% reduction in the VBE energy of the constituent layer as a result of rippling is an indication of a stronger interaction between the wavy layer and the H-diamond (100) surface, confirmed by the larger vdW energy and charge transfer in this system. Both heterostructures show a Type-II band alignment where CB edge of the hBN is below (above) the CB (VB) edge of the H-diamond (100) surface, which is consistent with the band alignment feature extracted by the XPS method in a recent experimental study [48].
Though the conduction band offsets (CBOs) are similar in both the heterostructures, the valence band offsets (VBOs) between the hBN and H-diamond (100) in the wavy heterostructure is relatively larger than the planner heterostructure. The VBO increase of 26% in the wavy hBN/diamond (100), as compared to the planar counterpart, indicates a higher hole confinement at the heterojunctions. The band alignment of the wavy heterostructure predicted using the PBE functional is consistent with our previous study using HSE functional [10]. LDOS of strained-planar and wavy graphene/H-diamond (100) systems are illustrated in Fig. 7b and d. The LDOS plots of graphene/Hdiamond (100) indicate that regardless of the type of the graphene layer, the system remains metallic due to the preserved Dirac cone characteristic.

Conclusion
In this paper, we presented a systematic DFT analysis to compare the structural stability and electronic properties of rippled (wavy) and planar hBN (graphene) layer and hBN(graphene)/H-diamond (100) heterostructures. We carefully selected and analyzed the structural alignments between the 2D layers and H-diamond (100) surface using chemical intuition and cluster analysis. We found that rippling affects the structural and electronic properties of 2D layers and 2D/ diamond heterostructures regardless of the 2D layer type. The rippled hBN layer shows a larger electron affinity and bandgap compared with the planar counterpart, with an indirect gap type in both structures. In graphene, remarkably, rippling does not open a band gap, but does cause the Dirac cone to shift away from the corner of the BZ. Rippling of the strained layers is an exothermic process, i.e., energetically favorable. Also, the heterostructures with a rippled layer were identified as the low-energy structures. Therefore, the strained layers are expected to develop ripples through the growth process, both as an isolated layer and as an interfacial layer on the H-diamond (100) surface. Thus, the planar heterostructures can only be achieved under applying additional constraints to the 2D layer, such as external forces (or clamping) to keep the layer flat. A larger vdW interaction and charge accumulation was observed at the interface between the rippled layers and H-diamond (100) surface. Both wavy and planar hBN/H-diamond (100) heterostructures indicate a Type-II band alignment, but a relatively larger VBO results in stronger hole confinement in the rippled system. Although rippling had a minimal effect on the electronic properties, its impact can be magnified depending on the 2D layer quality, vdW spacing and magnitude of rippling. Hence, understanding the impact of ripples, though small, is essential in determining the overall performance of the fabricated device design.

Methodology
The structural and electronic property analysis in this paper were performed through density functional theory (DFT) calculations using the Vienna ab-initio Simulation Package (VASP) [49,50] For these calculations, the plane-wave-pseudopotentials based on the Perdew-Burke-Ernzerhof (PBE) exchange correlation functionals and generalized gradient approximation (GGA) correction were employed [30]. For computations of electronic properties in 2D layers, we used the HSE06 functionals [31]. A maximum energy cutoff value of 520 eV was used for ionic relaxation of the structures. The Brillouin Zone (BZ) of graphene (hBN) unit cell was sampled using a 12 × 12 × 1 Gamma-centered mesh of k points which was found sufficient to optimize the structures to an accuracy of 10 -4 eV. For the heterostructures' geometric relaxation hence, a reduced 4 × 12 × 1 k-point mesh was used. We used the Grimme's D2 approximation [51], zero damping DFT-D3 [52], and DFT-D3 with the Becke-Jonson (BJ) damping [53] in order to include the effect of van der Waals interactions between the layer and the H-diamond surface. To include SOC in band structure calculations of graphene, we used the SOC implementation in VASP [54]. The convergence criteria for ionic relaxations were set to converge electronic forces on ions to better than 10 -5 eV/Å. For this study, we used a 2 × 1 Hydrogen terminated reconstructed diamond (100) surface with converged lattice parameters of a = 5.04 Å and b = 2.52 Å and a thickness of 10 Å, equivalent to 10 carbon layers, which was found sufficient to yield convergence as discussed in our previous article [10].

Data availability
The datasets generated and/or analyzed during the current study are available from the corresponding authors on reasonable request.

Conflict of interest
The authors declare they have no actual or potential competing financial interests.

Open Access
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.