Hydrogen trapping and diffusion in polycrystalline nickel: The spectrum of grain boundary segregation

Hydrogen as an interstitial solute at grain boundaries (GBs) can have a catastrophic impact on the mechanical properties of many metals. Despite the global research effort, the underlying hydrogen-GB interactions in polycrystals remain inadequately understood. In this study, using Voronoi tessellations and atomistic simulations, we elucidate the hydrogen segregation energy spectrum at the GBs of polycrystalline nickel by exploring all the topologically favorable segregation sites. Three distinct peaks in the energy spectrum are identiﬁed, corresponding to different structural ﬁngerprints. The ﬁrst peak ( − 0.205 eV) represents the most favorable segregation sites at GB core, while the second and third peaks account for the sites at GB surface. By incorporating a thermodynamic model, the spectrum enables the determination of the equilibrium hydrogen concentrations at GBs, unveiling a remarkable two to three orders of magnitude increase compared to the bulk hydrogen concentration reported in experimental studies. The identiﬁed structures from the GB spectrum exhibit vastly different behaviors in hydrogen segregation and diffusion, with the low-barrier channels inside GB core contributing to short-circuit diffusion, while the high energy gaps between GB and neighboring lattice serving as on-plane diffusion barriers. Mean square displacement analysis further conﬁrms the ﬁndings, and shows that the calculated GB diffusion coeﬃcient is three orders of magnitude greater than that of lattice. The present study has a signiﬁcant implication for practical applications since it offers a tool to bridge the gap between atomic-scale interactions and macroscopic behaviors in engineering materials. © 2023 Published by Elsevier Ltd on behalf of The editorial oﬃce of Journal of Materials Science & Technology. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Hydrogen, the smallest and most abundant element in the universe, is recognized as a highly efficient and clean energy carrier that can help achieve a zero-emissions future.As hydrogen-based energy systems become more prevalent, material challenges associated with hydrogen storage and transport remain a bottleneck needing to be addressed.One significant issue is hydrogen embrittlement (HE), where dissolved hydrogen in metals leads to dramatic degradation of mechanical properties and potentially catastrophic failure.Since first observed in 1875 [1] , numerous research effects have been contributed to elucidating the nature of HE and many mechanisms have been proposed [2][3][4][5][6] .To name some of them, there are hydrogen-enhanced decohesion (HEDE) [ 7 , 8 ], hydrogen-enhanced localized plasticity (HELP) [ 9 , 10 ], hydrogenenhanced strain-induced vacancies (HESIV) [ 11 , 12 ].All the proposed mechanisms are supported by experimental and theoretical evidence, however, due to the multifaceted nature of HE, it is challenging to apply any one paradigm universally.Researchers today tend to suggest that HE is a result of the synergistic interactions of multiple mechanisms [ 4 , 13 ], without enumerating the quantitative contribution from each mechanism or the critical condition for one dominating mechanism.One of the reasons for the current inability of quantifying the contribution of HE mechanisms is the lack of proper understanding of hydrogen-microstructure interactions.
Of all the HE-related phenomena, hydrogen-induced intergranular fracture represents a particularly severe degradation mode which is accompanied by flat and smooth fracture surface and decreased toughness [14][15][16][17][18][19] .The intergranular fracture is usually attributed to that the saturation of interstitial hydrogen atoms at GB reduces its cohesive energy and promotes cleavage-like failure, as proposed by the HEDE mechanism.This concept helps to interpret experimental observations and develop predictive models for quantitative assessment of hydrogen-related degradation [ 20 , 21 ].In contrast, HELP theory [ 10 , 22 ] tends to conclude that it is the hydrogen-influenced plasticity that leads to the attainment of critical hydrogen concentration at GB and results in final decohesion.Furthermore, recent studies have revealed the significant or even dominant role of GB vacancies in facilitating and controlling intergranular fracture [23][24][25] .All the evidence highlights the fundamental role of hydrogen-GB interactions in governing intergranular fracture.
To gain a better understanding of hydrogen-GB interactions and thus intergranular fracture, it is a prerequisite to have a precise description of hydrogen segregation and transport in polycrystalline materials.Multiple obstacles currently impede progress in this endeavor.Firstly, although a large experimental data reached agreement on hydrogen diffusivity in bulk [2] , the reported value of hydrogen diffusivity at GB, i.e., GB diffusion coefficient D GB shows a large dispersion spanning over several magnitudes [25][26][27][28][29][30][31][32][33] .The D GB is calculated by comparing diffusivities in single crystal and polycrystal, however, complex experimental setups, contributions from other crystal defects, and imprecise estimation of GB volume could all lead to incoherent results.Secondly, hydrogen diffusion and segregation behavior at typical GB based on the coincidence site lattice (CSL) model have been elucidated by a large number of atomistic simulations [34][35][36][37][38][39][40][41] .Nevertheless, those analyses are limited to one or several specific systems which are not necessarily representative of a real polycrystal network dominated by random GBs [25] .A collective description of hydrogen segregation and diffusion map in polycrystals is needed to bridge the gap between experiments and theories.Thirdly, most current hydrogen segregation models to describe equilibrium hydrogen concentration at GB originate from the one derived by Oriani [42] , which assumes single-value segregation energy by fitting experiments.In the 1980s, Kirchheim [6] considered the diversity of GB sites by assuming that hydrogen segregation energies at GBs fall into a single Gaussian distribution.However, the specific nature of this diversity was speculative at that time, and there is still a lack of direct demonstration at the atomic level.
Substitutional segregation, which also addresses the diversity of sites at GBs, has been a highly active research area in recent years.White et al. [ 43 , 44 ] introduced the concept of binding energy spectra to explain experimental segregation data in polycrystalline metals, going beyond a single segregation energy.More recently, Wagih and Schuh [45][46][47] demonstrated the anisotropic nature of GB segregation by extracting the spectra of GB segregation energies from atomistic simulations on various alloying systems.They derived a thermodynamic framework that incorporates these spectra to predict the equilibrium GB segregation state, successfully transitioning from segregation energy to GB enrichment.However, it is important to note that these studies focused exclusively on substitutional segregation and cannot be applied to cases of interstitial segregation.In hydrogen-metal systems, such as hydrogen-nickel systems, hydrogen acts as an interstitial solute at the GB network of polycrystals, enhancing the susceptibility to intergranular fracture.In such cases, a solely thermodynamic model is insufficient to account for hydrogen diffusion and trapping near GBs, kinetic analysis is needed to fully understand these phenomena [ 4 8 , 4 9 ].
In this study, we design a new computational route combining Voronoi tessellations and molecular statics to identify the interstitial sites and calculate the hydrogen segregation energy at all geometrically favorable interstitial sites in a polycrystal consisting of general GBs.Polycrystalline nickel is chosen as a representative host system because of its broad industry application and susceptibility to HE [50] .The hydrogen segregation spectrum in nickel is extracted as a collection of hydrogen segregation energy at millions of interstitial GB sites.GB core and GB surface sites are identified based on their distinctive energy profiles and unique geometrical arrangements.The spectrum was incorporated with a thermodynamic framework to predict the equilibrium hydrogen distribution in polycrystals.Hydrogen diffusion and trapping behavior is studied both statically by comparing neighbor energy gap and dynamically through mean square displacement analysis, with the GB diffusion coefficient D GB precisely calculated over a range of temperatures.

