Nb multilayer thin films

Computational Materials Science 108 (2015) 212–225

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Investigation of interfaces in Mg/Nb multilayer thin films A. Junkaew a,⇑, B. Ham c, X. Zhang b,c, R. Arróyave b a NANOTEC, National Science and Technology Development Agency (NSTDA), 111 Thailand Science Park, Thanon Phahonyothin, Tambon Khlong Nueng, Amphoe Khlong Luang, Pathum Thani 12120, Thailand b Department of Materials Science and Engineering, Texas A&M University, College Station, TX 77843-3003, USA c Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843-3123, USA

a r t i c l e

i n f o

Article history: Received 6 May 2015 Received in revised form 29 June 2015 Accepted 1 July 2015

Keywords: Mg/Nb thin films Interface bcc Mg Metastable DFT calculations

a b s t r a c t A density function theory (DFT) approach was used for understanding the stable and metastable interfacial structures observed experimentally in Mg/Nb multilayer films. Different types of interfaces, including the hcp Mg(0 0 0 1)/bcc Nb(1 1 0), the bcc Mg(1 1 0)/bcc Nb(1 1 0), and the hcp Mg(0 0 0 1)/hcp Nb(0 0 0 1), were investigated. The calculated interfacial energies help rationalize the stabilization of metastable interfaces involving bcc Mg-like arrangements in the experimentally characterized Mg/Nb multilayers. Electronic charge properties of the interfacial slabs are compared to corresponding bulk structures. The metallic bonding present in the thin film models has been explained by electron localized function and density of state analysis. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction Magnesium-based materials are some of the most promising light-weight hydrogen carriers investigated to date as Mg can store up to 7.6 wt.% of hydrogen [1]. However, practical use of conventional Mg-based bulk hydrogen storage materials is precluded by their rather high thermal stability as indicated by substantial hydride formation enthalpy, 76 kJ/mol, which results in dehydrogenation temperatures in excess of 573 K at hydrogen pressure of 1 bar [1]. In addition to the thermodynamic limitations, the kinetics of hydrogen sorption is rather slow. These hurdles make it difficult to use MgH2 as a hydrogen storage material for mobile applications. Noble metals and transition metals, such as Pd, Nb, Ti, V, Sc, Ni and others, have been used as catalysts for improving hydrogen sorption kinetics [2–9]. Moreover, improvement of hydrogen sorption in Mg has been investigated through thin film engineering. Higuchi et al. reported that capping Pd layer on top of Mg films decreases the dehydriding temperature to lower than 463 K. Reducing the degree of crystallization in the Mg/Pd layers has been shown to lower the dehydrogenation temperature, and hydrogen storage capacity reaches up to 5.6 wt.% [3]. In later work, the improvement in the dehydrogenation was found in Pd/Mg/Pd ⇑ Corresponding author. E-mail addresses: [email protected] (A. Junkaew), [email protected] (R. Arróyave). http://dx.doi.org/10.1016/j.commatsci.2015.07.003 0927-0256/Ó 2015 Elsevier B.V. All rights reserved.

trilayers. An increase in the thickness of Mg layer from 25 nm to 800 nm while keeping Pd at the 50 nm constant thickness is able to decrease their dehydriding temperatures from 465 K to 360 K with 5 wt.% hydrogen capacity [5]. Remarkable enhancement of the hydrogen sorption kinetics in the Mg/Pd thin films has been reported in other works, which used different thin films preparation and hydrogen loading conditions [10,11]. Moreover, other additives have been added to Mg-based thin films for improving their hydrogen sorption kinetic and thermodynamic properties, such as Cu [12], Ti [9], Mg–Ti–Fe [13], and Ni [14,15]. Nb has also been used to improve both the hydrogen sorption kinetics and thermodynamics in MgH2. The kinetics of hydrogen desorption in MgH2 with 5% NbH was reported by Pelletier et al. [16]. They proposed that an intermediate metastable NbHx hydride acts as catalyst for hydrogen sorption. MgH2 releases hydrogen rapidly at 573 K. de Castro et al. [17] suggested that the Mg–Nb nanocomposite can decrease the hydrogen desorption temperature to 543 K. The reduction in the desorption temperature indicates that the presence of Nb not only accelerates the kinetics of the hydrogen desorption process but also reduces the stability of the hydride. Hydrogenation in co-sputtered Mg–Nb thin films with various Nb compositions were determined by Mosaner et al. [18]. They found that these thin films provide enhanced kinetics in hydriding and lower activation energy in dehydriding when compared with pure Mg samples. More recently, the present authors [19,20] have shown that the stress associated with the epitaxial constraints in Mg/Nb

213

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

multi-layered thin films dramatically decrease the hydrogenation temperature MgH2 to 373 K. Detailed analysis of the structures of the hydrogenated films in combination with electronic structure calculations rationalized the dramatic de-stabilization of the magnesium hydride as a result of the formation of a metastable (orthorhombic) allotrope of MgH2, which is significantly less thermodynamically stable than the normal form of MgH2. In those and subsequent studies [21,22] the formation of metastable interfacial configurations involving bcc-like Mg and hcp-like Nb structures were observed depending on the thickness and composition (fraction of Nb/Mg) of the bilayer constituting the multilayer thin films. Metastable phases in Mg based nanocomposites films could play an important role on the hydrogen sorption properties of these films. This aspect has been largely overlooked in previous studies. In fact, metastable phases have been frequently observed in various metallic multilayer films, such as fcc Ti in the Ti/Ag multilayers [23], fcc Co and fcc Mn in Co/Mn multilayers [24], bcc Cr and hcp Co in Co/Cr multilayers [25], and fcc Co in Co/Cu multilayers [26]. The stabilization of the metastable polymorphs can be explained by favorable changes in surface and interfacial energies when considering classical bulk free energy differences between stable and metastable phases [27], the effect of coherency strains [28], and the effect of stacking faults in close-packed structures [29]. Ham et al. first observed metastable bcc Mg in Mg/Nb multilayers [30] and Junkaew et al. reported the stabilization of bcc Mg/bcc Nb and hcp Mg/hcp Nb interfaces [22]. In prior work we focused [22] on the investigation of the interface stability regimes for different thin film configurations and to do this we compared bulk and interfacial contributions to the interface system, relying on a Clausius–Clapeyron analysis to determine the coexisting lines in the so-called bi-phase diagram. This work focuses mostly on the structure, energetics and electronic properties of the stable and metastable interface structures in Mg/Nb multi-layers. The pseudomorphic growths of bcc Mg and hcp Nb in the Mg/Nb multilayers are investigated by performing the density function theory (DFT) calculations. Possible Mg/Nb multilayer models, which consist of hcp Mg(0 0 0 1)/bcc Nb(1 1 0), bcc Mg (1 1 0)/bcc Nb(1 1 0), and hcp Mg(0 0 0 1)/hcp Nb(0 0 0 1) multilayers and corresponding freestanding films are considered in this paper. Structural, energetic and electronic charge properties of these films are compared with those properties of bulk Mg and bulk Nb. 2. Methodology 2.1. Experimental method Mg/Nb multilayers were deposited on Si (1 0 0) substrates with 1 lm thermal oxides by DC magnetron sputtering at room temperature. The deposition rates for Mg and Nb were approximately 2 and 0.5 nm/s, respectively. The base pressure of the sputter chamber was 5 108 Torr prior to the deposition and argon partial pressure during deposition was under 2.5 mTorr. Mg and Nb targets, with 99.99% purity, were used. Transmission electron microscopy (TEM) specimens were prepared by grinding, polishing and finished by low energy (3.5 keV) Ar ion milling. TEM experiments were performed using a FEI Tecnai G2 F20 microscope operated at 200 kV. 2.2. Computational method The cohesive and structural properties of the interfaces in various systems have been explored through both computational and experimental work [31–35]. Two quantities used for understanding the mechanical and thermodynamic properties of the interfaces that can be directly extracted from DFT calculations [33,35]. The first quantity is the work of adhesion (W ad ), which is defined as the

