Atomistic insights into the anisotropic mechanical properties and role of ripples on the thermal expansion of h-BCN monolayers

Monolayer boron–carbon–nitrogen (h-BCN) has been studied in comparison with graphene and hexagonal boron nitride (h-BN) using classical molecular dynamics (MD) simulations with an aim to better understand the structural and thermal behaviors and the anisotropic mechanical properties. The structural features of the simulated sample were analyzed using the pair-correlation function and a full width at half maximum (FWHM). As a hetero-structure of h-BN and graphene, the C–C bond in the h-BCN is responsible for an improved FWHM compared to graphene. Consistent with graphene and h-BN, the in-plane lattice parameter of h-BCN shows thermal contraction over a wide range of temperatures and exhibits a system size dependence. The observed thermal contraction is explained by the presence of out-of-plane bending modes excited at finite temperatures. A tensile test has been performed as a suitable means of measuring the mechanical properties of the h-BCN sheet for zigzag and armchair orientations and found that it is mechanically anisotropic and stable under various strain directions and temperatures. The fracture strength of h-BCN is affected by loading direction and temperature. We found that the Young's modulus of h-BCN is smaller than that of graphene but is higher than that of an h-BN monolayer, suggesting that h-BCN has high mechanical stiffness. Our modeling-based findings provide a guide for future experiments concerning the physical properties of this advanced composite material.


Introduction
Due to the success of graphene, 1,2 alternative layered and nonlayered two-dimensional (2D) materials have become a center of interest for potent research due to their unique chemical and physical properties. Graphene-analogous materials have a wide range of applications for electrochemical storage devices and optoelectronic devices due to their unusual features. This 2D materials family has unique anisotropic mechanical, optical, electrical, and thermal properties which create promising opportunities for many applications, including the design of novel devices 3-6 and aerospace components. Hexagonal boron nitride (h-BN) is one of the well-studied graphene-analogous materials. 7,8 The atomic lattice of h-BN consists of a honeycomb structure with sp 2 hybridization and possesses high bond strength. Even though the physical properties of h-BN are comparable to graphene, it is a wide-bandgap insulator with an electronic band-gap of 5.97 eV, 7 whereas graphene is a zeroband-gap material. However, the application of graphene in nanoelectronics is limited because it does not exhibit a bandgap and is typically categorized as a semi-metallic material. The very high bandgap with insulating characteristics of h-BN also restrict its application in optoelectronics or nanoelectronics.
Inspired by the outstanding features of graphene and h-BN, efforts have been made to study additional hexagonal structures composed of nitrogen, boron, and carbon atoms. Most recently, Beniwal et al. reported the successful experimental synthesis of a 2D graphenic with a ternary monolayer containing boron, carbon, and nitrogen atoms; they named this material "monolayer boron-carbon-nitrogen" (h-BCN). 9 By using density functional theory (DFT) calculations, they predicted that the direct electronic band gap of h-BCN is 1.5 eV and is midway between that of gapless graphene and that of insulating h-BN. 9 In addition to the successful synthesis and characterization of h-BCN, the possible synthesis of other graphitic carbon nitride materials has been reported recently. 10,11 By employing rst-principles DFT calculations, Jiao et al. explored the structural stability, electronic properties, mechanical stability, and optical properties of different types of monolayer BC 2 N. 12 Mortazavi et al. performed extensive MD simulations on graphene like C 3 N and reported superior mechanical properties and thermal conductivity comparable to graphene. 13 Using MD simulations, Lin et al. studied the anisotropic thermal transport properties in monolayer BC 2 N. 14 Most recently, Zhang et al. examined the thermal conductivity in the h-BCN monolayer using MD simulations. They reported a signicantly lower thermal conductivity for h-BCN in comparison with graphene and h-BN. 15 Using classical MD simulations, we performed original studies on h-BCN to determine the following: the structural stability, the effect of ripples on the thermal expansion properties, and the anisotropic mechanical stability as a function of the temperature. To the best of our knowledge, there have been no previous theoretical or experimental investigations of these properties. Therefore, predictions based on our computation results can serve as a guide to researchers interested in the various properties of this novel 2D material.

