Quantum and isotope effects on hydrogen diffusion, trapping and escape in iron

We calculate the rate coef ﬁ cient as a function of temperature for lattice diffusion of hydrogen and its isotopes in a -iron; and also for trapping and escape from a vacancy. We employ Monte-Carlo and molecular dynamics methods based around the Feynman path integral formulation of the quantum partition function. We ﬁ nd large quantum effects including tunnelling at low temperature and recrossing at high temperature due to the ﬁ nite extent of the particle probability density. In particular these serve to in- crease the rate of trapping and to decrease the rate of escape at low temperature. Our results also show very clear non classical isotope effects. © Acta Materialia Inc. the CC


Introduction
One of the largest known diffusivities in the solid state is that of hydrogen in a-iron [1,2]. There are two reasons for this. One lies in the geometry of the body centered cubic lattice and its tetrahedral interstices; the other arises from the small mass of the proton leading to strong quantum effects, including large zero point energies and tunnelling [3]. It is typical that a rate coefficient may show Arrhenius behaviour at high temperature, T, and be essentially independent or weakly dependent on T otherwise [4] (see Fig. 4). A further complication arises in a-iron in that the transport of H is much attenuated by trapping of protons by crystal defects: dislocations, grain boundaries and vacancies, among others [5]. Hence the measurement of the lattice diffusivity presents many technical challenges [1] so that one would like to be able to separate out the effects of lattice diffusion and trapping by suitable theoretical calculations. We are concerned with the traps' capture probabilities [6]. In addition we are keenly interested in the mean residence time, t, of a proton trapped at a defect [6]; this is the inverse of the associated rate coefficient for jumping out of the trap. In addressing these matters we arrive at some rather startling conclusions concerning the roles of tunnelling through the barrier and recrossing at the saddle point in the potential energy surface.
We find that an interesting qualitative interpretation can be made from the behaviour of the "beads" in the path integral "necklace" in Feynman's picture.
The most severe approximation that we make is to assume that the hydrogen atom moves in a static lattice of iron atoms. This means that we cannot admit phonon assisted tunnelling [7]. However it allows us to work with a potential energy surface which provides a single degree of freedom in the classical transition state theory [4,8]. We do allow relaxation of the iron atoms, albeit in a rather stilted form: when the proton is in a reactant or product state (before or after a hop) it sees a lattice of iron atoms relaxed about the proton in its metastable position. We also hold the proton in a saddle point position and relax the iron atoms to provide an "activated complex" state. To locate saddle points we use a "nudged elastic band" (NEB) energy minimisation [9]. Interatomic forces that are required for these procedures are obtained from a magnetic tight binding (TB) model of H in iron [10]. This is not a severe approximation, since comparison with density functional theory (DFT) calculations shows good agreement in both concentrated and dilute limits [10]. In Fig. 1 we show a contour plot of the two potential energy surfaces.
Having established potential energy surfaces we calculate position probability densities (PPD) and quantum partition functions employing the Feynman path integral method [11] in a manner described earlier [12] using WangeLandau Monte Carlo [13]. To address trapping we consider the singly occupied vacancy as an archetypal trap for hydrogen. This is a much studied defect and regarded to have particular significance for the behaviour of H in airon [2,14]. The TB model has been shown to give a good account of the atomic structure and energetics of H binding to a vacancy compared to DFT calculations [10,15,16]. Fig. 2 shows a cartoon of the six possible hydrogen absorption sites. Fig. 3 shows the PPD of H trapped at an a-iron vacancy. It is very significant that at 50 K, although the centroid of the particle (in the language of path integral theory [17]) is constrained to remain at the dividing surface, the PPD clearly indicates that the proton has tunnelled through the barrier and largely escaped from the trap (indeed the proton has "split into two"). Note also, that at high T the proton is by no means localised and samples a considerable region of configuration space orthogonal to the reaction path, ie, in the region of the "dividing surface" having potential energies greater than at the saddle point. Therefore we would expect quantum effects effectively to lower the activation barrier at low T, but to raise it at high T.

Theory and results
The Feynman path integral method is a means to obtain quantum partition functions. We illustrate this for a single particle whose equation of motion is Schr€ odinger's equation in a potential energy V(x) (which becomes our potential energy surface in the configuration space of all the atoms in our simulations). The partition function of this particle is [11,19].
Here, E i are eigenvalues of the Schr€ odinger equation, r is the one particle density matrix, and the particle's probability density (see  (1) is taken are closed in the sense that the particle starts and finishes at the same point, x, in configuration space. The integral in (1) can be discretised for numerical purposes and this furnishes us with the approximate formula [19].

