Doping-Induced Enhancement of Hydrogen Evolution at MoS2 Electrodes
Sander Ø. Hanslin Hannes Jónsson Jaakko Akola

1em




[
Power-grid modelling via gradual
      improvement of parameters
    I. Papp
    September 4, 2024
================================================================



	
]

		













§ INTRODUCTION



MoS2 has emerged as a promising candidate among earth-abundant compounds to replace the precious metal catalysts traditionally used for the hydrogen evolution reaction (HER) <cit.>. Its two-dimensional layered nature allows for novel engineering on the nanoscale, for example by maximizing the presence of edge sites, which display higher activity than the basal plane. 

Transition metal doping has been successful in improving the reaction rates <cit.>, but different and sometimes conflicting results raise interesting questions regarding the underlying activation mechanism. For example, Deng et al. <cit.> found significant reduction of the HER overpotential with Pt-doping of few-layer MoS2 samples, while Co and Ni showed signs of weak activation and deactivation, respectively. Lau et al. <cit.> found that Co-doping lead to activation with respect to undoped MoS2, while Fe, Ni, and Ag were detrimental for the electrochemical current density. Wang et al. <cit.> found reduction in the HER overpotential for Fe, Co, Ni, Cu in vertically aligned MoS2, while Humphrey et al. <cit.> conversely found that Co-doped MoS2 displayed a larger overpotential on basal-oriented MoS2. 

In light of these rich properties of doped MoS2, it is of interest to determine which mechanisms are responsible for the experimentally observed activation or deactivation. The experimental situation is complex, and the doping, especially at high levels, may lead to significant changes in the sample morphology, or other large-scale modifications such as phase transitions, phase separation or formation of impurity particles or clusters on the MoS2 substrate. Such effects are relevant for catalyst performance, but constitute a regime which is challenging to assess theoretically due to the wide scope. In this work, we therefore limit our focus to the relatively low concentration doping regime, considering single- or few-atom impurities in the 2H phase of MoS2. We note that multi-elemental codoping is also a promising approach towards increasing MoS2 activity <cit.>, however such synergistic effects are not explored herein. 

Activation with respect to HER is often attributed to improved values of the H-adsorption free energy (Δ G_H) on both the basal plane and sulfur-terminated edges. However, theoretical works have shown that the nature of hydrogen evolution via Mo-sites and S-sites differs, as in the latter case a large Heyrovsky activation energy must be overcome <cit.>. This challenges the common descriptor-based view as metal-bound hydrogen contributes to evolution more readily despite the near-thermoneutral adsorption in both cases. Thus, improved HER via activation of sulfur-sites is an unlikely mechanism. A further implication is that calculations of activation energy are necessary to predict the HER performance at different dopant sites. 

Hereby, we investigate the consequences of transition/noble metal impurities (Co, Ni, Pt) on the HER kinetics through theoretical reaction modelling at the atomic scale. Our starting point is to consider the active sites of the undoped edge- and basal-oriented MoS2, which are respectively the sulfur-depleted Mo-edge (Mo_0) and S-vacancies, and their modification via doping. Overall, the results demonstrate how the presence of impurities can lead to reduction in the overpotential, and provide insight into the role of dihydrogen intermediates.





§ METHODS


The electronic structure calculations were performed at the level of density functional theory (DFT) within the plane-wave formulation implemented in the Vienna ab initio simulation package (VASP) <cit.>, with core electrons described by the projector-augmented wave approach <cit.>. All calculations were spin-polarized, with plane wave basis sets up to 400 eV and 3×3×1 Monkhorst-Pack sampling of the Brillouin zone. Valence electrons for transition metals include the outer s- and d-electrons. The revised Perdew-Burke-Ernzerhof exchange-correlation functional <cit.> was used, along with the D3 dispersion corrections <cit.>. Transition states for elementary reaction steps were approximated by first-order saddle points on the potential energy surface, obtained by the climbing-image nudged elastic band method <cit.>. Minima and saddle points were optimized with force convergence criteria of 0.02 eV/Å and 0.05 eV/Å, respectively. 

