Widely available active sites on Ni2P for electrochemical hydrogen evolution - insights from first principles calculations

We present insights into the mechanism and the active site for hydrogen evolution on nickel phosphide (Ni 2 P). Ni 2 P was recently discovered to be a very active non-precious hydrogen evolution catalyst. Current literature attributes the activity of Ni 2 P to a particular site on the (0001) facet. In the present study, using Density Functional Theory (DFT) calculations, we show that several widely available low index crystal facets on Ni 2 P have better properties for a high catalytic activity. DFT calculations were used to identify moderately bonding nickel bridge sites and nickel hollow sites for hydrogen adsorption and to calculate barriers for the Tafel pathway. The investigated surfaces in this study were the (10 % 10), ( % 1 % 120), (11 % 20), (11 % 21) and (0001) facets of the hexagonal Ni 2 P crystal. In addition to the DFT results, we present experiments on Ni 2 P nanowires growing along the h 0001 i direction, which are shown as eﬃcient hydrogen evolution catalysts. The experimental results add these nanowires to a variety of diﬀerent morphologies of Ni 2 P, which are all active for HER.


Introduction
Water splitting by electrolysis or photocatalysis is attracting attention as a prospective sustainable source of hydrogen for energy storage applications. [1][2][3] The hydrogen evolution reaction (HER) is the cathode half reaction: H + (aq) -1/2H 2 (g) The most active and stable catalysts in acid environments for HER are currently based on precious metals. 4 At low overpotentials, HER from Pt(111) is dominated by the Volmer-Tafel mechanism 5-7 which can be written as: This requires two H atoms to be adsorbed in proximity for fast diffusion and recombination. Alternatively the reaction may take place via a Volmer step (1), followed by the Heyrovski step: H + (aq) + H* + e À -H 2 (g) Reaction (3) is considered to be relevant at high over-potentials only. 7 In the last decade, several classes of non-precious materials have been found to be active catalysts for the HER. MoS 2 was proven as a promising non-precious HER catalyst material, which is stable in a wide pH range. However MoS x catalysts are not as active as platinum and they are only active at undercoordinated edge sites. [8][9][10][11] Hydrodesulfurization (HDS) catalysts such as Mo 2 C and MoB have recently attracted attention as hydrogen evolution catalysts with good stability in both acid and alkaline solution. 12,13 Ni 2 P has also previously been subjected to experimental and theoretical studies for the catalysis of hydrodesulfurisation [13][14][15][16] water-gas-shift, 17 and more recently for hydrogen evolution in acid. [18][19][20] The active sites and the details of the mechanism remain unknown for several of the newly discovered HER catalysts including Ni 2 P. The observed Tafel slopes of the Ni 2 P catalysts are similar to what is observed for MoS 2 edges, and the specific activity is one of the highest of the non-precious catalysts. Further experimental and theoretical studies can yield new insights for further design of electrocatalysts, which is the motivation of this study. In this paper, density functional theory (DFT) calculations are combined with experiments on high aspect ratio crystalline Ni 2 P nanowires to understand the mechanism of HER from Ni 2 P.
The trends in hydrogen evolution activity over various transition metals, 21,22 and various metal-and non-metal combinations have been investigated using DFT. The free energy of hydrogen adsorption, DG H* , has been established as a descriptor for predicting the exchange current density, [21][22][23][24] of transition metal catalysts. The best catalysts have free binding energies close to 0.0 eV, which is explained by the Sabatier principle; 25 stronger binding results in hydrogen poisoning, leaving no free sites for adsorption. Weaker binding results in a high overpotential needed to adsorb protons. Transition metal catalysts have binding energies that are slightly coverage dependent. On Pt(111) the binding energy calculated with DFT varies from À0.12 eV at low coverage to 0.04 eV at one monolayer coverage. 21 In addition, it has been shown by DFT calculations that platinum has no significant diffusion barrier between the adsorbing hollow sites on the (111) facet. 26 Liu and Rodriguez have published several studies on DFT calculations of hydrogen binding energies on the Ni 2 P(0001) surface. In 2005, they predicted the hydrogen evolution activity of Ni 2 P to be comparable to that of hydrogenase based on trends in adsorption energies. The ensemble of phosphor atoms available as proton acceptors next to moderately binding metal hollow sites and a weak binding Ni-P bridge is proposed to facilitate HER catalysis. 18 The binding energies can be compared with the trends in HER exchange current density calculated by Nørskov et al. 21 In these studies, the binding free energies are presented, which means that the calculated potential energy is corrected by +0.24 eV accounting for entropy and differences in zero point energy. When comparing the binding energies at the relevant coverage, 7,27 one observes that the metal hollow sites on Ni 2 P will be fully occupied and that the Ni-P sites will need an overpotential of at least 0.31 V. This does not agree with the very high activity observed in new studies of hollow and multifaceted Ni 2 P nanoparticle catalysts. 19,20,28,29 Another active site may therefore exist, and we investigate this using more detailed DFT calculations. In addition we report experiments showing that nano wires grown along the (0001) direction are highly active for HER, supporting the results of the present DFT calculations.

