Superior hydrogen permeation resistance via Ni – graphene nanocomposites: Insights from atomistic simulations

Designing hydrogen-resistant Ni-based alloys from the perspective of the Ni/graphene interface (NGI) provides the potential to increase hydrogen trapping away from potential fracture paths. Nonetheless, numerous essential mechanisms of hydrogen penetration behaviors in the Ni-graphene nanocomposites are presently not well understood. Here we investigate the influence of Ni/graphene interfaces (NGIs) on the behavior of hydrogen diffusion and trapping in their vicinity using atomistic simulations. Hydrogen diffusion is competitively affected by elevated temperatures and NGIs. The difference in the mean square displacement for hydrogen between the composites and pure Ni can be of two orders of magnitude, highlighting the sluggish diffusion in the composites. As NGIs reduce hydrogen formation energy and diffusion barrier, hydrogen prefers to migrate towards the interfaces. Hydrogen readily forms sp 3 C – H bonds with C atoms, thereby impeding its detachment from graphene and subsequent entry into a non-diffusible state. Results of the study will contribute to the use of Ni – graphene

Generally, there are two mainstream HE mechanisms proposed for high-strength metallic materials including Ni-based alloys, i.e., hydrogen-enhanced decohesion (HEDE) and hydrogen-enhanced localized plasticity (HELP) [14][15][16].The HEDE mechanism arises from that hydrogen reduces the cohesion strength along crystallographic planes, grain boundaries (GBs), or precipitate/matrix interfaces in metals, and thus degrades the fracture threshold of materials [16].The HELP mechanism postulates that the presence of hydrogen facilitates the nucleation and mobility of dislocations, ultimately rendering a localized softening and plastic failure of materials [17].Many studies have shown that the processes of hydrogen transport and trapping around a crack tip are critical to the activation of the specific HE mechanism in Ni-based alloys [18][19][20].The hydrogen states in alloys are typically divided into two categories: diffusible hydrogen, which can move and spread throughout the materials, and non-diffusible hydrogen, which is trapped and cannot easily move [21,22].Diffusible hydrogen is typically found in reversible traps, meaning it can be released and reabsorbed, while non-diffusible hydrogen is found in irreversible traps, meaning it cannot be easily released [21].Among them, the diffusible hydrogen is generally believed to be the primary factor contributing to HE fractures [23].Since it is the diffusible hydrogen that causes damage to Ni-based alloys, any approach capable of immobilizing it should effectively mitigate its deleterious effects [16].
To control diffusible hydrogen, all that is necessary is to introduce benign traps in metals or inhibit the ingress of hydrogen [23].Based on this principle, designing hydrogen-resistant materials from the standpoint of interface engineering opens the possibility of enhancing the hydrogen trapping away from potential fracture paths (e.g., trapping at heterointerfaces instead of dislocations) [23][24][25][26].More interestingly, it has recently been proposed that introducing graphene into metals is feasible for improving the diffusion resistance of hydrogen [17,27,28].This assumption can be supported by the following evidence.First, the unique atomic structure of graphene with a honeycomb-like pore size of only approximately 0.064 nm, is theoretically predicted to be able to block H atoms from penetrating across the carbon sheet [29].Second, H atoms could be adsorbed on the surface of graphene by forming sp 3 C-H bonds, inhibiting their diffusion [17,30].Finally, metal/graphene interfaces may also act as effective sinks for trapping H atoms, which can be inferred from previous irradiation-resistant studies [31][32][33].As a result, Ni-graphene nanocomposites (NGNC) may offer a viable solution to circumvent the HE dilemma of Ni-based alloys.For example, Fan et al. [17] reported experimentally that a graphene coating deposited onto metallic Ni can effectively protect the metal against HE.Zhou et al. [34,35] demonstrated that the addition of graphene into Ni coatings can further increase the impact of coating on the kinetics of hydrogen evolution reaction and hinder hydrogen permeation.Nevertheless, to the best of our knowledge, the study of enhancing mechanism of graphene and/or Ni/graphene interfaces (NGIs) on hydrogen permeation resistance in NGNC is still in its infancy, and even many fundamental mechanisms, especially atomic-level behaviors, have been not well understood at present.
To predict the hydrogen permeation resistance, a comprehensive understanding of the interactions between hydrogen and various traps such as GBs and heterointerfaces is necessary [24][25][26][36][37][38][39][40][41][42][43].As one of the most important approaches, atomistic simulations offer essential insights into comprehending the intricate atomic-level processes governing hydrogen diffusion behaviors in the vicinity of traps.For example, combined molecular statics (MS) and molecular dynamics (MD) simulations, Taketomi et al. [37] investigated the hydrogen distribution and diffusion around a {112}<111> edge dislocation in α-Fe, and they revealed that strong trap sites are widely distributed across the slip plane surrounding the dislocation core, resulting in an anisotropic diffusion behavior of hydrogen in this region.The MD observation of Foster et al. [38] indicated that tritium prefers to concentrate at various types of GBs in α-Zr and the concentration is approximately three times higher than that in bulk regions.The hydrogen segregation energy maps of Li et al. [39] based on MS calculations showed that the trapping capacity of hydrogen is highly contingent on the GB structure in fcc-Ni, as evidenced by the notable changes in the maximum excess hydrogen concentration with increasing misorientation angle.Zhou et al. [40] reported that fcc-Ni GBs can increase the in-plane hydrogen diffusion but reduce the out-of-plane diffusion due to a trapping effect.Foster et al. [38] demonstrated that planar GBs in hcp-Zr also can enhance tritium diffusivity on the boundary plane while exerting minimal influence on out-of-plane diffusion.Instead, Shim et al. [41] found that the hydrogen diffusivity at the Σ3 (and Σ5) tilt GBs of bcc-V is lower than that within the bulk in the temperature range between 473 and 1473 K, albeit the in-plane hydrogen diffusion is dominant at the Σ3 boundary with a relatively small displacement in the direction normal to the boundary plane.Moreover, Islam et al. [42] demonstrated a significant three orders of magnitude reduction in the diffusion coefficient of H atoms trapped at α-Fe/Fe 3 C interfaces compared to the α-Fe bulk, and they observed the absence of H atoms crossing the interface from either of the phases.Hodille et al. [43] showed that Be/BeO interfaces serve as effective trapping sites for deuterium atoms.Additionally, they identified the presence of additional kinetic barriers, such as deuterium trapping sites and deuterium molecule formation/dissociation in BeO, which significantly hinder the attainment of equilibrium.Hence, gaining a comprehensive understanding of how graphene and/or NGIs influence the diffusion and trapping of hydrogen in their proximity through atomistic simulations holds fundamental significance in the development of NGNCs that can effectively resist HE fractures in nuclear reactors.
In this study, MD simulations were performed to investigate the resistance of NGNC to hydrogen penetration, focusing on the observation of the behavior of hydrogen diffusion and trapping at early stages around an NGI.Noted that pure Ni was chosen as the model material for the matrix of NGNCs, which is consistent with previous studies in terms of simplifying the system [32,34,35,38,39].First, thermodynamic simulations were carried out at temperatures ranging from 500 to 1200 K, which allowed observation of the dynamic evolution mechanisms of hydrogen and quantifying the diffusion parameters.Subsequently, the hydrogen diffusion and trapping mechanisms of NGNC at the atomic scale were examined from both energetic and kinetic perspectives.The study shows that NGIs can reduce the diffusivity of hydrogen in their vicinity and promote hydrogen trapping in interfacial regions, resulting in increased durability of NGNC as a hydrogen-resistant material for nuclear reactors.

