Neutron spectroscopy study of the diffusivity of hydrogen in MoS₂

Abstract

The diffusion of hydrogen adsorbed inside layered MoS₂ crystals has been studied by means of quasi-elastic neutron scattering, neutron spin-echo spectroscopy, nuclear reaction analysis, and X-ray photoelectron spectroscopy. The neutron time-of-flight and neutron spin-echo measurements demonstrate fast diffusion of hydrogen molecules parallel to the basal planes of the two dimensional crystal planes. At room temperature and above, this intra-layer diffusion is of a similar speed to the surface diffusion that has been observed in earlier studies for hydrogen atoms on Pt surfaces. A significantly slower hydrogen diffusion was observed perpendicular to the basal planes using nuclear reaction analysis.

---

1 Introduction

Diffusion of molecules on surfaces is crucial for a range of chemical processes, such as catalysis or sensing, but it remains extremely difficult to observe experimentally on atomic length scales. This is particularly true for the case of light molecules that have few electrons, such as hydrogen. Substantial breakthroughs in this field have been achieved by quasi-elastic helium atom scattering (QHAS) and quasi-elastic neutron scattering (QENS) in recent years.^1–3 In the following we present a study that exploits QENS in this tradition.

With an increase in renewable energy production, the problem of energy storage becomes more and more significant. One of the most promising ideas is storage in form of hydrogen gas. It involves the production of hydrogen by electrolysis and later reconversion of hydrogen in fuel cells. It is of paramount importance for these processes to be carried out most efficiently with the most readily accessible materials. Therefore, there have been intense studies in the search for catalyst materials for the hydrogen and oxygen evolution reactions,^4,5 in order to replace the expensive platinum, which is widely used today.

Molybdenum disulphide, MoS₂, has shown promising behaviour as a catalyst in the hydrogen evolution reaction (HER),^6–10 which is completely in line with its known activity for hydrogenation reactions. However, the structural properties are rather complex. MoS₂ forms a crystal of hexagonal Mo and S layers, with an interlayer distance of 6.2 Å.¹¹ On the basis of this structure, three stacking variants have been observed, each with slightly different conductivity and catalytic activity.^6,12 In earlier studies catalytic activity had been attributed to amorphous particles,¹³ but recently it has been shown that although the active centres are probably sulphur edge atoms, the stoichiometry of the active surface is 1Mo:2S, the same as for the perfect crystal.^6,7,10,14

However, the amount of active centres may not be the only reason for chemical activity. For palladium, for example, during hydrogen evolution also the movement of hydrogen atoms to sites below the surface has been observed (HIR, hydrogen injection reaction).¹⁵ For platinum surfaces, the recombinative desorption of hydrogen was shown to be the rate limiting step in HER.¹⁶ For MoS₂ some theoretical studies suggest that intercalation of hydrogen is possible. Therefore, the diffusion of hydrogen can be expected to play a role in the overall reaction mechanism.^17,18

Within the scope of the present study we investigated the adsorption and diffusion of hydrogen on MoS₂ using quasi-elastic neutron scattering (QENS) and neutron spin-echo spectroscopy (NSE) together with nuclear reaction analysis (NRA) and X-ray photoelectron spectroscopy (XPS).

This paper is organised as follows: after this short introduction, we will give a short description of the theoretical model that we use to analyse the spectroscopic neutron scattering data. Then, we will describe the sample preparation and the experimental apparatus that we used. Finally we will describe our experimental findings and discuss them in the framework of earlier studies.

2 Theory

A general equation describing the incoherent neutron scattering off the adsorbed hydrogen consists of four parts:² the Debye–Waller factor, , the elastic scattering function, S_el(Q,ω), the quasi-elastic scattering function, S_qe(Q,ω), and a resolution function R(Q,ω):

S(Q,ω) = S_DW(Q)δ(ω) * [S_el(Q)δ(ω) + (1 − S_el(Q))S_qe(Q,ω)] × R(Q,ω). (1)

Q = k_i − k_f is the momentum transfer and ħω = E_f − E_i is the energy transfer of the neutron upon scattering with k_i and k_f being the initial and final momenta, and E_i and E_f being the initial and final energies of the neutron, respectively. ħ is the reduced Planck constant. The symbol * describes a convolution action. The term, which contains the most important information about the diffusive dynamics in the system, is S_qe(Q,ω). We will consider a model developed by Sears based on the first Born approximation for the case of a moving polyatomic molecule exhibiting both translational and rotational motion with negligible correlations between them.^19,20 The model has been applied successfully to a wide range of molecules in confinement, including hydrogen, water, benzene and neopentane.^21–24 For undercooled water the approximation of negligible correlations was tested by molecular dynamics simulations and shown to be applicable at temperatures where molecules are free to move and rotate.²⁵ For hydrogen molecules in MoS₂, due to the low hydrogen density, the translational–rotational decoupling can be assumed to be valid for diffusion within the layers of MoS₂. For the diffusion perpendicular to the layers, however, this approximation might be difficult to apply due to the spacial confinement. Since we are mainly interested in the hydrogen motion within the layers, however, we consider that the approximation of negligible correlations is reasonable. In this case S_qe(Q,ω) can be represented by a convolution of translational, S_T(Q,ω), and rotational, S_R(Q,ω), components, respectively. The rotational component is described by:¹⁹where j_l(x) is a spherical Bessel function of the first kind of order l, δ(ω) is the Dirac delta function, R is the radius of gyration of a molecule, is a normalised Lorentzian function and is the half-width at half-maximum (HWHM) of the corresponding rotational contribution with τ_rot being the characteristic time of rotation. We note that for describing the motion of hydrogen atoms, the rotational terms disappear (S_R(Q,ω) ≡ δ(ω)).

Depending on a system, translational motion can either have a free or a spatially constrained behaviour. In the former case it can be expressed as a simple normalised Lorentzian . In the latter case a part of the quasi-elastic scattering remains elastic and a delta-function term, A₀(Q)δ(ω), is added:²

S_T(Q,ω) = A₀(Q)δ(ω) + (1 − A₀(Q))L(Γ_T,ω), (3)

where A₀(Q) is called an elastic incoherent structure factor (EISF).²⁶ For confinement within a sphere of radius r_c it can be written as:²⁷In the following we will use the confined diffusion model since hydrogen motion in a layered MoS₂ is naturally constrained by the MoS₂ matrix. S_qe can then be expressed as:The convolution in eqn (5) can be performed explicitly, giving:

A quasi-elastic broadening caused by translational diffusion Γ_T depends heavily on the type of motion, which can be continuous or jump-like. A number of different models have been considered that describe each kind of diffusion as usually several of them need to be tested on a system in order to get a reliable result. Three models will be discussed in detail later in the Results and discussion section.

Here, we will consider a commonly occurring case, Γ^l_R ≪ Γ_T, which allows us to neglect Γ^l_R in the last series in (6). Thus, the fourth item becomes . The summation here starts from l = 1, however one can note that the third item [1 − A₀(Q)]j₀²(QR)L(Γ_T,ω) is actually a term of this sum with l = 0. Hence, withdrawing Γ_T from the sum as it does not depend on l, one can get:

