Exploring the electron density localization in single MoS2 monolayers by means of a localize-electrons detector and the quantum theory of atoms in molecules

The nature of the electron density localization in a MoS2 monolayer under 0 % to 11% tensile strain has been systematically studied by means of a localized electron detector function and the Quantum Theory of atoms in molecules. At 10% tensile strain, this monolayer become metallic. It was found that for less than 6.5% of applied stress, the same atomic structure of the equilibrium geometry (0% strain) is maintained; while over 6.5% strain induces a transformation to a structure where the sulfur atoms placed on the top and bottom layer form S2 groups. The localized electron detector function shows the presence of zones of highly electron delocalization extending throughout the Mo central layer. For less than 10% tensile strain, these zones comprise the BCPs and the remainder CPs in separates regions of the space; while for the structures beyond 10% strain, all the critical points are involved in a region of highly delocalized electrons that extends throughout the material. This dissimilar electron localiz...


INTRODUCTION
A special attention has been paid to the single layer of the molybdenum sulfide, MoS 2 , in the recent years. 1 Upon thinning from the bulk, the electronic structure of MoS 2 undergoes an interesting transition. [1][2][3][4][5][6] Single layers, can be extracted using micromechanical cleavage, 7 or lithium-based intercalation. 8 While the MoS 2 bulk is an indirect band-gap semiconductor with a band-gap of ∼1.2 eV; the MoS 2 monolayer, MM, is a direct band-gap semiconductor with a band-gap of 1.8-1.9 eV. 9 MM exhibit enhanced electrical, optical and mechanical properties associated to their band gap. This gap can be tuned by means of doping, 10 external electric field 11 or mechanical strain. 1,[12][13][14][15][16] Studies have confirmed that applying strain is one of the best possible strategies to tune the MM band gap, since it neither attenuates the properties nor is inefficacious for single layers. 1,12 In general, the band gap monotonically decreases as the tensile strain increases. A 4% of tensile stress reduces de band gap to 0.83 eV and a stress larger than 10% leads to be metallic. 13 This has raised enormous interest in exploring the extraordinary properties of the MoS 2 mono layers.
On the other hand, in a not well-known paper published in 1964, Walter Kohn 17 showed for the first time that the insulating nature of matter is a consequence of electron localization in the ground state. This theory remained largely unexplored until it was revisited and generalized by Resta,18 who emphasized that the cause of the insulating behavior of matter is electron localization, defined in a suitable way: It is the electron (de)localization in the ground state that renders a system conducting or insulating. A remarkable link between the so-called localization tensor, LT, and bond orders described by a class of delocalization indices, DI, of chemical bonding theory was recently reported. 19 The chemistry of ground states and the physics of conductivity was linked in that way. This was achieved through a real space partition of LT according to quantum chemical topology. into intra-and interatomic contributions. Each atom or functional group in a system may be classified as essential, if its contribution to the LT diverges, or a spectator, if it converges, as electrical conductivity is regarded. This allows for the exploitation of traditional chemical intuition to identify essential atomic groups in determining electrical conductivity. 19 Otherwise, a localized-electrons detector (LED) function defined in the quantum mechanics local moment representation and given exclusively in terms of the electron density and its gradient was recently proposed. 20,21 Using this LED function, an opposite pattern of the electron density localization for typical insulator cubic ceramic and good conductors transition metal materials, corroborating the Kohn conclusion, was been recently reported. 22 It was found that, for the elemental ceramic materials, the zone of low electron localization is very small and show areas containing the BCPs and the remainder CPs in separates regions of the space: the electrons are localized in separate areas and cannot move freely through the material. On the contrary, for the transition metals, there are low-valued localized electrons detector isocontours defining a zone of highly delocalized electrons that extends throughout the material.
To further understand the nature of the MM transition from semiconductor to metallic behavior, a MoS 2 monolayer from its equilibrium geometric at 0 % to 11% tensile strain by means of the LED function and the Quantum Theory of atoms in molecules, QTAIM 23 were systematically studied. At 10% tensile strain, this monolayer become metallic. It was found that for less than 6.5% of applied stress, the same structure of the equilibrium geometry is maintained; while beyond 6.5% TS a transformation is induced, as consequence of the change of a cage in a bond CPs, to a structure where the sulfur atoms placed on the top and bottom layer form S 2 groups. The LED function showed the presence of zones of highly electron delocalization extending throughout the Mo central layer. For less than 10% tensile strain the corresponding structures show zones of highly electron delocalization containing the BCPs and the remainder CPs in separates regions of the space, while for the structure beyond 10% strain, the corresponding zone involves all the CPs. This dissimilar electron localization pattern is like to that reported for semiconductors, such as Si, Ge and Sn bulks and metallic systems, such as transition metals bulk. 22

