First principles investigation of copper and silver intercalated molybdenum disulfide

We characterize the energetics and atomic structures involved in the intercalation of copper and silver into the van der Waals gap of molybdenum disulfide as well as the resulting ionic and electronic transport properties using first-principles density functional theory. The intercalation energy of systems with formula (Cu,Ag)xMoS2 decreases with ion concentration and ranges from 1.2 to 0.8 eV for Cu; Ag exhibits a stronger concentration dependence from 2.2 eV for x = 0.014 to 0.75 eV for x = 1 (using the fcc metal as a reference). Partial atomic charge analysis indicates that approximately half an electron is transferred per metallic ion in the case of Cu at low concentrations and the ionicity decreases only slightly with concentration. In contrast, while Ag is only slightly less ionic than Cu for low concentrations, charge transfer reduces significantly to approximately 0.1 e for x = 1. This difference in ionicity between Cu and Ag correlates with their intercalation energies. Importantly, the predicted...


I. INTRODUCTION
Transition metal dichalcogenides (TMDCs) are layered materials with formula MX 2 , where M is a group IV, V, or VI transition metal and X is a chalcogen atom. The TMDC layers are held together via van der Walls (vdW) interactions and the X-M-X intra-layer bonds have covalent-ionic character. The ability to exfoliate and synthesize single and few layer TMDCs and the range of properties achievable make them attractive for a wide range of applications. Quantum confinement and significant flexibility in terms of chemistry lead to a range of electronic properties, in contrast to the gap-less graphene. In addition to chemistry and size, their properties can be modulated via strain, doping, alloying, and combining dissimilar materials to create heterostructures. [1][2][3][4] These properties and progress in fabrication enabled the demonstration of several devices, including transistors, photovoltaic cells that require thin and transparent semiconductors, and flexible electronics. [5][6][7] Recently, the semiconductor-to-semimetal phase transition in two-dimensional materials via gate voltage has been investigated by ab initio first principles methods 8 suggesting potential applications in non-volatile memory devices 9 and other bi-stable electronics. 10 In addition to these applications, the presence of a vdW gap between layers makes TMDCs attractive as potential hosts for intercalation with electron donor or acceptor species. Recently, intercalation of group VI TMDCs with alkali metals has been investigated for nanoelectronics, energy storage, and catalysis applications. Lithium intercalation of MoS 2 , for example, has been reported to improve the optical transmission and electrical conductivity due to changes in the electronic structure and large injection of free carriers, 11 and potassium intercalated MoS 2 has been investigated as a promising catalyst for the synthesis of alcohols. 12 Reversible intercalation/deintercalation of Li, Na, and Mg in MoS 2 opens the potential for battery [13][14][15][16] and pseudocapacitor 17 applications.
From a basic science point of view, metal intercalation enabled the exploration of interesting material properties; for example, 3d transition metal species such as V, Cr, Mn, Fe, Co, and Ni intercalated into group V TMDCs form ordered superlattices that exhibit ferromagnetic and antiferromagnetic orderings depending on the concentration and temperature. 18 Cu, Cs, and Li intercalated TaS 2 and TaSe 2 exhibit a charge density wave phase at room temperature, 19 and Cu intercalated TiSe 2 has been recently explored for low temperature thermoelectric applications. 20  some TMDCs exhibit semiconductor-to-metal transition 21 and even the development of a superconducting state for copper into TiSe 2 and potassium into MoS 2 . [22][23][24] Despite significant interest in the intercalation of group VI TMDCs with metallic ions, several questions remain unanswered. For example, the energetics of the intercalation process and the resulting atomic structures as a function of intercalate concentration are not known, nor are the ionic mobility and the evolution of the electronic structure during metallic loading. In this paper, we report on density functional theory (DFT) calculations to answer these questions with the objective of helping assess the possible use of intercalated TMDCs in applications of interest in nanoelectronics and energy storage.
The calculations in this paper demonstrate that electrochemical intercalation is possible up to very high metal concentrations with an accompanying reduction in electrical resistance. The high-mobility predicted for the metallic ions within the vdW regions between the TMDC layers indicates that these materials are mixed ionic-electronic conductors and that copper ions could be made to migrate through the vdW gaps under relatively low external bias.
The remainder of the paper is organized as follows. Section II describes the metal intercalated MoS 2 models and details the energetics, diffusion, and transport simulation methods; in Section III, we discuss the intercalation energy as a function of metal concentration as well as the activation energy for diffusion of metal ions in MoS 2 . Section IV A presents the electronic band structure and density of states of the intercalated systems. The electronic transport properties are discussed in Section IV B and conclusions are provided in Sec. V.