Within the scope of our experiment (R = 0.37 Å for the hydrogen molecule, Q < 2 Å⁻¹) we generally have QR < 0.74 and in the first series we can neglect all terms higher than the 3rd order, as . The last sum equals unity as for any value of the variable x.²⁸ Finally, we obtain:

By substituting eqn (8) into eqn (1), adding a background term and regrouping elements we end up with the ultimate fitting expression for our quasi-elastic spectra:

Here, S₀(Q) = I₀·S_DW(Q) is an overall intensity factor, which accounts for signal loss/gain due to vibrations (Debye–Waller factor S_DW(Q)) and due to hydrogen desorption/adsorption (scaling factor I₀).

3 Experimental

3.1 Sample preparation

Molybdenum disulphide single crystals with a thickness of about 0.2 mm and edge lengths between 6 and 15 mm were obtained from PLANO (PLANO GmbH Wetzlar, Germany, distributor of SPI, USA). Hydrogen loading via electrolysis was performed several times, prior to every experiment. For nuclear reaction analysis (NRA) one MoS₂ single crystal was loaded with hydrogen shortly before the experiments, and for the XPS study four days in advance. Neutron studies require a larger volume of material and therefore all 10 single crystals were used. They were hydrogen loaded 4–6 days prior to the experiments. The loading procedure was the same for every experiment with slight variations in the duration of the process. Those samples that had been used in experiments before were degassed at temperatures >750 K. Prior to each loading all MoS₂ single crystals were cleaved by peeling off the surface layer using adhesive tape. Then the crystals were mounted onto a tantalum plate electrode and subsequently loaded by electrolysis in a 1 : 1 mixture of phosphoric acid (H₃PO₄ 85%) and glycerol, using a Pt counter electrode. The electrolysis lasted 2 hours at 4 V on one side of the crystal for the nuclear reaction analysis, 4 hours at 4 V on one side for XPS and 4 hours at 3 V on each side for neutron experiments.

The structure of MoS₂ shown in Fig. 5 and 12 was plotted using the VESTA 3 software.²⁹

3.2 Nuclear reaction analysis

Hydrogen depth profiles were measured by irradiating the sample with a ¹⁵N ion beam and counting the emitted gamma-rays on the NRRA beamline at RUBION, Bochum, Germany. The underlying reaction is:

¹⁵N + ¹H → ¹²C + α + γ,

with a narrow resonance at the energy of 6.4 MeV, which is used for the detection of hydrogen. In order to reach the resonance energy at different depths inside the sample, the beam energy was varied from 6.39 MeV to 6.65 MeV. The maximum beam current was set to 6 nA. Data were collected in two runs, shortly after loading the sample with hydrogen and after four weeks downtime. During the waiting time the sample was held in air under ambient conditions. The depth scale in nm was calculated on a basis of the stopping power of MoS₂.^30,31 The concentration scale was constructed taking into account the stopping power of both MoS₂ and hydrogen. More detailed information about the instrument and the analysis can be found elsewhere.^32,33

3.3 X-ray photoelectron spectroscopy

XPS surface characterisation was carried out on the HE-SGM station³⁴ at the BESSY II synchrotron, Helmholtz-Zentrum Berlin (HZB), Germany. The XPS signal was collected with a Scienta R3000 XPS analyser. The pass energy was set to 50 eV in order to decrease the background. Excitation energies used were 650 eV for the oxygen 1s line, and 385 eV for the other elements. The energy binning step was set to 50 meV. Data were collected posterior to the sample being heated to several temperatures: 300, 400, 450 and 500 K, respectively. During heating the pressure in the analysis chamber rose due to desorption from the sample. Therefore we had to wait several minutes until the pressure dropped to a value appropriate for the XPS experiment. The measurements thus show the state of the surface several minutes after the sample had been heated up to a certain temperature. For energy calibration the adventitious carbon peak at the corresponding temperature was used. The background was described by a Shirley step model.³⁵ The intensities were corrected for the atomic subshell photoionisation cross sections obtained by Yeh and Lindau.³⁶ Data treatment was performed using the XPST package within the Igor Pro software package (WaveMetrics, Portland, USA).

3.4 Quasi-elastic neutron scattering experiments

QENS data were collected at the neutron time-of-flight (TOF) spectrometer TOFTOF³⁷ at the research neutron source Heinz Mayer-Leibniz (FRM II), Garching, Germany, at a wavelength λ = 5 Å. The samples were held in a flat aluminium sample cell. Two sample orientations, an in-plane and an out-of-plane set-up, were used in order to obtain information about dynamics parallel and perpendicular to the MoS₂(0001) basal plane, respectively. For both set-ups the surface normal (crystallographic c-axis) was within the scattering plane. For the in-plane set-up, the c-axis was turned by ϕ = 20° away from the direct beam direction, i.e., Q was within the basal plane exactly for a scattering angle 2θ = 40° (Fig. 1). For the out-of-plane set-up we used ϕ = 110°, i.e., for 2θ = 40° the momentum transfer was perpendicular to the basal planes. The temperatures scanned in the out-of-plane set-up were 3, 100 and 300 K, whereas in the in-plane set-up data at 3, 50, 100, 200, 300, 400, and 500 K were measured. The spectra were measured consecutively with rising temperature up to 500 K and a further spectrum at 300 K in the in-plane set-up was measured after the last 500 K measurements to check for permanent changes of the sample upon heating to 500 K. Spectra at 3 K were used as a resolution function. Data reduction was performed using the Mantid software package.³⁸ Every spectrum was normalised to the beam monitor signal and, after deduction of an empty container spectrum, corrected by a vanadium spectrum to calibrate for variations in detector efficiency. Both empty container and vanadium spectra were measured at room temperature. All data were binned in energy transfer space from −10 to 3 meV with a step of 0.05 meV and in momentum transfer space from 0.05 to 2.35 Å⁻¹ with a step of 0.1 Å⁻¹. A representative QENS scan of MoS₂ at 400 K in the in-plane sample geometry is shown in Fig. 2. No Bragg reflections or significant data reduction artefacts are observed. Further in the text the geometry specification is dropped out, by default referring to the in-plane set-up. Every reference to the out-of-plane set-up is indicated explicitly.

Fig. 1 A scheme of the TOFTOF experiment. The incident beam k_i is shown as a horizontal line, c is a crystallographic c-axis of MoS₂ which is perpendicular to both sample surface and flat aluminium sample holder, k_f is a scattering vector for which the momentum transfer vector Q is strictly parallel to the MoS₂ basal planes. Here the so-called in-plane set-up with ϕ = 20° is depicted. For the out-of-plane geometry the sample was rotated 90° counter-clockwise, i.e. ϕ = 110°.

Fig. 2 Scattering function, S(Q,ω), of hydrogen loaded MoS₂ single crystals, measured at 400 K with the sample's basal planes oriented parallel to the momentum transfer vector Q. The plot uses a logarithmic intensity scale. Higher intensities are represented by red colour and lower intensities by blue. White area: no data, since it falls outside of the spectroscopic window.

3.5 Neutron spin echo measurement