reversible work needed to separate the interface into two free surfaces. Plastic deformation and diffusion are neglected for simplification. W ad represents the reversible work necessary to separate an interface and can be quantified through the following expression:

W ad ¼

1 ðEslðXÞ þ EslðYÞ  EslðX=YÞ Þ A

ð1Þ

where EslðXÞ and EslðYÞ are the total energy of isolated slab X and Y. The total energy of the interfacial slab is represented by EslðX=YÞ . A denominator (A) corresponds to the interface area. It is to be noted that in our calculations we considered isolated slabs corresponding to the ‘pure’ constituents as well as a composite slab corresponding to the interfacial structure. Effects due to the exposed surfaces on opposite sides of the interface are cancelled out as they appear in each of the states used in the expression above. The second quantity is the interfacial energy, ci , which is defined as the free energy difference between the interfacial atoms and interior atoms. This quantity describes the thermodymic stability of the interface. This term involves a more sophisticated calculation procedure. However, DFT calculations can be used for attaining this quantity as expressed in the following relation [33,35]:

c ¼ rX þ rY  W ad

ð2Þ

The surface energy by:

ri ¼

ri of the slab i (i = X or Y) can be calculated

    NslðiÞ 1 EbðiÞ EslðiÞ  2A NbðiÞ

ð3Þ

where EslðiÞ and EbðiÞ stand for the total energies of the surface slab and the bulk structure of component i, respectively. The number of atoms i contained in the slab- and the bulk structures are N slðiÞ and N bðiÞ , respectively. The factor of 2 is used because there are two free surfaces in the surface model separated by a vacuum region. The surface energy represents the stability of the surface. From the literature, these formulations have been used for understanding many interfacial systems such as metal–metal interfaces, and metal-ceramic interfaces [33,35,31,32,36,37]. In this work, all calculations were performed with a DFT approach implemented in the Vienna ab initio simulation package (VASP) [38–41]. The projector augmented method (PAW) [42] was used for a plane-wave basis. The local density approximation (LDA) [43] and the generalized gradient approximation (GGA) [44] were chosen for the exchange and correlation terms. In this work, 4

electrons in 4p6 4d 5s1 and 3p6 3s2 orbitals have been taken into account as the valence electrons of Mg and Nb elements. The plan-wave cutoff is 520 eV. As for bulk calculations, Monkhorst– Pack grids of 18  18  18 and 18  18  11 were employed for k-point sampling in the Brillouin zone of the bcc and hcp primitive cells, respectively. For the thin film calculations, the thickness of the simulated film is represented by the number of monolayers (ML). In all slab models, the slab is separated by a 20 Å vacuum layer in the direction normal to the surface. This vacuum prevents interaction between periodic images of the slabs. The k-mesh of 18  10  1 was used. Here, we considered three main types of interfacial models, which are the hcp Mg(0 0 0 1)/bcc Nb(1 1 0), the bcc Mg (1 1 0)/bcc Nb(1 1 0), and the hcp Mg(0 0 0 1)/hcp Nb(0 0 0 1) models. The interfaces were constructed based on the orientation relations observed by the experimental evidence [30], ½0 0 0 1hcp jj½1 1 0bcc ; 1  0 jj½1  1 1 , and ½0 1 1  0 jj½2 1 1  . The (1 1 0) stacking ½2 1 hcp

bcc

hcp

bcc

and (0 0 0 1) stacking are used in the bcc- and hcp-layers, respectively. In this paper, only coherent interfaces are considered. For example, the in-plane lattice parameters of the bcc Mg(1 1 0)/bcc

214

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

Nb(1 1 0) are laterally constrained to the lattice constants of bcc Nb. This simulated coherent interface relates to the HRTEM result in that the metastable bcc Mg has a lattice parameter similar to that of the bcc Nb layer [30]. Similarly, the in-plane lattice parameters of the hcp Mg(0 0 0 1)/hcp Nb(0 0 0 1) are based on the lattice parameters of hcp Mg. To account for the lattice mismatch in the hcp Mg(0 0 0 1)/bcc Nb(1 1 0) interface, coherent interfaces were constructed by assuming the matching of the in-plane lattice parameters of Mg or Nb layers to each other occurs. Then, there are two models that can be considered, which are the distorted in-plane hcp Mg (pseudo-hcp Mg) and the distorted in-plane bcc Nb (pseudo-bcc Nb). The Mg–Mg and the Nb–Nb interlayer distances are obtained from the spacing between two (1 1 0) planes in bcc Nb and the spacing between two (0 0 0 1) planes in hcp Mg, respectively. The interfacial spacing is an average value of the Mg–Mg and Nb–Nb interlayer distances in hcp Mg and bcc Nb. The interfacial bcc Mg/bcc Nb models were created by fixing 5ML-bcc Nb thickness and varying the thickness of the Mg layer, while the rest of the models were created by fixing 7ML-hcp Mg thickness and varying the thickness of the Nb layer. The adhesion of pseudomorphic bcc Mg on bcc Nb was monitored through calculations of W ad . Under the thin film condition, atoms are allowed to move along an out-of-plane direction, while they are constrained along the in-plane direction. Fig. 1(a)–(c) demonstrate the interface model of 5ML-bcc Mg(1 1 0)/5ML-bcc Nb(1 1 0) structure, 7ML-hcp Mg(0 0 0 1)/7ML-hcp Nb(0 0 0 1) structure and 7ML-hcp Mg(0 0 0 1)/5ML-bcc Nb(1 1 0) structures, respectively. Furthermore, the interactions between interfacial atoms were investigated by their electronic properties. The charge distribution, electron localization function (ELF) and density of states (DOS) were used to examine the interaction at the interfaces. The charge density and ELF contour of the unit cell were visualized by using VESTA [45].

3. Results 3.1. Experimental result Fig. 2 shows three different types of lattice structures of Mg/Nb multilayers. The formation of metastable phases in Mg/Nb multilayers can be observed in HRTEM micrographs and this is

Fig. 1. Slab models of (a) 5ML-bcc Mg(1 1 0)/5ML-bcc Nb(1 1 0) structure, (b) 7MLhcp Mg(0 0 0 1)/7ML-hcp Nb(0 0 0 1) structure and (c) 7ML-hcp Mg(0 0 0 1)/5ML-bcc Nb(1 1 0) structure. Orange and green balls represent Mg and Nb atoms, respectively. (For interpretation of the references to colours in this figure legend, the reader is referred to the web version of this paper.)