Computational route
We utilize large-scale atomistic simulations combining Voronoi diagram [ 36 , 51 ] to identify and calculate hydrogen segregation energy at all geometrically favorable interstitial sites in polycrystalline nickel.The atomic interactions are described by the embedded atom method (EAM) potential developed by Angelo et al. [ 52 , 53 ], which well captures hydrogen interaction with a wide range of defects.Molecular dynamic and molecular statics simulations are performed using the LAMMPS code [54] , Voronoi tessellations are constructed with the Voro ++ package [55] , and atomic structure is visualized in Ovito [56] .
Several nickel polycrystal models are built by filling cubes of varying size (14-30 nm) with randomly oriented grains, constructed by Voronoi tessellation using the Atomsk toolkit [57] .The default model in Fig. 1 (a) (20 nm) 3 has 733,655 atoms and 12 randomly oriented grains, with 122,688 atoms composing 68 individual GBs.[58] .The polycrystal is structurally optimized using the conjugate gradient algorithm, followed by thermal relaxation in the isothermal-isobaric (NPT) ensemble at 300 K for 300 ps.The equilibrium structure is then cooled down to 0 K with a rate of 1 K/ps and finally minimized with an energy tolerance of 10 −25 .
For each nickel atom in the polycrystal, we draw Voronoi cells consisting of all points of the plane closer to that central nickel atom than to any other neighbor nickel atoms.The vertices of all the cells are collected as potential hydrogen segregation sites, with a total of 4,019,870 sites identified in the default model (20 nm 3 ).These sites are divided into GB trapping sites and lattice trapping sites, by comparing their distance to the nearest GB atoms with a cutoff of 2.2 Å.Molecular statics are utilized to calculate the hydrogen solution energy at any site i , E sol i , defined as: where E i +H and E i are the total energies of the system before and after trapping of one hydrogen atom at site i, respectively, and E H 2 is the energy of one isolated H 2 molecule in vacuum.The hydrogen segregation energy at GB, , is defined as: where and E sol oct are the total energies of the system with a hydrogen atom trapped at the GB site i and an octahedral site in a perfect fcc lattice, respectively.Negative values of GB segregation energy indicate the GB sites are energetically more favorable for hydrogen absorption and trapping than the lattice.We screen the pair distance and energy gap between all trapping sites, pairs with the closest distance (cutoff 0.1 Å) and limited energy gap (cutoff 0.1 eV) are counted only once to avoid overlapping sites.The segregation energy of all GB sites after screening is plotted into a histogram to draw the final GB segregation spectrum.It should be noted that the computational framework presented in this study can be extended to other interstitial segregation systems.With the appropriate interatomic potential and parameters, the segregation spectrum can be obtained for additional systems using the procedures involving identification through geometric tessellation and calculation through molecular statics.

