Adsorption, intercalation, diffusion, and adhesion of Cu at the 2Hâ‹TMMoS2 (0001) surface from first-principles calculations

Study of the adsorption of a transition metal on the surface of a layered material and the possible subsequent intercalation into that layered material is of fundamental interest and potential technological importance. In the present work, we choose the transition metal Cu as the adsorbate or intercalant and 2H-MoS2 as the layered material. Energetics are calculated characterizing four of the most basic surface and interfacial phenomena: adsorption, intercalation, diffusion, and adhesion. Using first-principles density functional theory (DFT), we find that intercalating a Cu atom into the van der Waals (vdW) gap below the MoS2 (0001) surface is 0.665 eV more favorable than adsorbing the Cu atom on top of the surface, i.e., intercalation of single Cu atoms is strongly favored thermodynamically. Also, we find that the system with adsorbed Cu is magnetic, while the system becomes nonmagnetic after the Cu atom is intercalated into the vdW gap. We obtain the diffusion barriers of the Cu atom on the surface and in the vdW gap to be 0.23 and 0.32 eV, respectively. We also obtain an adhesion energy of 0.874 J/m2 for a Cu (111) slab bonding with a 2H-MoS2 (0001) slab. The DFT value of adhesion energy, as well as the gap width at the interface between the Cu and MoS2, depends strongly on the choice of the functional. From our analysis on bulk properties of both MoS2 and Cu, we suggest that our vdW-DF2-B86R results listed above are reliable for applications, e.g., interpretation of experimental results and physical modeling of this adsorption system. Disciplines Biological and Chemical Physics | Materials Chemistry | Physical Chemistry Comments This article is published as Han, Yong, Michael C. Tringides, James W. Evans, and Patricia A. Thiel. "Adsorption, intercalation, diffusion, and adhesion of Cu at the 2H−MoS2 (0001) surface from firstprinciples calculations." Physical Review Research 2, no. 1 (2020): 013182. DOI: 10.1103/ PhysRevResearch.2.013182. Posted with permission. Creative Commons License This work is licensed under a Creative Commons Attribution 4.0 License. This article is available at Iowa State University Digital Repository: https://lib.dr.iastate.edu/chem_pubs/1210 PHYSICAL REVIEW RESEARCH 2, 013182 (2020) Adsorption, intercalation, diffusion, and adhesion of Cu at the 2H-MoS2 (0001) surface from first-principles calculations Yong Han ,1,2,* Michael C. Tringides,1,2 James W. Evans ,1,2 and Patricia A. Thiel1,3,4 1Ames Laboratory, U.S. Department of Energy, Ames, Iowa 50011, USA 2Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA 3Department of Chemistry, Iowa State University, Ames, Iowa 50011, USA 4Department of Materials Science and Engineering, Iowa State University, Ames, Iowa 50011, USA (Received 22 November 2019; revised manuscript received 10 January 2020; accepted 22 January 2020; published 21 February 2020) Study of the adsorption of a transition metal on the surface of a layered material and the possible subsequent intercalation into that layered material is of fundamental interest and potential technological importance. In the present work, we choose the transition metal Cu as the adsorbate or intercalant and 2H -MoS2 as the layered material. Energetics are calculated characterizing four of the most basic surface and interfacial phenomena: adsorption, intercalation, diffusion, and adhesion. Using first-principles density functional theory (DFT), we find that intercalating a Cu atom into the van der Waals (vdW) gap below the MoS2 (0001) surface is 0.665 eV more favorable than adsorbing the Cu atom on top of the surface, i.e., intercalation of single Cu atoms is strongly favored thermodynamically. Also, we find that the system with adsorbed Cu is magnetic, while the system becomes nonmagnetic after the Cu atom is intercalated into the vdW gap. We obtain the diffusion barriers of the Cu atom on the surface and in the vdW gap to be 0.23 and 0.32 eV, respectively. We also obtain an adhesion energy of 0.874 J/m2 for a Cu (111) slab bonding with a 2H -MoS2 (0001) slab. The DFT value of adhesion energy, as well as the gap width at the interface between the Cu and MoS2, depends strongly on the choice of the functional. From our analysis on bulk properties of both MoS2 and Cu, we suggest that our vdW-DF2-B86R results listed above are reliable for applications, e.g., interpretation of experimental results and physical modeling of this adsorption system. DOI: 10.1103/PhysRevResearch.2.013182