Neutron spin echo (NSE) measurements were performed on the NSE spectrometer IN11 at Institut Laue-Langevin (ILL), Grenoble, France,³⁹ using the high signal set-up IN11C, which uses a 30° detector bank.⁴⁰ Here, we used a wavelength of λ = 5.5 Å for maximum signal. Following the same procedure as in the QENS study, data were collected in two sample orientations with the surface parallel and perpendicular to the scattering vector Q. The spectra were normalised using a standard procedure: all data were divided by a spectrum that was obtained in situ at 2 and 5 K for the in-plane and out-of-plane geometry, respectively. At these temperatures the system can safely be assumed to be static in the dynamic window of the spectrometer. The experiment covered a range of momentum transfers of Q = [0.1 to 0.7] Å⁻¹ and sample temperatures up to 600 K.

4 Results and discussion

4.1 Nuclear reaction analysis

We commence the description of our experiments with the structural and chemical characterisation of the samples: Fig. 3 shows the hydrogen concentration in dependence of the depth measured by NRA quickly after electrolysis loading of MoS₂ and after four weeks of ageing, respectively. For the freshly prepared sample the data clearly show a high concentration of hydrogen in the surface and sub-surface region and vanishing signal at higher depths (we note that the width of the surface peak is equal to the width dictated by the energy resolution of the incoming beam and, hence, the width of the surface peak gives only an upper limit for the width of the distribution function at the surface). After 4 weeks of ageing in ambient air environment, the surface peak decreased in intensity by half and a significant bulk signal was recorded that remained constant from a depth of 20 nm up to 100 nm, the maximum penetration depth recorded in our NRA measurements. Hence, the data show a diminution of hydrogen containing species in the surface layer and an increase of the concentration of hydrogen containing species in the bulk. After this first straightforward result, however, there remain two important further questions: (a) what hydrogen species do we have in this experiment and (b) can we say anything about the diffusivity?

Fig. 3 Hydrogen depth profile in H loaded MoS₂ single crystals, measured by NRA shortly after loading (red curve) and after four weeks of ageing (black curve). The concentration scale applied to the ordinate does not apply to the surface peak.

Considering the first question, one can imagine three hydrogen containing species that can be present in significant concentrations after the electrolytical loading: hydrogen atoms, hydrogen molecules and water. To distinguish between them, we have performed a test experiment with a sample as received and a sample that was exposed to pure water without electrolysis. In these cases, we also obtained a surface peak at the start of the experiment, but no diffusion into the bulk over time. Therefore, we can conclude that water may contribute significantly to the surface signal, but does not penetrate across the basal planes into the bulk in a significant amount. The hydrogen containing species that are found in the bulk of the sample after electrochemical loading are, hence, hydrogen atoms or hydrogen molecules. The concentration scale, in general, is not fully applicable to the surface peak due to several reasons, namely, hydrogen multilayer formation, surface contaminations and failure of the thick target approximation.³³ The first is not expected in the current experiment, however, light organic contaminations and only partial validity of the thick target approximation renders the concentration scale, to some extent, overestimated. Nevertheless, it is still possible to extract some information from the decay of the surface peak over time: We have demonstrated above, that only hydrogen molecules or hydrogen atoms diffuse into the bulk. This migration must lead to a loss in intensity of the hydrogen surface peak. Simultaneously, a gain in intensity may be observed due to a possible water adsorption from the atmosphere at ambient air conditions. Since during the ageing process the surface peak area decreased approximately by half, we conclude that several hours after loading process at least half of surface and subsurface hydrogen is in form of atomic and molecular hydrogen.

Now, we will try to analyse the diffusion in some more detail. The above considerations about the limited usability of the concentration scale in application to the surface peak concern mostly the depth profile up to ≈15 nm. Therefore the data beyond 15 nm can be used to extract some quantitative results. Fitting the data for the freshly loaded sample to the solution of Fick's second law of diffusion in the form of a Gaussian yields an order of magnitude for the diffusion coefficient D_⊥ of 10⁻²¹ m² s⁻¹. Thus, the observed diffusion is extremely slow. However, even such a slow diffusion would result in almost uniform hydrogen concentration across the bulk after a 4 weeks ageing period, which is clearly observed in the data.

4.2 X-ray photoelectron spectroscopy

Prior to the neutron spectroscopy study of hydrogen diffusivity, see below, the chemical state of the MoS₂ surface was monitored by X-ray photoelectron spectroscopy at temperatures up to 500 K (Fig. 4). The XPS spectra reveal dramatic changes in the electronic structure of the MoS₂ surface upon heating. The splitting of the molybdenum and sulphur lines into two peaks each with almost identical areas was observed starting from 400 K up to the highest measured temperature 500 K, indicating the separation of the surface into oxidised and slightly reduced sections. These sections may be assigned to the surface and subsurface regions respectively. We note, that the analysis depths of the XPS experiment is limited to less than 10 nm.⁴¹ In view of the results of the NRA study, this covers the region, where we expect to find water and hydrogen. The analysis of the Mo and S binding energies suggests that at room temperature the surface consists mostly of clean MoS₂, i.e. Mo⁴⁺ and S²⁻. Heating to 400 K partially reduces the surface, with Mo staying closer to the 4+ state than to 3+. Further increase in temperature prompts a new oxidation process, which results in a significant increase in the binding energy of the more oxidised state (labelled 2 in Fig. 4), causing Mo⁶⁺ to be formed at 500 K. The assignment of oxidation states remains tentative, however, since charging of the modified part of the surface may also have an impact on the observed binding energy. Keeping this reservation in mind we describe the occurring process as follows. During the loading procedure, apart from the production of –SH groups on the surface, a number of hydrogen atoms are incorporated into the bulk. At 400 K they become mobile and are subsequently bound also in the subsurface region, thus performing a reduction. At 450 K hydrogen atoms start to desorb from the surface, which leads to a surface oxidation. At the same time, as the H coverage on the surface decreases, subsurface hydrogen starts to migrate further to the surface, slowing down the surface oxidation and simultaneously oxidising the subsurface region. At 500 K the same process continues. It is also important to note that since MoS₂ is considered relatively active for the hydrogen evolution reaction, another source of subsurface H atoms apart from the bulk H exists, namely, catalytically split H₂ molecules. However in the XPS experiment little information can be drawn for hydrogen molecules, since they are not chemisorbed to the host material.

Fig. 4 XPS scans of various elemental lines for MoS₂ single crystals at different temperatures. From left to right: sulphur 2p, molybdenum 3d, carbon 1s and oxygen 1s. The temperature colouring is the same for all four plots: blue for 300 K, green for 400 K, red for 450 K and orange for 500 K. The grey dotted lines show the evolution of the binding energies of the two main oxidation states, labelled 1 and 2. For visual clarity a vertical offset is applied to the S and Mo data at elevated temperatures.