The electrolyte is represented by an explicit Eigen cation water cluster, as well as an implicit polarizable continuum model as implemented in VASPsol <cit.>. For the final kinetic evaluations, the grand-canonical reaction and activation energy is evaluated so as to account for the constant electrode potential during the electrochemical reaction. This is done by allowing the number of electrons in the system to vary, implicitly tuning the electrode potential <cit.>. The grand-canonical energy is then defined as Ω = E_n + δ n e Φ where E_n is the DFT energy with n electrons, δ n = n-n_0 is the difference in number of electrons from the neutral state, e is the elementary charge and Φ is the electrode potential, given by the effective work function. In the following, we refer the electrode potential to the standard hydrogen electrode (SHE), defining U=Φ-Φ_SHE with Φ_SHE = 4.43 V. During the grand-canonical evaluation, the saddle point structure representing the transition state is kept fixed.

The basal plane is modelled as a single monolayer, represented by a periodically repeating 5×5 MoS2 unit cell (75 atoms excluding adsorbates and solvent species). The edge model consists of alternating layers terminated with 50 % and 0 % S-coverage on the S- and Mo-edges, respectively, resembling a slab of vertically aligned 2H-MoS2 sheets. Each layer is 3×4 MoS2 unit cells for reaction calculations, and 5×4 for formation energy and adsorption calculations (64 and 112 atoms per simulation cell). Single transition metal atoms are introduced as Mo-substitutional impurities to these model systems, yielding a basal plane doping concentration of 4 % in terms of Mo-substitution, and 25 % or 100 % substitution of the Mo edge. The inset in Figure <ref> shows a top-view of the two simulation cells.

Formation energy calculations are performed via the cohesive energy

    E_c = E_tot - ∑_iE_i n_i,

where the sum goes over unique atomic species i with energy E_i and of quantity n_i in the simulation cell. 

The free energy of hydrogen adsorption in the gas phase is used for initial characterization of different adsorption sites, and is given by 

    Δ G = E_nH-E_(n-1)H - 1/2 E_H_2+ Δ E_zpe - TΔ S,

where E_nH is the DFT energy of a state with n adsorbed H-atoms and E_H_2 is the DFT energy of the H2 molecule. Δ E_zpe denotes the difference in zero-point energy upon adsorption, and is determined by vibrational analysis. The entropic contribution TΔ S is estimated by the standard entropy of the H2 molecule, effectively neglecting the entropy of the adsorbed state. 

The concentration of sulfur vacancies and hydrogen adatoms directly depends on the electrochemical conditions, specifically the pH, electrode potential U as well as the presence of sulfur species in solution. We assume that desulfurization occurs via electrochemical H2S production, so that the sulfur chemical potential is

    μ_S = μ_H_2S - 2(μ_H^+ + μ_e^-),

where the chemical potential of the proton-electron pair in solution is defined via the computational hydrogen electrode approach <cit.>:

    μ_H^+ + μ_e^- = 1/2μ_H_2 - eU + k_BT ln a_H^+,

where pH=-log_10 a_H^+. As for H2, μ_H_2S is given by the DFT energy, zero-point energy and standard entropy. Acidic conditions (pH=0) and a partial H2S pressure of p_H_2S/p^∘_H_2S = 10^-8 are assumed. These definitions of sulfur and proton-electron pair chemical potentials are used to evaluate relative energies of hydrogenated sulfur vacancies on the basal plane, as a function of electrode potential. 

A kinetic model is used to evaluate the relative reaction rates of the different configurations, and construct simulated polarization curves. We consider one site at a time, and separately evaluate the Volmer-Tafel and Volmer-Heyrovsky pathways. Assuming constant concentration of solvated species, the state space is given by the possible surface hydrogen configurations for the given site (0 or 1 H for the Heyrovsky pathway, 0, 1 or 2 H for the Volmer-Tafel pathway). Transitions between different microstates occur via the elementary Volmer, Heyrovsky and Tafel steps, and reverse reactions involving H2 are neglected, i.e. H2 is considered to be removed from the local system as soon as Tafel or Heyrovsky combination has occurred. A general transition rate from microstate a to b via path s is given by r^s_ab = θ_a k^s_ab, with θ_a being the occupation of state a. Within transition state theory, the associated rate constant takes the form

    k^s_ab = k_BT/h e^-Ω^(s)_ab / k_B T,