Computational details
Various properties of the h-BCN monolayer are studied using classical MD simulations that performed using the LAMMPS 16 package. Simulation boxes with a total number of 9600 atoms were used in most of the MD simulation. In this work, a Tersoff type empirical inter-atomic potential reported in Kinaci et al. 17 parameterised by Sevik et al. 18 and Lindsay et al. 19 has been used to model the interaction between the boron-carbon-nitrogen atoms in the h-BCN monolayer. We have studied the properties of the system by varying temperatures between 100 K and 1000 K. All the simulations were initiated by relaxing the atomistic conguration of the h-BCN system to a minimum potential energy state. In the rst phase of the simulations, we relaxed the h-BCN nanosheet under isothermal-isobaric ensemble (NPT) ensemble for 200 picoseconds (ps) with a time step of 0.001 ps. The simulation temperature and pressure were controlled using a Nose-Hoover thermostat algorithm and a barostat algorithm, respectively. We applied periodic boundary conditions in all three dimensions, and we kept the direction of thickness along the Z-direction as high to ensure that only a single sheet of h-BCN was involved in the simulation procedure. During the MD study of the h-BCN structure, we solved the equations of motion via the Velocity-Verlet scheme. Aer an NPT run of 100,000 steps (a time range of 100 ps), we performed an MD run with the canonical ensemble (NVT) for another 100 000 steps to further equilibrate the system at each temperature.
During the stress-strain calculations, in order to avoid spurious behavior when predicting the mechanical properties using optimized Tersoff interatomic potential, we modied the lower cut-off value of the potential to 2.10Å. Aer ensuring proper equilibration and relaxation using NPT ensemble for 200 ps, we applied uniaxial tension by pulling the sheet along the xdirection (zigzag) or along y-direction (armchair) at a strain rate of 0.001 per picoseconds, while in the direction perpendicular to the loading direction, the traction force was free. To maintain the desired boundary conditions, we enforced the NPT along all directions. We performed the simulations for 500 picoseconds. Stress calculations were performed with virial formulation using the following equation: where i and j denote indices in the Cartesian coordinate system; a and b are the atomic indices; m a and v a are the mass and velocity of atom a; r ab is the distance between atoms a and b, f is the atomistic volume of the system, and f ab i is the force along direction i on atom a due to atom b.