Simulation methodology
A sandwich-like NGI model with top-fcc configuration was first generated in this work, as shown in Fig. 1, and the creation and optimization details of the model have been described elsewhere [31,32].Periodic boundary conditions were imposed in all the coordinate directions.The NGI system comprises 32,500 atoms with a size of 6.21 × 5.59 × 10.10 nm 3 .The embedded atom method (EAM) potential parameterized by Angelo et al. [44] for Ni-H systems (including Ni-Ni, H-H, and Ni-H interactions) has been widely accepted [39,40], and it was also used in this work.The EAM potential can appropriately describe the dissolution and diffusion of hydrogen and accurately reproduce the migration and formation energies of hydrogen in fcc-Ni.The C-C and C-H interactions were described by the adaptive intermolecular reactive empirical bond order (AIREBO) potential [45].The Lennard-Jones potential was used to describe the Ni-C interactions [31,32].All simulations were implemented utilizing MD simulator LAMMPS code [46], and atomic visualizations were achieved with OVITO software [47].
Subsequently, a solute H atom was introduced into the system by inserting it at a random interstitial site at 2.61 nm away from the graphene (see Fig. 1).After energy minimization, the NGI-H system was thermally relaxed at a certain temperature for 5.0 ns with a timestep of 1.0 fs?The Nose-Hoover isobaric-isothermal ensemble (NPT) was employed to ensure no significant fluctuation of system pressure.The thermostat temperature was set within the range varying from 500 to 1200 K (with an interval of 100 K) to elucidate the hydrogen diffusion properties dependent on temperature.Because of the negligible effect of quantum tunneling on hydrogen diffusion at these temperatures [48], classical MD simulations can be justified in this work.During thermal relaxation, the positions of H atoms at any time t were recorded.Thus, the mean square displacement (MSD) of H atoms is calculated from the total diffusion distance according to the following equation where r i (t 0 ) and r i (t + t 0 ) represent the initial reference position of ith atom at time t 0 and its instantaneous position at time t 0 +t, respectively, N is the total number of atoms, and 〈…〉 indicates the ensemble average.Based on the Einstein relation, the self-diffusivity of hydrogen in NGNC for a given temperature can be quantified as shown below where d is the dimensionality of the system and describes one-, two-, or three-dimensional diffusion of H interstitials when d = 1, 2, or 3, respectively.To accurately calculate the diffusion coefficient, a "trajectory time decomposition (TTD)'' technique [49,50] was adopted, in which the single and long trajectory was decomposed into a series of shorter independent segments with equal duration, and then each segment can be calculated D i (i denotes ith segment) based on Eq. ( 2).The time interval of each segment was determined to be 50.0ps, followed by averaging D i over all segments.After obtaining the diffusion coefficients at different temperatures, the diffusion activation energy (E a ) of H atoms can be deduced from the Arrhenius equation where D 0 and k β represent the pre-exponential factor and the Boltzmann constant, respectively.Additionally, to serve as a control, the simulation settings and hydrogen diffusion calculations for pure Ni were maintained consistent with those applied to the NGNC.