Calculation methods
Calculations were carried out using the GPAW code with projector augmented wave functions on a real space grid and ultra soft pseudopotentials. 30 The RPBE 31 functional was used for the exchange-correlation contribution. The Ni 2 P crystal structure 32 was copied from The Materials Project, 33 and imported to The Atomic Simulation Environment (ASE). 34 Bulk lattice constants were converged to a = 5.151 Å and c = 3.408 Å using a third order polynomial fit to the energy versus the lattice constant. 35 A 1 Â 1 supercell was used for the (0001) and (11% 21) surfaces, a 3 Â 1 supercell was used for the (10% 10) surface and a 1 Â 2 supercell was used for the AE(11% 20) surfaces. The (0001) surface had five atomic layers, the (11% 21) and (10% 10) slabs had four atomic layers, and the AE(11% 20) surfaces had three atomic layers. The dipoles across the unit cells were all less than 0.10 V.
The slabs were relaxed using the Broyden-Fletcher-Goldfarb-Shanno LineSearch algorithm within ASE until the forces were below 0.01 eV Å À1 . A recalculation was later carried out with double the k-point sampling and a grid spacing of 0.12 Å for the most interesting slabs and adsorbate configurations. The resulting differences in adsorption energies did not exceed 0.02 eV, which is below the accuracy usually attributed to DFT.
The choice of facets and surface termination were based on calculations of the minimum energy configuration of all different ways to cut the crystal into the lowest index planes. Several STM and LEED studies [37][38][39][40] show that the Ni 2 P surfaces can have a stable phosphor termination, but under the hydro-thermal treatment before testing, it is expected that they lose the phosphor layer and expose active metal sites, as on the structures investigated in our calculations. The chosen structures are (10% 10), (11% 20), (% 1% 120), (11% 21) and the Ni 3 P terminated (0001), as shown in Fig. 1 along with their (hkjl) indices from a top view.
The free adsorption energies are found from calculated potential energies by correcting for the gas phase entropy DS = ÀS 0 (H 2 ) and the difference in zero point energy DZPE using the equation where E denotes the ground state energy at 0 K obtained from DFT. The correction for ZPE and entropy makes an addition of 0.24 eV to the adsorption energy of a hydrogen atom. 21,41 The calculation of adsorption free energies and adsorbate coverages was carried out using the self consistent scheme as in the work reported by Skulason et al. 7 The integral adsorption energy G int (n) is where G(N,n) is the free energy of a surface which includes N nickel atoms in the top layer and n adsorbed hydrogen atoms. Thus we use a definition of the coverage y = n/N, where n is the number of hydrogen adsorbates and N is the number of Ni atoms on the surface. The configuration which is relevant at a given chemical potential of protons and electrons is the one with the minimum integral adsorption energy (G int (n)). Using the computational hydrogen electrode approach 41 with self-consistent coverage and adsorption free energies, 27 the adsorption phase diagram was calculated for every facet. The required over-potential to favor an intermediate step is given through the differential binding energy, which is given by: We apply Z = DG diff (n)/e to find the minimum required overpotential to adsorb the intermediate, 41 which enables a surface recombination reaction (the Tafel step). The barriers were calculated using the nudged elastic band (NEB) method with a climbing image. 42