Zz mP
again over all closed paths having x iþP ¼ x i , ci [17]. But this is more than just a numerical convenience. It reveals that the partition function of a quantum particle is identical to that of a necklace of P classical particles, or beads, moving in a reduced potential energy V(x)/P and connected by springs of stiffness mP=b 2 Z 2 . Note that the stiffness is proportional to T 2 which means that at high T the necklace tightens towards a point particle (the classical limit) while  The six absorption sites of a hydrogen atom bound to a vacancy in bcc Fe. Up to six hydrogen atoms may be absorbed exothermically from bulk tetrahedral sites. The vacant site itself is not a trap site. Roughly, each proton is found near a tetrahedral site on the faces bounding the vacancy; however each is displaced slightly towards the vacancy and in some cases there are small lateral shifts. For precise locations, see Refs. [16,10]. at low T the necklace spreads out over the potential energy surface. The estimate converges to the exact result in the limit P / ∞. For our purposes we find P ¼ 20 is adequate. In what follows we will use the mapping of the quantum particle to the classical necklace to provide insight in the roles of temperature and isotope mass on trapping and escape of hydrogen at crystal defects in a-iron.
To describe transport and trapping we use reactive flux theory [20,4,21]. This is best introduced through classical transition state theory [22] (TST) beginning with the classical picture of a system that moves from reactant to product basins in a potential energy surface, F, passing through a saddle point. The rate coefficient in this reaction is [8], The term 1= ffiffiffiffiffiffiffiffiffiffiffiffi ffi 2pmb p may be interpreted as the "velocity" of the system as it approaches the saddle point. Since the iron lattice is fixed in our approximation we can take m to be the mass of the proton. The integral in the denominator is over all coordinate space A on the reactant side of the dividing surface and hence is essentially a reactant partition function Z r . Z * is the partition function for the system constrained to remain at a dividing surface S separating reactant and product basins and belonging to a configuration space of one fewer dimension. For comparison with experiment, even in the quantum case, one would like to cast this into Arrhenius form. To this end Vineyard imagines a third partition function Z Ã r (the "star" indicating reduced dimensionality) which belongs in a configuration space after the dividing surface has been translated from the saddle point down into the reactant basin [8]. Multiplying and dividing by this, the rate coefficient becomes   4. Path integral QTST rate coefficients k QTST (T) for lattice (bulk) diffusion and for trapping by, and escape from, a vacancy. Results are shown for the anti-muon, hydrogen, deuterium and tritium; and by scaling with ffiffiffiffi ffi m p the deviation from a classical isotope effect is revealed. The insets show the data at high T; for the bulk case we show the diffusivity, D ¼ a 2 k QTST /12 where a is the bcc lattice constant, and include a cross and circle at electrochemical permeation measurements of H from Refs. [1] and [18] respectively.
which is in the required Arrhenius form after recognising that DF is the reversible work, or free energy difference, in taking the constrained system from the reactant basin to the saddle point. This now passes rather neatly into a quantum transition state theory (QTST) by replacing classical partition functions with quantum mechanical ones [11,24]. In reactive flux theory the rate coefficient is [4,21], We can further calculate Z Ã r ðTÞ and so quite formally we arrive at The term k(t p ) is the transmission coefficient which we turn to presently. It is only necessary to find the two partition functions Z * and Z r to determine the QTST rate coefficient, k QTST , but if Z Ã r is also calculated then one can make contact with the Arrhenius form and enquire into how strong are the temperature dependences of the frequency prefactor and activation energies, which in the classical theory are T-independent. Fig. 4 shows QTST rate coefficients for lattice diffusion, and for trapping and release at a vacancy. In order to study isotope effects we show results for the anti-muon, proton, deuteron and tritium nucleus. In the larger plots we scale the rate coefficient with the square root of the mass: in a classical picture the data would thus fall into a single curve. Deviations from classical behaviour are very striking, except in the case of escape from the deep trap in which the barrier is too high to admit quantum effects. We should note that our theory is in reasonable agreement with the observed diffusivity of H in a-iron as indicated by data points in the inset of the left hand figure. The diffusivity shows a very large change in the Arrhenius slope below about 400 K so that an extrapolation from high temperature would lead to an error in the room temperature diffusivity of hydrogen of about a factor of two. The most striking quantum effect is seen in the case of trapping. The rate coefficient for jumping into a tetrahedral site neighbouring a vacancy from the nearest tetrahedral site in the next unit cell shows a very large non classical temperature dependence and isotope effect. This is evident in Fig. 3(a) which shows that at low T, with the centroid positioned at the dividing surface the proton has a small probability density at the reactant, trap site, although it has largely tunnelled out of the trap into the neighbouring tetrahedral site.
The transmission coefficient, k(t p ), enters as a factor in (2) consistent with the reactive flux theory of Miller et al. [20] who write the exact rate coefficient [21,25] k À t p Á k QTST ðTÞ ¼c fs À t p Á Z r ðTÞ in terms of a Kubo transformed fluxeside correlation function [23],c fs ðtÞ ¼ Tr e ÀbHF e iHt=Z qe ÀiHt=Z in which the Kubo transform is defined as [26].
Here H and F are the Hamiltonian and flux operators, the brackets indicating a commutator, and is a step function which is zero or one depending on whether the system is on the reactant or product sides of the dividing surface.
The transmission coefficient is calculated as the canonical ensemble average [21].
kðtÞ ¼ hdðq Ã À qÞðp=mÞqðqðtÞ À q Ã Þi hdðq Ã À qÞðp=mÞqðpÞi Here q is a reaction coordinate, the "star" denoting a constraint onto the dividing surface (as above) and a "bar" an average over the beads in the path integral necklace [17,21] At long times, t~t, k(t) tends exponentially to zero [4]; the hope is to identify a characteristic plateau time t p << t around which the transmission coefficient is fairly constant [4]. In classical TST k(t) is one: an assertion of the axiom that if the system reaches the saddle point with finite velocity in the direction of the reaction it will proceed into the product basin. The transmission coefficient thus accounts for recrossing of the dividing surface [4]. Ideally we'd like to find t À1 simply by multiplying the two terms in (2), however as we now show it's not always straightforward to identify the plateau. WangeLandau Monte Carlo has the benefit that it calculates directly the density of states so that thereafter the partition function can be obtained without further effort at any desired T [13]. On the other hand to find k(t) requires a lengthy set of computer simulations at each temperature of interest. We have calculated k(t) using ring polymer molecular dynamics (RPMD) [17,21], averaging over more than 10 6 trajectories. Fig. 5 shows our calculations at some representative temperatures for the cases of trapping and escape of H and D from the vacancy. In some cases it is not possible to identify a plateau region. We should point out that the RPMD is an approximation to the exact quantum dynamics only in the short time limit; as the simulation time increases, so does the possibility that the transmission coefficient calculated in RPMD deviates from the exact quantum transmission coefficient [25]. The product in (2) is independent of the choice of dividing surface [21]; however the two individual factors are not. Because of the difficulty in calculating k(t p ) unequivocally and as a smooth function of T, it is best to consider the two terms separately and examine their product qualitatively as we do here. Thereby we may observe clear and striking quantum and non Arrhenius properties of k QTST (T) as we do in Fig. 4 and match these with the rather startling conclusions from Fig. 5. We return to these shortly.

