A Theoretical Investigation on the Physical Properties of Zirconium Trichalcogenides, ZrS 3 , ZrSe 3 and ZrTe 3 Monolayers

: In a recent advance, zirconium triselenide (ZrSe 3 ) nanosheets with anisotropic and strain-tunable excitonic response were experimentally fabricated. Motivated by the aforementioned progress, we conduct ﬁrst-principle calculations to explore the structural, dynamic, Raman response, electronic, single-layer exfoliation energies, and mechanical features of the ZrX 3 (X = S, Se, Te) monolayers. Acquired phonon dispersion relations reveal the dynamical stability of the ZrX 3 (X = S, Se, Te) monolayers. In order to isolate single-layer crystals from bulk counterparts, exfoliation energies of 0.32, 0.37, and 0.4 J/m 2 are predicted for the isolation of ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers, which are comparable to those of graphene. ZrS 3 and ZrSe 3 monolayers are found to be indirect gap semiconductors, with HSE06 band gaps of 1.93 and 1.01 eV, whereas the ZrTe 3 monolayer yields a metallic character. It is shown that the ZrX 3 nanosheets are relatively strong, but with highly anisotropic mechanical responses. This work provides a useful vision concerning the critical physical properties of ZrX 3 (X = S, Se, Te) nanosheets.


