Atomistic understanding of hydrogen bubble-induced embrittlement in tungsten enabled by machine learning molecular dynamics

TL;DR

This paper develops a machine learning potential for tungsten-hydrogen systems, enabling atomistic simulations that reveal how hydrogen bubbles cause embrittlement and brittle fracture in tungsten.

Contribution

It introduces the first transferable machine learning model capable of capturing hydrogen trapping, H₂ formation, and bubble-induced fracture in tungsten at the atomistic level.

Findings

Hydrogen clusters nucleate along {100} planes near voids.

Hydrogen bubbles promote brittle fracture and crack formation.

The model achieves near first-principles accuracy with high efficiency.

Abstract

Hydrogen bubble formation within nanoscale voids is a critical mechanism underlying the embrittlement of metallic materials, yet its atomistic origins remains elusive. Here, we present an accurate and transferable machine-learned potential (MLP) for the tungsten-hydrogen binary system within the neuroevolution potential (NEP) framework, trained through active learning on extensive density functional theory data. The developed NEP-WH model reproduces a wide range of lattice and defect properties in tungsten systems, as well as hydrogen solubility, with near first-principles accuracy, while retaining the efficiency of empirical potentials. Crucially, it is the first MLP capable of capturing hydrogen trapping and H\textsubscript{2} formation in nanovoids, with quantitative fidelity. Large-scale machine-learning molecular dynamics simulations reveal a distinct aggregation pathway where planar hydrogen clusters nucleate and grow along \{100\} planes near voids, with hexagonal close-packed structures emerging at their intersections. Under uniaxial tension, these aggregates promote bubble fracture and the development of regular \{100\} cracks, suppressing dislocation activity and resulting in brittle fracture behavior. This work provides detailed atomistic insights into hydrogen bubble evolution and fracture in nanovoids, enables predictive modeling of structural degradation in extreme environments, and advances fundamental understanding of hydrogen-induced damage in structural metals.

Tables

Table 1. Table 1: Static properties of element W and H, and binary W-H systems predicted by BOPs [ 48 , 49 ] , EAM potentials [ 26 , 50 ] , DP -WH [ 38 ] , and NEP -WH models as compared to experimental (Exp.) and/or DFT reference values. For DFT results, both values from previous studies (Previous DFT ) and those obtained in the present work ( DFT ) are included. All lattice parameters a a and bond length L b L_{\rm b} are given in Å \rm\AA , elastic constants C i ​ j C_{ij} in GPa, and all other energies E E in eV, expect for surface formation energy E f surf E_{\rm f}^{\rm surf} , which is given in J m -2 . a bcc a_{\rm bcc} and a fcc a_{\rm fcc} denote the lattice parameters of BCC and FCC W, respectively. C i ​ j C_{ij} denote the elastic constants of BCC W. E f SIA E_{\rm f}^{\rm SIA} refers to the formation energy of interstitial dumbbells with different orientations. E f vac E_{\rm f}^{\rm vac} and E m vac E_{\rm m}^{\rm vac} represent the formation energy and migration energy of a mono-vacancy, respectively. E b divac − 1 ​ N ​ N E_{\rm b}^{\rm divac-1NN} and E b divac − 2 ​ N ​ N E_{\rm b}^{\rm divac-2NN} are the binding energies for the first- and second-nearest neighbor divacancy configurations, respectively. E f surf E_{\rm f}^{\rm surf} denote formation energies of low-index free surfaces with different orientations. E f tis E_{\rm f}^{\rm tis} and E f ois E_{\rm f}^{\rm ois} represent the solution energies of a hydrogen atom at TIS and OIS , respectively. E m t → t E_{\rm m}^{\rm t\rightarrow t} denotes H migration energy between two nearest-neighbor TISs . E b Tplanar E_{\rm b}^{\rm Tplanar} , E b Oplanar E_{\rm b}^{\rm Oplanar} , and E b rock − salt E_{\rm b}^{\rm rock-salt} denote the binding energies of H in the TIS and OIS planar self-clusters, and the rock-salt structures, respectively. E b surf E_{\rm b}^{\rm surf} are the binding energies of a hydrogen atom on the most stable adsorption site of low-index surfaces. E f H 2 E_{\rm f}^{\rm H_{2}} denotes the formation energy of a H 2 molecule in vacuum. L b H 2 L_{\rm b}^{\rm H_{\rm 2}} denotes the bond length of H 2 molecule. Values taken from previous studies are cited with references; all other values are obtained in the present work.