I. INTRODUCTION
The interaction of metals with quantum materials, including layered materials or layered heterostructures, is important for applications. For example, metals and particularly Cu would naturally serve as interconnects or heat sinks in electronic devices based on these materials [1]. As another example, intercalated metals can serve to modify the properties of layered materials, such as transport [2] or magnetism [3].
Several basic energetic parameters are very useful for understanding or predicting the interaction between a metal and a layered material for a given combination of metal and substrate. One is the adsorption energy, which reflects the strength with which a single metal atom binds to the substrate. Similarly, the adhesion energy reflects the strength with which a metal film or slab binds to the substrate. Together with other quantities such as surface energy or cohesive energy, * y27h@ameslab.gov either of these quantities can be used to predict the growth mode or morphology of the adsorbate [4]. The adsorption energy also determines the residence time of the adsorbate on the substrate, hence determining the thermal stability of the adsorbate relative to the gas phase. The adhesion energy, together with surface energies, can also be used to predict the detailed equilibrium crystal shape of an adsorbate cluster supported on the substrate [5,6]. A third important parameter is the intercalation energy for individual atoms, which can be compared with the adsorption energy to predict whether there is a thermodynamic preference for the adsorbate to remain on top of the substrate or to penetrate the spaces or galleries between layers. A final parameter is the diffusion barrier for individual atoms adsorbed on top of a substrate surface or intercalated in the gallery. Such barriers are important since diffusion plays a role in both the nucleation dynamics [7,8] and coarsening kinetics [9] of metal islands.
In this paper, we analyze all four of these types of parameters: adsorption energies for various sites, adhesion energy, intercalation energies for various sites, and diffusion barriers directly using density functional theory (DFT). The results are used to make predictions for the interaction of Cu with a MoS 2 substrate which could be tested by suitably crafted deposition experiments.
The choice of Cu as a metal is partly motivated by the fact that Cu is commonly used in electronic devices, owing to its high conductivity and relatively low cost. MoS 2 is chosen in part because it can be considered a prototypical transition-metal dichalcogenide (TMD). TMDs have excellent structural, electronic, and optical properties, thus having multiple promising applications, in particular, as materials for semiconducting devices or catalysis [10][11][12][13]. TMDs have the formula MX 2 , where M represents a transition metal (e.g., Mo or W) in group IV, V, or VI, and X denotes a chalcogen atom (e.g., S, Se, or Te). A bulk TMD crystal is composed of an array of parallel X -M-X triple layers (TLs) (often called "monolayers" in the literature) where adjacent TLs are held together via weak van der Waals (vdW) attractions. Thus, the space between two adjacent TLs has been commonly called the "vdW gap", whereas the bonds between X and M atoms within each X-M-X TL are covalently ionic in character.
The MoS 2 crystal can form several different phases [10][11][12], e.g., 2H-MoS 2 and 3R-MoS 2 , where the "H" and the "R" indicate hexagonal and rhombohedral symmetry, respectively, while "2" or "3" indicates that the corresponding unit cell of a MoS 2 phase includes 2 or 3 TLs. In this study, we focus on the phase that is thermodynamically stable at room temperature, 2H-MoS 2 . The 2H-MoS 2 unit cell (space group P6 3 /mmc) contains two Mo atoms and four S atoms, as illustrated in Fig. 1(a). The experimental lattice parameters [14] are a = 0.315 nm and c = 1.230 nm; the distance between two S monatomic layers (MLs) within one TL is d S = 0.317 nm, as indicated in Fig. 1(a). Then the thickness of the vdW gap is d vg = c/2 − d S = 0.298 nm for bulk 2H-MoS 2 . Figure  1(b) shows the vdW gaps for a 3-TL 2H-MoS 2 (0001) slab. Below, we sometimes omit the prefix "2H-" in the notation of 2H-MoS 2 for conciseness.
There exist no previous DFT analyses for the Cu-plusbulk-MoS 2 system. While a few DFT results are available for Cu interacting with one or two TLs of MoS 2 , their collective utility is limited because of issues related to accuracy and other factors discussed below. For example, Huang et al. [15] have performed DFT calculations for the adsorption energies of a Cu atom on a 1-TL MoS 2 slab (space group P6m2), but there have not yet been calculations for a Cu atom adsorbed on a bulklike MoS 2 (0001) slab. As another example, Guzman et al. [16] calculated the intercalation energies of a Cu atom at two high-symmetry sites in the vdW gap between 2 TLs, but they did not calculate the adsoption energies on top of the MoS 2 slab. Hence, their results cannot be used to predict whether intercalation is preferred. As a third example, in the above DFT calculations, Huang et al. [15] used the optPBE-vdW functional [17], whereas Guzman et al. [16] used DFT-D3 [18]. One does not know whether results from the two groups can or should be compared, since the choice of functional is particularly critical in systems that combine metals with vdW (layered) materials. For a pure metal, it is often true that the properties predicted from a functional with a dispersion correction are significantly worse relative to the results from a functional without the dispersion correction, as will be discussed for Cu in Sec. II. Therefore, a critical consideration for a transition-metal-plus-vdW-substrate system is that the use of a functional with a dispersion correction good for the substrate material (MoS 2 ) cannot significantly worsen the predicted properties for the transition metal (Cu).
The present work has three distinguishing features: (i) All quantities are calculated self-consistently so that comparisons can be made, e.g., between adsorption and intercalation energies, and thus the energies can be used collectively for modeling. (ii) Careful attention is paid to choosing a functional that is optimized for both the metal and the layered substrate (a feature which especially affects the adhesion energy). (iii) The MoS 2 substrate is modeled as multiple TLs, which makes our results applicable to systems with multilayer or bulk MoS 2 .
This paper is organized as follows. Section II describes the DFT method and benchmark analyses indicating the most appropriate choice of functional from among six different types of functionals. Section III provides analysis of the energetics and magnetism of a Cu atom at various sites on top of the MoS 2 (0001) slab and also in the vdW gap. The diffusion barrier in each case is also calculated. Section IV presents the adhesion energy between a Cu (111) slab and a MoS 2 (0001) slab. Section V provides a summary.