II. SIMULATION DETAILS
A. Density functional theory calculations All structural optimizations and total energy calculations were carried out using DFT as implemented in the Vienna ab initio simulation package (VASP). 25,26 For the electron-ioncore interactions, we use the projector-augmented wave (PAW) description, 27,28 and the electron exchangecorrelation potential is calculated within the generalized gradient approximation as proposed by Perdew, Burke, and Ernzerhof (GGA-PBE). 29 It is well known that London dispersion forces, originating from induced-dipole induceddipole interactions, are not captured by DFT based on the Born-Oppenheimer approximation. To account for such interactions, we use a non-local van der Waals correction to the total energy (DFT-D3) as proposed by Grimme et al. 30,31 The D3 flavor of dispersion-corrected DFT has been used in a wide variety of materials to predict their energetic, thermodynamic, and kinetic properties, 32,33 as well as intercalation chemistry in layered materials. [34][35][36] In DFT-D3, the total energy of the system is obtained as where E DFT is the Kohn-Sham energy obtained selfconsistently using GGA-PBE as the exchange-correlation functional, 37 and E non-local is the dispersion correction, which contains the sum of a two-and a three-body term. 30 The D3 method is more sophisticated than its predecessor D2 38 in the sense that the dispersion coefficients C AB n are dependent on the local environment of atoms A and B.

B. Metal intercalated MoS 2 models
We start with the hexagonal trigonal prismatic 2H-MoS 2 unit cell (space group P6 3 /mmc) with experimental lattice parameters a ¼ b ¼ 3.16 Å and c ¼ 12.294 Å . [39][40][41] A metastable phase of MoS 2 is known to exist with an octahedral coordination (1 T); 42 however, this phase shows a metallic character and its stability is contingent upon certain environments that require electron donation to the TMD layers. 43,44 Furthermore, we predicted the 1 T phase to be about 0.8 eV per formula unit higher in energy than the 2 H phase. This work focuses only on the stable 2H phase of MoS 2 .
The lattice parameters and atomic positions of the 2H-MoS 2 model are fully optimized using the DFT flavor mentioned above. The Brillouin zone integrals are approximated using a Monkhorst-Pack 45 k-mesh of 15 Â 15 Â 3 points for the structural optimization and 21 Â 21 Â 3 for the total energies and calculation of the electronic properties. The kinetic energy cutoff for the plane-wave basis is set to 500 eV, the energy convergence tolerance for the self-consistent field calculation is taken as <1 Â 10 À5 eV, and the force convergence tolerance for the conjugate gradient ionic positions and lattice optimization is set to <0.01 eV/Å .
The predicted lattice parameters for the pristine 2H-MoS 2 bulk crystal are a ¼ b ¼ 3.160 Å and c ¼ 12.334 Å , and a ¼ b ¼ 90 and c ¼ 120 . Our DFT-D3 calculations capture the experimental in-plane lattice constant and tend to slightly overestimate the out-of-plane lattice parameter by 0.36%. For completeness, we used a van Der Waals density functional (vdW-DF) [46][47][48] to compute the lattice parameters of bulk MoS 2 ; we obtained that the in-plane and out-of-plane lattice constants are overestimated with respect to the experiment by 1% and 0.8%, respectively. To check the influence of the vdW correction on the lattice parameters of bulk MoS 2 , we performed the calculation with GGA only and found that a remains unchanged while c increases by as much as 5% with respect to the experimental value. We concentrate on the DFT-D3 geometries and energies throughout this study.
To model the metal intercalation of MoS 2 , we use the optimized structure to construct a hexagonal 2 Â 2 Â 1 supercell where the desired number of copper or silver is placed in the vdW gap of the bulk MoS 2 . The metal ions can occupy either octahedral sites (H sites) or tetrahedral sites (T sites), as shown in Figure 1. For the supercell size we use, eight configurations of M x MoS 2 can be constructed with x ¼ 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, and 1.0. For each configuration, the metallic atoms are distributed between the two van der Waals gaps maximizing the distance between neighboring intercalates. To better probe the intercalation energy, we generated more diluted systems where only one metal atom was intercalated in 6 Â 6 Â 1, 5 Â 5 Â 1, 4 Â 4 Â 1, and 3 Â 3 Â 1 supercells corresponding to metallic concentrations of x ¼ 0.014, 0.02, 0.031, and 0.053, respectively. All configurations are fully optimized (atomic positions and lattice constants) using the same setup used for the optimization of pristine bulk MoS 2 . For all M x MoS 2 systems, the partial occupancies for each orbital are determined using the order-1 Methfessel-Paxton 49 scheme with an electronic temperature of r ¼ 0.2 eV.

