Electrocatalytic Study for Hydrogen Evolution Reaction on MoS$_2$/BP and MoSSe/BP in Acidic Media

Molecular hydrogen (H$_2$) production by electrochemical hydrogen evolution reaction (HER) is being actively explored for non-precious-metal based electrocatalysts that are earth-abundant and low cost like MoS$_2$. Although it is acid-stable, its applicability is limited by catalytically inactive basal plane, poor electrical transport and inefficient charge transfer at the interface. Therefore, the present work examines its bilayer van der Waals heterostructure (vdW HTS). The second constituent monolayer Boron Phosphide (BP) is advantageous as an electrode material owing to its chemical stability in both oxygen and water environments. Here, we have performed first-principles based calculations under the framework of density functional theory (DFT) for HER in an electrochemical double layer model with the BP monolayer, MoS$_2$/BP and MoSSe/BP vdW HTSs. The climbing image nudged elastic band method (CI-NEB) has been employed to determine the minimum energy pathways for Tafel and Heyrovsky reactions. The calculations yield that Tafel reaction shows no reaction barrier. Thereafter, for Heyrovsky reaction, we have obtained low reaction barrier in the vdW HTSs as compared to that in the BP monolayer. Subsequently, we have observed no significant difference in the reaction profile of MoS$_2$/BP and MoSSe/BP vdW HTSs in case of high coverage (25 %) and 1/3 H$^+$ concentration (conc.). However, in the case of small coverage (11 %) and 1/3 H$^+$ conc., MoSSe/BP shows feasible Heyrovsky reaction with no reaction barrier. Finally, on comparing the coverages with 1/4 H$^+$ conc., we deduce high coverage with low conc. and low coverage with high conc. to be apt for HER via Heyrovsky reaction path.


Introduction
The availability of clean and renewable energy source governs the tenable development. Innovation in systems like fuel cells, metal-air batteries and water electrolysis positively impacts the environment. 1 The cleanest alternative for the same is the molecular hydrogen (H 2 ) and hence, in the present context, we consider materials that support its production. 2,3 The electrochemical reactions that are in sync with the clean environment aim involve hydrogen oxidation reaction (HOR), oxygen reduction reaction (ORR), hydrogen evolution reaction (HER) and oxygen evolution reaction (OER). 4 The former two are associated with fuel cells, while the latter two are associated with water splitting or water electrolysis. There exists wide range of materials that can catalyze these electrochemical reactions by photocatalytic or electrocatalytic pathways. 5-7 The present paper focuses on HER by the electrocatalysts.
HER requires large overpotential to be initiated, and therefore catalysts are required to lower the overpotential. 8 In this respect, Pt has established itself to be an efficient catalyst. 9 However, its high cost and low abundance have urged the scientific community to find new materials for catalytic applications. 10 In fact, any heterogeneous catalysis under periodic boundary conditions faces the challenge of possessing an apt catalytic material that decreases the reaction barrier. 11 HER can occur in both acidic and alkaline media. In either of the media, the reaction steps follow (i) adsorption of H, (ii) its reduction and (iii) desorption as H 2 . 12 Now, HER has been reported to have sluggish kinetics in alkaline media with ambiguous active sites. 13 Since the electrolytic reactions at the electrode are acidic, we are focusing on the acidic media in the present study. The adsorption step is very fast and is termed as the Volmer step: 14 Volmer reaction (fast): H + + e − → H ad The subsequent steps take place either as Tafel or Heyrovsky paths.
As previously mentioned, the concept is to obtain material for the reaction that does not include precious metals like Pt. The literature has shown the transition metals (Fe, Ni, Co), carbides, metal oxides (RuO 2 , IrO 2 ), graphene, non-layered 2D materials, metal-organic frameworks (MOFs) and transition metal dichalcogenides (TMDs) as effective HER catalysts. 8, 12,[15][16][17][18][19] We restrict our study to the 2D materials that showcase quantum confinement effects with increased carrier mobility and large surface area. 20,21 This results in their increased catalytically active sites. The monolayer TMDs (in place of graphene) have established themselves as a potent material with optimal band gap suitable for optoelectronics, photocatalysis and electrocatalysis. [22][23][24][25][26][27] In addition, due to their flexibility, these are widely studied for flexible electronic devices. Literature has reported their use as catalysts for HER, especially on the surface of 1T -MoS 2 and edge sites of 2H-MoS 2 . 28,29 MoS 2 being acid-stable is an added advantage. 30 Furthermore, its heterojunctions have also shown promising HER catalytic behaviour. 13,[31][32][33] It is pertinent to mention here that the tunability of 2D materials for specific applications is prevalent by defect engineering, strain engineering, stacking order, external field implementation, alloying and forming heterojunctions. [34][35][36][37][38] Amongst them, formation of heterojunctions with van der Waals forces in between the constituent monolayers are classified under van der Waals heterostructures (vdW HTSs). These have proved a real boon to the field of work because the constituent monolayers retain their properties simultaneously with their combined vdW HTS properties. In addition, the electronegativity difference between the constituent monolayers actuates electron transfer, thereby affecting the HER. 39 Even if the constituent monolayers have inactive sites, the resulting vdW HTS can be obtained as an active electrocatalyst due to an inbuilt electric field at the interface. 40 Presently, we explore Boron Phosphide (BP) monolayer, MoS 2 /BP and MoSSe/BP vdW HTSs for HER. Recent works have reported vdW HTSs with BP instead of graphene as it has a similar single atomic layered hexagonal structure, however, along with a band gap. 41 BP monolayer has been reported with low carrier effective mass, high carrier mobility, good mechanical strength, and stability in water environments. 42,43 Since the lattice parameter of MoS 2 and BP is similar, the MoS 2 /BP vdW HTS becomes a plausible system with minimal lattice mismatch. 44 BP monolayer has also been synthesized experimentally . 45 In addition, since Janus (MoSSe) has established itself with more catalytically active sites than MoS 2 , we have also analyzed MoSSe/BP vdW HTS. Any prior investigations for HER on these systems are hitherto unknown; hence we have considered these systems for our work.
The aforementioned HER reaction path should be accounted for the proton and electron free energies. These are incorporated by the computational hydrogen electrode model as proposed by Norskov et al. 14 The model caters to the fundamental problem of large-scale calculation of a real system along with electrolyte by following the electrochemical double layer approach rather than external charge formation. The underlying approximation considers solvated proton upto first bilayer. Until now, no study has been reported with the analysis of vdW HTSs by computational hydrogen electrode model for HER in acidic media to the best of our knowledge. 46 We have initially discussed the stacking configuration and electronic structure.
Subsequently, the computational hydrogen electrode model is discussed. Thereafter, Tafel and Heyrovsky reaction paths are analyzed. Finally, we discuss the electrode potential, and the reaction and activation energies.
The generalized gradient approximation (GGA) that accounts for the exchange-correlation and it corroborates with the prior research. 44,66 Note that, initially, two stacking styles (see Fig. 2 (a) and (b)) between the constituent monolayers were considered, wherein the stacking corresponding to Fig. 2 (b) has minimum binding energy. 31 Therefore, we have proceeded with this stacking in our work.