Quantum theory of atoms in molecules
The quantum theory of atoms in molecules, QTAIM, 23 provides a quantitative link between the total one-electron density (regardless of how it is obtained: either through calculations or from experiment) and important physical properties of molecules or crystals. The topological properties of a molecule or crystal electron density distribution are summarized by their critical points (CPs). 23 These are points where the gradient vector field, ∇ρ(r) vanishes, and are classified by the ρ(r) curvatures, i.e., the three eigenvalues λ i (i=1, 2, and 3) of the Hessian matrix (H ij = ∂ 2 ρ(r)/∂x i ∂x j ). They are labeled by their rank (number of nonzero eigenvalues) and signatures (excess number of positive over negative eigenvalues) and correspond to maxima (3, 3) with λ 1 , λ 2 and λ 3 < 0, minima (3, +3) with λ 1 , λ 2 and λ 3 > 0, and saddle points (3, +1) with λ 1 > 0, λ 2 > 0 and λ 3 < 0 and (3, 1) with λ 1 < 0, λ 2 < 0 and λ 3 > 0. Every CP has a characteristic pattern of trajectories or gradient paths of ρ(r). Only the (3, 3) points are three-dimensional trajectory attractors: trajectories terminate at this CP. The (3, 3) points occur generally at the nuclear positions so that each nucleus is a 3D-point attractor in the vector field ∇ρ(r). The region traversed by the gradient paths which terminate at a given attractor, defines the basin of the 3D-attractor.
In general, a (3, 1) CP, is present between every pair of neighboring nuclei. It represents both local maxima (it is a 2D-attractor) in the directions of the two negative eigenvalues and a local minimum (1D-source) in the direction of the positive eigenvalue. The gradient paths associated with the negative eigenvalues terminate at this critical point, BCP, and define the interatomic surfaces, IAS, that partition the molecule into unique fragments (the atomic basins). Additionally, there is a unique pair of trajectories associated with the positive eigenvalue of BCP that originates at this point and terminate at the neighboring two nuclei (3D-attractors). These two unique gradient paths define a line linking the nuclei, the bond path. The network of bond paths defines the chemical structure and it is called the molecular graph. 23 The other CPs occur because of the geometrical arrangements of bond paths, 23 and they define the remaining elements of molecular structure: rings, and cages. If the bond paths are linked to form a ring of bonded atoms, then a (3, +1) critical point, RCP, is found in the interior of the ring. The (3, +3) CPs are called cage critical points, CCP, because they are located inside a cage nuclear arrangement. The electron density is a local minimum at a cage CP, and in crystals, it is the main source of trajectories of ∇ρ(r): trajectories only originate at such CPs and terminate at nuclei, bond or ring CPs.

