Vacancy-Induced Thermal Transport and Tensile Mechanical Behavior of Monolayer Honeycomb BeO.

Because of the rapid shrinking trend of integrated circuits, the performances of nanodevices and nanomechanical systems are greatly affected by the joule heating and mechanical failure dilemma. In addition, structural defects are inevitable during experimental synthesis of nanomaterials, which may alter their physical properties significantly. Investigation of the thermal transport and mechanical behavior of nanostructured materials with structural defects is thus a crucial requirement. In this study, the thermal conductivity (TC) and tensile mechanical behavior of monolayer honeycomb BeO are systematically explored using molecular dynamics simulations. An infinite length bulk TC of ∼277.77 ± 8.93 W/mK was found for the pristine monolayer BeO. However, the insertion of 1% single vacancy (SV) and double vacancy (DV) defects reduces the TC by ∼36.98 and ∼33.52%, respectively. On the other hand, the uniaxial tensile loading produces asymmetrical fracture stress, elastic modulus, and fracture strain behaviors in the armchair and zigzag directions. The elastic modulus was reduced by ∼4.7 and ∼6.6% for 1% SV defects along the armchair and zigzag directions, respectively, whereas the reduction was ∼2.7 and ∼ 5.1% for 1% DV defects. Moreover, because of the strong symmetry-breaking effect, both the TC and mechanical strength were significantly lower for the SV defects than those for the DV defects. The highly softening and decreasing trends of the phonon modes with increasing vacancy concentration and temperature, respectively, were noticed for both types of defects, resulting in a reduction of the TC of the defected structures. These findings will be helpful for the understanding of the heat transport and mechanical characteristics of monolayer BeO as well as provide guidance for the design and control of BeO-based nanoelectronic and nanoelectromechanical devices.

Introduction

With the advent of nanotechnology, the scale of transistor integration and the demand for miniaturized chips and circuits have increased considerably to achieve superior device performances.[1,2] However, significant Joule heating, which can lead to the early damage of the miniaturized chips and circuits, is a great disquiet in nanodevices.[3] In conjunction with the excessive Joule heating, the mechanical failure of nanostructures also significantly hampers the reliability of nanoelectronics and nanoelectromechanical systems.[4,5] Reliable nanostructured materials with superior thermal conductivity (TC) and extremely high mechanical strength are thus mandatory to address the challenges of future nanoscale systems.[6−9] Although the first exfoliated monolayer system, graphene,[10] finds substantial attention in solving the heat dissipation and structural failure problem, the zero bandgap aspect makes it implausible in realistic uses.[11−13] Beyond graphene, monolayer beryllium oxide (BeO), a compound of Be and O atoms with a honeycomb crystal structure, has attracted great attention in nanoscale device fabrication because of its prospective features of large electronic bandgap, extraordinary mechanical strength, superior piezoelectric behavior, high excitonic binding energy, and surprising optical properties.[14−17] Recently, combining auger electron spectroscopy, low-energy electron diffraction, and first-principles density functional theory (DFT) calculations, Afanasieva et al.[18,19] experimentally synthesized BeO nanosheets on a Mo (112) substrate. Moreover, using molecular beam epitaxy, Zhang et al.[20] has shown the first experimental growth of monolayer BeO on the Ag(111) surface. In addition, Luo et al.[21] uncovered that monolayer BeO possesses extraordinary dynamic, thermal, and mechanical stabilities. The advancement of BeO synthesis has sparked a lot of interest in the contemporary opto- and nanoelectronics applications.

During the large-scale growth of nanoscale materials, complex fabrication technology often produces different structural defects. The law of thermodynamics demonstrates that even minor structural flaws in nanoscale materials have a significant impact on their thermal, optical, electrical, and mechanical properties.[22] Among the different structural defects, single-vacancy (SV), and divacancy (DV) are crucially important,[23−26] and they can affect the thermal transport and tensile mechanical behavior of the nanomaterial greatly owing to their robust symmetry failure effect.[27,28] These defects are responsible for the localization of phonons, which are the foremost heat carriers in low-dimensional materials and cause strong phonon–phonon and phonon-vacancy scattering.[29] At elevated temperatures, when the defect concentration is increased, the low-frequency phonons produce Umklapp type phonon scattering, which can substantially decrease the TC of the material.[23] On the other hand, vacancy defects can significantly affect the mechanical performance of the nanostructures by altering their fracture behavior.[5,27] Several investigations suggest that a defect-induced structure breaches at a breaking stress that is substantially inferior to the ideal strength of the perfect sheet, and increasing the density of defects can significantly reduce the fracture stress, elastic modulus, and fracture strain.[30−32] To apply monolayer BeO in realistic applications, particularly in environments under which the system operates at high temperatures or at high loadings, it is important to discover the dependence of the mechanical behavior of monolayer BeO on vacancy defects.