Thermodynamic model for determining equilibrium hydrogen distribution
Given sufficient time to diffuse, all the hydrogen inside lattice and GB will redistribute and reach equilibrium.Following Mclean's isotherm [59] , Fermi-Dirac statistics incorporated with the single segregation energy value or segregation spectrum have been successfully applied in previous studies [ 45 , 47 , 60 , 61 ] to predict substitutional solute segregation at GBs. Regarding hydrogen as interstitial segregation at GB, Oriani [42] introduces a theoretical model to describe the relationship between equilibrium hydrogen concentration at GB and lattice: where c l and c GB are the atomic hydrogen occupation in the lattice and GB, respectively, E seg GB is the hydrogen-GB segregation energy, k is Boltzmann's constant and T is temperature.One thing to be noted is that the hydrogen occupation denotes the fraction of atoms (ratio between hydrogen atoms and nickel atoms) at a specific site, not the fraction of atoms in a unit volume.However, this model could be oversimplified to take E seg GB as one certain value by assuming there is just one type of GB segregation site for hydrogen.This single-value assumption is robust for comparing segregation feasibility between different solute systems at the same lattice concentration, however, it neglects the natural variation of segregation at different GBs and different sites which could lead to biased predictions when applied in a broader range.This issue is addressed by White et al. [44] and Kirchheim [6] , elaborated by Wagih and Schuh [45][46][47] , through introducing a set of segregation sites with different energy levels: where c i ( GB ) is the atomic hydrogen occupation at a specific GB site i and E seg i ( GB ) is corresponding segregation energy.By taking c i ( GB ) as the probability of hydrogen at site i , the total hydrogen occupation at GB could be viewed as c i ( GB ) contributed from all GB sites: where c GB is the total hydrogen occupation at GB and N is the total GB sites, the rightmost term in the equation is obtained by substituting Eq. ( 4) into the middle term.Thus, Eq. ( 5) describes the complete GB segregation by taking into account the energy difference and density distribution of all GB sites.Kirchheim et al. [6] have assumed a Gaussian distribution of hydrogen segregation energies at GBs to fit the experimental data, however, the full spectrum describing atomic hydrogen-GB segregation is still lacking which makes the application of Eq. ( 5) difficult.

Theoretical analysis for hydrogen diffusivity
Compared to the equilibrium hydrogen distribution, which can be directly evaluated by the energy differences at varying sites, hydrogen diffusion is a more complex and path-dependent process where analysis of dynamic hydrogen-GB interaction is needed.Experimentally, the hydrogen diffusion coefficient in single crystal or polycrystal could be measured by the permeation approaches [ 25 , 62 ].The diffusivity of hydrogen at GB could be derived linearly by subtracting the contribution by the lattice part from the whole polycrystal as follows: where D lattice , D GB , and D poly are the diffusion coefficients in the lattice, GBs, and the whole polycrystal, respectively, V lattice , V GB , and V poly are the volume of lattice, GBs, and the whole polycrystal, respectively, with V lattice + V GB = V poly .Alternatively, by regarding the polycrystal as a composite of grains embedded in an intergranular matrix composed of random GBs where the empirical Hashin-Shtrikman upper bound (HS + ) model can be applied [ 25 , 63 ], the GB diffusion coefficient could be derived from the following formula: where f GB = V GB /V poly is the volume fraction of GB.Actually, given the sufficient difference between D lattice and D GB , the value from Eq. ( 7) is equal to Eq. ( 6) in most cases.However, a significant amount of uncertainty in the material system has resulted in a substantial discrepancy in the calculated D GB as mentioned.Many factors such as pre-existing dislocations and vacancies, and inaccurate GB volume estimation could lead to inaccurate predictions.For example, a corrected estimation (5 Å in this paper) of GB thickness compared to Oudriss' work [25] will lead to a change of one or two magnitudes in the calculated D GB value.Theoretically, one typical method [ 37 , 64 ] to calculate D GB is the combination of nudged elastic band method (NEB) [65] and density functional theory (DFT), where the NEB method provides a way to find a minimum energy path given the initial and final states and DFT is used to calculate the energy of associated configuration.The maximum energy jump in the path is taken as the energy barrier, which can be used to estimate the ratio between GB and lattice diffusivities by assuming the same jump frequency for hydrogen in lattice and GB [37] : (8) where E bar GB and E bar lattice are the energy barrier (maximum energy minus minimum energy in the diffusion path) of the GB region and lattice region, respectively, k is Boltzmann's constant and T is temperature.However, for complex systems where the migration path is not known as a priori, the NEB is not applicable and mean squared displacement (MSD) analysis is needed [ 66 , 67 ].Molecular dynamics can provide the trajectory − → r i (t) of the particle useful for the MSD analysis, where the hydrogen diffusion coefficients D in a three-dimensional system can be determined as follows from Einstein [68] : where − → r i (0) and − → r i (t) are the position of the i th hydrogen atom at time 0 and time t , respectively, N is the total number of hydrogen atoms in the system, and MSD is equal to 1 By writing the D poly into the sum of MSD from lattice hydrogen and GB hydrogen we obtain: where N lattice , N GB , and N are the numbers of hydrogen atoms in the lattice, GBs, and the whole polycrystal.If the hydrogen atoms are randomly distributed in the whole polycrystal during the MSD collection period, the hydrogen proportion can be viewed as equal to their volume proportion: . And then Eq. ( 10) is transformed into Eq.( 6) .To be noted, most systems in the previous calculations are limited to CSL bicrystal models which could not address the collective diffusion behavior among all random GBs and along all directions.

