Mechanical and electronic properties of boron nitride nanosheets with graphene domains under strain

Abstract

Hybrid structures comprised of graphene domains embedded in larger hexagonal boron nitride (h-BN) nanosheets were first synthesized in 2013. However, the existing theoretical investigations on them have only considered relaxed structures. In this work, we use Density Functional Theory (DFT) and Molecular Dynamics (MD) simulations to investigate the mechanical and electronic properties of this type of nanosheet under strain. Our results reveal that the Young's modulus of the hybrid sheets depends only on the relative concentration of graphene and h-BN in the structure, showing little dependence on the shape of the domain or the size of the structure for a given concentration. Regarding the tensile strength, we obtained higher values using triangular graphene domains. We find that the studied systems can withstand large strain values (between 15% and 22%) before fracture, which always begins at the weaker C–B bonds located at the interface between the two materials. Concerning the electronic properties, we find that by combining composition and strain, we can produce hybrid sheets with band gaps spanning an extensive range of values (between 1.0 eV and 3.5 eV). Our results also show that the band gap depends more on the composition than on the external strain, particularly for structures with low carbon concentration. The combination of atomic-scale thickness, high ultimate strain, and adjustable band gap suggests applications of h-BN nanosheets with graphene domains in wearable electronics.

---

1. Introduction

Produced from graphite by using micromechanical cleavage technology,^1,2 graphene is a 2D material composed of sp²-hybridized atoms arranged in a honeycomb-like hexagonal lattice. Various applications have been proposed for graphene in modern technologies^3,4 due to its extraordinary properties, which include: high electron mobility at room temperature,⁵ high thermal conductivity,⁶ and broad light absorption spectrum.⁷ Moreover, graphene has motivated the search for other 2D materials arranged in a honeycomb-like hexagonal lattice. Well-known examples include hexagonal boron nitride (h-BN),^8,9 silicene,¹⁰ and many transition metal dichalcogenides.¹¹ In particular, h-BN, also known as white graphene, is an insulator material (energy gap >4 eV).^8,9,12 This 2D material also presents remarkable properties, such as high thermal stability,¹³ high resistance to oxidation,¹⁴ and strong optical absorption in the UV region.^15,16 Several synthesis methods have been used to achieve large area h-BN monolayers. For instance, h-BN has been grown through ultra-high vacuum chemical vapor deposition (CVD).^17,18 We can list applications for h-BN as atomically flat substrates,¹⁹ as 2D dielectric materials,²¹ and in lubrication.²⁰

We can certainly also include the mechanical properties of graphene and h-BN in their list of remarkable properties. Graphene is the strongest material ever synthesized, with high values of Young's modulus (≈1 TPa) and ultimate tensile strength (≈130 ± 10 GPa).^22–25 It is also able to withstand tensile strains as large as 25%.²⁶ h-BN also features high values of Young's modulus (≈0.865 TPa) and tensile strength (≈70.5 ± 5.5 GPa).^27–29 Extensive literature exists describing the mechanical properties of graphene and h-BN and suggesting applications. For example, these structures have been used as reinforcing materials in nanocomposites^30,31 and in highly flexible touch screens.³²

The contrasting properties of graphene and h-BN motivated investigations on hybrid 2D materials with intermediate properties. Arrangements of C, B, and N atoms were used to propose graphene-like hexagonal sheets with h-BN nanodomains and vice versa (graphene/h-BN sheets).^33–39 These structures were first synthetized by Ci et al. using a thermal catalytic CVD method.⁴⁰ Since then, other synthesis methods were proposed as well.^41–43 Regarding the mechanical, electronic, and optical properties of these hybrid structures, first-principles calculations provided values intermediate between those found in graphene and h-BN. Interestingly, these values could also be controlled by adjusting the size and geometry of the nanodomain.^33–38 Calculations also showed that the properties of the hybrid structures depend on the graphene/h-BN interface, where the less stable C–B and C–N bonds are found.^33,37,38 Overall, the current results indicate that graphene/h-BN sheets may be useful for applications in optoelectronic devices, as they exhibit variable band gaps (<2 eV) and an optical absorption spectrum in the visible region.

It is also possible to control the physicochemical properties of 2D materials (including graphene/h-BN hybrid sheets) through mechanical deformation.^{22–25,27–29,44–50} One way to achieve this is by applying a tensile strain, with its magnitude and direction determining whether the material undergoes bending, wrinkling, stretching, or breaking.^22,51–54 From an experimental point of view, strain can be applied to 2D structures through various methods.^55–58 For example, the material may be deposited in an elastic substrate, which is then elongated.⁵⁵ As a result of the applied strain, the atomic arrangement of the system is deformed, and many properties of the material are affected by this modification. For instance, changes in thermal conductivity,⁵⁹ electronic and optical properties,^{22–25,27–29,44–48} and chemical activity⁶⁰ have already been reported. In turn, these changes could improve or impair the performance of a 2D material for specific applications. These developments have given rise to a new field of research called “straintronics”, which aims to control electrical properties in 2D materials through mechanical deformation.⁶¹

Additionally, strain exerts a significant effect on the interactions between atomic orbitals in 2D materials, and hence it greatly affects their energy band gap. For example, Fujimoto et al. found that a tensile strain reduces the band gap of h-BN, whereas a small compressive strain increases it.⁶² Other studies found similar results.^22,27 In fact, each 2D material has a particular response to strain. For instance, in silicene a semimetal–metal transition was observed for biaxial strain values greater than 7.5%.⁶³ For transition metal dichalcogenides, calculations and experiments investigated the effects of strain on the electronic properties, revealing a semiconductor–metal transition for large enough tensile strain values.^64,65 And for graphene deposited on suitable substrates,²⁴ strain opens the energy gap.