Potential	BOP-Juslin	BOP-Li	EAM-Wang	EAM-Mason	DP-WH	Exp.	Previous DFT	DFT	NEP-WH
Element W:
$a_{\mathrm{bcc}}$	3.165 [48]	3.165 [49]	3.140 [26]	3.145 [50]	3.186 [38]	3.165 [51]	3.165 [52], 3.172 [53]	3.175	3.176
$a_{\mathrm{fcc}}$	4.005 [48]	4.002 [49]	4.053 [26]	3.975	4.040		4.054 [52], 4.023 [54]	4.028	4.029
$C_{11}$	542 [49]	515 [49]	544 [26]	513 [50]	516 [38]	522 [51]	523 [52], 517 [53]	-	534
$C_{12}$	191 [49]	203 [49]	208 [26]	201 [50]	201 [38]	204 [51]	203 [52], 198 [53]	-	199
$C_{44}$	162 [49]	162 [49]	160 [26]	163 [50]	146 [38]	161 [51]	160 [52], 142 [53]	-	156
$E_f^{\mathrm{SIA}-⟨111⟩}$	9.62 [49]	9.33 [49]	10.52 [26]	9.74 [50]	10.16 [38]		10.53 [52], 9.55 [55]	9.99	10.21
$E_f^{\mathrm{SIA}-⟨110⟩}$	8.77 [49]	9.53 [49]	10.82 [26]	9.62	10.35 [38]		10.82 [52], 9.84 [55]	10.27	10.49
$E_f^{\mathrm{SIA}-\mathrm{tetra}}$	8.60 [49]	10.75 [49]	11.94 [26]	10.60	11.81		12.27[52], 11.05 [55]	11.47	10.95
$E_f^{\mathrm{SIA}-⟨100⟩}$	8.93 [49]	12.01 [49]	12.86 [26]	9.97 [50]	12.03 [38]		12.87 [52], 11.49 [55]	11.96	11.07
$E_f^{\mathrm{SIA}-\mathrm{octa}}$	9.92 [49]	12.05 [49]	12.64 [26]	11.25	11.91		13.11 [52], 11.68 [55]	12.07	11.09
$E_f^{\mathrm{vac}}$	1.68 [48]	3.52 [49]	3.49 [26]	3.63 [50]	3.30 [38]	$3.67$ [56]	3.49 [52], 3.35 [57]	3.17	3.36
$E_m^{\mathrm{vac}}$	1.77 [49]	1.81 [49]	1.85 [58]	1.75 [50]	1.95	$1.78$ [56]	1.78 [55], 1.73 [57]	-	2.18
$E_b^{\mathrm{divac}-1NN}$	0.36 [49]	0.62 [49]	0.50 [26]	0.11	0.20		-0.05 [59], 0.01 [60]	0.02	0.35
$E_b^{\mathrm{divac}-2NN}$	0.08 [49]	0.15 [49]	0.39 [26]	0.03	-0.21		-0.39 [59], -0.35 [60]	-0.33	-0.09
$E_f^{\mathrm{surf}-\left(100\right)}$	1.446 [49]	3.157 [49]	2.721 [61]	3.973 [50]	3.900 [38]	2.990 [62]	4.021 [53], 4.149 [63]	3.984	4.007
$E_f^{\mathrm{surf}-\left(110\right)}$	0.931 [49]	2.319 [49]	2.306 [61]	3.524 [50]	3.302 [38]	3.220 [64]	3.268 [53], 3.396 [63]	3.241	3.180
$E_f^{\mathrm{surf}-\left(111\right)}$	1.720 [48]	3.222 [61]	2.963 [61]	4.261 [50]	3.563 [38]		3.556 [53], 3.829 [63]	3.457	3.491
Binary W-H:
$E_s^{\mathrm{tis}}$	1.04 [26]	0.86 [26]	1.05 [26]	0.80 [50]	0.95 [38]	$1.04$ [65]	0.88 [66], 0.85 [67]	0.92	0.92
$E_s^{\mathrm{ois}}$	1.40 [26]	1.18 [26]	1.40 [26]	1.16 [50]	1.35 [38]		1.26 [66], 1.33 [67]	1.30	1.29
$E_m^{t\to t}$	$0.34$ [48]	0.23 [26]	0.22 [26]	0.22 [50]	0.20 [38]	0.39 [65]	0.20 [66], 0.21 [67]	-	0.19
$E_b^{\mathrm{Tplanar}}$	0.01 [38]	0.97 [38]	1.04 [38]	1.27	0.65 [38]		0.64 [38]	0.61	0.56
$E_b^{\mathrm{Oplanar}}$	- [38]	0.67 [38]	0.65 [38]	0.95	0.87 [38]		0.87[38]	0.87	0.74
$E_b^{\mathrm{rock}-\mathrm{salt}}$	-0.36 [38]	0.85 [38]	0.36 [38]	1.14	0.50 [38]		0.50 [38]	0.48	0.52
$E_b^{\mathrm{surf}-\left(100\right)}$	0.39 [38]	-2.29 [38]	-0.18 [38]	-0.31 [50]	-0.90 [38]		-0.87 [68], -0.91 [69]	-	-0.80
$E_b^{\mathrm{surf}-\left(110\right)}$	0.09 [38]	-1.08 [38]	0.28 [38]	-0.49 [50]	-0.76 [38]		-0.76 [68] -0.75 [69]	-	-0.75
$E_b^{\mathrm{surf}-\left(111\right)}$	0.51 [38]	-0.47 [38]	-0.01 [38]	-0.39 [50]	-0.81 [38]		-0.62 [69]	-	-0.59
Element H:
$E_f^{H_2}$	-4.75 [26]	-4.75 [26]	-4.72 [26]	-4.24 [50]	-4.53 [38]	-4.52 [50]	-4.53 [38]	-4.56	-4.59
$L_b^{H_2}$	0.74 [48]	0.74 [49]	0.73 [26]	0.74 [50]	0.75 [38]		0.75 [38]	0.75	0.74

Equations

$$U_i=\sum _{\mu =1}^{N_{\mathrm{neu}}}{w}_{\mu }^{(1)}\tanh \left(\sum _{\nu =1}^{N_{\mathrm{des}}}{w}_{\mu \nu }^{(0)}{q}_{\nu }^i-b_{\mu }^{(0)}\right)-b^{(1)}.$$

Full text

INTRODUCTION Hydrogen-induced degradation has long been a critical challenge in materials science, often leading to premature or catastrophic failure of structural metals [1]. Over the past decades, multiple mechanisms have been proposed to explain hydrogen embrittlement [2], with the formation and evolution of internal high-pressure bubbles or voids recognized as a central factor [3, 4, 5, 6]. Tungsten (W) is a representative system for studying hydrogen-induced embrittlement, as it is widely considered a leading candidate for plasma-facing components in future fusion reactors owing to its exceptional physical properties [7, 8]. As a typical metal exposed to hydrogen-rich environments, tungsten has exhibited extensive hydrogen bubble and failure phenomena in experiments [9, 10, 11, 12, 13, 14]. Despite extensive experimental evidence of hydrogen bubble–induced embrittlement [11, 13, 14], the atomic-scale processes governing bubble nucleation and growth remain poorly understood, due to the challenges of in situ characterization of their internal structure and pressure.