II. DFT METHOD AND BENCHMARKING
We use the Vienna Ab Initio Simulation Package (VASP) [19] to perform first-principles DFT calculations with the projector-augmented-wave pseudopotentials [20] generated and released in 2013 by the VASP group for the electroncore interactions. Considering the significant contributions from vdW interactions for a dispersion-bonded system such as the layered material MoS 2 , we test six different types of functional with dispersion corrections for the exchangecorrelation energy: vdW-DF [21], DFT-D3 [18], optPBE-vdW [17], optB88-vdW [17], vdW-DF2-B86R [22], and SCAN-rVV10 [23], as described below. We choose the -centered k mesh according to the supercell size, as listed below for each specific system. Converged total energies are reached when the force on each relaxed atom is less than 0.1 eV/nm. In calculations for a surface system modeled as a periodic slab, the vacuum thickness between two adjacent slabs is not less than 2.0 nm. Dipole corrections and spin polarization are considered in all surface calculations.  Fig. 1). B 0 and B 0 are the bulk modulus and the first derivative of the bulk modulus with respect to pressure, respectively. The cohesive energy is E coh = E Cu − σ bulk , where E Cu is the energy of a single Cu atom in the gas phase, and σ bulk is the energy per primitive cell for fcc Cu.

System
Method For a benchmark analysis, we calculate the bulk properties of 2H-MoS 2 and fcc Cu. For MoS 2 , we use a unit cell containing two Mo atoms and four S atoms, as illustrated in Fig. 1(a); the k mesh is taken to be 51 × 51 × 13 and the cutoff energy is 550 eV. For fcc Cu, we use the corresponding primitive cell containing just one Cu atom; the k mesh is taken to be 51 × 51 × 51 and the cutoff energy is 600 eV. Table I lists the bulk properties of 2H-MoS 2 and fcc Cu from six different functionals. To compare with experimental values, we calculate the relative error δ = (Q DFT − Q exp )/Q exp × 100%, where Q DFT (Q exp ) represents a quantity such as a a lattice parameter or a cohesive energy, obtained from DFT calculation (measured from experiment). In Fig. 2, we show the δ values for five experimentally available quantities including the lattice parameters a, c, and d S for 2H-MoS 2 , lattice constant a for fcc Cu, and cohesive energy E coh for fcc Cu. As indicated by the dashed green curve with diamonds in Fig. 2, E coh (Cu) predicted from vdW-DF, DFT-D3, or SCAN-rVV10 is 15.5% smaller, 13.5% larger, and 15.6% larger, respectively, than the experimental value. Thus, these three functionals are not good choices for calculating a system with Cu, even though DFT-D3 or SCAN-rVV10 can indeed predict good lattice constants for both MoS 2 and Cu (see Fig. 2). The vdW-DF functional also predicts the worst lattice parameters, i.e., the δ values of the lattice parameter c for MoS 2 , lattice constant a for Cu, and lattice constant a for MoS 2 are 6.6%, 2.8%, and 2.7%, respectively. The optPBE-vdW functional is a significantly better choice than vdW-DF, DFT-D3, or SCAN-rVV10, with the biggest errors being δ = −3.4% and 3.4% for E coh (Cu) and c (MoS 2 ), respectively. Overall, as in Fig. 2, both optB88-vdW and vdW-DF2-B86R can reproduce all five experimental values well. The value of δ = 1.7% for E coh (Cu) from optB88-vdW is smaller than δ = 3.2% from vdW-DF2-B86R, but the relative errors for a (Cu), a (MoS 2 ), and c (MoS 2 ) from vdW-DF2-B86R are even smaller than those from optB88-vdW (Fig. 2). Thus, both optB88-vdW and vdW-DF2-B86R provide excellent choices if one only considers the five bulk properties (Fig. 2) which are experimentally available. In this work, we choose the vdW-DF2-B86R functional for the Cu-plus-MoS 2 systems.
In the above calculations for the lattice parameters of bulk MoS 2 , we fully relax all atoms in the unit cell and allow cell volume and shape to change. To confirm that the final  Table I. structure is at the global energy minimum E 0 with the volume V 0 , we calculate the energies of structures with a series of volumes fixed around the V 0 value by allowing cell shape to change. Then, we fit the Birch-Murnaghan equation [26,27] to the energy-volume data points. As an aside, the fitting curve also yields the bulk modulus B 0 of 2H-MoS 2 and the first derivative B 0 of the bulk modulus with respect to pressure, which are listed in Table I. No experimental values are available for B 0 and B 0 of bulk 2H-MoS 2 .
For the above calculations, all six functionals include dispersion corrections. A brief comparison with previous DFT calculations from the Perdew-Burke-Ernzerhof (PBE) functional without dispersion correction is worthwhile. Using the PBE functional, Li et al. [28] calculated the bulk properties of 2H-MoS 2 and obtained the lattice constants a = 0.3199 nm with δ = 1.56% and c = 1.2493 nm with δ = 1.57%, which is better than the values from DFT-DF and optPBE-vdW, but significantly worse than DFT-D3, optB88-vdW, and vdW-DF2-B86R. They also obtained a bulk modulus of 63.361 GPa [28], which is significantly larger than B 0 values that we obtain from all six functionals with dispersion corrections in Table I. For fcc Cu, the cohesive energy E coh is 3.486 eV with δ = −0.97% (the lattice constant a is 0.3635 nm with δ = 0.91%), obtained from our recent PBE calculations [29], while the δ values for E coh (for a) of fcc Cu from vdW-DF, DFT-D3, and optPBE-vdW are −15.46% (2.85%), 13.50% (−0.96%), and −3.41% (1.35%) (see Fig. 2), respectively. This indicates that a functional with the dispersion correction can sometimes give significantly worse predictions of the bulk properties than PBE, e.g., for the pure metal, Cu, as we mentioned in Sec. I. Thus, choosing an optimal functional for calculating a A-plus-B system needs to balance the properties of both A and B.
To summarize this section, we have carefully assessed six different functionals in terms of their ability to correctly predict five bulk parameters for both Cu and MoS 2 . Three are clearly better than the others, and from those we choose vdW-DF2-B86R for subsequent work in this paper.