Strain is also often applied to materials in order to study their mechanical properties. For the hybrid graphene/h-BN sheets, Zhao et al.⁴⁴ used molecular dynamics (MD) simulations to show that the fraction of h-BN contained in the structure determines its mechanical properties. In turn, this fraction depends on the geometry and/or the number of nanodomains in the sheet. They also found Young's modulus values intermediate between those of graphene and h-BN, with lower values for higher h-BN concentrations. More recently, Oliveira et al.⁴⁸ combined DFT calculations and MD simulations to investigate whether the mechanical properties of hybrid sheets depend on their size. The authors found that the mechanical properties are independent of scale, so long as the graphene/h-BN concentration remained constant as the sheet size increased or decreased. Additionally, B_xC_yN_z hexagonal sheets with atomic arrangements that included many B–C and N–C bonds were studied with MD⁴⁹ and DFT⁵⁰ simulations. The reported values of stiffness and tensile strength were lower than those reported for structures composed of h-BN nanodomains embedded in graphene sheets. Regarding the effect of strain on the electronic properties of graphene/h-BN sheets, Azevedo and Kaschny determined that the band gap of this hybrid material increased with strain, from 0.08 eV to 0.5 eV.⁴⁵ However, note that the investigated sheets were composed mainly of carbon atoms. Overall, most studies in the literature concentrate on hybrid structures composed mostly of carbon atoms.^44–50

On the other hand, hybrid structures with high h-BN concentration have received comparatively little attention, even though reports detailing their synthesis exist.^41–43 So far, theoretical investigations on BCN sheets with high h-BN concentration have focused on their electronic and magnetic properties.^38,39 In the present contribution, we combine DFT calculations and MD simulations to investigate the mechanical properties as well as the effect of composition and strain on the electronic properties of h-BN sheets with graphene domains. We find that the band gaps of the investigated systems depend strongly on the graphene/h-BN concentration and moderately or weakly on the external strain (depending on the composition of the hybrid sheet). Furthermore, they can withstand strain values above 10% before permanent deformation. The combination of elasticity with a controllable band gap suggests possible applications of the investigated nanosheets as semiconductor elements on wearable electronics.

II. Computational details and methods

In the present work, first principles calculations and MD simulations were employed to investigate the mechanical and electronic properties of h-BN nanosheets containing graphene nanodomains (h-BN/graphene sheets). Some of the studied structures are shown in Fig. 1 and 2. In our simulations, we considered square h-BN sheets (L_x = L_y = L), with side lengths (L) ranging from 2 to 50 nm (see Tables 1 and 2). In order to create the hybrid structures, we introduced graphene nanodomains in the center of the sheet by replacing B and N atoms with C atoms. Regarding the geometry of the nanodomains, for the L = 2 nm structure we considered hexagonal domains (see Fig. 1), whereas for the larger structures we considered circular, triangular, and star-shaped domains (see Fig. 2). We denote the number of C atoms in the graphene nanodomain by n_C, which ranges from 6 to 3534 atoms. To specify a structure, we need to know both the sheet length and the graphene domain size, and so we identify each structure using the symbol L_j–C_x. For example, L_2nm–C₆ stands for the hybrid sheet illustrated in Fig. 1(a), which has 2 nm in length and a graphene nanodomain with 6 carbon atoms. Additionally, for each structure we calculated the atomic fraction of graphene (γ), which is given by

γ = n_C/n_T, (1)

where n_C follows the previous definition and n_T is the total number of atoms. The values of γ for the investigated structures vary from 0.04 to 0.34, as shown in Table 1.

Fig. 1 Structures investigated using both DFT and MD simulations, composed of hexagonal graphene nanodomains embedded within square h-BN monolayers, with side length L = 2 nm. The values of the atomic fraction of graphene (γ) are (a) 0.04, (b) 0.15, and (c) 0.34. The directions of the applied strain are indicated in (d) and (e).

Fig. 2 Structures investigated using MD simulations, composed of (a) circular, (b) triangular, and (c) star-shaped graphene nanodomains (inscribed within a circle of diameter d) embedded within square h-BN sheets (with side length L). We considered structures for which we kept L constant and varied d (model-I), and also structures for which we kept the ratio d/L constant and varied L (model-II).

Table 1 Calculated mechanical properties of the structures illustrated in Fig. 1(a)–(c). L is the side length, γ is the atomic fraction of graphene, Y is the Young's modulus, σ is the tensile strength, and ε is the ultimate strain. Results for the structures displayed in Fig. 2(a)–(c) are presented in Table 2

DFT calculations	L (nm)	γ	Y (GPa)	Y (GPa)	σ (GPa) [ε]	σ (GPa) [ε]
			Armchair	Zigzag	Armchair	Zigzag
L2nm–C6	2	0.04	771.1	762.8	97.53 [0.22]	81.77 [0.19]
L2nm–C24	2	0.15	781.8	778.9	94.00 [0.19]	85.13 [0.16]
L2nm–C54	2	0.34	818.4	799.4	94.92 [0.19]	85.59 [0.15]

---

MD simulations	L (nm)	γ	Y (GPa)	Y (GPa)	σ (GPa) [ε]	σ (GPa) [ε]
			Armchair	Zigzag	Armchair	Zigzag
L2nm–C6	2	0.04	717.5	712.2	104.5 [0.19]	89.9 [0.17]
L2nm–C24	2	0.15	744.5	729.6	96.5 [0.16]	85.1 [0.16]
L2nm–C54	2	0.34	790.3	767.6	92.9 [0.17]	80.6 [0.15]

---

Table 2 Calculated mechanical properties of the structures illustrated in Fig. 2(a)–(c). L is the side length, γ is the atomic fraction of graphene, Y is the Young's modulus, σ is the tensile strength, and ε is the ultimate strain

Circle	L (nm)	γ	Y (GPa)	Y (GPa)	σ (GPa) [ε]	σ (GPa) [ε]
			Armchair	Zigzag	Armchair	Zigzag
