A computational study of energy barriers of structural transformations and hydrogen transfer in boehmite

The crystal structure of boehmite (γ-AlOOH) contains a large amount of hydrogen bonds that are joined into chains by sharing hydrogen-bond donor and acceptor oxygen atoms. The hydrogen ions in the hydrogen-bond chains are highly mobile and have complicated structural characterizations, and this feature may well be utilized for proton-conducting applications, but the mechanism is unknown without the dynamic parameters of the hydrogen-transfer processes. We propose probable hydrogen-transfer paths and compute their energy barriers using density functional theory with van der Waals density functionals, on both perfect and vacancy-containing crystal structures. It is revealed that the energy barriers are generally below 21 kJ mol−1 in a perfect crystal, and 14 kJ mol−1 in a vacancy-containing structure. The low energy barriers are indicators of the high proton conductivity of boehmite even at room temperature.


Introduction
Boehmite (g-AlOOH) is an aluminum oxyhydroxide, metastable relative to diaspore (a-AlOOH) under ambient conditions. Boehmite occurs naturally as one of the main components of bauxite ore, 1,2 and has been studied by various research groups for multiple applications, including catalysts, 3 adsorbents, [4][5][6][7] and thin porous membranes. 8 For these applications, boehmite is frequently used as a precursor for manufacturing alumina phases, such as g-Al 2 O 3 and a-Al 2 O 3 , usually through thermal dehydration. 9 Less explored, however, is that boehmite itself may serve as a proton conductor for electrochemical reactions in aqueous environments.
Recently, a particularly large amount of work has been focused on proton-conducting materials, 10 which have many applications in the electrochemical industry. For instance, proton exchange membrane (PEM) fuel cells (PEMFC) rely critically on proton conductors, in form of thin membranes, for charge and mass transport. The general requirements of a good electrolyte membrane for PEMFCs are: high proton conductivity (>10 À2 S cm À1 ), good chemical stability in highly acidic media and thermal stability at low to medium temperatures, low cost, and ability to be processed into thin lms. 11 Many PEMs have been developed, some being commercialized, but no membranes meet all the criteria. Commercially available products like Naon developed by Dupont de Nemours also have some limitations, among which the most serious one is that proton conductivity drops signicantly if operating temperature is above 90 C. 12,13 Therefore, nding alternative materials to improve or replace currently available PEMs is an important aim of numerous researches about PEMs.
Previous studies have revealed that hydrogen bonding plays an important role in proton conduction mechanisms. 14,15 Thus, those structures that contain a large amount of hydrogen bonds are likely good proton conductors. Boehmite has a similar hydrogen-bond network as water, which was deemed to be an excellent proton conductor. 16 It is of technological interest to investigate the proton-conducting ability of the hydrogen-bond network in boehmite. Since boehmite dehydrate at about 400 C (variable with samples and moisture), 17 much higher than that of boiling temperature of bulk water under ambient pressure, it is possible to prepare proton exchange membranes that can work at higher temperatures than current Naon membranes. Elevated working temperatures in the context of PEMFCs mean improved reaction kinetics of the electrochemical reactions, and reduced chance of catalyst poisoning. Other advantages of boehmite include its low cost and well-developed synthesis methods. With all these advantages, the proton conductivity of boehmite, however, has seldom been measured, partly because boehmite is highly anisotropic and protons may only diffuse along a specic direction.
In boehmite's crystal structure, hydrogen-bonds connect to one another and form chains along crystallographic [001] direction of the orthorhombic cell. These parallel hydrogenbond chains bind interlocking double layers of edge-sharing AlO 6 octahedra that are stacked along [010] direction. The locations of Al and O atoms in boehmite have been well measured by X-ray diffraction (XRD) and neutron diffraction experiments. 18,19 The locations of H atoms in boehmite, however, have been elusive for a long time in spite of multiple structural characterizations based on Raman, 20,21 infrared (IR) 22 and nuclear magnetic resonance (NMR) spectroscopy. 23 Because H atoms are weak X-ray scatters, it is difficult to locate the H atoms in boehmite from XRD experiments. The crystal structure of boehmite was deduced to be in the space group Cmcm from the denite positions of Al and O atoms. The neutron diffraction technique, which has much higher resolution in measuring H atoms' positions than XRD, conrms the structural characterization results by XRD. The previous researches on synthetic lepidocrocite (g-FeOOH) 24,25 also attributed the structure to the Cmcm group.
In the unit cell of boehmite, two probable positions of hydrogen atoms are found to be compatible with the space group Cmcm in boehmite. The hydrogen atoms may reside in Wyckoff site 8e, forming O-H bonds that are perpendicular to the double-sheet planes. This is highly unlikely because it doesn't conform to the oxygen's orthogonal bonding structure, and researchers have not found signicant electron density at this position. 26 Ewing 25 proposed that the hydrogen atom is located in Wyckoff site 4a, the midpoint tween two interlayer oxygen atoms. By this way, symmetric O-H-O hydrogen bonds are formed. This structure is, however, not acceptable because of the exceptionally long bond length of the symmetric O-H-O bonds, approximately 2.67Å, which is apparently larger than 2.55Å for most of conrmed symmetric O-H-O bonds. 27 Several space groups conforming to asymmetric O-H/O hydrogen bonds have been considered, such as Cmc2 1 , 28 Pca2 1 and P2 1 /c. 29 The only real difference among these structures is the relative positions of the H atoms. Most of previous studies focus on relative energetic stabilities of the probable crystal structures of boehmite. Some questions, like whether the above crystal structures can transform into each other and how defects affect the structural transformations, are unsolved. Hence, we address the structural transformations for both perfect crystal structure and crystal structures with defects. In the rst part, we investigate the Cmc2 1 and Pca2 1 space groups, as proposed by previous studies, and propose a new space group of Pmc2 1 for boehmite. For the second part, we focus on the defects of hydrogen vacancies for simplicity. The structural transformations involve breakage and formation of O-H bonds in the hydrogen-bond chains, so we calculate energy barriers for these processes to assess the easiness (or difficulty) of the transformations.
These short-distance displacements of hydrogen atoms in boehmite, when combined with the long hydrogen-bond chains, may result in long-distance displacement of the hydrogen atoms, following the Grotthuss hopping mechanism. 30 By calculating the energy barriers for the structural transformations, we are able to estimate proton conductivity as the two properties are tightly correlated. Proton exchange membranes are primarily characterized by proton conductivity. 31 Before extensive experimental measurement of the proton conductivity, it is very useful of computational studies to assess probable mechanisms of hydrogen transfer in this highly anisotropic material.

