Modeling and stabilities of Mg/MgH2 interfaces: A first-principles investigation

Modeling and stabilities of Mg/MgH2 interfaces: A first-principles investigation Jia-Jun Tang,1,2,3 Xiao-Bao Yang,1 Li-Juan Chen,1 and Yu-Jun Zhao1,2,a 1Department of Physics, South China University of Technology, Guangzhou 51064, China 2Key Laboratory of Advanced Energy Storage Materials of Guangdong Province, South China University of Technology, Guangzhou 510640, China 3School of Materials Science and Engineering, South China University of Technology, Guangzhou 510640, China


I. INTRODUCTION
Magnesium (Mg) is a promising candidate 1, 2 for application of mobile hydrogen storage material due to the high gravimetric hydrogen content 3 (7.6% in MgH 2 ), high volumetric density 3 (6.7 × 10 19 H atoms/L), and high energy density 4 (9MJ/kg) of its hydride MgH 2 . To make Mg applicable, there are still two technical obstacles 1, 2 to be solved: (1) Slow hydrogenation/dehydrogenation kinetics; (2) A too high working temperature of 350 -400 • C. In order to improve the material hydrogen sorption performance, extensive theoretical and experimental researches including strain engineering, 5,6 material surface/bulk ratio increasing, 7-9 transition metal catalyzing [10][11][12] and mixing 13,14 with other compounds have been done to mitigate the thermodynamic and kinetic drawbacks of Mg.
Hydrogen sorption in Mg concerns to the Mg/MgH 2 or other related interfaces. Recently, some experimental and theoretical efforts have been made to study the Mg/MgH 2 interface structure and improve MgH 2 dehydrogenation properties. Epitaxial growth of MgH 2 (110) on Mg(0001) is experimentally observed. 15 The significant increase of equilibrium hydrogen pressure in Pd/Mg systems was first believed to be caused by the clamping effect 16 but later found to be induced by the interfacial alloy hydride. 17 Mooij et al. 9 have found that the nanoconfined MgH 2 is more easily destabilized compared to its bulk counterpart because of an interface energy effect. Moreover, decreased dehydrogenation enthalpies could be also observed in Mg/Mg 2 Ni interface system 18  by the following formula where E tot (n) and E bulk represent the total energies of the slab model and the bulk with n formula units, A is the area of the surface slab and the factor 2 denotes two surfaces being created. Of note, the Fiorentini approach 36 is adopted to get a reasonable convergence of surface energies. 37