Synthesis of nickel phosphide nanowires
The synthesis was performed following a previous report. 43 A stock solution containing 0.75 mmol of Ni(acac) 2 , 1.8 mmol of oleic acid and 10 mL of OA was heated at 120 1C. Then the stock solution was very slowly injected (0.05 mL min À1 ) using a syringe pump into a stirred mixture containing 5 mL of OA and 2.4 mmol of TOP heated at 320 1C under reflux. The reaction mixture was kept at 320 1C under reflux until the stock solution was used up. Over the course of the reaction, the mixture color changes from transparent to dark yellow, orange and finally black. After cooling to room temperature, the products were washed using a mixture of hexane and ethanol, and separated by centrifugation. The supernatant was discarded and the washing step was repeated two more times. After the supernatant was discarded the nanowires were collected on a watch glass and dried in an oven at 50 1C overnight.
Physical characterization X-ray diffraction (XRD) was carried out on a X'Pert Philips diffractometer in Bragg-Brentano geometry with Cu K a1 radiation (l = 0.1540 nm) and a fast Si-PIN multi-strip detector. The tube source was operated at 45 kV and 40 mA. The Scherrer equation was used to calculate the average crystallite size of the powder.
The diffraction pattern was analyzed and compared with references in the international center of diffraction data (ICDD). Transmission electron microscope (TEM) images were taken on a Philips (FEI) CM12 with a LaB 6 source operated at 120 kV accelerating voltage. Specimens were prepared by ultrasonic dispersion of the samples in ethanol. The suspension was mixed with a micropipette by several suction-release cycles to ensure a representative and reproducible TEM sample is obtained. A few drops of the mixed suspension were deposited onto the carbon-coated grid.

Electrochemical measurements
Electrochemical measurements were recorded using a Gamry Instruments Reference 3000t potentiostat. A traditional threeelectrode configuration was used. For polarization and electrolysis measurements, a platinum wire was used as the auxiliary electrode and a double-junction Ag/AgCl (KCl saturated) electrode was used as the reference electrode. Both counter and reference electrodes were rinsed with distilled water and dried with compressed air prior to measurements. Potentials were referenced to a reversible hydrogen electrode (RHE) by adding a value of (0.2 + 0.059 Â pH) V. The current was normalised over the geometric surface area of the electrode. Ohmic drop was corrected using the current interrupt method. A total electrolyte volume of B50 mL was used to fill the glass cell. All potentials were converted and referred to the RHE unless stated otherwise. The electrolyte used throughout all electrochemical experiments was a 1 M H 2 SO 4 solution. During electrochemical experiments, the electrolyte was agitated using a magnetic stirrer rotating at 300 rpm. For all electrochemical experiments a glassy carbon electrode (B0.071 cm 2 ) was used as the working electrode. The cyclic voltammetry (CV) experiments were conducted in 1 M H 2 SO 4 at 25 1C using a scan-rate of 5 mV s À1 across a potential window of À0.3 to +0.1 V vs. RHE.

Pretreatment of the working electrode
Prior to loading of the catalyst, the working electrode was pretreated to achieve a better performance. 44 First, the electrode was manually polished using alpha alumina powder (CH Instruments, Inc.) with decreasing grain sizes (typically 0.3 and 0.05 mm) on a 73 mm diameter nylon polishing pad (MasterTex, Buehler). Between each polishing step the electrode was rinsed with deionized water and ultrasonicated in distilled water for 10 seconds. Then the electrode was dried using compressed air. This ultrasonication and drying cycle was repeated two more times: once in absolute ethanol and a second time with distilled water. This polishing process resulted in a shiny mirror finish. The bare working electrode was subjected to a constant potential of 2 V vs. RHE in 1 M H 2 SO 4 under vigorous stirring at 25 1C over 1 hour. Then, the pretreated electrode was rinsed with absolute ethanol and dried with compressed air prior to catalyst deposition.