consistent with our recent studies [30,21]. Coherent interfaces with bcc structure were observed in Mg 5.6 nm/Nb 10.4 nm multilayers (Fig. 2(a)). The atomic arrangement in this sample is obviously different from the regular hcp/bcc structure observed in Mg/Nb 5 nm multilayers (Fig. 2(e)). Mg has the bcc Nb lattice parameters and showed epitaxial growth as confirmed by the fast Fourier transform (FFT) of the image in Fig. 2(b). This result also corresponds to our previous study on bi-phase diagram predictions [22]. On the other hand, Mg 1.8 nm/Nb 0.2 nm multilayers in Fig. 2(c) shows coherent hcp/hcp type of interfaces and the corresponding FFT of the image in Fig. 2(d) confirms hcp crystal structure of the specimen. Although Mg/Nb multilayers show coherency in their structure, the chemically discernible interfaces reveals that Mg and Nb form an immiscible system. 3.2. Bulk Mg and bulk Nb 3.2.1. Energetic and structural properties The optimized cell parameters of bcc Nb, bcc Mg, and hcp Mg are displayed in Table 1. LDA calculations show smaller cell constants than those in GGA calculations for all models. The results obtained in this work are in good agreement with literature as it is usual for LDA calculations to overestimate cohesive energies in metals, resulting in smaller lattice parameters. GGA calculations, on the other hand, tend to underestimate cohesion and result in larger lattice parameters. In the case of bcc Mg, the lattice parameters from literatures were characterized under high pressure condition [46]. Additionally, electronic charge properties of these bulk structures were also carried out. 3.2.2. Bonding and electronic charge properties In this section, quantification of the ELF is combined with charge density analysis in order to gain a better understanding of the electronic structure of the interfaces under study. ELF was first proposed by Becke and Edgecombe [50] as a useful tool to distinguish chemical bonding types in solids, such as metallic, covalent, and ionic bonds [51,52]. In this work, the electronic structures of the bulk phases were evaluated through valence charge density and ELF analysis. ELF and charge density contours projected on the (0 0 0 1) plane of the conventional hcp structure and the (1 1 0) plane of the conventional bcc structure are discussed in this section. Those properties in bulk Mg and bulk Nb are demonstrated in Figs. 3(a)–(d) and 4(a)–(d), respectively. Charge density plots in Mg and Nb structures are displayed in Figs. 3(a)–(b) and 4(a)–(b), respectively. Charge densities at atomic positions show the denser charge densities than outer regions. Changes of charge densities around atomic sites can be indicated by contour lines. For Mg crystal, their charge and ELF contours of hcp form are more distorted than those in bcc form. Fig. 3(a) and (b) show electron distribution directions and interactions between Mg sites in hcp-and bcc-Mg lattices, respectively. Metallic characteristics can be observed in their ELF. At the core electron region, ELF presents very low value related to low charge density in the pseudopotential regime. Apparently, ELF maximum is 0.55 appearing between an atom and its nearest neighbor atoms, which identifies the free electron behavior in bcc Mg. The ELF maximum on (0 0 0 1) plane is 0.69, showing the covalent-like behavior in hcp Mg. For bulk Nb, the slightly distorted shapes of charge distributions, presented in Fig. 4(a) and (b), indicate interactions among valence charges of Nb atoms. For typical metallic bonding, ELF value exists in range of 0.3–0.6 in the interstitial regions [53]. Then, the characteristic of metallic bonding also appears in the bcc Nb crystal shown in Fig. 4(d). Its ELF contour is influenced by the contribution of its 4d states. The complexity of d orbitals is reflected on the features of the ELF contour plot. The

215

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

Fig. 2. (a) A high resolution TEM image of Mg 5.6 nm/Nb 10.4 nm multilayers showing the formation of metastable bcc Mg, which is coherent with Nb at interface. (b) FFT of  1 1 zone axis. (c) HRTEM image of Mg 1.8 nm/Nb 0.2 nm multilayers and (d) corresponding FFT. Both Mg and Nb interfaces showing bcc diffraction pattern along ½1 1  0 zone axis. (e) HRTEM micrograph of Mg/Nb 5 nm multilayers showing regular hcp Mg and bcc Nb crystal structure components show hcp crystal structures along ½2 1 confirmed by FFTs in (f and g).

Table 1 The optimized cell parameters (in Å) of bulk structures by DFT-GGA and DFT-LDA. Structure

hcp Mg bcc Mg bcc Nb hcp Nb

Calculation

Experiment

LDA

GGA

a = 3.140, c/a = 1.622 (a = 3.13, c/a = 1.615 [48]) a = 3.508 a = 3.264 a = 2.827, c/a = 1.831

a = 3.197, c/a = 1.622 (a = 3.20, c/a = 1.660 [48]) a = 3.571 a = 3.324 a = 2.879, c/a = 1.828

non-uniform occupation of particular d orbitals affects the non-spherical symmetry of ELF contour [54]. According to the electronic states of bcc crystals, there are two irreducible representations of the d orbitals, which are the eg and the t2g states. The eg and the t2g states orient along [1 1 1] and [1 0 0] directions to the nearest-neighbor atoms. The bcc structure is more stable when more electrons are filled into the eg states. Thus bulk Nb exists in bcc form. Meanwhile bcc Mg can form under high pressure, since compression can stabilize the bcc Mg structure by lowering of its 3d band toward the Fermi level (EF) [55,56,21]. Both partial density of states (PDOS) and total density of states (DOS) were calculated for bulk Nb and Mg. The results are demonstrated in Fig. 5. A red dash line represents the Fermi level (EF). In Fig. 5, all Nb and Mg crystals have metallic behavior as there is no significant band gap separating between conduction and valence bands. DOSs of hcp Mg reveal the same behavior with a literature [57]. For bcc Nb crystal, the most important contribution to the DOS comes from electrons in d band. There are three sets of peaks existing in DOSs of bcc Nb, which are 4 eV to 2 eV, 2 eV to 0.5 eV, and 2 eV to 6 eV. This result is in good agreement with literatures [58–60]. Both Mg forms provide similar DOSs, while there are some differences in the calculated DOS of two Nb forms. More continuous change can be observed around EF in hcp Nb.

a = 3.22, c/a = 1.624 [47] a = 2.9530  0.002 [46] a = 3.300 [49]

3.3. Surface and interface calculations of Mg/Nb films 3.3.1. Energetic properties In this work, the interfacial properties of three major forms of interfacial structures in the observed synthesized multi-layer films, which are the hcp Mg(0 0 0 1)/bcc Nb(1 1 0), the bcc Mg(1 1 0)/bcc Nb(1 1 0), and the hcp Mg(0 0 0 1)/hcp Nb(0 0 0 1), are determined by DFT calculations. By neglecting the incoherent interfaces, four coherent interfacial structures were constructed with the optimized lattice parameter of bulk Nb and bulk Mg. The number of Nb layers is varied from 1ML to 5ML, while the number of Mg layers is fixed at 7ML, in case of the hcp Mg(0 0 0 1)/hcp Nb(0 0 0 1), the hcp Mg(0 0 0 1)/pseudo-bcc Nb(1 1 0), and the pseudo-hcp Mg(0 0 0 1)/bcc Nb(1 1 0). For the bcc Mg(1 1 0)/bcc Nb(1 1 0) models, the number of Mg layers is varied from 1ML to 5ML, while the number of Nb layers is fixed at 5ML. The surface energies in isolated slabs, work of adhesion and interface energies in interfacial structures are shown in Figs. 6 and 7, respectively. Note that in Figs. 6 and 7 we do not include the results of calculations involving the 2ML-hcp Nb configurations. During our study, we found that the 2ML hcp Nb configuration was highly unstable, with a calculated surface energy in excess of 6 J/m2, which is approximately twice as large as the surface energy calculated for the other isolated hcp Nb slabs. In fact, during structural