The localized electrons detector, LED, function
Within the local quantum theory, any quantum observable A can be decomposed into two parts: a real partĀ(r) (the local value) and an imaginary partÃ(r) (the local spread). 24 The real part corresponds to the classical estimate to the observable A, whereas the imaginary part is its quantum fluctuation: For the momentum operatorP = ih − ∇, the local momentum representation 20 is: Thus, the quantum fluctuationP is proportional to ∇ρ(r)/ρ(r).P is an orbital-free direct threedimensional representation of localized pairs in atoms and molecules. 21 The single particle kinetic energy provides the electronic localization of electron pairs in molecules, and hence its associated momentumP is an ideal localized electron pair detector. 20,21,24 It is given in terms of the electron density and its gradient, The magnitude ofP depicts the regions where electrons are spatially confined in pairs, as in the case of bonding regions and atomic shells 21,24 and provides a direct three-dimensional representation of bonding interactions in molecules. LED isocontours show the same kind of symmetries 21 around the bond critical points that are given by the Laplacian of the electron density. Thus, it can show 3D-dimensionally the different kinds of bonding interactions identified by the topological analysis: the closed shell and the shared interactions. LED suggests local concentrations of charge (electron pairs) in regions around the ρ(r) critical points. 21 For the (3, -3) CPs which are 3D nuclear attractors (three negative curvatures) of the ρ(r) gradient field, LED shows high local charge concentration of ρ(r) on the inner core shells electrons. For the (3, -1) CP which are 2D attractor (two negative curvatures perpendicular to the bond) of ∇ρ(r), LED shows local charge concentration in bonding regions around the BCPs. The shape of the LED contour around these critical points are closed related to the ρ(r) curvatures. 20,21 Thus, in a typical covalent bond, |λ 3 | < |λ 1 | |λ 2 | and the associated LED contour has a cylindrical shape, whose main axis is aligned with the bond path. For a closed shell type bond with an appreciable ionic contribution, |λ 3 | > |λ 1 | ≥ |λ 2 | and the LED contours corresponds to discs perpendicular to the bond paths. 21 In general, LED contour shape can be uses as interaction identifiers and their value as an electron localization measure. 21 For example, Figure 1 shows the electron localization pattern for Au and Si bulks, a very good conductor and a semiconductor, respectively. As was found in ref. 22, for gold, most of the electron density is located inside big spheres (yellow zones in Fig. 1a with LED=0.6 au) around the nuclei core; while a low-valued (LED ≤ 0.3 au) localized-electron detector contour (Fig. 1b) defines a zone of highly delocalized electrons that extends throughout the material, occupying all the interstitial zone between the yellow isocontours (See Fig. 1c). An opposite pattern of the electron density localization for the semiconductor can be observed (See Fig 1d) where the electron density is mainly located around the nuclei (yellow areas) and the bond (green zones) CPs. For Si, a low-valued localizedelectron detector contour (LED=0.2 au) defining the most delocalized electron density is mainly located around two separated zones (Fig. 1e), one around the BCPs (green) and other around the RCP and CCPs (blue), occupying a very small space of the material. In general, for transition metals, TM, 22 most of the electron density, corresponding to LED>0.5 au, is located inside big spheres around the nuclei core. The interstitial zone between those spheres is occupied by highly delocalized electrons. For TM, the best conductors in nature, LED≤0.3 au define a zone of highly delocalized electrons that extends throughout the material and contain all the electron density bond, ring and cages critical points. For the semiconductors, 22 LED>0.5 au contours define separate zones around the nuclei and around the BCPs. For these materials, LED≤0.3 au contours defining the most delocalized electron density is also mainly located around two separated zones, one around the BCPs and other around the RCP and CCPs. Thus, for these materials, the electrons are in separate areas and cannot move freely through them and therefore it is not possible to form an electron current. LED suggests where the conduction channels in real space are.

CALCULATIONS
The lattices of bulk systems considered in this work were optimized (cell parameters and atomic coordinates) using Dmol 3 standard methods. For these periodic models, ρ(r) was calculated by means of the DMol 3 program 25,26 using the Kohn-Sham Hamiltonian with the gradient-corrected Perdew-Becke-Ernzerhof (PBE) exchange-correlation functional. 27 DMol 3 calculates variational self-consistent solutions to the DFT equations, expressing numerically the atomic orbital basis functions in an accurate spherical-polar mesh. The solutions to these DFT equations provide the molecular electron densities, which are used to evaluate the Laplacian and LED functions of the system. The numerical double-zeta plus polarization basis set DNP+ 25 was used in all calculations. The determination of the topology of ρ(r) for QTAIM analysis was performed using the AIMUC software. 28 The MARCHINCUBE, MC33 29,30 software was modified to build the LED contours.