III. ADSORPTION, DIFFUSION, AND INTERCALATION OF A Cu ATOM AT THE MoS 2 (0001) SURFACE
One key goal of our analysis is to assess the relative stability of a single Cu atom adsorbed on the surface of 2H-MoS 2 (0001) versus intercalated underneath the surface. To this end, we calculate the Cu chemical potential defined as [30,31] where E tot is the total energy of the Cu-plus-MoS 2 system, E MoS 2 is the energy of the MoS 2 substrate, and E Cu is the energy of one Cu atom in the gas phase. This quantity, μ Cu , reflects the stability of a Cu atom in any configuration, including adsorbed and intercalated. For adsorption, Eq. (1) is equivalent to the conventional adsorption energy, i.e., E ads ≡ μ Cu . In order to facilitate comparison with results in the literature, we also introduce an intercalation energy [16,15] where σ bulk is the energy per primitive cell (including one Cu atom) for bulk fcc Cu, instead of E Cu in Eq. (1). Then, E int quantifies the energy required for a Cu atom to become intercalated, starting from bulk fcc Cu. To simulate adsorption or intercalation for various Cu-plus-bulk-MoS 2 systems, we use a MoS 2 slab with thickness of 3 TLs [see Figs. 1(b) and 3(c)], and the bottommost TL is always fixed during relaxation. The lateral size of the supercell is taken to be 5 × 5 in units of a. The k mesh is taken to be 3 × 3 × 1.
To find the lowest μ Cu (or E ads ) of a Cu adatom on the MoS 2 slab surface, we initially place the Cu atom on five sites: fcc, hcp, top, bridge, and center, as shown in Fig. 3(a). The most favorable adsorption site is the hcp site with μ Cu = −1.485 eV. The Cu atom is 0.126 nm higher than the nearest three S atoms, and the Cu-S bond length is 0.226 nm. At the fcc site and top site, μ Cu = −1.342 eV and μ Cu = −0.970 eV, respectively. Adsorption of a Cu atom at the bridge or center site is unstable, since the adatom moves to a local hcp site after full relaxation. The difference μ Cu between μ Cu at the two most stable sites (fcc and hcp sites) is (−1.342 eV) − (−1.485 eV) = 0.143 eV.
In order to avoid excessively expensive computations [31] in estimating the diffusion barrier using a climbing nudgedelastic-band (CNEB) method [32], we consider a 1-TL MoS 2 (0001) slab with the bottommost S ML fixed during the relaxation. We first obtain the chemical potentials at hcp and fcc sites using the 1-TL slab and find almost-unchanged values (the errors are less than 0.5 meV) relative to those from the above 3-TL-slab calculations. This indicates that the effect of slab thickness on these adsorption energies can be neglected. We then use seven images plus the configurations with the Cu atom at hcp and fcc sites as two endpoints and obtain a diffusion barrier of 0.23 eV from the corresponding minimum energy path (MEP), as depicted in Fig. 4(a).
Using a 1 × 1 rectangular cell and the optPBE-vdW functional, Huang et al. [15] calculated the adsorption energies of a Cu atom on a 1-TL MoS 2 slab. Their results are very similar to ours. They identified the hcp site as the most stable, found μ Cu ≈ −1.36 eV at the hcp site (vs μ cu = −1.485 eV in the present work), and reported differences in adsorption energies at various sites comparable to ours. Notably, they reported that μ Cu ≈ 0.13 eV, close to our value of 0.143 eV. We also note that μ cu = −1.95 eV at the hcp site, which is the most favorable adsorption position from Li et al.'s PBE-D2 calculations for a 1-TL MoS 2 slab with a 4 × 4 supercell size [34]. The possible contributions to the discrepancy between their and our values of μ cu include differences in functional, cell shape and size, cutoff energy, and thicknesses of the MoS 2 slab, but our tests show that the main contribution comes from the use of different functionals.
A rough assessment of the growth mode of the metal film on MoS 2 comes from knowledge of the adsorption energy together with cohesive energy. The guideline is that threedimensional (3D) growth is possible when the strength of the adsorption energy of the metal on the substrate is less than the cohesive energy of the metal [35,36]. In the present case, the magnitude of adsorption energy of Cu is 1.49 eV (Table II) and the cohesive energy is 3.63 eV (Table I), so the criterion for 3D growth is met easily in this system. We emphasize that a more sophisticated assessment of the growth mode is discussed in Sec. IV.
To assess the propensity for intercalation of a Cu atom in the gallery or vdW gap underneath the top TL, we also calculate μ Cu for five sites in the gap: fcc-i, hcp-i, top-i, bridge-i, and center-i. These locations are assigned based on the position of the intercalated Cu atom relative to Mo and S atoms in the top TL. In each case we initially place the Cu atom underneath the top TL and then relax the structure. Values of μ Cu and E int , obtained after full relaxation, are listed in Table II Table II. The atom at the bridge-i or center-i site is unstable and moves to a nearby local energy minimum site, which is found to be a fcc-i site, after full relaxation. Therefore, for a Cu atom, intercalation into the vdW gap beneath the top TL is always more favorable than adsorption on the top TL. In addition, there is a small difference of only 0.020 eV between the μ Cu values at top-i and hcp-i sites.
To efficiently compute the diffusion barrier in the vdW gap, we use a 2-TL MoS 2 slab with an outermost S ML fixed during the relaxation to perform our CNEB calculations. We first obtain the chemical potentials of −2.035, −2.070, and −2.093 eV at fcc-i, hcp-i, and top-i sites using the 2-TL slab, and the differences are 0.019, 0.061, and 0.057 eV relative to those from the above 3-TL-slab calculations, respectively, indicating that the slab thickness effects on these chemical potentials are larger than on adsorption energies but not very significant (the errors are below about 60 meV). Our CNEB calculations show that the MEP in the vdW gap is from a top-i site to fcc-i site and then to a hcp-i site with two almost-equalheight saddle points, the highest of which is 0.32 eV, i.e., the diffusion barrier, as depicted in Fig. 4(b).
Using a 6 × 6 rectangular cell and DFT-D3 functional, Guzman et al. [16] calculated E int for a Cu atom at top-i and fcc-i sites in the vdW gap between 2 TLs and obtained values of around 1.2 eV with about 0.12 eV lower energy at the top-i site than at the fcc-i site. This result is in reasonable agreement TABLE II. DFT results for chemical potential μ Cu and magnetic moment m (in units of Bohr magneton μ B per cell) for adsorption and intercalation of a Cu atom after full relaxation at different initial sites. Sites fcc, hcp, top, bridge, and center are defined in Fig. 3(a) 4. (a) The MEP of a Cu atom absorbed on top of a 1-TL MoS 2 slab, diffusing from a hcp site to its nearest-neighbor fcc site. (c) The MEP of a Cu atom within the vdW gap of a 2-TL MoS 2 slab after intercalation, diffusing from a top-i site to its nearest-neighbor fcc-i site, and then to its nearest-neighbor hcp-i site. In (a) or (b), the inset shows the top-view trajectories of atoms for the corresponding MEP. The curves are generated from a modified Bézier method [33] by fitting the data points from our CNEB calculations versus reaction coordinates. Green dots on a curve represent CNEB images. with our result (in Table II), E int = 1.483 eV at the top-i site, which is is 0.097 eV lower than at the fcc-i site. The differences between their and our values of E int are, again, related to the uses of different DFT methods (including functional, cell size, cutoff energy, etc.) and different thicknesses of the MoS 2 slab, but mainly from the choices of different functionals.
We also analyze the bond lengths between the Cu atom and its nearby S atoms for the three (meta)stable intercalation sites: (i) Cu at top-i. As discussed above, the energy of the Cu atom at this site is lowest. From the symmetry of the system, this site is a top site relative to its upper TL, while it is a hcp site relative to its lower TL, as shown in Figs. 3(b) and 3(c). The vertical spacing (i.e., the height difference) from upper to lower S atoms is 0.326 nm. The Cu atom is 0.212 nm below the upper S atom (i.e., the Cu-S bond length is 0.212 nm) and is 0.114 nm higher than the lower three S atoms (for which the Cu-S bond length is 0.224 nm). Therefore, the intercalated Cu atom is located about 1/3 of the way up the vdW gap from the lower TL, i.e., in the center of a near-tetrahedron defined by four nearest S atoms.
(ii) Cu at hcp-i. Just considering the top and middle TL, the symmetry of this subsystem implies that the geometry of the hcp-i site is just the inversion of that for the top-i site. We indeed find that the bond lengths from our DFT calculations for the hcp-i site are correspondingly inverted relative to those for the top-i site. This indicates that the effect of the MoS 2 slab thickness (i.e., the number of TLs) on the geometry of intercalated Cu sites is negligible.
(iii) Cu at fcc-i. From the symmetry of the system, the fcci site is located at the hollow sites of both upper and lower TL. The vertical distance from upper to lower S atoms is 0.315 nm, and the Cu atom is 0.158 nm lower (0.157 nm higher) than the upper (lower) S atoms. Therefore, in this site the Cu atom is almost at the midpoint of the vdW gap. The Cu-S bond lengths with upper and lower S atoms are 0.241 nm.
Finally, we note that Cu-adsorbed MoS 2 systems are magnetic, but the magnetic moments are significantly reduced after intercalation (see Table II), while a pure MoS 2 slab is nonmagnetic due to the strong covalent bonding between Mo and S atoms. To understand the origin of magnetic properties for adsorption and intercalation, in Fig. 5 we plot the projected density of states (PDOS) for the most stable configurations for a Cu atom adsorbed at a hcp site on a 3-TL MoS 2 slab surface and intercalated at a top-i site of the 3-TL MoS 2 slab, as shown in Fig. 3. For Cu adatom adsorption [ Fig. 5(a)], near the Fermi level there are significant overlaps of 4s states of the Cu atom, 3s and 3p states of S atoms, and 4d states of Mo atoms, i.e., the unpaired 4s electron of the Cu atom induces the strong orbital hybridization corresponding to these overlapped states and also the notable separation between spin-up and spindown states of Cu, S, or Mo, therefore causing the relatively large magnetic moment. After intercalation at the top-i site, the magnetic moment decreases to almost zero, indicating that the intercalation can equate the spin-up and spin-down states, as revealed in Fig. 5(b).
To summarize the results in this section, the preferred adsorption site of a Cu atom on the top surface of MoS 2 is the hcp site, where it is at the center of three S atoms and directly atop a Mo atom. Note that the magnitude of μ Cu = −1.485 eV is well below that of the bulk cohesive energy E coh = 3.634 eV for Cu from our vdW-DF2-B86R calculation (Table I). Thus, the system energy is lowered by incorporating isolated adsorbed Cu atoms into a 3D island, suggesting a 3D growth of Cu on MoS 2 (in the absence of intercalation). A more rigorous analysis of growth mode is presented in Sec. IV. The difference between the μ Cu values for hcp and fcc sites is 0.143 eV, and the diffusion barrier of the Cu atom on top of the surface is estimated as 0.23 eV from our CNEB calculations.
Regarding intercalated Cu atoms, the preferred adsorption site is the top-i site, where it sits at the center of a neartetrahedron defined by four nearest S atoms and has μ Cu = −2.151 eV. The difference between the μ Cu values for this site and the hcp-i site is very small (about 0.02 eV) and is reasonable based on symmetry arguments. The diffusion barrier of the Cu atom in the vdW gap is estimated as 0.32 eV by using a 2-TL MoS 2 slab. Regardless of site, the value of μ Cu is much more negative for the intercalated Cu atom than for the adsorbed Cu atom on the top. The energy at the most favorable site in the vdW gap is μ Cu = 0.665 eV lower than at the most favorable adsorption site on top of the surface. Thus, there is a strong thermodynamic preference for intercalation of single atoms. Regarding kinetics, intercalation can only occur if there is a facile pathway for accessing the gallery or vdW gap, e.g., through surface defects which act as portals. Then, the above single-atom energetics suggests that deposited atoms may preferentially diffuse to the gallery, which could lead to 3D island formation in the gallery rather than on top of the surface (provided the strain penalty for such island growth is not too large). This observation motivates experiments involving Cu deposition on MoS 2 which has been sputtered prior to deposition with the aim of creating portal defects facilitating intercalation. This scenario for intercalation is realized for the Cu-plus-graphite system where μ Cu = 0.500 eV is smaller than the current system [30,31,37]. As a final comment, we also find that an adsorbed Cu on MoS 2 system is magnetic, while the system can become nonmagnetic after intercalation.