where k_B is the Boltzmann constant, h is the Planck constant, T is temperature, and Ω^(s)_ab is the grand-canonical activation energy, i.e. the energy of the transition complex along s connecting the a and b minima, referred to that of the state a. T=300 K is used in rate calculations. To account for rearrangement of the solvent beyond those included in the simulation, a minimum limit of 0.2 eV is enforced for the Volmer and Heyrovsky forward/reverse barriers. The total transition rate from a to b is then a sum over the available pathways, r_ab = ∑_s r^s_ab.  Solving for the steady state condition, where the occupation of each state is constant, the steady state current density due to electron transfer in the forward (+) and reverse (-) Volmer (v) and Heyrovsky (h) steps is given by 
    j = -e/A∑_ab( r_ab^v(+) - r_ab^v(-) + r_ab^h(+)),

with the elementary charge e and effective area per site A. j is voltage-dependent via the grand-canonical activation energies, and with this knowledge one can construct theoretical polarization curves. The area per site is here given by the number of sites (1 for basal plane, 1-4 for edges) and lateral dimensions of the chosen simulation cell. Our goal is not a direct comparison of the absolute magnitude of the current density with experiments, but nonetheless we should note that the edge content may vary considerably in experiments, depending on the sample morphology. The Mo_0 edge content of this simulation cell is ca. 0.8 nm/nm^2, which is of comparable magnitude to experimental values <cit.>. As the current density is proportional to the site density, the exponential nature of the polarization curve somewhat mitigates the importance of the chosen area on determining the overpotential. When discussing the theoretical overpotential η in the following, we refer to the magnitude of the negative electrode potential at which j exceeds the arbitrary threshold of 10 mA/cm^2. 





§ RESULTS AND DISCUSSION



Considering first the relative free energy of the impurity formation (Figure <ref>), we find that substitution on the edges is the most favored thermodynamically for all three metals, and especially for the Mo_0 edge. Thus, it is more likely to occur on edge-locating Mo. Further, interaction with sulfur-vacancies stabilizes the basal plane impurity. Defining chemical potentials from the respective bulk metal phases, the absolute basal plane substitution energies (zero-levels in Figure <ref>) are 3.8 (Co), 4.3 (Ni), and 4.5 eV (Pt). However, conclusions regarding absolute stability warrant a more thorough thermodynamic study, considering other possible reference phases in appropriate conditions. Hence we consider the relative formation energy to evaluate the different sites for each impurity metal. In this regard, Co, Ni, and Pt all show similar behavior. Based on these energies one might expect a significant edge-substitution at thermodynamic equilibrium, if no phase separation occurs. 

As indicated by horizontal lines in Figure <ref>, the formation energy per impurity atom remains low in the case of complete edge substitution, suggesting that a complete filling of the Mo_0 edge is more favorable than occupying other sites. Moving forward, we assume from this that the doped Mo_0 edge is present in the alternating vertical layer model. Note that the selection of the Mo_0 and S_50 edges qualitatively corresponds to sulfur-poor conditions, and that the selectivity of the substitution in general depends on the chemical potential of sulfur, as well as the impurity element <cit.>. 

The differences in substitution energy between basal plane and edges are in good agreement with those of Ref. , where it was also found that doping of the Mo_0 edge was significantly more favorable than with 50 % or 100 % terminating S-coverage. The presence of impurity atoms on the Mo edge would then stabilize the sulfur-depleted configuration to some degree.



 §.§ Doping of vacancies on basal plane


Looking further into the interaction between S-vacancies and the impurity metal, we note first that the single vacancy is stabilized by almost 2 eV if situated next to an impurity atom in comparison with the undoped basal plane. This is in good agreement with an earlier theoretical work <cit.>. In thermodynamic equilibrium the impurities will then be accompanied by S-vacancies. Considering higher levels of S-deficiency and interaction between vacancies, neighboring vacancy pairs are weakly favored over dispersed configurations by ca. 0.1 eV on the undoped MoS2 basal plane, see Figure <ref>a. 