L10nm–C3116	10	0.189	756.5	741.0	87.92 [0.144]	89.14 [0.152]
L20nm–C14316	20	0.193	755.2	740.8	88.66 [0.141]	77.73 [0.125]
L30nm–C27090	30	0.194	753.1	739.7	86.72 [0.139]	81.23 [0.132]
L40nm–C58794	40	0.012	714.6	700.4	86.27 [0.151]	86.04 [0.160]
L40nm–C56624	40	0.049	721.6	707.8	85.56 [0.147]	75.77 [0.131]
L40nm–C53010	40	0.109	733.7	719.5	85.37 [0.142]	79.97 [0.136]
L40nm–C47950	40	0.194	752.4	736.9	87.27 [0.141]	75.67 [0.122]
L40nm–C41444	40	0.304	776.2	763.2	81.16 [0.124]	81.68 [0.128]
L50nm–C74722	50	0.195	751.1	737.7	81.14 [0.127]	81.72 [0.133]

---

Triangle	L (nm)	γ	Y (GPa)	Y (GPa)	σ (GPa) [ε]	σ (GPa) [ε]
			Armchair	Zigzag	Armchair	Zigzag
L10nm–C3534	10	0.079	729.1	715.1	106.55 [0.219]	80.02 [0.141]
L20nm–C13850	20	0.079	727.4	714.6	106.27 [0.209]	82.95 [0.149]
L30nm–C30948	30	0.079	727.4	714.1	95.10 [0.179]	85.81 [0.154]
L40nm–C59214	40	0.005	712.7	698.5	100.11 [0.196]	81.43 [0.149]
L40nm–C58330	40	0.020	715.6	701.9	104.92 [0.211]	82.83 [0.152]
L40nm–C56868	40	0.045	720.9	706.4	101.16 [0.199]	84.53 [0.155]
L40nm–C54759	40	0.079	726.9	712.4	110.29 [0.225]	90.19 [0.166]
L40nm–C52038	40	0.126	736.5	720.9	97.43 [0.180]	85.43 [0.149]
L50nm–C85318	50	0.080	727.6	713.7	96.68 [0.181]	85.57 [0.153]

---

Star	L (nm)	γ	Y (GPa)	Y (GPa)	σ (GPa) [ε]	σ (GPa) [ε]
			Armchair	Zigzag	Armchair	Zigzag
L10nm–C3432	10	0.106	735.5	721.5	93.24 [0.169]	81.01 [0.142]
L20nm–C13454	20	0.105	731.9	719.6	93.24 [0.140]	83.70 [0.146]
L30nm–C30064	30	0.105	733.3	719.4	75.87 [0.122]	79.45 [0.135]
L40nm–C59112	40	0.007	713.5	699.2	77.25 [0.129]	79.34 [0.142]
L40nm–C57934	40	0.026	716.5	702.9	83.49 [0.143]	82.44 [0.151]
L40nm–C55984	40	0.059	723.8	709.3	77.19 [0.127]	78.16 [0.136]
L40nm–C53172	40	0.107	733.0	718.6	81.86 [0.134]	78.65 [0.132]
L40nm–C49544	40	0.168	744.7	729.6	76.07 [0.119]	77.77 [0.127]
L50nm–C82824	50	0.107	733.2	719.8	75.78 [0.121]	77.63 [0.131]

---

All first-principles calculations are based on density functional theory (DFT),⁶⁶ as implemented in the SIESTA code.^67,68 Norm-conserving Troullier–Martins pseudopotentials⁶⁹ were used, in the Kleinman–Bylander factorized form.⁷⁰ A double-ζ polarized basis set (DZP) was used, composed of numerical atomic orbitals of finite range. The exchange–correlation energy is expressed within the generalized gradient approximation (GGA),⁷¹ in the form of the Perdew–Burke–Ernzerhof (PBE) functional. All geometries were fully relaxed, with residual forces smaller than 0.1 eV Å⁻¹. We used a grid cutoff of 150 Ry for the grid integration, with the Brillouin zone sampled using a 10 × 10 × 1 k-point mesh within the Monkhorst and Pack scheme. We adopted a convergence criterion where self-consistency is achieved when the maximum difference between the output and the input of each element of the density matrix is less than 10⁻⁴ eV, in a self-consistent field cycle. We adopted a rectangular unit cell and imposed periodic boundary conditions. A vacuum region of 100 Å was added along the z direction to avoid artificial interactions between a monolayer and its periodic images.

The MD simulations were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code.⁷² The interactions between B, C, and N atoms were described with the Tersoff potential,⁷³ using the parameters adjusted by Kinaci et al.⁷⁴ Our MD simulations were performed using a reasonably small timestep of 0.1 fs, and proceeded in three steps:

(1) We first evolved the system for 2 × 10⁵ steps in the NPT ensemble, using Nose–Hoover thermostats and barostats^75-77 to set temperature and pressure values to 10 K and 0 Pa.

(2) We turned the thermostat off and then evolved the system for 2 × 10⁵ steps in the NPH ensemble, using the same barostat described above to set the pressure to 0 Pa. Note that the barostat was only applied to the planar direction that is not under strain from this step onward.

(3) In the final step we maintained the barostat configuration described above, and then elongated the system for 3 × 10⁶ steps using a strain rate of 10⁻⁶ fs⁻¹ (for a total strain of 30%).

We used the above method to obtain stress–strain curves for graphene and h-BN. The MD results were then compared against experimental results, as detailed in the ESI.‡ Overall, we find very good agreement between our results and the experimental ones for graphene.²⁶ In contrast, for h-BN, we find our results predict a lower Young's modulus and higher tensile strength when compared to the experiments.²³ However, we note that our results are consistent with other theoretical investigations³⁰