Methodology
In this study we investigate thermodynamic stability of boehmite in space groups Cmc2 1 (space group no. 36), 28 Pmc2 1 (space group no. 26), and Pca2 1 (space group no. 29) 29 by comparing their total energies. Energy barriers associated with proton-transfer paths are calculated using the climbing-image nudged elastic band (CI-NEB) method implemented in VTST code. 32,33 All the calculations are done using density functional theory (DFT) implemented in Vienna Ab initio Simulation Package (VASP). [34][35][36] The electronic exchange and correlation interactions are described by the PBE (Perdew-Burke-Ernzerhof) functionals, 37 and van der Waals density functionals (vdW-DF) [38][39][40] combined with exchange part of PBE, referred to as vdW-PBE hereaer. The effects of nuclei and core electrons are described using the Projector Augmented Wave (PAW) method. 41,42 The PAW datasets for elements Al, O, and H are from the PAW dataset library shipped with VASP. The core radii are 1.90 Bohr for Al, 1.52 Bohr for O, and 1.10 Bohr for H.
As for computational settings, the cut-off energy of plane wave functions is 500 eV, electronic state occupations are modeled using the Gaussian smearing method with a small smearing width of 0.005 eV, and total energies are extrapolated to zero electronic temperature. Sampling in the rst Brillouin zone is done on special k-points generated according to the Monkhorst-Pack method 43 using a regular grid of 8 Â 8 Â 6 for the primitive cell of the Cmc2 1 structure, 10 Â 6 Â 8 for the Pmc2 1 structure, and 6 Â 6 Â 8 for the Pca2 1 structure, according to tests by a convergence criterion of 0.001 eV per formula unit. The numbers of irreducible k-points are 60, 60, and 36 for the Cmc2 1 , Pmc2 1 , and Pca2 1 structures, respectively. The convergence criteria are set to 1.0 Â 10 À7 eV for energy differences in self-consistent-eld iterations, and 0.005 eVÅ À1 for residual Feynman-Hellmann forces in relaxations.
Structures with hydrogen vacancies are generated by removing one hydrogen atom from 2 Â 1 Â 2 and 3 Â 1 Â 2 supercells with respect to the conventional cell of Cmc2 1 . We use two supercells (instead of one) to check the interactions between vacancies of image cells in the periodic boundary setting. The grids for generating k-points in the Monkhorst-Pack scheme are 8 Â 4 Â 8 for the 2 Â 1 Â 2 supercell, and 4 Â 2 Â 8 for the 3 Â 1 Â 2 supercell, resulting in 64 and 16 irreducible k-points, respectively.
For the transition-state calculations, we create 4 images between the starting and end structures by linear interpolations using a utility program of the VTST scripts. Our tests with 8 and 16 images for the transition from Cmc2 1 to Pmc2 1 show that the calculated energy barriers are already converged within 0.005 eV for the calculations with 4 images. The force-based quick-min algorithm is used for nding minimum-energy path. Cell parameters are kept constant for all the transition states.

