Impact of host phonons on interstitial diffusion.

The net effect of host phonons on interstitial diffusion has remained as a fundamental knowledge gap in our current theories since the motions of the host atoms and interstitials were coupled in these theories. Here we study this effect through molecular dynamics simulations of hydrogen diffusion in palladium, in which the motions can be decoupled through pinning the host atoms. Mathematically this decoupling corresponds to expanding the total diffusion coefficient into a Taylor series, which separates the phonon contribution from the intrinsic interstitial jumping. Our results clearly show that palladium phonons significantly promote hydrogen diffusion. The phonon contribution, being linear with temperature at high temperatures and exponential at low temperatures, is fitted with Brownian motion model. The total diffusion of interstitials can be understood as the intrinsic interstitial jumping in a pinned host plus phonon-induced Brownian diffusion. The generality of our findings is validated by examining the motion of lithium in manganese oxide and carbon in iron.

Lattice vibrations or phonons play a major role in many physical properties, such as thermal and electrical conductivity, of condensed matter systems. Their effect on mass transport, i.e., atomic diffusion, however, is not well understood despite studies for decades. It is well known that the diffusion in solids occurs via atomic jumping and follows the Arrhenius equation. In particular, based on the established random walk theory[1,2] the diffusion coefficient for interstitials can be written aswhere is the activation enthalpy for the jump, is the Boltzmann constant, and is a constant virtually independent of temperature T and combines the effects of solid structures, interstitial vibration frequency, and an entropy change term associated with the jump. In the early versions[2,3] of the random walk theory, the effect of host phonons was not considered[4]. Later, a number of classical models[4-7] replaced the vibration frequency of the diffusing atoms with an effective frequency which involves the many body effect from the lattice phonons. The phonon effects on atomic diffusion were also discussed[8,9] based on the quantum theory[10]. In all these theoretical studies, however, the motions of the host atoms and the jumping atom are coupled and hence the net effect or importance of phonons is not clear. In other words, although one can compute the diffusion coefficients with lattice phonons included based on these theories[11], it is not easy to know qualitatively whether lattice vibrations positively or negatively contribute to the diffusion and quantitatively how significant the contribution is. This represents a knowledge gap for our fundamental understanding of solid diffusion, although the random walk theory has been established for decades.

The advance of our computational capacity and relevant algorithms allows us to study the above problem via simulations. Here we address the effect of lattice phonons on interstitial diffusion, using hydrogen (H) in face-centered-cubic (fcc) and amorphous palladium (Pd) as an example, based on molecular dynamics simulations. We consider the amorphous phase here because of its very different structural orders as compared to traditional crystalline solids. Hydrogen in solids has a long history in the interest of scientists[12,13], and metal hydrides for hydrogen storage have seen extensive research activities over the past decades[14]. Pd is one of the first metals of which the capacity for H absorption were discovered[15] and is still widely applied in the hydrogen energy sector[16,17], for example, as effective hydrogenation catalysts, hydrogen purification filters and hydrogen storage media. We validate the generality of our findings by examining two technically important interstitial systems, namely, carbon (C) in iron (Fe) and lithium (Li) in manganese oxide ().

