Controlled edge dependent stacking of WS2-WS2 Homo- and WS2-WSe2 Hetero-structures: A Computational Study

Transition Metal Dichalcogenides (TMDs) are one of the most studied two-dimensional materials in the last 5–10 years due to their extremely interesting layer dependent properties. Despite the presence of vast research work on TMDs, the complex relation between the electro-chemical and physical properties make them the subject of further research. Our main objective is to provide a better insight into the electronic structure of TMDs. This will help us better understand the stability of the bilayer post growth homo/hetero products based on the various edge-termination, and different stacking of the two layers. In this regard, two Tungsten (W) based non-periodic chalcogenide flakes (sulfides and selenides) were considered. An in-depth analysis of their different edge termination and stacking arrangement was performed via Density Functional Theory method using VASP software. Our finding indicates the preference of chalcogenide (c-) terminated structures over the metal (m-) terminated structures for both homo and heterobilayers, and thus strongly suggests the nonexistence of the m-terminated TMDs bilayer products.

the ML WS 2 synthesis. Among various chemical techniques, lithium intercalation in WS 2 powder, followed by the easy chemical exfoliation of the MLs, is also noteworthy 26,27 .
Apart from the exfoliation methods, various other growth techniques such as Physical Vapor Deposition (PVD) [28][29][30] and Chemical Vapor Deposition (CVD) 31,32 techniques were adopted. Various CVD techniques using different substrates such as, on Si/SiO 2 , Ti/TiO 2 , Graphite, Graphene oxides, Sapphire, Au foil, h-BN etc. were also reported [33][34][35][36][37][38][39] . Among several existing growth techniques, epitaxial CVD growth mechanism is the most preferred one due to its several technological advantages [40][41][42] such as reduction of defect density, consistency in the overall growth, growth products with sharper interfaces, the perseverance of in-plane electrical conductivity via the introduction of mirror twin grain boundaries, and simultaneous reduction of tilt grain boundaries [43][44][45] . However, these experimental growth techniques require further control over the layer thickness, edge sharpness, size, and the quality of the as-grown crystals in order to achieve specific applicability 33 . Along with the ML TMDs, synthesis of van der Waals 2D homo and heterostructures are also becoming an attractive field of research due to their layer dependent vast tunable properties 6,46,47 . Therefore, it is essential to study the stable layered growth products and perform a thorough computational analysis in order to better understand the electronic and associated properties. In this respect, Density Functional Theory (DFT) study can provide deep insights of the edge and stacking dependency of the stable post growth products.
In this work, we have performed DFT calculations to get an in-depth understanding of the structural and electronic properties of homo/hetero bilayers of triangular WS 2, and WSe 2 flakes. We investigated the nature and strength of the interlayer interactions, by testing various combinations of AA and AB stacking. The most stable and unstable combinations were sorted based on their stacking energy and structural changes along with their electronic distribution. This comparative study also represents the important differences based on homo (WS 2 -WS 2 ) and hetero bilayers (WS 2 -WSe 2 ). This study aims to provide a crucial insight based on the existence of a specific combination of growth products towards the growth dynamics. Most of the existing studies on TMD heterostructures are based on periodic structures [48][49][50][51] . This work considers non-periodic structures, which are most commonly observed in experiments [52][53][54] . Our results provide insight into the stability of the post-growth TMD homobilayers and heterobilayers.

computational Methodology and Models
Methodology. All the electronic structure calculations were performed using the DFT method as implemented in the Vienna Ab initio Simulation Package (VASP) 55 . Projector augmented wave (PAW) pseudopotential is taken for the inert core electrons, and valence electrons are represented by plane-wave basis set 56,57 . The generalized gradient approximation (GGA), with the Perdew-Becke-Ernzerhoff (PBE) exchange-correlation functionals, is taken into account 58 . In order to accurately estimate weak van der Waals interaction, a vdW-correction approach is used. The vdW-density functional (vdW-DF) combines nonlocal correlations directly within a DFT functional. All DFT based calculations in this study have been performed using the optPBE functional within the vdW-DFT family, implemented in the VASP package by Klimes and co-workers [59][60][61] . The plane wave basis was set up with a kinetic energy cutoff of 400 eV. Due to the non-periodic nature of the structure, the Brillouin zone is sampled using the gamma-k-point grid for each mono and bilayers. The box size was kept big enough in order to prevent the periodic image overlap. The distance between two periodic images was maintained to be more than 8 Å on each direction (details are present in the Section S1 of supporting information (SI)). The initial distances (W-W interlayer) between the two layers were kept in between the ~6.25-6.30 Å following the previous studies 62, 63 . The convergence criterion for the electronic relaxation was kept 10 −5 eV/cell, and the total energy was calculated with the linear tetrahedron method with Blochl corrections. All the internal coordinates were relaxed using conjugate gradient methods until the Hellmann-Feynman forces are less than 0.02 eV/Å. In this study, all the atoms present in the systems were fully relaxed, and no atomic positions were fixed during the optimization method.
The free energy for stacking (stacking energy) was calculated using the following formula where, ΔG f signifies the stacking energy of the bilayer TMD, and ′ G G or G , BL ML ML represent the total free energies of the bilayers and different monolayers. In case of homo-bilayers, ML and ML′ represent the same mono-layers, but they represent different monolayers for the hetero-bilayers.