The XPS data contain also carbon and oxygen lines. Both elements are likely to be adsorbed from the atmosphere in various forms. Prior measurements utilising deuterium induced gamma emission (DIGE) at RUBION, Bochum, Germany, showed a bulk oxygen content at a level of 0.5 at%. Hence, oxygen should mostly be located in the surface or subsurface region. The interpretation of the O 1s line requires a careful analysis: molybdenum oxides are expected to create an XPS signal around 530 eV^42,43 and, thus, are not observed in this study. The other oxygen containing species that are expected in our samples are adsorbed water and organic molecules containing oxygen also in the form of –CO and –CO₂. The binding energies of oxygen in these compounds are very similar and, therefore, cannot be resolved within the O 1s line. However, since the area under the peaks is proportional to the amount of corresponding C and O atoms in a sample, in first order approximation, some information can be drawn from the behaviour of the C and O lines under temperature treatment. If every oxygen atom was bound to carbon, then the desorption picture, i.e. the loss of area under a peak, would be the same for both the O 1s line and the peaks corresponding to –CO and –CO₂. In our case no simple correlation can be established. Therefore it can be concluded that a significant amount of the observed oxygen is bound in water. As we have shown in the NRA study, however, this water has to be in the surface or subsurface regions. Additional XPS measurements of two MoS₂ single crystals loaded with hydrogen under the same conditions and heated to 500 K ex situ showed a great difference between the areas of the corresponding O 1s lines. Given a very small measuring spot size, the most probable cause of this effect is water agglomeration within volume defects. In this case, the amount of oxygen seen in the XPS experiment would depend on whether the incoming X-ray beams hit the water agglomerate or not. The O 1s line splitting that appears at 500 K indicates that a large fraction of the water molecules transfer from a strongly adsorbed state to almost free H₂O molecules at this temperature.^44,45

To conclude, the XPS study shows that chemisorbed hydrogen shows a thermally activated mobility in MoS₂ when the temperature of the sample is raised to 400 K or above. This leads to a chemical transformation of the sample as reduction and oxidation of MoS₂ can be observed. Strongly adsorbed water is also present to some extent in the surface or subsurface region up to elevated temperatures and becomes free at 500 K.

4.3 Neutron time-of-flight spectroscopy – hydrogen dynamics

The first subject that we aimed to address in the neutron scattering study was the impact of the layered structure of MoS₂ on the hydrogen dynamics. This should be observable through a difference between the quasi-elastic spectra obtained at different sample orientations, i.e., with the momentum transfer being parallel or perpendicular to the basal planes of the 2D samples (Fig. 5). In QENS measurements, the diffusive dynamics of sample atoms and molecules lead to a broadening of the elastic scattering line (around energy transfer ħω = 0).² In our QENS spectra no systematic difference between the two sample orientations was found in the normalised data. This means that the properties of hydrogen diffusion parallel and perpendicular to the MoS₂ basal planes are mostly indistinguishable within the spectroscopic window of the present TOF QENS measurements. However, since the lowest exploitable Q values were relatively high, this finding does not put strict constraints on the model of hydrogen diffusion. Based on the structure of molybdenum disulphide we assume two possible diffusion patterns. In the first case hydrogen moves identically between the layers of MoS₂ and through the hexagonal wells (Fig. 5 inset), thus the length scale is of the order of the c unit cell parameter, i.e., 6.15 Å. The second possibility is a motion completely restricted by the sulphur planes, i.e., with no diffusion through the wells. In this case, the characteristic length of diffusion should be no larger than the gap between the sulphur planes of 3 Å. Fig. 5 displays the spectra at Q = 0.85 Å⁻¹, the momentum transfer for which the corresponding scattering angle best satisfies the in-plane and out-of-plane condition.

Fig. 5 Neutron quasi-elastic data obtained at T = 300 K, Q = 0.85 Å⁻¹ and different orientations of MoS₂ single crystals versus momentum transfer vector. The in-plane geometry is represented by the black line and the out-of-plane geometry is represented by the red line. Error bars are displayed as a shading with a corresponding colour. The inset shows the structure of MoS₂, side-view (bc plane) on the left-hand-side and top-view (ab plane) on the right-hand-side.

Now we turn to the in-plane data at various temperatures. Below 200 K any hydrogen motion appeared to be too slow to fit into the spectroscopic window of the instrument, thus yielding no detectable quasi-elastic broadening (Fig. 6a). At 200 K a clearly visible difference between the data and a resolution function occurs, however it is still too small to allow meaningful quantitative fits. A major gain in quasi-elastic signal takes place when heating the sample up to 300 K, and already at this temperature the QENS broadening exhibits clear signs of two different types of motion being present in the sample (Fig. 6b). Taking into account their dependence on the Q vector, which will be discussed later, we interpret these two modes preliminarily as translational and rotational diffusion of hydrogen. The assignment of these modes to specific hydrogen containing species, i.e., hydrogen atoms, hydrogen molecules or water, is difficult. The NRA experiments demonstrate that water does not penetrate into the sample volume and that at least half of the hydrogen is in the form of atoms and molecules. The XPS measurements also show that water changes its binding energy only just below 500 K. Since in the QENS data no radical alterations are observed when heating from 400 to 500 K and since the neutron signal comes from the whole sample simultaneously, we conclude that water may not play a major role in the QENS signal. Now we consider if we can distinguish between the mobility of H₂ and H. Obviously, hydrogen molecules may exhibit both translational and rotational diffusion, whereas atoms are expected to show mostly translational motion with only some degree of frustrated rotational modes. In addition, the XPS experiment indicates that H atoms start moving at ≈400 K, whereas in the QENS data a quasi-elastic broadening can already be observed at room temperature. This signal does not change qualitatively at elevated temperatures, so we conclude that most of the signal in the QENS experiment is caused by H₂ mobility.

Fig. 6 (a) Quasi-elastic neutron scattering scans of MoS₂ single crystals loaded with hydrogen, measured at a range of temperatures from 3 to 500 K at momentum transfer Q = 1.35 Å⁻¹. Statistical error bars are omitted for a clearer visualisation. (b) Closer look at the 300 K data, which are represented by green crosses with error bars, and their approximation by means of a rotational–translational model (red line). Rotational and translational components are shown as purple and blue dashed lines, respectively, the light orange dotted line is the resolution function (hydrogen containing sample at 3 K) and a background term is shown as a black dash-dotted line. The residuals function (difference between measured and calculated points weighted for experimental errors) is displayed on top of the graph.

In view of these considerations the whole set of data at temperatures higher than 200 K was treated by means of the rotational–translational model described in the Theory section above (eqn (9)). However, several spectra at low and high Q were not used in further analysis as the lack of data points in the QENS region (Fig. 2) made the corresponding analysis statistically unstable.

Based on a χ² analysis and on the behaviour of the residuals function (difference between observed and calculated points weighted for statistical errors) we conclude that our model is sufficient for the description of the data and no extra fitting parameters are required. The χ² values of the treated data set are distributed around 0.5 with a standard deviation of 0.1. The residuals functions for all spectra exhibit no systematic deviations from zero. The resulting fitting parameters are summarised in Table 1.

Table 1 Parameters of hydrogen diffusion inside MoS₂ single crystals at various temperatures. The asterisk (*) represents data from the second heating run. Mean square displacement 〈u²〉 was calculated via the Debye–Waller factor. τ_rot is the characteristic time of rotational diffusion. The radius of a confining sphere r_c was obtained by fitting eqn (4) to the A₀(Q) factor. Jump distance and residence time between jumps τ were estimated with the help of several models: Chudley–Elliot (CE), Singwi–Sjölander (SS) and Hall–Ross (HR) model. The subscripts denote, which model was applied for the corresponding numbers

