A study of size-dependent properties of MoS2 monolayer nanoflakes using density-functional theory

Novel physical phenomena emerge in ultra-small sized nanomaterials. We study the limiting small-size-dependent properties of MoS2 monolayer rhombic nanoflakes using density-functional theory on structures of size up to Mo35S70 (1.74 nm). We investigate the structural and electronic properties as functions of the lateral size of the nanoflakes, finding zigzag is the most stable edge configuration, and that increasing size is accompanied by greater stability. We also investigate passivation of the structures to explore realistic settings, finding increased HOMO-LUMO gaps and energetic stability. Understanding the size-dependent properties will inform efforts to engineer electronic structures at the nano-scale.

1 which contains the discrete symmetries: C 3 trigonal rotation, σ h reflection by the xy plane, σ v reflection by the yz plane, and all of their products 11 .
There have been significant efforts to understand the size-and edge-dependent, structural and electronic properties of MoS 2 monolayer nanoflakes. For example, quantum confinement effects in TMDC nanoflakes have been investigated by Miró et al., both experimentally and through density-functional theory (DFT) 7 . Wendumu et al. have presented the size-dependent optical properties of 1.6 to 10.4 nm MoS 2 nanoflakes 12 using the density-functional tight-binding (DFTB) method. An extensive DFT edge-dependence study on MoS 2 monolayer nanoribbons has been reported by Pan et al. 13 . Ellis et al. have studied the band gap tranistion in multilayered MoS 2 using DFT in gaussian09 with periodic boundary conditions 14 . Recently Nguyen et al. have experimentally studied the size-dependent properties of few-layer MoS 2 nanosheets and nanodots 15 but a complete study of the structural and electronic properties of very small single-layer MoS 2 nanoflakes has not yet been presented.
Here we report a DFT study of the 0 K size-dependent properties of 1H MoS 2 monolayers of size smaller than 2 nm. We begin our discussion by studying the relative stability of the armchair and zigzag configurations shown in Fig. 1. We present the geometries of the relaxed structures for different nanoflake sizes to thoroughly understand the structural response as a function of lateral size. We report the electronic properties: binding energy, flake formation energy, HOMO-LUMO (highest-occupied molecular orbital to lowest-unoccupied molecular orbital) gap, charge densities; and the passivation of the flakes. We are particularly interested in exploring how the HOMO-LUMO gap changes with the nanoflake size, leading to applications in HOMO-LUMO gap engineering. This is especially important as the HOMO-LUMO gap is the first step in determining the tunable fluorescent properties of nanoflakes, and as MoS 2 is known to be biocompatible 16 , nanoflakes of known size could be useful for biolabelling applications 17 .
This paper is organized as follows: first we discuss all the required methods and techniques. Then we study two different edge configurations for MoS 2 monolayers and find the most stable one, following with a discussion of structural stability as a function of size, the electronic properties and the properties of the passivated structures.