Recently, the thermal and mechanical properties of two-dimensional (2D) BeO have been investigated by utilizing various theoretical approaches.[17,33] Xia et al.[17] calculated the in-plane TC of 2D-BeO as 266 W/mK with the phonon Boltzmann transport equation. They found that the TC of 2D-BeO reduces significantly with the increasing temperature, and the TC of 2D-BeO is lower than that of graphene and h-BN due to the smaller phonon group velocity and phonon relaxation time. They also concluded that 75% of the TC was caused by phonons with frequencies in the range from 0 to 5.4 THz and that the phonons exhibit a saturated nature after 20 THz. Using the equilibrium molecular dynamics (MD) method with a modified embedded atom (MEAM) potential, Wei et al.[33] also calculated the TC of 2D-BeO at different temperatures. They found the TC as ∼300 W/mK at 300 K and due to the dominant effects of Umklapp phonon scattering, the TC shows a decreasing nature with the increasing temperature. Though these works have encompassed temperature-dependent thermal properties of monolayer BeO, it is yet unclear how the vacancy defects will influence the TC of monolayer BeO. The combined effect of vacancy defects and temperature on the TC of 2D-BeO is also unknown. On the other hand, using MD simulations, Zarghami Dehaghani et al.[14,34] explored the effects of temperature, cracks (circular and square shaped notches), and grain size on the mechanical properties of monolayer BeO. Although few works have been performed on the mechanical properties, the different structural vacancy effects such as the DV and SV with different concentrations in tailoring the mechanical properties of BeO for different chiral orientations are still missing. It is therefore crucial to explore the vacancy-defected thermal transport and mechanical behavior of monolayer BeO using some reliable techniques in order to include its applications in next-generation nanoelectronics and nanoelectromechanical systems.

In this paper, we have performed a thorough and systematic investigation of the thermal transport and tensile mechanical performance of monolayer BeO through MD simulation. Classical Tersoff potentials developed by Byggmästar et al.[35] are utilized for the exploration of the TC and fracture behavior. The different vacancy defects such as SV and DV with a concentration ranging from 1 to 5% are used for the investigation of TC and tensile mechanical strength. For efficient heat conduction, the combined effects of defects and temperature as well as for the validation of our employed simulation, the infinite length TC has also been investigated. Among the two considered defects, the SV defects show a more substantial decrease in TC and fracture strength than the DV defect due to its greater regularity failure effect. We have calculated the phonon density of states to explain the difference of the calculated TC at different lengths, vacancies, and temperatures. In addition, we explore that mechanical behaviors such as elastic modulus, fracture strength, and fracture strain show a considerable anisotropic nature for this monolayer BeO.

Computational Details

In order to calculate the TC and tensile mechanical properties of monolayer BeO, we have employed classical MD simulations. To describe the atomistic interactions concerning Be–Be, O–O, and Be–O atoms in the MD simulation, classical Tersoff empirical bond order potentials developed by Byggmästar et al.[35] have been utilized. This potential has been developed solely for the nanosheet and nanotube structure of BeO. Furthermore, the Tersoff potential comprises two body as well as three body terms, which are very much necessary for exposing the bond breaking and bond forming of monolayer BeO. According to Byggmästar et al., the Tersoff potential correctly produces the experimental cohesive energies, lattice parameters, elastic constants, phonon scattering rates, acoustic phonon branches, and defect formation energies of monolayer BeO systems, as verified by first principles calculations and experiment.[35]

First, for the estimation of the TC, we have used the reverse non-equilibrium molecular dynamics (RNEMD) technique[36] via the LAMMPS package.[37] The schematic arrangement of the monolayer BeO structure utilized for the TC calculation is depicted in Figure a. The zoomed-in view shows the typical DV and SV defects that can be formed experimentally when monolayer BeO is synthesized from the bulk structure.[38−40] The whole monolayer structure was divided into N slabs, each one with a thickness δ. We also applied periodic boundary conditions (PBC) along the X and Y directions to remove the edge as well as the size effect in RNEMD techniques. In addition, a large vacuum space of nearly 30 Å was taken along the out-of-plane direction to remove interlayer interactions. The conjugate gradient algorithm (CGA) was utilized for energy optimization of the structures, and all the calculations are performed with a time step of 0.5 fs. Finally, the TC along the sheet length (X) is assessed by using the “heat flux” as the “simulation input” and determining the “temperature gradient” as the “simulation output” with Fourier’s law:where J denotes the heat flux, dT/dx denotes the temperature gradient, and K designates the TC along the length (X) direction.

Figure 1

(a) Structural view of monolayer BeO for the computation of TC using the RNEMD process. The zoomed-in view represents the realization of SV and DV defects into the BeO sheet. The whole nanosheet (length) was allocated into N slabs, each with thickness δ. At a distance of X = “L/4”, a “hot slab” and a distance of X = “3L/4”, a “cold slab” is formed by injecting and removing a small amount of heat ΔQ, respectively, using a time step of 0.5 fs. (b) Atomistic view of vacancy-defected monolayer BeO for the calculation of uniaxial tensile mechanical behavior along the two chiral directions (armchair and zigzag).