To investigate the amorphous phase, a Pd-H system with 4000 Pd atoms and 200 H atoms was well liquefied at 2100 K and quenched to 300 K at 10 K/s, and the resulting structure was confirmed to be amorphous based on its Pd-Pd pair distribution function. For the crystalline system, we randomly put 200 H atoms at the octahedral interstitial sites of an fcc structure of 4000 Pd atoms. We relaxed the obtained crystalline structures at various temperatures to obtain the equilibrium system volume. These computations were carried out using the NPT (constant atom number, pressure, and temperature) ensemble by setting the normal pressures to zero. Using the NVT (constant atom number, volume, and temperature) ensemble at various temperatures with corresponding system volume obtained from the NPT ensemble, we then computed the mean squared displacements (MSDs) of H and Pd atoms in the obtained structures for various times after an initial relaxation for 50 ps or so. The MSD sampling time ranges from 0.1 to 200 ns, depending on temperature, Pd phase, and whether Pd atoms were pinned. The varying sampling time is to ensure the motion of H in a given situation becomes stable. Each MSD data point in this work represents the average over about ten independent computations. More details of MSD sampling can be found elsewhere[19]. The MSD data were used to compute the diffusion coefficient based on the Einstein relationwhere t is the time and is the ensemble average of MSDs of diffusing atoms. For H diffusion, only MSD data higher than 10 Å were used to make sure the system is in the diffusive regime. Pd diffusion is much slower, and hence we only sampled the MSD at a fixed time length (100 ps) for implicating its mobility. The computations for carbon in crystalline iron are similar to the above. The atomic interactions of the Pd-H system were described by an embedded-atom-method-based potential[20]. Simulations[21] based on this potential has confirmed the concentration dependence of hydrogen diffusion in metals as observed experimentally. All the simulations in this work were performed using the code LAMMPS (ref.[22]) with Nos-Hoover thermostat and barostat, and the timestep was set as 1 fs.

We started our study with the thermodynamics of Pd matrix. For fcc Pd, its simulated melting point was found to be 1100 K. For the amorphous phase below the glass transition temperature, , Pd atoms essentially only vibrate about their local favorable positions within the time scale of interest to this study (Fig. 1a). Above , the MSD of Pd increases monotonically with temperature. During the simulations, crystallization of Pd at elevated temperatures was avoided by selecting a reasonably short time scale (but long enough for observing stable H diffusion). We note the existence of a critical temperature () in the supercooled region, as also reported for Pd in multiple-component amorphous alloys[23], which corresponds to the turning point of the Arrhenius curve and signifies the change in the diffusion mode from liquid-like motion to solid-like hopping upon cooling[23]. For the purpose of this work, we limit our discussion on H diffusion in the amorphous phase below .

Figure 1

Diffusion of H in Pd matrix, with Pd fixed and not fixed, respectively. (a) MSD of amorphous Pd atoms with respect to temperature at a given sampling time t = 100 ps. ( 660 K) is the estimated mode coupling temperature for Pd diffusion. (b) Hydrogen diffusion in amorphous Pd. The Arrhenius fitting gives = 1.52 and = 0.091 (in cm/s and eV, same for below), and = 1.43 and = 0.13 for Pd fixed. ( 400 K) is the estimated glass transition temperature, and ( 1100 K) is the melting point of fcc Pd. (c) Hydrogen diffusion in fcc Pd. The Arrhenius fitting gives = 3.76 and = 0.16 , and = 3.20 and = 0.20 if Pd is fixed. The experimental diffusion data are from reference[18]. (d) The ratio of H diffusion coefficient in fixed Pd () to that in mobile Pd (D). The errorbars are for standard deviation of 10 independent simulations (same for other figures unless otherwise specified). The data for H diffusion are in Supplementary Information.

We study the net effect of host phonons by comparing H diffusion with the vibrations of Pd atoms enabled and disabled, respectively. While impossible in experiments, this strategy is convenient in simulations and useful for tackling complex problems[24]. We first simulated the diffusion of H atoms with the Pd atoms fixed (by setting the force on them to zero) after the system was annealed for some time. Note that in this case, the kinetic energy of H is still connected to an external thermostat in simulations. From a theoretical perspective, in this case the collision between H and Pd is elastic since the collision does not transfer kinetic energy to Pd. As indicated by the open circles in Fig. 1b–c, H diffusion in this case fits well into the Arrhenius equation. As can be seen, like the random walk theory, our simulations based on the static Pd matrices predict a well-behaved Arrhenius diffusion of H. We note that in the random walk theory lattice relaxation during the jumping of interstitials is possible. Here the host matrix is relaxed for some random time before it is fixed.