Discussion
In the following discussion, for brevity and clarity we will focus only on hydrogen and deuterium. For the case of bulk, lattice diffusivity we find that the transmission coefficient reaches a plateau value of one for both H and D at 200 K and 300 K; at 500 K k takes values of about 0.6 in D and 0.7 in H at 20 ps but has not yet reached a plateau. If we take it that at all temperatures k(t p ) / 1 in the case of bulk diffusion we can then use the QTST to study lattice diffusivity. Moreover the factorisation of the rate coefficient into Arrhenius form (3) is unambiguous since there is no arbitrariness in the choice of dividing surface in the case of a symmetric diffusion hop. To this end we show in Fig. 6 the activation energies and frequency prefactors as written in (3) as a formal rewriting of (2) with k(t p ) ¼ 1. We can draw a number of conclusions about lattice diffusivity from Figs. 4 and 6.
(i) The isotope effect is far from being described by the classical ffiffiffiffi ffi m p factor.
(ii) Both n and DF show a marked temperature dependence, at least at low temperature. (iii) Even at ambient temperatures and considerably above, quantum effects are non negligible. (iv) While the factorisation into frequency prefactor and activation terms is unambiguous in the case of lattice diffusion, the resulting terms are not necessarily those that are measured [1], so that whereas we make contact with experimental estimates of the activation energy, this is only an illustration. It is particularly striking that the actual proton frequency belonging to the potential energy surface is more than a factor of two smaller than the high temperature prefactor deduced after the factorisation.
We can now discuss Figs. 4 and 5 in relation to trapping and escape. The rate coefficient for trapping shows characteristic T-independence at low T to a greater extent than the bulk case because at low T we estimate that about 90% of the PPD resides in the reactant well, compared to exactly 50% in the bulk case. This serves effectively to reduce the low T activation barrier. Conversely the rate coefficient for escape is largely classical in its T-dependence. It is important to appreciate that the potential energy surface near the saddle point has negative curvature in the reaction coordinate direction, but otherwise is steeply rising. A classical particle may exist exactly at the saddle point but a quantum particle has a PPD that occupies a region of configuration space demanding that it acquires a larger potential energy, having amplitude to exist in the neighbourhood of the saddle point. This circumstance persists up to the highest temperatures that we have considered here; and it is useful to imagine that at low temperature the PPD may indeed be more closely confined near the saddle pointdthe PPD changes from "disk-shape" to "spherical" as the temperature is raised and this makes it harder to "squeeze" over the barrier. There is a wealth of curious phenomena contained in our calculated transmission coefficients. In the case of trapping, transmission coefficients for H are rather straightforward: at 200 K and 300 K, k(t p ) ¼ 1. This is due to the asymmetry of the PES which drops off steeply on the product side of the saddle pointdat low T the springs of the necklace are weak, the beads are well separated and once the leading bead falls into the trap it pulls the others after it. At 500 K the plateau is not yet reached after 20ps of simulation time but the transmission coefficient is expected to be less than one because the springs are become stiffer and recrossing is easier. The trapping of D shows similar behaviour. At low T the behaviour is as just described for H; again as T increases k falls, since tunnelling becomes less easy, and indeed reaches a minimum at around 1000 K after which we expect k to approach the classical value of one. In fact we have calculated k in this case for a single-bead, classical particle and it is indeed one. The escape of D is inhibited at and below 200 K by a vanishing transmission coefficientdthis is not the expected exponential decay as t / t, since in the QTST t z 200 ps (Fig. 4), this being a lower limit since k < 1. There is a deeper reason for the vanishing of  iron. The data drawn as error bars show estimates of the activation energy from H 2 gas equilibration (high T) and electrochemical permeation (low T) [1]. Estimates of the prefactor from the same data are 15e37 ps À1 and 11 ps À1 [1] which are not necessarily in reasonable accord with our theoretical calculations (see the text). The curvature of the tight binding, nudged elastic band, potential energy surface gives n¼13 ps À1 [15].
the transmission coefficient here. Escape is the reverse of the case for trapping: even if the leading bead reaches over the saddle point, it is dragged back into the trap by its weakly connected images which are still climbing the steep PES. At 500 K and above the springs are stiff enough to prevent this recrossing mechanism. Finally we turn to the escape of H from the trap. This is also "frozen" out at 200 K and 300 K as for D, but is reversed at very low T since at 100 K and 50 K k is small but non zerodan effect we attribute to tunnelling. The minimum in k for escape of H is at about room temperature, when tunnelling is weak, but the particle is not yet classical. We should repeat that in two situations, H escape and D trapping, Fig. 5 shows minima in the T-dependence of the transmission coefficient. For H escape k has a minimum of about 0.001 just below room temperature, and for D trapping k has minimum value of 0.2 at about 1000 K. In both these cases we attribute the increase in k as T / 0 to tunnelling and the increase at high T as due to the approach to classical behaviour. From that point of view neither H nor D are behaving classically at room temperature. Again, we point out that even up to 1500 K D trapping is not classical, since k ¼ 0.6 while a single bead simulation results in the classical result k ¼ 1.

Conclusions
In conclusion, we have calculated QTST rate coefficients and transmission coefficients for hydrogen and its isotopes undergoing lattice diffusion, trapping and escape from a vacancy in a-iron, using a tight binding formalism for interatomic forces and reactive flux theory using Wang-Landau Monte-Carlo and RPMD within the Feynman path integral formalism. At temperatures where tunnelling becomes competitive with classical over barrier reaction pathways (below 350 K) we observe opposing effects on escape and trapping rate coefficients. Transmission coefficients calculated by RPMD indicate that moderate tunnelling impedes escape from the vacancy while trapping into the defect is assisted. In the deep tunnelling regime (below 100 K) the trapping rate coefficient shows the largest deviation from Arrhenius form. The transmission coefficient for escape at temperatures around and below 100 K increases by two orders of magnitude which also is attributed to deep tunnelling. These discoveries will have a bearing on the interpretation of permeation and thermal desorption experiments.
At present it is necessary to examine k QTST and k independently, nonetheless the method is promising and shows a way forward to study defects of more immediate practical importance, viz. dislocations, and grain and interphase boundaries in steel.