In the presence of impurity atoms, however, the neighboring configuration is strongly favored by ∼1 eV for all dopants. Thus, the impurities enable neighboring vacancy configurations to a larger degree. Further, for the third vacancy a cluster-like configuration is favored for Ni- and Pt-doping, while Co-doped and undoped MoS2 favor linear vacancy distributions, see Figure <ref>b. The cluster triple vacancy configuration fully exposes a metal atom and is likely of interest for HER. Overall, the vacancies and impurity atoms are co-confining. 

On undoped MoS2, hydrogen atoms can adsorb in the single vacancy site with a near-zero change in free energy. The multiple-vacancy configurations and impurity atoms modify the adsorption, as shown in Figure <ref>. Notably, Ni and Pt do not modify adsorption onto the single vacancy (V_S) significantly, while Co leads to a more favorable adsorption. In the doped vacancy, the second adsorption is significantly less favorable, but the adsorption configuration consists of H2* adsorbed onto a Mo-atom, which suggests a possible Volmer-Tafel pathway. The importance of dihydrogen intermediates has been established for single-atom catalysts <cit.>, and similar trends may be relevant on the local doping-induced structures considered in this work. In this specific case the intermediate seems too unstable to support an efficient pathway.
In the double-vacancy (V_2S), adsorption becomes stronger by doping, and these configurations are likely not of interest for hydrogen evolution. In the clustered triple vacancy configuration (V_3S), undoped MoS2 has a near-thermoneutral adsorption onto the top-site once the more favorable vacancy sites are filled. The presence of an impurity atom makes the top site less favorable. A bridge-configuration is preferred instead, but is still high in energy. Notably, a fifth hydrogen atom may adsorb (unfavorably) onto the undoped top site to form a dihydrogen complex. On the doped structures such a complex would dissociate, but even then the adsorption is too endothermic. Overall, the doping seems to bring the vacancy sites out of the thermoneutral regime for hydrogen adsorption. 





Once the MoS2 system is subjected to electrochemical HER conditions, the potentially stabilizing effect of hydrogen adsorption is crucial in determining which vacancy configurations are formed. Under desulfurizing conditions (pH=0 and p_H_2S/p^∘_H_2S = 10^-8), Figure <ref> shows the free energy of hydrogenated vacancy configurations in doped and undoped MoS2. Single and double vacancy formation is favored at small applied negative voltages for the doped systems, while the triple vacancy is preferred below ca. U=-0.35 to -0.40 V. In undoped MoS2, the fivefold hydrogenated triple vacancy configuration surpasses the pristine basal plane around -0.55 V. At larger negative potentials, higher vacancy concentrations would be preferred, but importantly we note from these results that due to favorable hydrogen adsorptions, the clustered triple vacancy is favored over the linear or zig-zag variants in all systems. In the two latter cases, the adsorptions are higher in energy, more similar to those on the double vacancy. 
Given that the system can equilibrate under these conditions, we thus expect the clustered triple vacancy to be present (and favored) even in undoped MoS2, and going forward we will consider the cluster vacancy for all systems. Tsai et al. <cit.> observed clustered vacancy configurations after electrochemical desulfurization of basal-oriented (undoped) MoS2. Further, their observed onset in reduction of the S:Mo ratio between U=-0.5 and U=-0.6 V vs. RHE  coincides well with the stable region of V_3S for the undoped system in Figure <ref>.

Next we consider the single, double and triple (cluster) vacancies, and obtain the grand-canonical reaction- and activation energies via reaction modelling as outlined in the Methods section.  For the single-vacancy, a Volmer-Heyrovsky mechanism via the first hydrogen is preferred, despite the Tafel-relevant geometry upon double adsorption. The energy diagram for this mechanism (see Figure <ref>a), shows that the Heyrovsky barrier is significantly reduced by doping. The simulated polarization curves are presented in Figure <ref>b. As expected from the energy diagram, all dopants activate the single-vacancy with respect to the undoped case, reducing the required overpotential by 0.3-0.4 V. 