Our calculations were performed in two stages. In the first stage, we performed DFT calculations and MD simulations to determine the mechanical and electronic properties of the small hybrid h-BN/graphene sheets (L = 2 nm) that are shown in Fig. 1(a)–(c). In this stage we observed that both methods yielded similar stress–strain curves. So, in the second stage, we only performed MD simulations to investigate the mechanical properties of larger hybrid sheets (10 ≤ L ≤ 50 nm), which are shown in Fig. 2. In turn, we divided the calculations with the larger structures in two models. In model-I, the side length L is constant and the graphene/BN concentration is variable, while in model-II the side length L is variable and the graphene/BN concentration is constant. In both models, we considered different geometries (circular, triangular, and star-shaped). In order to determine the mechanical properties of these structures, we applied a tensile strain along one direction and then calculated the resulting tensile stress to obtain a stress–strain curve. Next, we calculated the slope of the linear region to determine the Young's modulus (Y). The tensile strength (σ) and the ultimate strain (ε) were taken at the point where the stress reaches its maximum value. The calculated values of Y, σ, and ε are presented in Tables 1 and 2.

One final remark is that the methodology used to obtain the stress is different for DFT and MD simulations. For the DFT approach, we first increased one side of the simulation box by 1% and then relaxed both the other side and the atomic positions until the calculation converged. After that, we increased the simulation box by 1% again and repeated the process. For the MD simulations, we increased one direction of the simulation box continuously at a fixed strain rate, and used a barostat to keep the other direction relaxed. Meanwhile, we allowed all atoms to freely evolve. Strain was applied along the armchair (x) and zigzag (y) directions, as illustrated in Fig. 1(d), (e) and 2(a).

III. Results and discussion

We have used DFT calculations and MD simulations to analyze the mechanical properties of the structures displayed in Fig. 1(a)–(c) (L_2nm–C₆, L_2nm–C₂₄, and L_2nm–C₅₄). For each structure, we obtained stress–strain curves for both the armchair and zigzag directions which are shown in Fig. 3(a)–(f). In Fig. 3 we highlight in green the linear elastic region, which occurs for strain values between 0 and 0.05 (or 5%). Observe the good agreement between the DFT and MD results in this region. This agreement is likely related to the fact that the parameters used for the MD simulations were parameterized to reproduce the results of DFT calculations.⁷⁴ As a consequence, the difference between Young's modulus values obtained through DFT and MD is only ∼10% (see Fig. 4(a)). We found a non-linear relationship between stress and strain for higher strain values, even though no atomic rearrangement occurred, which would indicate permanent deformation. Also, notice that the agreement between the DFT and MD results is not as good in this non-linear region, although the disparity is still reasonably small. We attribute these differences to distinct DFT and MD methodologies (for instance, we use non-zero temperatures in the MD simulations). As the strain continues to increase, eventually the material fractures for strain values between 15% and 22%. Both DFT and MD results predict that the maximum tensile strain is higher along the armchair direction than along the zigzag direction. The stress–strain curves obtained here are similar to other curves reported in the literature.^44–50

Fig. 3 Stress–strain curves obtained through DFT and MD simulations for the structures shown in Fig. 1. In plots (a), (b), and (c) strain is applied to the armchair direction; while in plots (d), (e), and (f) strain is applied to the zigzag direction. For DFT results, the red solid circles indicate data points and the line is only a guide to the eye. For MD results, the number of data points is large, and a blue line is used to connect adjacent points. The green region corresponds to the elastic region.

Fig. 4 Plot (a) shows Young's modulus results against graphene concentration for the structures displayed in Fig. 1. Plots (b) and (c) show Young's modulus results against graphene concentration and side length, respectively, for the structures shown in Fig. 2.

Let us discuss in more detail the DFT predictions for the mechanical properties of the h-BN nanosheets containing graphene domains. The Young's modulus (Y) and tensile strength (σ) values obtained for L_2nm–C₆, L_2nm–C₂₄, and L_2nm–C₅₄ are summarized in Table 1. Comparing the results obtained for the zigzag and armchair directions, we find that the mechanical properties of the hybrid sheets are weakly anisotropic, with differences of ∼2% for Y and ∼10% for σ. Note that the same behavior was observed for h-BN, graphene, and other h-BN/graphene sheets, since this weak anisotropy is related to the hexagonal structure of the unit cell.^{22–25,27–29,44–50} Our results also show that the Young's modulus values increase with increasing concentrations of graphene, regardless of the direction of applied strain. In contrast, we find no conclusive relation between the tensile strength and ultimate strain with the number of C atoms. These results are in agreement with others previously reported for graphene/BN sheets, which found higher Young's modulus for higher concentrations of graphene.^44,45,48 After the stress reaches its maximum value, for strain values of 21% for the armchair and 15% for the zigzag direction, the h-BN/graphene sheets fracture (see Fig. 5). Regarding the obtained fracture patterns, we find similar results for the zigzag and armchair directions. Bonds initially break in the vicinity of the nanodomain, and then the fracture propagates along a direction that is perpendicular to the applied strain. In this process, C–B, C–N, and B–N bonds are broken, whereas the graphene domains remain intact. After fracture occurs, the DFT calculations no longer converge.

Fig. 5 Snapshots from DFT calculations displaying the fracture patterns obtained for L_2nm–C₆, L_2nm–C₂₄, and L_2nm–C₅₄. The arrows in this figure indicate the directions in which we applied strain.

Let us now discuss the mechanical properties of the larger hybrid sheets (10 ≤ L ≤ 50 nm), which are shown in Fig. 2. Due to their large size, only MD simulations were used to investigate their mechanical properties. Note, however, that the comparison between DFT and MD results for the smaller structures indicates that the MD results are reliable, particularly in the linear elastic region. We consider two models (model-I and model-II) for the larger graphene/BN sheets. In model-I we keep the size of the hybrid sheet constant (L = 40 nm) while we vary the domain size (γ ranges from 0.07 to 0.304, see Table 2). The goal of model-I is to investigate the dependence of the mechanical properties on composition. In model-II we keep the fraction of graphene constant (γ = 0.19) while we vary the sheet size (10 ≤ L ≤ 50 nm). The goal of model-II is to investigate possible size effects (the ratio of hybrid C–B and C–N bonds to the total number of bonds decreases as the sheet size increases). We also consider whether the shape of the domain affects the mechanical properties for both model-I and model-II. The mechanical properties of the larger sheets are summarized in Table 2. These results were extracted from various stress–strain curves, which can be found in Fig. S2 of the ESI.‡ Our results indicate that the linear region for the larger structures extends from 0 to 3%.