Next we computed D of H without pinning the Pd atoms. As indicated by the solid circles in Fig. 1b–c, H diffusion in this scenario also follows the Arrhenius relationship. Figure 1c also shows the experimental[18] D of hydrogen from to K, which differs from our calculations within one order of magnitude. Other measurements[25,26] reported the values of D near 300 K to be cm/s, which are close to the extrapolation of the previous measurement[18] and differ from our calculations by more than one order of magnitude. Such a difference mainly results from the higher activation energy reported by experiments, which is not surprising in view of the trapping of hydrogen by lattice defects in experiments[27]. The accuracy of the EAM potential used in this work may also contribute to the difference, but overall the calculations resemble the trend of experiments.

Figure 1b–c clearly shows a general and significant difference in the diffusion coefficients for fixed and mobile Pd. We label the diffusion coefficients for the case of fixed Pd as since here only H jumping effect contributes to the diffusion. Figure 1d shows that the ratio decreases as temperature decreases and reaches below 30% at room temperature. To investigate the effect of Pd phonons on H diffusion, we split the total diffusion coefficient D aswhere by itself has the Arrhenius form , and represents the contribution of Pd phonons. The philosophy behind this treatment is as follows. The Pd phonon effect is related to the mass of Pd (). Imagine H diffuses in a hypothetical Pd isotope with infinite mass (). In this case, the vibration amplitude of Pd becomes zero and D reduces to . As w increases, the phonon effect emerges and contributes to H diffusion. This led us to formally write D as a Taylor serieswhere and the rest items can be wrapped as representing the phonon effect.

In the following we study from the perspective of Brownian motion. Brownian motion is generally regarded as a fluid phenomenon, but similar phenomena, such as thermal noise in electric conductors[28], exist in solids and all of them follow the fluctuation-dissipation theorem[29]. The splitting of D may raise concerns since the Brownian motion model is space/time continuous while the interstitials move in a space-discrete force field determined by the host structure. We note equation 3 implicitly treats the force field created by the host atoms aswhere corresponds to the host atoms at their pinned positions and the fluctuating augmentation arises from host phonons. The field results in a discrete diffusion component () and the Brownian contribution comes only from which is continuous in space/time. Such a treatment is indeed consistent with the solid-liquid phase transition: as temperature increases, the component increases and eventually the force field becomes liquid-like.

The diffusion coefficient of Brownian particles suspended in a fluid can be written as (see reference[30,31] or the Supplementary Information)where and a are the average kinetic energy and size of the Brownian particles, respectively, and is the liquid viscosity. A more familiar form for equation (6) is since equals based on the equipartition theorem. In view that H atoms are in thermal equilibrium with Pd phonons, we can directly apply Equation (6) to if we know the average energy of Pd phonons . In contrast to liquid of which is sensitive to temperature, here we can assume constant for Pd phonons. In other words, we assume a constant friction coefficient for the motion of H in the phonon fluid. At low temperatures[32,33], is approximately proportional to , where n is some constant. Substituting this into , we can writeAt high temperatures, the equipartition theorem gives . Note the actual (=) is smaller than by some roughly constant number because the heat capacity per phonon, C, at low temperatures is lower than its high temperature limit . In this case, we can write the phonon contribution aswhere c is a negative constant dependent on the material.

As shown in Fig. 2a, at high temperatures of both the amorphous and the fcc phases follows a linear behaviour with a negative intercept at . We note that for the amorphous phase the linearity holds for a range of temperatures above , which we attribute to the cancellation effect as at both and D deviate from the Arrhenius relationship for solids. At low temperatures, a function proportional to (n = 4 for amorphous and 5 for fcc) seems to fit well. The fitting parameters n = 4 and 5 are based on the Debye approximation[32] and the finite temperature field theory[33], respectively. However, as shown in the log-scale inset of Fig. 2a, a constant n actually does not fit well with the data, which contrasts with the established models. To build a theory for this observation is beyond the scope of this work, but we note the results are logically not surprising since, if n can increase from 1 at high temperatures to 4 or 5 at low temperatures, it is reasonable to assume higher n at even lower temperatures. Inspired by this, we found the data at low temperatures fit well with an exponential function, as shown in Fig. 2b. We note that can be approximated as D as because, with a higher activation energy, is an infinitesimal of higher order than D. This explains the exponential fitting of from the mathematical perspective.