Before calculating the steady-state temperature gradient, all the BeO sheets were relaxed with an isothermal-isobaric ensemble (NPT) simulation for 2 ns and then switched to a microcanonical ensemble (NVE) simulation for 2 ns. After the structural relaxation, a kinetic heat energy of amount ΔQ was constantly added into the heat source region (at a distance of L/4 of the sheet), and correspondingly a similar energy of −ΔQ was removed from the heat sink (at a distance of 3L/4 of the sheet) at each time step of 0.5 fs to produce a heat flux along the length direction (X) of the simulation system. The microcanonical ensemble (NVE) simulation for 2.5 ns is used to carry heat from the heat source to the heat sink. Throughout this evolution, the total energy of the structure stays constant and the accumulation (withdrawal) of energy can be accomplished by dividing (multiplying) the velocity of atoms. The created heat flux J along the X direction is then expressed aswhere t, M, and A represent the time step, frequency of heat energy alteration, and sheet area (equal to wh), respectively. The Δε represents the transferring energy difference between the coldest and hottest atoms in each time step of M. In steady-state circumstances, Δε can be stated aswhere vcold, vhot, and m denote the highest velocity atoms in the cold slab, lowest velocity atoms in the hot slab, and atomic mass, respectively. The 1/2 factor is incorporated owing to the structure’s periodicity, i.e., heat energy can transfer in both routes.

After calculating the heat flux, we again run the system for another 15 ns to determine the time-averaged temperature profiles; those are lastly required for the calculation of the dT/dx. Finally, we treat the linear temperature zone with the least-squares method to achieve the temperature gradient dT/dx so that the TC can be calculated using , by the Fourier’s law of heat transfer. When computing the in-plane TC, usually, the distance among the neighboring sheets in the bulk structure is chosen as the sheet thickness. Here, the thickness of monolayer BeO is used as 3.05 Å, as stated in an earlier work,[41] and the heat alteration frequency(M) is selected as 1000.[23,42]

Second, the mechanical stress of monolayer BeO under uniaxial tensile deformation was investigated by virial stress, which can be stated as[43]where the total sum is taken for all the atoms in the volume, m signifies the mass of atom i, u signifies the time derivative of the interchange, r signifies the location vector, and f signifies the interatomic force applied on atom i by atom j.

The atomic arrangements of monolayer BeO from its top view is shown in Figure b. Here, the two fundamental chiral orientations, namely, the armchair and the zigzag directions, are shown along the X and Y directions, respectively. For both chirality, the total atoms considered in the monolayer BeO were 13,376 and the dimension of the whole sheet was 419.5 nm2 (20.13 nm in the X direction and 20.08 nm in the Y direction). In order to remove the size effect, PBC were employed along the two (X and Y) directions. A large vacuum space of 15 Å was considered on both sides (along the Z direction) of the monolayer so that the particles do not interact along the Z direction, and only monolayer BeO was involved in the calculation. The CG method was employed for the energy minimization of the system. The ending tolerance for force and energy were taken as 10–8 eV/Å and 10–8, respectively. In the minimization process, the highest force/energy evaluations and the maximum iterations of the minimizer were considered as 10,000 and 5000, respectively. After completing the energy minimization, an equilibration simulation for 50 ps is applied by means of NVE microcanonical ensemble. An NPT ensemble for 40 ps was further employed to equilibrate the system pressure to 1 bar. During the NPT ensemble simulation, the system gets reconfigured and all the inner stress developed was eliminated. Lastly, we employed the NVT ensemble for 20 ps for thermal relaxation. A temperature of 300 K was used to achieve thermal equilibration. In each of the equilibration processes, we confirmed appropriate convergence of the energy (potential) and preferred thermodynamic extents. The equation of motion was solved through the velocity Verlet integration technique using a time step of 1 fs. We have used a stable strain rate of 109 s–1 along the tensile deformation orientations (armchair and zigzag). Though the applied strain rate is relatively large compared to realistic situations, this high strain rate is very much appropriate in atomistic calculations to explore the material fracture behavior with an affordable computational cost.[27,44] For performing the post-processing task and documenting the calculated strain, stress, and the trajectories of the considered structure, the open-source visualization software OVITO was used.[45] Furthermore, it is well recognized that the trajectory of the classical MD calculations is extremely stochastic. Therefore, all the MD calculations were performed for four different initial settings (i.e., molecular velocity) for both perfect and vacancy-defected structures to incorporate the influences of stochasticity and ambiguity.[46] We took an average from the four outcomes of these four distinct computations to accomplish the error assessment.

Results and Discussion

The temperature distribution profile has a crucial effect on the accuracy of the TC calculation in the RNEMD approach. We have thus considered a long simulation time to estimate the temperature gradient (dT/dx) of the system. Figure a–c shows the steady-state temperature profile along the length (X) of a 100 nm × 10 nm sheet for three different simulation times: 4 ns, 8 ns, and 12 ns. The temperature in the hot and cold regions is maintained as T + ΔT and T – ΔT, respectively, where ΔT is ∼25 K. Due to the high phonon scattering, a strong non-linearity was found near the hot and cold regions. A linear temperature profile is generated between the hot and cold zones as a result of the quick energy transfer. This linear region is taken into account while calculating the dT/dx. Moreover, it has been found that with the increase in the simulation time, the linear region of the temperature profiles shows greater smoothness and linearity with length, making it appropriate for the precise calculation of the TC. As the slope of the temperature profiles for different simulation times is different, thus, the average of the TC is taken. The heat flux needed to determine the TC using Fourier’s law is estimated from the energy profile of the hot and cold regions shown in Figure S1 (Supporting Information). The heat flux was estimated to be 1.94589 eV/ps.

Figure 2