Atomistic simulations provide a powerful complement to experiments. First-principles calculations at the density-functional theory (DFT) level can accurately describe hydrogen interactions with vacancies, grain boundarys (GBs), surfaces, and dislocations at zero temperature [15, 16, 17, 18, 19, 20, 21, 22]. However, their prohibitive cost in both system size and timescale prevents direct resolution of the mechanisms responsible for hydrogen embrittlement. In contrast, molecular dynamics (MD) simulations based on empirical interatomic potentials can access much larger systems and extended timescales, and have yielded valuable insights into hydrogen embrittlement[23, 24, 25]. However, their reliability is fundamentally limited by the quality of the underlying potentials, which show large discrepancies in predicting H trapping at defect[26, 27, 28] and fail to capture H-nanovoids interactions or reproduce experimentally observed gigapascal-level internal H2 [29, 26]. These deficiencies cast serious doubt on their ability to predict key processes such as hydrogen bubble nucleation and evolution.

By learning directly from quantum-mechanical electronic structure data, machine-learned potentialss (MLPs) [30, 31, 32, 33, 34] have emerged over the past two decades as a powerful framework that overcomes long-standing accuracy-efficiency trade-off of DFT and empirical potentials, enabling large-scale, high-fidelity simulations of complex materials at tractable computational cost [35, 36, 37]. For the W-H system, Wang et al. developed a deep potential (DP)-based model (denoted as DP-WH hereafter) [38] that successfully reproduces the pure W properties and solute H behavior. As the training dataset did not include configurations involving multiple vacancies decorated with numerous H atoms, the model’s transferability to H bubble dynamics remains uncertain [38].

In this study, we develop a unified MLP model for the W-H binary system based on the neuroevolution potential (NEP) framework (hereafter denoted as NEP-WH) [34], using a feedforward neural network (NN) trained against DFT data via an active learning scheme. The NEP framework is selected for its combination of near-DFT accuracy and exceptional computational efficiency [39], enabling simulations of more than two million atoms over time spans exceeding $10\text{\,}\mathrm{ns}$. Such capability is essential for resolving atomistic processes of hydrogen bubble-induced embrittlement. The developed NEP-WH model surpasses all previously reported W-H potentials, accurately reproducing fundamental lattice and defect properties of W, solute H behavior in both perfect and defect matrices, and the formation energy of the H${}_{\textrm{2}}$. Crucially, it provides substantial improvements in describing void-hydrogen interactions. Leveraging this model, we perform extensive MD simulations across nanovoids of varying radii, which reveal the formation of {100} planar hydrogen clusters and H-rich hexagonal close-packed (HCP) regions driven by planar hydrogen agglomeration. This evolution gives rise to H concentration-dependent a ductile-to-brittle transition under uniaxial tension.

RESULTS

NEP-WH model training

We employed an active learning strategy based on committee uncertainty estimates [42, 39] from an ensemble of NEP models to systematically explore the configuration space for training the NEP-WH model. The initial dataset comprised diverse configurations, including bulk W with and without interstitial H, H-decorated vacancy clusters, manually perturbed structures, and structures sampled from ab initio molecular dynamics (AIMD) simulations. Single-point DFT calculations using the Perdew-Wang (PW91) functional [43] (see Methods for details) were performed to obtain total energies and atomic forces for these structures, which were then used to train an ensemble of four NEP models. One model from this ensemble was then used to carry out MD simulations to further sample the configuration space. Model uncertainty was estimated from the variance in predicted forces across the ensemble, and configurations with large discrepancies were identified as outliers. These outliers were recalculated with DFT calculations to obtain energies, forces, and virial tensors, and subsequently added to the training set for the next iteration. The iterative process continued until the prediction errors across all sampled configurations converged within acceptable limits. The active learning procedure spanned a wide range of thermodynamic conditions, with temperatures from $50\text{\,}\mathrm{K}$ to $5100\text{\,}\mathrm{K}$ and pressures from $-2\text{\,}\mathrm{MPa}$ to $20\text{\,}\mathrm{GPa}$. To further enhance robustness, we manually incorporated configurations of melted and strained structures. The final training set included pristine, defective, liquid, and short-range interaction configurations for both W and W-H systems as well as $\mathrm{H}_{2}$ molecule, totaling 57,400 structures (8,267,784 atoms) (see Figure 1a). Details of the iterative process for different structural types are provided in Supplementary Note S1 and Table S1.

In each iteration, the fourth-generation of NEP framework [34, 39] was employed, with a radial cutoff of $6\text{\,}\mathrm{\text{Å}}$ and an angular cutoff of $5\text{\,}\mathrm{\text{Å}}$. To prevent unphysical clustering due to insufficient short-range repulsion, the Ziegler-Biersack-Littmark (ZBL) potential [44, 45] was integrated into the trained NEP-WH model. The input file used for training the NEP-WH model is provided in Supplementary Note S3, and additional details are described in the Methods section.

The parity plots for energy and force demonstrate the high accuracy of the NEP-WH model (Figure 1b-c). Despite the wide ranges in both quantities, the achieved root-mean-square errors (RMSEs) are $4\text{\,}\mathrm{meV}\text{\,}{\mathrm{\text{atom}}}^{-1}$ for energy and $122\text{\,}\mathrm{meV}\text{\,}{\mathrm{\text{Å}}}^{-1}$ for force. In addition to its accuracy, the NEP-WH model, as implemented in GPUMD [46], exhibits excellent computational performance. It achieves a computational speed of $4.5\times 10^{6}$ atom step s*-1*, approximately two orders of magnitude faster than the DP-WH model [38] implemented in LAMMPS (version 2Aug2023) [47], and only several-fold slower than the EAM potential by Wang et al. (denoted as EAM-Wang hereafter) [26], also implemented in LAMMPS (version 29Aug2024) [47] (Figure 1(d)). This performance benchmark was obtained by running MD simulations in the micro-canonical ensemble for 1,000 steps on bulk BCC W containing 18,522 atoms, using a single NVIDIA V100 GPU with 32 GB of memory. Notably, this is the maximum system size that the DP-WH model can handle under these conditions without encountering memory errors. In contrast, the NEP-WH model can scale up to five million atoms on the same hardware (see Supplementary Fig. S1). The superior computational efficiency of NEP-WH is crucial for subsequent hydrogen bubble simulations, which require extensive sampling to insert H${}_{\textrm{2}}$ molecules until steady-state bubble pressure is reached (see Figure 4(a)).