Properties of the spectrum
By exploring all the interstitial segregation sites through molecular statics, the map of all potential hydrogen trapping sites with their segregation energy is shown in Fig. 1 (b).As indicated by the color, octahedral (green) and tetrahedral (red) interstitial sites in the lattice display a layer-by-layer pattern (visualization effect), while the GB network provides a number of energy-favorable sites (cyan and blue) in the region associated to those interstitial sites.The distance between interstitial sites and the nearest GB atom is defined as the interstitial-GB distance, and its relationship with segregation energy is plotted in Fig. 1 (c).Only the closest cluster (blue) within a cutoff distance of 2.2 Å is taken as GB segregation sites, and the remaining interstitial sites (pink) are regarded as lattice trapping sites.It is clear both qualitatively in Fig. 1 (b) and quantitatively in Fig. 1 (c) that GB sites show a broader region of energy distribution compared to the lattice, which is attributed to its diverse nature and complex local structures.In contrast, the energy distribution of lattice sites is monotonous and only oscillates around 0 eV and 0.4 eV in line with typical octahedral and tetrahedral sites, respectively.
After a screening procedure to avoid overlapping sites, a total number of 491,799 GB interstitial segregation energies in the default model is extracted to draw their probability density histogram in Fig. 2 (a), i.e., the energy spectrum of hydrogen-GB segregation in nickel.The calculated spectrum shows a well-fitted Gaussian mixture distribution, three peaks with distinct energy levels are identified after deconvolution.The separated Peak1, Peak2, Peak3, and their sum are denoted by red, green, blue, and black dash curves, respectively.The fitted probability density distribution of segregation energy f ( E seg GB ) can be expressed as: where u i , σ i , and w i are the mean, standard deviation, and weight of each peak, respectively, and the integral of probability den- ( u 3 = 0 .287 eV , w 3 = 0 .17 ) have a more similar energy level to lattice octahedral (green vertical line) and tetrahedral (blue vertical line) segregation energies.The segregation energies at octahedral and tetrahedral sites are directly obtained from the two crests in Fig. 1 (c).Thus, the trapping behavior in the GB region is dominated by the Peak1 sites due to their more energy-favorable nature.The cumulative probability density is shown in Fig. 2 (b).Approximately 80% of sites are favorable for hydrogen GB segregation with a neg-ative value, with the lowest ( u 1 − 3 σ 1 ) segregation energy located at −0.376 eV.In comparison, a bunch of segregation energies based on typical coincidence site lattice (CSL) [69] GB models are listed in the literature [ 37 , 40 , 52 , 35 , 70 , 71 ].Most of the reported energies from CSL models belong to the top 20% of low-energy sites, which account for exactly the left half of Peak1.Only some special GBs are located in the relatively high-energy area.On the one hand, it reflects the CSL model can capture the key feature of those most favorable trapping sites at GBs and highlight the importance of Peak1 sites during hydrogen occupation.On the other hand, it indicates that any single or combination of the CSL GBs cannot be representative of hydrogen interstitial segregation in a full polycrystal.Thus, to verify the generality of this spectrum, we repeat the computational route for several polycrystals with varying sample sizes and grain sizes.The resulted spectrums are plotted in Fig. 2 (a, b) with edge lines of different colors, which agrees well with the original spectrum extracted from the default model in both probability density and cumulative probability density.These spectrums are further fitted into the Gaussian mixture model, with the parameters listed in Table 1 .The location ( u i ) and shape ( σ i ) of each peak from different polycrystal samples match perfectly, proving the universality of the hydrogen-GB spectrum.
The small deviation in weights ( w i ) of each peak is acceptable given that small errors could be introduced during the classification between lattice and GB interstitial sites.Thus, this spectrum is a unique representative of the statistical hydrogen-GB interaction in polycrystals which provides a robust tool to predict equilibrium hydrogen concentration, trapping, and diffusion behavior in a polycrystal.We will elucidate the usage of the spectrum in the subsequent Sections 3.3 and 3.4 .We have noticed that Waigh and Schuh [72] reported a similar two-Gaussian mixture distribution in the hydrogen-palladium system.Its similarity with the hydrogen-nickel system discussed in this manuscript suggests that the Gaussian mixture distribution may provide a universal representation of interstitial hydrogen segregation at GBs in fcc metals.