RESULTS AND DISCUSSION
MoS 2 bulk is a layered type crystal whose lattice is described by the hexagonal space group P63/mmc with a = b = 3.160 Å and c = 12.294 Å. 31 Its crystal structure belongs to a family of polytypic structures with close-packed triangular double layers of S with Mo atoms arranged in the trigonal-prismatic holes of the S double layers. Each Mo atom is bonded to six S atoms in a trigonalprismatic arrangement. While the intralayer M-S bonding is strong, the inter-layer S-S bonding is weak as it arises from van der Waals forces. The space group of MM reduces to P6m2 due to loss of inversion symmetry. We have built a 3x3 cell of this monolayer separating a single layer from the bulk geometry. Vacuum layers thicker than 10 Å were used to ensure that there were no interactions between adjacent slabs. The geometry of this model was optimized using algorithms included in the Dmol3 program. 25 For this system, all the CPs were located and the parameters that characterize them are given in Table I (  the critical point). Additionally, systematic tensile strain on the cell parameters a and b in a range from 0% to 11% strain was studied and the CPs for these strained structures were also calculated. The CPs characteristic parameter for 2.5%, 6.5 %, 9.0% and 11% tensile strain are also reported in Table I. The data given in this Table shows that, like the elemental transition metals, 22,32 the Mo-S BCPs present characteristics associated with closed-shell type interactions: electron density at BCP, ρ b , small and the λ 3 curvature dominates, consequently ∇ 2 b > 0. These bonds contrast with the shared type interactions, typical of covalent bonds, in which ρ b is large, the λ 1 and λ 2 curvatures dominate and ∇ 2 b < 0. Typical ρ b values in shared interactions are 0.722, 0.551, and 0.252 au for N 2 , O 2 , and CC bond in ethane molecules, respectively. 23 Figure 2 shows the dependence of the band gap energies of MM with the tensile stress. As it can be seen in this Figure and in Table I, it has found that for less than 6.5% of applied stress, the same structure of the equilibrium geometry (that we call type A) is maintained: same number and type of CPs. The bond paths network (the monolayer graph) defining the monolayer geometry at 0% strain is shown in Figs. 3a and 3b. Almost the same type of CPs 30,33 as that of the bulk cases were found: Mo-S bond critical points (gray spheres in Fig. 3), four-membered ring CPs (blue spheres in Fig. 3), and trigonal prism like cage (red spheres in Fig. 3). Each cage is defined just by three nonplanar curved faces with a four-membered ring CP just at the center of each face. Additionally, for the monolayer case, six-membered ring CPs (light pink spheres in Fig. 3) were also found. Note in Figure 3 that the Mo-S bond paths, i.e., the trajectories joining the atoms, are not straight lines. The bond path is usually curved from the perimeter of strained rings. 23 The deviation of the length of the bond path from the internuclear distance is thus considered an indicator of the strain in a bond. 23,34 For the MoS 2 monolayer each bond path matches the internuclear line (See Fig. 3b) indicating no strain in the bond, however, the bond path is not a straight line. Curved bond paths have been reported in closed-shell type bonds such as hydrogen bonds, 35,36 agostic bond involving transition metals, 37,38 and metal alloys. [39][40][41] The basin, i.e., the volume spanned by the paths ending at each Mo atom is defined by six curved rhomboidal Mo-S IAS (See Fig. 3c), three with S atoms on the top layer and three with S atoms located on the bottom layer. Thus, the Mo atoms form a layer sandwiched (See Fig. 3d) between two layer of S atoms. Interestingly, the Mo basins reach the border of the monolayer. LED isocontours for this monolayer at 0 % strain are shown in Fig. 4. An isocontour with a high value of 1.0 au, meaning zones of high electron localization, is shown in Figs. 4a and 4b. This contour shows regions enclosed by light blue spherical and red contours contain the core and lone pairs of the S atoms, respectively. Yellow distorted big spheres, where most of the Mo electron density is located, can be also observed. Figs 4c and 4d, display a LED low-valued green isocontour equal to 0.3 au defining a zone of small electron localization (high electron delocalization). This green isocontour shows as separated regions oblate disks containing the BCPs of the Mo-S bonds and a shell-like region containing all the ring and cage critical points forming a shell of six-membered ring centered at the Mo nuclei (See Fig. 4c). This pattern of separated regions containing the BCPs is like  of that reported for the elemental semiconductors Si (See Fig. 1e), Ge and Sn bulks. 22 However, for the MoS 2 monolayer, the shell of the most delocalized electron is mainly found in the central layer of Mo atoms (See Fig. 4d). We have found that this structure (type A) is maintained under a tensile strain less than 6.5%, i.e., the monolayer graphs are isomorphic (same number and type of CPs) and the pattern of the most delocalized electrons for these structures is like to that reported for elemental semiconductor materials. The effect of the 10% tensile strain on the structure of the monolayer is shown in Figure 5. Clearly it can be seen in Figs. 5a and 5b, that this strain induces the transformation of the CCPs in BCPs, forming S-S bonds between the top and bottom sulfur layers. Thus, the sulfur atoms form week bonded S-S groups which assemble about the metals Mo atoms: three S 2 groups by each Mo atom. A similar structure (that we call type B) was found for tensile strain larger than 6.5% (See Fig. 2 and Table I). Similar change from CCPs to BCPs has been reported for equilibrium structures in metal alloys in which the enhanced of the alloys elastic constants was related with the BCPs formation in the alloy phase. [39][40][41] Despite this change in structure, the LED=0.3 au contour pattern for the monolayer with structure type B is like that of the type A (See Fig. 6), up to the TS reaches 10%, where the metallic behavior starts. Note in Fig. 6, the systematic manner that the LED zones around the BCPs (blue areas) coalesce toward the shell-like region containing all the ring and cage critical points (green areas). Observe that, even for a tensile strain of 9% involving a small gap of 0.2 eV, the blue and green areas are still separates. Figures 7 show the effect of the 10% tensile strain on the LED=0.3 au. Note that this zone of highly delocalized electrons now encompasses all CPs. Thus, 10% tensile strain induces a change in the zone of high electron delocalization of the structures containing in separates regions the BCPs and the remainder critical points (Fig. 6a), to a zone that involves all CPs (Fig. 6d), which is similar to that of good metallic conductor such as transition metals. 22 In MM this zone of high electron delocalization extends throughout the Mo central layer. Thus, the present results show that the metallic character of the MoS 2 monolayer produce by the tensile strain bigger than 10%, is associated with the formation of a region of highly delocalized electrons that extend along the central plane of Mo atoms, forming a kind of channel in which the electrons can easily move across the material.  A recent quantum conductance calculation work 5,42 on a mechanically deformed MM, using the nonequilibrium Green's functions method combined with the Landauer-Buttiker approach for ballistic transport, showed that at a critical tensile strain level close to 11%, the band gaps vanish to zero; the transport channels are completely open and current flow is enabled. Our results agree with the above findings; thus, Fig. 6 shows that the BCPs areas (blue) systematically merge to the shell-like region containing all the ring and cage critical points, as the TS increases. Just around of 11% TS, when the suggested 42 transport channels are open, the BCPs and the remainder CPs form a zone that involves all CPs (Fig. 6d), which is like that of good metallic conductor.