Crystal structures without defects
We can clearly see that boehmite has a layered structure (see Fig. 1) that contains AlO 6 octahedron layers and hydrogen-bond layers. Each hydrogen-bond layer is composed of parallel hydrogen-bond chains along the [001] crystallographic direction. The boehmite crystal structures of Cmc2 1 , Pmc2 1 and Pca2 1 space groups are all orthorhombic, and the main difference is orientations of the hydroxyl groups. Since around each oxygen atom in the hydrogen-bonding layers there are two positions at which a hydrogen atom may reside, corresponding to two possible congurations for a hydrogen-bond chain. For the Cmc2 1 structure, all the hydroxyl groups (O-H) tilt towards the [00 1] direction; for the Pmc2 1 structure, the hydroxyl groups switch directions in adjacent hydrogen-bond layers; and for the Pca2 1 structure, the hydroxyl groups in adjacent hydrogen-bond chains in a same layer have opposite directions.
The calculated lattice parameters of the Cmc2 1 structure are listed in Table 1, together with experimental values. By comparing the lattice parameters from plain PBE and vdW-PBE calculations, we nd that the values for a are almost the same (less than 0.01Å or 0.3% in differences) in both settings, the values for c are also almost identical, and both are larger than the experimental values by about 1.1%. As for the lattice parameter b that is pertinent to interlayer distances, the value calculated by vdW-PBE is slightly smaller than the experimental value by about 0.01Å, while the value by PBE differ from experimental value by about 0.1Å. In general, vdW-PBE performs better in reproducing the lattice parameters of boehmite.
It is noted that the orthorhombic cell of Pca2 1 is twice that of Cmc2 1 and Pmc2 1 . The total energies of the three structures, i.e., Cmc2 1 , Pmc2 1 , and Pca2 1 , are compared to sort their relative thermodynamic stability. From Table 2 we can see that the Pmc2 1 structure has almost identical total energy as Cmc2 1 in both PBE and vdW-PBE calculations, but Pca2 1 has consistently lower total energy than Cmc2 1 and Pmc2 1 , by À0.27 kJ mol À1 in PBE or À0.26 kJ mol À1 in vdW-PBE calculations. Such small energy differences are close to resolution limits of experimental measurements and computations. The difference should be treated carefully because even though it is not negligible, the value is one order of magnitude smaller than the thermal uctuation energy at room temperature (about 2.5 kJ mol À1 ). Thus, we show that the thermodynamic stability of Cmc2 1 and Pmc2 1 are essentially the same, and Pca2 1 is slightly more stable than Cmc2 1 and Pmc2 1 .
The order of thermodynamic stability of the three structures is related to the orientation of hydroxyl groups in hydrogenbond layers. In Cmc2 1 and Pmc2 1 , all hydroxyl groups in a hydrogen-bond layer point to a same direction (either [0 0 1] or [0 0 1]), and have the same polarity. In Pca2 1 , hydroxyl groups in a hydrogen-bond layer point to opposite directions in alternating chains, and collectively have zero polarity. Interactions between electric dipoles of a same polarity increases the system's energy, while interactions between electric dipoles of opposite polarities are associated with a negative energy. It appears that dipolar interactions are negligible across neighboring hydrogen-bond layers, since the energies of Cmc2 1 and Pmc2 1 are almost identical. The interactions of dipole moments in adjacent chains in a same layer should be very weak, at the order of 0.1 kJ mol À1 that is readily overridden by thermal uctuations at room temperature. We conclude that the three structures have the same thermodynamic stability within an error of 0.5 kJ mol À1 , and tentatively propose that hydroxyl polarity contributes slightly to their relative thermodynamic stability. In a same layer, the orientation of the hydroxyl groups in different zigzag-shaped chains should be independent of each other at normal conditions. This conclusion is consistent with a previous report in the literature. 26 Because the three structures have almost the same thermodynamic stability, a boehmite sample may have domains of the three probable structures. The next question is whether they can transform into one another. If the structural transformations are easy at certain conditions, then the domains should be dynamic or volatile, otherwise the domains are static. We can assess energy barriers associated with the structural transformations, and compare the energy barriers with possible driving forces. Since the dipolar interactions between hydrogenbond chains are weak, it is convenient for us to assume that only one hydrogen-bond chain changes and all other chains are constant. We study the transformation between the Cmc2 1 and Pmc2 1 structures, and propose two paths for hydrogen atoms to cause this structural transformation: the stretch mode and the   swing mode (see Fig. 2). In the stretch mode, all the hydroxyl bonds break and hydrogen atoms move to nearby oxygen atoms to form new hydroxyl bonds (Fig. 2a). In the swing mode, a hydrogen atom remains around a same oxygen atom, and the hydroxyl bond does not break but only changes its direction (Fig. 2b). Each mode requires cooperative motion of all hydrogen atoms in a whole hydrogen-bond chain, otherwise a high-energy state may occur that either has two hydrogen atoms bonded to one oxygen atom, or has a very small distance between two hydrogen atoms boned to two neighboring oxygen atoms. We point out that both modes can independently transform the structure between Cmc2 1 and Pmc2 1 , but longrange transfer of hydrogen atoms is only possible with both modes cooperating in concert. We use the CI-NEB method to calculate the energy barrier for the two hydrogen-transfer modes. As shown in Fig. 3, the PBE results indicate that the stretch mode has lower energy barrier (13.52 kJ mol À1 ) than the swing mode (20.68 kJ mol À1 ), while the vdW-PBE results show that the stretch mode and swing mode have almost the same energy barrier (18.70 and 18.26 kJ mol À1 , respectively). The energy-barrier difference between these two modes is 7.16 kJ mol À1 for PBE and 0.44 kJ mol À1 for vdW-PBE, showing the inuence the van der Waals forces on the transition states that involve hydrogen bonds in boehmite. Especially for the stretch mode, the energy barrier estimated by vdW-PBE is 5.18 kJ mol À1 higher than that by PBE. The large difference in energy barriers by the two methods can be understood by the arrangement of the atoms along the transition paths. In the stretch mode, the hydrogen atoms that are attached to one AlO 6 -octahedron layer move to another AlO 6 layer, and this geometric change for the intermediate states is accompanied by energy increase by van der Waals interactions. In the swing mode, the hydrogen atoms are invariantly attached to one AlO 6 -octahedron layer along the whole transition path, for which the geometry change is smaller than that in the stretch mode.
Since the energy barriers are around 20 kJ mol À1 , the energy required for transformation between the crystal structures can be provided by thermal uctuations. Energy uctuations at temperature T is roughly estimated to be RT (where R is gas constant, 8.31 J K À1 mol À1 ). Even though it is thermodynamically metastable relative to diaspore (a-AlOOH) under ambient conditions, boehmite can be stable at temperatures below its dehydration temperature, about 700-750 K, before transforming into g-Al 2 O 3 . 45-47 At T ¼ 300 K, RT z 2.5 kJ mol À1 , and at the dehydration temperature (around 750 K), RT z 6.2 kJ mol À1 . In this temperature range (300-750 K), the energy barrier of the transition between Cmc2 1 and Pmc2 1 space groups is 3-8 times the RT value. It turns out that this structural transition is highly probable even at room temperature.
We can use the Arrhenius equation to calculate structural transition rate as 48 where n stands for the frequency of successful structural transformation, A is a prefactor that is assumed to be a systemdependent constant, DE is the energy barrier for the transition, R is the gas constant, T is the temperature, and n 0 is the attempt rate, which corresponds to the frequency of a characteristic phonon mode. The proton-transfer modes involve coordinative motion of atoms, and certain phonon vibration modes may assist the processes if the motion of atoms in the phonon modes coincides with that in the proton-transfer modes. Tunega et al.
have obtained boehmite crystal's phonon spectra, 48 where two OH vibration modes match the stretch and swing modes of proton transfer in the present study. The value of n 0 for the stretch mode (2947 cm À1 ) is about three times that of the swing mode (1283 cm À1 ). Assuming that the prefactor A is the same for the two-transition modes, we can estimate the relative rate of the two transition modes. At room temperature (approximately 300 K), for the energy barriers calculated using PBE, n stretch n swing z 52:6  and for the vdW-PBE case: n stretch n swing z 3:6 In both PBE and vdW-PBE calculations, the rate of stretchmode transitions is higher than swing-mode transitions. But, the rate differences and the causes are quite different in PBE and vdW-PBE calculations. In PBE, the rate of stretch-mode transitions is several tens of times of swing-mode transitions, and the difference is mainly caused by its smaller energy barrier; while in vdW-PBE, the rate of stretch-mode transitions is only about four times that of swing-mode transitions, and the difference is mainly caused by the characteristic phonon frequency.