Models.
Our study considers only zigzag edged TMD structures due to their well-known energetic and as well as various electronic preferences over the armchair ones [64][65][66][67][68][69][70][71][72] . In this regard, triangular MX 2 flakes are known to be zigzag edged from each side (see Fig. 1c), and they can either be transition metal terminated (m-) or chalcogen terminated (c-). More importantly, triangular flakes are one of the well-known geometries formed by the MLs 73-75 . Fig. 1a represents the top view of the hexagonal honeycomb geometry of an MX 2 flake with both armchair and zigzag edges. Top and bottom side of the figure represents the armchair edges, while, the left and right sides are the representative of the differently terminated zigzag edges. To understand the influence of the different edge termination towards the bilayer formation, one has to begin the analysis by considering either fully metal or chalcogen terminated MX 2 monolayer. In this study, we have investigated the stacking of triangular ML flakes and have considered both the homo and the hetero-bilayers cases. In the case of homo-bilayers, we have considered WS 2 -WS 2 . For WS 2 -WSe 2 hetero-bilayer, we have changed the chalcogen atoms (X) of one of the layers to Se from S. Apart from the edge termination, we have also considered the possible contribution of different stacking of this monolayer MX 2 flakes towards the formation of the bilayers. In order to do so, we have only focused on the two most important stacking possibilities vis a vis. AA and AB (see Fig. 2). Our investigation is focused on the zigzag edged WS 2 and WSe 2 MLs with two different edge terminations -(a) metal (m-) erminated ( Fig. 1c1) and (b) chalcogen (c-) terminated (Fig. 1c2). Each of these MLs gets stacked together, forming bilayers (BLs) via two different stacking orientations named AA stacked (m-m and c-c), and AB stacked (m-m and c-c) as shown in Fig. 2. Here Fig. 2 is a general representation of all the homo/hetero layered (WS 2 -WS 2 (homo) and WS 2 -WSe 2 (hetero)) structures. AA is when the M and M (similarly X and X) sit on top of each other (Fig. 2a). On the other hand, AB represents the situation when M and X (similarly X and M) sit on top of each other (see Fig. 2b). It is important to mention here that due to the high computational cost and non-periodic nature of the models, small triangular flakes were considered to gain insight about the stacking phenomenon. Furthermore, the stacking of the triangular ML flakes was considered in order to fully understand the influence of edge termination of both layers and experimental reports on the stacking of triangular MLs were also present in the literature 76,77 . Since the main purpose is to understand the edge effect, and the second layer nucleates at the edge of the first layer instead of the basal plane 78 , we haven't considered the edge passivation.