Methods
We investigated the structural and electronic properties of neutral MoS 2 monolayer nanoflakes with stoichiometry Mo n S 2n using DFT in gaussian09 18 . In experiments, usually triangular shaped islands of MoS 2 have been reported but it has been theoretically speculated that MoS 2 islands can exist in various shapes, such as trigonal, hexagonal, truncated hexagonal and rhombohedral [19][20][21][22] . We used rhombic flakes to maintain the neutrality and Mo n S 2n stoichiometry of the flakes. Also, we experienced convergence issues with triangular flakes.
To choose an appropriate functional for our modelling, we conducted an in-depth analysis of the functionals listed in Table 1. We picked a relaxed 72-atom flake as this was the largest size we could model with the B3LYP functional. We compared the relative atomic positions of each atom in the central zone of the 72-atom flake with the bulk structure (infinitely large and regular structure in all three dimensions) 23 . The displacement ΔR i of each atom from the bulk position is defined as where i indexes the atoms in the central zone of the 72-atom flake. The mean value of ΔR i , i.e., ΔR for each functional is given in Table 1. All functionals except B3LYP 24-26 result in less than 5% variation from the bulk atomic positions. This indicates that the three functionals, BHandHLYP 27 , PBE1PBE 28 , and M052X 29 predict similar structures at similar levels of accuracy. We also calculated the HOMO-LUMO gap as function of flake size for all these functionals as shown in Fig. 2. We expect the HOMO-LUMO gap to decrease with increasing flake size, approaching the experimental monolayer MoS 2 gap for larger flakes, as reported by Gan et al. 30 (Fig. 2). Hence, we do not consider these two functionals further. For smaller flakes, BH and HLYP and M052X both produce HOMO-LUMO gaps well above the monolayer experimental value 9 and we can expect the band gap with these functionals to converge close to the experimental monolayer band gap for larger flakes. Cramer and Truhlar report that M052X is not a recommended functional for transition metal chemistry 31 . Considering this, we therefore used the BHandHLYP functional for this article, although we have also performed all the calculations with M052X functional and did not find any major difference in the results. A  32 , ∆E x Becke88 is Becke's 1988 24 gradient correction to the local-spin density approximation (LSDA) for the exchange term, and E c LYP is the Lee-Yang-Parr correlation term 25 .
The basis set used was an effective-core potential basis set of double-zeta quality, the Los Alamos National Laboratory basis set also known as LANL2DZ 33 and developed by Hay and Wadt [34][35][36] . These basis sets are widely used in the study of quantum chemistry, particularly for heavy elements 33 .
gaussian09 optimization criteria: calculations were converged to less than 4.5 × 10 −3 Hartree/Bohr maximum force, 3 × 10 −4 Hartree/Bohr RMS force, 1.8 × 10 −3 Hartree maximum displacement, and 1.2 × 10 −3 Hartree RMS displacement. All the flakes were converged to the default SCF (self-consistent field) limit of <10 −8 RMS change in the density matrix except those specified in the next section. The charge multiplicity (net charge) was 0 and the spin multiplicity was 1 (singlet; spin neutral).
In the geometry optimization process, the geometry was modified until a stationary point on the potential surface was found. Analytic gradients were used and the optimization algorithm was the Berny algorithm using GEDIIS 37 . We calculated the electronic properties of the optimized structures. The charge densities were plotted in avogadro 38, 39 from a compatible gaussian09 checkpoint file. Data availability. The datasets generated and/or analysed during the current study are available from the corresponding author on reasonable request.
Size-dependent structural properties. The properties of MoS 2 monolayers are often investigated under the assumption of an infinite slab and real effects arising due to the confinement and boundaries are ignored. A nanoflake is a monolayer with spatial dimensions less than 100 nm. The structural, electronic and optical properties of such nanoflakes will be strongly influenced by varying their lateral size.
We study the MoS 2 monolayer nanoflakes for two commonly known edge structures, zigzag and armchair, to investigate the stable edge structure for smaller nanoflakes. Structures before geometry relaxation without any edge termination are shown in Fig. 1. Zigzag structures have double-coordinated, bridge-like S or Mo atoms on the edges [ Fig. 1(a)], whilst armchair have single-coordinated, antenna-like S or Mo atoms [ Fig. 1(b)]. We relaxed both of these types of structure, encountering convergence issues for the two larger structures (72 atoms and 105 atoms). We succeeded in getting convergence of <10 −7 RMS change in the density matrix for the 72-atom structures in both zigzag and armchair edge configurations. For the 105-atom structure, we obtained convergence of <10 −5 RMS change in the density matrix in zigzag edge configuration, but could not converge the 105-atom armchair edge configuration at all. This therefore, sets the maximum structure size in our calculations. In gaussian09, the energy change is not a criterion for convergence, however, the worst level of convergence for the largest structure, i.e., <10 −5 RMS change in density matrix, typically corresponds to <10 −10 Ha change in energy 18 . For the larger structures, we are more confident of the trends instead of the absolute values of energy.
The ground-state energies as functions of the size of the nanoflakes are shown in Fig. 3. Assuming that the edge width remains constant for any flake size, as the flakes get larger the ratio of number of edge atoms to core atoms decreases significantly because the number of core atoms increases more rapidly. (A quick circular approximation shows the core area ∝L 2 , whilst treating the edge as an annulus gives area ∝L, where L is the radius of the core.) The structure becomes more stable as it becomes larger. Figure 3 shows that zigzag-edged structure is more stable than the armchair configuration for nanoflakes of size less than 2 nm. All further properties are discussed for the zigzag edge configuration only because it is the stablest.
The relaxed structures of MoS 2 monolayer nanoflakes are shown in Fig. 4. We compared the atomic positions in the relaxed structures with their unrelaxed positions in the bulk structure 23 Fig. 4(d,e)]. Both of these structures show identical geometric behaviour and the maximum distortions are on the Mo-Mo lengths shown by red-arrowed lines d5 = d6 = 2.60 Å. These two structures show a well-established core whose mean structural parameters approach the bulk structure values 23 .
We have done an analysis of the Mo-S bond lengths in the central zones of our relaxed structures and compared them with the bulk Mo-S bond lengths of 2.41 ± 0.06 Å reported in 23 . Figure 4  n n flake 2 is the energy of the flake having n Mo atoms and 2n S atoms. As defined, FFE < 0 indicates that the flake is more stable than its constituent atoms. Figure 5 shows that with the increase in nanoflake size the FFE decreases sharply, so more energy is released by adding atoms in the  larger flakes indicating that the flakes tend to grow energetically. Conversely, more energy is required to break the larger flakes into their constituents.
We calculated the binding energies for all flake sizes and present them as a function of size in Fig. 6. We removed a Mo or S atom from as close as possible to the centre of the core or the edge as possible. The binding energy for the Mo atoms is given by Similarly, the binding energy for S atoms is given by B n n n n 2 2 1 S Negative values of the binding energy indicate that energy is required to remove an atom from a nanoflake. The negative dependence with size means that the cost rises with flake size. For example, removing a Mo atom from the core of a 45-atom flake requires ~1.2 eV more energy than removing it from the core of a 24-atom flake.
= − E E B D form , where E D form is the defect-formation energy so we can also calculate the energy required to create a Mo or S vacancy in the core or on the edge of the nanoflakes. From Fig. 6, significantly more energy is required to create a Mo vacancy as compared to a S vacancy. Also there is no major difference in the energy required to create a Mo vacancy in the core or in the edge in smaller flakes but as the size of the flakes increases, comparatively it becomes easier for defects to form on the edges. In case of S atoms, approximately the same energy is required to create a S vacancy in the core or in the edge as shown in Fig. 6(b).
To predict the electronic properties of ultra-small MoS 2 monolayer nanoflakes, we calculated their HOMO-LUMO gaps and charge densities of their HOMO and the LUMO (Fig. 7). With an increase in flake size, the HOMO-LUMO gap decreases for both unrelaxed and relaxed structures which is in keeping with intuition around the increase in the HOMO-LUMO gap with decreasing particle size as discussed in the methods section. Mak et al. 9 measured the band gap of 1.88 eV for a large MoS 2 monolayer sheet as shown by the dashed line in Fig. 7. For larger flakes, we have not observed the band gap converging to this value. One possible cause could be dangling bonds in the nanoflakes. To address this, we study passivated structures in the next section.
To get deeper insight into the HOMO-LUMO behaviour as a function of nanoflake size, we calculated charge-density plots (Fig. 7) for structures before and after the geometry relaxation. We can see that the majority of the HOMO and the LUMO charge densities are lying on the corners and edges in all of these structures except the 9-atom nanoflake where they are scattered over the whole structure. No single, stand-out trend is observed across all the structures. In short, the charge density is highly sensitive to the structural size for these small sized nanoflakes.   Table 2. A comparison of the mean lengths in the relaxed, passivated structures with and without H dimers on the Mo corner rings, encircled by green on all the structures in Fig. 8. There is a maximum mismatch of 2% in the S-S length in the 9 atom structure. To understand the effects of dangling bonds on the properties of the structures, we passivated both Mo and S edges with H atoms. We passivated each edge Mo atom with two H atoms as we expect Mo atoms to be bonded with six atoms in this particular MoS 2 stoichiometry. We also tested single H-termination of all edge Mo atoms and could not obtain converged, relaxed structures. We suspect this means that such structures are energetically unfavourable. We terminated each edge S atom with one H atom as all the central S atoms form three bonds with their neighbouring Mo atoms. We relaxed these passivated structures and observed that on the acute-Mo corner of all the nanoflakes, the H atoms are pushed away and they do not appear to bond to Mo atoms (Fig. 8). We investigated this non-bonding of corner Mo atoms with H atoms by checking their bond lengths. The average Mo-H bond length for all the edge Mo atoms is 1.665 ± 0.005 Å while on the corner it is 1.94 Å. The two H atoms on the  We removed the acute-Mo corner H atoms, relaxed the structures again and observed almost the same structural parameters on the corner as with the corner H atoms. We compared the mean Mo-Mo, S-S, and Mo-S lengths of the acute-Mo corner ring (encircled by green in Fig. 8) in Table 2 for the relaxed structures with and without the H dimer on the corner Mo atom. For all the structures, there is a minimal change in the bond lengths between 0-2%. All the S atoms bond well to one H atom each with an average S-H bond length of 1.365 ± 0.005 Å. We could not obtain a relaxed, converged 105-atom (we are not counting the number of H atoms to keep the number of atoms in each flake consistent with the previous discussion) passivated structure.
To calculate the stability, we have compared the energies of the passivated structures with the corresponding unpassivated ones. We found that the passivated structures are significantly more stable than the unpassivated ones by 4.33, 5.9, 6.96, and 9.66 eV for 9, 24, 45, and 72 atoms respectively as shown in Fig. 9 where the relative formation energy (RFE) is: n n m n n 2 2 2 where m is the number of H atoms in the passivated structures. Passivation of the dangling bonds modifies the electronic structure, charge densities, and hence the HOMO-LUMO gap. In Fig. 10, the HOMO-LUMO gap of the passivated structures is contrasted against the unpassivated ones. We find that the HOMO-LUMO gap widens with passivation. We suspect this is because of the removal of dangling bonds. This effect is significant in smaller nanoflakes but as the size increases, the ratio of edge to core atoms decreases. Hence, due to fewer edge states in the larger structures, the HOMO-LUMO gap difference (both relative and absolute) between the passivated and the unpassivated structures becomes smaller. The energy level diagram for the unpassivated and passivated flakes is shown in the Supplementary information.
The charge densities of the passivated structures are shown in Fig. 11. These are much more distributed states in contrast to the charge density plots for unpassivated, relaxed structures [ Fig. 7(b)]. Thus passivation makes HOMO/LUMO states in these small-sized flakes more like the expected infinite monolayer.

Conclusions
In summary, we have investigated the size-dependent structural and electronic properties of MoS 2 monolayer nanoflakes of sizes up to 2 nm using DFT. Our main focus has been to explore the small-sized nanoflakes. We provide more-detailed information for engineering small-sized nanoflakes by reporting the energetically favourable edge configuration and size of the nanoflakes. We predicted the trends in the energetics as functions of size. We passivated the structures to explore the effects of passivation on small-sized nanoflakes. We found the passivated structures to be more stable, with wider HOMO-LUMO gaps than unpassivated ones. We observe several strong size dependencies of various properties.
The size-dependence of the HOMO-LUMO gap of these small-sized nanoflakes holds promise for opto-electronic applications. However, due to the size-dependent energetics involved, one must take care in the manufacture/selection of these flakes. Due to limited computational resources, we were able to model only small-sized nanoflakes and can predict trends for larger flakes only by extending the fit functions. However, an extension of the current work to nanoflakes larger than 2 nm would be a good benchmark for the DFTB size-dependent HOMO-LUMO gaps reported by Wendumu et al. 12 .