Calculated temperature gradient profiles of monolayer BeO at 300 K for simulation times of (a) 4 ns, (b) 8 ns, and (c) 12 ns. The BeO nanosheet under consideration has a length of 100 nm and a width of 10 nm. The linear temperature division between hot and cold slabs is represented by the blue, lime, and aqua colors, respectively, which was further fixed with the least-squares procedure to obtain dT/dx for the TC calculations.

First, we have calculated the length-dependent TC of pristine monolayer BeO by varying the length of the sheet from 10 to 400 nm with a fixed width of 10 nm. In the RNEMD method, the considered monolayer structure width of 10 nm was sufficient to eradicate edge effects, and we identified no undesirable temperature profiles in our simulation. The length-dependent TC at room temperature is shown in Figure a. With the increase in sheet length, the TC increases. For a sheet length of 100 nm, we have found a TC of ∼116.46 ± 6.38 W/mK and when the sheet length is increased to 400 nm, the TC is estimated as ∼152.99 ± 7.13 W/mK, which is nearly 31.36% greater compared to the 100 nm length. The infinite length bulk TC of monolayer BeO can be found by extrapolating the TC of the considered structures (i.e., 10–400 nm), using the following relation[23,29]where K, L, and Λ symbolize the infinite length bulk TC, sheet length, and the phonon mean free path, respectively. By plotting the inverse sheet length and corresponding inverse TC of the monolayer BeO (shown in Figure b), the infinite length bulk TC was found to be ∼277.77 ± 8.93 W/mK. Our computed infinite length bulk TC of monolayer BeO is nearly equal to the calculation of Xia et al.[17] (266 W/mK) using the first-principles phonon Boltzmann transport equation and that of Wei et al.[33] (300 W/mK) using equilibrium MD simulations with the MEAM potential. In addition, our estimated TC shows a very good agreement with the first-principles results calculated by Mortazavi et al.[15] By means of full iterative solution of the Boltzmann transport equation, they estimated a superior TC value of 385 W/mK for this monolayer BeO, which is significantly higher than our obtained TC value of ∼277.77 ± 8.93 W/mK. The underestimation of the TC calculation may be endorsed due to the consideration of higher phonon scattering rates that are not counted in DFT-based results. Hence, our employed RNEMD approach and potential show very good agreement with the earlier thermal behavior prediction of the BeO nanosheet.

Figure 3

(a) Length-dependent TC at 300 K. (b) Inverse correlation between TC and length utilized for the determination of the infinite length bulk TC.

Structural defects such as vacancy can affect the thermal transport behavior of the nanoscale material in an unexpected way, and linking the effect of vacancy defects may be an appropriate process to regulate the phononic TC of monolayer BeO. The missing atoms in the structure of a monolayer system are normally termed as vacancy defects. The vacancy defects can be again classified into single, double, and multi-vacancy defects based on the number of atoms missing from the pristine configuration.[23] The loss of a single atom produces an SV and creates three dangling-like bonds due to the three under coordinated edge atoms. In contrast, the loss of two atoms from the pure monolayer produces a DV defect. Although it contains under-coordinated atoms in the SV configuration, the DV-induced configuration does not produce any under-coordinated atoms. Therefore, a DV causes a lesser effect on the phonon properties compared to the SV defect.[27] Here, we have investigated the SV- and DV-induced TC of monolayer BeO using vacancy concentrations ranging from 1.0 to 5.0%. To explore the impact of vacancies on the TC of monolayer BeO, the dimension of the 2D-BeO was taken as 100 nm × 10 nm. The TC of monolayer BeO as a function of vacancy concentration for both SV and DV is depicted in Figure a. As the figure suggests, the introduction of vacancies causes the TC of monolayer BeO to exhibit a decreasing behavior for both types of defects. For a vacancy concentration of 1.0%, the estimated TC of a 100 nm × 10 nm monolayer BeO for the DV and SV is 77.42 ± 5.23 W/mK and 73.38 ± 6.09 W/mK, respectively. On the other hand, when the concentration is increased to 5.0%, the DV- and SV-induced TC is found to be 26.35 ± 6.56 W/mK and 21.96 ± 5.21 W/mK, respectively. Because of the DV and SV defects, a comparable decreasing trend of the TC was also obtained for graphene, ZnS, MoS2, and SiC.[47,48] The percentage (%) decrease in the TC of the monolayer BeO owing to the different vacancy concentrations can be estimated by the following formula:

Figure 4

DV- and SV-induced (a) TC and (b) percentage decrease in the TC of monolayer BeO for different defect concentrations at 300 K. The nanosheet used for the defect-induced TC calculation was 100 nm × 10 nm.

