Theoretical study on hydrogen storage capacity of expanded h-BN systems

In this work, the hydrogen storage capacity of the expanded hexagonal Boron Nitride (eh-BN) systems has been presented. We have employed a new equation of state (EOS) for hydrogen gas to ﬁgure out the hydrogen density distribution proﬁles in the eh-BN systems. In this regard, the environmental conditions (i.e., temperature and pressure) are considered in the prediction procedure using DFT single point calculations. The eh-BN systems with different layer spacings are studied by PBE method with consideration of the long range dispersion corrections. On account of the in-plane polar bonds, a series of adsorption positions are considered. Additionally, the adsorption energy and hydrogen density proﬁles are reported. Furthermore, the relation between uptakes and the interlayer spacings with the effects of the environmental conditions are also provided. The limit of the physical hydrogen storage capacity in a perfect eh-BN system at 243 K and 10 MPa is founded to be 2.96 wt.%.


Introduction
As one of the most promising clean energy source, the compact storage and safe delivery are the main challenges for the application of the hydrogen energy.The short-term goal of the designed novel hydrogen storage materials is to reach the gravimetric density limit of hydrogen, i.e., 5.5 wt.% and 0.04 g/cm 3 at ambient temperature and pressure by the year 2017, which is targeted by the U. S. Department of Energy.In the automobile industry, to fuel a typical vehicle, the hydrogen storage system is required to contain at least 7.5 wt.% and 0.07 g/cm 3 hydrogen fluid at the acceptable limit T ¼ 243 K and P ¼ 10 MPa [1].Another notable criterion includes the system recyclability.The hydrogen molecule should also be able to be desorbed from the storage system under ambient conditions.Hence, the adsorption energy of a single hydrogen molecule in an effective storage system should lie between the physisorbed and chemisorbed states.More precisely, this adsorption energy needs to be around 0.1-0.3eV per H 2 molecule [2][3][4].
Initially, the perfect graphene was considered as the most promising hydrogen storage material owing to its large surface area and light weight.It has already been confirmed that the perfect monolayer graphene do not have strong hydrogen storage capability ascribe to its surface inertia.However, graphene is still a rational model to test and verify the theoretical methods for gas storage prediction due to its intrinsically simple structure and high symmetry.The monolayer h-BN has nearly identical geometrical configuration as graphene.On the other hand, the B-N bond is a sort of in-plane polar bond.In this situation, the distribution of the hydrogen molecules on the h-BN surface is more complicated than that on graphene.Therefore, a purposeful test for theoretical method with monolayer h-BN model is worthwhile.Moreover, the polar B-N bond was expected to bring greater hydrogen storage capacity than the C-C band in graphene.For this reason, many nanostructures based on the h-BN have also been studied as molecule adsorption materials [5][6][7].Zhao et al. [5] claimed that a 7.69 wt.% hydrogen storage capacity was found in the graphene/boron nitride heterobilayer system through the theoretical studies at LDA level.Another even higher theoretical uptake value, around 8.7 wt.%, was found in the Li-decorated hybrid boron nitride and graphene domains by Hu et al. [6], also at LDA level.Oxygen-doped BN nanosheets were studied by Lei and co-workers [8], experimentally and theoretically.The system had 5.7 wt.% capacity for hydrogen storage under 5 MPa at room temperature.Recently, Kumar et al. [9] presented the importance of the bond exchange spillover mechanism in the h-BN nanostructures as hydrogen storage materials.A 7.4 wt.% uptake was predicted in the h-BN nanostructures at GGA level with dispersion correction in their work.At once, Zhang et al. [10] studied the Li-doped porous BN material for hydrogen generation and storage.The authors reported a 7.5 wt.% theoretical uptake derived from the around 0.16 eV adsorption energy in such systems.However, the theoretical predictions via DFT calculations always seem to be overly optimistic because of the missing thermodynamic information.
The most theoretical studies on the hydrogen material design are based either on the DFT calculations at ground state, or on the statistical mechanic methods, e.g. the molecular dynamics (MD) or Monte-Carlo (MC) simulations [11].The advanced QM/ MD-lVT simulation [12] is until now still expensive and timeconsuming.Morris' group, following the works of Patchkovskii [13] and Cabria [14], combined the DFT single point calculations and the thermomechanical boundary conditions by using an effective continuum model together to predict the hydrogen and methane uptakes in nanoporous carbons [15][16][17][18].This method can provide the hydrogen storage capacity with full thermodynamic information only via a series of DFT single point calculations.The time for the thermodynamic simulations can be saved in this way.A new EOS for hydrogen fluid is introduced in our previous work [19] to improve the prediction accuracy of this combination method.The determination of the fugacity coefficients of the hydrogen fluid is the key point in a similar prediction research.We believe that with the help of our EOS, the fugacity coefficients of the hydrogen fluid can be estimated with full physical meanings, and show very good agreements with the experimental results at reasonable temperatures and pressures.
The eh-BN is a simple and effective theoretical model of the nanoporous h-BN material.In the present work, the relation between the storage capacity and the interlayer spacing (i.e., the pore size) is discussed to find the storage limit of the eh-BN systems at certain temperature and pressure.This work is organized as follows: In Section 2, the models of the hydrogen molecule and eh-BN systems will be introduced.The simulation method and parameter arrangement will be presented there, too; In Section 3, the main simulation results will be presented and discussed.The adsorption energy map and the pressure profiles inside the eh-BN systems will be shown.The relation of the hydrogen uptakes to the interlayer spacings of the eh-BN will be discussed; A short conclusion will be given in Section 4.