Crystal structure with defects of hydrogen vacancy
Real boehmite crystals always contain various forms of defects, among which we concentrate on hydrogen vacancies. The normally continuous hydrogen-bond chains in a perfect lattice break at hydrogen vacancies, splitting into shorter chains. The hydrogen vacancies are expected to interrupt transition of nearby hydrogen atoms. In this study, we create hydrogen vacancies by removing a hydrogen atom from a supercell of boehmite, then design paths for hydrogen transfer to or from this vacancy position, and calculate energy barriers associated with these transition paths. Along the transition paths there are intermediate states for which the hydrogen atom in motion resides at one of the two equilibrium positions around an oxygen atom. In Fig. 4 we present ve intermediate states, namely VacH, VacH_s1, VacH_s1t2, VacH_t1 and VacH_t1s2. Two paths can be composed by connecting subsets of these states for transferring a hydrogen atom along the [001] or [00 1] direction. The rst path is from VacH to VacH_s1 via a swing-mode motion of a nearby hydrogen atom, and to VacH_s1t2 via another stretch-mode motion of the hydrogen atom. By repeating this path, the hydrogen-vacancy moves continuously towards the [00 1] direction and, equivalently, the hydrogen atoms in this hydrogenbond chain move toward the [001] direction. The second path is from VacH to VacH_t1 via a stretch-mode motion of a nearby hydrogen atom, and to VacH_t1s2 via another swing-mode motion. The motion of the hydrogen vacancy and hydrogen atoms are both opposite to those in the rst path.
These four steps involved in the two paths have different energy-barriers. We have used the CI-NEB method to calculate these energy-barriers and present the results in Table 3 and Fig.  5. We can clearly see that the two paths have very similar highest energy barriers in both PBE and vdW-PBE calculations. From our PBE calculations, the swing-mode steps, i.e. from VacH to VacH_s1 and from VacH_t1 to VacH_t1s2, have higher energy barriers than the stretch-mode steps, i.e. from VacH to VacH_t1 and from VacH_s1 to VacH_s1t2. The differences in these energy barriers are 1-2 kJ mol À1 . For vdW-BPE calculations, the stretch-mode steps have higher energy-barriers than the swingmode steps by approximately 6-8 kJ mol À1 . The reasoning on the effects of vdW interactions on the stretch-and swing-modes has been discussed in the previous section for the perfect crystal, and also applies for the structures with hydrogen vacancies. The effects of hydrogen vacancies on the energybarriers will be discussed in the next paragraph.
By comparing the energy-barriers in Fig. 3 and Table 3, one nds that the hydrogen vacancy lowers energy-barriers of both  stretch-and swing-modes, but to different extents. The energybarrier of the stretch mode is lowered by 4-6 kJ mol À1 , and that of the swing mode is lowered by 10-12 kJ mol À1 , both consistently in PBE and vdW-PBE calculations. The lowering in the energy barriers can be explained as follows. In a perfect crystal, all hydrogen atoms in a hydrogen-bond chain can only move coordinately to avoid the otherwise high-energy states. In the proximity of a hydrogen-vacancy, a hydrogen atom can move independently of other hydrogen atoms in either stretch or swing mode, requiring less energy than the motion of many hydrogen atoms in a perfect crystal. As for the difference in the vacancy-induced decreases of energy-barriers between the two modes (stretch and swing), it can tentatively be explained from the change in electrostatic interactions between the hydrogen ions. We assume the hydrogen ions carry the same amount of positive charge in a perfect crystal and vacancy-containing structures, thus the electrostatic interactions depend solely on distances between the ions. For simplicity, we only consider the electrostatic interaction between a hydrogen ion and its nearest  neighbor. Along each transition path the state of the highest energy, or the saddle point, can be identied, and some representative ones are shown in Fig. 6. For the saddle-point structures in the two stretch-mode transition paths, the distance between the nearest neighboring hydrogen ions changes from 1.866Å in a perfect crystal (Fig. 6a) to 2.125Å in the vacancycontaining structure (Fig. 6b). For the two swing-mode transition paths, the distance between the nearest neighboring hydrogen ions changes from 1.874Å in a perfect crystal (Fig. 6c) to 2.505Å in the vacancy-containing structure (Fig. 6d). Starting from similar distances between nearest neighboring hydrogen ions, the larger increase in the H-H distances in the swingmode steps leads to larger decrease in the associated energy barriers.
We have tested two sizes of supercells, namely 2 Â 1 Â 2 and 3 Â 1 Â 2, in calculating the energy barriers in the vacancycontaining structures that, similar to many other defectcontaining structures, are expected to be sensitive to the supercells' sizes in the periodic boundary setting. For both PBE and vdW-PBE calculations, the differences between the energybarriers in the two supercells are below 1 kJ mol À1 , as can be derived from the numbers in Table 3. The trends of the energy barriers with respect to the steps are also consistent (see Fig. 6). We are condent that the calculated energy barriers are converged with respect to the supercell's size at 2 Â 1 Â 2.