On the double vacancy configuration, the Volmer-Heyrovsky pathway is preferred, but only Co shows any significant activity. The Heyrovsky barriers are larger than in the single vacancy due to the more favorable adsorption, leading to the current onset at a more negative potential. Interestingly, the Heyrovsky barrier becomes smaller than the Volmer barrier at larger negative potentials, and continues to scale while the Volmer barrier stagnates. This leads to the flat shape of the Co-MoS2 polarization curve, as the (more) rate-determining step is not being significantly reduced. For reference, calculations of the Tafel mechanism on these doped double-vacancy configurations showed that the barrier for Tafel combination is very large (2 eV upwards) due to repulsion from the impurity atom, rendering this mechanism irrelevant. The double vacancy thus strictly inhibits evolution in both doped and undoped cases.

For the triple vacancy, different mechanisms are supported in doped and undoped systems. A Volmer-Volmer-Tafel mechanism via a dihydrogen intermediate is possible in the undoped case, while the endothermic H-adsorption of the doped systems favors a direct Volmer-Heyrovsky mechanism (see Figure <ref>a). The corresponding energy diagrams are shown in Figure <ref>b, where the required activation energy is significantly smaller in the undoped case. This leads to the low overpotential in Figure <ref>b. The Co- and Pt-systems display very similar overpotentials compared to the single vacancy case, while the Ni-doped system is limited by a large Volmer-barrier and is deactivated in this configuration.
Importantly, formation of a dihydrogen complex enabled a more efficient reaction path, which is further explored in the following section on edges. 



The improved kinetics of the clustered triple vacancy on undoped MoS2 has some important implications. Due to the large overpotential required on the undoped single vacancy, triple vacancies would form before this point according to the analysis in Figure <ref>, indicating that the origin of basal plane activity is not only any S-vacancy, but specifically the undercoordinated Mo-atoms in threefold vacancy configurations. This is in contrast with the typical understanding that the single vacancy itself enables evolution due to its near-thermoneutral hydrogen binding. Comparing again with work by Tsai et al. <cit.>, the onset of reduced S:Mo atomic ratio around U=-0.6 V vs. RHE was accompanied by significantly increased HER activity. This aligns well with the understanding provided by our results, namely that the fully exposed Mo atom is the active site on the basal plane. This configuration is generated by the same conditions that drive the hydrogen evolution, and connection with the initial vacancy concentration and distribution in the sample may be elusive. 

Doping deactivates the triple vacancy (relatively for Co and Pt, and completely for Ni), which suggests that the basal plane is deactivated by doping overall. Doped single-vacancies are activated by a reduction of 0.3-0.4 V in the overpotential, and are present at lower voltages than in the undoped case. However, they are only stable at very small voltages, before they are replaced by doped double- and triple vacancies, both of which inhibit hydrogen evolution compared to the undoped triple vacancy. Experimental work by Humphrey et al. <cit.> found reduction in HER activity upon doping basal-oriented MoS2 with Co, in agreement with these conclusions. Rather than deactivating the single vacancy, our results suggest that this is due to deactivation of the triple vacancy.



 §.§ Doping of edges


We have chosen to focus on the Mo_0 edge, as this is both the most kinetically active edge <cit.> as well as the most thermodynamically stable doping site. Two cases of edge doping levels are considered, 25 % and 100 %. At 25 % edge doping, the equilibrium H-coverage is reached at 1.25 monolayers where H is bound to Mo-atoms, avoiding the slightly higher energy impurity site. From there on, the next adsorption can occur onto the impurity atom or onto one of the Mo atoms, in most cases with a weakly positive free energy cost. The exception is adsorption onto the Ni-atom, which is unfavorable by ca. 0.5 eV, see Figure <ref>. 