Methodology
As a premise, the gas storage system in this research is assumed to be in equilibrium.As known, one of the most important equilibrium condition is that the chemical potential should be a continual function at the spatial boundary, i.e., the chemical potentials in the external free gaseous phase and in the adsorbed phase should have the identical values on the boundary line.Therefore, the following equation should be ensured in the adsorbed phase and the free gaseous phase: where the subscripts in and ext denote the adsorbed phase inside the material system and the external free gaseous phase, respectively.b is the inverse of k B T, where k B is the Boltzmann constant.
The adsorption energy is marked as E ads in in Eq. ( 1).The fugacity, f ¼ /P, is introduced in the calculation as an effective pressure to describe the chemical thermodynamic process.The fugacity coefficient / can be fixed at particular temperature T and pressure P by ln /ðT; PÞ ¼ 1 RT where V m denotes the molar volume of the gas, and R is the Avogadro constant.The molar volume V m can be numerically determined with the aid of the virial expansion of the EOS for a real gas: where BðTÞ; CðTÞ and DðTÞ are the second, third and fourth virial coefficients, respectively.According to our previous work [19], until the fourth virial term is taken into account, the accurate fugacity coefficients of the normal hydrogen gas can be obtained compared with the experimental values.Details for how to determinate these virial coefficients can be found in our previous work [19].
A Gaussian type effective potential is introduced by Feynman and Hibbs (FH) in their book [20] to get rid of the quantum derivation between the experimental results and the Lennard-Johns (LJ) simulations.The Taylor expansion until second order for FHeffective potential [21] is employed to describe the body-body interaction between the hydrogen moleculs, as follows: where m r is the reduced mass of the LJ pairs.In this work, m r is set to 1:6744 Â 10 À27 kg [22] for the diatomic hydrogen molecule, which is treated as an entity.The LJ-parameters used here are set exactly to be the same as in the previous work [19].Naturally, the adsorbed and free hydrogen molecules are in two different states.However, the difference is very small owing to the essence of the physical adsorption.The electron density maps of an adsorbed H 2 molecule at certain typical positions on a monolayer h-BN are illustrated in Fig. 1.Note that the atomic configurations of the adsorption systems are fully relaxed during the calculations.According to our calculations, the largest deviation of the Mulliken charges between the adsorbed and free hydrogen molecules will not be larger than 4.7%.The change of the Coulomb interaction owing to this charge perturbation will be smaller than a percent.Considering the tiny variance, our simulations are performed with the free hydrogen model.We assume also that the hydrogen fluid keeps to be the normal form under our study conditions.The effect of the spin isomers ratio can be covered in the estimation errors.
If the hydrogen storage system is considered as a continuum model, the adsorption energy of a single hydrogen molecule in the eh-BN system can be introduced in Eq. ( 1) as E ads in , which is defined as follows: where the X denotes the hydrogen storage material system, and E tot ðXÞ represents the total energy of the material system, E free ðH 2 Þ points out the energy of a free hydrogen molecule, E tot ðX þ H 2 Þ is the total energy of the material system X with an adsorbed hydrogen molecule.Note that the basis set superposition error (BSSE) is already taken into account in our calculations.Here, the energy quantities are obtained from the calculations performed by the CRYSTAL14 computer code [23].Three all electron Gaussian basis sets of triple-zeta valence with polarization quality (TZVP) are applied for boron, nitrogen and hydrogen atoms, which are developed by Peintinger et al. [24].It is worth mention that the adsorbed hydrogen molecules may have different spatial forms, e.g., standing or lying on the h-BN surface.Nonetheless, the physical adsorption energies of the hydrogen molecule with different spatial forms have the differences less than 10 meV.Therefore, the spatial forms of the hydrogen molecule will not affect our prediction results.On the other hand, the change of the molecule orientation in a dynamic system is caused by the dynamic energy (in other words, the tem-perature of the system) and the influence from the surrounding molecules.Precisely, the most part of the orientation effect has already been taken into account in the fugacity calculations.A 3 Â 3 supercell is created for the eh-BN systems.A 32 Â 32 Pack-Monkhorst net [25] is set to sample the two-dimensional Brillouin zone in the reciprocal space.A thick vacuum layer is automatically builded via setting up c ¼ 500 Å by the code for all calculations to eliminate the effects from the periodic images.The dispersion correction method, named D2 introduced by Grimme [26], is also employed to get a reasonable adsorption energy.The dispersion coefficients C 6 are set to 3.13, 1.23 and 0.14 J nm 6 mol À1 for B, N and H atoms, respectively.The atomic van der Waals radii are set to 1.92, 1.55 and 1.20 Å for B, N and H atoms, respectively.Under this arrangement, our calculated lattice parameters of h-BN are a ¼ 2:497 and c ¼ 6:659 Å after full geometry optimization at PBE-D2 level.Compared to the experimental measurement for the lattice constants [27], the related deviations from our simulation are less than 1%.Five different inter layer spacings: 5, 6, 7, 8 and 9 Å between the expanded h-BN layers are considered.The interactions between the h-BN layers are very weak owing to the large distances.Therefore, the prediction process can be simplified through stacking the h-BN sheets as hexagonal (AA) arragement, as shown in Fig. 2(a).Actually, it is not only a simplification in our simulations, but also part of practice stacking forms in porous h-BN materials.Since h-BN contains two different species, the distributions of the electron and the electrostatic potential on h-BN surface are heterogeneous.Thus, the adsorption energy of the hydrogen molecule is also nonuniform in the fixed plane.However, taking full symmetry advantage of the B-N six-ring, the distribution map of the adsorption energy can emerge via sampling an equilateral triangle area as shown in Fig. 2(b).A net with 21 points is chosen for the sampling in an 1/6 BN-ring.For the entire area fenced by the ring, one only needs to roll over this sampled image another five times along its side.

