ARTICLE pubs.acs.org/JPCC
Vacancy Diffusion in NaAlH4 and Na3AlH6 Kyle Jay Michel and Vidvuds Ozolin-s* Department of Materials Science and Engineering, University of California, Los Angeles, P.O. Box 951595, Los Angeles, California 90095-1595, United States ABSTRACT: In recent years there has been considerable interest in sodium alanate as a prototypical complex hydride for solid-state hydrogen storage. Much effort has gone into understanding the rate-limiting processes in its hydrogen release and absorption reactions. The diffusion of metal species has been suggested as a possible kinetic bottleneck in Ti-doped materials. In this paper, we outline an approach for calculating the diffusivity of defects in complicated lattices using a combination of first-principles density-functional theory calculations and stochastic kinetic Monte Carlo methods. We apply this methodology to the diffusion of metal defects in NaAlH4 and Na3AlH6 that have been predicted to exist in large concentrations. We find that of the metal defects that exist in the largest concentrations, a neutral AlH3 vacancy is the most mobile in NaAlH4 (ΔHmig = 0.34 eV, D0 = 1.30 102 cm2/s, and DT=400K = 7.55 107 cm2/s) and that a negatively charged Na vacancy is the most mobile in Na3AlH6 (ΔHmig = 0.33 eV, D0 = 6.67 103 cm2/s, and DT=400K = 4.96 107 cm2/s). At T = 400 K, the calculated diffusion rates are an order of magnitude lower for charged AlH4 vacancies in NaAlH4 (ΔH = 0.44 eV, D0 = 1.41 102 cm2/s, and DT=400K = 3.92 108 cm2/s) and charged Na vacancies in NaAlH4 (ΔH = 0.43 eV, D0 = 2.96 103 cm2/s, and DT=400K = 1.19 108 cm2/s). This information is necessary for understanding the kinetics of mass transport during the hydrogen release and absorption reactions of NaAlH4.
’ INTRODUCTION Sodium alanate has been studied extensively since Bogdanovic and Schwickardi1 discovered that doping with Ti greatly increases its reversible hydrogen capacity. However, a clear picture of the atomic-level processes involved in this reaction has yet to be presented. Recently, studies have focused on determining if the creation and diffusion of native metal-containing defects through bulk NaAlH4 may be the rate-limiting step in the hydrogen release and absorption reactions of this material,2,3 which proceed according to NaAlH4 T
1 2 Na3 AlH6 þ Al þ H2 3 3
ð1Þ
Many processes must occur for reaction 1 to go to completion including nucleation, growth, and diffusion of hydrogen to the surface. Particular to the growth process is that there must be a long-range segregation of metal (Na and Al) when NaAlH4 decomposes into Na3AlH6 and Al. Many authors have suggested that a flux of native, metal-containing defects through bulk phases may provide the segregation mechanism and that this may be the rate-limiting step for the entire reaction.26 This is the third paper in a series dealing with the thermodynamic and kinetic properties of defects in NaAlH4 and Na3AlH6. In the first two papers, defect concentrations are calculated as functions of temperature with7 and without8 Ti substitution in the bulk. Having predicted which defects exist in the largest concentrations, and thus those that are likely involved in the fastest masstransport paths in reaction 1, the diffusivities of these defects are calculated in this paper. Together with the results of refs 7 and 8, these results can be used to quantitatively evaluate the fluxes of bulk defects and rates of mass and charge transport in reaction 1 r 2011 American Chemical Society
assuming that the rate of defect creation is fast enough to establish local equilibria at the relevant solidsolid interfaces (see ref7 for a discussion of our model). This work extends earlier studies of defect kinetics in NaAlH4 which either considered only the migration energies for simple charged defects (Na, AlH4, and H vacancies, and H interstitials),3 or obtained free energy barriers for complex neutral defects (AlH3 and NaH vacancies) at one temperature.2 Wilson-Short et al.3 also gave estimated lower bounds on the migration energies of complex defects such as AlH3 vacancies, which were approximated from combinations of simple defects. Neither of these studies attempted to calculate the rates or actual diffusivities as functions of temperature, nor was diffusion in Na3AlH6 considered, one of the key product phases of reaction 1 for which high concentrations of Na vacancies have been predicted.7,8 Therefore, open questions remain regarding the identification of the fastest mass transport mode of metal species that must accompany hydrogen release and absorption in NaAlH4. Here, we use ab initio molecular dynamics simulations to predict metal diffusion mechanisms via charged and complex neutral defects, followed by nudged elastic band9,10 refinement of the transition states. First-principles calculations of vibrational frequencies are used to obtain effective jump attempt frequencies, following which diffusivities are calculated using stochastic kinetic Monte Carlo methods. We apply these methods to the diffusion of Na, AlH3, and AlH4 vacancies in NaAlH4 and to the diffusion of Na vacancies in Na3AlH6; these have been found to Received: April 20, 2011 Revised: August 15, 2011 Published: August 23, 2011 21465
dx.doi.org/10.1021/jp203675e | J. Phys. Chem. C 2011, 115, 21465–21472
The Journal of Physical Chemistry C
ARTICLE
exist in the highest concentrations of all metal-containing defects7 and are therefore expected to dominate the kinetics of metal transport in reaction 1. Surprisingly, we find that the diffusivity of neutral AlH3 vacancies (which can be thought of as bound complexes of positively charged AlH4 vacancies and negatively charged H interstitials) is an order of magnitude higher than the diffusivities of charged Na and AlH4 vacancies in NaAlH4. This shows that in comparison with the AlH4 vacancy, the additional hydrogen in the AlH3 vacancy acts to reduce the migration barrier, contrary to what is commonly assumed in the literature.3 Our results provide a comprehensive picture of the temperature dependence of metal diffusion kinetics in a prototypical complexhydride-based hydrogen storage material.
’ METHODS Diffusivity. Diffusivities were obtained from kinetic Monte Carlo (KMC) simulations parametrized by attempt frequencies and migration energy barriers from first-principles DFT calculations. Given that all of the defects considered in ref 7 were predicted to exist in low concentrations (a few defects per 104 f.u. or lower at T = 400 K), interaction effects between them are negligible and KMC simulations can be carried out using single defects. Periodic boundary conditions were employed to remove edge effects. The residence time algorithm was used to extend the KMC simulations to time scales that are consistent with the calculated diffusion rates.11 At each Monte Carlo step, the total rate of a defect jumping from the current site i to any other site j was given by
Ritot
¼
Niend
∑ Ri, j
ð2Þ
j¼1
where Ri,j is the rate of jumps between each pair of sites and Nend i is the number of possible sites to which a defect at site i can jump. A random number, w1, was generated between 0 and Rtot i . The different ending states were placed in series with widths equal to their jump rates, and the system was advanced to the state j for which Ri,j1 e w1 < Ri,j. A second random number, w2, was generated between 0 and 1, and the simulation time was advanced by an amount equal to Δt = ln(w2)/Rtot i . We used the transition-state theory (TST)1214 to calculate the rate of a defect jumping from site i to site j mig
ΔHi, j Ri, j ¼ veff i, j e
=kB T
ð3Þ
where νeff i,j is the effective frequency at which a defect attempts to jump from i to j, ΔHmig i,j is the electronic energy of migration between these sites at absolute zero, kB is the Boltzmann constant, and T is the temperature. The effective jump frequency, as proposed TS GS where Q is the by Eyring,12 was described by νeff i,j = Qi,j /Qi partition function in the ground (GS) or transition state (TS). For a quantum mechanical system in the harmonic approximation this can be expressed as 2 !3 3N Q hνGS m 2 sinh 7 6 6 2kB T 7 kB T 6m ¼ 4 7 eff !7 6 ð4Þ νi, j ¼ TS 7 h 6 3N Q hνm 5 4 2 sinh 2kB T m¼5 th Here h is the Planck constant, νGS m is the m vibrational frequency of the structure with a defect in its ground state, and similarly
th νTS m is the m frequency of the structure with a defect in the transition state (frequencies were arranged in ascending order). The product of frequencies in the ground state was taken over all 3N 3 optical modes where N is the total number of atoms in the supercell. In the transition state this product also omitted the unstable optical mode at the saddle point since the translational partition function for this mode is included in the factor of kBT/h. DFT methods were used to calculate the enthalpies of migration and vibrational frequencies, as described in the next subsection. Due to the stochastic nature of diffusion processes, many simulations were performed at each temperature in order to obtain meaningful statistical results. One thousand KMC simulations were carried out for each defect at each temperature for a given length of time (Δt), after which the distances traveled in the x, y, and z directions (Δx, Δy, and Δz, respectively) were obtained. Total simulation times of Δt between 0.1 and 10 s were used depending on the calculated jump rates for each defect type. The diffusivity along a direction X (equal to x, y, or z) was obtained from
DX ¼
ÆðΔXÞ2 æ 2Δt
ð5Þ
where the average, Æ(ΔX)2æ, was taken over all simulations at a given temperature. The values obtained from eq 5 formed the diagonal elements of the diffusivity tensor. By symmetry, the offdiagonal elements of the diffusivity tensor are equal to zero for the tetragonal crystal structure of NaAlH4. In the monoclinic structure of Na3AlH6, deviation of the cell angles from orthorhombic symmetry are small (β = 89.491 from experimental results15), and it is reasonable to assume that the off-diagonal components can be neglected (i.e., we assumed that fluxes are only produced along the direction of a composition gradient and that the perpendicular elements of the flux can be neglected). From these results, the error in the diffusivity in the X direction, ΔDX, was obtained from sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Nsim 1 ΔDX ¼ ððΔXÞ2m ÆðΔXÞ2 æÞ2 ð6Þ 2Nsim Δt m ¼ 1
∑
where Nsim is the number of simulations at one temperature and (ΔX)2m is the square of the distance traveled in the mth simulation. The migration barrier for diffusion is often obtained in experimental studies by assuming that the diffusivity can be approximated by the Arrhenius equation so that mig
DðTÞ ¼ D0 eΔHfit
=kB T
ð7Þ
Here, D0 is the Arrhenius prefactor that has units cm2/s. In all figures showing the diffusivities of defects, a least-squares fit to the data was added and the migration energy for diffusion, ΔHmig fit , was extracted. This incorporated any temperature dependence that may be in νeff i,j as would be done in experiments. Ab Initio Calculations. Migration energies and attempt frequencies were studied using density-functional theory (DFT) as implemented in the Vienna ab initio simulation package16 (VASP). For all calculations we used ultrasoft pseudopotentials17 (USPP) with a cutoff energy of 250 eV and the generalized gradient approximation18 (GGA) for the electronic exchangecorrelation functional. The structures of NaAlH4 and Na3AlH6 were discussed in ref 7. Defects were introduced in supercells of NaAlH4 {Na3AlH6} where the unit cell lattice parameters and cell shape were fixed. Supercells of NaAlH4 {Na3AlH6} were 21466
dx.doi.org/10.1021/jp203675e |J. Phys. Chem. C 2011, 115, 21465–21472
The Journal of Physical Chemistry C used with dimensions 2a 2a c so that in the defect-free supercell there were 96 {80} atoms. Atomic positions15 were relaxed while the cell shape and volume were fixed to their values in the defect-free cell. Charge states for the various defects were achieved through the manual adjustment of the total number of electrons with the addition of a uniform compensating background charge. During relaxations, energies were found to converge to better than 0.1 meV/atom with a 2 2 2 Γ-centered k-point mesh when compared to a denser, 4 4 4 mesh. In order to predict the diffusion pathways of defects, we used ab initio molecular dynamics (AIMD) simulations starting from the fully relaxed structures. The temperature was set above 400 K (close to the temperature needed for significant hydrogen release according to reaction 11) and controlled using the NoseHoover thermostat.19,20 When diffusion events were observed during these simulations, snapshots of atomic configurations were taken before and after they occurred. Atomic positions from these snapshots were fully relaxed and served as the ground-state end points for subsequent DFT calculations of transition states and migration barriers. This allowed for the exploration of diffusion events that may have been difficult to intuit, especially for more complex defects such as AlH3 vacancies in NaAlH4. Diffusion pathways were found using the nudged elastic band (NEB) method9 where starting guesses for the intermediate structures were generated by interpolating images between the relaxed initial and final states obtained as just described. Once these paths had been determined, the transition state for each was further refined using the climbing image nudged elastic band (cNEB) method.10 The migration energy (ΔHmig i,j in eq 3) was calculated as the difference between the ground- and transition-state total energies. In order to include vibrational effects of solid phases, we used an implementation of the supercell force constant method21 to find the vibrational modes within the harmonic approximation. Here we constructed the dynamical matrix by calculating in VASP the forces that resulted from all symmetry-inequivalent atomic displacements in the 96- {80-}atom NaAlH4 {Na3AlH6} supercells. This matrix was diagonalized to obtain the characteristic frequencies (νm in eq 4) that were used to find the effective attempt frequency for the diffusion of each defect.
’ DIFFUSION IN NaAlH4 The crystal structure of NaAlH4 (see Figure 1) is comprised of 3+ ion negatively charged AlH 4 groups (nominally with one Al and four H ions) and positively charged Na+ ions. Here, the diffusion of various metalHx vacancies in NaAlH4 is discussed. We consider the native vacancies [AlH3], [AlH4]+, and [Na] in NaAlH4 where the naming convention is such that [AlH4]+ represents the removal of a negatively charged AlH 4 unit and [Na] represents the removal of a positively charged Na+ ion. These have been found to exist in the highest concentrations of all considered native metal defects7 and are therefore likely to lead to the highest mass transport rates in NaAlH4. For several reasons, we choose to omit the diffusivities of charged hydrogen defects in this paper even though some types have been predicted to exist in concentrations that are comparable to the metal defects.7 First, we are mostly concerned with metal-containing defects since there is experimental evidence that they are the rate-limiting species for mass transport.46 Second, in the section on defects in Na3AlH6 we present results strongly indicating that diffusion of hydrogen vacancies is faster
ARTICLE
Figure 1. Ninety-six-atom supercell of NaAlH4. Na ions are represented by yellow spheres, Al ions by blue spheres, and H ions by pink spheres. Tetrahedra are shown at AlH4 sites for visual clarity. Black lines show the allowed jumps between Na sites with equivalent vectors assumed for jumps between AlHx sites.
than diffusion of metal defects. Finally, diffusion of charged hydrogen vacancies, [H]+, and hydrogen interstitials, iH, in NaAlH4 has been studied by Wilson-Short et al.3 where it was shown that the migration barriers are 0.26 and 0.16 eV, respectively. These barriers are lower than the migration barriers for metal diffusion, and therefore hydrogen diffusion is unlikely to represent the rate-limiting processes in either Na3AlH6 or NaAlH4. Due to charge neutrality, hydrogen defects form in compensating pairs with metal defects, such as [Na][H]+ and [AlH4]+iH, and diffusional fluxes of metal species will not be limited by the diffusion of hydrogen vacancies or interstitials. We note here that previous studies3,22 proposed that charge transport via hydrogen defects is rate-limiting in undoped NaAlH4, and the catalytic effect of Ti doping was attributed to the lowering of the formation energy of [H]+ and iH due to Ti-induced pinning of the Fermi level. Our results presented in the second paper of this series8 show that the assumptions of Peles and Van de Walle22 about the structure of Ti dopants in NaAlH4 are not correct, and therefore the effect of Ti on the Fermi level is very different from that surmised in ref 22. During KMC simulations, only nearest-neighbor jumps in NaAlH4 are allowed and all symmetry-equivalent jumps between similar sites (i.e., Na to Na and AlHx to AlHx) are assigned the same rate. Jumps between second-nearest-neighbor sites were also considered but not included in the KMC simulations (in Figure 1 these correspond to jumps in the x and y directions). For [Na] diffusion, the rate of jumps to second-nearest-neighbor sites was calculated to be lower than that to the nearest-neighbor sites by a factor of roughly 104105 at T = 400 K. Therefore, the probability that these jumps will be sampled during the simulation is very low, and their omission will have a negligible effect on the results. During AIMD simulations, only jumps of AlH3 vacancies between nearest-neighbor AlHx sites were observed. 21467
dx.doi.org/10.1021/jp203675e |J. Phys. Chem. C 2011, 115, 21465–21472
The Journal of Physical Chemistry C
ARTICLE
Table 1. Migration Enthalpy (eV) and Jump Rate (jumps/s) for All Defects Considered in NaAlH4. Also Shown Are the Migration Enthalpy (eV) and Arrhenius Prefactor (cm2/s) from Least-Squares Fits to the Diffusivities in Figure 2 Averaged over the x, y, and z Directions defect [AlH3] +
[AlH4] [Na]
ΔHmig
R(T = 400 K)
ΔHmig fit
D0
0.329
8.127 108
0.336
1.301 102
0.418 0.377
4.265 10 1.195 107
0.441 0.430
1.410 102 2.963 103
7
Therefore, we do not have a prediction for the path that this defect would take to a second-nearest neighbor site. However, by analogy to the [Na] defect, the omission of this longer jump from KMC simulations likely does not have a significant effect on the results. Calculated migration enthalpies and jump rates at T = 400 K are given in Table 1 for defects in NaAlH4, and the calculated diffusivities are shown as functions of temperature in Figure 2 where linear least-squares fits to the diffusivities have been added. Also shown in Table 1 are the migration barriers and diffusivity prefactors from fits of the Arrhenius equation to the calculated diffusivities. The processes involved in the diffusion of these defects are discussed in the following subsections. [Na]. As shown in Table 1, the calculated migration barrier for [Na] vacancies between nearest-neighbor sites is 0.377 eV. The diffusivity of [Na] in NaAlH4 is shown as a function of temperature in Figure 2c. At all temperatures, the diffusivities along the x and y directions are lower than the diffusivity along the z direction roughly by a factor of 2. This is attributed to the local geometry at each Na site where the coordination number between nearest-neighbor Na sites is equal to four. Of the vectors connecting these four allowed jumps, only two have an x component or y component where all four have a z component. Since the rates for all jumps are equivalent, on average only onehalf of all jumps result in displacement along the x or y directions, while all jumps result in displacement along the z direction. This same argument applies to the AlHx vacancies that also have higher diffusivities along the z direction than along the x and y directions. The diffusion of [Na] in NaAlH4 is a relatively simple process and can be used to test the accuracy of the present KMC implementation. This defect is characterized by a single diffusion rate and jump distance, and negligible correlation effects. Thus, the diffusivity in one direction can be explicitly written as13 D ¼ nc f λ2 R
ð8Þ
Here, nc is the coordination number, f is the probability that a given jump is in the forward direction, λ is the distance traveled in a single jump, and R is the rate of jumps between sites. For the diffusion of [Na] along the x axis of NaAlH4 (as in Figure 1), this expression reduces to Dx = λ2xR where λx is the magnitude of the jump distance projected onto the x axis. This analytic expression is plotted along with the values obtained from KMC simulations in Figure 3 which shows that the two are in very good agreement. [AlH3]. The diffusion path of an AlH3 vacancy serves as an example of how AIMD simulations can be used to find processes that are not completely intuitive (this path is shown in Figure 4). Here, the diffusion event can be viewed in three stages. In the first stage, an AlH2 5 unit, which is created when the AlH3 vacancy is
Figure 2. Diffusivity of defects in NaAlH4 with temperature. Error bars are based on 1000 KMC simulations at each point. x, y, and z labelings are based on Figure 1.
introduced, moves toward the vacant site. In the second stage, a hydrogen ion is passed from the diffusing unit to a neighboring AlH 4 anion, and in the third and final stage, the diffusing unit, now AlH 4 , moves into the previously vacant site. At the transition 3 state two AlH 4 anions share an extra H ion, forming an Al2H9 complex where the migration barrier is 0.33 eV. This process results in the complete diffusion of [AlH3] as the structures in the initial and final states are equivalent by symmetry. A second path was also observed during AIMD simulations but was found to have a migration barrier that is higher in energy by 114 meV. This second path is similar to the low-energy one except that the diffusing AlH2 5 unit does not approach the AlH 4 unit as closely when transferring the H ion. As a result, the AlH bond must be extended further before the H is transferred which results in a higher energy barrier. Using the 21468
dx.doi.org/10.1021/jp203675e |J. Phys. Chem. C 2011, 115, 21465–21472
The Journal of Physical Chemistry C lowest-energy path as input for KMC simulations, the diffusivity for this defect is given as a function of temperature in Figure 2a.
Figure 3. Comparison of the explicit diffusivity from eq 8 and the values obtained from KMC simulations for the diffusion of [Na] in NaAlH4 along the x direction of Figure 1.
ARTICLE
By itself, the path shown in Figure 4 allows for net diffusion of an AlH3 vacancy along only the x or y axes of Figure 1 and does not allow for a defect moving along the x axis to make a jump along the y axis, or vice versa. However, we find that during AIMD simulations, the extra hydrogen left by the introduction of an AlH3 vacancy (which is bound in the AlH2 5 unit) is highly mobile. AIMD simulations of NaAlH4 with an AlH3 vacancy were carried out for a total of 250 ps during which three jumps of the vacancy were observed. However, during this same time more than 10 jumps were observed exchanging the extra hydrogen between the AlH 4 groups neighboring the vacancy. Due to hydrogen exchange, the vacancy and AlH2 5 are able to reorient and net diffusion is possible in all directions. [AlH4]+. During AIMD simulations, the lowest-energy path that was found for [AlH4]+ diffusion in NaAlH4 involves a rotation of the diffusing AlH4 unit in addition to the necessary translation. The migration barrier for this path is calculated to be 0.486 eV. A second path was also investigated in which the starting and ending configurations are manually adjusted in such a way that no rotation occurs during the translation of the AlH4 unit. The migration barrier for this path is found to be 0.418 eV, lower than the previous path by 68 meV. The harmonic
Figure 4. Diffusion of [AlH3] in NaAlH4. In all figures, hydrogen ions in the diffusing AlH unit are assigned distinct colors to emphasize the rotation and transfer of an H ion that occurs. Figures show the structure in the initial (a), transition (b), and final (c) states. 21469
dx.doi.org/10.1021/jp203675e |J. Phys. Chem. C 2011, 115, 21465–21472
The Journal of Physical Chemistry C
ARTICLE
Table 2. Migration Barriers (eV) and Jump Rates (jumps/s) of Na Vacancies in Na3AlH6a jump
ΔHmig
R(T = 400 K)
1
0.313
1.931 108
2 3
0.330 0.341
2.179 108 9.659 107
4
0.311
1.831 108
5
0.073
2.365 1011
6
0.747
4.275 102
7
0.100
1.172 1011
8
0.071
1.973 1011
9
0.635
1.322 104
10 11
0.620 0.090
8.288 103 2.634 1011
ΔHmig fit
D0
0.328
6.670 103
See Figure 5 for visualization of these jumps. The final row gives the migration barrier (eV) and Arrhenius prefactor (cm2/s) from a least-squares fit to the diffusivities in Figure 6 averaged over the x, y, and z directions. a
Figure 5. Na3AlH6 structure where yellow spheres represent Na ions, blue spheres Al ions, and pink spheres H ions. Octahedra are added at AlH6 sites to emphasize their orientations. Black lines and labeling show the allowed jumps during KMC simulations. Numbering of jumps corresponds to Table 2, and coloring of these numbers is simply added for visual clarity.
transition-state theory (TST) expressions in eqs 3 and 4 predict that the lower-migration-energy process also has a higher rate (assuming that the jump attempt rates are the same), which seemingly contradicts the AIMD simulation data. It remains an open question whether this is a failure of the harmonic TST and the process seen in AIMD simulations really corresponds to a lower free energy barrier or whether this is a statistical artifact due to short simulation time. Indeed, the [AlH4]+ diffusion events are very rare (we observed 1 event in 100 ps), and it is practically unfeasible to quantify the statistical significance of observing the higher-energy pathway instead of the lower-energy one. We have used the lower-energy jump in the KMC simulations as it represents a higher-rate process within the harmonic TST. Corrections beyond the harmonic TST are likely to increase the calculated jump rates, and therefore the calculated diffusivity values should be viewed as lower bounds. The diffusivity obtained from KMC simulations parametrized using this second path is shown in Figure 2b.
’ DIFFUSION IN Na3AlH6 Previous calculations involving native defects in Na3AlH6 have shown that negatively charged Na vacancies have the largest concentrations of any defect that contains a metal species.7 Therefore, we limit the evaluation of diffusivities in Na3AlH6 to this type. The Na3AlH6 structure is shown in Figure 5. There are two symmetry-inequivalent Na sites in this structure that
Figure 6. Diffusivity of [Na] in Na3AlH6 with temperature. Error bars are based on 1000 KMC simulations at each point. x, y, and z labelings are based on Figure 5.
are often referred to as large and small sites. The large site is in the octahedral void formed by AlH6 units and the small site in the tetrahedral void. From DFT calculations, the electronic energy to form a Na vacancy at the large site is lower than at the small site by 0.24 eV. In addition, there are also different orientations of AlH6 units, and combined, this results in many inequivalent jumps between Na sites. During KMC simulations, all jumps shorter than 4 Å between Na sites are allowed. Within this constraint we find 11 inequivalent jumps in the structure. Numbering in Figure 5 labels these different jumps, and Table 2 lists migration energies and diffusion rates for each one of these paths. Also shown in this table are the migration barrier and prefactor for diffusion from a least-squares fit of the Arrhenius equation to KMC data. Using the energies and frequencies for these paths as input for KMC simulations, the diffusivity of [Na] in Na3AlH6 is shown as a function of temperature in Figure 6. We find that positively charged hydrogen vacancies in Na3AlH6 charge-balance [Na] at all temperatures of interest for hydrogen storage.7 Due to the large number of possible jumps for 21470
dx.doi.org/10.1021/jp203675e |J. Phys. Chem. C 2011, 115, 21465–21472
The Journal of Physical Chemistry C hydrogen vacancies in Na3AlH6, the calculation of all migration barriers and frequencies needed for a KMC simulation is extremely cumbersome and has not been attempted here. Instead, we have calculated these values for only a small number of jumps of [H]+ between AlH6 units. Comparing these values to those obtained for [Na] jumps, we find that the magnitudes of migration barriers are comparable between these types. However, the effective frequencies of jump attempts (νeff i,j in eq 4) are much higher for [H]+ than for [Na] diffusion by factors of 510. Therefore, the diffusivity of [Na] is expected to be significantly lower that that of [H]+ in Na3AlH6. In order to preserve charge neutrality, diffusion of negatively charged [Na] will be accompanied by the diffusion of positively charged [H]+, and quantitative calculation of the limiting mass transport rates does not require the precise value of the higher-rate hydrogen vacancy process.
’ DISCUSSION AND CONCLUSIONS The migration barriers of defects in NaAlH4 obtained in this work can be compared with those obtained in a previous study by Wilson-Short et al.3 (current study vs ref 3): the barrier is 0.38 vs 0.41 eV for [Na], 0.42 vs 0.46 eV for [AlH4]+, and 0.33 vs 0.46 eV for [AlH3]. Small differences in the calculated values for [Na] and [AlH4]+ are attributed to different supercell sizes, numerical parameters, and/or pseudopotentials used in the DFT calculations. However, there are significant quantitative and qualitative differences in the migration mechanism predicted for [AlH3]. In the work by Wilson-Short et al., the diffusion of [AlH3] is regarded as a combination of the diffusion of [AlH4]+ and charged, interstitial hydrogen. The highest migration energy of these two defects is taken as a lower bound for the migration energy of the composite defect. However, the path pictured in Figure 4 with a migration barrier of 0.329 eV shows that this approximation is not valid as the diffusion of [AlH3] can follow more complicated paths than those that may be suspected by viewing the neutral vacancy as a superposition of charged [AlH4]+ and interstitial hydrogen. In comparison with the migration of an AlH 4 anion, we find that the extra hydrogen near an AlH3 vacancy facilitates the jump of the AlH2 5 unit due to the special structure of the transition state where this hydrogen ion is shared with a neighboring AlH 4 anion, cf. Figure 4. Therefore, we propose that although AIMD simulations are computationally expensive, they are an important tool for determining the often nonintuitive diffusion paths of complex defects. In Gunaydin et al.2 it was found that the free energy of migration for AlH3 vacancies in NaAlH4 is equal to 12 kJ/mol (0.125 eV/defect) at T = 400 K. We can use the harmonic frequencies in the ground and transition states already obtained to calculate the vibrational free energies as functions of temperature and evaluate entropic contributions to the migration barrier. We find that vibrational entropy lowers the free energy barrier at T = 400 K by a few tens of millielectronvolts from the 0.33 eV value obtained including only the static DFT energy. The former value is significantly higher than the results obtained from umbrella sampling.2 A possible explanation for this disparity is that umbrella-sampling methods use AIMD simulations that include anharmonic effects, while the harmonic TST only includes harmonic effects. Indeed the melting temperature of NaAlH4 is near 450 K,23 and anharmonic effects are expected to be significant at typical temperatures of hydrogen release (approximately 400 K). Therefore, the differences between the harmonic TST and AIMD umbrella-sampling results are not surprising.
ARTICLE
Comparison of the diffusivities of defects in NaAlH4 shows that the most mobile metal defect type is the neutral AlH3 vacancy. At T = 400 K, the average diffusivity of this defect is on the order of 106 cm2/s; the diffusivity of [AlH4]+ at this temperature is an order of magnitude lower, and lower by still another order of magnitude for [Na]. In Na3AlH6, the diffusivity of [Na] is slightly lower than 106 cm2/s. Therefore, in NaAlH4 metal ions are likely transported in the form of AlHx vacancies, while in Na3AlH6 metal ions are transported by Na vacancies. Typically, ball-milled grains of NaAlH4 are a few micrometers in diameter.23 In order for the two fastest diffusers, AlH3 vacancies in NaAlH4 and Na vacancies in Na3AlH6, to travel 10 μm at T = 400 K, it will take on average 1.3 and 2.0 s, respectively. This however relates only to a single defect where in actuality an enormous number of defects are needed to carry reaction 1 to completion. Therefore, concentrations of these defects at local equilibria must also be included when estimating the true time that this reaction will take to complete via the diffusion of defects. This further analysis is reserved for a forthcoming publication. Since the goal of this study is to identify the rate-limiting process in reaction 1 when carried out in bulk samples, defect diffusivities have only been studied in bulk phases. It is known that the activation energy for dehydrogenation of NaAlH4 decreases with particle size.24 Therefore, it is not appropriate to compare the results of this work to experiments carried out on nanoscale particles since the rate-limiting processes may be changed when the particle is made small and the characteristic diffusion times may fall below the average time required for defect creation.7 Furthermore, for sufficiently small particles, Na3AlH6 is not formed and NaAlH4 decomposes directly into NaH, Al, and 3/2 H2.25,26 In this regime, diffusion of defects in Na3AlH6 is not relevant. In this paper we have outlined an approach for the determination of diffusivities of defects in complicated structures using a combination of stochastic and first-principles calculations. Where not completely apparent, transition paths for vacancies have been investigated using AIMD simulations. Using this combination of methods, we have calculated the diffusivity of metal-containing defects in NaAlH4 and Na3AlH6. We find that of all the metal defects that exist in high concentrations in NaAlH4, AlH3 and positively charged AlH4 vacancies are the most mobile. In Na3AlH6, negatively charged Na vacancies exist in the highest concentration of all considered metal-containing defects and have diffusivity values that are comparable to those of AlH3 vacancies in NaAlH4. Therefore, the fastest metal transport in NaAlH4 occurs via diffusion of AlHx vacancies, while in Na3AlH6 metal mass flow is carried by Na vacancies. In conjunction with the results for defect concentrations in refs 7 and 8, these results suggest that diffusion of Na vacancies in Na3AlH6 provides the fastest mode of metal mass transport for reaction 1. With this and the two previous papers in the series, a comprehensive study of the thermodynamic and kinetic properties of native defects in NaAlH4 and Na3AlH6 has been presented. In particular, the concentrations of defects under local equilibrium conditions, as well as their bulk diffusivities, have been calculated as functions of temperature. In a forthcoming paper, these results will be combined within a mechanistic model for the diffusional mass transport in the (de)hydrogenation of NaAlH4, providing theoretical means to unambiguously identify the activation energies and rate-limiting diffusion processes from first-principles theory. 21471
dx.doi.org/10.1021/jp203675e |J. Phys. Chem. C 2011, 115, 21465–21472
The Journal of Physical Chemistry C
ARTICLE
’ AUTHOR INFORMATION Corresponding Author
*E-mail
[email protected].
’ ACKNOWLEDGMENT The authors gratefully acknowledge financial support from U.S. Department of Energy, Office of Science, Basic Energy Sciences, under grant No. DE-FG02-05ER46253. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. ’ REFERENCES (1) Bogdanovic, B.; Schwickardi, M. J. Alloys Compd. 1997, 253254, 1–9. (2) Gunaydin, H.; Houk, K. N.; Ozolin-s, V. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 3673–3677. (3) Wilson-Short, G. B.; Janotti, A.; Hoang, K.; Peles, A.; Van de Walle, C. G. Phys. Rev. B 2009, 80, 224102. (4) Lohstroh, W.; Fichtner, M. Phys. Rev. B 2007, 75, 184106. (5) Borgschulte, A.; Zuettel, A.; Hug, P.; Barkhordarian, G.; Eigen, N.; Dornheim, M.; Bormann, R.; Ramirez-Cuesta, A. J. Phys. Chem. Chem. Phys. 2008, 10, 4045–4055. (6) Ivancic, T. M.; Hwang, S.-J.; Bowman, R. C.; Birkmire, D. S.; Jensen, C. M.; Udovic, T. J.; Conradi, M. S. J. Phys. Chem. Lett. 2010, 1, 2412–2416. (7) Michel, K.; Ozolin-s, V. J. Phys. Chem. C 2011, DOI: 10.1021/ jp203673s. (8) Michel, K.; Ozolin-s, V. J. Phys. Chem. C 2011, DOI: 10.1021/ jp203672u. (9) Jonsson, H.; Mills, G.; Jacobsen, K. W. Classical and Quantum Dynamics in Condensed Phase Simulations; World Scientific: River Edge, NJ, 1998. (10) Henkelman, G.; Uberuaga, B. P.; Jonsson, H. J. Chem. Phys. 2000, 113, 9901–9904. (11) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. J. Comput. Phys. 1975, 17, 10–18. (12) Eyring, H. J. Chem. Phys. 1935, 3, 107–115. (13) Wert, C.; Zener, C. Phys. Rev. 1949, 76, 1169–1175. (14) Vineyard, G. H. J. Phys. Chem. Sol. 1957, 3, 121–127. (15) Ozolins, V.; Majzoub, E. H.; Udovic, T. J. J. Alloys Compd. 2004, 375, 1–10. (16) Kresse, G.; Furthmuller, J. Phys. Rev. B 1996, 54, 11169–11186. (17) Vanderbilt, D. Phys. Rev. B 1990, 41, 7892–7895. (18) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. Rev. B 1992, 46, 6671–6687. (19) Nose, S. J. Chem. Phys. 1984, 81, 511–519. (20) Hoover, W. G. Phys. Rev. A 1985, 31, 1695–1697. (21) Wolverton, C.; Ozolins, V.; Asta, M. Phys. Rev. B 2004, 69, 144109. (22) Peles, A.; Walle, C. G. V. D. Phys. Rev. B 2007, 76, 214101. (23) Bogdanovic, B.; Brand, R.; Marjanovic, A.; Schwickardi, M.; Tolle, J. J. Alloys Compd. 2000, 302, 36–58. (24) Balde, C. P.; Hereijgers, B. P. C.; Bitter, J. H.; Jong, K. P. D. J. Am. Chem. Soc. 2008, 130, 6761–6765. (25) Mueller, T.; Ceder, G. ACS Nano 2010, 4, 5647–5656. (26) Majzoub, E. H.; Zhou, F.; Ozolins, V. J. Phys. Chem. C 2011, 115, 2636–2643.
21472
dx.doi.org/10.1021/jp203675e |J. Phys. Chem. C 2011, 115, 21465–21472