216

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

Fig. 3. Electronic charge structures projected on the hcp(0 0 0 1) and bcc(1 1 0) planes in Mg structures (a) the charge density distribution in hcp Mg, (b) the charge density distribution in bcc Mg, (c) the ELF contour in hcp Mg, and (d) the ELF contour in bcc Mg. For charge density distributions, the contour lines are drawn from 0 to 0.5 at 0.0015 e/ Å3 intervals. The contour lines are drawn from 0 to 0.80 at 0.05 intervals in ELF contours.

Fig. 4. Electronic charge structures projected on the hcp(0 0 0 1) and bcc(1 1 0) planes in Nb structures (a) the charge density distribution in hcp Nb, (b) the charge density distribution in bcc Nb, (c) the ELF contour in hcp Nb, and (d) the ELF contour in bcc Nb. For charge density distributions, the contour lines are drawn from 0 to 0.5 at 0.0015 e/ Å3 intervals. The contour lines are drawn from 0 to 0.80 at 0.05 intervals in ELF contours.

optimization, it was observed that this structure tended to relax into two isolated hcp-Nb monolayers. It is unclear whether this observation was due to a deeper physical reason or was more related to convergence and relaxation behavior of the code. In any case, following the trends for the calculations (work of adhesion and surface energies) corresponding to the 1, 3, 4 and 5ML-hcp Nb, it is clear that the interfacial energy corresponding

to the hcp Mg/2ML-hcp Nb system should be on the order of 0.5 J/m2, although these estimates are not rigorous due to the issues already discussed. Regardless of the uncertainty in this single point, the remaining calculations show systematic trends, as evidenced in Figs. 6 and 7. In Fig. 6, calculated r results of nML-bcc Mg, nML-pseudo-bcc Nb, and nML-bcc Nb are in the range of 0.78–0.99 J/m2, 2.67–

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

217

Fig. 5. Electronic DOS of (a) hcp Mg, (b) bcc Mg, (c) hcp Nb, and (d) bcc Nb.

Fig. 6. Surface energies,

r, in simulated surface models.

3.38 J/m2, and 2.36–2.96 J/m2, respectively. The bcc Nb(1 1 0) surface has been claimed to be the most stable one among other low-index crystalline orientations, such as (1 0 0) and (2 1 0) surfaces [61,62]. From the literature, the surface energy of Nb (1 1 0) is in the range of 2.4–2.9 J/m2 [61–66]. The 2.655 J/m2 of surface energy of this (1 1 0) plane was measured by experiment [67]. A slight increment of surface energies of the isolated bcc Mg(1 1 0) occurs as slab thickness increases. While a declining trend happens in a few isolated bcc Nb slabs, although the oscillations tend to stabilize as the number of layers increase. In case of the hcp Nb(0 0 0 1) and the distorted bcc Nb(1 1 0), their surface energies tend to increase as the number of Nb layers are greater than three layers. Typically, the surface energy is related to the stability of one particular surface structure, so the Nb(1 1 0) facet is more stable as its thickness increases. Differently, there is less stability in the Mg(1 1 0) facet with increment of its thickness, which may be related to the fact that Mg is not stable in the bcc structure under normal conditions. Similar to the Mg(1 1 0) facet, the metastable hcp Nb and the distorted bcc Nb facets are less stable as the number of layers increases. Despite the fact that in some of

Fig. 7. (a) Work of adhesion W ad , and (b) interface energies c in simulated interface models.

the systems studied the energies were not stablized in a rigorous manner, we believe that the variance in surface energy as a function of number of layers will not be significant and we thus take the values observed with the largest number of layers as a good approximation to the large ML number limit. Next, the interface models were explored. Plots of W ad and c as a function of varied n layers, which specified in each interfacial model are shown in Fig. 7. Results for W ad of the bcc Mg(1 1 0)/bcc Nb(1 1 0), the hcp Mg(0 0 0 1)/pseudo-bcc Nb(1 1 0), and the pseudo-hcp

218

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

Table 2 Selected data of percentage of interlayer spacing change (%Ddi;j ) of isolated Mg and Nb slabs. Structure

7ML-hcp Mg

7ML-pseudo-hcp Mg 5ML bcc Mg 5ML bcc Nb

5ML pseudo-bcc Nb 5ML hcp Nb a b c d e f g

%Ddi;j %Dd1;2

%Dd2;3

%Dd3;4

+0.89 +1.9  0.3a +1.76b +1.18c +1.13d +4.20 +4.41 3.89 3.60e 3.70f 4.3g 10.78 24.42

0.17 +0.8  0.4a +0.0b +0.36c +0.31d +3.05 +3.45 0.00 0.5e

0.37 0.4  0.5a +0.0b 0.73c +0.21d +3.18

4.93 16.42

Experimental data (LEED 100 K) [68]. Experimental data (LEED 100 K) [69]. Theoretical values: [70]. Theoretical values: [48]. Theoretical values: [71]. Theoretical values: [63]. Theoretical values: [60].

Mg(0 0 0 1)/bcc Nb(1 1 0) models are in range of 2.71–2.81 J/m2, 2.26–3.03 J/m2, and 2.71–3.25 J/m2, respectively. In the hcp Mg(0 0 0 1)/hcp Nb(0 0 0 1) models, W ad values are in range of 2.55–

3.31 J/m2, excluding the 7ML-hcp Mg/2ML-hcp Nb interface. Interestingly, all W ad values are greater than zero. Significantly, these positive W ad values demonstrate the attractive interactions between two layers at the interfaces in all simulated models, whether they are formed by metastable bcc Mg or metastable hcp Nb. However, W ad values of the hcp Mg(0 0 0 1)/hcp Nb(0 0 0 1) models decrease when the number of hcp Nb layers increases. In case of metastable bcc Mg on bcc Nb, their W ad are similar when the number of bcc Mg layers increases from 1ML to 5ML. Fig. 7(b) shows that the energy of the interface of the bcc Mg(1 1 0)/bcc Nb(1 1 0) increases with increasing of bcc Mg layers. The hcp Mg(0 0 0 1)/pseudo-bcc Nb(1 1 0), and the pseudo-hcp Mg(0 0 0 1)/bcc Nb(1 1 0) also have the similar trend when the number of Nb layers increases. While pseudo-hcp Mg(0 0 0 1)/bcc Nb(1 1 0) models have the similar c when the number of Nb increases. As a result, the interface energies in the pseudo-hcp Mg(0 0 0 1)/bcc Nb(1 1 0) models reveal similar values even if the number of bcc Nb layers change. Differently, the reductions of c with lowering n values are presented in the rest interfacial models.

3.3.2. Microstructures of the surface and the interface models In order to connect with the real (experimental) thin films conditions, the lattice parameters of the interface models are fixed so they correspond to the bulk forms, but they are allowed to relax along the normal direction of (1 1 0) and (0 0 0 1) planes for bcc and hcp layers, respectively. The relaxed interlayer distances di;j ,

Fig. 8. A percentage of interlayer spacing change %Ddi;j in (a) nML-bcc Mg/5ML-bcc Nb, (b) 7ML-hcp Mg/nML-hcp Nb, (c) 7ML-hcp Mg/nML-pseudo-bcc Nb, and (d) 7MLpseudo-hcp Mg/nML-bcc Nb, where n = 1–5.

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