HER Study
Now, we advance on HER study, for which we have constituted 2×2 and 3×3 supercells. The former being smaller restricts the proton concentration (conc.) variability, therefore we need larger supercells. In view of this, we are analyzing 2 × 2 supercell along with 3 × 3, because unlike monolayer, the vdW HTS with further large supercell size becomes computationally demanding. The subsequent paragraphs discuss the concepts of coverage and proton conc.
for clarity.
The first step is to obtain the coverage that gives ∆G H 0 for our study. The number of adsorbed hydrogen (H ads ) per surface atoms is defined as the coverage. ∆G H is the free energy of atomic hydrogen adsorption and is expressed as: where, In the aforementioned equations, ∆E H is the hydrogen binding energy on the surface of vdW  Fig. 2 (f)). We now discuss the optimized systems consisting of water layer (water-solid interface with 3 Å thick water layer) without and with solvated protons (i.e., H + ). Fig. 2  are essential because the electrostatic potential, as seen from the solid surface, also depends on the same. The stability of the vdW HTS along with water layer orientation is established by the similar profile of radial distribution plot at 0K and 300K (see Fig. S6 in SI).   instead, a slight hump is present. This corresponds to the buckling in the BP monolayer.

Tafel Reaction Step
As previously mentioned, the site of H ads is buckled with respect to other sites, and during the H 2 evolution process, the corresponding BP site adjusts itself to the planar configuration (see Fig. 5). Note that BP surface is considered for the reaction analysis as basal plane of MoS 2 is not catalytically active.

