I. INTRODUCTION
Metals used to store tritium are subjected to extensive interactions with helium atoms that are produced by transmutation reactions from energetic tritium fusion neutrons.Reference Yu, Zhao, Zhang, Yang and Chu1 Because of extremely low solubility of helium in materials, helium atoms tend to be trapped at defects, such as vacancies, dislocations, and grain boundaries, and can consequently form helium nanobubbles.Reference Yang, Zu, Xiao, Gao, Heinisch, Kurtz and Liu2–Reference Vassen, Trinkaus and Jung5 The formation of helium bubbles will lead to high-temperature intergranular embrittlement,Reference Barnes6 which can significantly degrade the mechanical properties of materials. At the same time, the release of helium atoms from tritium-storing materials will decrease the purity of tritium and will also influence the plasma performance in a fusion environment. Therefore, the retention and release behavior of helium is an important issue for tritium-storing materials. As such, the helium release from Ti, Pd, and Ni et al.Reference Lasser7–Reference Shanahan and Hölder13 has received great attention. By analysis similar to that for dislocation loop punching, it was found that further bubble growth ceases by one of two mechanisms, (i) the equilibration of the chemical potentials for gas atoms in the bubble and in interstitial solution or (ii) interbubble fracture when the average ligament stress exceeds the theoretical tensile strength of the materialReference Wolfer14, Reference Wolfer15 as the helium content in the host metal reaches critical dimensions.Reference Bowman and Attalla16 Previously, two theoretical models were suggested to explain the release behavior of helium atoms from metal materials. One is the percolation theoryReference Camp17 proposed by Camp, in which helium release is correlative with the percolation of helium in the tritide. As the helium concentration reaches to a critical value, the percolation network forms and helium can release along this network. The other is the model proposed by Schober et al.Reference Schober, Trinkaus and Lasser18 in which dislocation loops and networks can accelerate helium release. These studies provide important insights into helium release from materials.
Recently, Snow and BrewerReference Snow and Brewer11 have studied the helium release and microstructural changes in Er(D,T)2−x3Hex films, and it was found that helium is released only from the near-surface region. There exists a helium-bubble-free zone, or denuded zone, with a thickness of approximately 15 nm. Using a set of rate equations, a continuum-scale evolutionary modelReference Cowgill19 successfully predicted the formation of a helium-bubble-free zone, but the thickness of the denuded zone is different from the experimental measurements. Whether or not atomic helium exists in the denuded zone is unclear both experimentally and theoretically. A detailed investigation of the surface influence on helium release has also not been performed, and the atomic-level release mechanisms of helium nanobubbles are largely unknown, because the helium release process cannot be easily attained by experiments. Therefore, a computational approach has been developed as an important method to study the release of helium atoms and to provide the atomic-level mechanisms associated with the motion, diffusion, and release of helium atoms.
It is known that the helium retention in metal tritides depends on stoichiometry, helium distribution, and aging effects. In the present work, we have chosen to investigate the evolution of helium nanobubbles in a palladium crystal in terms of the variation of their depth from the surface and the release mechanisms involved using molecular dynamic (MD) simulations. The tendency of 3He atoms to form nanobubbles in Pd tritide has been clearly demonstrated by TEM observations, even for low 3He contents,Reference Thiébaut, Douilly, Contreras, Limacher, Paul-Boncour, Décamps and Percheron-Guégan20 and Pd is an excellent material for tritium storage that retains the helium generated in the matrix by radioactive decay for at least 19 years.Reference Emig, Garza, Christensen, Coronado and Souers21 Therefore, the present simulation results are used to explain the physical origins of the helium nanobubble release in Pd and to understand the formation of the denuded zone.
II. MODEL AND METHODS
The modified analytic embedded atom method (MAEAM)Reference Hu, Naka and Yamane22–Reference Xiao and Hu24 is used, employing the MorseReference Baskes and Melius25 potential to describe the interactions between Pd–Pd and Pd–He as well as the Lennard–JonesReference Johnson26 potential for He–He interactions. A Pd slab consisting of 12a 0 × 12a 0 × 25a 0 unit cells is modeled, having periodic boundary conditions along the (100) and (010) directions and a free-boundary condition along the (001) direction, where a 0 is the lattice parameter (0.386 nm for palladium). Within this model a spherical void with a 2 nm diameter is formed at a given depth below the surface and filled with helium atoms in the void region. It is well known that the increase of helium atoms in a bubble can promote the expansion of the bubble in the void region,Reference Lasser7 and thus models containing 321, 587, 723, and 795 He atoms in the bubble are considered in the simulations. For each case, the different geometry and configuration of the bubble at different depths are studied at 300 K. Two different algorithms were used to check their sensitivities on the helium motions within the bubble: (i) a constant pressure and constant temperature (NPT) simulation is performed in a bulk system to calculate the lattice parameter at 300 K, and then the model is simulated under the constant volume and constant temperature (NVT) condition and (ii) all the simulations are directly performed with the NVT ensemble. It was determined that there are no obvious differences in the helium motions within the bubbles in the two models.
III. RESULTS AND DISCUSSION
A. Correlation between bubble’s shape and depth
It is critical to understand the evolution and release behavior of helium nanobubbles from the crystalline surface and the corresponding mechanisms. To study the effect of crystalline surface on helium bubbles, we simulated several Pd–He systems containing one helium nanobubble located at different depths, h, from the Pd surface. Here h is defined as the distance between the free surface and the center of the helium bubble in the host metal. The thickness of a thin film can be described by t, the distance between the surface and the top (the part of the bubble closest to the surface) of the helium bubble, and thus, we have h = t + r He, where r He is the radius of the helium bubble considered. The transformation of the bubble’s shape as a function of its depth h is shown in Fig. 1, with 723 helium atoms in the bubble. When the helium bubble is located at depths h below the surface ranging from 4.63 to 2.31 nm, the shape of the bubble is approximately spherical, and it maintains a symmetric structure for a long time at the initial temperature, similar to its initial configuration. The interface between the bubble and palladium atoms is pliable, as shown in Fig. 1(a). It is of interest to note that, in Figs. 1(b) and 1(c), the shape of the helium bubble starts to distort, and it is obviously elongated toward the surface, which is due to the influence of the free surface. Once the balance between the internal pressure of the bubble and the surface strain is obtained, a longer simulation does not result in further deformation of the bubble. When the depth h decreases, i.e., when the bubble is closer to the surface, helium atoms suddenly release from the palladium materials, along with strong surface deformation, as shown in Figs. 1(d) and 1(e).
The atomic configurations of the bubble with 723 helium atoms at different depth h: (a) from 2.31 nm to 4.63 nm, (b) 1.93 nm, (c) 1.54 nm; (d, e) show the helium eruption from the surface and the release behavior of helium atoms at 10 and 600 ps at a depth of 1.35 nm. Here the helium atoms are described by light spheres and the Pd atoms by the dark spheres.
B. Dynamic process of helium release
The release behavior of the helium atoms is described as following. During the release process, the helium atoms rapidly propagate toward the surface, pushing the palladium atoms out of their pathway, and creating a channel for further helium release. This process induces the crystalline surface to deform and eventually rupture. As a result, a cavity is formed with a channel connecting the bubble to the surface, and the helium atoms can easily escape from the material along the channel one by one through the cavity. Thus, the formation of the surface cavity accelerates the helium release from the near-surface region.
It is of interest to note that only helium atoms are released during this process, without any palladium release being observed. Thus, the possibility for using tritium-storing materials in a fusion environment that can release He without releasing other impurities may be an important design issue affecting plasma performance.
Figure 2(a) shows the number of released helium atoms from the surface as a function of time. It can be seen clearly that the release is a nearly instantaneous process, with almost all of the helium atoms being released from the surface within a time of 75 ps. To further describe the helium release process effectively, we have determined the release rate of helium as a function of time, which is shown in Fig. 2(b). The release rate appears to have three stages: (i) an accelerating release rate from 0 to about 10 ps, which corresponds to the high-pressure release of helium bubble, (ii) a decelerating release rate from 10 to 35 ps, which may be associated with the decrease of the pressure in the bubble, and (iii) a constant, relatively low release rate for the time longer than 35 ps (the pressure in the bubble becomes small, but is large enough to push helium atoms releasing along the channel, as discussed previously). The helium release also results in the reconstruction of the crystalline surface [Fig. 1(e)].
(a) The evolution of helium amounts released as a function of time, where the small figure shows the change trend of release amounts from 0 to 80 ps and (b) the evolution of release speed, where the small figure indicates the trend from 0 to 40 ps.
The different configurations of the surface layers with time are shown in Fig. 3(a). A small cavity at the crystalline surface, accompanied by the helium release, initially appears at 2 ps, and then it increases gradually with time [Fig. 3(b)]. After it reaches a maximum size, at about 10 ps, the cavity region starts to shrink and reaches a stable size as the time advances beyond about 100 ps. The cavity opening can be considered approximately a circle with a specific diameter that is defined as an average diameter of the opening at the surface and can be estimated from the atomic configurations. The variation of the cavity diameter with time is shown in Fig. 3(b). Initially, the cavity diameter grows with time, then decreases with increasing time from 10 to about 30 ps, and finally remains a constant for the further simulations. Moreover, the helium release also induces the palladium atoms on the surface to move away from the cavity center, forming some protuberances around the cavity as the time increases, as shown in Fig. 4(a). The atomic configurations clearly demonstrate that the protuberances become more evident from 0 to 10 ps, but the protuberance is gradually reduced from 10 to 20 ps, which is consistent with the increase and decrease in the cavity size with increasing time, as observed previously. For a more quantitative description of the protuberance variation, we have defined a parameter ΔD:

(a) The transformation of cavity and (b) the evolution of cavity size as a function of time, where the small figure shows the trend from 0 to 17 ps.
(a) The transformation, deformation and reconstruction of the surface atoms as a function of time, where the base surface is indicated by lines and (b) the evolution of the ΔD, where the small figure shows the trend from 0 to 18 ps. Here, ΔD describes the degree of surface protuberances.
where Δdi is the vertical distance between the palladium atom of the protuberance and the base-surface, and M is the total number of palladium atoms in the protuberance region on the surface plane. Physically, ΔD describes the extent of protuberances, and the variation of ΔD with time indicates the evolution of the protuberance as time advances. The time variation of ΔD is shown in Fig. 4(b), and its trend is generally similar to that of the cavity size, i.e., rapidly increasing during the first few picoseconds, and then decreasing with time to a constant value. The helium atoms escape one by one from the surface, rather than in the small cluster forms.
C. Helium release mechanism
To further investigate the surface effects on helium release, the volume of the helium bubble was calculated at different depths from 4.63 to 1.54 nm below the surface using an algorithm for Convex HullReference Barber, Dobkin and Huhdanpaa27 that determines a set of facets and a set adjacency lists of an irregular object, from which its volume can be estimated. Figure 5 shows the evolution of the bubble volume as a function of depth, with different amounts of helium atoms in the bubble. From these curves, it is clear that the transformation of the bubble shows three stages as the bubble approaches the surface. At first, when the bubble is far away from the surface its volume is approximately constant, and then it starts expanding as it gets closer to the surface. There is a critical depth, h c, at which the bubble begins to expand and deform that depends on the amount of helium in the bubble, and increases with increasing bubble size, as indicated by a dotted line in Fig. 5. Finally, the bubble cracks, and most of the helium atoms are released from the Pd surface. If the depth h of the helium bubble is larger than the critical depth h c, the helium atoms remain as a bubble within the matrix, even if the simulation has been extended for a much longer time.
Variation of the volume of He bubbles as a function of their depth from the surface, h, for helium atom clusters containing (a) 587, (b) 723, and (c) 795 He atoms, respectively. The dashed line intersecting curves (a), (b), and (c) indicates the “cracking” of helium bubbles.
To release helium atoms, it is expected that the metal film above the bubble will be ruptured when the gas pressure P of the helium bubbles exceeds its tensile strength. The virial functionReference Selyutin28–30 is used to calculate the pressure P of the helium bubbles in a Pd crystal:
where N is the number of the helium atoms, k B is the Boltzmann parameter, T is the ambient temperature, V is the corresponding volume of the helium bubble, and ri and
are the distance and the mutual force of the helium atoms, respectively. Figure 6 shows the gas pressures of the bubbles (with 321, 587, 723, and 795 helium atoms) at different depths h, where the solid line indicates the strength of metal film. For each bubble, the gas pressure remains essentially constant when it is far away from the surface and the strength of metal film above the bubble is high enough to retain the helium bubble in the host lattice. For example, the strength of the metal film at h = 3 nm is larger than the pressure of the helium bubble with 795 helium atoms (d) in Fig. 6. Therefore, no helium atoms are released, and the gas pressure and volume of the bubble change little. In general, the strength of the film decreases with increasing its thickness. When the helium bubble is located at the critical depth h c near the surface, the film cannot burden the stress from the helium bubble. At this time, the initial gas pressure is strong enough to induce some fractures, or even break the film, which leads to the deformation and swelling of the bubble. Then a cavity initiates on the surface, microcracks appear near the bubble surface, and helium atoms release instantaneously and erupt like a volcano, which represents the microscopic mechanisms of helium nanobubble release from the Pd crystalline surface. The critical gas pressure P c at the depth h c increases with the increase of the helium atoms (or helium-to-vacancy ratio) in the bubble.
The relationship between the gas pressure of He bubbles and their depth h, with (a) 795, (b) 723, (c) 587, and (d) 321 helium atoms, respectively, where the black solid line shows the host-metal film tensile strength as a function of depth h.
The present simulations indicate that the film tensile strength σ is an important parameter for determining helium release. The critical gas pressure P c of the helium bubble is equal to the maximum tensile strength σ of the host-metal film with a thickness of t. The tensile strength σ of the host-metal film can be evaluated using the gas pressure P c of the helium bubbles and the film thickness. An exponential function can be used to describe the relation between the film strength σ and thickness t,
where σb and t 0 are two constants, σm is the theoretical (ideal) tensile strength when the film thickness t is unlimited. The three parameters can be obtained by fitting Eq. (3) to the strength σ and the thickness t of the host-materials film, which are equal to the P c and h c–r He of the helium bubbles in the present simulations, respectively, and they are determined to be σm = 34.5 GPa, σb = 42.46 GPa, and t 0 = 1.03 nm. The fitted results are also shown in Fig. 6, as indicated by a solid line. Using first-principles calculation, Černý and Pokluda calculated the theoretical tensile strength of 10 cubic metal crystals,Reference Černý and Pokluda31 and the theoretical tensile strengths for Ni, Pt, and Au were determined to be 35.2, 34.1, and 44.5 GPa, respectively. The theoretical tensile strength of Pd obtained from the present simulations is in good agreement with their prediction.
As described in the introduction, both experimental observationReference Snow and Brewer11 and theoretical predictionReference Cowgill19 have found a helium-bubble-free zone near the crystalline surface. The similar surface helium release phenomenon has also been observed in our atomic-scale simulations. There is a near-surface zone where all the bubbles will release their helium instantly, the so-called denuded zone. We have found that the helium bubble is retained in the host metal when it is located outside the denuded zone, whereas the bubble will expand rapidly and the helium release occurs spontaneously as long as the bubble is within the zone. This suggests that all the helium bubbles are unstable in the denuded zone, which determines a helium-bubble-free zone. The denuded zone thickness, λd, is calculated from the bottom of the helium bubble to the surface. Mathematically, λd is equal to the value of the critical depth hc plus the helium bubble’s radius r He, as shown in Fig. 5. The thickness λd of the denuded zone correlates with the bubble size, or the amount of helium atoms in the bubble. With an increase in the number of helium atoms or the bubble size, the critical depth of the bubble deformation will increase, and the thickness of the denuded zone becomes larger, as seen in Fig. 5. Considering a bubble with 723 helium atoms as an example, the radius r He of the helium bubble is about 1.0 nm and is almost a constant before it releases its helium during the simulation. In this case, the critical depth h c is 2.3 nm in the host metal and thus the thickness λd of the denuded zone is determined to be 3.3 nm. In general, the thickness of the denuded zone depends mainly on the size of the helium bubble.
In our molecular dynamics simulations, the thickness of the denuded zone is several nanometers in size (about 3 nm typically, as mentioned previously), and it is different for different bubble sizes. As predicted by Cowgill, using a rate equation model,Reference Cowgill19 the effective helium bubble escape depth in Pd–T0.65 is about 2 to 5 nm, depending on the age of the material. In Er(D,T)2−x3Hex films, Snow and BrewerReference Snow and Brewer11 have found experimentally that the thickness of the denuded zone is approximately 15 nm. It is clear that there are some differences in the thicknesses of denuded zones obtained with different methods, which may be correlated with the properties of the host materials used. The helium bubbles are difficult to release from the host materials with excellent helium retention ability, which indicates that the thickness of the bubble denuded zone will be thinner, and vice versa. The present study is aimed at understanding the helium bubble release in a palladium lattice, while this experiment was carried out in Er(D,T)2−x3Hex films. Although the bubble denuded zones are both very thin, they are different quantitatively. Because of the different materials used, the effects of lattice defects, such as vacancies, dislocations, and grain boundaries, are important in these host materials. Furthermore, the effects of oxidation near the surface may affect the ability to determine the thickness of the bubble denuded zone experimentally. However, both Cowgill’s model and the present atomic simulations have successfully predicted that there exists a denuded zone near the surface, which provides a physical explanation for the experimental observation.
As for Cowgill’s continuum-scale evolutionary model,Reference Cowgill19 a set of rate equations were fitted to many experimental data and were successfully applied to study the helium nanobubble release for aging metal tritides, but atomistic processes have not been considered. Snow and BrewerReference Snow and Brewer11 assumed that it is possible that the bubble-free region has a high concentration of atomic helium, but in their experiments it may be rather difficult to measure the residual helium atoms. We also find that there are still some atomic helium residues in the denuded zone after the helium releasing, as clearly shown in Fig. 1(e), but the concentration of atomic helium is very low. It is clear that the present simulations have successfully explained the major features of the experimental phenomena.
IV. CONCLUSION
The microscopic mechanisms of helium release from a palladium crystal are investigated using MD methods. A denuded zone has been observed in the modeling results, which supports the recent experimental observation. Outside the denuded region, further away from the surface, the helium bubbles are retained in the host metal, but within the denuded zone a bubble can expand rapidly and finally “crack,” creating a channel for helium to propagate to the crystal surface. Almost all the helium atoms in the bubble escape the material, and just a few remain in the denuded zone. The present simulations provide a physical explanation for why the helium is released only from the near-surface region rather than through bulk diffusion. The helium release is an instantaneous process, with the helium erupting from the surface, creating surface deformation and the nucleation of a cavity on the surface. The cavity formation can further accelerate the helium release from the surface. The present studies demonstrate that atomic-level simulations provide an important method to study the dynamic process of helium release as well as determining the effects of the actual material microstructures on helium behavior in metals. The relation between the tensile strength and thickness of metal film has also been established.
ACKNOWLEDGMENTS
This work is financially supported by the National Natural Science Foundation of China (50871038), China Institute of Atomic Energy, and the High Performance Computing Center of the Hunan University. F. Gao and H.L. Heinisch are grateful to the support from the United States Department of Fusion Energy Science under Contract No. DE-AC06-76RLO 1830.