Electrode preparation
Before being deposited on the working electrode, the nanowires were annealed at 450 1C for 4 hours under 5% H 2 /N 2 gas in order to remove any surfactant present on the nanowire surface. 19,43 10826 | Phys. Chem. Chem. Phys., 2015, 17, 10823--10829 This journal is © the Owner Societies 2015 The catalyst was loaded on a pretreated working electrode via drop-casting of 10 mL of the catalyst ink, equivalent to a loading of 1.42 mg cm À2 . The catalyst suspension was a 500 mL solution consisting of 5 mg of the catalyst, 400 mL of distilled water, and 100 mL of absolute ethanol. The slurry was ultrasonicated for 5 hours and mixed with a micropipette by several suctionrelease cycles prior to deposition to ensure that a representative and reproducible catalyst sample is obtained. The temperature of the ultrasonic bath was kept below 45 1C at all times to avoid any undesired heating effect. The nanowires are ultra-sonicated for 5 h to minimize aggregation. Once the catalyst was deposited, the electrode was dried in an oven at 50 1C for 10 minutes. 5 mL of 0.2% Nafion s solution in absolute ethanol was then dropcasted on the glassy carbon surface and the electrode was dried in air for 5 minutes prior to electrochemical measurements.

Results and discussion
XRD measurements were conducted on the annealed Ni 2 P nanowires (Fig. 2). According to Scherrer's equation, the average crystallite size of the nanowires is 11 nm. Ni 2 P is the predominant species in the sample, while small amounts of Ni 12 P 5 and NiO are also present. Recently, Zhang et al. reported that pure Ni 12 P 5 has a similar activity to that of Ni 2 P for HER in acid. 45 We showed earlier that NiO is less active than Ni 2 P; furthermore, NiO tends to dissolve at a reductive potential in acid. Thus, the HER activity of the nanowires (see below) can be mostly attributed to Ni 2 P.
The transmission electron microscope (TEM) image and selected area electron diffraction (SAED) patterns of the nanowires are shown in Fig. 2a-d. The nanowires are rather uniform and the width of the nanowires observed by TEM is in agreement with the crystallite size calculated from the XRD pattern.
The focused SAED pattern (Fig. 2b) was indexed using the JEMS software program (P. Stadelmann, JEMS, EPFL). The line linking the (000n) spots is the (0001) direction. The nanowire can be observed on the SAED pattern when the focus is reduced ( Fig. 2c and d). Comparison between the indexed (Fig. 2b) and unfocused SAED patterns ( Fig. 2c and d) allows the determination of the growth direction. The results indicate that the nanowires grow along the c-axis, i.e., the h0001i direction.
Thus, the nanowires do not expose many (0001) facets, since the (000n) planes exist only in their cross section (see Fig. 3). The low index facets AE(11% 20) and (10% 10) are assumed to be widely available on the samples. These morphologies would be expected to exhibit a less than optimal performance if the ensemble on the Ni 3 P terminated (0001) facet was the only active site. Fig. 4 shows the HER activity of the Ni 2 P nanowires in 1 M H 2 SO 4 . The nanowires are excellent HER catalysts. This result combined with the calculations (see below) does not indicate that the (0001) ensemble is the only active site of HER from Ni 2 P, although it is not possible to separate contributions from the various facets and sites in the experiments. The overpotential to drive a current density of 10 mA cm À2 is 133 mV (Fig. 4a). Two Tafel slopes are observed. At Z o 125 mV, the Tafel slope is about 60 mV dec À1 , while at Z E 125-275 mV, the Tafel slope is 126 mV dec À1 (Fig. 4b). Despite the difference in sample preparation and morphology, the catalytic activity of the Ni 2 P nanowires is very similar to the activity of hollow and multifaceted Ni 2 P nanoparticles reported by Schaak et al., 19 that of polydispersed Ni 2 P nanoparticles reported by our group 20 and that of high surface area Ni 2 P nanoflakes reported by Han et al. 29 This suggests that widely available facets or sites of Ni 2 P are active. This is contrary to MoS x , where only specific morphologies are efficient because the edge sites need to be exposed. 8 The facets with (ab À(a + b) 0) indices are propagating along the growth direction of the nanowires (see Fig. 3), and are thus expected to be abundant in the sample.
The DFT calculations show that the most strongly binding sites tend to be nickel bridge sites or nickel 3-fold hollow sites. The nickel bridge sites are found in continuous rows on the (10% 10) and (% 1% 120) facets and adjacent to nickel hollow sites on the (11% 21) facet. (see ESI, † for geometries and adsorption  energies.) Adsorption on Ni-P sites and on top P atoms was generally 0.2 eV to 0.4 eV less stable than the metal sites. Thus, adsorption on Ni-P sites and on top P atoms would require a higher over-potential.
It takes a high over-potential to favor a high coverage of atomic hydrogen on all the facets. Thus, Ni 2 P is on the weak binding side of the Sabatier volcano. Adsorbates on Ni 2 P interact strongly leading to a steeper increase in binding energy with increasing coverage (see ESI, † for the surface phase diagrams), compared to most transition metals including platinum.
The (0001) facet has the hollow site occupied at equilibrium, but as shown in Fig. 6, it requires a large over-potential of 0.41 V to stabilize the initial Ni-P bridge state for the Tafel step. It is also improbable that adsorbates from separate hollow sites on (0001) can combine since they are not neighbouring. This was confirmed by calculating the diffusion barriers on the (0001) and the (10% 10) facets (see ESI †). The (10% 10), (% 1% 120) and (11% 21) facets only require an over-potential of 0.0 V, 0.06 V and 0.19 V, respectively, to stabilise neighbouring H* or H*, which are mobile in the row of nickel atoms. These results suggest that H* in Ni hollow sites on the (0001) facet do not contribute to the HER rate at low over-potentials. To investigate further, the barriers for the Tafel step were calculated as described in the following.
Previous calculations on Pt(111) suggest that the Tafel pathway is faster than the Heyrovski pathway at low over-potentials. If proton transfer is faster to H*(Ni 2 P) than to H*(Pt(111)), it is possible that isolated metal hollow sites on the (0001) facet play a role. In the following, it is assumed that the surface coverage is in equilibrium with protons in solution, which means that the Tafel step is rate limiting. Calculation of the barrier for proton transfer (Volmer or Heyrovski step) is out of the scope of this paper, since it requires very precise information on the interfacial structure. 46 As shown in Fig. 5, the rate limiting HER barrier is G TS , where G TS is the energy of the transition state relative to the 2(H + + e À ) state. The initial states were the most stable configurations at the lowest possible over-potential according to the surface phase diagrams. As shown from calculations presented in Fig. 5, G TS can be lowered by further increasing the over-potential until the free energy of protons in solution are at the Tafel transition state level. This agrees well with the exchange current density being a good indication for the activity at higher overpotentials.
The results for a Tafel pathway are summarized in Fig. 6, with comparison between the studied facets. The adsorption energies, which are easy to calculate, are usually a good descriptor, since they are expected to scale with the transition state energies. In the case of a Tafel pathway, it is more accurate to use calculated G TS to compare the activity from the different facets, since it is the highest barrier in the reaction pathway, which limits the rate.
Observing the calculated G TS in Fig. 6 for the Tafel steps, it is clear that the (% 1% 120) facet and the (11% 21) facet should have the highest exchange current density and thus be the most active facets. The reason why these sites have the lowest combination barrier may be found in the geometry of the Ni atom binding  10) facet. The configuration of every point is shown by the insets, where Ni is shown in green, P in yellow and H in white. The energy of the transition states is shown by the peak of the splines. These were calculated using the nudged elastic band method with a climbing image. 41 This shows how G TS is reduced by applying a higher over-potential. Fig. 6 Free energy diagram assuming a Tafel mechanism. Calculated barriers for the Tafel step are shown by the peak of the splines for the configurations, which are most stable at 0 V vs. RHE. The insets show transition states of the two facets with the lowest Tafel barriers. Ni atoms are shown in green, P in yellow and H in white.