The interfacial energy
The interfacial energy E i at the grain boundary between two phases is defined as the formation energy of the interface system in a particular context, where the chemical potentials at which two phases are at equilibrium could be determined. The calculated formula of interfacial energy proposed by Arras et al. 38 is given by where E tot is the total energy of interface model, N Mg and N H are the total number of atoms of magnesium and hydrogen respectively in the interface system, μ 0 Mg and μ 0 H are the chemical potentials of magnesium and hydrogen when two phases, namely Mg and MgH 2 , are at equilibrium. In our model in which the two phases represent the atom reservoirs, 39 μ 0 Mg and μ 0 H are calculated by and where ε is the total energy per formula unit and x is the atomic concentration of an element (subscript) in a phase (superscript are equal to 1, 0, 1/3 and 2/3, respectively.

The H defect formation energy
The H defect formation energy E f is defined as the formation energy of the H vacancy in the MgH 2 phase, or the H interstitial in Mg phase of the interface model. The calculated formula of the defect formation energy is given by where E tot and E perfect are the total energies of interface models with and without H defect, N defect indicates the number of H defects, and the signs + and -represent H vacancy and interstitial, respectively. Here, in order to stimulate the H-rich context, μ H is equal to 1/2E H2 .

Surface energies
The surface energies strongly reflect the possibilities of formation of a particular surface. It should be noted that our first-principles calculation of surface energies is conducted at 0K, where the bond cutting and the surface reconstruction are the two major factors that dictate the surface energies. Previous papers 37,40,41 showed that both of these two effects played an important role on surface energies. Moreover, the contribution of vibration entropy, another factor in determining the surface energies, goes up when the temperature rises so that the two major factors may not be the only determining ones and the surfaces with high surface energies at 0 K may appear. 42 As a matter of fact, the hydrogen desorption in magnesium hydride requires a temperature much higher than the room temperature. Therefore, the surfaces with high surface energies at 0 K also have to be considered in dictating the combination of two surfaces of two phases.

The matching of shape and lattice constants of the cells
The periodic interface model which consists of two phases requires that the two phases share the same shape and close lattice constants. The perfectness of matching leads incoherent interfaces to coherent or semi-coherent interfaces.
The matching of shape takes priority over the matching of lattice constants. Specifically, the cells of rectangle shape (NB: the hexagonal shape can be transformed into rectangle shape) could always be matched easily after some unit repeating operations, which is in pursuit of existence of a coincidence site lattices 43,44 (CSL). Comparatively, the cells of parallelogram with unequal neighboring sides and a non-right angle are difficult to be matched with the others. For instance, it is found that the Mg(1011) surface cell is in its difficulty to be matched due to its shape of parallelogram where a Mg(1011) = 6.11 Å, b Mg (1011) = 6.38 Å and the included angle γ = 74.87 • . The cell of this shape, unlike the rhombus, is also impossible to be transformed into the periodic cell of rectangle shape due to its specified atomic arrangement.
The unit repeating operations follow two geometrical criteria: 45,46 (1) The lattice misfit strain between two phases should be as small as possible.
(2) The size of repeated cells, manifested by the period of the CSL, is as small as possible. The first criterion is easily understandable. The second one is based on the observation 43, 45 that the grain boundaries between two identical materials are often found to have a CSL with small unit cell volume. The information about the size of our repeated cells is specified in Table II. In summary, the rectangle cell (original or transformed) is always preferred in the process of choosing the surface of each phase. In practice, we choose the Mg(0001), Mg(1010), MgH 2 (101) and MgH 2 (210) surfaces whose information is specifically listed in Table I, which will be further discussed in following section. Fig. 1 provides an intuitive illustration of interface modeling and interface pairs in which we are interested.  Table II.

Determination of the mutual lattice constants
For an interface system with lattice mismatch, both of the two phases would be possibly under elastically tensile or compressive effect. Therefore, the elastic cost subject to the mutual lattice constants is required when the interface is created. Here, the mutual constants are defined as the lattice constants of the interface model where two phases coexist with minimized elastic energy cost. The minimum of this elastic energy cost is expected to be found in our interface models. Basing on the assumption that both of two phases would compromise energetically to reach the mutual constant, we propose that the energy cost on elastic effect is calculated by where E EEC represents the total elastic energy cost of two phases, E Phase.1 and E Phase.2 denote the elastic energy costs of two phases with respect to the most energetically favored structure. The search of the mutual lattice constants is required to determine the mutual constants. In practice, the calculations of E Phase.1 and E Phase.2 are calculated respectively in bulk forms of two phases by fitting the E-c curve at given a and b. Thus, the mutual constants could be determined by the minimized elastic energy cost. It is in line with widely adopted theoretical models, although it is relatively simple to the real interface where the lattice constants of both phases will be restored beyond corresponding critical thickness. 47

Determination of relative position between two phases
The stability of interface, after the determination of the mutual lattice constants, is highly influenced by the bonding nature in interfacial region, namely the interfacial structure. Naturally, the factors that lead to a change of interfacial structure should be focused on. One may easily come up with several factors including the relative position of two phases, the concentration gradient between two phases and the effect of defect or dopants in interfacial region.
The authentic hydrogen sorption process always exhibits migration of grain boundaries, which is affected not only by the thermodynamic stability but also by the kinetic one. From the point of view of theoretical calculation, it is not easy to describe all these factors simultaneously. However, when it comes to the discussion of the interfacial bonding, we simply assume that the relative position of two phases should be the primary factor in affecting the interfacial structure in our calculations since the variation of the relative position would lead to a dramatic change of interfacial bonding, whose effect is similar as the sliding shear acted on the system.

A. The bulk properties of magnesium and magnesium hydride
The crystal structure of magnesium is hexagonal close packing and the calculated Mg bulk lattice constants are a Mg = b Mg = 3.190Å with the c/a value of 1.623, which are in good agreement with the experimental values 48,49 and previous theoretical results. [50][51][52] Magnesium hydride, whose crystal structure is simple tetragonal, is a highly stable hydride with strong Mg-H bond, which is a mixture of ionic and covalent bonding. 53 The calculated bulk lattice constants are a MgH2 = b MgH2 = 4.452 Å, c MgH2 = 2.992 Å with the c/a value of 0.672, which agree well with the previous experimental [54][55][56] and theoretical 57-59 results.

B. The surface energies of Mg and MgH 2
Here, surface energies of two phases are calculated and compared to find the relative stabilities of surfaces. It is known that the low-index surfaces are often more stable than the high-index ones since the number of breaking bones 60 that are required to create a low-index surface is relatively less. However, it has been reported 37 that the low-index Mg(1010) has a high surface energy thanks to its massive basal bond cuttings. MgH 2 (210) is of the highest surface energy in Table I, while the other MgH 2 surfaces share low surface energies.
In a theoretical calculation, the adopted model should be periodic. Consequently, the shapes and lattice constants of two phases are identical, which are widely accepted in previous literatures 24,26,38,[61][62][63] concerning interface. It means that the shapes and the lattice constants should be taken into account in combining the two surfaces of two phases.
Considering these factors mentioned in the above two paragraphs, i.e. the surface energies, the shape and lattice constants of a surface-slab cell, we choose the combinations, as listed in Table I, of which two surfaces of two phases are identical in shape and close in lattice constants (repeated along the axes if necessary). In practice, the low-energy Mg(0001) and MgH 2 (101) are chosen while the high-energy Mg(1010) and MgH 2 (210) are selected to exemplify our approaches.

Elastic effect
Formation of interfaces always exhibits energy compromise of two phases, which can be reflected by elastic effect. We define here, for an ideal condition, the mutual constant as a constant which leads to minimized energy cost on elastic effect. Therefore, the search of mutual constants for interface models could be interpreted as searching the minimum E EEC among a wide range of lattice constants. In practice, the search of mutual lattice constant for interface could be divided into two steps: (1) The elastic costs on each phase are estimated. (2) Take into account the whole contribution from two phases caused by elastic effect. It should be pointed out here that the interface-induced strain effect on both phases may lead to different strains along axes a and b. Therefore the consideration of nonequivalent strain on axes a and b makes sense. Since the total energies of Mg and MgH 2 per formula unit are −1.56eV and −8.91eV respectively, the contributions from two phases, caused by elastic effect, also differ. It is clear that the thickness of the two phases significantly affects the minimization of the elastic energy in our model. Practically, the critical thickness, including the strain release effect of the layers within the critical thickness, should be included in the minimization of the elastic energy. Considering the difference in contribution of two phases, the energies of E Mg and E MgH 2 are calculated with respect to the most stable structure within the searching range. Moreover, the number of Mg layers is nearly twice as that of MgH 2 layers so that the contributions from two phases could be competing. Here, we assume that our choice of thicknesses, made by the two above mentioned strategies, is capable of investigating the elastic effect, i.e. determining the mutual constants. Similar ratio between numbers of two phases, i.e. Ge and GeO 2 , is used in previous literature. 64 Fig . 2 shows the case of Mg(1010)/MgH 2 (210) interface. We apply various structures, which are distinguishable from their lattice constants, to two phases in their bulk form cell. The elastic effects are found to be obviously influencing the phase stabilities. The relative stabilities of varying lattice constants could be reflected by the energy difference of the respective structures with respect to the most energetically favored one (as shown in Figs. 2(a) and 2(b)). Considering the Possion effect, we fixed a and b but varied c, and then fitted the E-c curve to determine the energy for each structure. Afterwards, the whole elastic contribution from two phases are estimated by simply calculating E EEC = E Phase.1 + E Phase.2 [see Fig. 2(c)]. The most energetically favored structure could be found among many structures that we considered. In this Mg(1010)/MgH 2 (210) interface case, the mutual constants are found to be a = 3 Å and b = 10.2 Å. Here, we have made a cell repeating operation for Mg(1010) in order to match the interface (see Table II).
Similarly, the mutual constants searching proceeds on the cases of Mg(0001)/MgH 2 (101) and Mg(1010)/MgH 2 (101) interfaces as shown in Fig. 3. The details of three interface combinations are specifically listed in Table II. With the obtained lattice constants, the interfaces could be modeled.

The search of relative position of two phases
The determination of mutual constants mainly focuses on elastic effect but does not involve the interfacial bonding, which primarily affects the interface stabilities. Recent papers have shown that interface properties highly depend on the lattice matching 65 or interface chemistry. 66 At Mg/MgH 2 interfaces, the Mg-H bonding is the most important bonding that leads to interface stability.
Authentically, there are interfaces consisted of two phases with distinctive symmetries in many cases, especially for the cases with high-index surfaces. For example, correlation between interface trap densities and interface anisotropy are found 67 on the interfaces consisted of high-index surfaces. Another case 68 is that the interface of Ag/high-index substrate enhances resonance by enhancing the scattering cross-sections. From a point of view of theoretical calculation, it is rather difficult to find the most energetically-favored structure for interfaces where two phases are with different symmetries. Therefore, the full consideration of relative position is essential.
We propose that the relative position of two phases could be seen as the vital reason causing change of interfacial bonding in theoretical calculations. For theoretical modeling, considerable interfacial structures could be obtained by adjusting the relative position of two phases. Therefore, it is possible to compare the stabilities of relaxed interfaces with various structures.
We have calculated the effect of relative position in Fig. 4. It could be obviously observed the varying of relative position could cause a dramatic change in interfacial energies (the largest difference is 73.2 meV/Å 2 (1.171 J/m 2 ), up to 22.73% of the maximum). In every structure (at particular x and y), an energy dependence on d could be seen. Here, d is defined as the distance between two phases of initial configuration along z axis, as shown in the inset of Fig. 4(c). This is easy to understand since the increasing of d lengthens the distance of two phases and thus greatly alters the Mg-H bond lengths at interface. It would lead to the dramatic change of Mg-H bond energies and thus influence the interfacial stabilities. We find there are local energy minimums in every structure within the proposed d range. It means that there is a preference for d, which dovetails with strengthening of interfacial bonds, in interface models. Logically, it could be interpreted that the optimization of d is always required for interface models. Moreover, it is found the varying of relative position on x-y plane (at fixed d) also induce a change of interfacial energies although not as dramatic as that induced by d varying. Specifically, the varying of x-y relative position would cause a change of local structure. But this effect would be obviously offset by relaxations of interfacial atoms and thus just induce a slight effect on the change of interfacial energies. Comparably, the varying of d would be more effective on tuning interfacial bonding. Higher effectiveness is mainly attributed to the fact that all interfacial bonds have to respond as d varies. Compared with x-y relative position varying, it is more global, which relaxations of interfacial atoms are not easy to offset. Similar findings could be observed at interfaces of Mg(0001)/MgH 2 (101) and Mg(1010)/MgH 2 (101) interfaces (see Fig. 5). In details, every structure experiences quadratic curves. Thus, the relative position of two phases, in theoretical interface models and calculations, should be always taken into account, especially for the influence of d varying on interfaces where no obvious symmetries are adopted.
Notably, it is also observed that the interfacial energies of three interface pairs in Fig. 4 and  Comparatively, the Mg-Mg bonding is much less contributive to interfacial stability. Therefore we argue that the Mg-H bonding is the primary cause leading to higher interfacial stabilities. The most energetically favored interfacial bonding structure could be obtained by varying the relative position, i.e. the search of relative position of two phases. Interestingly, we found that the interfacial stabilities are influenced by surface energies of two phases. The Mg(1010)/MgH 2 (210), the interfacial pair of surfaces with high surface energy, are expected to have low interfacial energy due to great number of dangling bonds that each phase has. However, the results show that Mg(1010)/MgH 2 (210) are owning a high interfacial energy around 4 meV/Å 2 (0.064 J/m 2 ) while the Mg(0001)/MgH 2 (101) pair consisted of two low energy surfaces has a much lower one of nearly 2.1 meV/Å 2 (0.034 J/m 2 ). Meanwhile, the interfacial energy of Mg(1010)/MgH 2 (101) pair, which is with one high energy and one low energy surfaces, is about 2.8 meV/Å 2 (0.045 J/m 2 ). This difference of interfacial energies could be explained by the fact that the partial saturation of interfacial dangling bonds is not capable enough to offset the instabilities of high energy surfaces. Experimentally, it is observed 44 that the grain boundaries comprised of low energy surfaces have relatively low energies. Our finding is in line with the experimental ones, and indicates that consideration of surface energies should be taken into when choosing interface combinations.

D. Defect formation at interface
The properties of Mg/MgH 2 interface, especially the hydrogen sorption property, are intriguing. As hydrogen sorption dovetails with the migration of Mg/MgH 2 interface, both of the hydrogen adsorption or desorption are considered. We here report comparative formation energies of hydrogen vacancies and interstitials, representing interfacial effect on hydrogen desorption and adsorption respectively. Fig. 6(a) shows the interfacial effects on hydrogen vacancies. Clearly, it could be observed that hydrogen vacancy formations are highly facilitated in near-interface region. Overall, the formation of hydrogen vacancies at interface is much easier than that in MgH 2 bulk. It could be well understood as the formation of interface dovetails with some dangling (under-coordinated atoms) or floating (overcoordinated atoms) bonds. The existence of dangling or floating bonds means the interfacial H atoms are less bound by Mg atoms and weaker stabilities than those in MgH 2 bulk and thus the interface dramatically lowers the formation energies of hydrogen vacancies.  effect affects the H vacancy formation energy as expected. As the interface lattice constants (tuned by strain) are close to those of MgH 2 bulk, the V H formation energies go up due to increased stability of MgH 2 phase of interface thus parabolic curves appear. Interestingly, the MgH 2 bulk data is found to be monotonic, which is not parabolic. Actually, it is raised from difference between interface and bulk data. Specifically, the C axis of bulk is allowed to relax when strains are applied. It means all the atoms are able to relax in response to the strain. Two MgH 2 atom layers in the interface model are fixed in order to stimulate the stable MgH 2 substrate therefore the atom relaxations in far-interface regions are not enough so that similar parabolic curves could be observed. In authentic interface environments, the latter situation is more likely to happen since inadequate responses happen to a substrate that are 'thick' enough, which we usually deem to be strains. Based on the above analyses, we conclude that hydrogen vacancies are more energetically preferred in interface regions, especially the near-interface one, compared with the bulk region. The interface also facilitates hydrogen interstitials. It is extremely easy for H inte 1st to form, which is reflected by its negative value of formation energy. H interstitial formations are much more approachable to the bulk region situation as it is more distance from interface. The reason for this phenomenon is that Mg bulk context would facilitate hydrogen uptake thanks to the strong Mg-H bonding. The strong bonding strongly increases stability of system. When a tensile strain is applied, the space for H interstitials is enlarged and thus resulting in a monotonic trend of formation energies. It is observed that the lengths of Mg-H bond under strains of −6% and +6% are 1.92 Å and 1.96 Å for H 4th inte . The altering of Mg-H bonding rationalizes the monotonic trend of formation energies.

IV. CONCLUSIONS
Three combinations of Mg/MgH 2 interface models are built to investigate the interface stabilities. The criterions including surface energies, interface matching, elastic effects and relative positions are discussed for Mg/MgH 2 interface modeling. It is found that the elastic effects and relative positions strongly affect the incoherent interface stabilities by affecting incoherent interfacial structure. Thus the searches of mutual constants and relative positions provide feasible solutions and make prediction of interfacial structure possible. Lastly, the defect formations are highly facilitated at Mg/MgH 2 interfaces.