Results and discussion
A Structural stability First, we examined the possible h-BCN monolayer stoichiometry as reported in the experimental work of Beniwal et al. 9 In our atomistic studies herein, we also noticed that h-BCN has three possible structures (h-BCN_v1, h-BCN_v2 and h-BCN_v3) with different coordination environments around the B, N, and C atoms of the networks within the same B 2 C 2 N 2 stoichiometry. Our cohesive energy calculation shows that the h-BCN_v3 model possesses the highest structural stability, followed by h-BCN_v1 and h-BCN_v2, as shown in Fig. 1. The three h-BCN monolayers are composed of B 2 C 2 N 2 units, but the orientation and arrangement of B 2 C 2 N 2 units are different for each model. Therefore, they have different connectivity between the B 2 C 2 N 2 units and different types of hexagonal units: B 2 NC 3 , BN 2 C 3 , B 2 C 2 N 2 , B 3 N 3 , N 3 C 3 , and B 3 C 3 . In contrast, while the h-BCN_v1 and h-BCN_v3 models have C-C and B-N connecting bonds between the B 2 C 2 N 2 units, the h-BCN_v2 model has a B-C connecting bond. Comparing the bond length of the connecting bonds, we found that the B-C (1.52Å) bond is longer than both the C-C (1.42Å) and B-N (1.45Å) bonds. The weak B-C connecting bond makes the h-BCN_v2 model unstable. In addition, h-BNC_v2 has a B 3 C 3 unit composed of only a weak B-C bond.
Although both h-BCN_v1 and h-BCN_v3 models have the same C-C and B-N connecting bonds, the different arrangement of their B 2 C 2 N 2 units leads to varying contributions of C-C and B-N bonds to the material's structural stability. The h-BCN_v1 model is more stable than the h-BCN_v2 model due to the greater contribution of the C-C bond. Since the h-BCN_v3 model possesses the most energetic stability of the three models, we focus the remaining discussion on the MD simulation with the h-BCN_v3 conguration. Hereaer, h-BCN_v3 is denoted simply as h-BCN. Fig. 2a shows a representative structure of the h-BCN monolayer. Analysis of the structural properties of h-BCN was carried out at various temperatures. The radial distribution function (RDF) (or equivalently, the pair correlation function), is a measure that can be used to forecast the stability of the atomistic structure. To evaluate the structural stability of h-BCN nanosheet, the RDF was used. The RDF is a Dirac delta function at zero kelvin temperature due to the fact that there exists a well-dened separation between atom pairs in any perfect crystal. But the delta functions broaden into smooth peaks at high temperatures due to thermal vibrations. The width of the RDF peaks is proportional to the root mean squared displacement of the atoms from their equilibrium position. As the temperature increases, the height of peaks in the RDF decreases due to thermal broadening. The RDFs of h-BCN for all the possible atomic distances at 300 K are shown in Fig. 2b and c. During MD simulation, the full width half maximum (FWHM) of RDF is considered to be intimately related to the integrity of the atomistic structure. We predicted that the structures with sharp peaks and relatively smaller values of FWHM for RDF would possess better structural integrity. We also studied the stability and structural integrity of h-BCN nanosheet at different temperatures with the help of the FWHM values estimated using the RDF, as shown in Fig. S1. † The calculated FWHM values for various atomic distances of h-BCN at a temperature of 300 K is shown in Table 1. Compared to the C-N bond (1.42Å), the C-B bond (1.52Å) is weaker and is highly sensitive to temperature. From Fig. S1, † it is clear that B-B and B-C show a large increase in FWHM around 1000 K. Moreover, the absence of a chemical bond in B-B is responsible for the bond weakening. Compared to the C-C atomic distance of graphene, a lower FWHM value of h-BN and h-BCN implies that the structure of h-BCN is robust and stable possessing better structural integrity with respect to the empirically observed interatomic potential. It is also worth noting that the FWHM of h-BCN shows a similar trend to that of h-BN. Even though h-BCN, h-BN, and graphene have similar hexagonal lattice structures, the different atom constituents lead to visible changes in the structural and energetic stability. In addition to this, the C-C bond in the h-BCN hetero-structure is responsible for an improved FWHM compared to graphene. A similar pattern of structural robustness was reported for hydrogenated h-BN nanotubes 20 and h-BN nanosheet. 21,22 B Linear thermal expansion coefficient The structural stability of 2D atomic crystals and surfaces have been a point of disagreement in the theory of condensed matter. Based on the Mermin-Wagner theorem, 23 which states that long wavelength uctuations destroy the long-range order for 2D crystals. Similarly, 2D membranes embedded in a 3D space tend to be rippled. However, these dangerous uctuations can result in a two-dimensional membrane, and exhibit strong height uctuations. But in the case of graphene and 2D h-BN, the repression of long wavelength ripples takes place via a strong anharmonic coupling between the in-plane stretching and the out-of-plane bending modes. The height uctuations on the surface are known as ripples, 24 and they essentially stabilize these 2D membranes. 25, 26 We expected similar out-of-plane excursions in h-BCN sheet. One advantage of classical MD simulations is that they can include the full anharmonicity of the inter-atomic potential. Another advantage is that the simulations can include millions of atoms. We studied the  Fig. 3a and b display the temperature dependence of the in-plane lattice parameter and the LTEC. We found that the lattice parameter decreases as temperature increases. The LTEC is computed by the direct numerical differentiation of the lattice constant a using the equation: The LTEC is negative in the entire computed temperature interval up to 1000 K. This is consistent with the case of graphene and h-BN, as tabulated in Table 2. We also found that the temperature evolution of the a-lattice is system size-dependent (Fig. 2). This system size dependence of the a-lattice is due to the existence of large-scale ripples 24,32 in h-BCN, which is similar to graphene and h-BN and in the actual physical context, which cannot be tted adequately inside the small simulation cells. 33 The high thermal stability of h-BCN can be effectively utilized in the manufacturing of h-BCN composite materials. Knowledge of how h-BCN composite materials react to high temperatures is vital to optimizing the manufacturing process. The negative thermal expansion coefficient of h-BCN composites may reduce the thermal stresses and weaken the bonds as temperature rises.