Introduction
After the experimental isolation of graphene [1][2][3] reported in 2004, two-dimensional (2D) nanomaterials have been extending continuously, owing to their exceptional physical and chemical properties. High surface-to-volume ratios in 2D systems can not only evolve to exceptional electronic and optical features, but are also highly appealing for practical chemistry-related applications, such as sensing and energy storage. In recent years, several 2D crystals with interesting physical properties have been fabricated, such as MoSi 2 N 4 [4,5], fluorinated diamane [6], penta-palladium phosphide selenide (PdPSe) [7], niobium oxide diiodide (NbOI 2 ) [8], penta-palladium phosphide sulfide (PdPS) [9], graphene-like BC 2 N [10], and penta-nickel diazenide (NiN 2 ) [11] nanosheets. Highly bright prospects for the application of 2D nanomaterials in critical technologies and their outstanding physical and chemical features act as a continuous driving force for experimental endeavors to design and fabricate novel crystals. In line with continuous experimental accomplishments in the field of 2D nanomaterials, most recently Li et al. [12] succeeded in the exfoliation of the zirconium triselenide (ZrSe 3 ) nanosheets. Experimental observations and theoretical calculations confirm highly anisotropic and strain-tunable semiconducting excitonic effects in ZrSe 3 nanosheets [12]. This experimental advance also highlights the appealing possibility of the exfoliation of zirconium trisulfide and tritelluride (ZrS 3 and ZrTe 3 ) Figure 1 depicts different views of the crystal structure of the ZrX 3 (X = S, Se, Te) monolayers. According to Figure 1, a ZrX 3 monolayer includes a rectangular primitive unit cell with a P21/M (No. 11) space group and highly anisotropic atomic arrangement along the x and y directions. A ZrX 3 monolayer can basically be considered laterally aligned (along the x direction) and alternatively inverted triangular prismatic ZrX 3 chains, which are connected through Zr-X bonds in the zx plane (l a Zr−X ). The inter-chain Zr-X bonds (l a Zr−X ) were found to be longer than the Zr-X bonds within the chain (l b Zr−X ). This observation was expected to lead to anisotropic mechanical, optical, and electronic properties in the ZrX 3 monolayers. Table 1 summarizes the optimized lattice constants, bond lengths, and calculated band gap properties of each ZrX 3 monolayer. Our predicted lattice constants for ZrX 3 monolayers agree well with previous data: a = 5.138 Å and b = 3.619 Å for ZrS 3 , a = 5.423 Å and b = 3.745 Å for ZrSe 3 , and a = 5.942 Å and b = 3.909 Å for ZrTe 3 [12]. For the bulk structures the box size along the out-of-plane direction for the ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers were found to be 9.01, 9.46, and 10.14 Å, respectively, with excellent agreement with corresponding experimental values of 8.98 [24], 9.43 [24], and 10.10 Å [24], respectively. In the Supplementary Materials document, the complete crystal data for the single-layer and bulk ZrX 3 systems are given. According to spin-polarized calculations, we found that these systems are not magnetic. Looking carefully to the crystal structure of a ZrX 3 monolayer, two different types of X atoms can be identified: surface X atoms, which exist in pairs (X 2 ), and the internal ones, which contribute to the inter-chain binding and only coordinate with the Zr atom. Our calculated bond lengths for X 2 were 2.07, 2.38, and 2.83 Å, which are only slightly different than the typical values for single covalent S-S, Se-Se, and Te-Te bonds [25,26]. According to Bader charge analysis, Zr atoms are positively charged and transfer electrons to the surface and internal chalcogen X atoms, indicating that electrostatic interactions play a key role in keeping the structural integrity of ZrX 3 monolayers. The amount of charge exchange decreased as the electronegativity of chalcogen  3 . Although all chalcogen atoms were negatively charged, it was found that surface X atoms yielded almost twice the charge as those of the internal counterparts. The electron localization pattern along the Zr-X bonds shown in Figure 1b also supports the idea of the dominance of ionic interaction within the ZrX 3 monolayers. The electron localization maps also show the gradual increase in electron localization in between adjacent surface X 2 moieties from ZrS 3 to ZrTe 3 monolayers. For the case of ZrTe 3 , an appreciable electron localization was found in between adjacent surface Te 2 moieties (find Figure 1d), indicating interactions between them. It is worth noting that our calculated distance between surface Te 2 moieties in the ZrTe 3 monolayer was 3.05 Å, which is only slightly larger than Te-Te bonds calculated for experimentally synthesized 2D αand β-Te: 3.02 Å [27]. Therefore, surface Te 2 moieties in the ZrTe 3 monolayer can actually be considered infinite Te 2 chains extended along the x direction. In other words, by comparing the ELF contours for the over-surface nonbonded X-X bonds, around the center of the Te-Te pairs (Figure 1d), the ELF values were considerably higher than those in the corresponding S-S ( Figure 1b) and Se-Se (Figure 1c) pairs, which reveals the formation of stronger Te-Te interactions in the ZrTe 3 monolayer. that electrostatic interactions play a key role in keeping the structural integrity of ZrX3 monolayers. The amount of charge exchange decreased as the electronegativity of chalcogen X atoms decreased: 1.77 e for ZrS3, 1.61 e for ZrSe3, and 1.33 e for ZrTes3. Although all chalcogen atoms were negatively charged, it was found that surface X atoms yielded almost twice the charge as those of the internal counterparts. The electron localization pattern along the Zr-X bonds shown in Figure 1b also supports the idea of the dominance of ionic interaction within the ZrX3 monolayers. The electron localization maps also show the gradual increase in electron localization in between adjacent surface X2 moieties from ZrS3 to ZrTe3 monolayers. For the case of ZrTe3, an appreciable electron localization was found in between adjacent surface Te2 moieties (find Figure 1d), indicating interactions between them. It is worth noting that our calculated distance between surface Te2 moieties in the ZrTe3 monolayer was 3.05 Å , which is only slightly larger than Te-Te bonds calculated for experimentally synthesized 2D α-and β-Te: 3.02 Å [27]. Therefore, surface Te2 moieties in the ZrTe3 monolayer can actually be considered infinite Te2 chains extended along the x direction. In other words, by comparing the ELF contours for the over-surface non-bonded X-X bonds, around the center of the Te-Te pairs (Figure 1d), the ELF values were considerably higher than those in the corresponding S-S ( Figure 1b) and Se-Se (Figure 1c) pairs, which reveals the formation of stronger Te-Te interactions in the ZrTe3 monolayer.  After an effective analysis of the structural and bonding characteristics of the ZrX 3 monolayers, we next examined their dynamic stability by evaluating phonon dispersion relations. The predicted phonon dispersion along highly symmetrical points for the ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers are illustrated in Figure 2. As the first important finding, the phonon modes were free of imaginary frequencies, confirming the dynamic stability of these systems. By increasing the atomic number of chalcogen atoms in the ZrX 3 nanosheets, it was clear that although preserving the general form of dispersions, the phonon modes showed narrower frequency ranges, which indicates lower group velocity. It is noticeable that all acoustic and optical modes appeared with considerable intersections, stimulating Energies 2022, 15, 5479 4 of 10 the scattering and reducing the lifetime for modes with higher frequencies. With generally lower group velocity and lifetime, it was expected that by increasing the atomic number of chalcogen atoms in ZrX 3 nanosheets, they would show lower lattice thermal conductivity, which is consistent with recent theoretical results for the ZrS 3 [13] and ZrTe 3 [14] systems. Table 1. Structural and electronic properties of ZrX 3 (X = S, Se, Te) monolayers.  Figure 1). 2 Bond lengths of those Zr-X bonds that are in the xz and yz planes, respectively ( Figure 1). 3 Bond length of covalently bonded X 2 moieties ( Figure 1). 4 Distance between adjacent X 2 moieties. 5 Average Bader charges on each of the Zr atoms, X atoms of X 2 moieties (X surface ), and internal X atoms (X internal ), respectively ( Figure 1). 6 Band gaps using PBE functional   Figure 1). 2 Bond lengths of those Zr-X bonds that are in the xz and yz planes, respectively ( Figure 1). 3 Bond length of covalently bonded X2 moieties ( Figure 1). 4 Distance between adjacent X2 moieties. 5 Average Bader charges on each of the Zr atoms, X atoms of X2 moieties (Xsurface), and internal X atoms (Xinternal), respectively ( Figure 1). 6 Band gaps using PBE functional ( )/PBE functional with inclusion of spin-orbit coupling (SOC) effect/ HSE06 functional ( 06 )/HSE06 functional with inclusion of SOC effects ( 06+ ). For ZrS3 and ZrSe3, K-points at which valance band maximum and conduction band minimum occurs are shown in parentheses.
After an effective analysis of the structural and bonding characteristics of the ZrX3 monolayers, we next examined their dynamic stability by evaluating phonon dispersion relations. The predicted phonon dispersion along highly symmetrical points for the ZrS3, ZrSe3, and ZrTe3 monolayers are illustrated in Figure 2. As the first important finding, the phonon modes were free of imaginary frequencies, confirming the dynamic stability of these systems. By increasing the atomic number of chalcogen atoms in the ZrX3 nanosheets, it was clear that although preserving the general form of dispersions, the phonon modes showed narrower frequency ranges, which indicates lower group velocity. It is noticeable that all acoustic and optical modes appeared with considerable intersections, stimulating the scattering and reducing the lifetime for modes with higher frequencies. With generally lower group velocity and lifetime, it was expected that by increasing the atomic number of chalcogen atoms in ZrX3 nanosheets, they would show lower lattice thermal conductivity, which is consistent with recent theoretical results for the ZrS3 [13] and ZrTe3 [14] systems. The Raman spectrum of each ZrX3 monolayer is shown in Figure 3. Apparently, due to the same crystal symmetry of each monolayer ZrX3, the number of Raman active modes was the same in all structures. It is seen that there existed three prominent Raman active phonon modes for each monolayer structure. The three prominent peaks are labeled I, II, and III, as shown in the figure. As the chalcogenide atom changed from S to Te, which means the atomic radius increased, the frequency of each prominent peak displayed red shift due to the larger atomic mass. In addition, the motion of individual atoms showing the vibration of the corresponding phonon mode is also given in the right panel of the figure. The phonon mode I was calculated to be at the frequencies 522, 292, and 151 cm −1 for ZrS3, ZrSe3, and ZrTe3 monolayers, respectively. Apparently, mode I arose from the in- The Raman spectrum of each ZrX 3 monolayer is shown in Figure 3. Apparently, due to the same crystal symmetry of each monolayer ZrX 3 , the number of Raman active modes was the same in all structures. It is seen that there existed three prominent Raman active phonon modes for each monolayer structure. The three prominent peaks are labeled I, II, and III, as shown in the figure. As the chalcogenide atom changed from S to Te, which means the atomic radius increased, the frequency of each prominent peak displayed red shift due to the larger atomic mass. In addition, the motion of individual atoms showing the vibration of the corresponding phonon mode is also given in the right panel of the figure. The phonon mode I was calculated to be at the frequencies 522, 292, and 151 cm −1 for ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers, respectively. Apparently, mode I arose from the in-plane vibration of the outermost chalcogenide atoms against each other. On the other hand, modes II and III stemmed from the out-of-plane vibration of the Zr and chalcogenide atoms. The frequencies of modes II/III were calculated to be 321/274, 220/171, and 74/62 cm −1 for the ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers, respectively. The mode II phonon represents the out-of-plane vibration of the internal chalcogenide atoms against Zr atoms, whereas the outermost chalcogenide atoms made little contribution to the vibration. Finally, the mode III phonon arose from the out-of-phase vibration of each atomic level with respect to each other, that is, the outermost chalcogenide atoms vibrated against each other, and Zr atoms also moved out of phase. The three characteristic phonon peaks were quite important for the detection of the ZrX 3 monolayers. represents the out-of-plane vibration of the internal chalcogenide atoms against Zr atoms, whereas the outermost chalcogenide atoms made little contribution to the vibration. Finally, the mode III phonon arose from the out-of-phase vibration of each atomic level with respect to each other, that is, the outermost chalcogenide atoms vibrated against each other, and Zr atoms also moved out of phase. The three characteristic phonon peaks were quite important for the detection of the ZrX3 monolayers. Worth mentioning is that to prepare 2D materials, the common approaches include mechanical exfoliation, chemical vapor deposition (CVD), and liquid-phase exfoliation. The mechanical exfoliation of bulk 2D materials into single or multiple layers can be achieved by using external driving forces [28][29][30]. Before analyzing the electronic properties, it is thus very useful to investigate the mechanical exfoliation energy required for the isolation of the ZrS3, ZrSe3, and ZrTe3 monolayers from their native bulk structures. For this purpose, we first acquired the energy-minimized six-layer slabs of the ZrX3 nanosheets, with the same stacking pattern as that of their bulk systems. In the next step, the last layer was steadily separated toward the out-of-plane vacuum direction, with a small step of 0.25 Å . The change in the energy of the systems was subsequently calculated and the cleavage energy was recorded. As shown in Figure 4, the relative energies showed sharp initial increases and later reach converged values. According to our DFT-D3 simulations, the exfoliation energies of 0.32, 0.37, and 0.40 J/m 2 were predicted for the isolation of the ZrS3, ZrSe3, and ZrTe3 monolayers, which are comparable to that of graphene: 0.37 J/m 2 [31]. These findings reveal that the separate layers in these systems showed relatively weak interactions and moreover highlight that by increasing the atomic number of chalcogen atoms, the exfoliation energy increased, consistent with earlier studies [32]. We remind that Li et al. [12] synthesized ZrSe3 nanosheets via the mechanical exfoliation method. Taking into account our predictions for the exfoliation energies and the aforementioned experimental achievement, the experimental isolation of ZrS3 and ZrTe3 monolayers from their bulk structures is very bright. Worth mentioning is that to prepare 2D materials, the common approaches include mechanical exfoliation, chemical vapor deposition (CVD), and liquid-phase exfoliation. The mechanical exfoliation of bulk 2D materials into single or multiple layers can be achieved by using external driving forces [28][29][30]. Before analyzing the electronic properties, it is thus very useful to investigate the mechanical exfoliation energy required for the isolation of the ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers from their native bulk structures. For this purpose, we first acquired the energy-minimized six-layer slabs of the ZrX 3 nanosheets, with the same stacking pattern as that of their bulk systems. In the next step, the last layer was steadily separated toward the out-of-plane vacuum direction, with a small step of 0.25 Å. The change in the energy of the systems was subsequently calculated and the cleavage energy was recorded. As shown in Figure 4, the relative energies showed sharp initial increases and later reach converged values. According to our DFT-D3 simulations, the exfoliation energies of 0.32, 0.37, and 0.40 J/m 2 were predicted for the isolation of the ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers, which are comparable to that of graphene: 0.37 J/m 2 [31]. These findings reveal that the separate layers in these systems showed relatively weak interactions and moreover highlight that by increasing the atomic number of chalcogen atoms, the exfoliation energy increased, consistent with earlier studies [32]. We remind that Li et al. [12] synthesized ZrSe 3 nanosheets via the mechanical exfoliation method. Taking into account our predictions for the exfoliation energies and the aforementioned experimental achievement, the experimental isolation of ZrS 3 and ZrTe 3 monolayers from their bulk structures is very bright.
To explore the electronic characteristics of the ZrX 3 monolayers, we calculated electronic band structures using PBE and the more accurate HSE06 functional. The effect of spin-orbit coupling (SOC) on band gap properties of ZrX 3 monolayers was also examined. Figure 5 depicts the HSE06 band structures of the ZrX 3 monolayers without and with the inclusion of SOC effects. The corresponding PBE results are also given in Figure  S1 in the Supplementary Materials. Band-gap values as well as transition k-points for HSE06 results are listed in Table 1. According to results shown in Figure 5 and Table 1, the ZrS 3 and ZrSe 3 monolayers were indirect gap semiconductors with HSE06 band gaps of 1.93 and 1.01 eV. In the ZrS 3 monolayer, VBM was located at a k-point on the Γ-X path (0.184210 k 1 + 0 k 2 , k 1 and k 2 are lattice vectors of reciprocal space) and CBM was located at the Γ point. The direct gap at the Γ point in ZrS 3 was only 0.08 eV larger than the indirect Energies 2022, 15, 5479 6 of 10 gap, indicating that ZrS 3 may behave as a quasi-direct gap semiconductor. VBM and CBM in the ZrSe 3 monolayers , were, however, located at the Γ and X points, respectively. Both band gap values and transition k-points are in good agreement with previous data (1.92 eV for ZrS 3 and 0.92 eV for ZrSe 3 ) [33]. It can be seen that for each of the ZrS 3 and ZrSe 3 monolayers, PBE and HSE06 band structures looked similar, except the PBE band gaps, as expected, were underestimated. It was also observed that the inclusion of SOC did not yield a detectable effect on the band structure of the ZrSe 3 , but it reduced the HSE06 (PBE) band gap of the ZrSe 3 monolayer by 0.08 (0.06) eV (find Figures 5 and S1). Unlike the ZrS 3 and ZrSe 3 monolayers, the ZrTe 3 monolayer exhibited a metallic character irrespective of the functional used, which is also in agreement with a previous report [33].   Figure S1 in the Supple Materials. Band-gap values as well as transition k-points for HSE06 results are listed 1. According to results shown in Figure 5 and Table 1, the ZrS3 and ZrSe3 monolay indirect gap semiconductors with HSE06 band gaps of 1.93 and 1.01 eV. In the ZrS3 m VBM was located at a k-point on the Γ-X path (0.184210 1 + 0 2 , 1 and 2 are la tors of reciprocal space) and CBM was located at the Γ point. The direct gap at the Γ ZrS3 was only 0.08 eV larger than the indirect gap, indicating that ZrS3 may behave a direct gap semiconductor. VBM and CBM in the ZrSe3 monolayers, were, however, l the Γ and X points, respectively. Both band gap values and transition k-points are agreement with previous data (1.92 eV for ZrS3 and 0.92 eV for ZrSe3) [33]. It can be for each of the ZrS3 and ZrSe3 monolayers, PBE and HSE06 band structures looked except the PBE band gaps, as expected, were underestimated. It was also observed inclusion of SOC did not yield a detectable effect on the band structure of the ZrS reduced the HSE06 (PBE) band gap of the ZrSe3 monolayer by 0.08 (0.06) eV (find and S1). Unlike the ZrS3 and ZrSe3 monolayers, the ZrTe3 monolayer exhibited a meta acter irrespective of the functional used, which is also in agreement with a previo In order to further rationalize the band-gap change trend in ZrX 3 monolayers, for each one, we calculated atom-type projected density of states (PDOS) and charge density distributions of valance band maximum (VBM) and conduction band minimum (CBM), as shown in Figure 5. For the three monolayers, VBM was mainly contributed by internal X-(p y ,p z ) with a minor contribution from Zr-(p y ,d yz ,d xy ), representing a shallow bonding s(Zr-X) state. Considering the fact that the energy position of p orbitals of chalcogen X atom increased from S to Te, VBM in the ZrX 3 monolayers was expected to move upward in energy from the ZrS 3 to ZrTe 3 monolayer. This expectation was confirmed by the absolute energy positions of VBMs of the ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers, calculated to be −6.57, −5.68, and −5.21 eV, respectively. CBM in the ZrS 3 lattice was made of Zr-(s,d z 2 ,d x 2 −y 2 ), representing a bonding (Zr-Zr) state propagating along the y direction. The charge densities of CBMs in the ZrSe 3 and ZrTe 3 monolayers were, however, almost exclusively distributed over surface X 2 moieties, and they both represented bonding (X-X) states made of X-(s,p x ,p z ). As our electron localization analysis revealed, the strengths of interaction between adjacent X 2 moiety increased from ZrS 3 to ZrTe 3 . The stronger the interaction between X 2 moieties, the lower the absolute energy position. Putting the conclusions together, from ZrS 3 to ZrTe 3 , VBM shifted upward in energy whereas CBM shifted downward, leading to smaller band gaps (E g (ZrS 3 ) > E g (ZrSe 3 ) > E g (ZrTe 3 )). In order to further rationalize the band-gap change trend in ZrX3 monolayers, for each one, we calculated atom-type projected density of states (PDOS) and charge density distributions of valance band maximum (VBM) and conduction band minimum (CBM), as shown in Figure 5. For the three monolayers, VBM was mainly contributed by internal X-(py,pz) with a minor contribution from Zr-(py,dyz,dxy), representing a shallow bonding s(Zr-X) state. Considering the fact that the energy position of p orbitals of chalcogen X atom increased from S to Te, VBM in the ZrX3 monolayers was expected to move upward in energy from the ZrS3 to ZrTe3 monolayer. This expectation was confirmed by the absolute energy positions of VBMs of the ZrS3, ZrSe3, and ZrTe3 monolayers, calculated to be −6.57, −5.68, and −5.21 eV, respectively. CBM in the ZrS3 lattice was made of Zr-(s, 2 , 2 − 2 ), representing a bonding (Zr-Zr) state propagating along the y direction. The charge densities of CBMs in the ZrSe3 and ZrTe3 monolayers were, however, almost exclusively distributed over surface X2 moieties, and they both represented bonding (X-X) states made of X-(s,px,pz). As our electron localization analysis revealed, the strengths of interaction between adjacent X2 moiety increased from ZrS3 to ZrTe3. The stronger the interaction between X2 moieties, the lower the absolute energy position. Putting the conclusions together, from ZrS3 to ZrTe3, VBM shifted upward in energy whereas CBM shifted downward, leading to smaller band gaps (Eg(ZrS3)> Eg(ZrSe3)> Eg(ZrTe3)).
Finally, we examined the mechanical responses by performing uniaxial tensile simulations along the x and y directions, as distinguished in Figure 1. Uniaxial stress-strain responses of the ZrX3 monolayers along the x and y directions are illustrated in Figure 6. In these results, real volumes of the deformed lattices were considered in the conversion of the stress values to the standard GPa unit [34][35][36]. The real area of the deformed Finally, we examined the mechanical responses by performing uniaxial tensile simulations along the x and y directions, as distinguished in Figure 1. Uniaxial stress-strain responses of the ZrX 3 monolayers along the x and y directions are illustrated in Figure 6. In these results, real volumes of the deformed lattices were considered in the conversion of the stress values to the standard GPa unit [34][35][36]. The real area of the deformed nanosheets can be obtained using the simulation box sizes along the in-plane directions. The real volume at every strain was calculated by finding the normal distance between boundary chalcogen atoms plus their effective vdW diameter. According to our geometry-optimized bulk lattices, the thicknesses of the stress-free ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers were predicted to be 9.01, 9.46, and 10.14 Å, respectively, equivalent to the effective vdW dimeters of 3.05, 3.20, and 3.26 Å for the S, Se, and Te atoms in the ZrX 3 monolayers, respectively. The stress-strain curves plotted in Figure 6 are uniaxial, which means that during the complete deformation and after the geometry minimization, these kagome monolayers exhibited a stress component only along the loading direction and showed negligible values along the two other perpendicular directions. As expected, and stemming from the anisotropic structure, the stress-strain curves along the x and y directions were considerably different, confirming highly anisotropic mechanical features. The elastic modulus of ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers along the x (y) directions were predicted to be 93 (142), 90 (118), and 120(56) GPa, respectively. The ultimate tensile strength of the ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers along the x (y) directions were predicted to be 6.2 (16.5), 4.9 (17.4), and 4.6 (11.8) GPa, respectively. As expected, due to the existence of more connecting bonds along the y direction than the x counterpart, these systems showed considerably higher tensile strength along this direction. The same observation is also consistent for the elastic modulus of the ZrS 3 and ZrSe 3 monolayers. The ZrSe 3 monolayers unexpectedly showed a higher elastic modulus along the x direction, which, as discussed earlier with the ELF results, was due to the formation of continuous over-surface Te-Te interactions along this direction in this system. As is clear, due to the multifaceted structural and bonding effects, the ZrX 3 nanosheets showed highly anisotropic and complex mechanical behavior. Worthy to mention that complex material properties can be explored using MTPs with high accuracy and accelerated computational costs [37][38][39][40][41].
were considerably different, confirming highly anisotropic mechanical features. The elastic modulus of ZrS3, ZrSe3, and ZrTe3 monolayers along the x (y) directions were predicted to be 93 (142), 90 (118), and 120(56) GPa, respectively. The ultimate tensile strength of the ZrS3, ZrSe3, and ZrTe3 monolayers along the x (y) directions were predicted to be 6.2 (16.5), 4.9 (17.4), and 4.6 (11.8) GPa, respectively. As expected, due to the existence of more connecting bonds along the y direction than the x counterpart, these systems showed considerably higher tensile strength along this direction. The same observation is also consistent for the elastic modulus of the ZrS3 and ZrSe3 monolayers. The ZrSe3 monolayers unexpectedly showed a higher elastic modulus along the x direction, which, as discussed earlier with the ELF results, was due to the formation of continuous over-surface Te-Te interactions along this direction in this system. As is clear, due to the multifaceted structural and bonding effects, the ZrX3 nanosheets showed highly anisotropic and complex mechanical behavior. Worthy to mention that complex material properties can be explored using MTPs with high accuracy and accelerated computational costs [37][38][39][40][41].