Results and Discussion
Differently Stacked Homo/Hetero BLs Formation of WX 2 . After performing a rigorous scaling process (details are given in the section S1 of the SI), we have optimized the following formula for the c-and m-terminated MLs: a) W 10   www.nature.com/scientificreports www.nature.com/scientificreports/ and hetero-BL-c-AB), W 30 S 40 (homo-BL-m-AA and homo-BL-m-AB) and, W 30 S 20 Se 20 (hetero-BL-m-AA and hetero-BL-m-AB). We have relaxed the systems mentioned above to find stable orientations. Figure 3 represents the optimized structures for all the BLs, and Fig. S1 of the S1 represents the optimized structures of the above mentioned MLs. It is visually significant from Fig. 3 that for both the AA and AB stacked homo/hetero BLs, c-terminated structures are less distorted as compared to their m-terminated counterparts. However, it is quite evident from Fig. 3g,h that both the m-terminated AA/AB stacked heterobilayers are almost distorted from their original trigonal geometry.
To predict the stability of the afore-mentioned BL geometries, we have calculated the stacking energies Δ .
Moreover, a proper structural analysis was also conducted before studying their electronic properties. Tables 1 and 2 depicts all the stacking energies values for the homo and hetero BLs respectively. High negative stacking energies of all the c-terminated BLs are also indicative of their extreme stability. On the other hand, low negative stacking energies of the m-terminated homo/hetero BLs strongly suggests the lower probability of their formation during the growth process. Most importantly, a moderately high positive stacking energy (Table 2) of the m-terminated AA stacked WS 2 -WSe 2 BL implies the non-existence of the same during the growth process. We   www.nature.com/scientificreports www.nature.com/scientificreports/ also found that the c-terminated heterobilayers are more stable from their homo counterparts since they have higher negative stacking energies.
Structural analyses of these BLs and a comparative study with MLs is the first step to address the differences in the ΔG f values. Table 3 (c-terminated) and 4 (m-terminated) represents the average bond distances of all the BL and ML geometries. We have first considered the average interlayer W-W distances to get an insight of the same towards the stacking energy. In the case of c-terminated AA stacked BLs, both the homo and hetero W-W interlayer distances fall under the same range (~6.25-6.28 Å), which is comparable to the previously reported studies 62,63 as well (see Table 3). In the case of AB stacked structures, the W-W interlayer distances increase up to ~6.45 Å since the interlayer W's are not on top of each other due to their stacking arrangement. On the other hand, all the m-terminated geometries show increased W-W distances, and they are in the range of 6.56-6.84 Å (see Table 4).
It is important to note here that the increased interlayer W-W distances also correspond to their lower stacking energies as compared to the c-terminated geometries. Most importantly, the W-W distances are highest in the case of the AA stacked m-terminated heterobilayers. The increase in W-W distance weakens the van der Waals bond, which further results in lesser stability. The average intralayer W-W, W-S, and W-Se distances also play important roles in their corresponding energies. It is evident from Table 3 that all the intralayer bond distances (W-S/W-Se/W-W) are almost similar in the c-terminated BLs and MLs. Moreover, BLs are achieving the extra van der Waals attraction between the two layers and thus providing extreme stability. On the contrary, the intralayer W-S/W-Se distances of the m-terminated BLs vary a lot from their corresponding MLs and increased up to ~0.2-0.4 Å (see Table 4). The increase in the intralayer distance indicates the weakening of the bonds, which is in accordance with our findings. Furthermore, it is also evident from the optimized m-terminated MLs, that the intralayer bond length distribution varies a lot. A few bonds are shortened, and a few are elongated, resulting in a non-uniform distribution for the same. However, the interlayer bond distances are uniform in cases of all the c-terminated MLs and BLs.
Bader charge analysis of the BLs. Bader Charge Analysis 79 is one of the well-known approximation schemes to calculate the total electronic charge around each atom within the Bader volume. We have used post analysis scripts by Henkelman group 80 to calculate the total valence electron distribution around each atom in the system. In our case, the system contains W, S, and Se, and according to our used pseudopotential, each of these atoms contain 6 valence electrons. Thus, the charge is calculated based on a total of 6 electrons. All the c-and m-terminated bilayers were considered. Electron partitioning is tabulated with their corresponding numbers in section S3 of the SI. The net charge on atom can be calculated as follows: a) if an atom contains less than 6 electron and the number of electrons is x (x can be a fraction too), then that atom is positively charged and the corresponding charge is +(6-x) and, b) if an atom contains more than 6 electrons and the number of electrons is y (y can be a fraction too), then that atom is negatively charged and the corresponding charge is −(y-6). Figs. S2, S3 of the SI represent the side views of the c-and m-terminated homo/heterobilayers respectively. All atoms are numbered (this number corresponds to the numbered atoms in the section S3 of the SI) and differently oriented  www.nature.com/scientificreports www.nature.com/scientificreports/ such that most of them are visible. From the tabulated electron distribution (Section S3 of the SI), it is clear that all the W atoms are positively charged and all the X (S or Se) atoms are negatively charged. However, the distribution of electrons varies in a large extent for each of these systems. The variation is significant in between the m-and c-terminated systems for the homo and heterobilayers.
The highest and lowest charges on the W, S, and Se were calculated and marked with blue color inside the bracket (Section S3 of SI). The AA and AB stacked homobilayers show a uniform electron distribution for each W and S atoms, and it ranges in between ~3-4 electrons for W and ~6-7 electrons for S. Average electron partitions on W and S are in the range of ~3.7 and, ~6.8 respectively. On the other hand, the electron distributions show some amount of non-uniformity for the c-terminated AA and AB stacked heterobilayers. W attached to the S atoms show an average electron distribution in the range of ~3.7 and the W attached to the Se atoms show an average electron distribution in the range of ~4.1. As a result, the difference in electron distribution also significant in the S (higher electron distribution) and Se (lower electron distribution) atoms respectively. From our observation, it is clear that the difference in the electron distribution in heterobilayers (presence of two types of W atoms attached to different chalcogen (S/Se) atoms) as compared to the homobilayers acts as a stabilizing factor for the same. Our results suggest that the difference in electron distribution (charge separation) in between the W and X (S/Se) atoms is lower for the cases of m-termination and thus account for their lower negative stacking energies.