Heyrovsky Reaction Step
Unlike Tafel, the Heyrovsky reaction step is not a pure surface reaction. It involves charge transfer, thereby affecting the reaction barrier and Fig. 4 (e)-(h) substantiates the same.
We have observed reduction in reaction barrier in the vdW HTSs as compared to that in BP monolayer (refer Fig. 3 and Fig. 4). The 3 × 3 supercell configuration puts forth decreased reaction barrier than in the case of 2 × 2 supercell configuration. The MoS 2 /BP and MoSSe/BP demonstrates this reduction from 0.48 eV (Fig. 4 (e)) to 0.43 eV (Fig. 4 (g)) and 0.59 eV (Fig. 4 (f)) to 0 eV (Fig. 4 (h)), respectively. The reduced coverage and hence reduced charge redistribution on the surface can be attributed to the reduced reaction barrier in the 3×3 supercell. Further, we observe significant change in case of MoSSe and this may be attributed to the combined effect of the coverage and the electronegativity difference within the MoSSe layer that affects the charge transfer at the interface.
Apart from the factors that are discussed above, there are structural parameters that affect the reaction steps. optimize the H 2 molecule in the water layer. The reaction barrier, therefore, depends on the buckling in the monolayer, the water molecule's orientation, and the coexisting water molecules with H + (see Fig. 5).
Finally, we discuss the Heyrovsky reaction in MoS 2 /BP for 1/4 H + conc. both in case of 2 × 2 (i.e., 1 H + /4H 2 O) and 3 × 3 (i.e., 2 H + /8H 2 O) supercells. We observed reaction barrier decreases from 0.48 eV (Fig. 4 (e)) to 0.09 eV ( Fig. 6 (a)) and 0.59 eV (Fig. 4 (g)) to 0 eV ( Fig. 6 (b)) in 2 × 2 and 3 × 3 supercells, respectively. This indicates that high coverage prefers low H + conc. and vice versa for reduction in reaction barrier. We correlate this with the overpotential of the reaction, as discussed in the following section. Overpotential is the difference between the experimentally obtained reaction potential and the electrode potential. The electrode potential is analyzed only in the Heyrovsky reaction as it involves proton transfer. Therefore, this affects the work function and the potential at which the reaction takes place.  The electrode potential (U) of the slab is reported relative to the normal hydrogen electrode (NHE):

Electrode Potential
Here φ (E vac -E fermi ) is the work function, and φ NHE is taken to be 4.44 eV. 14,30,67 The work The potential drops are evident in Fig. 7, with a significant drop at the interface of BP and the water layer. The values of U corresponding to water layer with and without H + are reported in Table 1, which are in the range of -2.5 V to 1.3 V. We have incorporated dipole corrections as the vdW HTSs with two different surfaces maintain two potentials. Moreover, the H ads and, therefore, the coverage affects the dipole-dipole interactions. As a result, we report the two values of U, i.e., U 1 and U 2 , corresponding to two vacuum levels of E vac_1 and E vac_2 , respectively.
As previously discussed the dependence of φ on water orientation, we have explicitly opti-  Now we progress towards the extrapolation approach to cater the problem of potential change from initial to final in case of proton transfer Heyrovsky reaction. In this approach, we obtain reaction energies and activation energies of system with different supercell sizes and H + conc. The former represents the energy difference between initial and final states, while the latter is the amount of energy required to overcome the reaction barrier. Thereafter we obtain ∆E R and E a vs ∆U plot. The ∆U signifies change in electrode potential from initial to final. Moreover, the change in U 1 (corresponding to water layer potential) is significant as compared to the change in U 2 (corresponding to MoS 2 layer potential). The potential drop and charge transfer would accordingly affect the U 1 and U 2 . Hence, the reaction taking place at BP layer surface is crucial and we must consider U 1 for our analysis of electrode potential.
Therefore, the ∆U represented in the plot is corresponding to the U 1 . On extrapolating the ∆E R to ∆U = 0, we obtain -1.33 eV. The negative value indicates the spontaneity of the Heyrovsky reaction step. In case of positive ∆E R , the Heyrovsky reaction would have been the rate determining step. The corresponding E a is obtained at 0.11 eV. Hence, on comparing vdW HTS with monolayer, the synergistic effect of the two layers play role in affecting the overpotential and hence the reaction mechanism.

Conclusion
In summary, we have modelled a dynamically stable MoS 2 /BP and MoSSe/BP vdW HTSs.
They have been studied for HER by computational hydrogen electrode model. The optimized structure with the water layer showed a significant potential drop at the surface-water interface. The electrostatic potential is further affected by the proton solvated in water layer and the H ads constituting coverage over the surface.     The optical properties have also been calculated by the GW approach (Many Body Perturbation Theory (MBPT)). Fig. S3 shows the optical response using GW@BSE method that calculates the dielectric function. This is a complex function where the expression for interband process is the imaginary part thereby giving the absorption spectra. The real part (Re( )) is deduced from the Kramers-Kronig relation. We observe both the vdW HTSs