Results
At a given temperature, a solute H atom may leave its original position to walk randomly, and during the atomic-scale dynamic process, the diffusion trajectory can be tracked and logged.Fig. 2 shows the diffusion trajectories of single H interstitial near the NGI at different temperatures, in which displacement vectors between the atomic positions at 0 and 5.0 ns are plotted per 10.0 ps and connected one by one.At 500 K, it is so difficult for the solute H to migrate each step that it cannot move to a further position.Suddenly, this restraining behavior began to be lifted at 600 K, and the solute H can achieve a long-distance migration.Since then, with the increasing temperature, both the diffusion frequency (i.e., the reciprocal of dwelling time) and diffusion distance of hydrogen increase accordingly.As the temperature exceeds 900 K, the hydrogen diffusion behaviors change again: the individual diffusion distance increases further while the total diffusion distance and frequency gradually reduce.In this situation, the solute H can cross the box boundary and reach the other side of the bulk.Above 900 K, the solute H for each case will eventually be trapped in the NGI.Especially, the solute H entered the NGI without being disturbed too much at 1200 K.In addition, the diffusion trajectory of a solute H atom initially located at 0.60 nm from the graphene is also displayed in Fig. S1.Evidently, the solute H at each temperature shows the tendency of directly migrating toward the interface.Even when the diffusion trajectory is changed by the temperature-induced stochastic lattice disturbance (e.g., at 700 or 1000 K), the solute H is still quickly trapped by the interface.These suggest that the closer a solute H is to the interface, the more it will be affected by the interface.Nevertheless, no H atoms are observed to penetrate the graphene of NGI during the trapping processes at different temperatures, unlike the case that pure graphene membrane allows thermal protons to permeate it [51].This agrees with the speculation about the synergistic effects of sp 3 C-H bonds and NGI sink character in inhibiting hydrogen diffusion.To further highlight the effects of NGI, the hydrogen diffusion trajectory in pure Ni is shown in Fig. S2 for comparison.Different from the NGNC, an increase in temperature will always intensify both hydrogen diffusion frequency and diffusion distance in pure Ni, suggesting that H atoms in pure Ni are difficult to fix and thus evolve into diffusible hydrogen.This may give preliminary evidence for the excellent resistance ability of NGNC to HE fractures.
The diffusion behavior of solute H near the NGI can be further identified by analysis of the above computer-generated trajectories.As a benchmark reference, the details of H trajectories in pure Ni are first studied and the corresponding MSDs are displayed in Fig. 3(a).The solute H exhibits random walk all the time, which makes the MSDs difficult to stabilize for long periods of time.The MSDs increase generally with increasing time at different temperatures but they exist significant fluctuations, which are also observed in other similar studies [52].As the temperature rises, the peak value of MSD augments first, reaching its ceiling at 1100 K, and then reduces, implying that the diffusion of solute H in a perfect fcc crystal is limited by the level of its neighboring stochastic lattice disturbance due to temperature.In Fig. 3 (b), the MSDs of the single H interstitial at 2.61 nm away from the NGI at different temperatures are exhibited.At 500 K, several steps appear in the diffusion profile, but the MSD is always lower than 0.3 nm 2 with time.This suggests that although the solute H has jumped, the extent of the jump is not obvious, explaining the phenomena in Fig. 2(a).In addition, the diffusion profile looks very rough due to many burrs at 500 K, illustrating that the solute H significantly oscillates around each dwelling site but tends to next jump.With the temperature elevating, the number of solute H jumps gradually increases and its corresponding dwelling time decreases.However, when the temperature reaches 900 K, the solute H will no longer jump after a certain period but remain to dwell.After this, the final dwelling duration is positively associated with temperature, in other words, the total jump time gradually shortens.Especially, the hydrogen diffusion quickly becomes stable within 0.046 ns at 1200 K. Combined with Fig. 2, the long dwelling behavior stems from the fact that the solute H is trapped by the NGI and then may form a sp 3 C-H bond with a C atom on the surface of graphene, making it difficult for the solute H to diffuse further.This behavior differs markedly from that of large-angle GBs with a large excess of free volume, which can accelerate the diffusion of solute H along them, but resembles that of small-angle GBs with high-density dislocations and that of other the MSD rises to fluctuate distinctly around its platform value at 1200 K like that at 500 K, suggesting that the energy for the next jump of solute H is being gathered.It can be expected that if the temperature elevates further, the solute H should be able to gain enough energy to get rid of the bound of the C-H bond and continue to migrate along the NGI [54].
In fact, this phenomenon is indeed observed above 1800 K but has not been shown herein, since it is beyond the scope of this work.The trend of the peak value of MSD with temperature is roughly in line with that in pure Ni, while its ceiling occurs at 900 K, a lower temperature than that of pure Ni, which may be because the sink role of NGI in trapping solute H starts to work at the temperature.At the same temperature, the MSD for NGNC is significantly less than that of pure Ni, even in the case that the solute H is not yet captured by the NGI below 900 K.This may be due to the interfacial attraction inhibiting the in-plane hydrogen diffusion in the bulk while enhancing the out-of-plane diffusion, like a magnet effect.The MSD difference between NGNC and pure Ni can be of two orders of magnitude, highlighting the sluggish diffusion of solute H in NGNC.
A single H interstitial often migrates three-dimensionally along different crystallographic directions.Hence, the MSD can be further decomposed along the x, y, and z directions, and the three components are presented in Fig. 4. As for pure Ni, the solute H diffuses mainly along the [111] direction at 500 K, the [‾110] direction at 600 K, or the [‾1‾12] direction at 700 K, in the form of one-dimensionally forward or backward movement.The diffusion is changed to two-dimensional mode within the temperature range from 800 to 900 K.A similar diffusion pattern for single He interstitials has also been found in α-Fe [49,50].However, the main contributions to the MSDs return to one-dimensional diffusion above 900 K. Changes in the diffusion pattern of H-interstitials in bulk Ni with temperature may be synergistically regulated by their forward-backward motion mode and by their adjacent lattice perturbations.Overall, the motion of the solute H in pure Ni is preferred along the [111] or [‾110] direction, which may be because it has a relatively lower activation energy barrier in this direction [49,55,56].This can also be confirmed in the subsequent results.As for NGNC, the diffusion behavior of solute H below 700 K almost remains consistent with that in pure Ni, i.e., switching from the [111] direction at 500 K to the [‾110] direction at 600 K.As the temperature rises to the range of 700-900 K, the solute H mainly diffuses two-dimensionally along the xy-crystallographic plane.These suggest that the diffusion pattern of solute H is not much affected by the interface below 900 K.However, at any temperature above 900 K, the diffusion gradually shifts from a two-dimensional to a one-dimensional mode with time, like the behavior of single He interstitials near a symmetric tilt GB in bcc W at 1000 K [57], and the higher the temperature, the faster the mode shift.Unlike below 700 K, one-dimensional diffusion along the [111] direction is dominated above 900 K.The shift in diffusion mode may originate from the following reasons.The previous studies [58,59] have demonstrated that an interface has a narrower range of influence on defects in its vicinity at low temperatures than at high temperatures.Below 700 K, the NGI hardly affects the hydrogen diffusion because the initial NGI− H distance is likely greater than the range of NGI influence on a solute H, and thus the solute H can retain its original diffusion behavior like that in pure Ni.However, with the temperature elevating, the vibrations of lattice atoms intensify, causing the gradual expansion of interface width and the concomitant increase of the range of NGI influence.On the other hand, the increase in the individual diffusion distance of hydrogen at high temperatures raises the probability of hydrogen diffusion towards the widening interfacial region.Consequently, after a period of two-dimensional walk, H atoms tend to move toward the NGI and are then rapidly captured by the interface at high temperatures, especially above 1200 K.
The MSD analysis described above can be further used to extract the solute H diffusion coefficients in pure Ni and the NGNC system at different temperatures, and then the diffusion activation energies based on the Arrhenius equation, all summarized in Fig. 5.In the pure Ni system, all MD data points (see Fig. 5 The prefactor and activation energy of hydrogen diffusion in pure Ni is calculated to be 1.880 nm 2 ps − 1 and 0.451 eV, respectively, which are approximately comparable with the results (i.e., 3.835 nm 2 ps − 1 and 0.508 eV) obtained by Zhou et al. [40] using the same EAM potential.These can also confirm the rationality of the calculation method used in this work.In the NGNC system, a similar linear law also appears in the MD data points below 900 K, while the Arrhenius equation is not satisfied with them above 900 K, which shows a gradual downtrend (see Fig. 5(b)).In particular, this nonlinear behavior above 900 K in the NGNC system is reflected in all dimensions, due to the remarkable interference of the interface on the natural diffusion of solute H.The hydrogen diffusivity at 1200 K is even close to that at 500 K, meaning that once H atoms reach the interface, there is almost no diffusion out of it.This behavior is distinct from that of small-angle GBs, where the in-plane diffusion of solute H is still possible although the diffusivity is reduced due to high-density dislocations [40,41,52].Thus, the MD data points displayed in Fig. 5(b) can be divided into the linear Arrhenius regime at low temperatures and the nonlinear regime at high temperatures to investigate.Within the linear Arrhenius regime, the prefactors of hydrogen diffusion in all dimensions decrease concerning those in pure Ni, especially in the [111] direction by 0.81 nm 2 ps − 1 .Furthermore, compared to pure Ni, the corresponding activation energies also tend to decrease in NGNC, while an opposite phenomenon occurs in the [‾1‾12] direction, improving by 0.017 eV.These changes in hydrogen diffusion parameters due to NGIs are approximately consistent with those due to fcc-Ni GBs found by Zhou et al. [40], highlighting the role of the interface in attempting to influence the diffusion behavior of nearby solute H at low temperatures.On the other hand, either in the NGNC or in pure Ni system, the prefactor of hydrogen diffusion shows the largest value in the [‾1‾12] direction, followed by that in the [111] direction, and the smallest value in the [‾110] direction.There is a difference of about 3.5 (or 2.2) times in the prefactor of hydrogen diffusion between the [‾1‾12] and [‾110] directions in the NGNC (or pure Ni) system.This implies that the [‾1‾12] direction facilitates faster hydrogen transport compared to the [111] and [‾110] directions due to the least dense arrangement of atoms among them.Therefore, the diffusion of solute H may be synergistically regulated by the atomic density of the bulk lattice and the trapping character of NGIs before being completely captured by the interface.Within the nonlinear regime, the hydrogen diffusion coefficients in all dimensions exhibit much lower values than those in pure Ni, especially at 1200 K.In addition, because the solute H becomes swiftly trapped and exhibits limited movement, the diffusion activation energies for H atoms within the NGI at temperatures exceeding 900 K fail to accurately describe the actual diffusion barrier.Taken together, the presence of NGIs in the composites can reduce the diffusivity of hydrogen in their vicinity compared to single crystal Ni, especially at high temperatures, resulting in slower hydrogen diffusion and thus easier trapping by the interfaces.

Discussion
The presence of NGIs leads to significant changes in the hydrogen diffusion behavior, especially causing the migration of solute H toward the interfaces.An increase in temperature can accelerate this process.Moreover, once solute H atoms enter the NGIs, they will be instantaneously adsorbed on the surface of the graphene.Unless stimulated by an extremely high temperature (exceeding 1800 K), it is difficult for the solute H to pull against the restraint of graphene and migrate further along the interfaces.To uncover the underlying physical mechanisms driving these phenomena, it becomes essential to analyze how the diffusion and trapping behaviors of the solute H atoms at different positions around the NGIs are affected by their local environments.This analysis should be approached from both energetic and kinetic perspectives.

Energetics of hydrogen diffusion and trapping
In a perfect fcc crystallographic texture, there are two most possible interstitial configurations for solute H: octahedral and tetrahedral interstitial sites (OIS and TIS) [36,60].The formation energy of solute H in a Ni matrix can be calculated by where E H t and E t represent the total energy of the observed system with and without the solute H, and 1  2 E H2 denotes the chemical potential of a single H atom separated from the dissociation of the H 2 molecule at 0 K, respectively.Here the hydrogen-chemical potential was determined from the EAM potential to be − 2.37 eV [61], close to the experimental value of − 2.26 eV [62].Further, the capture efficiency of a given site α of NGIs on H atoms from the bulk can be quantified by the segregation energy, defined as where is the formation energy of the solute H in the bulk or interface region of the NGI system, respectively.Utilizing these provided definitions, formation energies displaying negativity signify the thermodynamic preference for the dissociation and incorporation of hydrogen within the Ni matrix.Conversely, segregation energies characterized by positivity suggest a propensity for hydrogen to preferentially segregate at the NGI.
The formation energy of solute H related to its distance from the NGI is shown in Fig. 6 (a), in which the data at each distance was sampled from 20 different interstitial sites.The formation energy of an isolated H atom at the OIS and TIS in the bulk is 0.189 and 0.60 eV, respectively, which is approximately in agreement with the values (i.e., 0.17 and 0.58 eV) of Huang et al. [61] using the same EAM potential.This result also supports the consensus that H atoms remain more stable at the OIS than at the TIS in a Ni matrix [36,60], while exhibiting opposite characteristics relative to He atoms [63].Therefore, the H atoms at the OIS have been focused here and used for further calculations of the energetics.Intuitively, the energetics of solute H is mirror symmetry centered on the graphene.When the solute H approaches the interface, the formation energy deviates significantly from its value in the bulk and presents a gradual decrease followed by a sharp increase, and its minimum value is − 2.297 eV, which means that the solute H is more stable in the middle between the graphene and the Ni layer within the interface than in the region closer to the graphene or closer to the Ni layer.Conversely, the segregation energy of the solute H near the NGI shows a tendency to increase and then decrease, but always remains positive and reaches its maximum at 2.484 eV.The endothermic characteristic indicates an attractive interaction between the NGI and the solute H near it.One very striking thing is that the minimum value of segregation energy (i.e., 0.247 eV) near the interface is even comparable to the maximum value (i.e., 0.25 eV) near the Ni GBs [60], highlighting the trapping capability of NGI to solute H.The range of effective interaction of NGI with an isolated H atom reaches 15.5 Å, which is less than the initial NGI− H distance in the MD simulations but larger than that of the interface interacting with an isolated He atom [31].In fact, according to the theory of hydrogen segregation energy [64], it is also possible to give the hydrogen occupation probability at a particular site α within the NGI, derived as where θ bulk denotes the hydrogen concentration in atomic ratio in the bulk.Based on the above calculations of energetics, once the hydrogen concentration of bulk in actual scenarios is identified, the content of solute H atoms that are trapped by the NGI can be predicted quantitatively, and vice versa.Taking the temperature of 1200 K as an example, the hydrogen concentration of bulk is approximately 0.344 appm in the above MD simulations, and thus the θ α can be determined to 1.27 × 10 − 15 when E α seg = 2.484 eV.Considering that hydrogen segregation energy is not uniform near the NGI, this probability may also screen all favorable trapping sites of solute H within the NGI.Hence, the presence of solute H within the NGI is characterized by favorable energetics, while an elevated temperature milieu can expedite their migration towards this interface.Such a phenomenon contributes to mitigating the susceptibility of Ni-based alloys to embrittlement by employing the principles of "heterointerface engineering" [23].
In addition, different from the self-interstitial atoms as well as He atoms [31,32], solute H within the interface does not have any barrier-free layers.Instead, the gradual change in the formation energy suggests the presence of several metastable locations for solute H within the interface, affected by both local crystal-structural compositions and the surrounding stress conditions.To elucidate this matter, numerous H atoms were poured into the interfacial region randomly and then allowed to relax statically to their nearby stabilization sites based on the θ α .The stable sites of H atoms within the core of the NGI and the corresponding formation energies are shown in Fig. 6(b).Note that Fig. 6(b) displays only part of H atoms at the occupiable sites within the core of the NGI rather than the maximum excess, which is sufficient for unveiling the general hydrogen behaviors within the interface and beneficial to save computational resources.By identifying the difference in formation energy, it is easy to figure out that there are three different types of stable sites for H atoms within the core of the NGI, which are highlighted near the virtual hexagons in Fig. 6(b) and named sites S I , S II , and S III , respectively.Among them, all the S I sites are distributed in the center of each six-membered ring of graphene when viewed along the z-axis direction, and the S II and S III sites overlap each other at the C atoms.The formation energies of solute H atoms at the three types of stable sites are approximately − 0.232, − 0.348, and − 2.298 eV, and the corresponding segregation energies are 0.425, 0.542, and 2.484 eV, respectively, suggesting that the S III sites demonstrate the highest stability for accommodating H atoms within the NGI.The three types of stable configurations of solute H atoms are further extracted from Fig. 6 (b) and exhibited in detail in Fig. 6(c).The solute H atoms occupying the S I sites exhibit a preference for locating themselves within the central position of a hexagonal pyramid formed by one Ni atom and six C atoms.Conversely, other variants of H atoms position themselves at the center of tetrahedra composed of three Ni atoms and one C atom.As the solute H switches from site S I to S II to S III , the Ni-H distance gradually increases from 1.676 to 1.986 Å, whereas the C-H distance starts almost unchanged (i.e., 1.798 Å) and then decreases to 1.103 Å, implying that the formation volume of solute H tends to expand and then contract.Hallil et al. [60] have demonstrated that the segregation energy of the solute H near a fcc-Ni GB relates to its atomic volume (Ω H ), exhibiting a linear increase with the increasing Ω H while remaining constant (i.e., 0.25 eV) after a critical volume (Ω H ≈ 6.30 Å 3 ).Based on the Voronoi tessellation method [31,46], the atomic volume of the solute H at the sites S I , S II , and S III is calculated to be 7.23, 7.45, and 7.14 Å 3 , respectively.Clearly, although the Ω H at different sites within the NGI exceeds the critical volume within the Ni GBs, their segregation energies are not constrained but have further increased.Moreover, an extended linear relationship exists between the segregation energy and the Voronoi volume of the solute H at the sites S I and S II , while at the site S III there is an anomaly that does not follow any of the previously found laws [60].Generally, there is a typical range for the sp 3 C-H bond length from 1.056 to 1.115 Å, and the corresponding bond energy is approximately 413 kJ mol − 1 or 4.302 eV [65].The distance of 1.103 Å indicates that the solute H forms a sp 3 C-H bond at the S III sites, which can result in a large release of heat in the system as well as the local configuration change of graphene.The high covalent bond energy ensures that the solute H at the S III sites has the most stable configuration, making it difficult to break away from the graphene and thus fall into a non-diffusible state.This also confirms the above MD observation and highlights the role of NGIs as benign traps.On the other hand, the bond energy released at the S III sites is nearly twice the hydrogen segregation energy, and this extra energy may be mainly dissipated in the lattice displacement and twisting of C atoms, which leads to the fact that the hydrogen properties within the NGI are clearly different from those of the Ni GBs.

Kinetics of hydrogen diffusion and trapping
Because solute H is more stable at the OIS of fcc Ni as mentioned above, it preferably migrates between two adjacent OIS.However, direct jumps along the path connecting two adjacent OIS are extremely difficult, since the energy landscape does not have a saddle point present along this path but has a local maximum [36].To achieve a smooth jump, the intermediate path image is usually positioned at the TIS in the middle between the adjacent OIS.This is exactly the mirror image of the jumps of the He atoms along the TIS-OIS-TIS path [31].Hence, only migrations along an OIS-TIS-OIS path are considered for the solute H in fcc Ni.Following this most probable path, the kinetic details of solute H jumping from the bulk to the NGI can be unveiled by using the climbing-image nudged elastic band (CI-NEB) method [66], including the transition state, optimum migration path, and diffusion energy barrier, etc. Six different migration paths related to the initial spatial coordinates of solute H are investigated here, and the corresponding diffusion barriers toward the interface are shown in Fig. 7(a).The solute H diffusion barrier is 0.470 ± 0.006 eV in the bulk, almost equal to the result (i.e., 0.46 eV) in pure Ni provided by Angelo et al. [44] and close to the diffusion activation energy obtained from the above MD simulations.As solute H approaches the NGI, the diffusion barrier gradually decreases until 0.122 eV, confirming that within a specific range of the NGI, the solute H exhibits a distinct preference for migrating towards the interface.Moreover, to observe the atomic transition of solute H in detail, one representative diffusion path is displayed in Fig. 7(b).
H. Huang et al.Significantly, the solute H at the OIS does need to temporarily reside in its adjacent TIS after crossing a high barrier to jump to the next OIS.Once the solute H enters the interface, the energy of the system will decrease, with a maximum decrease of 0.537 eV.After overcoming a small barrier of 0.122 eV (i.e., from site K to L in the inset of Fig. 7(b)), the solute H tends to fall into the S II sites within the core of the NGI and then jump between two adjacent S II sites, with a diffusion barrier of 0.08 eV (e.g., from site L to M in the inset of Fig. 7(b)), which has also been observed in other migration paths.In addition, following comprehensive testing, it has been determined that for the transition of solute H from site S II to S III (see Fig. S3), a diffusion barrier of approximately 0.139 eV must be surpassed, which is slightly higher than that between two adjacent S II sites.These suggest that solute H exhibits enhanced migratory ease within the NGI compared to its mobility in the bulk, with a notable theoretical preference for transitioning between two adjacent S II sites.However, due to the attraction of the vibrational C atoms activated by ambient temperature, the solute H at the S II sites is too late to migrate along the most favorable path within the NGI, but instead crosses a slightly higher barrier and is then trapped in the S III sites to form sp 3 C-H bonds.In other words, the direct jumps along the S II -S II path are not always feasible, especially at high temperatures, as observed in the above MD simulations.But in turn, a solute H located at the S III site must surpass a saddle point that has a barrier height of at least 2.0 eV to move back to the neighboring S II site.Moreover, H atoms require overcoming a much higher barrier while attempting to migrate from site S II (or S III ) to S I , which is difficult to achieve under the conditions of nuclear reactors.On the other hand, H atoms can hardly penetrate graphene through the six-membered carbon ring due to the difficulty of jumping to the S I sites, which confirms the findings of previous studies indirectly [29].Therefore, the NGI is only possible to offer two types of migration paths for a solute H, i.e., S II -S II or S II -S III -S II .However, in practice, the likelihood of both being activated is almost negligible.Overall, the kinetic observations reveal a fascinating phenomenon: solute H atoms near the NGI exhibit a preference for migrating towards it due to an exceptionally low diffusion energy barrier, serving as an alternative driving force distinct from the above-mentioned energetics.Meanwhile, H atoms captured by the interface are almost non-diffusible, thus preventing them from forming hydrogen molecules or bubbles at the interface.
It has also been observed that the entry of H atoms from the bulk into the NGI causes pronounced lattice distortions.For example, the interatomic distances around the solute H atom change significantly as it moves from site C to K in the inset of Fig. 7(b).This suggests that the formation volume of solute H may also be one of the physical reasons that affect its kinetic behavior near the NGI.Moreover, the formation volume of solute H can be regulated by the thermal expansion or contraction of the lattice due to ambient temperature [58].Consequently, when using the CI-NEB method to predict the kinetic behavior of solute H in the vicinity of the NGI, it becomes essential to consider the temperature-dependent nature of the diffusion energy barrier, which contributes to a comprehensive comparison with the preceding MD results.This issue needs to be tackled in the future with a more advanced model that accurately depicts the features of the NGI− H system under high temperatures.Additionally, hydrogen diffusion in traditional polycrystalline metals and alloys is affected by various parameters, such as crystal grain size, GB energy, and excess free volume [53].For instance, alterations in the crystal grain size can have a dual effect on hydrogen diffusion [67].On one hand, GBs boost hydrogen mobility by increasing surface area per unit volume as the grain size decreases, thus facilitating hydrogen diffusion along them.On the contrary, a microstructure with smaller grains has a significantly increased density of triple junctions.These junctions can serve as trap sites for H atoms and lead to a decrease in hydrogen mobility.These competing effects indicate the presence of a particular grain size that maximizes hydrogen diffusivity.Typically, the elevated hydrogen diffusivity along GBs could enhance the transport of solute H, ultimately leading to its segregation into critical zones [53].This phenomenon subsequently induces the creation of new vacancies and a subsequent decrease in the cohesive strength of GBs, thereby initiating HE failures [53].Consequently, there are still substantial challenges in utilizing the "GB engineering" strategy to mitigate the intergranular embrittlement of metallic materials in the presence of solute H.However, following the definition of GBs [25] and prior studies [31,32,68,69], it can be deduced that NGIs are a class of interfaces featuring low energy levels which act as reservoirs for trapping of solute H. Based on the preceding analysis, solute H should be capable to be uniformly distributed alongside these NGIs, successfully hindering its mobility and segregation towards vulnerable locations within the composites.Therefore, the strategy of "heterointerface engineering" from an NGI perspective may be more expedient and effective for designing HE-resistant metallic nickel.

Conclusions
In summary, this study has investigated the effect of graphene and/ or NGIs on the behavior of hydrogen diffusion and trapping in their vicinity, using atomistic simulations.The thermodynamic results showed that as the temperature increases, the individual diffusion distance of solute H exhibits an increase, while its total diffusion distance and frequency experience a decrease.When the temperature exceeds 900 K, the solute H ceases to exhibit intermittent jumps after a certain duration and is instead trapped by the NGI.This circumstance leads to a gradual transition of the diffusion mechanism from two-dimensional to one-dimensional over time, manifested specifically as diffusion along the [111] direction.There is a noticeable difference in the MSD for solute H diffusion between NGNC and pure Ni.This difference can be substantial, spanning two orders of magnitude, highlighting the hindered rate of solute H diffusion in the NGNC.The presence of NGI in the composites leads to reduced hydrogen diffusivity in their vicinity compared to single-crystal Ni.This effect is particularly pronounced at elevated temperatures, leading to slower hydrogen diffusion, and consequently facilitating its entrapment in the interfacial regions.Energetic and kinetic calculations have been employed to comprehend the dynamic behaviors of solute H further.The segregation energy of solute H near the NGI demonstrates a propensity to rise and subsequently decline, peaking at 2.484 eV.Notably, the minimal segregation energy near the interface is comparable to the maximum value near the Ni GBs.Discerning the disparity in formation energy clarifies the existence of three distinct and stable sites for H atoms within the core of NGI.At the S III sites, solute H readily forms a sp 3 C-H bond with the C atoms on the surface of graphene as the distance between C-H falls within the typical bond length range.This bond formation produces significant heat release in the system and causes local configuration changes in graphene.The robust covalent bond energy guarantees the utmost stability of the soluble H at the S III sites, thereby impeding its detachment from graphene and subsequent entry into a non-diffusible state.Upon approaching the NGI, the diffusion barrier gradually diminishes to 0.122 eV, affirming the solute H's distinct inclination to migrate toward the interface within a specified NGI range.The NGI has the potential to provide two distinct migration paths for solute H. Nonetheless, the practical activation likelihood for both paths is extremely minimal.The present study advances understanding of the behaviors of hydrogen permeation in NGNC.Achieving an appropriate dispersion of graphene within the matrix can considerably reduce the diffusion coefficient of solute H in the material.This, therefore, guides the active H atoms towards non-diffusible H atoms.Such a process can successfully provide an efficient strategy to impede hydrogen penetration in the material and establish a persistent hydrogen barrier.This approach may enable the successful avoidance of HE failures in high-strength Ni-based alloys, ultimately ensuring the reliability of nuclear reactor operation.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Atomic model of the NGNC with a solute H atom at 2.61 nm away from the graphene.

Fig. 3 .
Fig. 3. MSDs of a single H atom in pure Ni (a) and the NGNC system (b) during 5.0 ns in the temperature range of 500-1200 K.

Fig. 4 .
Fig. 4. MSD components along different directions obtained from a single H atom in pure Ni (a) and the NGNC system (b) in the temperature range of 500-1200 K.

Fig. 5 .
Fig. 5. Diffusion coefficients for a single H atom in pure Ni (a) and the NGNC system (b) at various temperatures, along with their components in different directions.The parameters obtained from linear fitting using the Arrhenius formula are presented.

Fig. 6 .
Fig. 6.Distributions of the formation energies of solute H. (a) Formation energies of solute H varying with its distance from the NGI.The segregation energy and effective interaction range of NGI towards the solute H are showcased.(b) Formation energies of solute H for different stable sites within the core of NGI.(c) Three types of stable atomic configurations of solute H.

Fig. 7 .
Fig. 7. Diffusion of solute H toward the NGI.(a) Diffusion energy barriers of solute H varying with its distance from the NGI along various paths.Stable H atoms located in the OIS and TIS are presented.(b) Transition energy profile for hydrogen diffusion towards the interface.The inset exhibits the atomic configurations of hydrogen at various sites along the diffusion path.