Evaluation of the basic properties of W-H systems

To assess the reliability of the NEP-WH model beyond RMSE metrics, we compare its predictions for a broad range of static properties of W-H systems with those from other models, including BOP models [48, 49], EAM models [26, 50], and the DP-WH model [38], as well as with previous experimental and DFT results, and additional DFT calculations conducted in this work (Table 1). The evaluated properties span elemental W and H, as well as binary W-H system. We performed static calculation for mainly using the NEP-WH model, and also other models when the relevant data are missing from the literature. The definitions and calculation details of all benchmark properties are provided in Supplementary Note S4.

For elemental W, all models reasonably predict lattice parameters and elastic constants. The formation energy of interstitial dumbbells ($E_{\rm f}^{\rm SIA}$) is correctly ordered only by the BOP-Li [49] and NEP-WH models. Monovacancy formation energies ($E_{\rm f}^{\rm vac}$) and vacancy migration barriers ($E_{\rm m}^{\rm vac}$) are well captured by all models except BOP-Juslin [48], which significantly underestimates $E_{\rm f}^{\rm vac}$. Among all models, only DP-WH and NEP-WH correctly predict the signs of the binding energies for both first- ($E_{\rm b}^{\rm divac-1NN}$) and second-nearest neighbor divacancies ($E_{\rm b}^{\rm divac-2NN}$). For surface formation energies ($E_{\rm f}^{\rm surf}$), only DP-WH and NEP-WH reproduce the reconstruction of the (100) free surface and both the correct magnitude and relative ordering for low-index facets, while NEP-WH additionally shows strong agreement with DFT for high-index surfaces, outperforming DP-WH (Supplementary Figure S3). Similarly, for grain boundary formation energies across various types, both MLPs perform significantly better than empirical potentials (Supplementary Fig. S4). To further evaluate the performance of the various models in reproducing static properties of elemental W, we calculated formation energies of interstitial dislocation loops, vacancy clusters, and the generalized stacking fault (GSF) energies along the $\left\langle 111\right\rangle$ direction of the (1 $\overline{1}$ 2) plane and (1 $\overline{1}$ 0) plane (Supplementary Figure S5). Both NEP-WH and DP-WH exhibit good accuracy, while the BOP and EAM models are less accurate.

In the binary W-H system, all models reasonably reproduce H solution energies at both TISs ($E_{\rm s}^{\rm tis}$) and OISs ($E_{\rm s}^{\rm ois}$) in BCC W. They also correctly identify the minimum-energy diffusion path as hops between adjacent TISs, with predicted migration barriers ($E_{\rm m}^{t\rightarrow t}$) in good agreement with DFT and experimental results. Furthermore, the solution energy of hydrogen at a TIS in hydrostatically strained BCC W lattice, and the migration barrier along the minimum energy path as a function of hydrostatical lattice strain (Supplementary Figure S5 f-g) calculated by both NEP-WH and DP-WH exhibit good accuracy. Recent DFT studies reveal that hydrogen prefers energetically to form planar clustered structures, where a planar hydrogen cluster comprising TIS H or OIS H aligns along a (001) plane, rather than existing as individual interstitials. Additionally, the rock-salt mono-hydride structure, with hydrogen occupying all [001] type OIS, is more energetically stable than both TIS and OIS planar self-clusters [70, 25]. Both the DP-WH and NEP-WH models successfully reproduce the DFT-predicted energetic trend: $E_{\rm b}^{\rm rock\text{-}salt}<E_{\rm b}^{\rm Tplanar}<E_{\rm b}^{\rm Oplanar}<E_{\rm s}^{\rm tis}<E_{\rm s}^{\rm ois}$, while all empirical potential fail to capture this order (Supplementary Fig. S6). Similarly, for the binding energy of a hydrogen atom at the most stable adsorption site on low-index surfaces ($E_{\rm b}^{\rm surf}$), MLP predictions closely match DFT reference values, while all empirical potentials exhibit either qualitatively incorrect trends or large quantitative deviations. Finally, for the H2 molecule, all models well accurately reproduce the formation energy $E_{\rm f}^{\rm H_{2}}$ and bond length $L_{\rm b}^{\rm H_{2}}$ of the H2 molecule in vacuum.

Across the benchmarks above, both the NEP-WH and DP-WH models clearly outperform all considered empirical potentials in capturing W properties, and W-H interactions in vacancy-free environments, demonstrating strong consistency with DFT reference values. However, this consistency breaks down in the presence of hydrogen-vacancy complexes. The DP-WH model has been shown to inaccurately predict the maximum number of H atoms that can be accommodated in monovacancies due to an insufficient diverse training set [38], raising concerns on its reliability for simulating hydrogen bubble nucleation and growth [15, 6]. To further assess these limitations, we extended our benchmarks to evaluate hydrogen trapping energetics in larger vacancy clusters ($\rm V_{4}H_{n}$, $n>0$) at finite temperatures. The $\rm V_{4}H_{n}$ structures were obtained by sequentially introducing hydrogen atoms at $600\text{\,}\mathrm{K}$, following the protocol described in Ref. [6] (see Supplementary Note S3 for details). While NEP-WH continues to deliver accurate predictions, DP-WH fails to reproduce correct hydrogen trapping energies when H2 molecules form in V4 clusters (Figure 2). Other empirical potentials exhibit even larger deviations at all trapping stages. Currently, NEP-WH is the only W-H potential capable of reliably describing both hydrogen trapping and $\rm H_{2}$ nucleation in nanovoids.