Regarding model-I, we find that the Young's modulus increases as the concentration of graphene increases, whereas we find no clear relationship between tensile strength and composition. Young's modulus results for model-I are summarized in Fig. 4(b). Observe that Y increases linearly with the concentration of graphene for both armchair and zigzag directions. Comparing results for both directions, we again notice that the Young's modulus is slightly higher for the armchair direction. Regarding the domain shape, we find that it has no noticeable effect on Young's modulus values, other than the fact that different shapes have distinct concentrations of graphene. This conclusion is supported by the data presented in Fig. 4(b). Notice that all Young's modulus results obtained fall within the same line, irrespective of the domain shape. We discuss possible reasons for the lack of relationship between Young's modulus values and domain shape in the ESI.‡

Moving on to the tensile strength results for model-I, they are summarized in Table 2. Observe that σ values vary as the domain size increases, but with no clear trend, sometimes increasing and sometimes decreasing as the graphene concentration increases. This occurs because the tensile strength of the hybrid sheets is limited by the weaker B–C and C–N bonds (the fracture process is described ahead). Finally, note that, unlike the Young's modulus, the tensile strength has a weak dependence on the domain shape, being overall higher for triangular shapes when comparing structures with similar γ values. Furthermore, while we did not study the effect of the shape on the electronic properties, previous investigations on this topic found that shape affects these properties.³³

Regarding model-II, we find that the mechanical properties of the BN nanosheets with graphene domains do not depend on the sheet size, as long as the graphene concentration remains constant as the sheet size varies. Young's modulus results for model-II are presented in Fig. 4(c). For constant γ, notice that the Young's modulus remains constant as the size of the nanosheet increases, for all pulling directions and domain shapes. Inspection of Fig. 4(c) also reveals that the Young's modulus is yet again slightly higher along the armchair direction. Also, notice that Y is higher for circular and lower for triangular domains. However, it is important to remark that although γ is constant for structures with the same domain shape, it does vary for structures with different domain shapes. For model-II, we used γ values of 0.19, 0.08, and 0.11 for circular, triangular, and star-shaped domains, respectively. Hence, model-II structures with circular domains are stiffer simply because they contain more graphene. With regard to the tensile strength, we find no clear relationship between this quantity and the nanosheet size in the results presented in Table 2. Finally, the tensile strength appears to be higher for structures with triangular domains, particularly when strain is applied along the armchair direction.

We have also considered systems with lower symmetry than those discussed so far. We used MD simulations to study h-BN monolayers with (i) asymmetrical graphene domains and (ii) rotated triangular graphene domains. The results obtained for structures with lower symmetry are presented and discussed in the ESI.‡ In general, we continue to find Young's modulus values intermediate between h-BN and graphene. However, we found a higher degree of anisotropy for the mechanical properties of the asymmetrical structure. We also found that the orientation of the triangular domain does not affect Young's modulus results but does affect tensile strength and ultimate strain values.

Fig. 6 shows the stress and fracture patterns obtained in an MD simulation where the strain was applied along the armchair direction. Fig. 6(a) shows atomic configurations while Fig. 6(b) shows the corresponding stress distributions. Similar results for the zigzag direction are presented in Fig. S4 of the ESI.‡ Since the overall behavior is the same for both directions, we only detail armchair results. Before fracture the stress is higher (blue) in the middle region of the hybrid sheet, whereas the stress is lower (white) in the regions above and below the graphene domain. To understand this result, notice that if we think of the sheet as composed of multiple thin segments spanning its length along the x-direction, we find that each thin segment is under the same strain (as the entire sheet is under the same strain). Hence, as the regions above and below the graphene domain are composed entirely of the less stiff h-BN, lower stress values produce the same strain. Consequently, fracture starts in the middle region, where stress is higher on average (t = 0 ps). However, note that strain is not particularly high at the graphene–BN interface, where the fracture always begins. This region is where we find the weaker hybrid bonds, and the C–B bonds were always the first to break. To understand this result, we note a previous DFT study that calculated first-neighbor bond energies for C–C, B–N, C–N, and C–B bonds, for hybrid B_xC_yN_z structures.⁷⁸ The authors found that the C–B bond had the lowest bond energy. Fracture then propagates perpendicular to the applied strain (t = 0.2 ps). This result is in agreement with the previously discussed DFT results. As the process continues, stress decreases (t = 0.4 ps).

Fig. 6 Fracture results, for strain applied along the armchair direction for the hybrid sheet with L_x = L_y = 10 nm and a circular graphene domain (d = 5 nm). (a) Snapshots from MD simulations detailing the time evolution of the monolayer fracture. (b) Corresponding stress distribution for the structures presented in (a).

Let us now discuss the DFT predictions for the electronic properties of the structures shown in Fig. 1(a)–(c). Fig. 7(a) and (b) present the energy band gap E_g against the strain for the armchair and zigzag directions. These values were calculated from the electronic band structures with 100 k-points along the Γ–X–Z direction. For zero strain, each structure has the same value of E_g in both directions, as expected. The gap values for L_2nm–C₆, L_2nm–C₂₄, and L_2nm–C₅₄ are 3.52 eV, 2.27 eV, and 1.56 eV, respectively. Comparing the smaller graphene domain with the others, we find a decline in the band gap of 35.5% for L_2nm–C₂₄ and 55.7% for L_2nm–C₅₄. This energy gap reduction for increased concentrations of graphene is in agreement with the results found by Manna et al.³³ and Azevedo et al.³⁷ Fig. 8(a) shows the projected density of states (PDOS) for selected h-BN nanosheets containing graphene nanodomains. We can infer from the PDOS that the electronic states close to the Fermi energy E_f are in general associated with C atoms. On the other hand, B and N atoms contribute to the electronic states in the conduction and valence bands, respectively. In agreement with the PDOS, the localized density of states (LDOS) results presented in Fig. 8(b) reveal that the bottom of the conduction band and the top of the valence band are associated with p_z orbitals from C atoms. These results are common features of hybrid nanostructures with stoichiometry B_xC_yN_z.^12,36,37 Finally, note that we did not calculate the band gap of the larger structures presented in Fig. 2 because they were only investigated using classical molecular dynamics methods, which cannot be used to calculate electronic properties.