T, K	〈u^2〉, Å^2	τ rot, ps	r c, Å	d CE, Å	d SS, Å	d HR, Å	τ CE, ps	τ SS, ps	τ HR, ps
300	0.23 (1)	2.1 (4)	0.68 (2)	2.8 (2)			0.20 (1)
400	0.40 (1)	2.1 (3)	0.63 (1)	2.1 (2)	2.6 (6)	2.2 (2)	0.29 (3)	0.19 (3)	0.26 (2)
500	0.51 (1)	1.5 (3)	0.74 (1)	2.3 (2)	2.9 (6)	2.5 (2)	0.32 (2)	0.22 (3)	0.29 (2)
300*	0.32 (1)	1.7 (5)	0.71 (3)	1.9 (7)	2.4 (16)	2.1 (7)	0.21 (8)	0.15 (6)	0.20 (5)

---

In the following we discuss the parameters of hydrogen diffusion in MoS₂ obtained by fitting eqn (9) to the data. The change in the overall intensity S₀(Q) at different temperatures that is shown in Fig. 7a reveals that a significant fraction of hydrogen was desorbed during the first heating from 300 to 400 K (of the order of 50%) and then from 400 to 500 K (an additional 10–20%). This desorption happened rather fast during the heating and temperature equilibration phases and did not impact the QENS profile. The hydrogen was not re-adsorbed upon cooling back to 300 K because the sample holder was not tightly closed and the environment around the sample holder was at low pressure during experiments. In the XPS data only minor changes are spotted during heating from room temperature to 400 K. Therefore, the desorbed species are interpreted as surface hydrogen molecules or water, which do not significantly impact on the chemistry and are, thus, not observable in the XPS study. The Q dependence of S₀(Q) allows us to estimate the level of hydrogen vibrations 〈u²〉. The calculated values are moderate for the whole range of temperatures (Table 1).

Fig. 7 Total neutron count characteristics of H/MoS₂ at different temperatures. The asterisk represents second heating run. (a) The overall intensity factor S₀(Q) (scale parameter times Debye–Waller factor) and its fit, represented by solid lines with corresponding colour. (b) The ratio of elastic scattering to the normalised total neutron count S_el(Q).

The factor S_el, displayed in Fig. 7b, shows the ratio of the scattering intensity off the immobile hydrogen and the host material, to the total scattering intensity S_el = (I^im_H + I_MoS2)/(I^total_H + I_MoS2). The term “immobile” refers rather to atoms or molecules, whose movement is very slow and falls into the energy window of the resolution function, than to completely stationary particles. Although the incoherent scattering cross section of MoS₂ is very small as compared to hydrogen, it still may contribute to the S_el factor due to a larger number of scatterers, and thus it is taken into account in the further considerations. Hence, the decrease of S_el with temperature indicates that a part of the immobile hydrogen either becomes mobile, i.e., becomes sufficiently fast to enter the spectroscopic window of the spectrometer, and thus contributes to 1 − S_el, or desorbs. Both room temperature runs, before and after heating to 500 K, give almost identical values of S_el. As no significant changes of S_el are observed between the two room temperature runs, it can be concluded that the balance between mobile H₂ and H atoms that are immobile at room temperature is preserved, which is in line with the above observations of H and H₂ desorption at elevated temperatures. During heating from 300 to 400 K the S_el factor drops and then remains similar upon heating to 500 K, indicating that heating from to 400 to 500 K does not induce any additional hydrogen to become mobile.

The decay of A₀(Q) with Q (Fig. 8a) shows the presence of a certain confinement effect beyond the natural confinement of the rotational diffusion. The approximation of A₀(Q) with a model of a diffusion confined in a sphere gives a good fit with a radius of confinement r_c ≈ 0.7 Å. This length may indicate that hydrogen molecules are, to some extent, caged by four neighbouring sulphur atoms (the corresponding distance calculated on the basis of the MoS₂ structure is around 0.6 Å). In order to better understand the confinement effect two further EISF models were tested, namely, jump reorientations between two and three equivalent sites, respectively, equally spaced on a circle.² The fits of these models are of the same quality as for a spherical confinement model and yield a circle radius of around 0.6 Å, very close to what was obtained with the previous model. This result is not surprising as for such a small confining length all the sensitive information should be contained in a high Q region, which could not be assessed in the current experiment. We would like to stress that the presence of this additional confinement length, in addition to the radius of the hydrogen molecule, does not exclude the existence of long range diffusion, and can be interpreted by a motion in a confined space, occurring between long jumps. For all three models extrapolations of the fit functions to Q = 0 Å⁻¹ yield values slightly higher than 1, despite the fact that the expected value is less or equal to unity. This is an indication of the existence of various hydrogen containing species inside the sample. Indeed, if slow motion of atomic hydrogen falls into the energy window of the rotational mode of H₂, i.e. very close to the resolution function, it would increase the A₀(Q) factor since the observed area of the rotational peak increases (see eqn (9)). As expected the effect is minor as it is caused by the presence of atomic hydrogen. The characteristic time of rotation τ_rot, as theoretically expected, does not depend on the momentum transfer Q (Fig. 8b).

Fig. 8 Diffusion parameters of hydrogen in MoS₂ single crystals at various temperature. The asterisk represents second heating run. (a) A₀(Q) factor and its approximation using a model of translational diffusion confined within an impermeable sphere (eqn (4)). (b) Characteristic time of rotational diffusion τ_rot.

The part of the quasi-elastic broadening that is caused by translational diffusion of molecular hydrogen exhibits a typical behaviour for a jump-diffusion process (Fig. 9) following a DQ² law at low values of momentum transfer and reaching a plateau at high Q.² As was observed earlier,⁴⁶ diffusion models with a continuous distribution of jump lengths, like the Singwi–Sjölander⁴⁷ (SS) and the Hall–Ross⁴⁸ (HR) model, are well suited for treating such cases. Both were tested in this work, together with a general discrete-jump-diffusion model by Chudley and Elliot⁴⁹ (CE). The Singwi–Sjölander and Hall–Ross models consider different distributions of jump distances. The Singwi–Sjölander assumes an exponentially decaying distribution

ρ(r) = r/r₀²·exp^−r/r₀, (10)

while in the Hall–Ross model the jump distances are distributed normally:

Fig. 9 Quasi-elastic broadening due to translational motion of hydrogen in MoS₂Γ_T(Q) plotted versus momentum transfer Q, represented by marks. And its approximated by different translational diffusion models, represented by solid lines: Hall–Ross, Singwi–Sjölander and Chudley–Elliot model for a hexagonal array of jump vectors. Data were collected at various sample temperatures, indicated by different colours in the plot. The asterisk represents the second heating run.

The Chudley–Elliot model, on the other hand, considers jump diffusion between periodic discrete jump sites. In the case of a Bravais lattice, the half-width-at-half-maximum (HWHM) of a translational Lorentzian can be derived from the following general expression:

where n is a number of neighbouring jump sites and the summation goes over jump vectors l_i in a positive hemisphere. For the hexagonal MoS₂ structure it is reasonable to choose jump sites that are distributed with a hexagonal symmetry and a unique jump length l, since it was already shown that the diffusion should be isotropic. Hence, a jump site array was chosen in the form of a right hexagonal pyramid with a jump origin in the centre of the base and jump vectors pointing at the vertices.