where i and j represent two adjacent layers, were measured in both of surface and interface structures. The relaxation behavior are explained in term of the percentage of interlayer spacing change comparing with an initial value %Ddi;j , which is ðdi;j  d0 Þ 100=d0 . The d0 values are 2.547 Å and 2.308 Å for the hcp and bcc slabs, respectively. Calculated %Ddi;j of selected Mg and Nb slabs are reported in Table 2. Expansion and contraction of the interlayer spacings referenced to the bulk values are indicated by the positive and negative signs of %Ddi;j . Firstly, the microstructure of surface structures, which are isolated Mg and Nb slabs, are considered. It should be noted that these isolated films were simulated by placing slabs in the middle of cell and they are separated by 20 Å vacuum on both sides. Thus, both sides of slabs are symmetric. For instance, the properties of the first and the second layers of 5ML-Nb are identical to its fourth and fifth layers. Their structural changes after relaxation are reported in Table 2. The results in Table 2 indicate that both the surface and the in-plane constraint affect the relaxation behavior of the Mg and Nb slabs. The topmost interlayer distances of the Mg slabs are expanded comparing with the bulk values, and larger change can be examined in the pseudo-hcp Mg and the bcc Mg than that in the hcp Mg. There are both experimental and computational works investigating the interlayer relaxation in clean Mg(0 0 0 1) surfaces [48,72,68,70,73–75]. An expansion of the topmost interlayer spacing in Mg(0 0 0 1) was observed in all works, although in slightly different magnitude, while the next interlayer spacings show expansion or contraction behaviors in different works. From the

219

experiments, the relaxations of the first three interlayer spacing of the clean Mg(0 0 0 1)-(1  1) surface had been observed by using a dynamical low-energy electron diffraction (LEED) I–V analysis at 100 K [68]. They reported that Dd1;2 ¼ þ1:9  0:3, Dd2;3 ¼ þ0:8  0:4, Dd3;4 ¼ 0:4  0:5, Dd4;5 ¼ 0:0  0:5, where d0 value is 2.61 Å. For the relaxations of the next interlayer separations, which are %Dd2;3 and %Dd3;4 , are different from other reported values in Table 2. However, these contractions of %Dd2;3 and %Dd3;4 were present in theoretical investigations of the quantum-size effect (QSE) on Mg(0 0 0 1) surface by Li et al. [75] and Zhang et al. [66]. In contrast with the relaxations in the Mg(0 0 0 1) surface, a contraction of the topmost interlayer is found in the 5ML bcc Nb(1 1 0), which is comparable with the results from references [71]. Nevertheless, interlayer spacing of inner layers seems to be bulk-like for the 5ML bcc Nb(1 1 0). Additionally, compressions of their interlayer spacing can be monitored in all Nb slabs. The distances between a surface and the first subsurface layer %Dd1;2 are more compressed than next interlayer spaces %Dd2;3 . The most contraction is shown in the topmost interlayer spacing of 5ML hcp Nb. This shortening of interlayer spacing between the surface and subsurface layer will be discussed by their charge properties in next section. Fig. 8 shows the comparison of %Ddi;j as a function of varied layers. In Fig. 8(a), %Ddi;j of nML-bcc Mg/5ML-bcc Nb is presented. Not only the surface effect, but also the interface effect influences their structure. The last five layers are Nb layers in this case. After relaxation, all interlayer distances in Nb region are contracted, excluding

Fig. 9. Charge distributions and ELF contours projected on (1 0 0) planes in freestanding slabs (a) charge distributions of the 7ML-hcp Mg model, (b) charge distributions of the 5ML-bcc Mg model, (c) charge distributions of the 5ML-hcp Nb model, (d) charge distributions of the 5ML-bcc Nb model, (e) ELF contour of the 7ML-hcp Mg model, (f) ELF contour of the 5ML-bcc Mg model, (g) ELF contour of the 5ML-hcp Nb model, and (h) ELF contour of the 5ML-bcc Nb model. For charge density distributions, the contour lines are drawn from 0 to 0.5 at 0.0015 e/Å3 intervals. The contour lines are drawn from 0 to 0.80 at 0.05 intervals in ELF contours.

220

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

d2;3 which is expanded slightly. It can be observed that the interlayer distances in the first four Nb layers are constant, even if the number of Mg layers increases. However, the nearest Nb layer beneath an interface is influenced by the interface effect, with d4;5 increasing as the number of Mg layers increases. In contrast to the Nb region, distances of all upper Mg layers, including the interface region, are greater than the initial spacing. This can be explained by the lattice mismatch strain effect of the interface [76]. Elastic distortion of the epitaxial layers, i.e. Mg layers, can be observed when their in-plane lattice parameter is constrained to the lateral (in-plane) Nb parameters. Then, the epilayer lattice parameter normal to the interface relaxes along the out-of-plane direction. According the bulk calculations, the optimized lattice parameters in bcc Mg and bcc Nb are 3.508 Å (3.571 Å) and 3.264 Å (3.324 Å) from LDA (GGA) calculations, respectively. In this case, the lattice parameter of epilayer is greater than that of substrate, and then the epilayer lattice parameter expands along the interface normal. In addition, the topmost interlayer spacing exposes the most extension due to the free surface effect. All Mg interlayer spacings, including the interface, tend to decrease when the number of Mg layers increase. In Fig. 8(b)–(d), the first seven layers represent the Mg region and the next five layers are Nb layers. %Ddi;j in the 7ML-hcp Mg/nML-hcp Nb, the 7ML-hcp Mg/nML-pseudo-bcc Nb, and the 7ML-pseudo-hcp Mg/nML-bcc Nb are shown in Fig. 8(b)–(d), respectively. Most of Nb interlayer spacings in Nb regions are

contracted. Considering the Nb regions in these three models, %Ddi;j as a function of Nb layers change in the same way, but in different values or magnitudes. For instance, the topmost Nb interlayers have the most contraction, and then next sublayers are less contracted than those in the topmost Nb interlayers. In these three major models, the black filled circles are the interfacial interlayers, which their reference values (d0 ) are the average of the hcp Mg(0 0 0 1) and the bcc Nb(1 1 0) spacings. Comparing Fig. 8(b) and (c), their %Ddi;j as a function of Nb monolayers alter as the same trends, but the interfacial distances in the 7ML-hcp Mg/nML-hcp Nb models are closer than those in the 7ML-hcp Mg/nML-pseudo-bcc Nb. For the Mg regions in Fig. 8(b)–(d), %Ddi;j as a function of Nb layers of the 7ML-pseudo-hcp Mg/nML-bcc Nb models have most change than other models. Obviously, the interfaces dominate the change in Mg regions in these three models. Also, the lattice mismatch strain effect can be used for explaining these results. The interaction between two layers at the interface will be discussed in further section. 3.3.3. Bonding and electronic properties in surface and interface structures In order to identify the type of atomic bonding across the interfaces studied, the electronic charge properties of the interfacial models are presented in this part. In this section, the electronic charge density distribution, ELF, and DOSs analysis were evaluated to clarify interactions in slab structures.