IV. ADHESION ENERGY BETWEEN Cu (111) SLAB AND MoS 2 (0001) SLAB
The adhesion energy between two slabs bonded at an interface is calculated as [29] where E s1 , E s2 , and E s1s2 are the energies of the slab s1, slab s2, and the s1-s2 system, respectively; A is the s1-s2 interface area.
To minimize possible artificial contributions from strain effects due to lattice mismatch between a Cu (111) slab and MoS 2 (0001) slab, we use a (5 × 5)Cu/(4 × 4)MoS 2 unit cell, noting also that there exists experimental evidence for this interface configuration [38]. See Fig. 6 for a 5 × 5 Cu (111) ML on a 4 × 4 MoS 2 (0001) substrate. In this work, we always use the unstrained MoS 2 (0001) slab, i.e., with a lateral supercell size of 4 × 4 (in units of a MoS 2 , which is taken to be the vdW-DF2-B86R value of 0.3168 nm, as listed in Table I). The lattice of the Cu (111) slab with a lateral supercell size of 5 × 5 matching the unstrained 4 × 4 MoS 2 (0001) slab has a lateral lattice constant of a str = 4a MoS 2 /5. Calculations of adhesion energies use a 3-TL MoS 2 (0001) slab plus a 5-ML Cu (111) slab. During energy minimization, atoms in the bottommost TL are fixed at their bulk lattice points and other atoms relax, as indicated in Fig. 7(a). In the determination of β from Eq. (3), the lateral lattice constant of the separated Cu slab is retained at the value a str .  For a Cu (111) slab with a thickness of 5 ML on the unstrained 4 × 4 MoS 2 (0001) slab, there is a tensile lateral strain (a str − a 5 )/a 5 × 100% = 0.4%, where a 5 = 0.2523 nm is the lateral equilibrium lattice constant of the 5-ML Cu (111) slab from our vdW-DF2-B86R calculations. Note that for a bulk Cu (111) slab (i.e., with infinite thickness), the lateral equilibrium lattice constant is a ∞ = a Cu / √ 2 = 0.2550 nm, which is larger than a 5 . Here a Cu = 0.3606 nm is our vdW-DF2-B86R value in Table I. Then there is a compressive lateral strain (a str − a ∞ )/a ∞ × 100% = −0.6% relative to a ∞ . If one uses the energy of the Cu slab with lattice constant a 5 (a ∞ ) in Eq. (3), the adhesion energy is denoted as β * (β * * ) [29]. Table III lists the β, β * , and β * * values of the 5ML_Cu/3TL_MoS 2 system [ Fig. 7(a)]. The values corresponding to two k meshes of 3 × 3 × 1 and 7 × 7 × 1 in Table III indicate the good convergence of the adhesion energies.
In Table III, we also list the values of adhesion energies for s1 and s2 slabs, which are "frozen" or "relaxed" when their energies are calculated separately by cutting the fully relaxed geometry of their corresponding combined s1s2 system into s1 and s2 slabs.
The appropriate definition of adhesion (β, β * , or β * * ) depends on the specific application. As an example, for the construction of equilibrium crystal shape of epitaxial Cu nanoclusters on MoS 2 , we would suggest using the adhesion energy β = 0.874 J/m 2 , given that the Cu cluster is strained to match the bulk MoS 2 substrate. Let us return to the issue discussed in Sec. III for the growth mode of Cu films on MoS 2 . The precise criterion for 3D Cu growth is that the adhesion energy is less than twice the Cu (111) surface energy γ Cu(111) , where γ Cu(111) ≈ 1.3 J/m 2 [29]. Thus, again it is clear that the growth mode for Cu films on MoS 2 is 3D. TABLE III. Cohesion energies for the 5ML_Cu/3TL_MoS 2 and 5ML_Cu/1TL_MoS 2 systems from our DFT calculations using the vdW-DF2-B86R functional compared with the vdW-DF results for the 5ML_Cu/1TL_MoS 2 system. The structures are shown in Fig. 7. Cases where the s1 and s2 slabs are frozen or relaxed, as well as the k mesh used in each calculation, are noted. The unit of eV/cell for β is eV per unit cell of 2H -MoS 2 (0001) surface. For the definitions of β, β * , and β * * , see text. Previously, Le et al. [38] calculated the adhesion energy (which they called negative binding energy) of a 1-TL MoS 2 slab on 5-ML Cu (111) using the vdW-DF functional [21] and obtained a value of 0.27 eV/cell, which is much smaller than our value β = 0.474 eV/cell (= 0.874 J/m 2 ) for the 5ML_Cu/3TL_MoS 2 system. To determine the reason for this significant discrepancy, we first calculate β * and β * * values, considering that (i) a different lateral strain in the Cu slab due to the use of a different lattice constant of Cu slab can result in some error, and (ii) information about the lattice constant of the Cu slab in Ref. [38] is not available. However, as listed in Table III, both the β * = 0.461 eV/cell and β * * = 0.478 eV/cell from our vdW-DF2-B86R calculations do not significantly differ from the β = 0.474 eV/cell. Therefore, the significant discrepancy is not from the definition of adhesion energy.
Next we use a 5ML_Cu/1TL_MoS 2 system, matching the analysis of Le et al. [38], and obtain a β value of 0.463 eV/cell (see Table III), which also does not significantly differ from β = 0.474 eV/cell for the 5ML_Cu/3TL_MoS 2 system. This result indicates that the thickness (in TLs) of the MoS 2 (0001) slab has little effect on the β value. For this calculation, the relaxed structure is shown in Fig. 7(b).
Finally, we calculate the same 5ML_Cu/1TL_MoS 2 system as above but use the vdW-DF functional instead of vdW-DF2-B86R, and the relaxed structure is shown in Fig. 7(c). This result reveals that β is reduced to a much smaller value of 0.200 eV/cell, which is close to the value of 0.24 eV/cell that Le et al. [38] obtained using vdW-DF. The difference between these two values, 0.200 and 0.24 eV/cell, might be from the difference in other input parameters, e.g., for convenience, we still use a lateral lattice constant a MoS 2 = 0.3168 nm from vdW-DF2-B86R. In fact, the lateral lattice constant value of 0.3168 nm that we use for the 5ML_Cu/1TL_MoS 2 system is close to 0.3255 nm from Le et al.'s vdW-DF calculation [38]. Thus the resulting small lateral strain should not markedly influence the β value [29]. Consequently, the large discrepancy in the β value is from the use of the different functionals, vdW-DF2-B86R versus vdW-DF. The latter is the old version of the former and predicts much worse bulk properties for both MoS 2 and Cu, as discussed in Sec. II.
Along with the large discrepancy in the β value from the two different functionals, vdW-DF predicts larger interlayer spacings than vdW-DF2-B86R, especially for the interface spacing between Cu and MoS 2 slabs: 0.280 nm (vdW-DF) vs 0.237 nm (vdW-DF2-B86R), as indicated in Figs. 7(b) and 7(c). The significant increase in the interface spacing from Fig. 7(b) to Fig. 7(c) corresponds to weaker adhesion between the two slabs. Notably, the corrugation of each ML from the two functionals is similar (Fig. 7), and therefore the corrugation does not play a role in the different β values. In principle, this provides a way to experimentally test which of the two functionals is better: A measurement of the spacing at the Cu-MoS 2 interface, e.g., via transmission electron microscopy, of 0.28 nm would favor vdW-DF, whereas 0.24 nm would favor vdW-DF2-B86R.
To summarize this section, values of the adhesion energy fall between 0.46 and 0.48 eV/cell, depending upon the lattice constant used for the Cu slab. This variation in adhesion energy is small compared with the difference by a factor of about 2 that can arise from using different functionals. Smaller adhesion energy corresponds, physically, to a larger gap at the interface between the Cu slab and the MoS 2 slab. With the functional that we employ here, a gap width of 0.237 nm is predicted.