All properties discussed above were evaluated as static, zero-pressure configurations to establish the fundamental material characteristics. However, experimental observations have estimated hydrogen bubble pressures to reach approximately $10\text{\,}\mathrm{GPa}$ at $300\text{\,}\mathrm{K}$ in H${}_{\textrm{2}}$ bubbles [71, 72], underscoring the importance of accounting for pressure effects at room temperature. Figure 3 presents the H-H RDFs at $300\text{\,}\mathrm{K}$ and $5\text{\,}\mathrm{GPa}$, computed from MD simulations using various potential models and compared against DFT-driven MD results. The EAM model and DP-WH fail to reproduce the average intramolecular H-H distance associated with the first RDF peak, while the BOP model fails to capture the average intermolecular H-H distance reflected in the second peak. In contrast, NEP-WH model accurately captures both features, making it the only potential capable of reliably simulating bubble behaviour under intrinsic H2 molecule gas pressure. Its predictive accuracy is further validated at higher pressures of $10\text{\,}\mathrm{GPa}$ and $20\text{\,}\mathrm{GPa}$, as shown in Supplementary Figure S7.

Hydrogen bubble formation and evolution

After validating the reliability of the developed NEP-WH in capturing the lattice and defect properties of both W and W-H systems, and intrinsic H${}_{\textrm{2}}$ bubble behavior, we investigated the behavior of hydrogen dissolution in W nanovoids of varying sizes and the size dependence of the resulting hydrogen-bubble pressure. The dissolution process was modeled via an iterative procedure, where H${}_{\textrm{2}}$ molecules were sequentially inserted into the void core followed by structural relaxation at finite temperature (see Methods for details), effectively capturing the dynamic evolution. The bubble pressure was calculated as the hydrostatic pressure in the core sphere with a radius of $2\text{\,}\mathrm{\text{Å}}$ smaller than the nanovoid radius.

Figure 4a shows the evolution of hydrogen bubble pressure in nanovoids of varying radii over simulation time. For all nanovoids, H${}_{\textrm{2}}$ molecules initially fill the voids rapidly, resulting in a swift increase in pressure. During this process, some H${}_{\textrm{2}}$ molecules dissociate into H atoms that adsorb onto the nanovoid surface, while others remain trapped inside (see Supplementary Movie 1). A representative atomic snapshot of the void region at $0.109\text{\,}\mathrm{ns}$ for a $1\text{\,}\mathrm{nm}$ nanovoid is shown in Figure 4b, consistent with the stable dissolution configuration previously predicted by DFT calculations [6]. As additional H${}_{\textrm{2}}$ molecules are introduced, the system reaches equilibrium, maintaining a steady population of H${}_{\textrm{2}}$ molecules within the nanovoid and consequently stabilizing the hydrogen bubble pressure. Due to the higher solubility of hydrogen atoms in the tungsten matrix compared to helium [73], and the limited void sizes considered [6], we did not observe the continuous pressure buildup reported for helium bubbles that leads to loop punching [74]. Instead, excess hydrogen molecules dissociate at the void surface and diffuse into the surrounding matrix, resulting in a dynamic equilibrium of hydrogen bubble pressure within the void.

The equilibrium hydrogen bubble pressures, averaged over the final $1\text{\,}\mathrm{ns}$, are $18.73\text{\,}\mathrm{GPa}$, $17.58\text{\,}\mathrm{GPa}$, $16.95\text{\,}\mathrm{GPa}$, and $16.01\text{\,}\mathrm{GPa}$ for nanovoids with radii of $0.75\text{\,}\mathrm{nm}$, $1\text{\,}\mathrm{nm}$, $1.25\text{\,}\mathrm{nm}$, and $1.5\text{\,}\mathrm{nm}$, respectively. For larger nanovoids with radii of $1.75\text{\,}\mathrm{nm}$-$5\text{\,}\mathrm{nm}$, hydrogen atoms were initially placed inside the voids and iterative procedure was applied to accelerate the simulation. The resulting hydrogen bubble pressure exhibits an inverse proportionality to the void radius, consistent with the functional form predicted by established empirical model [6] (see Supplementary Note S5 and Supplementary Figure S8). While quantitative differences are present, the predicted pressure range qualitatively agrees with available experimental values [71, 72]. These deviations may originate from neglecting temperature and void surface energy effects. According to the Young-Laplace theorem [75], hydrogen bubble pressure in nanovoids should be proportional to $\gamma/R$, where $\gamma$ is the surface energy and $R$ is the void radius. This relationship has long lacked direct experimental validation. Our study provides quantitative data of hydrogen pressures in nanovoids of varying sizes, demonstrating a strong linear relationship between pressure and $\gamma/R$ (see inset in Figure 4a). This predictive relationship establishes a fundamental reference for future investigations of hydrogen bubble thermodynamics.

As the hydrogen bubble approaches equilibrium, additional H${}_{\textrm{2}}$ molecules dissociate into H atoms that diffuse into the surrounding matrix. Notably, hydrogen diffusion around all studied voids is anisotropic, with a clear preference for distribution along the {100} crystallographic planes. Figure 4c presents an MD snapshot at $1.417\text{\,}\mathrm{ns}$, illustrating the early-stage formation of platelet-like self-clusters along the {100} planes surrounding a $1\text{\,}\mathrm{nm}$ void. Similar visualizations along the [010] and [001] directions are provided in Supplementary Figure S9. These planar hydrogen self-clusters stochastically form and grow along the (100), (010), and (001) planes. In each orientation, the hydrogen atoms are confined to a $4\text{\,}\mathrm{\text{Å}}$-thick slab spanning two adjacent tungsten atomic planes and the interstitial plane between them. During the simulations, hydrogen atoms frequently hop between these planes but rarely escape from the planar cluster. These {100} planar hydrogen clusters extend outward from the void surface, and the hydrogen atoms within them occupy TISs. Their presence increases the interplanar spacing between adjacent W layers and induces BCC-to-FCC transitions in small, localized regions of the lattice; a minor HCP phase also emerges at the intersections of the (100) and (001) planar hydrogen clusters.