Fig. 10. Charge distributions and ELF contours at the interfaces, projected on (1 0 0) planes in the interface structures (a) charge distributions of the 5ML-bcc Mg/5ML-bcc Nb model, (b) charge distributions of the 7ML-hcp Mg/5ML-hcp Nb model, (c) charge distributions of the 7ML-hcp Mg/5ML-pseudo-bcc Nb model, (d) charge distributions of the 7ML-pseudo-hcp Mg/5ML-bcc Nb model, (e) ELF contour of the 5ML-bcc Mg/5ML-bcc Nb model, (f) ELF contour of the 7ML-hcp Mg/5ML-hcp Nb model, (g) ELF contour of the 7ML-hcp Mg/5ML-pseudo-bcc Nb model, and (h) ELF contour of the 7ML-pseudo-hcp Mg/5ML-bcc Nb model. For charge density distributions, the contour lines are drawn from 0 to 0.5 at 0.0015 e/Å3 intervals. The contour lines are drawn from 0 to 0.80 at 0.05 intervals in ELF contours.

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

3.3.3.1. Charge distribution and ELF analysis of 5ML-Mg and 5ML-Nb structures. According the symmetry of both surfaces, the contour plots of selected surface slabs are presented from the surface region to the middle layer of structures. The charge distributions and ELF contours of selected freestanding slabs are illustrated in Figs. 9(a)–(d), and 9(c)–(h), respectively. In case of the Mg slab, the charge density distributions are distorted along an in-plane direction, which are shown in Fig. 9(a) and (b). Fig. 9(e) and (f) show the ELF contours corresponding with those charge distributions of 7ML-hcp Mg and 5ML-bcc, respectively. For all ELF contours, ELF value is zero in the vacuum region. The ELF maximum reaches 0.8 at the surface, and ELF values in an interior region also comprise higher magnitude than ELF values in bulk bcc Mg. This high ELF value is thought to be related to the covalent-like feature of metallic bonding which involves mixing of the s and p atomic valence orbitals. This result is supported by DOS analysis in a latter section. For surface structures, interactions between atoms in Nb structures can be observed in both of charge density distributions and ELF contours. Their charge distributions and ELF contours of

221

5ML-hcp Nb and 5ML-bcc Nb illustrate in Figs. 9(c)–(d) and 9(g)– (h). The distributions of charge densities among those atoms can be notified obviously. In 5ML-bcc Nb, interactions between atoms in the second Nb layer have and atoms in the surface layer are also observed in Fig. 9(d). Additionally, these valence charge features of Nb slab can explain the shortening of interlayer spacing between its surface and subsurface layers. Fu et al. [77] proposed effects of localized d electrons and delocalized sp electrons on the relaxation of a transition-metal surface. Considering bonding between surface atoms and subsurface atoms, the contraction of their interlayer distance happens as the d–d bonding is increased and the delocalized sp electrons are extended outward to the vacuum region. In Fig. 9(g) and (h), the highest ELF value, which is 0.5, can be observed at a surface region. This high ELF value notifies the delocalized sp electrons on the topmost layer.

3.3.3.2. Charge distribution and ELF analysis of Mg/Nb interfaces. Contour plots of valence charge density and ELF on (100) cuts of selected interfaces are illustrated in Figs. 10(a)–(d)

Fig. 11. The layer-projected DOS (LPDOS) of the 5ML-bcc Mg/5ML-bcc Nb model (a)–(e) the first Mg layer to the fifth Mglayer (f)–(j) the sixth Nb layer to the tenth Nb.

222

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

and 10(e)–(h). There are charge density gradient at the interface corresponding with the interaction across the interface. The valence electron density distributions and ELF contour in the 5ML-bcc Mg/5ML-bcc Nb model are illustrated in Fig. 10(a) and (e). Its valence electron densities are distorted along an in-plane in each layer. Similar with the surface models, Mg region has higher ELF than Nb region. The ELF maximum is 0.8 at the surface Mg plane, which is not shown here. For Mg, ELF decreases along the direction from surface to interface area. For Nb, ELF values less than 0.5 are found in inner layers, and its ELF maximum is shown at the surface. At the interface, the ELF values are close to 0.5 which is associated to free electron-like behavior in metallic bonding. The other three interface models (7ML-hcp Mg/5ML-hcp Nb, the 7ML-hcp Mg/5ML-pseudo-bcc Nb, and 7ML-pseudo-hcp Mg/5ML-bcc Nb), reveal similar characteristics in their charge structures. The interactions across the interface and metallic features can be elucidated as same as in the 5ML-bcc Mg/5ML-bcc Nb model.

3.3.3.3. Density of states (DOS). The DOS in interface structures were calculated and compared to the data for bulk Mg and Nb. Figs. 11–14 demonstrate the layer-projected DOS (LPDOS) of selected interface models, which are the 5ML-bcc Mg/5ML-bcc Nb, the 7ML-hcp Mg/5ML-hcp Nb, the 7ML-hcp Mg/5MLpseudo-bcc Nb, and the 7ML-pseudo-hcp Mg/5ML-bcc Nb models, respectively. LPDOS at the interfaces present different behavior from that observed in the bulk forms. In all interface structures, there are more filled electrons in electronic states than in their bulk states. In all Nb layers, LPDOS peaks are broadened near the Fermi level (EF ). For example, the LPDOS of each layer of the 5ML-bcc Mg/5ML-bcc Nb is illustrated in Fig. 11. Fig. 11(g)–(i) notify that interior layers pose in similar aspect. LPDOS in multilayers also change from those in bulk forms. Likewise, the electronic behaviors in the other models are similar with the 5ML-Mg/5ML-Nb model. In the multilayers, there are more states occupying the valence bands than that in their bulk forms. In part of Nb layers, their

Fig. 12. The layer-projected DOS (LPDOS) of the 7ML-hcp Mg/5ML-hcp Nb modell (a)–(g) the first Mg layer to the seventh Mg layer (h)–(l) the eighth Nb layer to the twelfth Nb.

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

223

Fig. 13. The layer-projected DOS (LPDOS) of the 7ML-hcp Mg/5ML-pseudo-bcc Nb model (a)–(g) the first Mg layer to the seventh Mg layer (h)–(l) the eighth Nb layer to the twelfth Nb.

LPDOS are broad around EF region compared to DOS of bulk Nb. The hybridization of sp orbitals are verified by overlapping of LPDOS in Mg layers. A result of this interface structure is similar to its isolated Nb and Mg slabs in the region of the surface to interior layers. A larger number of states in both of s- and p-bands participates in valence band than in conduction band. These aspects imply that forming Mg/Nb multilayers affects on electronic properties of its components. Both of surface and interface parts dominate on their electronic characteristics. 3.4. Final discussion on relative stability of metastable interface structures As mentioned in the introduction to this manuscript, the major motivation to carry out this work was the observation (see Fig. 2) of different interfacial structures in multi-layered Mg–Nb thin films with different fractions of Mg and Nb (f Mg ; f Nb ) and different bi-layer thicknesses (k). In prior work [22] the stabilization of these different interfaces was rationalized in terms of a so-called biphase

equilibrium analysis. Using thermodynamic coexistence arguments, it is possible to define the change in free energy of a bi-layer system as the phases composing the interface undergo phase transitions [22]:

Dg ¼ 2Dc þ ½DGA f A þ DGB f B k

ð4Þ