The calculated percentage decrease in TC with different concentrations of SV and DV is shown in Figure b. We have found that for a vacancy concentration of 1.0%, the TC of monolayer BeO experienced a percentage decrease of 36.98% for the SV and 33.52% for the DV. When the concentration is increased to 5.0%, a decrease of ∼81.14% was obtained for SV and a decrease of 77.37% was found for DV. As the SV was formed by removing one Be or O atom, giving rise to three two-coordinated Be or O atoms, it can effectively damage the symmetry of the sp2 configuration and produce a less stable structure. On the other hand, a DV defect was created by removing two bonded atoms (Be and O); hence, the native configuration can reorganize to reconstruct to the three-coordinated sp2 connection by producing two octagon and pentagon arrangements.[49] In the SV-induced structure, a larger degree of phonon-defect scattering thus occurs compared to the highly stable three-coordinated sp2-bonded DV defect and makes a greater degree of reduction in TC compared to the DV. Moreover, for a similar concentration of vacancy, the SV-induced structure is highly anomalous compared to the DV-induced structure and causes a greater decrease in the phonon mean path owing to the scattering of phonons from the different distributed vacancy centers.[50] As Figure a,b shows, a lower TC reduction is thus found for the DV defect than for the SV defect. These types of findings are analogous and in good agreement with another study.[51] Moreover, to understand the dependency of the TC behavior on the types of vacancies, we have further considered two types of SV: SVBe (by removing only Be atoms) and SVO (by removing only O atoms). The obtained results are compared to the earlier SV made by removing Be and O atoms from different locations (SVBe&O). It is found that the SVBe shows higher formation energy and stability than the SVO and SVBe&O. The SVBe&O defect shows the greatest reduction in the TC for the same vacancy concentration. Conversely, the SVBe defect shows the smallest decreasing trend in the TC at different vacancy concentrations. The formation energy and calculated TC value for SVBe, SVO, and SVBe&O defects at different concentrations are presented in Tables S1 and S2 (Supporting Information).

Effective thermal management at different temperatures is a critical issue for the stability and reliability of nanoelectronics. Because of the increased workload, there occurs significant temperature rises in the electronic devices.

Considering the temperature as well as the combined effect of temperature and vacancy defects on the TC of monolayer BeO is thus an important issue and must be investigated before applying BeO in realistic applications. Figure a shows the variation of the TC of pristine monolayer BeO as a function of temperature, ranging from 100 to 700 K. As the figure suggests, with the increase in temperature, the TC of monolayer BeO will be significantly reduced. At 100 K, a TC of ∼253.22 ± 7.55 W/mK is found for the considered pristine monolayer BeO of the dimension 100 nm × 10 nm. However, when the system temperature is increased to 700 K, the reduced value of TC is obtained as ∼53.79 ± 6.27 W/mK, which is about ∼78.75% lower compared to the TC at 100 K. Moreover, by considering a 1% vacancy concentration, as well as temperature ranging from 100 to700 K, the combined effect of temperature and vacancy on the TC of monolayer BeO is depicted in Figure b. As illustrated in Figure b, with the increase in temperature from 100 to 700 K, the TC shows a monotonically decreasing trend for both the SV- and DV-induced BeO. At all the temperatures, the SV induced BeO shows a greater reduction of the TC than the DV induced system. This trend of defect-induced TC at different temperatures is comparable to those of earlier investigations on 2D materials.[52−54] From Figure b, it can be seen that at room temperature for 1% SV and DV defects, the TC of monolayer BeO is estimated as ∼73.38 ± 6.09 W/mK and ∼77.42 ± 5.23 W/mK, respectively. When the temperature is increased to 700 K, the TC for SV- and DV-induced BeO is decreased to ∼36.60 ± 5.21 W/mK and ∼40.83 ± 5.05 W/mK, respectively. Hence, the percentage reduction of TC for SV defect is much greater than that of the DV defects. The temperature-dependent reduction of the TC of vacancy-defected BeO can be explained in terms of phonon scattering due to vacancies as well as Umklapp phonon scattering due to temperature.[55,56] With the increase in temperature, the density of high-frequency phonons increases. Hence, the Umklapp phonon scattering also increases. Consequently, the integrated effect of the phonon scattering produced by vacancies as well as the Umklapp phonon scattering produced by high temperatures lessens the TC a significant amount. Moreover, from Figure b, it can be noticed that the temperature-dependent decreasing trend of the TC of SV- and DV-induced BeO is not the same. The TC is influenced by both types of scatterings (Umklapp and phonon-vacancy scattering) at temperatures below 350 K, resulting in a very high TC reduction rate. On the other hand, the TC reduction due to temperature does not change rapidly above 350 K, and the predicted TC for both types of defects is nearly the same. The TC is only governed by the Umklapp phonon scattering at higher temperatures, and the phonon scattering by vacancy defects is less important than the Umklapp phonon scattering.[55,56] Therefore, the TC decreases slowly with temperature, and the difference between SV- and DV-induced TC is negligible.

Figure 5

TC variation of (a) pristine and (b) vacancy-defected BeO nanosheet as a function of temperature.

Phonons are primarily responsible for heat conduction in semiconducting materials. With the variation of different parameters such as constituent atoms, system dimension, vacancy concentration, and system temperature, the phonon modes show a varied nature and explain the variation of the estimated TC.[11] The phonon density of states (DOS) calculation is a valuable property that can be used for the quantitative interpretation of the heat transportation behavior of nanomaterials under different circumstances. The phonon DOS of monolayer BeO can be computed using the Fourier transformation of the VACF[23] with the following relation:where v(0) and v(t) represent the velocity of the ith component of species α (Be or O) at time 0 and time t, respectively, and the bracket denotes an average over time and atoms species. Lastly, the phonon DOS is estimated with the following formula:where F(ω) represents the complete phonon DOS, Fα(ω) indicates the partial phonon DOS of atoms Be or O, and ω signifies the wavenumber. The partial phonon DOS can be quantified as . The attained VACF for both the Be and O atoms are shown in Figure S2a,b (Supporting Information), respectively. The VACF shows an enormous extent of oscillation at small correlation period. However, the oscillating character diminishes with the increasing correlation period because of the interactions between the atoms and the classical forces from nearby atoms. The obtained VACF is then utilized to calculate the partial phonon DOS and finally the total phonon DOS of monolayer BeO.