Interstitial sites at GB core and surface
In order to comprehend the three-peak spectrum, it is important to explore the fingerprint microstructure of each peak and their underlying structure-property correlation.To track this, we first divide the spectrum into two parts using a vertical line with the mean value of Peak1, u 1 in Fig. 3 (a).Only the left part ( E seg GB < −0 .205 eV ) of the spectrum is visualized in Fig. 3 (a1) according to the energy level, while the whole spectrum is visualized in Fig. 3 (a2), but the sites belonging to the right part ( E seg GB > −0 .205 eV ) are all colored red.The primary observation obtained from the visualization is that the low-energy sites are completely surrounded by sites with relatively high energy levels, suggesting that the geometric arrangement of GB sites is organized based on their respective energy values.It agrees with previous CSL model studies where core structures of GB behave as the deepest trapping sites and are mostly located at the left half of peak1 ( E seg GB < −0 .205 eV ), as shown in Fig. 2 (b).The GB sites are further divided into core and surface sites according to their local atomic environment as described below.The blue peak in Fig. 3 (b) shows the density distribution of GB core sites and the remaining olive  part corresponds to GB surface sites.Only GB core sites are shown in Fig. 3 (b1) colored by their energy level, while GB surface sites are all colored olive in Fig. 3 (b2).The significantly overlapping area (90% of the Peak1 area) between GB core sites and Peak1 (red dash curve) in the spectrum obviously shows their statistical similarities, which is also demonstrated by the similar hierarchical structure between GB core and GB surface in Fig. 3

(b2).
To distinguish between the GB core and surface sites structurally, merely using a distance parameter between interstitial sites and the nearest GB atoms is not sufficient, a more sophisticated local environment descriptor is needed.The number of neighboring fcc atoms n and their proportion fn among all the neighbors are chosen to describe the local structure.n denotes the number of atoms in the vicinity of an interstitial site belonging to the fcc (lattice) structure, where the non-fcc atoms are treated as GB atoms using common neighbor analysis (CNA) [73] .fn corresponds to the ratio between n and the coordination number of the interstitial site.Typically, the coordination number denotes the number of nearest neighbor atoms.For a perfect fcc lattice, the coordination number of an octahedral site is 6 setting the cutoff as 1.76 Å ( 1 / 2 a 0 , where a 0 is the lattice constant).To take into account the effect of the second nearest neighbor, the coordination number could be increased to 14 by setting the cutoff as 3.05 Å ( √ 3 / 2 a 0 ).For an interstitial site in GB, the coordination number can vary dramatically (4 to 15) due to its complex nature.The n and fn are taken as the measure of similarity of the structures between tested GB sites and perfect fcc lattice sites, with low values indicating the GB core sites with less similarity and high values indicating the GB surface with more similarity.Setting a proper cutoff for n and fn could then help to classify the GB core ( n < n cutoff or fn < fn cutoff ) and surface ( n > = n cutoff or fn > = fn cutoff ).The performance of the classification is measured by the nonoverlapping area between Peak1 and the identified GB core in the spectrum, with a lower value indicating better performance.The results from combinations of different neighbor distance cutoffs (1.6-3.2Å) and n cutoff (0 −15) or fn cutoff (0 −1) are shown in Fig. 4 (a, b), with the corresponding GB core proportion indicated by the blue area in the spectrum in Fig. 4 (a1 −a3) and (b1-b3).The GB core area starts to fill Peak1 with increased n cutoff and fn cutoff , while discrepancy also generates around the right corner of Peak1.The balance between increased overlapped area and generated discrepancy is reached at the minimum nonoverlapping area of 0.15 given n cutoff = 3 or fn cutoff = 0.3 with a 3.05 Å neighbor cutoff.At this point, the selected GB core sites show the largest coincidence with the Peak1 area in Fig. 4 (a2, b2).Thus, by selecting appropriate parameters, we can quantitatively connect the low-energy sites in Peak1 to a group of sites with specific local geometrical arrangements, i.e., the GB core sites.