Fig. 7 Results of energy band gap versus strain for the structures presented in Fig. 1. (a and b) Show results for the case where the strain was applied to the armchair and zigzag directions, respectively.

Fig. 8 Plots (a) and (b) show the projected density of states (PDOS) for structure L_2nm–C₂₄ for different strain values. In (a), the strain was applied to the armchair direction and in (b) to the zigzag direction. The dotted vertical line indicates the Fermi energy E_f. (c) Shows the local density of states (LDOS) associated with the bottom of the conduction band and the top of the valence band for the L_2nm–C₂₄ monolayer without strain. (d and e) Show the LDOS of L_2nm–C₂₄ for strains applied to the armchair and zigzag directions, respectively.

In general, we observe that the energy gap decreases gradually as the applied strain increases, except in the case of L_2nm–C₆ for strain values >16%, where we observed a sharp gap decline. The following discussion on the band gap variation with strain is concentrated on results for values below 15%, since fracture becomes an issue for some structures for larger strain values. In general, the energy gap decreases slightly more when the strain is applied to the armchair direction. Furthermore, we observe that the band gap of the considered structures varies weakly or moderately with deformation. These results can be seen more easily in Table 3, which provides band gap values and percentage variations for selected strain values. Examining these results, we notice that the gap values become more sensitive to the external strain as we increase the fraction of graphene in the hybrid sheet. The relatively small gap variation for structures with a high concentration of h-BN is expected, since the band gap of this material presents small variation with the strain for a uniaxial strain below 15%.⁷⁹ In contrast, the band gap of other 2D materials can vary by more than 50% for uniaxial strain values around 10%. For example, this is the case for many transition metal dichalcogenides⁸⁰ (the smallest variation was observed for MoTe₂, for which the band gap decreased around 35% for 10% strain). Finally, when comparing the effects of composition and strain, we find that the former has a more substantial influence on the band gap (except for the case of L_2nm–C₆ deformed along the zigzag direction). We can also combine both effects to decrease the band gap from 3.5 eV to values below 1.0 eV.

Table 3 Electronic band gap (E_g) of the structures illustrated in Fig. 2(a)–(c) for selected strain values. The numbers in parenthesis correspond to the percentage variation relative to the initial value

DFT calculations	Armchair 5%	Armchair 10%	Zigzag 5%	Zigzag 10%
	Eg (eV)	Eg (eV)	Eg (eV)	Eg (eV)
L2nm–C6	3.25 (7.67%)	3.00 (14.8%)	3.28 (6.82%)	3.02 (14.2%)
L2nm–C24	2.00 (11.9%)	1.72 (24.2%)	2.01 (11.5%)	1.76 (22.5%)
L2nm–C54	1.32 (15.4%)	1.09 (30.1%)	1.35 (13.5%)	1.12 (28.2%)

---

An important note is that the DFT method underestimates the band gap of semiconductors,⁸¹ and consequently, the numbers presented here are underestimated. However, we remark that other articles calculated band gaps of strained structures using different methods and found that, while the band gap values depend on the method employed, the general trends do not.^79,82 The size of the structures considered here in the DFT calculations (160 atoms) hampered the use of more accurate methods.

Altogether, the calculated properties of the h-BN nanosheets with graphene domains indicate possible uses in wearable electronics. For such applications, a material needs to (i) be stretchable, (ii) be flexible, (iii) operate under large deformation.⁸³ Regarding requirement (i), our calculations show that the investigated monolayers can be extended by at least 10% before undergoing permanent deformation. Concerning condition (ii), although the hybrid sheets are very stiff, they are also very thin, requiring small forces to deform (for example, flexible silicon devices have been produced by reducing their thickness⁸⁴). Finally, regarding requirement (iii), our results show that the considered structures remain semiconductors over a wide range of strain values. They also reveal that it is possible to control the band gap of the considered monolayers by changing their composition. Semiconductors have been used in wearable electronics as photodetectors,⁸³ as transistors,⁸⁵ and as components in sensors.⁸⁶ Therefore, h-BN nanosheets with graphene domains may play an important role in wearable electronics.

IV. Conclusions

In summary, we combined DFT and MD simulations to investigate the mechanical and electronic properties of h-BN nanosheets with graphene domains under strain. Regarding the mechanical properties, we found good agreement between the DFT calculations and MD simulations – the former was used to study smaller systems and the latter to study both smaller and larger structures. We found that Young's modulus values depend only on the relative concentration of graphene and h-BN. It does not depend on either the shape of the domain or the size of the structure (provided the h-BN and graphene domain are increased/decreased proportionately). Our results also revealed that the tensile strength is slightly higher for triangular graphene domains, particularly when applying strain to the armchair direction. In the case of DFT calculations, the system fractures for strain values between 15% and 22% at the junction between h-BN and graphene. In the case of MD calculations, the ultimate strain was slightly lower, but fracture also starts in the same region, even though stress is not higher at the interface. Fracture begins at the junction between h-BN and graphene because it is the location of the weaker C–B bonds.

Regarding the electronic properties, we found that the band gap depends more on the composition of the hybrid sheets than on the external strain. Our results also show that the band gap of structures with more carbon atoms is more sensitive to the applied strain. However, the variation of the band gap with the external strain is still relatively small compared to other 2D materials. Furthermore, when composition and strain are combined to modify the band gap, we can obtain a broad range of values (from 1.0 eV to 3.5 eV). The proposed hybrid monolayers combine elasticity (the nanosheets can withstand considerable stress before plastic deformation), atomic-scale thickness (allowing the monolayers to bend and deform easily), and adjustable band gap (allowing the customization of the material for different applications). Those properties indicate that the h-BN nanosheets with graphene domains could be promising materials for applications in wearable electronics.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We would like to thank the Brazilian Research Agencies CNPq and CAPES for financial support. LDM acknowledges the computational support provided by the High Performance Computing Center at UFRN (NPAD/UFRN). SA acknowledges the PRONEX – FAPESQ/PB – MCT/CNPq Grants no. 006/2018 and 151/2018.