Results and discussions
The EOS from virial expansion (VE) is used in this work to pave the way for determining the fugacity coefficients of the hydrogen fluid at arbitrary temperature and pressure.According to Eq. ( 1), the pressure inside the system can be estimated with the aid of the adsorption energy of the hydrogen molecule, when the fugacity coefficients of the hydrogen gas are already on hand through Eq. ( 2).Energy distribution information for the sample adsorption in the 1/6 BN-ring has been investigated.The contour map of the hydrogen adsorption energy in the middle of two h-BN layers with 7 Å interlayer spacing is presented in Fig. 3 as an example.It clearly demonstrates, that the adsorption energy in the center of the sixthring is the lowest as compared with the other points in the same intermediate plane.Furthermore, the nitrogen atoms are more attractive to the hydrogen molecule than the boron atoms on this plane.Even so, the scan along other planes (which approaches close enough to one of the eh-BN layers) shows that the adsorption energy in the center of the sixth-ring appears to be the lowest, but the boron atoms turn out to be more attractive than nitrogen, as shown in Fig. 4. Note that the calculated results until now are from the PBE-D2 calculations without any thermodynamic information.According to the definition of adsorption energy by Eq. ( 5), the negative energy points out the attraction and the positive one describes the repulsion.Fig. 4 shows that the values at different sampling points on the planes between À0.6 and 0.6 Å are almost the same.This region takes up the most of the integral area.Therefore, the average values of the adsorption energies are employed in Eq. ( 1) to estimate the pressure distributions inside the systems.Herein, the thermodynamic parameters (i.e., temperature and pressure) attend into our estimations.
Once the pressure information inside the eh-BN systems is on hand, the density profiles of the physisorbed hydrogen gas can be obtained from our virial EOS as defined by Eq. (3).As examples, the density profiles for the systems with 6, 7 and 8 Å interlayer spacings are presented in Fig. 5.It has been observed that the adsorption capacity on the middle plane becomes lower and lower with an increase in interlayer spacing.A double wall form is found in the system with large interlayer spacing (i.e., d P 8 Å), which illustrates that the large interlayer spacing leads to the gradual separation of the eh-BN system into two individual monolayer h-BN systems.After the integration of the density profile, the uptake of corresponding eh-BN system at certain thermodynamic condition can be obtained.The results are plotted in Fig. 6 for three meaningful temperatures 298, 273 and 243 K, respectively.The uptakes of the systems with d = 8 and 9 Å are not the highest ones in the low pressure region.With the increase of the external pressure, the uptakes of these two systems will exceed the other ones with smaller interlayer spacings.The temperature is another key parameter.The lower the temperature, the earlier the cross points appear.The estimated uptakes are also fitted with the aid of exponential function model to predict the uptake limit of the system at certain temperature.However, the predicted limits are based on the assumption that the physisorbed hydrogen gas will not change its phase inside the system under high pressures beyond 20 MPa.The uptake limits predicted by extrapolating the fitting curves are listed in Table 1.It is obvious from the density profiles that the uptake is also a function of the interlayer spacing within certain rang.The relation of uptake to interlayer spacing at T ¼ 243 K is presented in Fig. 7.The extrapolations of the fitting curves show that the perfect eh-BN system can only provide the maximal uptakes 2.24, 2.96 and 3.86 wt.% at 243 K under external pressures 5, 10 and 20 MPa, respectively.Our results have established that hydrogen storage capacity of the perfect eh-BN is no more better than the corresponding expanded bilayer graphene system.It might be the reason that the attraction of the in-plane polar BN-bond can not fully compensate the absence of the delocalized p-electrons on the surface.The predicted uptakes in present work are reasonable compared with the experimental result of the oxygen doped BN nanosheets reported in Ref. [8].At any rate, our prediction method for the hydrogen storage capacity is proved to be effective for some complex systems, e.g., eh-BN systems.The doped h-BN or graphene systems can also be similarity treated without any difficulty.It is worth mentioning that Li and co-workers experimentally showed that the hydrogen uptake capacity of a high-quality microporous/mesoporous BN material (named HBBN-1) with an exciting weight adsorption up to 5.6%