C. Ion mobility
To compute the activation energy for diffusion and minimum energy paths (MEPs) of metal ions diffusing through the vdW gap of MoS 2 , we use the nudged elastic band (NEB) method as implemented in VASP. 50 The NEB calculations start from a set of geometries interpolating between initial and final states; then, the ionic positions of all the geometries are iteratively optimized using only the ionicforce components perpendicular to the hypertangent. The minimum energy path is then determined from spline interpolation of the total energy of the individual images and the tangential projection of the 3N force components on each image. The initial and final geometries are obtained by translating the metal ion between equivalent octahedral or tetrahedral sites within the MoS 2 hexagonal lattice. Three possible paths were optimized using NEB, namely, octahe- , refer to Section III.

D. Electronic transport
The electronic transport through Cu and Ag intercalated MoS 2 is treated in the ballistic regime within the phasecoherent approximation of self-consistent non-equilibrium Green's function (NEGF) method. All electronic transport computations are carried out using the TranSiesta code, 51 implemented in the SIESTA package. [52][53][54] Double zeta plus polarization (DZP) numerical orbital basis sets are used for all atomic species and the exchange-correlation potential is calculated within the GGA-PBE 29 functional. To account for non-local long range interactions, we include a dispersion correction to the total energy as prescribed by Grimme,38 known as DFT-D2. All structures are relaxed at the DFT-D2 level of theory using a maximum force tolerance on each atom of <0.01 eV/Å with a Monkhorst-Pack 45 k-mesh of 15 Â 15 Â 3 points for the ionic relaxation.
The device model for the electronic transport adopts a two probe configuration where the central region is modeled as a rectangular 2 Â 2 ffiffi ffi 3 p Â 1 supercell built from the previously fully relaxed M x MoS 2 geometries (Sec. II B). The left and right electrodes have the same composition as the central region, as discussed in detail in Section IV B. A k-mesh of 4 Â 4 Â 100 was used for the computation of the electrodes self-energy while the transmission spectroscopy was carried out on an 8 Â 8 k-grid.