V. CONCLUSIONS
We have performed first-principles DFT calculations for adsorption, intercalation, and adhesion energies of Cu at a 2H-MoS 2 slab surface using vdW-DF2-B86R, which is a reliable functional in terms of predicting bulk properties of both 2H-MoS 2 and fcc Cu. We find that the intercalation of a Cu atom underneath the top TL, i.e., into the first vdW gap, is 0.665 eV more favorable than adsorption of the Cu atom on top of the slab surface. On top of the slab surface, the most favorable position of a Cu atom is at a hcp site with chemical potential of −1.485 eV, while in the first vdW gap, the most favorable position is at a top-i site with a chemical potential of −2.151 eV. DFT results from our CNEB calculations show a diffusion barrier of 0.23 eV for the Cu atom on top of the surface and 0.32 eV within the vdW gap. We also find that the magnetic moment of 1 μ B for adsorption of a Cu atom on the MoS 2 (0001) surface can reduce to almost zero after intercalating into the vdW gap and therefore provides useful information for potential technological applications in, e.g., spintronics.
We also obtain a value of 0.874 J/m 2 for the adhesion energy of a Cu (111) slab bonding with a 2H-MoS 2 (0001) slab. The value of 0.874 J/m 2 (0.474 eV/cell) is about twice the value (0.24 eV/cell) from Le et al.'s vdW-DF calculation [38]. Our analysis shows that this discrepancy is due to the different functionals. Along with the difference in adhesion energies, there is a corresponding difference (0.043 nm) in the interlayer spacing between Cu and MoS 2 slabs. The DFT energy parameters obtained from our vdW-DF2-B86R calculations in this work can be directly applied for analyzing the experimental results for Cu deposition on 2H-MoS 2 substrates as well as the energetic parameters for modeling such systems. In particular, comparison of our result for adhesion energy with the Cu (111) surface energy implies a Volmer-Weber growth mode for Cu on 2H-MoS 2 (0001).