Equilibirum hydrogen concentration at GB
By combining the atomic spectrum presented in this study with Eq. ( 5) , we are able to construct a comprehensive hydrogen segregation map for polycrystals.The full spectrum prediction is derived directly from the discrete spectrum, which is based on all 491,799 segregation sites without any approximation.By assigning weights to each site based on their probability density, we can convert the prediction into a continuous description using the following equation: where c GB is the total hydrogen occupation at GB, c i ( GB ) ( E seg GB ) is atomic hydrogen occupation at a specific GB site i with segregation energy E seg GB , f ( E seg GB ) is probability density distribution of E seg GB .By substituting Eq. ( 11) into the middle term, we obtain the Gaussian mixture model prediction in the rightmost of Eq. ( 12) .The predicted hydrogen occupations at GB c GB as a function of lattice hydrogen occupations c l at different temperatures (300 K, 450 K, 543 K) are shown in Fig. 5 (a) with the detailed segregation states illustrated in Fig. 5 (c1 −c9).The perfect agreement between discrete full spectrum prediction (solid) and continuously fitted spec-trum prediction (dash) at three temperatures further proves that the spectrum has been well captured by the Gaussian mixture model.Taking an example of a lattice occupation c l = 0 .001 and temperature T = 300 K, the predicted GB occupation is 0.303 which means 30.3% of GB interstitial sites have been filled with hydrogen atoms at equilibrium.Among those filled interstitial sites, the contribution from the 1st Gaussian (dotted), i.e., Peak1, is 0.273, indicating that 90% of GB trapped hydrogen is concentrated at the GB core region.It is worth noting that in cases where the hydrogen concentration exceeds the dilute limit, the repulsive interac-  13) , the prediction based on Gaussian mixture model in Eq. ( 14) , and the contribution of Peak1 are illustrated by the solid, dash, and dot lines, respectively.(b) Effective single-value segregation energy E seg GB by substituting c GB and c l from (a) into Eq.( 11) .(c1 −c9) Detailed equilibrium hydrogen segregation state at GB with varying lattice hydrogen occupations and temperatures.The GB hydrogen occupations c GB and contribution from Peak1 is denoted by solid and dot curves, respectively, with green, yellow, and orange colors.The spectrum is filled with blue and the Peak1 is indicated by the red dash curve.tion between hydrogen atoms becomes significant and can lead to a lower real hydrogen occupation compared to the theoretically predicted value here.As c l keeps increasing to 0.005, the c GB increases to 0.423 and the contribution from 1st Gaussian increases to 0.35.While the newly inserted hydrogen keeps pumping into the GB core region, the proportion of GB core hydrogen decreased to 82.6%.This is because the GB core sites are closed to their full segregation value (0.424, w 1 , the weight of Peak1) and GB surface starts to take the responsibility to accommodate hydrogen atoms.However, in most experimental charging levels [ 14 , 74 , 75 ], the hy-drogen occupation is below 0.005 (5000appm) and it is still the GB core sites that dominate the nickel GB trapping behavior in this region.As the temperature increases, the equilibrium hydrogen occupation at GB decreases but the GB core still plays the main role in trapping hydrogen in the displayed region.
To further compare with experiment results based on the single-value GB segregation assumption as in Eq. ( 3) , the concept of effective segregation energy is employed.As previously derived by Steigerwald et al. [76] and Wagih et al. [45] , it can be written with a discrete form in the middle or a continuous form in the rightmost: where the c GB and c i ( GB ) are the total GB hydrogen occupation and atomic hydrogen occupation at a specific GB site i , respectively, E seg GB is the effective segregation energy by deriving the spectrum back into one equivalent value in Eq. ( 3) at a specific temperature and lattice concentration.The result E seg GB as a function of c l is shown in Fig. 5 (b).As the c l increases, the absolute value of E seg GB decreases because the GB core hydrogen proportion decreases and newly occupied GB surface trapping sites weaken the overall GB trapping capacity.While at a higher temperature, less hydrogen is trapped at GB and E seg GB is represented by fewer amount of sites but with lower average segregation energy (more energy favorable), which makes the absolute value of E seg GB higher.We also attach an example comparing recent thermal desorption spectroscopy (TDS) results from Wada et al. [74] , the experimentally calculated singlevalue segregation energy of about 0.2 eV at 543 K and c l = 1 ‰ agrees well with our prediction.