At 100 % edge doping, the first adsorption is (by default) endothermic for Co and Ni, and 0 eV for Pt. However, interestingly, another hydrogen atom can be adsorbed onto the same dopant site, forming an adsorbed H2* complex. The adsorbed complex is stable for Co, and for Ni the second adsorption is nearly thermoneutral with respect to the first. On the contrary, the second adsorption is strongly unfavorable for Pt. The same trend can be seen on the impurity atom at the 25 % doping level. The affinity towards H* and H2* are thus qualitatively different for these metal impurities. We note that the trend of exothermic adsorption to form stable H2* following an initial endothermic H-adsorption is similar to that observed in single transition metal atoms supported on the MoS2 basal plane <cit.>. The H2* complex can also form on Mo-atoms, and in this case the adsorption is closer to thermoneutral also for the Ni- and Pt-systems. The formation of such H2* species is relevant when considering the possible reaction pathways for HER. Such species with near-thermoneutral binding are prime candidates for efficient intermediates. This also illustrates the importance of considering coverage effects and exploring the adsorption configuration space prior to modelling HER itself. 



The possibility of further H-adsorption leading to H2* complexes warrants revisiting the undoped Mo_0 edge, where the evolution was found to proceed through a Volmer-Heyrovsky pathway <cit.>. Interestingly, the incoming proton interacts favorably with the surface in the corresponding Heyrovsky transition state (TS). This attractive interaction stabilizes the transition complex, and leads to geometrically similar TS for the Volmer and Heyrovsky steps which are also close on the potential energy surface (Figure <ref>a). From the Volmer perspective, TS is stabilized by favorable H-H interaction in addition to the surface attraction. This also means that something resembling a H2* complex near the surface is part of the reaction for both the Heyrovsky and Volmer processes. Depending on whether the energetics favor adsorption of the H2* complex, H2 may release directly into solution or stay on the surface, (nominally a Heyrovsky or Volmer step). In the latter case, H2* may then desorb in a Tafel-like process at a later point, completing HER. In the following we refer to these mechanisms as Volmer-Heyrovsky and Volmer-Volmer-Tafel, although the Tafel process only involves desorption, not H-H combination. 

Figure <ref>b displays the grand canonical energy diagrams at U=-0.2 V for the undoped Mo_0 edge at a hydrogen coverage θ_H = 1.25,  starting after the first Volmer step (identical for all cases). The H2* intermediate is seen to proceed with a very small Tafel barrier after adsorption. However, the direct H2-release via a Heyrovsky process is slightly preferred due to the smaller barrier, and proceeds at a lower overpotential, as shown in Figure <ref>c. Note that this distinction may be sensitive to computational specifics, e.g. choice of exchange-correlation functional and dispersion corrections. The Volmer-Tafel process through the H2* intermediate (A) is also compared to that of the alternative 2H* intermediate (B). The 2H* configuration is more stable, but also requires a large Tafel barrier to combine, leading to detrimentally worse overall kinetics. The difference between these two Tafel paths illustrates that a more thermoneutral intermediate is not always associated with better kinetics, as possible pathways and barriers must be considered. The presence of these similar transition- and intermediate complexes suggests that the potential energy surface is quite flat in this immediate neighborhood of the configuration space, encompassing configurations of H2* to H2 via water-surface complexation with comparable energies.



Based on the findings above and considering the adsorption energetics, we map out the Volmer-Heyrovsky and Volmer-Volmer-Tafel pathways through the Mo- and dopant-sites. The polarization curves are presented in Figure <ref>, where we consider hydrogen coverages immediately below and above the equilibrium coverage. The H atoms partaking in evolution are colored blue in the right panels. For both mechanisms, the kinetics on the Mo-site farthest from the impurity atom (Mo^1) are rather similar, which is to be expected. Closer to the dopant atom the behavior is more different. The Volmer-Heyrovsky mechanism is deactivated on both Mo-sites, but activated on the dopant site for 25 % of Co and Pt impurities. The Ni-site is severely Volmer-limited due to the large adsorption energy. At 100 %, all systems are moderately deactivated for Volmer-Heyrovsky, compared to the pristine Mo edge. The kinetics at 25 % match the expectations set by adsorption energies, but this is not the case at 100 %, where Pt performs worse despite Δ G = 0 eV. The Mo-sites also perform worse despite near-zero adsorption energies, demonstrating that the interaction between the protonated water cluster and electronic states of the hydrogenated surface at the transition state can not be inferred from the hydrogen binding energy alone. The Volmer-Volmer-Tafel process is seen to proceed at a reduced overpotential via the dopant site for Co (25 % and 100 %), and via the neighboring Mo-site for Ni and Pt (25 %), where the H2* complex is not favored on the impurity site. 



