{'fraction_non_alphanumeric': 0.03383242153757558, 'fraction_numerical': 0.018283904405413186, 'mean_word_length': 4.990944372574385, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hydrogen Segregation in Palladium and the Combined Effects of Temperature and Defects on Mechanical Properties

Atomistic calculations were carried out to investigate the mechanical properties of Pd crystals as a combined function of structural defects, hydrogen concentration and high temperature. These factors are found to individually induce degradation in the mechanical strength of Pd in a monotonous manner. In addition, defects such as vacancies and grain boundaries could provide a driving force for hydrogen segregation, thus enhance the tendency for their trapping. The simulations show that hydrogen maintains the highest localization at grain boundaries at ambient temperatures. This finding correlates well with the experimental observation that hydrogen embrittlement is more frequently observed around room temperature. The strength-limiting mechanism of mechanical failures induced by hydrogen is also discussed, which supports the hydrogen-enhanced localized plasticity theorem.


INTRODUCTION
High activities of hydrogen in palladium have provided practical implications in many areas of materials science and chemistry, including catalysis 1-4 , composite membranes [5][6][7][8] , nanosensors [9][10][11] , hydrogen storage and fuel cells [12][13][14] . In fact, a unique solubility of hydrogen in palladium was first reported as early as the mid-nineteenth century 15,16 . However, it is believed that undesirable hydrogen absorption in crystals could substantially degrade their mechanical endurance, thus challenges the stability and performance of devices 17,18 . Therefore, a fundamental understanding of hydrogen effect is necessary to help improve the materials' design and production.
It was accepted that elastic response of hydrogenated metals and alloys changes due to the lattice relaxation and local hydrogen excitation [19][20][21][22] . Several major hydrogen embrittlement mechanisms have been established: (i) hydrogen-enhanced localized plasticity 23 ; (ii) stressinduced hydride formation 24 ; and (iii) lattice, grain boundary (GB) and interface decohesion 25 . In general, diffusion of hydrogen atoms in crystals is enabled by their occupations of interstitial sites until they are trapped in open volumes such as vacancies and GBs. Additionally, GBs could also act as the transport channel of impurity species. In fact, presence of vacancies and GBs in metals and alloys could substantially alter their mechanical response, with or without the localization of impurity elements [26][27][28] . However, various studies have reported that hydrogen propagates GB cracks and reduces the ductility of metals and alloys [29][30][31][32][33][34][35][36][37][38][39][40][41] . It was suggested that hydrogen causes embrittlement by promoting the formation of microvoids and decreasing the cohesive strength of GBs.
Lately, atomistic simulations became convenient techniques for investigation of electronic structure and geometry of nanomaterials. In particular, Daw and Baskes developed the embedded-atom method (EAM), which was shown to successfully describe the properties of metal and alloy systems 42 . For instance, the method has been used by Zhong et al 43 to study the hydrogen-loaded palladium under tensile stress and they reported about the pre-melting disorder caused by hydrogen. In addition, first-principles calculations based on density functional theory have also been applied to examine the electronic origin of hydrogen-vacancy and hydrogen-GB complexes in several metal and alloys [44][45][46][47] . In this work, we use the classical molecular dynamics (MD) with embedded-atom method to direct theoretical investigations on hydrogen segregation in Pd using different crystal models (Table I)

COMPUTATIONAL METHODS
We used the molecular dynamics method to compute equilibrium and mechanical properties of the H-Pd system. The nuclear motion of particles in MD simulation technique is based on classical mechanics 48 , in which the interactions between atoms are characterized by a force field.
In this work, we used the interaction potential from the EAM, which was developed by Daw and Baskes 42,49 . This semi-empirical, many-body model describes the total energy of a metal atom as the energy obtained by embedding that atom into the local electron density of its atomic neighbors. Unlike pair potential models, in which the energy is just the sum over pair bonds, the total energy in EAM consists of an embedding energy term and the electrostatic interaction term.
where the first term stands for the embedding energy and ρ is the background atomic electron density. The second term is related to the contribution from all cross-pair atomic interactions.
The EAM approach was reportedly capable of solving many problems of interests in metals and metal alloys, including defects, impurities, surface, fracture, etc. 50 . More specifically, it was successfully applied to metal hydride problems and hydrogen embrittlement phenomena [51][52][53][54] . We used in this work the EAM hydrogen-palladium interatomic potential developed by Zhou et al. 55 A supercell model of 10,000 atoms was used to simulate the face-centered cubic (FCC) palladium (single crystal). In case of a bicrystal (i.e. FCC Pd bulk containing a grain boundary), we studied specifically a ∑5(210) grain boundary that separates the supercell into two differentlyoriented grains of 5000 Pd atoms on each ( Figure 1). This supercell size (10,000 Pd atoms) yields the convergence of GB surface energy as low as 0.1 meV/A 2 . As a high-angle grain boundary, this GB plane offers considerable free volume for impurity transports and segregations and therefore could be potentially vulnerable to impurity-induced embrittlement.
Initial H-Pd configurations were prepared by randomly inserting H into interstitial or empty space inside Pd cells (single crystal, vacancy-containing and bicrystal supercells).
Afterwards, the modulus and strength were studied by means of tensile test simulations. The MD time-step is chosen as 0.5 fs and the periodic boundary conditions are imposed in all x, y and z directions. At each temperature T and hydrogen content x H , the reference zero-strain state was firstly obtained from constant temperature and stress ensembles by optimizing the volume of the supercell. The uniaxial tension is then applied in the z-direction (which is [210] for all supercell models, and normal to the GB plane in case of the bicrystal), and is increased gradually in steps of 50 MPa, followed by adequate relaxation to optimum volume (at that external pressure). The tensile strength σ TS is defined as the maximum tensile stress at which the system is sustained as stable. The tensile modulus in the applied tension direction is calculated by fitting the linear stress-strain relation 56 within a small range of strain (up to 0.5%): where σ z and ε z correspond to stress and resulting strain along applied tension direction (z-axis). Figure 2 shows the stress-strain relationship for a bicrystal at x H =0.1 and T=300K, from which tensile modulus and tensile strength can be derived. When the applied stress is increased in the z-direction, boundaries of simulation box are allowed to shrink or dilate in order to maintain zero stress in x and y dimensions. In other words, in this method, we control the stress and calculate the subsequent strains (positive along z and negative along x and y), induced by that uniaxial stress. Having knowledge of these domain changes, we were also able to conduct a separate simulation, in which the strains are totally controlled and then stress could be obtained based on kinetic energy and interatomic interactions, using the canonical ensemble (NVT). This simulation approach (strain control) is actually a more precise reproduction of the tensile test experiment. However, the stress-strain curves obtained by theses two approaches (stress control and strain control) are virtually similar to each other ( Figure 2). In the later method, the stress drops when the material is pulled up to a certain stretch, which indicates the breaking of some interatomic bonds. The ultimate tensile strength, or the maximum tensile stress at which the material can sustain stable deformation, is identical in both approaches. In this work, we adopted the stress-control procedure.

A. Hydrogen segregation at grain boundary
In the single crystal, the equilibrium distribution of hydrogen is expectedly uniform throughout the supercell. However, at the vicinity of grain boundary, the hydrogen segregation could be more complicated. At low temperature, hydrogen atoms have low diffusion coefficient.
Upon temperature elevation, the diffusion would become active after H gains sufficient kinetic energy. Our simulation showed that GB might provide a great gradient for hydrogen absorption.
At high temperatures (above 400K), the equilibrium distribution of hydrogen is reached quickly in the bicrystal and tension also facilitates the hydrogen segregation at grain boundary.
However, we found that, around room temperature, the saturation concentration of H at GB is relatively high. In addition, the H diffusion is quite active and far from the steady state, even after the equilibrium volume has been reached. Therefore, in order to obtain final H-Pd bicrystal configurations, we firstly heat the initial H-Pd systems up to 300K and maintain at this condition until the saturation of hydrogen at grain boundary was fully reached (approximately after lengthy 10 7 time-steps). The structures used for tensile test then will be obtained from this configuration by heating (or cooling) and then sufficiently equilibrating at desired temperatures. Figure 3a shows the final hydrogen distribution along z dimension at different temperatures (for the case of x H = 0.06). Hydrogen distribution is monitored by dividing simulation box into 10 layers along z-direction, in which layers number one and six contain the GBs (showing high peaks in H content). We saw that between 100 and 300K, two GB-layers are able to trap totally up to 90% of total hydrogen in the bicrystal cell. Interestingly, H concentration at GBs reduces significantly after transition from 300 to 400K. The H contents at GB layers are also displayed in Figure 3b. H distribution is high and consistent in 100-300K and decreases above 300K (visualization of H distribution with x H =0.06 is presented in Figure 4). In Table II and Table III we present the tensile modulus E [210] and tensile strength σ [210] at different T and x H for the single crystal; and the data is visualized in Figure 5. The maps reveal a monotonous decrease for both E and σ with respect to elevated T, which agrees with common trends that metallic materials soften at high temperature or reduced pressure 59,60 . Also, the same behavior was observed for the increase in hydrogen concentration.  56 . This discrepancy between ∑5(100) and ∑5(210) could be explained as the geometry factor of each GB type, in which a portion of atoms situated at ∑5(100) GB in FCC crystals has short interatomic distance and overweighed the atoms with expanded distances in their contribution to elastic response, due to the anharmonicity of interatomic potential 56 .
Meanwhile, the large grain mismatch and big free volume in the ∑5(210) GB, instead, can result in material softening.
In Figure 6 and this suggests one to fit the data using σ T, x = σ ! T + σ ! x + σ ! Tx + σ ! and E T, x = e ! T + e ! x + e ! T + e ! for tensile strength and tensile modulus, respectively, as functions of temperature T and hydrogen content x H . These fitting coefficients were tabulated in Table IV. To make a specific comparison, in Figure 7 we presented their stress-strain behaviors at 300K and x H = 0.1, in which the curves for single crystal and vacancy crystal are quite compatible. The bicrystal has a lower tensile modulus (57.2 GPa), compared to the single (73.5 GPa) and vacancy ones (70.5 GPa). There is also a big degradation (34%) from that of the H-clean single crystal (87 GPa). Interestingly, the bicrystal has the same tensile strength with the other two, regardless of being softer. In other words, the bicrystal is more ductile, but still comes to unstable state at the same stress compared to the others, under high H absorption (x H = 0.1 cause a 40% reduction in tensile strength). The explanation for a higher ductility in bicrystal is that the mechanical degradation in a single crystal occurs throughout the whole cell, as opposed to one specific location (GB in the bicrystal cell).
In Figure 8 we display the atomic configuration of defective Pd atoms in different crystals, respectively, with x H = 0.1 at zero tress (Figure 8a-c) and corresponding maximum tensile strength (Figure 8d-f). The concentration of "defective" Pd induced by vacancy and hydrogen binding is pretty high and this can be the reason why the vacancy-containing crystal has a lower tensile strength compared to the bicrystal and single crystal. We also observed the formation of dislocations (stacking fault) at grain boundary at high external tensile stress (Figure 8f). This specific atomic plane was detected to be (1 1 -1) and belongs to the general family of closepacked slip planes in FCC metals. This emission of dislocations could infer about hydrogenenhanced localized plasticity (HELP) mechanism 23 and that the failures at stress above tensile strength are of plastic nature rather than brittle. Also, it aligns with several experimental observations that fracture process occurs in the vicinity of grain boundary by highly localized plasticity, as opposed to embrittlement along the grain boundary 23,61 .

IV. CONCLUSIONS
Grain boundaries are easily exposed to hydrogen absorption due to high segregation energy, which could accommodate a significant portion of hydrogen in nanocrystalline materials.