A. Energetics and structure of intercalation
To assess the energetics associated with the intercalation of metallic atoms into the vdW gap of MoS 2 , we computed the intercalation energy per metal atom where E M x MoS 2 is the total energy of the intercalated system with metal concentration x, n x is number of intercalants in the simulation cell, E MoS 2 is the total energy of the pristine host with the same number of formula units as the intercalated system, and E metal is the total energy per atom of the metal in its ground state crystal structure (fcc). The upper panels in Fig. 2 show the intercalation energy of metal intercalated MoS 2 as a function of concentration for Cu [ Fig. 2(a)] and Ag [ Fig. 2(c)]. Intercalation of a single Cu atom in a 6 Â 6 Â 1 MoS 2 supercell from the fcc metal requires an energy 1.2 eV; additional atoms require less energy and the average intercalation energy for x ¼ 1 is 0.8 eV. In the case of Ag, the concentration dependence is significantly stronger; intercalating the first Ag atom requires approximately 2.2 eV; while for x ¼ 1, the average energy is just above 0.5 eV.
To better understand the nature of the interaction between the ions and the TMDC, we characterized how the intercalation modifies the lattice parameter of the TMDC and charge transfer between the metal and host. As the metal concentration increases, the in-plane lattice parameters increase slightly, but, as expected, the largest expansion occurs in the out-of-plane direction, as shown in Figs. 2(b) and 2(d). Atomic charge of the intercalated species were calculated via Bader Atoms in Molecules (AIM) analysis; 55,56 and are shown in Table I. In the case of Cu, the charge at the lowest concentration studied is þ0.57 e and þ0.51 e for H and T site occupancy, respectively. We observe a weak decrease in ionicity with concentration (to þ0.48 e in the H site and þ0.44 in the T site for x ¼ 1). This is consistent with the relatively weak dependence of the c lattice parameter with concentration, see Figure 2  This increase in metallicity and, thus, the increase in atomic size result in a significant expansion of the c lattice parameter with concentration. The reduction in atomic charge as more Ag atoms are intercalated also reduces the Coulomb repulsion (unlike the case of Cu); this is consistent with the significant decrease in Ag intercalation energy with concentration.
In order to obtain additional insight into the concentration dependence of intercalation energies, we computed the ideal binding between the intercalants in the absence of the TMDs. These calculations were performed by removing the Mo and S atoms from the relaxed simulation cells and calculating the total energy of the Cu/Ag atoms in their intercalated positions (no relaxation). We find that increasing concentration from x ¼ 0.125 to x ¼ 1 results in a reduction in energy (that is, more binding) of approximately 2 eV per atom for both metals; tables with the calculated energies are included in the supplementary material. This shows that interaction of Cu with the TMD layers, including charge transfer, has a strong effect on the interaction between intercalants.
The energetics and partial charges are consistent with the fact that Cu is more electroactive than Ag and indicate that these metals could be electrochemically intercalated into MoS 2 . For comparison, electrochemical dissolution of Cu in amorphous silica is known to occur in resistive switching electrochemical metallization (ECM) cells [57][58][59] and the energetics associated with the process involve a dissolution energy of 2.7 6 2.4 eV. 60 Given the similarity of the intercalation energies of Cu/Ag at H and T sites, we would expect these competing metastable sites to be both occupied. This is supported by the calculation of the ion migration energy barrier through H and T, which differs by only a few meV, see Section II C.

B. Ion mobility
We now switch our attention to the mobility of copper and silver ions within the vdW gaps of MoS 2 . With knowledge of the stable geometries of the ions in the MoS 2 host, the activation energy for diffusion was computed by finding the minimum energy path (MEP) between two equivalent positions within the lattice, as explained in Section II C.
As shown in Figure 3(a), the diffusion of Cu between T sites involves a curved path that approaches the H sites and involves two transition states with an activation energy of 0.32 eV. Note that the H site in this path is not a transition state but rather a metastable configuration; once a copper ion falls in this state, it requires about 0.18 eV to complete the full T-to-T hop. For completeness, we studied a direct path between H sites involving the bridge state (denoted B) located on the projection of a Mo-S bond, see Figure 3(b). We find that this path yields a higher energy barrier than the zig-zag path alternating T and H sites and this should not play a significant role in ion transport.
The minimum energy path found by the nudged elastic band is consistent with the computed intercalation energy profiles shown in Figure 2 where we observed that intercalates seating in T sites are more stable. Ag diffusion involves similar paths but higher activation energies, 0.38 eV between T sites and 0.61 eV between H sites as shown in Figure 4.
The calculated activation energies for diffusion of copper in MoS 2 (0.38 eV) are comparable to that of ionic conductors used as solid electrolytes. For example, the activation energy for diffusion of Li þ in the LMPX (M ¼ Ge, Si, Sn, and Al and X ¼ O, S, and Se) family of super ionic conductors has been reported to be as low as 0.18 eV for M ¼ Ge and X ¼ S and as high as 0.36 eV for M ¼ Ge and X ¼ O. [61][62][63] Other DFT studies 64