In this equation, Dc is the change in the interfacial energy upon transformation, DGi is the allotropic free energy change (normalized per volume) of the reference component, i. Note that DGi may include not only the difference in the thermodynamic bulk free energies of the phases but also include other volume effects, such as those associated to strain energy. The condition for coexistence is arrived at when dg ¼ 0. By comparing all possible interfaces (in this case hcp Mg/bcc Nb, hcp Mg/hcp Nb and bcc Mg/bcc Nb), using the coexistence condition (Clausius–Clapeyron) one can then find the slopes for the boundaries separating each interfacial regime. What was concluded from that study was that the competition between interfacial and bulk contributions to the energetics of the multi-layered system could be greatly affected by the structure

224

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

Fig. 14. The layer-projected DOS (LPDOS) of the 7ML-pseudo-hcp Mg/5ML-bcc Nb model (a)–(g) the first Mg layer to the seventh Mg layer (h)–(l) the eighth Nb layer to the twelfth Nb.

Fig. 15. Comparison of the predicted diagram of Mg/Nb multilayers and experimental results. The bi-phase diagram are plotted as a function of k and f Nb . The experimental results are represented by points [22].

and thickness of the bi-layers constituting the thin film. As per that model [22], it was indeed possible to stabilize a bcc Mg/bcc Nb interface, for example, provided the relative fraction of Mg to Nb was small enough and the total interlayer thickness was sufficiently small. From Fig. 7, it is seen that the energy of hcp Mg/bcc Nb interface, despite the fact of being the stable interface structure, is significantly higher than that of hcp Mg/hcp Nb and bcc Mg/bcc Nb interfaces, according (roughly) to the relation: chcp Mg=bcc Nb > chcp Mg=hcp Nb > cbcc Mg=bcc Nb . This same qualitative relation was found in the CALPHAD-based analysis presented by the present authors in [22]. The stabilization of hcp Mg/hcp Nb and bcc Mg/bcc Nb can be rationalized in simple energetic terms: as long as the system size is small (small total bi-layer thickness, k), bulk-like effects may not be as important as interfacial effects and in this case, the interfacial system with the lowest interfacial energy would dominate—for a more detailed discussion please refer to prior work [22]. This is observed in the experimentally determined biphase diagram for the Mg/Nb system as first reported in [22] and reproduced here in Fig. 15.

A. Junkaew et al. / Computational Materials Science 108 (2015) 212–225

4. Conclusion DFT methods were employed to study the stable and metastable interfacial structures in Mg/Nb multilayers. Calculated W ad values correspond to favorable interactions in all simulated interface models. The calculated interfacial energies suggest the sequence chcp Mg=bcc Nb > chcp Mg=hcp Nb > cbcc Mg=bcc Nb , which tends to favor the metastable interfacial structures bcc Mg/bcc Nb and hcp Mg/hcp Nb at small system sizes. These results support the experimental observation of metastable bcc Mg and hcp Nb phases in Mg/Nb multilayers. Moreover, the relaxed interlayer distances of surface and interface models were measured by comparing them with interlayer spacing of bulk forms. For an isolated Nb film, the electronic properties dominate on the contraction of an interlayer spacing between its surface and subsurface layers. In case of Mg film, the interlayer expansions arise from the lattice mismatch strain according to in-plane constraint. The interface structures considered in this study reveal a similar behavior to surface models, and the interface effect also dominates at interface region. The expansion values of interfacial distance tend to decrease when the number of Mg layers increases. Furthermore, the electronic charge properties were used as a tool for indicating bonding and interactions in structures. The metallic characteristics are present in both Mg and Nb regions. Comparison with electronic charge properties in bulk forms, more states are involved in valence bands in all slab models. The d orbitals provide the most contributions on interactions among Nb atoms. The partial covalent feature in metallic bonding presented in Mg region is noticed by DOS and ELF analysis. Acknowledgments Most of the calculations were carried out on the CAT cluster of the department of Chemical Engineering of Texas A&M University, the HYDRA and EOS clusters of Texas A&M Supercomputing facility, and the Texas Advanced Computing Center at the University of Texas at Austin. This research was supported by the United States National Science Foundation under NSF Grant Nos. CBET-0932249 and CMMI-0953984. References [1] L. Schapbach, A. Züttel, Nature 414 (2001) 353–358. [2] A. Krozer, B. Kasemo, J. Less-Common Met. 160 (2) (1990) 323–342. [3] K. Higuchi, H. Kajioka, K. Toiyama, H. Fujii, S. Orimo, Y. Kikuchi, J. Alloys Compd. 293–295 (0) (1999) 484–489. [4] G. Liang, J. Huot, S. Boily, A.V. Neste, R. Schulz, J. Alloys Compd. 292 (1–2) (1999) 247–252. [5] K. Higuchi, K. Yamamoto, H. Kajioka, K. Toiyama, M. Honda, S. Orimo, et al., J. Alloys Compd. 330–332 (0) (2002) 526–530. [6] E. Callini, L. Pasquini, L.H. Rude, T.K. Nielsen, T.R. Jensen, E. Bonetti, J. Appl. Phys. 108 (7) (2010) 073513. [7] N. Hanada, T. Ichikawa, H. Fujii, J. Phys. Chem. B 109 (15) (2005) 7188–7194. [8] M. Tsuda, W. Dino, H. Kasai, H. Nakanishi, H. Aikawa, Thin Solid Films 509 (1– 2) (2006) 157–159. [9] A. Baldi, G.K. Pálsson, M. Gonzalez-Silveira, H. Schreuders, M. Slaman, J.H. Rector, et al., Phys. Rev. B 81 (2010) 224203. [10] J. Qu, B. Sun, R. Yang, W. Zhao, Y. Wang, X. Li, Scripta Mater. 62 (5) (2010) 317– 320. [11] S. Singh, S.W.H. Eijt, M.W. Zandbergen, W.J. Legerstee, V.L. Svetchnikov, J. Alloys Compd. 441 (1–2) (2007) 344351. [12] H. Akyilidz, M. Özenbasß, T. Öztürk, Int. J. Hydrogen Energy 31 (10) (2006) 1379–1383. [13] B. Zahiri, C.T. Harrower, B.S. Amirkhiz, D. Mitlin, Appl. Phys. Lett. 95 (10) (2009) 103114. [14] S.E. Stoltz, D. Popovic´, Surf. Sci. 601 (6) (2007) 1507–1512. [15] L.Z. Ouyang, H. Wang, C.Y. Chung, J.H. Ahn, M. Zhu, J. Alloys Compd. 422 (1–2) (2006) 58–61. [16] J.F. Pelletier, J. Huot, M. Sutton, R. Schulz, A.R. Sandy, L.B. Lurio, et al., Phys. Rev. B 63 (2001) 052103.

225