Hydrogen diffusion and trapping at GB
To study hydrogen diffusivity in polycrystals, we would like to revisit the coexistence of short-circuit diffusion and GB trapping phenomenon from the energy perspective first.By treating hydrogen diffusion as a random walk from one site toward its neighbor sites, the concept of mean absolute migration energy is introduced as the average absolute segregation energy gap between sites in microstructure A and their neighbor sites in microstructure B: where A and B are the identified microstructures in Section 3.2 that could be GB core, GB surface, or lattice.i and j are all interstitial sites belonging to microstructures A and B, respectively, E seg i and E seg j are the segregation energy of sites i and j , respectively.Only the neighbor pairs ( i and j) with a distance less than the cutoff distance d c of 1.55 Å ( √ 3 / 4 a 0 , which is the distance between the closest octahedral site and tetrahedral site), will contribute to the mean absolute migration energy E mig A , B , n is the total number of qualified pairs.Therefore, we only focus on the migration behavior in the nearest neighbor.For example, given A = lattice, B = lattice, the cutoff distance, and the fact that any octahedral site in the perfect lattice is surrounded by eight tetrahedral sites with an energy gap of 0.4 eV, only the nearest octahedral-tetrahedral pairs will be counted as qualified neighbor pairs.Thus, the E mig lattice , lattice is equal to 0.4 eV.The mean absolute migration energies inside or between each structure are displayed in Table 2 , with a lower value indicating a lower energy barrier for hydrogen diffusion.For a hydrogen atom in the lattice, it must overcome at least E mig lattice , lattice of 0.4 eV to diffuse.In contrast, hydrogen atoms at GB core only need to overcome an average energy E mig GB core , GB core of 0.072 eV, providing a fast channel to facilitate hydrogen transport.The value of E mig GB core , GB core could be also compared to the standard deviation of Peak1 ( σ 1 = 0 .057 eV) given the large overlap between Peak1 and GB core sites.The diffusion inside GB surface or between surface and core is also easier compared to the lattice.These low energy barrier inside GB provides direct evidence of GB short-circuit diffusion with the GB core making the most significant contribution.However, the transport between GB and lattice is characterized by relatively high energy gaps of 0.324 and 0.44 eV.In particular, diffusion between the GB core and lattice can be more challenging than in the lattice.This indicates that GB could provide a short-circuit diffusion channel along its core region while simultaneously preventing hydrogen diffusion perpendicular to it by jumping into the lattice.The collective diffusion behavior in the vicinity of GB is characterized by the energy gap between GB sites and all their neighboring sites E mig GB , all of 0.233 eV.It should be noted that the mean absolute migration energy serves as an indicator to characterize the energy difference between the initial state (at sites belonging to structure A) and the final state (at sites belonging to structure B).Its value reflects the minimum energy required for crossing the energy barrier, in contrast, the energy barrier is the energy difference between the initial state and the saddle point.Typically, the saddle point exhibits a higher energy state than the final state.However, in the nickel-hydrogen system, the energy difference between the saddle point and the final state can be relatively small compared to their absolute values [ 37 , 52 ].Given the large statistical base, we assume that the mean absolute migration energy, even if not linearly scaled, still characterizes the degree of difficulty for a hydrogen atom at structure A to diffuse to structure B. Therefore, by replacing the energy barrier E bar GB and E bar lattice in Eq. ( 8) with E mig GB , all and E mig lattice , lattice , respectively, we can make a rough estimation of the GB diffusion coefficient D GB , which is approximately 640 times larger than the lattice diffusion coefficient D lattice at 300 K.It is worth noting that linear averaging may introduce inaccuracies when applied to Arrhenius behavior.In this context, the aim is to provide a rough qualitative comparison of the overall diffusion properties among different structures.Therefore, further kinetic analysis is necessary to obtain a more accurate prediction.
Thus, molecular dynamics are employed to calculate the GB diffusion coefficient D GB through the MSD analysis in Eq. (10) .10,0 0 0 hydrogen atoms are randomly inserted into the default polycrystal model or single crystal model with the same size, and the systems are maintained at 300, 350, 400, 450, 500, 550, and 600 K to track the average MSD for all hydrogen atoms.100 ns is sufficient to get the converged MSD (shown in Fig. A1 in Appendix) without global redistribution where Eq. ( 10) is equal to Eq. ( 6) .The diffusion coefficients for single crystal D lattice and polycrystal D poly are displayed in Fig. 6 , with the corresponding GB diffusion coefficient D GB derived from linear model in Eq. ( 6) or HS + model in Eq. ( 7) .The polycrystal shows higher effective diffusivities about 170 times at 300 K and 4 times at 600 K than the lattice, and the calcu-Fig.6. Hydrogen diffusion coefficients as a function of temperatures.The diffusivity in single crystal and polycrystal are directly extracted from MSD analysis, while the GB diffusivity is calculated through linear model in Eq. ( 6) or HS + model in Eq. ( 7) .The Arrhenius equation is obtained by linear fitting of the diffusion coefficients at varying temperature.Experimental results for lattice diffusivity [ 25 , 26 ] at varying temperatures and GB diffusivity [25][26][27][28][29][30][31][32][33] at room temperatures are labeled by stars.lated GB diffusivity could be 970 (for linear) and 690 (for HS + ) times at 300 K and 19 (for linear) and 13 (for HS + ) times higher.It proves that the collective effect of GB is accelerating hydrogen diffusion, and the enhancement of three orders of magnitude at 300 K agrees well with the prediction from the mean absolute migration energy analysis above.At the examined region, the linear model and HS + model predict similar values of D .The experimentally obtained D lattice at varying temperatures [ 25 , 26 ] and D GB at 300 K [25][26][27][28][29][30][31][32][33] are displayed as stars.Compared to the consistent results of D lattice in the literature which also agree with our calculation, the D GB show a large fluctuation range from 1 to 10 5 times the lattice diffusivity.By eliminating the effects of other defects in the material system and identifying the GB volume fraction directly from the atomic structure, we have an accurate prediction of GB diffusivity to be 10 3 times at 300 K. We note that the polycrystal system examined in this study consists of grains with an average grain size of 10 nm, which is significantly smaller than the grain size typically observed in coarse-grain or even nanocrystalline materials in experimental studies.The high volume proportion of GBs and triple junctions in such a polycrystal, resulting from the small grain size, could lead to a high calculated hydrogen solubility and diffusion coefficient D poly .The amount of soluted hydrogen can be estimated to scale with the GB volume, and it is important to note that as the grain size increases, a substantial decrease in D poly is expected.In fact, Kirchheim et al. [29] have previously reported the hydrogen diffusivity to be only two orders of magnitude higher in nanocrystalline nickel with a grain size of 100 nm.As the grain size continues to increase, the volume of GBs decreases, and the short-circuit diffusion effect diminishes.With coarse-grain sizes on the micrometer scale, the enhancement of diffusivity from GBs can be largely mitigated, resulting in similar diffusivity between coarse-grain and single crystal materials [25] .Thus, D poly can decrease from two orders of magnitude higher to the same order of magnitude as D lattice , within the range of 100 nm to micrometer-sized grains.The diffusivities at varying temperatures are further fitted using the  where the activation energy barriers E bar , and the pre-exponential factors D 0 for different structures are shown in Table 3 .The activation energy barrier for lattice diffusion 0.484 eV is similar to the value of 0.4 eV [2] from previous experiments in a large temperature range, the slightly higher value could be explained by the ideally perfect lattice in the MD simulations compared to experiments.Actually, by introducing point defects such as vacancy, interstitial, and Frenkel pair into the single crystal, Zhou et al. [67] found the E bar for lattice diffusion can be decreased 0.05 −0.1 eV, which will make our results completely match the experimental value.As for GB diffusion, the linear model and HS + model output a similar value of energy barrier of 0.286 eV.The gap of 0.2 eV between diffusion barriers in the lattice and GB show conformity with the gap of 0.167 eV between E mig GB , all and E mig lattice , lattice in Table 2 .This further proves that the properties derived form spectrum could well capture the hydrogen diffusion behaviors at GBs.