A. Band structure and density of states
In this section, we study the effect of copper and silver intercalation on the electronic structure of the material Figure 5(a) shows the electronic band structure of pristine MoS 2 along the high symmetry points shown in the Brillouin zone in Figure 5(b). The shaded regions highlight the electronic band structure along the out-of-plane directions C-A, L-M, and K-H. As is expected, in pristine MoS 2 , we observe very few bands close to the Fermi level along these out-ofplane (vertical) directions due to the weak interaction between layers. Figures 6 and 7 show the electron band dispersions and densities of states of Cu-and Ag-intercalated TMDC at different concentrations, respectively. In all cases, we observe a semiconductor-to-metal transition upon intercalation. The general features of the electronic band structures of the intercalated systems are similar to that of neat MoS 2 . The main effect of intercalation is bringing the conduction band closer to the Fermi level, eventually crossing it at several points in the Brillouin zone. As the metallic concentration increases, more bands cross the Fermi energy, resulting in an increase in the density of states at the Fermi energy and conduction as will be shown in Section IV B. These effects are relatively independent of the location of the intercalant (T or H sites).
Additionally, we note that a significant amount of bands crosses the Fermi energy at multiple points in momentum space along out-of-plane directions.
In order to obtain additional insight into the effects of the intercalation on the electronic structure, we projected the electronic density of states on the various elements and the results are shown as green, red, and blue curves in the DoS plots for Cu/Ag, Mo, and S, respectively. Interestingly, the states near the Fermi level are contributed by molybdenum d and sulfur p bands, while the intercalated metal derived states are well below the Fermi energy.
Despite the qualitative similarities in the electronic band structures between H and T intercalated systems (upper and lower panels in Figures 6 and 7), the different atomic arrangements lead to differences in the density of states at the Fermi energy with consequences in the electron transport properties of the intercalated systems. We find a higher density of states at the Fermi energy for ions in the H sites than in the T sites. For example, for x ¼ 1 in H sites the density of states at E F is 20 and 10 states/eV for Cu and Ag, respectively, while for the T we find 10 and 2 states/eV. Therefore, the degree of metallization is higher for Cu/Ag ions at the H sites than for the T sites, which correlates with higher charge transfer, see Table  I. Also worth noting is the relatively low dispersion of the bands in the out-of-plane direction when Cu is intercalated in T sites, this will result in lower vertical conductance.

B. Electronic transport
To characterize the semiconductor to metal transition due to intercalation, we studied the electronic transport through the intercalated MoS 2 within the ballistic regime using the DFT non-equilibrium Green's function (NEGF) 51 formalism as implemented in TranSIESTA. 65 Figure 8 shows the two-probe configuration used for the transport calculations. The left and right leads of the systems have the same composition as the central region to avoid the effect of interfaces. The size of the device corresponds to a 1 Â 1 Â 3 rectangular supercell of the intercalated systems described in   1. The supercell is oriented in such a way that the c lattice vector matches the direction of transport. In order for the left and right leads to have the same composition as the channel, we use 1 Â 1 Â 1 rectangular supercells of the metal intercalated MoS 2 as extended electrodes. The self-energies of the electrodes are extracted from a standard DFT calculation and the non-equilibrium charge density matrix is computed using NEGF. Figures 9 and 10 show the equilibrium conductance as a function of copper and silver concentration for the T and H cases, respectively. In all cases, we show horizontal (inplane) and vertical transport. We find that conductance increases with concentration of the intercalant except in the case of vertical transport for Cu on T sites where we observe a saturation. These results correlate with the observed increase in DoS at the Fermi energy with increasing concentration. In most cases, Cu intercalation results in higher conductance than Ag, except in the case of vertical transport across the Cu intercalated at T sites which is consistent, as discussed above, with the almost dispersion-less bands near the Fermi energy in the out-of-plane directions C-A, L-M, and K-H, as shown in the upper panel of Figure 6. In contrast, the band dispersion in the in-plane directions of Cu(Ag) intercalated in H(T) sites is comparable to that of the out-of-plane directions leading to a similar behavior of equilibrium conductance. The fact that the cross-plane electronic transport is comparable to that of the in-plane transport opens the possibility of designing top electrodes for transistor applications.