With ongoing hydrogen diffusion, the planar clusters along {100} planes expand outward from the void surface. Although atomic positions of hydrogen fluctuate slightly, the cluster thickness and orientation remain unchanged. This expansion enhances the extent of the FCC phase transformation. However, these FCC domains are dynamically unstable, repeatedly forming and annihilating during the simulation. This instability likely arises from the severe lattice distortions required for the BCC –FCC transformation in W, which leads to the absence of the rock-salt structure in our MD simulations and is in agreement with prior results [38]. At the same time, the intersections of planar hydrogen clusters along the {100} planes promote the formation of HCP phases that act as new hydrogen traps (see Figure 4d). Two wedge-shaped hydrogen-rich HCP regions gradually develop near the void, reaching local hydrogen concentrations of 1.3-1.5 H atoms per W atom. Both the planar clusters and the hydrogen-rich HCP regions continue to grow until the end of the simulation (Figure 4(d)).

Hydrogen self-clustering has previously been reported under conditions of high hydrogen concentration or in the presence of local stress fields [76, 25, 38]. Hydrogen-rich HCP regions have also been observed near cracks under strained Fe [23]. Here, we directly observe pronounced hydrogen aggregation along {100} planes and hydrogen-rich HCP regions near bubble surfaces. We attribute this behavior to the persistent internal H${}_{\textrm{2}}$ bubble pressure and the H dissociation around bubble surface. This directional aggregation may significantly influence the mechanical response of W, underpinning the experimentally observed intragranular blisters that preferentially develop along the {100} planes [11, 14, 13]. In the following section, we investigate the mechanical response of these hydrogen-enriched configurations.

Uniaxial tension

To assess the impact of hydrogen bubble evolution on the mechanical properties of tungsten, we analyzed the uniaxial tensile response of representive configurations sampled from the hydrogen bubble simulation. Specifically, we examined four cases: the initial void (hydrogen concentration $C_{\rm H}$ = 0 at.%) and hydrogen bubbles with $C_{\rm H}$ = 1.45 at.% (Figure 4b), $C_{\rm H}$ = 2.62 at.% featuring {100} planar clusters (Figure 4c), and $C_{\rm H}$ = 4.96 at.% containing both planar clusters and hydrogen-rich HCP regions (Figure 4d). Uniaxial tensile simulations were conducted along the [100] orientation at $600\text{\,}\mathrm{K}$ with a strain rate of $1\text{\times}{10}^{8}\text{\,}{\mathrm{s}}^{-1}$. As shown in Figure 5, the pristine void without hydrogen exhibits the highest tensile strength of $27.46(1)\text{\,}\mathrm{GPa}$. Introducing hydrogen significantly reduces the strength, decreasing to $24.73(7)\text{\,}\mathrm{GPa}$ at $C_{\rm H}$ = 1.45 at.%, $21.44(4)\text{\,}\mathrm{GPa}$ at $C_{\rm H}$ = 2.62 at.%, and $17.73(6)\text{\,}\mathrm{GPa}$ at $C_{\rm H}$ = 4.96 at.%. Meanwhile, the fracture behavior transitions from ductile to brittle, with the ultimate strain dropping from $\sim$0.12 in the hydrogen-free void to $\sim$0.06 at the highest hydrogen concentration. These results highlight the pronounced embrittlement of tungsten induced by hydrogen bubbles.

Pristine voids without hydrogen exhibit distinct dislocation emission and slip behavior under tensile strain (Figure 6a-b). During yielding, a pie-shaped stress concentration forms around the void center, oriented perpendicular to the tensile axis, accompanied by partial BCC to FCC phase transitions at its tips (Supplementary Figure S10a). Simultaneously, two 1/2 $\left\langle 111\right\rangle$ dislocation loop fragments are emitted from the nanovoid. As the trailing ends of these dislocations exit the void, vacancy clusters form on the void surface (see Figure 6a and Supplementary Figure S10b). The cross-slip and recombination of dislocations initiate an additional hardening stage at a strain range of $\sim$0.0760 - $\sim$0.1002 (Figure 5) and leave behind some vacancies along the slip path (Supplementary Figure S10b-c). Continued dislocation emission leads to coalescence of these vacancy clusters, culminating in fracture (Figure 6b). The resulting fracture surface varies slightly among independent simulations (Supplementary Figure S11b-f), shaped by dislocation slip and vacancy distribution.

The presence of hydrogen significantly suppresses dislocation emission and promotes premature fracture. Under uniaxial tension, hydrogen bubbles undergo pronounced brittle cleavage (Figure 6c-h). As the hydrogen concentration increases, the fracture crack evolves from a symmetric (Figure 6d) to increasingly asymmetric (Figure 6f,h) about the vertical axis of nanovoid, which is perpendicular to the tensile direction. This hydrogen concentration-dependent crack morphology is closely linked to the formation and evolution of planar hydrogen clusters during tensile loading.

In the nanovoid containing 1.45 at.% hydrogen, no planar hydrogen clusters are present initially. However, under uniaxial tension, hydrogen atoms migrate toward pie-shaped stress concentration region and progressively assemble into planar clusters along the (100) crystallographic planes, as shown at a representative strain of 0.0534 prior to fracture (Figure 6c and Supplementary Figure S10e). As strain increases, the clusters expand and further separate adjacent tungsten layers. As a result, cracks nucleate at the nanovoid surface and propagate outward. Hydrogen atoms preferentially segregate to the crack surfaces, and some hydrogen molecules diffuse into the crack, collectively promoting brittle cleavage. The resulting fracture surface aligns along the (100) plane, exhibiting a symmetric fracture crack about the vertical axis of the nanovoid center (Figure 6d and Supplementary Figure S11h-i).

