Hydrogen Transport Between Layers of Transition Metal-Dichalcogenides

Hydrogen is a crucial source of green energy and has been extensively studied for its potential usage in fuel cells. The advent of two-dimensional crystals (2DCs) has taken hydrogen research to new heights, enabling it to tunnel through layers of 2DCs or be transported within voids between the layers, as demonstrated in recent experiments by Geim's group. In this study, we investigate how the composition and stacking of transition-metal dichalcogenide (TMDC) layers influence the transport and self-diffusion coefficients (D) of hydrogen atoms using well-tempered metadynamics simulations. Our findings show that modifying either the transition metal or the chalcogen atoms significantly affects the free energy barriers (Delta F) and, consequently, the self-diffusion of hydrogen atoms between the 2DC layers. In the Hh polytype (2H stacking), MoSe2 exhibits the lowest Delta F, while WS2 has the highest, resulting in the largest D for the former system. Additionally, hydrogen atoms inside the RhM (or 3R) polytype encounter more than twice lower energy barriers and, thus, much higher diffusivity compared to those within the most stable Hh stacking. These findings are particularly significant when investigating twisted layers or homo- or heterostructures, as different stacking areas may dominate over others, potentially leading to directional transport and interesting materials for ion or atom sieving.

∆F and higher self-diffusion coefficients (D) compared to their S-or W-based counterparts, making MoSe 2 systems ideal for H atom transport.Among different high-symmetry stackings, the R types have up to twice lower ∆F than the H types.These results are of significance in the exploration of 2D TMDC materials as potential transport channels for energy devices, such as fuel cells.In particular, twisted homo-and heterostructures will be of great interest, as differently sized domains of the mentioned high-symmetry stackings, with varying ∆F , will lead to local transport channels.

Effect of Elemental Composition
Initially, we examined the diffusion of H-atoms between layers of various Group-6 TMDC bulk materials in their H h h (2H) forms as depicted in Figure 4.In addition, we also investigated metallic NbS 2 (in H X h as the most stable polytype) and semiconductor main-group β-InSe for comparison.Our primary objective was to comprehend the influence of different elements in the layers on ∆F for hydrogen transfer and its diffusivity.
The optimized lattice parameters and the interatomic distances of all systems investigated in this section are given in Table I.The in-plane lattice parameters agree within a maximum of 0.6% with the experimental data, while somewhat larger errors were obtained for the out-of-plane lattice parameters, with maximum overestimation of 2.4% for MoSe 2 .[13] All systems have approximately 3 Å of space between the layers (d I ; see Table I), providing sufficient room to accommodate H atoms.The H atoms bind to the chalcogen atoms (X) with bond lengths of around 1.4 Å for S and approximately 1.5-1.6Å for Se.In all cases, the H atoms bound to X atoms point towards the interlayer void and the middle of hexagonal rings.