While considering these different possibilities in the context of experiment, the first remark is that the Mo_0 edge is significantly more active than the basal plane, with ca. 0.25 V lower overpotential in the undoped cases. Therefore, if the edge-content in experiment is not negligible (even small amounts of edge-sites will affect the effective overpotential significantly), e.g. in polycrystalline samples, it is reasonable to assume that the observed doping effect is due to modification of the edge sites. Within this regime, moderate doping levels of Co, Ni and Pt will enhance the edge activity. At high doping levels, Ni- and Pt-doping leads to deactivation while Co-doping still turns out beneficial.

The activation due to Pt by Deng et al. <cit.> and Co by Lau et al. <cit.> leads to shifts in the overpotential of similar magnitude as those found here for the 25 % Mo_0 edge. 
The deactivation due to Ni by Lau et al. is also similar to what we find for the 100 % edge substitution. Deng et al. find also weak activation by Co, yielding a slightly larger overpotential than for Pt-MoS2, while we find the opposite trend in the theoretical activation. The reason for this is not clear, but we should note that fixed wt %-doping constitutes different degrees of atomic substitution for Co/Ni and Pt. Due to the difference in atomic mass, the distribution of Pt-impurities will be more than three times as dilute. Deng et al. use a lower doping level (1.7 wt %) than Lau et al. (3.0 wt %) and see a less pronounced Ni-deactivation. We speculate that this difference is related to the degree of edge-substitution, as it is consistent with the theoretical understanding developed in this work (activation for partial edge substitution, deactivation upon complete edge substitution). The experimental work that most resembles the edge model used here is that of Wang et al. <cit.> where vertically aligned MoS2 layers were doped with 3d metals Fe through Cu. The Co-doped sample contained ca. 22 at % Co on the edge, decaying with depth, and all dopants resulted in overpotential reduction on the order of 0.1 V, in excellent agreement with our findings for the 25 % edge substitution. Rather than activation of the S-edge, our results indicate that this activation can be explained by improved kinetics on the Mo_0-edge. 

It should be noted that any systematic errors due to choices in the theoretical model, e.g. the implicit solvation scheme, explicit interface description, exchange-correlation functional or other computational parameters, may manifest as shifts in predicted barriers and reaction rates. We assume that these errors are similar across the studied systems, so that relative comparison is valid. Comparison with experiment suggests that the absolute rates are also of reasonable magnitude, but we must acknowledge the possibility of errors due to neglected effects, as well as partial cancellation of these.

Unlike the edges, our results show that the 2H-MoS2 basal plane is deactivated by Mo-substitutional doping. This doping mechanism cannot therefore account for the very low overpotentials observed in experiments. However, we can not disregard the possibility that the general basal plane may be activated by means of other mechanisms, notable examples being transition to the 1T phase or anchored single-atom impurities. 

Both the basal plane S-vacancies and the Mo_0 edge involve Mo-bound hydrogen with Δ G_H≈ 0, but the much faster kinetics on the edge (and the significant improvement from single to triple vacancy) emphasizes the significance of exposing undercoordinated metal atoms, and suggests that the resulting hydrogen binding is of different nature on the respective sites. The low coordination of the edge enables formation of the metal-bound dihydrogen complex. Such affinity towards dihydrogen chemisorption is typically associated with transition metal complexes <cit.>, as H2 on metal surfaces tends to either dissociate or physisorb <cit.> (dihydrogen may however bind to defects such as ridges <cit.> or form as transient states during Tafel combination <cit.>). This indicates that the Mo_0 edge represents a middle ground between undercoordination, enabling chemistry resembling a coordination compound, and maintaining the good electron transport and stability of a metallic surface, as can be a problem with e.g. supported single atom catalysts. The hapticity of H2 ligands on transition metal complexes is largely governed by π-backdonation from metal d-orbitals into the antibonding σ^*_H-H orbital, and equivalent mechanisms likely determine the H2* binding energy on the doped edges.  From the few dopants studied herein one can at least note that dihaptic binding on the edge is less favored for dopants with a higher number of valence d-electrons (Co, Ni and Pt being respectively 3d^7, 3d^8 and 5d^9), though more in-depth analysis would be insightful in this regard. As the Kubas interaction fundamentally defines the energy landscape of the H2*↔2H* transition, tuning this interaction may be an important tool also in optimizing electrocatalyst performance.