C Evaluation of the anisotropy in the elastic constants
This section details the results of the elastic constants of the h-BCN sheet at zero kelvin temperature, which we computed utilizing molecular statics simulations. We have carried out   calculations pertaining to an h-BCN sheet of innite spatial extent along the x-and y-axes, which correspond to the zigzag and armchair directions, as shown in Fig. 2a. A detailed description of the calculation procedure can be found in our previous work. 34 We evaluated the elastic behavior of h-BCN sheet by employing the energy-strain approach. The mathematical explanation of the in-plane interactions of the matrix of a second-order elastic constant is denoted as: where, E, A 0 , d 0 , and 3 are the equilibrium energy, effective area and the effective thickness (the van der Waals distance), and the strain tensor associated with the calculation. The elastic energy per unit area E(3) can be expressed using the polynomial form: We calculated the four non-zero 2D elastic constants: C 11 , C 22 , C 12 and C 66 (1-xx, 2-yy, 6-xy in the standard Voigt notation) of the h-BCN sheet by tting the strain-relative energy curve, as shown in Fig. 4 GPa), suggesting that h-BCN has good mechanical properties with high anisotropy. The highly anisotropic characteristic of h-BCN is attributed to its structural morphology. It has also been observed that h-BCN possess higher stiffness in comparison to other 2D materials such as monolayer SiC, 36 C 3 N, 13 and pentagraphene. 37 This suggests that h-BCN has potential use in the elds of transportation, aerospace, power generation and energy storage. In our study, the anisotropic behavior is clearly observed from the direction dependent elastic constants, and the h-BCN sheet also satises Born's mechanical stability criterion.

D Mechanical fracture behavior
In order to study the fracture behavior of h-BCN nanosheet, simulations were carried out using an optimized cut-off function of the interatomic potential. During the initial stages of stressstrain analysis, we carried out the relaxation for 200 ps in the NPT ensemble. The change in temperature and potential energy of the h-BCN system during the relaxation is shown in Fig. S2. † Even though we performed relaxation for 200 ps, we observed that the system reaches a properly relaxed state aer 50 ps, as shown in Fig. S2a and b. † Aer the relaxation, we applied deformation by pulling the h-BCN sheet along the x-direction (zigzag) or the ydirection (armchair) with temperatures ranging from 100 to 1000 K, in conjunction with three different strain rates: 0.0001, 0.001 and 0.005 ps À1 . We observed that this material possesses strong anisotropic characteristics that depend on numerous factors, including temperature, strain rate, and loading direction. MD simulations were performed to investigate the effect of strain rate and temperature on the mechanical responses of h-BCN sheet. We  . 4 The relative energy versus strain response graph of h-BCN sheet using molecular statics simulations for the analysis of elastic constants. The relative energy is plotted as a function of various specific strains necessary for calculating the elastic tensors.
monitored the engineering stress-strain response of h-BCN within a temperature range of 100-1000 K, under tension load along both armchair and zigzag directions. The observed anisotropic stressstrain response along the armchair and zigzag directions are shown in Fig. 5a and b, respectively. The fracture strength (engineering stress) and the strain decrease over the studied temperature range along both the armchair and zigzag directions. At 300 K, the observed fracture strength is 111 GPa along the armchair direction and 100 GPa along the zigzag direction with a strain rate of 0.001 ps À1 . The corresponding fracture strains are 0.17 and 0.15 respectively along the armchair and zigzag direction. The lower ultimate tensile strain and stress for the zigzag orientation is likely due to the deformation of bonds, which are directly subjected to the load at the beginning of the elongation process, whereas in the armchair situation the external loads are weakened through both bond elongation and bond angle variation. We also observed that the armchair deformations exhibit a mechanically strengthening behavior in the stress-strain evolution process. This can be explained on the basis of the spatial bond distributions and bond breakage mechanism, which suggest that more ductile materials may exhibit such features during the initial stages of the stressstrain evolution procedure. 38 Fig. 6 illustrates the stress-strain responses of the tensile properties of the h-BCN sheet at different strain rates and temperatures. We performed the simulations with temperatures ranging from 100 to 1000 K, in conjunction with three different strain rates (0.0001, 0.001, and 0.005 ps À1 ) along the armchair direction ( Fig. 6a and b) and the zigzag direction ( Fig. 6c and d). The maximum stress in both armchair and zigzag directions occurs at a strain rate of 0.001 ps À1 . The lack of proper experimental evidence limits any direct comparisons of the obtained stress-strain response of h-BCN with that of other materials. However, our approach is qualitatively in agreement with both the experimental and theoretical investigations of graphene and h-BN. [39][40][41][42][43] This qualitative agreement substantiates the validity of our numerical model. From our stress-strain investigations, we concluded that the failure  morphology of h-BCN sheet is dependent on the strain rate. 44 We also studied the failure morphology of h-BCN upon tensile loading. In most materials, the fracture mechanism is affected by temperature and strain rate. Due to h-BCN's atomically thin nature, it easily undergoes various types of fracture mechanisms and deformations, despite its high stiffness. Bond length and bond angle of both the pristine h-BCN sheet and the strained h-BCN sheet were studied, with varying orientations. Fig. 7a shows the calculated bond length and bond angles of an unstrained pristine h-BCN sheet. We observed considerable changes in bond length and bond angle aer applying a deformation in the armchair and zigzag directions, as shown in Fig. 7b and c respectively.
The structural details of an individual atomic network determine the material's strength in any particular direction. As we subjected the h-BCN network to uniaxial strain, its response was manifested by changes in the bond angles as well as the bond lengths. At a completely relaxed zero strain state, we observed a single bond angle around 120 in the angle distributions of one hexagonal ring in the h-BCN (Fig. 7a). This conguration corresponds to the most unrestrained geometry of the sp 2 hybridization of the atoms in the BCN network. As soon as we applied strain to the system, it deformed by changing the bond lengths and angles. The strained congurations along the armchair and zigzag directions are shown in Fig. 7b and c, respectively. For each case, we observed a notable difference between the armchair and zigzag directions in terms of their bond and angle distributions. The widened angle is higher in the armchair case than in the zigzag case. A similar pattern of bond angles and lengths has been reported in the case of other 2D materials such as biphenylene, phagraphene and nitrogenated holey graphene. 45 We attribute the high stretchability along the armchair direction to the signicant widening of the bond lengths and angles in the h-BCN system. We also found that certain stages of deformation in h-BCN involve the formation of defects, cracks, and the breaking of bonds and that the edges of the h-BCN material are easily affected by the deformation, as shown in Fig. S3 and S4. † In Fig. 8, the failure process of h-BCN sheet along armchair and zigzag directions at 300 K are depicted. During the initial stages of tensile loading, bond elongation followed by bond breaking happening in both armchair and zigzag orientations, as shown in Fig. S3 and S4. † Further increasing the strain level along the armchair and zigzag directions lead to more bond breakage. The formation of chain type rupture in the sheet indicates the ductile nature of mechanical failure for h-BCN sheet. Unlike pristine graphene, the bond breaking and failure process happen at distinctly different strain levels in h- BCN sheet, which also infers a ductile failure process in the material. The ductile failure behavior predicted by the Tersoff potential is in good agreement with the earlier atomistic studies of pristine and defective graphene, 39,40 graphene-like C 3 N, 13 boron nitride nanosheets, 46,47 and phagraphene 48 using the Tersoff-like potential. Based on our calculations we found that h-BCN monolayer shows a higher ductile failure along the armchair direction when compared to the zigzag direction. This can be conrmed by the obtained failure strain values and is found to be 0.17 and 0.15 respectively along the armchair and zigzag direction at 300 K. In the present work, we have focused on the mechanical stability analysis of pristine monolayer BCN. We are of the opinion that incorporating defects in the computational studies could provide more realistic picture prior to experimental analysis. In view of this, a systematic study of the thermomechanical properties of defective monolayer BCN would also be worthwhile.