V. DISCUSSION AND CONCLUSIONS
We presented an ab initio study of metal intercalation of Cu and Ag into MoS 2 based on density functional theory. Intercalating Cu at low concentrations requires approximately 1.0 eV per ion, while Ag requires over 2 eV; these values indicate the possibility of electrochemical intercalation. We also note a decrease in the ionic character of the intercalated metal as its concentration increases. Even though the charge transfer of copper and silver at low concentration is comparable, silver becomes practically metallic (neutral) at high concentration (x ¼ 1.0).
Nudged elastic band calculations indicate relatively small activation energy for diffusion of copper and silver in the vdW gap between TMDC layers, ranging from 0.3 eV for Cu to 0.4 eV for Ag. Finally, the simulations show that intercalation of Cu or Ag in MoS 2 induces a semiconductor to metal transition. The origin of this transition is that charge transfer from the metallic ions pulls the conduction band of the TMDC lower in energy and across the Fermi level. The conductance depends on concentration and the site of the intercalant.
The ability to electrochemically intercalate metallic ions into TMDCs could be of interest for several potential applications such as electrochemical metallization cells (ECMs) and contacts for TMDC based electronic devices. ECM cells are resistance switching devices which are of interest for memory and logic. [66][67][68][69][70] These devices switch between low and high resistance states via formation of a metallic filament as an active electrode, which is electrochemically dissolved into a dielectric. Other DFT studies have found that the filament formation is attributed to the increase of metal coordination and strong binding energy between metal ions. 71,72 Current ECMs use either an oxide (like SiO 2 , Al 2 O 3 , and WO 3 ) 73,74 or a chalcogenide solid electrolyte (such as GeS, GeSe, and Ag 2 S) [75][76][77] as the dielectric and layered TMDCs could represent an interesting alternative. Also related to resistance switching memories, materials whose resistance decreases sharply at a threshold voltage (so called threshold switching) are of significant interest as access devices in high-density memory. 78,79 We note that in this study we focused on configurations with evenly distributed intercalants, see Figure 1. We believe a study of the possible aggregation of intercalants, including dimerization and the possible formation of one dimensional paths would be very valuable Metal-doped TMDCs also have the potential to address a major challenge in TMDC-based electronics, namely, the fact that they tend to form high Schottky barriers (SBs) with most of the commonly used metal contacts. Several solutions to this challenge have been explored including doping the contact regions to reduce the SB or selecting contact materials with appropriate work functions. 80 Several metals have been proposed as contact materials to monolayer TMDC as well as different contact configurations leading to the conclusion that edge contacts are more desirable over top contacts due to an enhancement of the electron injection into the TMDC monolayers. 81,82 The development of techniques and novel materials to accurately control contact resistance is key for the realization of TMDC applications in electronics 83 and spintronics. 84,85 The intercalation of metal ions in the contact regions could address some of these issues.

SUPPLEMENTARY MATERIAL
In the supplementary material, the phonon dispersion of copper intercalated systems with concentration x ¼ 0.125 is shown in Sec. S1 along with the partial phonon density of states (Fig. S1). In Sec. S2, the binding energy of the intercalants (Cu and Ag) is presented in Table S1.

ACKNOWLEDGMENTS
This work was supported by the FAME Center, one of the six centers of STARnet, a Semiconductor Research Corporation program sponsored by MARCO and DARPA. The computational resources from nanoHUB.org and Purdue University are also acknowledged.