Conclusions and remarks
In this paper, our computational results are mainly about thermodynamic stability of three probable structures of boehmite, and energy barriers for hydrogen transfer in boehmite's lattice with and without hydrogen vacancies. The results show that the Pca2 1 structure has slightly lower total energy than the Cmc2 1 and Pmc2 1 structures due to different polarities of the OH bonds. Transfer paths of hydrogen atoms along the hydrogen-bond chains can be divided into steps of swing-and stretch-mode transitions, and are associated with different energy barriers that depend on computational settings. For a perfect crystal the energy barriers are below approximately 21 kJ mol À1 , and for the structures with hydrogen vacancies, the energy barriers are below approximately 14 kJ mol À1 . The computational settings of PBE and vdW-PBE lead to different orders of the energy barriers for swing-mode transitions and for stretch-mode transitions, in both perfect crystals and vacancy-containing structures. It has been observed van der Waals forces are important for hydrogen-bonds, 49 thus the results of vdW-PBE should be more reliable and closer to reality.
The energy barriers are generally low in comparison with thermal uctuation energies at room temperature and above, suggesting high mobility of hydrogen atoms in boehmite's lattice. Thus, the crystal structure of boehmite may comprise of all the three structures that are frequently changing into one another. It is difficult to identify domains of different polarity of the OH bonds without lowering the temperature. The average structure belongs to the Cmcm space group that is repeatedly conrmed by multiple experiments.
We suggest, based on the low energy barriers, that boehmite has a high proton conductivity that may be utilized for protonconducting applications. The proton transfer mechanism in boehmite do not rely on liquid-phase water, and the working temperature can be elevated to above 100 C. It has been stated in ref. 50 and 51 that a working temperature of 120 C and 50% or lower relative humidity are targets for the development for automotive use of hydrogen/air fuel cells. The only drawback is that the proton conductivity of boehmite is highly anisotropic. Synthesizing and assembling highly ordered boehmite crystals are challenging, but it should be much easier to make boehmite ller materials for organic membranes to enhance their hightemperature proton conductivity.

Conflicts of interest
There are no conicts to declare.