Conclusion
In this study, we performed classical MD simulations based on the optimized Tersoff potential to investigate the thermomechanical properties of an h-BCN monolayer. We examined the effects of temperature, system size, and tensile properties on h-BCN. Our major ndings are as follows: (1) the narrow FWHM of h-BCN as compared to graphene is attributed to the better structural integrity of h-BCN with respect to the used empirical interatomic potential which leads to a lesser spread in the interatomic distances. The C-C bond in the h-BCN heterostructure is also responsible for an improved FWHM compared to graphene. (2) At nite temperatures, a thermal contraction is observed and is due to the presence of out-ofplane bending modes excited, along with thermal energy which reduce the thermal stresses and weakening the bonds as temperature rises which guarantee the durability of h-BCN based composite materials for better practical applications.
(3) Akin to graphene and h-BN in the hexagonal lattice, the calculated Young's modulus of h-BCN is smaller than that of graphene, but is higher than that of h-BN monolayer, suggesting that h-BCN has good mechanical properties with high anisotropy. (4) The elastic moduli of h-BCN decreases with temperature and is affected by the loading direction. (5) Finally, the lower ultimate tensile strain and stress for the zigzag orientation is likely due to the bonds deformations, which are directly subjected to the load at the beginning of the elongation process, whereas in the armchair situation the external loads can are more spread out via both bond elongation and bond angle variation. We believe that the results we obtained may help in the design of h-BCN-based composite materials for large-scale applications.

Conflicts of interest
There are no conicts to declare.