Concluding Remarks
We studied the structural, phononic, electronic, and single-layer exfoliation energies and mechanical properties of ZrX3 (X = S, Se, Te) monolayers. The acquired phonon dispersion relations revealed the dynamical stability of the aforementioned 2D systems. Exfoliation energies of 0.32, 0.37, and 0.40 J/m 2 were predicted for the ZrS3, ZrSe3, and ZrTe3 monolayers' isolation, which confirms bright prospects for the mechanical isolation of ZrS3 and ZrTe3 monolayers from their native bulk structures. ZrS3 and ZrSe3 monolayers were found to be indirect gap semiconductors, with HSE06 band gaps of 1.93 and 1.01 eV, whereas the ZrTe3 monolayer yielded a metallic character. The elastic modulus of ZrS3, ZrSe3, and ZrTe3 monolayers along the x (y) directions were predicted to be 93 (142), 90

Concluding Remarks
We studied the structural, phononic, electronic, and single-layer exfoliation energies and mechanical properties of ZrX 3 (X = S, Se, Te) monolayers. The acquired phonon dispersion relations revealed the dynamical stability of the aforementioned 2D systems. Exfoliation energies of 0.32, 0.37, and 0.40 J/m 2 were predicted for the ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers' isolation, which confirms bright prospects for the mechanical isolation of ZrS 3 and ZrTe 3 monolayers from their native bulk structures. ZrS 3 and ZrSe 3 monolayers were found to be indirect gap semiconductors, with HSE06 band gaps of 1.93 and 1.01 eV, whereas the ZrTe 3 monolayer yielded a metallic character. The elastic modulus of ZrS 3 , ZrSe 3 , and ZrTe 3 monolayers along the x (y) directions were predicted to be 93 (142), 90 (118), and 120(56) GPa, respectively, and the corresponding ultimate tensile strength values were found to be 6.2 (16.5), 4.9 (17.4), and 4.6 (11.8) GPa, respectively. It is shown that because of multifaceted structural and bonding effects, ZrX 3 nanosheets showed highly anisotropic and complex mechanical behavior. The presented DFT results provide an effective overview of the key physical properties of the ZrX 3 (X = S, Se, Te) nanosheets, which can serve as valuable information for their practical application in nanodevices.