A more pronounced hydrogen-induced ductile-to-brittle transition is observed in nanovoids with higher hydrogen concentrations (Figure 6e-h). In both the 2.62 at.% and 4.96 at.% cases, hydrogen platelets are initially distributed along all three principal crystallographic directions (see Figure 4c-d and Supplementary Figure S9). Under uniaxial tension along the [100] direction, cracks preferentially propagate along the (100) hydrogen cluster planes rather than through the void center. The formation of hydrogen platelets generates localized stress concentrations (Figure 7a,c), which further accelerate brittle cleavage and reduce the ultimate tensile strength of the material (Figure 7b,d). During this process, hydrogen atoms dissolved from other hydrogen-rich regions—such as (001) and (010) planar clusters—migrate toward and accumulate along the newly formed crack surfaces (see Supplementary Movie 2 and Movie 3), further promoting fracture. Specifically, in the 2.62 at.% case, this leads to an asymmetric fracture crack about the vertical axis of the nanovoid center (Figure 6f and Supplementary Figure S11n-r).

Compared to the 2.62 at.% case, where only planar hydrogen clusters are present, the 4.96 at.% case exhibits brittle cleavage in voids where both hydrogen-rich HCP region and planar hydrogen clusters coexist (Figure 6g). In this case, the initial (100) hydrogen clusters are larger than those in the 2.62 at.% case (see Figure 4c-d), resulting in a lower tensile strength (Figure 5). Under uniaxial tension, the HCP phase gradually depletes, generating vacancies and passivating the crack in that region (Figure 6h). During the depletion of the HCP phase, $\left\langle 100\right\rangle$ dislocation fragments are emitted, which absorb hydrogen atoms from the hydrogen-rich HCP region and relieve local stress (Figure 7c,d). Notably, this dislocation emission is not associated with ductile behavior. Once new planar hydrogen clusters form within the hydrogen-rich regions, brittle cleavage is reinitiated. The interplay between dislocation emission and fracture along hydrogen platelets ultimately leads to a highly asymmetric fracture crack with respect to the vertical axis of the nanovoid center, with the crack front propagating predominantly downward, as shown in Figure 6(h) and Supplementary Figure S11t-x.

DISCUSSION In summary, we present a machine-learned NEP model for the W-H system, enabling accurate and efficient atomistic simulations of hydrogen trapping and bubble evolution in tungsten. The model reproduces DFT results across a wide spectrum of properties, including short-range interactions, lattice and defect characteristics, hydrogen energetics in both pristine and defected environments, and H${}_{\textrm{2}}$ properties under pressure. It also captures hydrogen trapping energies in nanovoids with near-DFT fidelity, even in cases where H${}_{\textrm{2}}$ molecules form inside, thereby enabling predictive simulations of hydrogen bubble nucleation and evolution dynamics. Such breadth of accuracy, unmatched by existing empirical or MLPs, establishes the NEP-WH model as a significant advance. Beyond accuracy, the model delivers exceptional computational efficiency, supporting million-atom simulations on a single desktop GPU. This unique combination of fidelity and scalability offers a powerful framework for investigating hydrogen bubble dynamics with near first-principles precision.

Our hydrogen bubble simulations driven by the NEP-WH model reveal the atomic-scale mechanisms of bubble formation and evolution in nanovoids. H${}_{\textrm{2}}$ molecules preferentially occupy the void core, while $\rm H$ accumulates on the void surface. H${}_{\textrm{2}}$ molecules aggregation induces a rapid pressure increase that equilibrates under continuous hydrogen injection. Excess hydrogen diffuses to form planar clusters at TISs along the {100} planes, and their intersections yield hydrogen-rich HCP regions at high concentrations.

Uniaxial tensile simulations reveal a ductile-to-brittle transition, driven by planar hydrogen clusters that suppress dislocation emission from void surface and promote brittle cleavage. The fracture response depends strongly on hydrogen concentration: (i) at low concentrations, symmetric fracture surfaces pass through the nanovoid center, as planar hydrogen clusters form dynamically at the crack front; (ii) at medium concentrations, pre-existing planar hydrogen clusters guide fracture along their planes, producing crack surfaces aligned with the cluster orientation; (iii) at high concentrations, the interplay between dislocation emission from H-rich HCP regions and fracture of planar hydrogen clusters generates a highly asymmetric fracture morphology, with the crack front propagating unidirectionally.

Recent experiments on tungsten exposed to low-energy hydrogen plasma have reported intragranular blisters preferentially developing on 100 planes [11, 14, 13], consistent with our MD observations. This correlation suggests that the experimentally observed crack-shaped blisters may originate from 100 planar hydrogen clusters surrounding voids. We propose that once the internal pressure of hydrogen bubbles exceeds a critical threshold, the combined effect of bubble pressure and planar cluster-induced stress concentration initiates cracks. These cracks subsequently propagate by incorporating nearby H, ultimately forming cavities observed experimentally. Together, these mechanistic insights establish a direct connection between atomistic bubble evolution and macroscopic blistering, providing guidance for the design of plasma-facing tungsten components in fusion energy applications.

METHODS

Density functional theory calculations

All DFT calculations were performed using vasp code [77, 78] with exchange-correlation described by the generalized gradient approximation proposed by Perdew-Wang (PW91) [43]. The Brillouin zone is sampled by the Monkhorst-Pack scheme with a grid spacing of $0.2\text{\,}{\mathrm{\text{Å}}}^{-1}$. The projector-augmented-wave method was applied, and a plane-wave cutoff energy of $500\text{\,}\mathrm{eV}$ was used. The 5$d$ and 6$s$ electrons of W and the 1$s$ electron of H were treated as valence electrons. Convergence thresholds were set to $1\text{\times}{10}^{-5}\text{\,}\mathrm{eV}$ and $1\text{\times}{10}^{-2}\text{\,}\mathrm{eV}\text{\,}{\mathrm{\text{Å}}}^{-1}$ respectively. The INCAR input file for vasp package is provided in Supplementary Note S2.

Machine learning potentials

We employ the NEP approach [34, 39] to construct the NEP-WH model. The NEP approach provides a framework for generating highly efficient MLPs and has been successfully applied to model structural, chemical, thermal, mechanical, and transport properties of complex materials [37]. Specifically, we adopt the latest NEP4 implementation, which has been applied to develop a general-purpose potential for many metals and their alloys [39].