The effect of length on the TC of monolayer BeO can be easily understood from the length-dependent phonon DOS behavior shown in Figure . As the figure indicates, because of the enhanced availability of the phonons, a monotonic, increasing trend of phonon DOS peak intensity with an increasing sheet length is found in this study. At room temperature, it is expected that all the acoustic mode phonons (ZA, TA, and LA phonons) can participate in the heat transport and their involvement would be intensified when the sheet length is increased.[29,57] Moreover, at room temperature, the effect of optical phonons on the heat transport of binary compounds is assumed to be negligible and can be understood from the sharp increasing trend of the phonon DOS peaks (they do not show any redshift or blueshift behavior). Thus, this character of the phonon DOS can control the increasing nature of the TC of monolayer BeO at an increased sheet length.

Figure 6

Length-dependent phonon DOS curves of monolayer BeO at a 300 K temperature.

We have also calculated the phonon DOS of monolayer BeO at room temperature by considering the SV and DV with different concentrations. To quantitatively verify the difference between the calculated TC of both pristine and vacancy-induced structure, we have also shown a comparison of the phonon DOS considering the pristine and SV- and DV-induced structure. Figure a,b demonstrates the variation of phonon DOS with the vacancy concentrations 1.0, 2.0, 3.0, 4.0, and 5.0% for the DV- and SV-induced BeO, respectively. In addition, the variation of the calculated phonon DOS for 4% DV- and SV-induced BeO compared to pristine BeO is also shown in Figure c. Due to the increase in the defect concentration, there is a reduction of the peak strengths of the phonon DOS for both SV- and DV-induced BeO. This type of reduction in PDOS peaks with defect concentration was also found in other investigations.[23,58] This decreasing trend of strength in phonon peaks with defect concentration occurs mainly due to the strong phonon scattering near the vacancy centers, which decreases the phonon group velocities as well as mean free path, creating the reduction of the TC. Moreover, from the comparison of the SV and DV defect phonon DOS with pristine BeO shown in Figure c, it clearly validates our TC calculation at room temperature.

Figure 7

Effects of (a) DV and (b) SV defect concentrations on the phonon DOS at 300 K. (c) Calculated phonon DOS variation of the pristine and DV (4%)- and SV (4%)-induced BeO nanosheet at 300 K.

In addition, to identifying the decreasing nature of the calculated TC of both pristine and vacancy (1%)-defected BeO with increasing temperature, here, we study the phonon DOS of monolayer BeO at different temperatures. Usually, when the temperature of a nanosheet is increased, the energy of the acoustic phonon modes increases considerably and produces strong phonon scattering phenomena. Conversely, for a binary system, when the system temperature increases, the high-frequency optical phonons participate in the thermal conduction process in parallel to acoustic phonons.[57,59] Therefore, with the increase in temperature, a blueshift occurs in the peaks of the low-frequency range and a redshift nature occurs in the peaks of the high-frequency range of the phonon DOS of the pristine monolayer BeO, as found in our results (Figure a). This type of behavior of the acoustic and optical phonons finally produces phonon modes shrinking (due to blueshift and redshift) as well as peak intensity lowering (due to strong phonon scattering) and is responsible for the reduction of the TC of the pristine sheet with the increasing temperature. Figure b,c also displays the calculated phonon DOS of DV- and SV-induced BeO for three different temperatures. For both types of vacancy, the phonon DOS also exhibits a similar decrease in the individual peak intensity as for the pristine case when the temperature is increased from 100 to 700 K. The first and second peaks of the phonon DOS exhibit a blueshift (changing of the phonon energy from lower to higher), while the third, fourth, and fifth peaks display a redshift (changing of the phonon energy from higher to lower) with the increasing temperature, similar to pristine monolayer BeO. Therefore, with the increasing temperature, a frequency-reducing trend is observed in the phonon DOS of vacancy-induced monolayer BeO. These declines, as well as the frequency contracting behavior of the phonon DOS, are responsible for the reduction of TC with the increasing temperature.

Figure 8

(a) Effect of temperature variation on the phonon DOS of pristine monolayer BeO. (b) Phonon DOS behavior of (b) DV- and (c) SV-defected BeO monolayer at three different temperatures.