§ CONCLUSIONS


	
The effect of transition/noble metal doping (Co, Ni and Pt) on the activity of MoS2 towards HER was studied by performing grand-canonical DFT simulations (incl. solvent description) and theoretical reaction modelling which enabled construction of theoretical polarization curves. For the undoped basal plane, the active site responsible for intrinsic activity was found to be the central Mo atom in clustered threefold S-vacancies. Despite the single- and double-vacancies allowing near-thermoneutral H-adsorption, hydrogen evolution via a H2* complex occurs through a Volmer-Volmer-Tafel pathway on the triple vacancy with a much lower overpotential (η>1.1 V vs. η=0.52 V). The clustered triple vacancy is also energetically favored by HER conditions. Impurity atoms interact favorably with S-vacancies in mutually stabilizing configurations, enabling increased vacancy generation. However, H-intermediates on the impurity atom are high in energy, and the resulting Volmer-Heyrovsky pathway requires a larger activation energy. Therefore, all dopants seemingly deactivate the basal plane. 

Like the undoped triple vacancy, the edges display moderate affinity towards surface-bound H2* complexes in both doped and undoped cases. This enables an efficient Volmer-Volmer-Tafel process of evolution in which the H* combination occurs directly in the second Volmer step, and also influences the Volmer-Heyrovsky process by stabilizing the transition complex via attractive surface interaction. At moderate doping levels (25 % edge substitution), evolution proceeds with a reduced overpotential (0.1-0.2 V) for all dopants. At full edge substitution, Ni and Pt are deactivated relative to the undoped Mo_0 reference (η = 0.27 V), and only the Volmer-Volmer-Tafel process on the Co site remains more active, seemingly due to the more stable dihydrogen intermediate, which may further be related to the valence d-electrons via the Kubas interaction. 

The large discrepancy in performance between sites on the basal plane and edges despite similar thermodynamics of the H* intermediate illustrates the importance of explicitly including activation energy in kinetic estimations, as the free energy descriptor does not account for this difference. In experimental context, it suggests that observations of low overpotential in doped 2H-MoS2 are unlikely to be due to Mo-substitutional doping of the basal plane, and rather due to modification of the Mo edge. 

In summary, edge doping can greatly increase the HER activity of MoS2, and the best results are achieved by partial (Co, Ni, Pt) or full substitution (Co) of the Mo_0 edge. The different trends on MoS2 basal planes, low-, and high-level doped edges may help explain experimental observations of differing effect of transition and noble metal impurities.
Further, the role of dihydrogen intermediates was identified, contributing towards understanding the chemical picture of hydrogen evolution in these systems. The behavior of transition metal dichalcogenide edges presents an intersection between surface- and coordination chemistry which seems essential for their role as efficient catalysts. Further understanding of this regime may guide the design of active yet stable catalysts.





§ ACKNOWLEDGEMENTS


The calculations were performed on resources provided by Sigma2 - the National Infrastructure for High Performance Computing and Data Storage in Norway, project No. NN9497K. JA acknowledges financial support from the Academy of Finland, project No. 322832 ‘‘NANOIONICS’’. 
HJ thanks Prof. Jens Nørskov for numerous inspiring discussions and collaborations over the past 30 years, and acknowledges financial support from the Icelandic Research Fund, project No. 207283-053.





§ CONFLICT OF INTEREST


There are no conflicts to declare.





 








Wiley-chemistry