At this point we should recall the geometry of the experiment. The in-plane set-up implies that a momentum transfer vector Q is almost parallel to the basal planes of the MoS₂ single crystals, i.e., Q_‖ ≫ Q_⊥. The SS and HR models assume isotropic diffusion, which implies that only the magnitude of the momentum transfer vector is relevant and the geometry of the experiment is insignificant. However, the CE model requires an averaging over all possible directions of Q, which in our case requires an averaging over Q_‖. By substituting the pyramid jump array into (eqn (12)) and averaging over Q_‖ one gets an expression for the width of the translational broadening as displayed in the table below (see Appendix A for details).

Singwi–Sjölander	Hall–Ross	Chudley–Elliott

Here, D is the diffusion coefficient, τ is the residence time between jumps and J₀(x) is a zeroth order Bessel function of the first kind. The mean squared jump length and the diffusion coefficient are related by the Einstein equation². The projection of Q in the current experiment can be evaluated as Q_‖ = Q·sin(ϕ + arcsin(Qλ/4π)).

All three models describe most of the data equally well (in terms of a χ² test). However, the Singwi–Sjölander and Hall–Ross models did not succeed in treating the data obtained during the first 300 K run. The calculated jump distances and residence times τ are listed in Table 1. The first room temperature scan indicates a faster diffusion than all the others. This is probably due to the weakly bound surface hydrogen and/or water, which can be expected to exhibit faster motion and which desorbs at higher temperatures. Therefore, it would be appropriate to exclude the 300 K results for the first temperature run from further consideration. That being so, one can note a slight, but consistent correlation of the jump distance with temperature, which holds for every model. An averaged value of the jump length for all temperatures and models is 2.3 Å. This value may be compared to certain characteristic distances in the MoS₂ structure and to the confinement length r_c, that was obtained from the A₀(Q) factor. The most important distances in the MoS₂ structure are the minimum distance between the surfaces of two neighbouring S atoms of 1.1 Å and a unit cell parameter of 3.15 Å. The obtained jump distance lies between those two values. The characteristic confinement length of about 1.4 Å appears to be smaller than the calculated average jump distance, therefore the observed jumps happen between neighbouring confinement regions and not within them. It is also important to note that this jump distance is much less than the 6.15 Å interval between two layers of MoS₂. Therefore, what is observed in the QENS experiment is likely a movement completely restricted by the MoS₂ layers, rather than a diffusion through the hexagonal wells.

Since molybdenum sulphide is a promising catalyst for the hydrogen evolution reaction (HER), we would like to compare the parameters of hydrogen molecule diffusion in MoS₂ to the currently most active catalyst, platinum. Unfortunately spectroscopic information about hydrogen diffusion on Pt is only available for hydrogen atoms. It was obtained by Graham et al.⁵⁰ for a (111) surface of Pt by means of quasi-elastic helium atom scattering. We also add to the comparison the diffusion of hydrogen molecules on exfoliated graphite, studied by Bahn et al.⁴⁶ with quasi-elastic neutron scattering. Exfoliated graphite has not shown any significant activity in the HER reaction⁴ but its layered structure is similar to MoS₂, hence it can be used as a reference non catalytic material. An Arrhenius plot of the temperature dependence of the diffusion coefficients of these three materials is displayed in Fig. 10. The diffusion coefficient for platinum was extracted via the Chudley–Elliot model, whereas the data for MoS₂ were averaged between the three models studied in this work (Table 2). The 300 K results from the first run were not used for this evaluation also because it was not possible to average over the models. The data for exfoliated graphite were extracted using the Singwi–Sjölander model. Interestingly, the diffusion coefficient of hydrogen atoms on Pt extrapolates quite well to the diffusion coefficient that was observed in the present study for hydrogen molecules at high temperatures. The diffusion coefficient of hydrogen on graphite, on the other hand, is substantially higher than the diffusion coefficient of hydrogen on MoS₂ if extrapolated to a similar temperature. We also note that the activation barrier, E_B, for the hydrogen molecule diffusion in MoS₂ appears to be very small as indicated by the small inclination of the Arrhenius plot. In fact, E_B appears to be similar to the activation barrier in graphite. However, we are unable to assert this finding with confidence and do not want to give a value here, since the number of temperatures that we were able to measure was limited and the error bars of the diffusion coefficients are quite significant.

Fig. 10 Comparison of diffusion coefficients calculated in the current work for hydrogen molecules on MoS₂ (red circles) with the results obtained earlier for hydrogen atoms on platinum⁵⁰ (blue line) and hydrogen molecules on exfoliated graphite⁴⁶ (black triangles). Lines are calculated based on the activation energies and D₀ factors stated in the corresponding papers. Vertical lines for Pt show the range of measured temperatures. The diffusion coefficient at room temperature shown here for MoS₂ represents the second run.

Table 2 Diffusion coefficient of hydrogen in MoS₂ obtained by means of three different diffusion models: Chudley–Elliot D_CE, Singwi–Sjölander D_SS and Hall–Ross D_HR; and an average over the models 〈D〉. Here the 300 K results from the second heating run are shown, as an averaging over the models for the first run is not possible

T, K	D CE × 10^−8, m^2 s^−1	D SS × 10^−8, m^2 s^−1	D HR × 10^−8, m^2 s^−1	〈D〉 × 10^−8, m^2 s^−1
300	2.8 ± 1.5	6.4 ± 5.0	3.6 ± 1.6	4.3 ± 2.7
400	2.5 ± 0.3	5.8 ± 1.5	3.2 ± 0.4	3.8 ± 0.8
500	2.7 ± 0.2	6.3 ± 1.5	3.5 ± 0.4	4.2 ± 0.7

---

4.4 Neutron spin-echo study

Neutron spin-echo (NSE) experiments provide a direct measurement of the time evolution of the space correlation function of the sample (also called the intermediate scattering function S(Q,t)).^2,51 NSE data are complementary to TOF spectroscopy data, because they provide the Fourier transform with respect to energy and time of the TOF data. The spectroscopic window of the technique, in addition, allows the observation of dynamics with a characteristic time from ≈5 ps to several ns and therefore allows us to investigate the diffusion phenomena at lower Q values than the ones studied in the QENS experiment, i.e., testing substantially larger real space distances. NSE spectra are shown in Fig. 11. Starting from 400 K the intermediate scattering function S(Q,t) of the in-plane set-up shows a clearly visible decay that can be fit by an exponential law at Q = 0.15 Å⁻¹, which is an indication of diffusive motion parallel to the surface. We note, that an exponential decay found with NSE is equivalent to a Lorentzian line shape in the TOF experiment, because S(Q,t) is the Fourier transform of S(Q,ω). Within the experimental error bar, we find the same diffusion coefficients through our NSE measurements as for the translational diffusion derived from the QENS data, i.e., (5 ± 3) × 10⁻⁸ m² s⁻¹ at 500 K. While we find a clear indication of diffusion parallel to the basal planes we observe that the signal in the out-of-plane orientation does not change substantially with temperature indicating that no significant diffusion is found perpendicular to the basal planes of MoS₂ within the spectroscopic window of QENS spectroscopy. We would like to stress that this observation is not in contradiction to the TOF findings (where we did not observe a difference between spectra parallel or perpendicular to the basal planes) because the smallest Q values that were exploitable in the TOF spectra represent distances that are in the range of the interlayer distance in MoS₂ of 6.2 Å. On the other hand, the distance corresponding to Q = 0.15 Å⁻¹, probed in the NSE experiment, is about 40 Å. It was not possible to measure the dynamics at higher Q values using NSE, because the neutron beam polarisation was destroyed by the incoherent scattering of the hydrogen for higher Q values: coherent scatterers leave the neutron spin unaltered upon scattering, whereas incoherent scattering atoms cause spin-flips upon interaction with the scattering neutron leading to a final beam polarisation after scattering of on average.⁵¹ We note, that rotational diffusion will not be visible at the low Q measured by our NSE measurement, since the intensity of the rotational diffusion vanishes at small Q.²⁵