The experimental construction of pristine nanostructures is practically unachievable. Numerous vacancy defects are produced throughout the synthesis and can have unexpected consequences on the tensile mechanical behaviors of nanoscale materials. Moreover, vacancy-like defects can be used as an exceptional tool to regulate the physical behavior of monolayer BeO. In this section, attempts are made to determine the effect of vacancies (DV and SV) on the mechanical strength of monolayer BeO under uniaxial tensile loading at room temperature. Figures a,b and 10a,b demonstrate the stress–strain behavior of DV- and SV-induced BeO sheet with vacancy concentrations varying from 1.0 to 5.0% in the armchair and zigzag directions, respectively. For both considered vacancies and directions, the stress–strain relationship represents a linear trend nearly up to the fracture stress and ultimately terminates with a fast deterioration resulting in a brittle fracture behavior of the monolayer BeO. As the figures (Figures a,b and 10a,b) imply, for both directions, the ultimate fracture stress exhibits a decreasing trend when the vacancy concentration is increased. For pristine BeO, the evaluated fracture stress, fracture strain, and elastic modulus along the armchair orientation are found as ∼65.61 ± 1.98 GPa, ∼0.1519 ± 0.023, and ∼388.28 ± 3.98 GPa, respectively, while these values are ∼78.84 ± 2.17 GPa, ∼0.2215 ± 0.025, and ∼385.22 ± 4.45 GPa, respectively, for the zigzag orientation. The obtained tensile mechanical behavior of the pristine monolayer BeO agrees well with previous MD simulation studies by Zarghami Dehaghani et al.[14,34] Mortazavi et al.[15] also calculated the tensile mechanical behavior of monolayer BeO using first-principles DFT calculations. They estimated an elastic modulus of 408 GPa and tensile strength of 53.3 GPa for the pristine BeO nanosheet using a sheet thickness of 3.06 Å. Our estimated tensile mechanical behavior shows a very good agreement with all of these studies.

Figure 9

(a) DV- and (b) SV-defected stress–strain relationship of monolayer BeO at 300 K along the armchair direction. Vacancy-defected (c) fracture stress, (d) fracture strain, and (e) elastic modulus behavior of monolayer BeO along the armchair direction.

Figure 10

(a) DV- and (b) SV-defected stress–strain relationship of monolayer BeO at 300 K along the zigzag direction. Vacancy-defected (c) fracture stress, (d) fracture strain, and (e) elastic modulus behavior of monolayer BeO along the zigzag direction.

The variations of mechanical behavior with different concentrations of vacancy are shown in Figures c,d and 10c,d. Due to the increase in vacancy concentration, both the fracture stress and strain decay monotonically. Along the armchair orientation, for a 1.0% DV concentration, the fracture stress is found as ∼52.1 ± 2.1 GPa, while for a 1.0% SV, it is ∼53.8 ± 1.8 GPa. On the other hand, in the zigzag orientation, a fracture stress of ∼56.76 ± 1.88 GPa and ∼52.74 ± 2.05 GPa is found for 1.0% DV and 1.0% SV concentrations, respectively. For both armchair and zigzag orientations, the variation of the elastic modulus with different considered vacancy concentrations is shown in Figures e and 10e. It is observed that for both orientations, the decreasing trend of elastic modulus is more substantial for the SV defect. At room temperature, due to the incorporation of 1% SV into the structure of monolayer BeO, the elastic modulus has decreased to ∼4.76 and ∼6.58% in the armchair and zigzag orientation, respectively, which is approximately 1.3–1.7 times greater than the decreasing trend of DV-induced BeO. Islam et al.[27] and Jing et al.[60] also reported an analogous trend of elastic modulus reduction for graphene and 2D-SiC with DV and SV defects. The impacts of DV and SV concentration on the % reduction of various mechanical properties calculated for monolayer BeO along the armchair and zigzag direction are listed in Tables and 2, respectively. Furthermore, we have also determined the SVBe- and SVO-induced fracture stress and elastic modulus at different concentrations and compared the results to the SVBe&O defect-based outcomes. Similar to the TC results, the SVBe&O defect shows the highest reduction of the fracture stress and elastic modulus for the same vacancy concentration. In contrast, the SVBe defect shows the smallest decreasing trend in calculated mechanical characteristics at different concentrations. The SVBe, SVO, and SVBe&O defect-induced tensile mechanical behavior at various concentrations are provided in Table S3 (for the armchair direction) and Table S4 (for the zigzag direction) (Supporting Information).

Table 1

Effect of DV and SV Concentration on the Percentage (%) Decrease of Different Tensile Mechanical Behavior of Monolayer BeO along the Armchair Direction

	armchair direction
	double vacancy			single vacancy
vacancy concentration (%)	fracture strength (%)	fracture strain (%)	elastic modulus (%)	fracture strength (%)	fracture strain (%)	elastic modulus (%)
1.0	20.6	30.1	2.7	18.1	23.6	4.7
2.0	25.3	30.8	5.6	33.7	38.9	9.6
3.0	26.1	32.1	8.1	36.3	39.1	14.4
4.0	26.9	33.2	10.9	39.1	39.8	18.7
5.0	30.3	33.4	14.2	43.8	41.6	22.2

Table 2

Effect of DV and SV Concentration on the Percentage (%) Decrease of Different Tensile Mechanical Behavior of Monolayer BeO along the Zigzag Direction

	zigzag direction
	double vacancy			single vacancy
vacancy concentration (%)	fracture strength (%)	fracture strain (%)	elastic modulus (%)	fracture strength (%)	fracture strain (%)	elastic modulus (%)
1.0	28.1	33.6	5.1	33.1	46.8	6.6
1.5	29.4	39.7	9.5	37.5	46.5	12.8
2.0	30.2	45.3	13.1	39.2	43.2	18.1
2.5	38.8	45.2	17.3	47.1	49.5	24.1
3.0	43.5	47.2	21.6	49.6	48.2	29.6