TABLE I. Structural properties of all investigated bulk materials in the H h
h forms (H X h for NbS2): lattice vectors,a, b, and c (in Å); bond lengths between metal and chalcogen atoms, dM−X , and between chalcogen atoms within a single layer, dX−X , (in Å); interlayer distances measured between metal atoms, dM−M , and between inner chalcogen atoms in the neighbouring layers, dI , (in Å); and the bond lengths between chalcogen and hydrogen atoms, dX−H , (in Å), cf. Figure 5b.For InSe, dM−M was measured between inner In atoms.Free energy barriers, ∆F , (in meV) and the self-diffusion coefficients, D, (in cm 2 s −1 × 10 −3 ) of all studied systems.We employed WTMetaD simulations, with two collective variables (CVs) defined for each system (see Section for detailed information), to calculate the free-energy (F ) landscape.The landscape exhibited two minima and an energy barrier (∆F ) separating them.As an example, Figure 2 displays the F landscape for bulk WS 2 along with its corresponding ∆F value.The energy landscapes for all the other systems share similar characteristics and primarily vary in terms of their ∆F values.For completeness, these energy landscapes are provided in Figure S1 in the Supporting Information.The values of ∆F for all systems are given in Table I.
Among the four Group-6 TMDC bulk materials, the highest ∆F belongs to WS 2 and the lowest to MoSe 2 .We noticed that Se-based TMDC have lower ∆F than their S-based counterparts and Mo-based systems have lower ∆F than the corresponding W-based materials.From Equation 2 in Section we know that lower the energy barrier the larger the self-diffusion coefficient.This allows a conclusion that H atoms diffuse easier in the Se-and Mo-based TMDCs.The significance of interlayer spacing cannot be understated, as it plays a pivotal role in self-diffusion.Our simulations, however, tend to overestimate the interlayer distance in Se-based systems to a greater extent compared to S-based materials.To address this, we extended our investigation on MoSe 2 bulk and adjusted the c lattice vector to match the experimental value (12.910Å). [13] By doing so, the H-binding sites are brought into closer proximity, reducing the ∆F even further.Consequently, the self-diffusion of MoSe 2 increases, supporting our conclusions.Our results agree well with previously published work, [11] where the authors obtained ∆F = 90 meV and D = 1.50 × 10 −4 cm 2 s −1 for MoS 2 bulk in H h h form.The small differences that are observed might be due to the 0.4% difference in the interlayer distance.
It is important to note that the defined collective variables (CVs) in WTMetaD simulations do not allow for the investigation of diffusion direction; they only result in a Brownian motion of H atoms. Nevertheless, these CVs enable the estimation of D (self-diffusion coefficients).
In agreement with our previous report, [11] the H atoms are transported between the layers of TMDC materials in a zigzag manner, binding to the inner chalcogen atoms of neighbouring layers (see Figure S2).The transfer between layers is also supported by low-energy shear modes.Additionally, our previous results shows that energy barrier for proton transport between the layers is only slightly higher than that for H atom, which should be similar for other TMDCs.This might be important when investigating proton transport for potential applications, such as fuel cells.
For comparison, we investigated two other layered materials, NbS 2 and β-InSe.NbS 2 is a metallic TMDC from Group-5.Nb atoms have one electron less than Mo atoms, what results in a ferromagnetic material.This system was FIG. 2. Exemplary free-energy (F ) landscape obtained from WTMetaD simulations of a single H atom moving between layers of H h h bulk WS2.The collective variables (CVs) were defined as coordination number (CN) of H atoms with the inner chalcogen atoms X1-H (CV1) and X2-H (CV2; cf. Figure 5), which correspond to the x and y axes.The color legend corresponds to F in eV.investigated in the H X h polytype, which is more stable than the H h h one.[14] The interlayer distance in this material is smaller than in the Group-6 TMDCs, however, ∆F is by almost 20 meV higher than that of WS 2 , resulting in a very low D (only 9 × 10 −5 cm 2 s −1 ).The second material, β-InSe (H h h stacking), is a well-studied layered system, which is a quasi-direct band-gap semiconductor.[15] It has the highest ∆F of all studied systems, over 300 meV, which is more than twice larger than that of NbS 2 , resulting in a negligible self diffusion (1.1 × 10 −7 cm 2 s −1 ).This might be due to the fact that H atom forms much shorter bonds with Se in this material than in the Group-6 TMDCs, which should result in different frequency of vibration (v 0 ).

Effect of Stacking
In this section, we investigated only bulk MoS 2 in different high-symmetry structural stacking, namely , and R M h (which is the same as R X h in homo-layered materials; see Figure 4).[16,17] Additionally, we considered the 3R stacking, which requires three layers in the bulk unit cell.
Among the five high-symmetry stackings, there are three that are low-energy systems: [18] the most stable H h h , and the R M h and H X h , which are only about 2 meV per unit cell less stable.The other two stackings are about 53 meV per unit cell less stable than H h h .It, therefore, was not surprising that, during the WTMetaD simulations with H atoms between the layers, the high-energy stackings (R h h and H M h ) reverted back to their corresponding R M h and H X h forms.Consequently, Table II presents the results for the most stable structural stackings of MoS 2 .The corresponding free-energy landscapes are shown in Figure S1.
While the phononic properties of the H h h and 3R MoS 2 polytypes show very similar frequency for shear modes (31.85 cm −1 and 32.85 cm −1 for 3R-and H h h polytypes, respectively), [19] these materials differ significantly in the diffusivity of H atoms.The H h h polytype has the highest ∆F and the lowest D among all the stackings.On the other hand, the lowest energy barrier and the highest diffusion belong to the R forms.We expect similar results for the other investigated TMDCs.
This finding holds significant importance when considering twisted bilayer TMDCs with moiré structures and high- ; bond lengths between metal and chalcogen atoms, dM−X , and between chalcogen atoms within a single layer, dX−X , (in Å); interlayer distances measured between metal atoms, dM−M , and between inner chalcogen atoms in the neighbouring layers, dI , (in Å); and the bond lengths between chalcogen and hydrogen atoms, dX−H , (in Å), cf. Figure 5b.Free energy barriers, ∆F , (in meV) and the self-diffusion coefficients, D, (in cm 2 s −1 × 10 −3 ) of all studied systems.. symmetry stacking domains [18] for ion-or atom-transporting materials in energy technologies.Depending on the twist angle between layers, different sizes of the mentioned high-symmetry stackings are formed, which may lead to localized transport channels and directional diffusion (a subject of our ongoing investigations).

CONCLUSIONS
In summary, our study focused on investigating the diffusivity of H atoms between layers of TMDC bulk materials, utilizing the interstitial voids that provide ample space for H-atom or -ion transport.We observed that Group-6 TMDCs containing Se or Mo atoms exhibit lower free-energy barriers for H-atom transport compared to their S-or W-based counterparts.Additionally, we explored the impact of layer stacking in MoS 2 on the self-diffusion coefficients.Our findings indicate that the R polytypes (3R and R M h ), present when the twist angle between layers θ → 0 • , demonstrate higher diffusivity than the slightly more stable H forms (present when θ → 60 • ).
These discoveries hold significant importance for future studies concerning transport properties in materials with twisted layers, leading to moiré structures and the formation of differently sized domains with various high-symmetry stackings.Such arrangements have the potential to form localized transport channels in these systems.
Furthermore, our results demonstrate the potential for optimizing intriguing materials based on TMDC homo-or hetero-layers with different twist angles for applications involving ion-or atom-transport, such as in fuel cells.
While our study focused on the diffusion of protiums, it is reasonable to expect that the results will be quite similar for protons as well.[11].