Charge density difference (CDD) analysis of the BLs.
To study the underlying reason behind the difference in stability of the considered BLs, charge density distribution (CDD) analyses were performed, and the spatial distribution of the electron cloud (isosurface) was plotted. Figures 4 and 5 describe the top and side views of the CDD of the homo and hetero bilayers respectively. The homogeneous charge density distribution for the c-terminated homobilayers (see Fig. 4a1,b1) and considerable charge density overlap between the bilayers (see Fig. 4a2,b2) are in accordance with their negative stacking energy. On the other hand, the diffused electron cloud distribution in the triangular edges of the m-terminated homobilayers (see Fig. 4c1,d1) and lesser overlap between the bilayers (see Fig. 4c2,b2) strongly support their lesser negative stacking energies as compared to their c-terminated analogues. However, the c-terminated heterobilayers show a strong electron cloud overlap between the two MLs (Fig. 5a2,b2) and the extent of charge distribution in both the AA and AB stacked heterobilayers are stronger as compared to the c-terminated AA and AB stacked homobilayers. These findings also corroborate their corresponding stacking energy from Section 3.1 as well. Nonetheless, Fig. 5d2 is suggestive of very less electron cloud overlap in between m-terminated AB stacked heterobilayer and therefore, validates its lower stacking energy. Most importantly, the m-terminated AA stacked heterobilayer (see Fig. 5c2) does not have any charge density overlap in between the two layers, and finally turns out to be an unstable system (positive stacking energy).

Density of states (DOS) analysis of the BLs.
To understand the changes in the electronic properties of these differently terminated bilayers, we have included the DOS analyses. The atomic DOS analyses give us an insight into the electron density distribution in or around the Fermi level, and the extent of atomic orbital overlap for each of these systems. Figures 6 and 7 describe the partial density of states (PDOS) plot of the homo/ heterobilayers. The AA/AB stacked c-terminated homo/heterobilayers (Figs. 6a,b, 7a,b) show a little discontinuity   www.nature.com/scientificreports www.nature.com/scientificreports/ and thus strongly suggests their extensive stability. In the case of heterobilayers, the extent of p-d overlap is much more significant and thus advocate towards their high stability.
From our analyses of the structural and electronic property, and the overall energetics, it clear that both c-terminated homo/hetero BLs are more likely to exist as bilayer 2D TMD growth products. The absence of structural deformity, homogeneous charge, and electron cloud distribution combine to form the most stable c-terminated geometries. Among m-terminated homo BLs, both the AA and AB stacking seems to be stable and are likely to exist, although, they are less stable than their c-terminated analogues. On the contrary, the AA stacked m-terminated BL is most likely to be non-existent. The instability or the lesser stability of the m-terminated cases can be due to the presence of dangling bonds (unsatisfied valence) of the heavy transition metal W. Due to W-termination, their valency is not satisfied fully and thus leading towards the structural instability. Generally, the passivation technique for the dangling bonds is used for multipurpose such as improving the quantum yield, preventing oxidation, increasing mobility, photoluminescence etc 81,82 . In this case, if we passivate and physically stabilize a TMD monolayer regardless of the type of termination, it will be very difficult for the bilayers to form/ exist. The reason is that the dangling bonds are required to grow the second TMD layer since the dangling bonds are known to act as the active sites for the epitaxial growth process 83 . The second layer grows from the edge of the first layer in the two-step growth method, as reported by the Ajayan group 78 . The edge passivation will inhibit or suppress the second layer formation. It should be noted that our work focuses on the existence of differently stacked bilayer post-growth products, and thus passivation will be far from the scope of this present study. Our study explores three different aspects such as (a) the effect of edge termination, (b) the effect of different stacking, and (c) the comparison of homo and hetero bilayers and draws an insightful picture for the probability of their existence as a bilayer TMD growth product. This study gives a motivation towards the controlled growth process of the homo and hetero layered TMDs and aims to provide a clear path towards the application-dependent TMD synthesis.

conclusion
In this DFT based study, electronic insights of the controlled bilayer TMD growth products with the changing parameters were studied. Three main parameters were considered such as different edge termination (c-and m-terminated), different bilayer stacking (AA and AB), and the change in the chalcogen atoms in one of the layers (homo and hetero). Tungsten (W) based chalcogenides (S and Se) were chosen due to their increasing demand in electronic applications over the well-studied Mo based systems. Any applications of these 2D Tungsten Chalcogenides (WX 2 ) require a well-established controlled synthesis path due to their structural dependent applications. Formation of bilayer flakes of these chalcogenides is a major research area and our calculation suggests the extreme stability of the c-terminated cases over the m-terminated bilayer growth products for both the