Fig. 11 Normalised intermediate scattering function as a function of Fourier time for sample temperatures between 300 and 600 K for a momentum transfer of Q = 0.15 Å⁻¹. In the main graph Q is parallel to the basal planes of MoS₂, whereas in the inset Q is perpendicular to the basal planes. Data have been normalised to the scans at the cryostat base temperature of 2 K.

5 Conclusions

We have performed the analysis of hydrogen diffusion in the bulk of H-loaded MoS₂ single crystals by means of quasi-elastic neutron scattering and neutron spin-echo spectroscopy. The results were compared to the currently most active HER catalyst platinum and a non-catalytic exfoliated graphite. A supplementary study of the surface chemistry of the MoS₂ under temperature treatment and ageing was carried out with the help of X-ray photoelectron spectroscopy and nuclear reaction analysis.

First of all, the nuclear reaction analysis has shown clear signs that hydrogen is able to penetrate the layered structure of MoS₂ and migrate from the surface to the bulk across the basal planes of this layered material, though rather slowly with a diffusion coefficient of the order of 10⁻²¹ m² s⁻¹. The NRA and XPS data also showed that water does not penetrate into the bulk, but can be found adsorbed at the surface or subsurface layer. The QENS study showed a clear translational and rotational diffusion of hydrogen molecules and only a negligible contribution of hydrogen atom diffusion within the spectroscopic window of the present experiments. The low Q data from NSE spectroscopy demonstrated that the hydrogen diffusion between 300 and 600 K is limited to a motion constrained between the layers of MoS₂. The QENS data could be fitted by models that consider a jump diffusion of rotating hydrogen molecules, that remain in a limited space of a dimension of about 1.4 Å between jumps. The calculated jump distance is of the order of 2.5 Å, supporting the idea that the diffusion observed here comes from hydrogen enclosed between the layers.

The calculated diffusion coefficients of hydrogen molecules in molybdenum disulphide were compared to the ones of hydrogen atoms in a catalytic platinum and of hydrogen molecules in non-catalytic exfoliated graphite. According to the Arrhenius plot (Fig. 10) the diffusion rate of H₂ in MoS₂ at room temperature or above is similar to the diffusion coefficient expected for hydrogen atoms on Pt at the same temperatures. This could have interesting technical implications even though the diffusion mechanisms are likely to be very different for hydrogen molecules and atoms. The Arrhenius plot indicates a significantly lower diffusivity for H₂ in MoS₂ than for hydrogen molecules in exfoliated graphite, however it also suggest a low activation barrier, which would be in accordance with the two-dimensional nature of the MoS₂ samples.

Conflicts of interest

There are no conflicts to declare.

Appendix

A Chudley–Elliot model for a pyramid array of jump sites

We start with a convenient representation of the Chudley–Elliot model for jump diffusion on a Bravais lattice:where n is the number of jump vectors, however the sum goes over only half of them, as the model is centrosymmetric. For layered materials it is safe to divide the diffusion into parallel and perpendicular constituents, which would represent the diffusion between and through the layers, respectively. We now consider a hexagonal array of jump sites in a plane, parallel to the basal planes of the material with equal jump distances l, and vertical jumps in a perpendicular direction with a jump distance l_c. Such a distribution can be imagined as two right hexagonal pyramids with their bases glued together, with an origin of a jump in the centre of the base. A visual representation of the jump array applied to the MoS₂ structure is shown in the Fig. 12.

Fig. 12 The structure of MoS₂ and a hexagonal jump array applied to it. Molybdenum and sulphur are displayed as purple and yellow spheres, respectively. Jump sites are shown as black balls and jump vectors as black arrows.

Substituting this jump array into eqn (13) one gets:

Setting a geometry of an experiment such that a momentum transfer vector Q is very close to the basal plane of the pyramid XY, i.e. θ ≈ π/2, allows to rewrite eqn (14) explicitly:

where Q_‖ = ‖Q_‖‖ is a projection of Q on the XY plane. Strictly speaking the current approximation cos(Ql_c cos(θ)) ≈ 1 depends not only on θ, but also on l_c. It remains valid even for smaller θ if l_c ≪ l and fails if l_c ≫ l. In the latter case the l_c containing cosine should be directly calculated.

In the case of an oriented polycrystalline sample, i.e. a set of single crystals whose basal planes are aligned, but may be rotated with respect to each other, the model requires further averaging over all the possible directions of Q_‖, i.e. over ϕ:

Converting cos(ϕ) to sin(π/2 − ϕ) and switching from ϕ to ϕ′ = π/2 − ϕ we come up with an integral definition of a Bessel function of the first kind:

This is the final expression for an in-plane jump diffusion over a hexagonal array of jump sites with a unique jump distance l and residence time τ.

Acknowledgements

This work is based upon experiments performed at the TOFTOF instrument operated by Technische Universität München at the Heinz Maier-Leibnitz Zentrum (MLZ), Garching, Germany. This project has received support from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 730872 CALYPSO plus for the work at the HZB BESSY II synchrotron, Berlin, Germany, and under grant agreement no. 731096 FILL2030 for covering the open access fee. The authors want to thank HZB for the allocation of synchrotron radiation beam time at BESSY II via project 192-08548 CR. Neutron spin-echo experiments were carried out at Institut Laue-Langevin (ILL), Grenoble, France. The authors want to thank HZB, MLZ and ILL for the allocation of synchrotron and neutron beam time. The help provided by Dimitri Renzy and Simon Baudoin at ILL is gratefully acknowledged. The authors would like to thank Laurène Gallot for her contributions in the framework of an Erasmus internship. The work was financed within the scope of the joint PhD project ref. ILL-1709.1 between Institut Laue – Langevin and Westfälische Hochschule.

Notes and references