Table 1
The uptake limits of five expanded bilayer h-BN systems from the extrapolation of the fitting curves in Fig. 6   at room temperature and at the relatively low pressure of 3 MPa [28].This result exhibits a bright future of the development of hydrogen storage materials.After a close examination, we find out that the adsorption energy of a hydrogen molecule on the top of a nitrogen vacancy in a monolayer h-BN system is around 469 meV which is above five times stronger than that on the perfect h-BN surface.This result confirms one of the deduction for the high hydrogen uptake capacity in Ref. [28] that the high density of structural defects provides many more active sites for the adsorption which consequently enhance the hydrogen storage in the HBBN-1 system.

Conclusion
Employing the continuum model and our virial expanded EOS, we have calculated the hydrogen uptakes in the perfect eh-BN systems with different interlayer spacings at a series of temperature and pressure.Our method has been proved to be valid for the prediction of the hydrogen storage capability of h-BN based nano-material.Our findings clearly evidenced that the previously reported hydrogen storage capacities of the h-BN models (around or even higher than 7 wt.%)[5,9] are far away too high.The maximal hydrogen uptake is around 2.96 wt.% at the thermodynamic condition T ¼ 243 K and P ¼ 10 MPa, which only meets 2/5 of the applicable standard 7.5 wt.% in automobile industry [1].However, impurity doping in the layer of h-BN, alkali metal atom decorating on the surface of h-BN and the porous BN materials with high defect density may be the best solutions to attain the goal of the practice hydrogen storage in the material systems based on h-BN.

Fig. 1 .
Fig. 1.The charge density maps of the adsorbed H 2 molecule at four typical positions on a monolayer h-BN.The contours are varied from 0 to 0.2 e=bohr 3 with increment 0.005 e=bohr 3 .

Fig. 2 .
Fig. 2. (a) A sketch for the AA stacking of the bilayer h-BN model with certain interlayer spacing.The energy differences between different stacking models are very tiny owing to the large interlayer spacings involved in this work; (b) the arrangement of the sampling points in a single B-N hexagonal ring.

Fig. 3 .
Fig. 3. Adsorption energy distribution for single H 2 molecule in the eh-BN model with 7 Å interlayer spacing.The sampling is performed in the middle plane between two h-BN layers by using PBE-D2 method.The adsorption energy changes from À156.3 meV at the top of B-atoms to À159.5 meV at the top of the six-ring center with increment 0.2 meV.

Fig. 4 .
Fig. 4. Adsorption energy profiles for single H 2 molecule adsorbed at four typical sampling points in the eh-BN model with an interlayer spacing 7 Å.The squares, circles, triangles and stars represent the sampling points on the tops of boron atom, nitrogen atom, center of a six-ring and the B-N bridge, respectively.The change of the distance between the adsorbed molecule and one of the h-BN layers is set as abscissa.The zero position is defined on the middle plane of the bilayer system.

Fig. 7 .
Fig.7.The dependence of the hydrogen uptakes on the interlayer spacings at 243 K and 5 MPa (squares), 10 MPa (circles), 20 MPa (triangles).The solid lines are fitting curves by using exponential functions to find the limits of the hydrogen uptakes.
at three typical temperatures 298, 273 and 243 K.