The name NEP originates from its two key components. On the one hand, it is a NN potential model [30]; on the other hand, the free parameters in the model are trained using an evolutionary algorithm [34]. The input layer of the NN consists of a number of descriptors $\mathbf{q}_{i}$, which represent the local chemical environments of a central atom $i$. There are two types of descriptors used in NEP4: radial descriptors $q_{n}^{i}$ that depend only on interatomic distances, and angular descriptors $q_{nl}^{i}$ that depend on both atom distances and bond angles. The index $n$ in $q_{n}^{i}$ ranges from 0 to $n_{\rm max}^{\rm R}$ and the index $n$ in $q_{nl}^{i}$ ranges from 0 to $n_{\rm max}^{\rm A}$. The index $l$ ranges from 0 to $l_{\rm max}$. To balance accuracy with efficiency, we choose $n_{\rm max}^{\rm R}=9$, $n_{\rm max}^{\rm A}=8$, and $l_{\rm max}=4$. The radial and angular descriptors are smoothly cut at $6\text{\,}\mathrm{\text{Å}}$ and $5\text{\,}\mathrm{\text{Å}}$, respectively. At short distances, the NEP model is smoothly connected to the ZBL potential [44] to account for the strong repulsive forces encountered under extreme conditions [45].

The NN model takes the descriptors as the input layer neurons and output the site energy $U_{i}$ of the atom $i$. With a hidden layer with $N_{\rm neu}$ neurons, the site energy can be explicitly written as

[TABLE]

Here, $\tanh(x)$ is the activation function, $w^{(0)}$ are the weight parameters connecting the input layer (with dimension $N_{\rm des}$) and the hidden layer, $w^{(1)}$ represents the weight parameters connecting the hidden layer and the output layer (the site energy), $b^{(0)}$ represent the bias parameters in the hidden layer, and $b^{(1)}$ is the bias parameter in the output layer. We use a relatively large value of $N_{\rm neu}=120$ in our model. All these parameters are optimized by minimizing a loss function that is a weighted sum of the RMSEs of energy, force, and virial. Details of the nep.in and zbl.in input files are provided in Supplementary Note S3.

Molecular dynamics

The NEP-WH model was used to conduct large-scale MD simulations to investigate hydrogen solubility in nanovoid and the influence of H on the mechanical properties of W. All simulations were performed using gpumd package [46], and visualizations were carried out with the ovito software [79].

We first constructed a cubic simulation cell with dimensions of $10\times 10\times 10$ nm3 along the [100], [010], and [001] directions. Periodic boundary conditions were applied in all three directions. Nanovoids with radii of $0.75\text{\,}\mathrm{nm}$, $1\text{\,}\mathrm{nm}$, $1.25\text{\,}\mathrm{nm}$, $1.5\text{\,}\mathrm{nm}$, $1.75\text{\,}\mathrm{nm}$ were created at cell center. For larger void radii of $2\text{\,}\mathrm{nm}$, $3\text{\,}\mathrm{nm}$, $4\text{\,}\mathrm{nm}$, and $5\text{\,}\mathrm{nm}$, the domain size was increased while maintaining a constant void volume fraction of 0.014, in order to minimize periodic boundary effects. An initial energy minimization was performed to relax the atomic configurations and eliminate any residual forces. The W system was then equilibrated for $20\text{\,}\mathrm{ps}$ in the isothermal–isobaric (NPT) ensemble at $600\text{\,}\mathrm{K}$, with a time step of $2.0\text{\,}\mathrm{fs}$. H${}_{\textrm{2}}$ molecules were randomly introduced into the void cores, with each pair of H atoms seperated by more than $1.3\text{\,}\mathrm{\text{Å}}$. The void core was defined as a sphere with a radius of $2\text{\,}\mathrm{\text{Å}}$ smaller than the actual void radius, to avoid interaction between H and the W surface. After every 10 H${}_{\textrm{2}}$ molecules were inserted, the system was relaxed for $1\text{\,}\mathrm{ps}$ in the canonical (NVT) ensemble, with a time step of $0.4\text{\,}\mathrm{fs}$. If an insertion attempt failed, it was skipped. Additionally, a $10\text{\,}\mathrm{ps}$ relaxation was performed after every 100 insertions to ensure full structural equilibration.

To study the effect of H on mechanical behavior of W, uniaxial tension was applied along the [100] direction for systems containing a $1\text{\,}\mathrm{ns}$ void with varying H concentrations, with a strain rate of $1\text{\times}{10}^{-8}\text{\,}{\mathrm{s}}^{-1}$. Simulations were carried out in the canonical (NVT) ensemble with a time steps of $0.4\text{\,}\mathrm{fs}$ W-H systems and $2\text{\,}\mathrm{fs}$ for pure W systems.

DATA AVAILABILITY

The NEP-WH potential files are freely available at the Zenodo repository https://doi.org/10.5281/zenodo.16973981. The corresponding reference datasets will be made available upon acceptance of this work. Additional data that support the findings of this study are available from the corresponding author upon reasonable request.

Author contributions P.Q. initiated and conceived the study. Y.B. and K.S. prepared the training and test structures. Y.B. performed the DFT calculations. Y.B. and J.L. tested the various hyperparameters and trained the NEP model. P.Y. analyzed the descriptor space. Y.B. and Y.N. evaluated the NEP model’s performance on basic properties. Y.B., K.S., and Y.W. performed MD simulations. P.Q. and P.Y. supervised the project. All authors discussed the results and contributed to the writing of the paper.

Declaration of competing interest The authors declare that they have no competing interests.

Acknowledgement This research was funded by the Advanced Materials-National Science and Technology Major Project (No.2024ZD0607800), the National Key Research and Development Program of China (2023YFB3506704), the National Natural Science Foundation of China (52371003), and the University of Science and Technology Beijing High Performance Computing university-level public platform. P.Y. is supported by the Israel Academy of Sciences and Humanities & Council for Higher Education Excellence Fellowship Program for International Postdoctoral Researchers. Y.W. is supported by the Academy of Finland’s Grant No. 353298 under the European Union – NextGenerationEU instrument.

References