1. A. P. Graham and J. P. Toennies, Collective Diffusion on Surfaces: Correlation Effects and Adatom Interactions, Dordrecht, 2001, pp. 1–10.
2. M. Bée, Quasielastic neutron scattering: principles and applications in solid state chemistry, biology and materials science, Adam Hilger, Bristol and Philadelphia, 1988.
3. R. Hempelmann, Quasielastic neutron scattering and solid state diffusion, Oxford University Press, Oxford, 2000.
4. M. Carmo, D. L. Fritz, J. Mergel and D. Stolten, Int. J. Hydrogen Energy, 2013, 38, 4901–4934.
5. Y. N. Regmi, G. R. Waetzig, K. D. Duffee, S. M. Schmuecker, J. M. Thode and B. M. Leonard, J. Mater. Chem. A, 2015, 3, 10085–10091.
6. Y. Yan, B. Xia, Z. Xu and X. Wang, ACS Catal., 2014, 4, 1693–1705.
7. J. D. Benck, T. R. Hellstern, J. Kibsgaard, P. Chakthranont and T. F. Jaramillo, ACS Catal., 2014, 4, 3957–3971.
8. Z. Chen, D. Cummins, B. N. Reinecke, E. Clark, M. K. Sunkara and T. F. Jaramillo, Nano Lett., 2011, 11, 4168–4175.
9. X. Dai, K. Du, Z. Li, H. Sun, Y. Yang, W. Zhang and X. Zhang, Int. J. Hydrogen Energy, 2015, 40, 8877–8888.
10. X. Zou and Y. Zhang, Chem. Soc. Rev., 2015, 44, 5148–5180.
11. P. A. Bertand, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 44, 5745–5749.
12. D. Voiry, M. Salehi, R. Silva, T. Fujita, M. Chen, T. Asefa, V. B. Shenoy, G. Eda and M. Chhowalla, Nano Lett., 2013, 13, 6222–6227.
13. D. Merki, S. Fierro, H. Vrubel and X. Hu, Chem. Sci., 2011, 2, 1262–1267.
14. J. D. Benck, Z. B. Chen, L. Y. Kuritzky, A. J. Forman and T. F. Jaramillo, ACS Catal., 2012, 2, 1916–1923.
15. C. Rozain and P. Millet, Electrochim. Acta, 2014, 131, 160–167.
16. J. O. Bockris, Modern electrochemistry 2B: electrodics in chemistry, engineering, biology and environmental science, Springer Science & Business Media, 1998, vol. 2.
17. I. N. Yakovkin and N. Petrova, Chem. Phys., 2014, 434, 20–24.
18. K. M. S. Takagi and T. Gotoh, Appl. Surf. Sci., 1999, 144, 278–282.
19. V. F. Sears, Can. J. Phys., 1967, 45, 237–254.
20. V. F. Sears, Can. J. Phys., 1966, 44, 1299–1311.
21. T. R. Prisk, M. S. Bryan and P. E. Sokol, Phys. Chem. Chem. Phys., 2014, 16, 17960–17974.
22. A. J. O'Malley and C. R. A. Catlow, in Neutron scattering – applications in biology, chemistry, and materials science, Experimental Methods in the Physical Sciences, ed. F. Fernandez-Alonso and D. Price, Elsevier, San Diego, 2017, vol. 49, pp. 349–401.
23. I. Calvo-Almazán, S. Miret-Artés and P. Fouquet, J. Phys.: Condens. Matter, 2012, 24, 104007.
24. H. Jobic, Phys. Chem. Chem. Phys., 2016, 18, 17190–17195.
25. S.-H. Chen, P. Gallo, F. Sciortino and P. Tartaglia, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1997, 56, 4231–4243.
26. M. Bée, Phys. B, 1992, 182, 323–336.
27. F. Volino and A. Dianoux, Mol. Phys., 1980, 41, 271–279.
28. A. N. Vavreck and W. Thompson, Q. Appl. Math., 1984, 42, 321–324.
29. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
30. J. F. Ziegler, M. D. Ziegler and J. P. Biersack, Nucl. Instrum. Methods Phys. Res., Sect. B, 2010, 268, 1818–1823.
31. J. F. Ziegler, Program SRIM (Stopping and Range of Ions in Matter), http://www.srim.org.
32. F. Traeger, M. Kauer, C. Wöll, D. Rogalla and H.-W. Becker, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 075462.
33. H.-W. Becker and D. Rogalla, in Nuclear Reaction Analysis, ed. H. Fritzsche, J. Huot and D. Fruchart, Springer International Publishing, Cham, 2016, pp. 315–336.
34. A. Nefedov and C. Wöll, in Advanced Applications of NEXAFS Spectroscopy for Functionalized Surfaces, ed. G. Bracco and B. Holst, Springer Berlin Heidelberg, Berlin, Heidelberg, 2013, pp. 277–303.
35. D. A. Shirley, Phys. Rev. B: Solid State, 1972, 5, 4709–4714.
36. J. Yeh and I. Lindau, At. Data Nucl. Data Tables, 1985, 32, 1–155.
37. H. M.-L. Zentrum, J. Large-Scale Res. Facil., 2015, 1, A15.
38. O. Arnold, J. Bilheux, J. Borreguero, A. Buts, S. Campbell, L. Chapon, M. Doucet, N. Draper, R. F. Leal, M. Gigg, V. Lynch, A. Markvardsen, D. Mikkelson, R. Mikkelson, R. Miller, K. Palmen, P. Parker, G. Passos, T. Perring, P. Peterson, S. Ren, M. Reuter, A. Savici, J. Taylor, R. Taylor, R. Tolchenov, W. Zhou and J. Zikovsky, Nucl. Instrum. Methods Phys. Res., Sect. A, 2014, 764, 156–166.
39. F. Traeger and P. Fouquet, Hydrogen diffusion in MoS ₂ , Institut Laue-Langevin (ILL) , 2016 DOI:10.5291/ILL-DATA.7-05-460.
40. B. Farago, Phys. B, 1998, 241–243, 113–116.
41. M. P. Seah and W. A. Dench, Surf. Interface Anal., 1979, 1, 2.
42. J.-G. Choi and L. Thompson, Appl. Surf. Sci., 1996, 93, 143–149.
43. D. Ganta, S. Sinha and R. T. Haasch, Surf. Sci. Spectra, 2014, 21, 19–27.
44. L. Trotochaud, A. R. Head, S. Pletincx, O. Karsloğlu, Y. Yu, A. Waldner, L. Kyhl, T. Hauffman, H. Terryn, B. Eichhorn and H. Bluhm, J. Phys. Chem. B, 2018, 122, 1000–1008.
45. D. I. Patel, D. Shah, S. Bahr, P. Dietrich, M. Meyer, A. Thißen and M. R. Linford, Surf. Sci. Spectra, 2019, 26, 014026.
46. E. Bahn, L. A. H. Giraldo, V. Kuznetsov, I. Calvo-Almazán, M. Zbiri, M. M. Koza, T. C. Hansen, P. F. Henry, A. Lapp, S. Pouget, M. Mesa and P. Fouquet, Carbon, 2020, 166, 307–315.
47. K. S. Singwi and A. Sjölander, Phys. Rev., 1960, 119, 863–871.
48. P. L. Hall and D. Ross, Mol. Phys., 1981, 42, 673–682.
49. C. T. Chudley and R. J. Elliott, Proc. Phys. Soc., 1961, 77, 353–361.
50. A. P. Graham, A. Menzel and J. P. Toennies, J. Chem. Phys., 1999, 111, 1676–1685.
51. Neutron Spin Echo, ed. F. Mezei, Springer, Berlin, 1980, vol. 128.

---