Computational Details
In this work, we investigated different semiconducting TMDCs from Group 6 (see Figure 4), with a general formula MX 2 (M -Mo or W; X -S or Se).MoS 2 was also simulated in the following high-symmetry stackings with two layers in the unit cell forming a bulk material: H h h (also known as 2H stacking), R M h (also known as 3R stacking if three layers are in the unit cell), R h h , H X h , and H M h .For comparison, we have investigated two other layered materials, namely, metallic NbS 2 and semiconducting beta-InSe (see Figure 4).All systems were represented as bulk materials and periodic boundary conditions were used accordingly.To investigate the diffusion of H atom between layers of these materials, we used the 6×6×1 supercels for all studied systems.
All systems were fully relaxed (atomic positions and lattice vectors) using density functional theory (DFT) with Perdew-Burke-Ernzerhof exchange-correlation functional (PBE), [20] Gaussian augmented plane wave (GAPW) basis sets, and Grimme's D3 London dispersion correction, [21] as implemented in the CP2K package.[22] The Quickstep method was employed with Goedecker-Teter-Hutter (GTH) [23] pseudopotentials together with DZVP-MOLOPT-GTH-SR basis set for all elements except hydrogen, which was represented using the full potential.We considered NbS 2 as a ferromagnetic metallic system.In all calculations performed using supercells, we used Γ-point approach.
Since we were interested only in the free-energy barriers and the corresponding self-diffusion coefficients, we omitted using spin-orbit coupling, which is important when investigating electronic properties.
We performed well-tempered metadynamics (WTMetaD) [24] simulations, as implemented in the CP2K package, to obtain the free-energy barriers (∆F ) for H-atom diffusion between the layers of studied materials.Here, the canonical NVT ensemble was employed with the temperature set to 300 K and CSVR (canonical sampling velocity rescaling) thermostat [25] with 0.5 fs timestep.Each trajectory was at least 45 ps long, to ensure convergence.Each WTMetaD simulation was preceded by standard Born-Oppenheimer molecular dynamics (BOMD) simulations (MD; NVT, 300 K, time step of 0.5 fs) with duration of at least 5 ps to ensure thermal equilibrium.We defined two collective variables (CVs; see Figure 5) to trigger the process of H-atom transfer between the layers: all chalcogen atoms in a given system were divided into three different kinds: X1 (inner layer with bound H atom), X2 (inner layer neighboring X1), and X (outermost layers).One of the collective variables, CV1, was defined as the coordination number (CN) of H to X1, CN H−X1 , and the other, CV2, as CN of H to X2, CN H−X2 .The form of the coordination function, CN H−X , is defined in CP2K as follows: where r ij is the distance between two sites that bind the H atom. R 0 , the cutoff distance in the CN , is set to 3.0 a 0 for X1-H and X2-H distances in all considered materials.[11] Note that we only considered transfer of H atoms between X1 and X2 because ∆F for transfer within the same layer is much larger.[11] I).
Gaussian hills were applied every 150 steps to all TMDCs and every 100 steps to InSe during WTMetaD calculations.From WTMetaD simulations, we obtained the free-energy barriers for H atom transfer between layers, ∆F , which was used to calculate the self-diffusion coefficients, D: [11,26] where k B is the Boltzmann's constant, T is the temperature (here set to 300 K), D 0 is a pre-factor, [11,27] ∆r is the distance between H atoms in two closest binding sites, v 0 is the frequency of the stretching modes between X and H atoms, and q i can be 2, 4 or 6 for the diffusion dimensionality of 1D, 2D or 3D, respectively (here set to q i = 4).The frequency, v 0 , of S-H and Se-H were taken from literature to be 2570 and 2300 cm −1 , respectively, assuming that the second neighbours do not affect it.[28]

FIG. 1 .
FIG. 1. Top and side views of all considered structural high-symmetry stackings.Unit cells are marked with black lines.Blue -metal atom (Mo, W, Nb), yellow -chalcogen atom (S or Se), violet -In atom.

FIG. 3 .
FIG. 3. (a) An exemplary representation of collective variables X1-H and X2-H (see details in Methods section) within a TMDC material, shown from the side view.(b) Definition of selected bond lengths and interatomic distances (cf.TableI).

FIG. 5 .
FIG. 5. (a) Schematics of the H-atom transfers path between layers of a TMDC material.(b) Exemplary X-H bond distance change during WTMetaD simulations in WS2 bulk.

TABLE II .
Structural properties of MoS2 bulk systems with different high-symmetry stackings: lattice vectors,a, b, and c (in Å)