Conclusion
Through a novel calculation route combining atomistic simulations and Voronoi diagram, we were able to identify all the possible interstitial sites in polycrystals consisting of realistic, random GB networks and derive the spectrum of hydrogen segregation at GBs.This spectrum, based on millions of local atomic configurations, revealed a comprehensive picture of hydrogen-GB interactions.Three peaks with distinct energy levels were identified using the Gaussian mixture model and further connected to groups of in-terstitial sites with characteristic geometrical fingerprints, i.e., the GB core sites and GB surface sites.The hierarchical arrangement between the GB core, GB surface, and lattice was found to facilitate both short-circuit diffusion and trapping, with the low barrier inside the same structure facilitating hydrogen transport while the high barrier between different structures inhabiting hydrogen diffusion.Through kinetics analysis, we further extracted the hydrogen diffusivity in polycrystals at varying temperatures, where the hydrogen diffusion coefficient at GB was found to be three orders of magnitude higher than in the lattice.Finally, we incorporated the spectrum with a thermodynamics framework to predict the equilibrium hydrogen concentration at GB without neglecting the diverse nature of GBs.It is important to acknowledge that in a real polycrystal, multiple defects such as vacancies, dislocations, solute atoms, and strain fields are present.These defects can have both individual and collective effects on the energy spectrum, thereby influencing the segregation tendencies and enrichment of hydrogen.For example, the formation of vacancies at GBs can significantly decrease the segregation energy and facilitate the clustering of hydrogen.Furthermore, the presence of strain fields can lead to a global shift in the energy spectrum, and the extent of this effect highly depends on the loading mode.Therefore, it is crucial to consider the influence of these multiple defects in order to obtain a comprehensive understanding of hydrogen segregation and enrichment in polycrystalline materials.In our forthcoming study, we will systematically elucidate the effect of strain fields on the energy spectrum, as well as on the diffusion and trapping behaviors of hydrogen.

Fig. 1 .
Fig. 1.(a) Default (20 nm) 3 nickel polycrystal.(b) Map of potential interstitial hydrogen trapping sites, colored by its corresponding segregation energy.(c) Relationship between hydrogen segregation energy E seg GB and distance from GB.Only the closest sites within a cutoff distance of 2.2 Å (blue) are taken as interstitial sites at GB, the remaining sites (pink) are regarded as in lattice.

Fig. 2 .
Fig. 2. Segregation energy spectrum of interstitial hydrogen at GB. (a) Probability density distribution of hydrogen segregation energies at 491,799 (model 1), 1,347,330 (model 2), 428,936 (model 3), and 203,052 (model 4) interstitial GB sites.The decomposed three Gaussian distributions and their sum are denoted by red, green, blue, and black dash curves, respectively.The segregation energies of octahedral and tetrahedral sites in the lattice are indicated by the vertical dash line.(b)Cumulative probability density distribution of segregation energies in the four models.The red symbols denote hydrogen segregation energy in several typical CSL GB models from previous DFT calculations[ 37 , 52 , 35 , 70 , 71 ]  or EAM potentials[ 40 , 52 ].

Fig. 3 .
Fig. 3. (a) Divided spectrum using the mean value of Peak1.Only the left part ( E seg GB < −0 .205 eV , group A1) of the spectrum is visualized according to their energy level in (a1), while the whole spectrum is visualized in (a2) but sites belonging to the right part ( E seg GB > −0 .205 eV , group A2) are all colored red.(b) Divided spectrum into GB core sites (blue) and GB surface sites (olive).Only GB core sites (group B1) are visualized in (b1), while the GB surface sites (group B2) are all colored olive in (b2).

Fig. 5 .
Fig. 5. (a) Equilibrium atomic hydrogen occupation at GB c GB as a function of lattice hydrogen occupation c l at T = 300, 450, and 543 K.The prediction based on the full discrete spectrum in Eq. (13) , the prediction based on Gaussian mixture model in Eq. (14) , and the contribution of Peak1 are illustrated by the solid, dash, and dot lines, respectively.(b) Effective single-value segregation energy E seg GB by substituting c GB and c l from (a) into Eq.(11) .(c1 −c9) Detailed equilibrium hydrogen segregation state at GB with varying lattice hydrogen occupations and temperatures.The GB hydrogen occupations c GB and contribution from Peak1 is denoted by solid and dot curves, respectively, with green, yellow, and orange colors.The spectrum is filled with blue and the Peak1 is indicated by the red dash curve.

Table 1
Parameters of Gaussian mixture distribution for the different polycrystal models, where u i , σ i , and w i are the mean, standard deviation, and weight of each peak.

Table 2
Mean absolute migration energy between different structures.

Table 3
Activation energy barrier E bar and the pre-exponential factor D 0 for different structures.