References

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science, 2004, 306, 666.
2. A. K. Geim, Science, 2009, 324, 1530.
3. M. J. Allen, V. C. Tung and R. B. Kaner, Chem. Rev., 2010, 110, 132.
4. A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov and A. K. Geim, Rev. Mod. Phys., 2009, 81, 109.
5. S. D. Sarma, A. Shaffique, E. H. Hwang and R. Enrico, Rev. Mod. Phys., 2011, 83, 407.
6. X. Xu, L. F. C. Pereira, Y. Wang, J. Wu, K. Zhang, X. Zhao, S. Bae, C. T. Bui, R. Xie, J. T. L. Thong, B. H. Hong, K. P. Loh, D. Donadio, B. Li and B. Ozyilmaz, Nat. Commun., 2014, 5, 3689.
7. C. Guo, J. Zhang, W. Xu, K. Liu, X. Yuan, S. Qin and Z. Zhu, Nanomaterials, 2018, 8, 1033.
8. L. H. Li and Y. Chen, Adv. Funct. Mater., 2016, 26, 2594.
9. A. Nag, K. Raidongia, K. P. S. S. Hembram, R. Datta, U. V. Waghmare and C. N. R. Rao, ACS Nano, 2010, 4, 1539.
10. Silicene: Structure, Properties and Applications, ed. M. J. Spencer and T. Morishita, Springer, New York, 2016.
11. B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti and A. Kis, Nat. Nanotechnol., 2011, 6, 147.
12. S. Azevedo, J. Kaschny, M. Caio and F. de Brito Mota, Eur. Phys. J. B, 2009, 67, 507.
13. J. Eichler and C. Lesniak, J. Eur. Ceram. Soc., 2008, 28, 1105.
14. W. Han, L. Wu, Y. Zhu, K. Watanabe and T. Taniguchi, Appl. Phys. Lett., 2008, 93, 223103.
15. Y. Kubota, K. Watanabe, O. Tsuda and T. Taniguchi, Science, 2007, 317, 932.
16. R. Beiranvand and S. Valedbagi, Diamond Relat. Mater., 2015, 58, 190.
17. G. Kim, A. Jang, H. Y. Jeong, Z. Lee, D. J. Kang and H. Suk Shin, Nano Lett., 2013, 13, 1834.
18. K. K. Kim, A. Hsu, X. Jia, S. M. Kim, Y. Shi, M. Hofmann, D. Nezich, J. F. Rodriguez-Nieva, M. Dresselhaus, T. Palacios and J. Kong, Nano Lett., 2012, 12, 161.
19. C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard and J. Hone, Nat. Nanotechnol., 2010, 5, 722.
20. Z. Pawlak, T. Kaldonski, R. Pai, E. Bayraktar and A. Oloyede, Wear, 2009, 267, 1198.
21. Y. Sasama, K. Komatsu, S. Moriyama, M. Imura, T. Teraji, K. Watanabe, T. Taniguchi, T. Uchihashi and Y. Takahide, APL Mater., 2018, 6, 111105.
22. S. Deng, A. V. Sumant and V. Berry, Nano Today, 2018, 22, 14.
23. A. Falin, Q. Cai, E. J. G. Santos, D. Scullion, D. Qian, R. Zhang, Z. Yang, S. Huang, K. Watanabe and T. Taniguchi, et al., Nat. Commun., 2017, 8, 15815.
24. C. Si, Z. Sun and F. Liu, Nanoscale, 2016, 8, 3207.
25. S. Choi, S. Jhi and Y. Son, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 081407.
26. C. Lee, X. Wei, J. W. Kysar and J. Hone, Science, 2008, 321, 385.
27. Q. Peng, Nanotechnology, 2018, 29, 405201.
28. Q. Peng, W. Ji and S. De, Comput. Mater. Sci., 2012, 56, 11.
29. A. Falin, Q. Cai, E. J. G. Santos, D. Scullion, D. Qian, R. Zhang, Z. Yang, S. Huang, K. Watanabe, T. Taniguchi, M. R. Barnett, Y. Chen, R. S. Ruoff and L. H. Li, Nat. Commun., 2017, 8, 15815.
30. D. Akinwande, et al., Extreme Mechanics Letters, 2017, 13, 42.
31. J. R. Young, I. A. Kinloch, L. Gong and K. S. Novoselov, Compos. Sci. Technol., 2012, 72, 1459.
32. A. I. Vlasov, D. S. Terent’ev and V. A. Shakhnov, Russ. Microelectron., 2017, 46, 192.
33. A. K. Manna and S. K. Pati, J. Phys. Chem. C, 2011, 115, 10842.
34. B. Xu, Y. H. Lu, Y. P. Feng and J. Y. Lin, J. Appl. Phys., 2010, 108, 073711.
35. X. Fan, Z. Shen, A. Q. Liuc and J. Kuod, Nanoscale, 2012, 4, 2157.
36. L. Leite, S. Azevedo and B. L. Bernardo, Superlattices Microstruct., 2017, 111, 1137.
37. S. Azevedo, F. Moraes and B. L. Bernardo, Appl. Phys. A, 2014, 117, 2095.
38. A. C. Bevilacqua, M. H. Kohler, S. Azevedo and R. J. Baierle, Phys. Chem. Chem. Phys., 2017, 19, 5629.
39. R. M. Santos, A. L. Aguiar, L. A. R. Junior and J. R. Martins, Chem. Phys. Lett., 2019, 724, 103.
40. L. Ci, L. Song, C. Jin, D. Jariwala, D. Wu, Y. Li, A. Srivastava, Z. F. Wang, K. Storr, L. Balicas, F. Liu and P. M. Ajayan, Nat. Mater., 2010, 9, 430.
41. Z. Liu, L. Ma, G. Shi, W. Zhou, Y. Gong, S. Lei, X. Yang, J. Zhang, J. Yu, K. P. Hackenberg, A. Babakhani, J. Idrobo, R. Vajtai, J. Lou and P. M. Ajayan, Nat. Nanotechnol., 2013, 8, 119.
42. C. Huang, C. Chen, M. Zhang, L. Lin, X. Ye, S. Lin, M. Antonietti and X. Wang, Nat. Commun., 2015, 6, 7698.
43. L. Cheng, J. Meng, X. Pan, Y. Lu, X. Zhang, M. Gao, Z. Yin, D. Wang, Y. Wang, J. You, J. Zhang and E. Xie, Nanoscale, 2019, 11, 10454.
44. S. Zhao and J. Xue, J. Phys. D: Appl. Phys., 2013, 46, 135303.
45. S. Azevedo and J. R. Kaschny, Solid State Commun., 2019, 290, 60.
46. Y. Zhang, Q. Pei, Z. Sha and Y. Zhang, Phys. Lett. A, 2019, 383, 2821.
47. L. Jiao, M. Hu, Y. Peng, Y. Luo, C. Li and Z. Chen, J. Solid State Chem., 2016, 244, 120.
48. I. S. Oliveira, J. S. Lima, A. Freitas, C. G. Bezerra, S. Azevedo and L. D. Machado, RSC Adv., 2021, 11, 12595.
49. Y. Zhang, Q. Pei, Z. Sha and Y. Zhang, Phys. Lett. A, 2019, 383, 2821.
50. L. Jiao, M. Hu, Y. Peng, Y. Luo, C. Li and Z. Chen, J. Solid State Chem., 2016, 244, 120.
51. S. Deng, A. V. Sumant and V. Berry, Nano Today, 2018, 22, 14.
52. C. Si, Z. Sun and F. Liu, Nanoscale, 2016, 8, 3207.
53. S. Choi, S. Jhi and Y. Son, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 081407R.
54. Z. H. Ni, T. Yu, Y. H. Lu, Y. Y. Wang, Y. P. Feng and Z. X. Shen, ACS Nano, 2008, 2, 2301.
55. C. Lee, X. Wei, J. W. Kysar and J. Hone, Science, 2008, 321, 385.
56. M. Huang, H. Yan, C. Chen, D. Song, T. F. Heinz and J. Hone, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 7304.
57. T. Yu, Z. Ni, C. Du, Y. You, Y. Wang and Z. Shen, J. Phys. Chem. C, 2008, 112, 12602.
58. F. Ding, H. Ji, Y. Chen, A. Herklotz, K. Dorr, Y. Mei, A. Rastelli and O. G. Schmidt, Nano Lett., 2010, 10, 3453.
59. M. Guo, Y. Qian, H. Qi, B. Kedong and Y. Chen, Carbon, 2020, 157, 185.
60. P. L. de Andres and J. A. Verges, Appl. Phys. Lett., 2008, 93, 171915.
61. A. A. Bukharaev, A. K. Zvezdin, A. P. Pyatakov and Y. K. Fetisov, Phys.-Usp., 2018, 61, 1175.
62. Y. Fujimoto, T. Koretsune and S. Saito, J. Ceram. Soc. Jpn., 2014, 122, 346.
63. G. Liu, M. S. Wu, C. Y. Ouyang and B. Xu, Europhys. Lett., 2012, 99, 17010.
64. T. Cheiwchanchamnangij and W. R. L. Lambrecht, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 205302.
65. J. Liang, J. Zhang, Z. Li, H. Hong, J. Wang and Z. Zhang, et al., Nano Lett., 2017, 17, 7539.
66. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133.
67. D. Sánchez-Portal, P. Ordejón, E. Artacho and J. M. Soler, Int. J. Quantum Chem., 1997, 65, 453.
68. J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón and D. Sánchez-Portal, J. Phys.: Condens. Matter, 2002, 14, 2745.
69. N. Troullier and J. L. Martins, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 1993.
70. L. Kleinman and D. M. Bylander, Phys. Rev. Lett., 1982, 48, 1425.
71. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
72. S. Plimpton, J. Comput. Phys., 1995, 117, 1.
73. J. Tersoff, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 39, 5566(R).
74. A. Kinaci, J. B. Haskins, C. Sevik and T. Cagin, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 115410.
75. G. J. Martyna, D. J. Tobias and M. L. Klein, J. Phys. Chem., 1994, 101, 4177.
76. W. Shinoda, M. Shiga and M. Mikami, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 134103.
77. M. Tuckerman, J. Alejandre, R. Lopez-Rendon, A. L. Jochim and G. J. Martyna, J. Phys. A: Math. Gen., 2006, 39, 5629.
78. M. S. C. Mazzoni, R. W. Nunes, S. Azevedo and H. Chacham, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 073108.
79. J. Wu, B. Wang, Y. Wei, R. Yang and M. Dresselhaus, Mater. Res. Lett., 2013, 1, 200.
80. P. Johari and V. B. Shenoy, ACS Nano, 2012, 6, 5449.
81. J. P. Perdew, Int. J. Quantum Chem., 1985, 28, 497.
82. D. Cakır, H. Sahin and F. M. Peeters, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 205421.
83. P. C. Y. Chow and T. Someya, Adv. Mater., 2020, 32, 1902045.
84. K. Dae-Hyeong, et al., Science, 2011, 333, 838.
85. S. Wang, et al., Nature, 2018, 555, 83.
86. M. Barbaro, A. Caboni, P. Cosseddu, G. Mattana and A. Bonfiglio, IEEE Trans. Inf. Technol. Biomed., 2010, 14, 758.

---

Footnotes

† PACS numbers: 73.20.At, 73.20.Hb and 71.15.Mb.

‡ Electronic supplementary information (ESI) available. See DOI: 10.1039/d1ra05831b

---