[17] J.F.R. de Castro, S.F. Santos, A.L.M. Costa, A.R. Yavari, F.W.J. Botta, T.T. Ishikawa, J. Alloys Compd. 376 (1–2) (2004) 251–256. [18] P. Mosaner, N. Bazzanella, M. Bonelli, R. Checchetto, A. Miotello, Mater. Sci. Eng. B 108 (1–2) (2004) 33–37. [19] B. Ham, A. Junkaew, R. Arroyave, J. Chen, H. Wang, P. Wang, et al., Int. J. Hydrogen Energy 38 (20) (2013) 8328–8341. [20] B. Ham, A. Junkaew, R. Arroyave, J. Park, H.C. Zhou, D. Foley, et al., Int. J. Hydrogen Energy 39 (6) (2014) 2597–2607. [21] A. Junkaew, B. Ham, X. Zhang, A. Talapatra, R. Arróyave, Mater. Res. Lett. 1 (2013) 161–167. [22] A. Junkaew, B. Ham, X. Zhang, R. Arróyave, Calphad 45 (2014) 145–150. [23] J.Q. Zheng, J.B. Ketterson, G.P. Felcher, J. Appl. Phys. 53 (5) (1982) 3624–3628. [24] K. Ounadjela, P. Vennegues, Y. Henry, A. Michel, V. Pierron-Bohnes, J. Arabski, Phys. Rev. B 49 (1994) 8561–8573. [25] P. Boher, F. Giron, P. Houdy, P. Beauvillain, C. Chappert, P. Veillet, J. Appl. Phys. 70 (10) (1991) 5507–5511. [26] F.J. Lamelas, C.H. Lee, H. He, W. Vavra, R. Clarke, Phys. Rev. B 40 (1989) 5837– 5840. [27] E.S. Machlin, An Introduction to Aspects of Thermodynamics and Kinetics Relevant to Materials Science, third ed., Elsevier Science Ltd., Oxford, 2007. [28] R. Bruinsma, A. Zangwill, J. Phys. (J. Phys. France) 47 (12) (1986) 2055–2073. [29] A.C. Redfield, A.M. Zangwill, Phys. Rev. B 34 (1986) 1378–1380. [30] B. Ham, X. Zhang, Mater. Sci. Eng. A 528 (4–5) (2011) 2028–2033. [31] D.F. Johnson, D.E. Jiang, E.A. Carter, Surf. Sci. 601 (3) (2007) 699–705. [32] H.R. Gong, Y. Nishi, K. Cho, Appl. Phys. Lett. 91 (24) (2007) 242105. [33] W. Liu, J.C. Li, W.T. Zheng, Q. Jiang, Phys. Rev. B 73 (2006) 205421. [34] A. Arya, E.A. Carter, Surf. Sci. 560 (2004) 103–120. [35] M. Christensen, S. Dudiy, G. Wahnström, Phys. Rev. B 65 (2002) 045408. [36] S. Wang, H. Ye, Curr. Opin. Solid State Mater. Sci. 10 (1) (2006) 26–32. [37] S.J. Lee, Y.K. Lee, A. Soon, Appl. Surf. Sci. 258 (24) (2012) 9977–9981. [38] G. Kresse, J. Hafner, Phys. Rev. B 47 (1993) 558–561. [39] G. Kresse, J. Hafner, Phys. Rev. B 49 (1994) 14251. [40] G. Kresse, J. Furthmüller, Comput. Mater. Sci. 6 (1) (1996) 15–50. [41] G. Kresse, J. Furthmüller, Phys. Rev. B 54 (1996) 11169–11186. [42] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758–1775. [43] J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048–5079. [44] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. [45] K. Momma, F. Izumi, J. Appl. Crystallogr 41 (3) (2008) 653–658. [46] H. Olijnyk, W.B. Holzapfel, Phys. Rev. B 31 (1985) 4682–4683. [47] A.W. Hull, Phys. Rev. 10 (1917) 661–696. [48] E. Wachowicz, A. Kiejna, J. Phys.: Condens Matter 13 (48) (2001) 10767. [49] W.F. Gale, T.C. Totemeier, Smithells Metals Reference Book, in: W.F. Gale, T.C. Totemeier (Eds.), eighth ed., Elsevier/Butterworth-Heinemann, Amsterdam, 2004 (u.a.). [50] A.D. Becke, K.E. Edgecombe, J. Chem. Phys. 92 (9) (1990) 5397–5403. [51] P. Ravindran, P. Vajeeston, R. Vidya, H. Fjellvag, A. Kjekshus, J. Power Sources 159 (1) (2006) 88–99. [52] A. Savin, J. Mol. Struct. THEOCHEM 727 (1–3) (2005) 127–131. [53] J.L. Han, L.Z. Sun, X.D. Qu, Y.P. Chen, J.X. Zhong, Physica B 404 (1) (2009) 131– 137. [54] M. Kohout, F.R. Wagner, Y. Grin, Theor. Chem. Acc. 108 (2002) 150–156. [55] A.K. McMahan, J.A. Moriarty, Phys. Rev. B 27 (1983) 3235–3251. [56] J.A. Moriarty, A.K. McMahan, Phys. Rev. Lett. 48 (1982) 809–812. [57] N. Novakovic´, L. Matovic´, J.G. Novakovic´, M. Manasijevic´, N. Ivanovic´, Mat. Sci. Eng. B 165 (3) (2009) 235–238. [58] A.R. Jani, N.E. Brener, J. Callaway, Phys. Rev. B 38 (1988) 9425–9433. [59] H. van Leuken, A. Lodder, R.A. de Groot, J. Phys.: Condens Matter 3 (22) (1991) 3945. [60] K. Shein, I. Shein, N. Medvedeva, E. Shalaeva, M. Kuznetsov, A. Ivanovskii, Phys. Met. Metallogr 102 (6) (2006) 604–610. [61] H.R. Gong, K. Cho, Microelectron. Eng. 86 (3) (2009) 240–243. [62] B.Q. Fu, W. Liu, Z.L. Li, Appl. Surf. Sci. 255 (20) (2009) 8511–8519. [63] M. Methfessel, D. Hennig, M. Scheffler, Phys. Rev. B 46 (1992) 4816–4829. [64] M. Weinert, R.E. Watson, J.W. Davenport, G.W. Fernando, Phys. Rev. B 39 (1989) 12585–12597. [65] L. Vitos, A.V. Ruban, H.L. Skriver, J. Kollár, Surf. Sci. 411 (1–2) (1998) 186–202. [66] J.M. Zhang, D.D. Wang, K.W. Xu, Appl. Surf. Sci. 252 (23) (2006) 8217–8222. [67] W.R. Tyson, W.A. Miller, Surf. Sci. 62 (1) (1977) 267–276. [68] P.T. Sprunger, K. Pohl, H.L. Davis, E.W. Plummer, Surf. Sci. 297 (1) (1993) L48– L54. [69] P.H. Ismail, A.P. Baddorf, E.W. Plummer, Phys. Rev. B 66 (2002) 245414. [70] J.L.F.D. Silva, C. Stampfl, M. Scheffler, Surf. Sci. 600 (3) (2006) 703–715. [71] J.S. Luo, B. Legrand, Phys. Rev. B 38 (1988) 1728–1733. [72] R.M. Wentzcovitch, M.L. Cohen, Phys. Rev. B 37 (1988) 5571–5576. [73] A.F. Wright, P.J. Feibelman, S.R. Atlas, Surf. Sci. 302 (1–2) (1994) 215–222. [74] P. Staikov, T.S. Rahman, Phys. Rev. B 60 (1999) 15613–15616. [75] X.G. Li, P. Zhang, C.K. Chan, Physica B 390 (1–2) (2007) 225–230. [76] W.K. Liu, M.B. Santos, Thin Films: Heteroepitaxial Systems, Series on Directions in Condensed Matter Physics, vol. 15, World Scientific Publ., 1999 (ISBN 9789810233907). [77] C.L. Fu, S. Ohnishi, E. Wimmer, A.J. Freeman, Phys. Rev. Lett. 53 (1984) 675– 678.