SUMMARY AND CONCLUSIONS
The structure of the MoS 2 monolayer defined by the electron density CPs, under tensile strain in a range from 0 to 10% strain was studied using QTAIM. For less than 6.5% of applied stress the same structure of the equilibrium geometry (0% strain) is maintained. In this structure, the Mo atoms form a layer sandwiched between two layer of S atoms, where each Mo atom is bonded by means of curved bond paths to three S atoms on the top layer and three S atoms located on the bottom layer. Strain larger than 6.5% induces the transformation of this structure to one where the sulfur second neighbors placed on the top and bottom layer form S-S bonds as consequence of the change of cage CPs in BCPs. In this new structure, the sulfur atoms form weak bonded S-S groups which assemble about the metals Mo atoms: three S 2 groups by each Mo atom. The electron density localization of those structures was studied using the LED function. For less than 10% tensile strain the corresponding structures show zones of highly electron delocalization containing the BCPs and the remainder CPs in separates regions of the space. This pattern is like semiconductors such as Si, Ge and Sn bulks. 20 On the contrary, the structure under 10 % tensile strain show a zone of highly delocalized electrons that encompass all CPs. This zone is like to that reported by metallic systems such as transition metals bulk. 20 Thus, 10% tensile strain induces the formation of a region of highly delocalized electrons that extend along the plane of Mo atoms. suggesting the formation of a kind of channel where the electrons can easily move across the monolayer.