Figure 2

Contribution of Pd phonons to H interstitial diffusion. The data points for the amorphous phase above are shown only for reference. (a) The linear parts can be fitted with = (4.06 ) cm/s and (3.62 ) cm/s for the fcc and amorphous phases, respectively. Power functions = 6.42 cm/s and 2.27 cm/s were used to fit the fcc phase, and = 4.74 cm/s was used for the amorphous phase. The plotted data for the fcc phase are shifted downwards for clarity. Inset: Zoom-in for low temperatures in log scales, same units as in the main. The linear lines are for eye guide only. The triangle indicates the upper limit of (by setting = 0) of fcc at 175 K (triangle slightly shifted to the left for clarity). (b) Arrhenius plot of .

For some comments to the above results, firstly we ignored the zero point vibrations in the above derivations, which only adds a small constant that is important as T approaches zero. Secondly, the quantum tunneling effect for H diffusion[10], which is important for temperatures below 200 K[34], was not included in this work. Thirdly, we expect the phonon contribution is a general effect and hence also affect other interstitials. Indeed, a recent molecular dynamics study showed that the permeability of He atoms and molecules in amorphous silica decreases if the thermal motion of silica is forbidden[35]. To extend our findings to larger atoms, we studied the movement of carbon in fcc iron at 1500 K and that of Li ions in at 1000 K, and the results (in Supplementary Information) confirmed the impacts of host phonons on interstitial motion. Finally, although atom-pinning decouples the motions of the host and interstitial atoms, it is a theoretical approach without an experimental counterpart. Hence, we also examined the phonon effects by only changing the mass of Pd atoms (without pinning Pd). This is equivalent to using Pd isotopes in experiments although in simulations one can arbitrarily set the mass of Pd. We used the fcc phase at 500 K as an example. Figure 3 shows that the diffusion coefficient of H decreases as the mass of Pd increases (or, equivalently, the phonon energy decreases), which supports our findings based on Pd-pinning.

Figure 3

Calculated hydrogen diffusion coefficients in fcc Pd lattice at 500 K as a function of the ratio of a hypothetical Pd mass () to the true Pd mass ().

In summary, based on molecular dynamics simulations we studied the effects of host phonons on interstitial diffusion, using hydrogen (H) within face centered cubic and amorphous palladium (Pd) as an example. Compared previous theoretical studies that coupled the phonons with the interstitial motions, this work decouples the two by pinning the host atoms and hence clearly reveals the net effects of host phonons. It was found that Pd phonons significantly promote H diffusion by causing Brownian-like diffusion of H atoms and dominate H diffusion below /2 ( being the melting point of Pd). Similar effects were also found for other important interstitials, such as lithium in manganese oxide and carbon in iron. This work establishes a new and improved physical picture for the general diffusion of interstitial atoms in solids and provides new perspectives for us to understand and design diffusion-related material behaviours. For example, the relaxation in metallic glasses[36] was found to correlate with the diffusion of the smallest constituting atoms within some loosely packed/bonded regions. Such diffusion was hypothesized to relate with the vibrations of atomic strings nearby[36], which is consistent with the current work. It is worth noting that the Brownian motion of nanosized liquid lead (Pb) inclusions within a solid aluminium (Al) was reported[37]. Although the authors found the motion is controlled by the shape of the inclusions and the diffusion of Al along the liquid/solid interface, the agitation force resulting from Al phonons seems to be an important reason for the observed motion.

Supplementary Information.