By means of MD simulations, the variation of the calculated mechanical properties of the monolayer BeO containing different types of defects can be easily and qualitatively understood from their initial calculated potential energy per atom.[5,61−63] At zero strain, the potential energy per atom shows a variation for the different types of structure. Generally, for a pristine structure, it should be the minimum. However, when the structure contains defects, the potential energy per atom increases compared to the pristine case. For an identical concentration, as the SV produces a stronger symmetry breakdown effect than the DV defect, the potential energy per atom at zero strain should show a higher value for the SV-defected structure. The evaluated potential energy per atom for pristine, 1% DV, and 1% SV defects along the armchair and zigzag orientations of monolayer BeO is shown in Figure a,b, respectively. From the results, it can be noticed that for both chiral directions, the pristine monolayer BeO shows a smaller initial potential energy value per atom compared to the vacancy-induced BeO. In addition, of the two types of vacancies considered in this study, the SV-induced BeO structure exhibits a greater initial potential energy per atom compared to the DV-induced BeO because of its greater regularity failure effect for a similar concentration.

Figure 11

Estimated potential energy per atom of pristine and vacancy-defected monolayer BeO along the (a) armchair and (b) zigzag directions. The zoomed-in view represents the initial energy variation of different sheets considered in this study.

Due to the unlikeness in bond orientations,[64] the bond elongating and bond breaking behaviors change considerably under uniaxial tensile force along the two chiral directions, producing a noticeably anisotropic fracture behavior both for pristine and for vacancy-induced systems. As revealed in Figure a, there are six force carrying bonds along the armchair alignment, among them four bonds have a relative angle of 60° with the force direction and the remaining two bonds are parallel to the force (armchair) direction. Oppositely, along the zigzag alignment, there are four force carrying bonds in every unit cell, and all the bonds have an angle of 30° with the force (zigzag) direction, as displayed in Figure b. For better understanding, the anisotropic lattice deformation owing to the uniaxial tensile force, we have computed the variation of bond length and bond angle as a function of strain along the two directions. The OVITO Pro bond analysis tool[45] is utilized to accomplish this analysis. Here, two bond lengths of types M and N and two bond angles of types δ and θ are used. When a uniaxial tensile force is employed along the armchair alignment, the modification in bond angles and bond lengths is indicated as armchair (δ or θ or M or N). In the same way, for tensile force along the zigzag alignment, the change in bond angles or bond lengths is specified as zigzag (δ or θ or M or N). We have used the room temperature tensile mechanical behavior to calculate the bond angles and bond length changes. From Figures and 10, the computed fracture stress (fracture strain) along the armchair and zigzag orientations, respectively, was 0.1519 (65.61 GPa) and 0.2215 (78.84 GPa) at 300 K. As shown in Figure d, the zigzag-directed bond (type N) stretches more than the armchair-directed bond (type M) with the enlargement of tensile force, and it is more likely to break down with higher fracture stress and strain along the zigzag orientation. In addition, the amplitude of the bond angle (shown in Figure c) along the zigzag orientation (zigzag (delta)) shows a significantly greater change than the armchair orientation (armchair (theta)) due to the application of tensile force. Thus, anisotropic fracture stress and elastic modulus are observed in armchair and zigzag orientations.

Figure 12

Pictorial view of bond angles (δ and θ) and bond lengths (M and N) of honeycomb BeO upon application of tensile force along the (a) armchair and (b) zigzag alignments. Change in (c) bond angle and (d) bond length as a function of strain when uniaxial tensile force is employed on monolayer BeO along the two chiral directions.

Conclusions

In conclusion, we investigated the phononic thermal transport and uniaxial tensile mechanical behavior of BeO nanosheets by utilizing MD simulations. Two types of vacancy defects, namely, the DV and SV with concentrations ranging from 1 to 5%, were considered to comprehensively analyze the TC and tensile deformation behavior of monolayer BeO. Moreover, the combined effect of vacancy (1%) and temperature (100 to 700 K) on the TC behavior of BeO nanosheets was also studied. Of the two considered vacancies, because of its remarkable regularity breakdown phenomena, the SV defect causes a greater reduction of the TC than the DV for all considered concentration. In order to qualitatively explain the attained TC and the phonon DOS for different lengths, different vacancies with varying concentrations were also explored. The phonon DOS peaks exhibit a strengthening nature with increasing length, quantifying the increasing nature of TC at an increased length. On the other hand, a softening nature of the phonon DOS peak intensities with an increasing concentration was found for both types of vacancies, and a substantial reduction of the peak intensities was found for the SV defect. Moreover, a significant phonon DOS reduction with an increasing temperature was obtained for both types of vacancy-defected BeO. During uniaxial tensile loading, owing to the bonding dissimilarity, the armchair-directed fracture stress, elastic modulus, and fracture strain show a greater value than the zigzag direction. For 1% SV, the reduction of the elastic modulus is 4.7 and 6.6% along the armchair and zigzag direction, respectively, whereas for 1% DV, the reduction was 2.7 and 5.1%. The greater reduction of mechanical strength due to the SV defect is elucidated by calculation of the potential energy per atom. Overall, the present results are beneficial for the understanding of the thermal transport and tensile mechanical behavior of monolayer BeO under the influence of vacancy defects and provide a potential way to tune the TC and mechanical strength for next-generation energy efficient nanodevices.