Volumetric Strain Dependence of Quantum Diffusion of Hydrogen in bcc Iron

Abstract

To predict hydrogen embrittlement in steels and clarify its mechanism, it is necessary to understand the time variation of hydrogen distribution. Since the quantum mechanical effect is remarkably observed in hydrogen diffusion, even at room temperature, there is a need to take into account it for analyses. In this study, we evaluated the diffusion coefficients of hydrogen in body-centered cubic (bcc) iron via density functional theory and small-polaron theory calculations. The analyses were carried out under various magnitudes of volumetric strains to investigate its effect on the diffusion coefficient. The temperature dependence of the diffusion coefficient was found to change at about 400 K. This is attributed to the fact that the tunneling between ground states dominantly contributes to the diffusion at lower temperatures, whereas at high temperatures, that between low excited states contributes dominantly. The diffusion coefficient was also found to increase with compressive volumetric strain and decrease with tensile volumetric strain. The volumetric strain dependence of the diffusion coefficient was clarified by the volumetric strain dependence of the tunneling matrix elements and that of the activation enthalpy.

1. Introduction

To predict hydrogen embrittlement in steels and to understand its mechanism, the time variation of hydrogen distribution, i.e., the diffusion behavior of hydrogen, needs to be studied. It has been found that in iron, the diffusion behavior changes remarkably at about 500 K due to the quantum mechanical effect, which emerges at lower temperatures.¹⁾ Thus, for iron, though the classical molecular dynamics (MD) calculation is effective in analyzing the phenomenon at higher temperatures, analyses that adopt the quantum mechanical effect are required at lower temperatures, including room temperature. As a hydrogen atom has the lightest mass, among all atoms, quantum mechanical effects, such as tunneling effect, remarkably affect its diffusion behavior in solids.²⁾ The quantum mechanical effect tends to strongly affect hydrogen diffusion in body-centered cubic (bcc) metals due to the short jump distance. The effect is remarkable in iron, which has the smallest lattice constant among the bcc transition metals.

In the prediction and evaluation of hydrogen diffusion in metals, a diffusion equation that incorporates a term dependent on the gradient of isotropic stress is employed.^3,4,5) On the other hand, it is also important to clarify the stress or strain dependence of the diffusion coefficient. Taketomi and Matsumoto carried out an analysis via the MD method, and found that the hydrogen diffusion behavior in bcc iron exhibits strong anisotropy around the edge dislocation. They attributed the behavior to the contribution of anisotropic stress.⁶⁾ Nagase and Matsumoto also evaluated the stress dependence of the diffusion coefficient of hydrogen in bcc iron via the MD method.^7,8) However, the effect of stress on the quantum diffusion of hydrogen is yet to be clarified.

The small-polaron theory^9,10) by Holstein, which was basically developed to describe the tunneling process of electrons in semiconductors and insulators, was extended to studying the diffusion of hydrogen in metals by Flynn and Stoneham.¹¹⁾ The hydrogen diffusion behavior of various bcc metals, such as vanadium,^12,13) niobium,^{12,14,15,16,17)} and tantalum,^12,15,16,17) has been analyzed via the small-polaron theory. In addition, the small-polaron theory has been employed in the analysis of iron;¹⁸⁾ however, the calculations considering higher-order excited states and the evaluation under stress or strain are yet to be carried out.

This study aims at clarifying the quantum diffusion behavior of hydrogen in bcc iron under volumetric strain via density functional theory (DFT) and small-polaron theory calculations.

2. Analysis Methods and Calculation Conditions

2.1. Self-trapped and Coincidence States

Hydrogen in bcc iron occupies a tetrahedral interstitial site (T-site).^19,20,21) The repulsive force from the hydrogen atom pushes the surrounding iron atoms outward, and thus, the potential well is lowered. This makes the hydrogen atom more stable at the T-site, and this is called self-trapped state.²⁾ For a hydrogen atom in the self-trapped state at a T-site to move through a potential wall to an adjacent T-site via tunneling, thermal activation is required to bring the energy levels of the initial and final sites of the migration to the same height.²²⁾ When the energy levels of the two sites are at the same height, they are referred to as a coincidence state. Though diffusion via tunneling can occur between excited states, further activation enthalpy is required.

2.2. Calculation Method of Diffusion Coefficient by Tunneling Process

In this paper, the following equation used by Klamt and Teichler as the jump frequency is employed to calculate the diffusion coefficient considering the quantum mechanical effect:¹⁵⁾   

$$w_{ij}=\frac{J_{ij}{}^2\pi }{h}\sqrt{\frac{\pi }{E_{ij}^{\text{a}}{k}_{\text{B}}T}}exp\left(-\frac{E_{ij}^{\text{a}}+{\epsilon }_i^0}{k_{\text{B}}T}\right)/Z$$(1)

where w_ij, J_ij, and E^a_ij indicate the jump frequency, tunneling matrix element, and activation enthalpy, respectively, when the initial state is the i-th excited state and the final state is the j-th excited state. h is the Planck constant, k_B is the Boltzmann constant, T is the temperature, ε⁰_i is the eigenenergy of i-th excited state for a hydrogen atom in the self-trapped state at the initial T-site, and Z is the partition function, expressed as   

$$Z={\displaystyle \sum }_{i}exp(-{\epsilon }_i^0/k_{\text{B}}T)$$(2)

Diffusion via tunneling requires some overlap between the wave functions of the initial and final states. The degree of overlap is quantitatively expressed by the tunneling matrix element J_ij.²⁾ Generally, J_ij changes remarkably with the distance between the adjacent sites; it increases as the distance decreases. Therefore, the quantum effect is anticipated to be remarkable in iron, which has the smallest lattice constant among the bcc transition metals.

There are two types of coincidence states: symmetric coincidence states, where the excited states of the initial and final states are the same, i.e., i = j, and asymmetric coincidence states, where they are different, i.e., i ≠ j. The contribution of the tunneling process through asymmetric coincidence states has been found to be small,^13,15) and hence, we considered only symmetric coincidence states in this study.

To evaluate w_ij using Eq. (1), ε⁰_i, J_ii, and E^a_ii are required; herein, they are determined using Yoshinari’s procedure for the vanadium–hydrogen system¹³⁾ (see also the flowchart shown in Fig. 1):

Fig. 1.

Flowchart of the calculation of diffusion coefficient.

1) The atomic configurations of iron are obtained when the hydrogen atom is in the self-trapped state at a T-site and when it is in the saddle point (SP) configuration during classical diffusion between adjacent T-sites. The former is obtained by inserting a hydrogen atom into the T-site of the iron lattice and performing structural relaxation, while the latter is obtained via the nudged elastic band (NEB) method.

2) The potential energy surface (PESs) of hydrogen for the atomic configurations of iron are calculated when the hydrogen atoms are at the T-site and at the SP. Mapping of 3D PES is carried out by assigning grid points to the rectangular parallelepiped region centered on the midpoint of the two adjacent T-sites (Fig. 2), and repeating DFT calculations where the iron atoms are fixed to the positions obtained in step 1) and the hydrogen atom is fixed to each grid point. In this study, to reduce the calculation cost, the following two contrivances are conducted. First, we consider the symmetry of the calculation region of the PES; it is symmetrical along the (001) plane containing two T-sites, as well as along a plane parallel to the xz-plane and including the initial T-site when the hydrogen atom is in the self-trapped state at the T site. It is as well-symmetrical along the (110) plane when the hydrogen is at the SP. Second, as a very high ground-state energy is obtained when the iron and hydrogen atoms are close to each other, the DFT calculation is not conducted at the grid points where the distance to the closest iron atom is 0.4 a₀ or less, where a₀ is the lattice constant of iron in the no-strain state. The ground-state energy for these configurations is set to 0 eV. With these two measures, the number of grid points for calculating the PES can be reduced to about 1/6 or 1/5.

Fig. 2.

Calculation region for PES. ●, ▼, and ■ indicate Fe atoms, T-sites, and SP, respectively. The 3D PES was calculated for the shaded region.

3) The eigenenergies of the hydrogen atom are calculated by solving the Schrödinger equation for the PES of the hydrogen at the T-site and SP configurations obtained from the spline interpolation. The Fourier grid Hamiltonian (FGH) method²³⁾ is employed to solve the Schrödinger equation. ε⁰_i is calculated from the initial configuration. Two coincidence states are obtained¹³⁾ at the SP configuration where the corresponding wave functions have the same and opposite signs at the two T sites, and those eigenenergies are ${\epsilon }_{ii}^{\text{SP}}$ and ${\epsilon }_{ii}^{\text{SP*}}$.

4) From the eigenenergies of the hydrogen atom at the T-site and at the SP, the parameters J_ii and E^a_ii are calculated, using which the jump frequency w_ii is obtained. J_ii is obtained as $\left|{\epsilon }_{ii}^{\text{SP}}-{\epsilon }_{ii}^{\text{SP*}}\right|/2$ because the energy difference between these two states is 2J_ii.^13,17,24) The mean energy of the two states is the energy of the coincidence state E^c_ii, and the difference between E^c_ii and ε⁰_i is E^a_ii, i.e., E^a_ii = E^c_ii − ε⁰_i.

5) On the basis of the random walk theory,²⁵⁾ the diffusion coefficient D is obtained from the jump frequency w_ii. The diffusion coefficient is given as   

$$D=\frac{1}{12}{\mathrm{\Gamma }}_0{a}^2$$(3)

where a is the lattice constant of iron and Γ₀ is the total jump frequency, which is obtained by summing the jump frequencies w_ii through each excited state; thus,   

$${\mathrm{\Gamma }}_0={\displaystyle \sum }_i{w}_{ii}$$(4)

2.3. Analysis Models and Conditions

In this study, the Vienna Ab-initio Simulation Package^26,27,28) is employed for the DFT calculations. In all analyses, the generalized gradient approximation exchange-correlation energy function of Perdew–Burke–Ernzerhof²⁹⁾ is employed for the approximation of the exchange-correlation energy. The projector augmented wave method³⁰⁾ is used for the treatment of the core electron, and the calculation considering the electron spin is carried out. The Monkhorst-Pack scheme³¹⁾ is used to define the k-points.

A 2 × 2 × 2 bcc lattice supercell, which contains 16 iron atoms and 1 hydrogen atom (Fe16H1), is used for all calculations. A k-point mesh of 9 × 9 × 9 is used. The cutoff energy is set to 500 eV for the structural relaxation and the NEB calculations and 475 eV for the PES calculation. These calculation conditions are set so as to obtain sufficient calculation accuracy after carrying out comparative calculations. To evaluate the volumetric strain dependence, the calculation is carried out under the conditions of 0.0, ±0.6, ±1.5, and ±3.0% volumetric strains. Under the no-strain condition, both cell size and atomic positions are relaxed in the structural relaxation and NEB calculations. For the with-strain conditions, the predetermined volumetric strain is applied to each atomic configuration obtained from the no-strain condition. During the structural relaxation and NEB calculations under strain conditions, the cell sizes are fixed, and thus, only the atomic positions are relaxed.

The number of intermediate atomic configurations in the NEB calculation is three, and the initial path is obtained from the linear interpolation between the initial and final states. The PES calculations are performed on a rectangular parallelepiped region of 0.75 a_x × 0.75 a_y × 0.8 a_z, where a_x, a_y, and a_z are the lattice constants in the x, y, and z directions, respectively, giving rise to PES grid points of 13 × 13 × 13. The FGH calculation is performed for 25 × 25 × 25 PES grid points obtained by applying the spline interpolation to 13 × 13 × 13 PES grid points.

3. Analysis Results

3.1. Results under No Strain Condition

Figure 3 shows the minimum energy path (MEP) for the classical diffusion of hydrogen atoms obtained via NEB under the no-strain condition. The result is in good agreement with that of the DFT calculation performed for the model of Fe₁₂₈H₁ by Jiang et al.²¹⁾ The SP energy, which is based on the energy of the self-trapped state, is determined as 86 meV in this analysis, which is also in good agreement with the 88 meV reported by Jiang et al.²¹⁾ Thus, it can be said that this analysis for the Fe₁₆H₁ model can appropriately reproduce the hydrogen diffusion in larger lattices. Here, the volumetric strain at the SP was obtained as 2.6 × 10⁻⁴ for the current model. We neglect its influence on the activation enthalpy in Section 3.2, because the contribution is very small, that is < 1 meV under 3% of volumetric strain condition.

Fig. 3.

Initial path and MEP obtained via NEB calculation.

Figure 4 shows the PES of hydrogen in the cross-section including the initial and final T-sites. Figure 4(a) shows the PES for the initial state. The minimum value of energy when the hydrogen atom is near the T-site is taken to be zero of PES. When the hydrogen atom is near the T-site located at (0.25, 0.5, 0.5), the energy is 50–100 meV lower than that in the case where the hydrogen is at other T-sites; that is, the hydrogen atom is in a self-trapped state at the T-site. The PES for the atomic configuration of iron when the hydrogen atom is present at the SP is shown in Fig. 4(b).

Fig. 4.

PES at the cross-section including initial and final T sites. (a) Self-trapped configuration at a T site. (b) SP configuration: The coordinates are expressed on the basis of the lattice constant in each axial direction. (Online version in color.)

The eigenenergies of the self-trapped state of the hydrogen atom obtained by solving the Schrödinger equation using the FGH method are listed in Table 1. The energy value is based on the ground state energy of the system when the hydrogen atom is at the self-trapped configuration. The wave functions of the various eigenstates are shown in the upper panel of Fig. 5. The wave functions are visualized using a 3D visualization code, VESTA.³²⁾ Here, each eigenstate is classified according to the shape of the corresponding wave function. First, the ground state is termed the g-state. The excited states are classified by the number and direction of nodal planes of the wave function. For example, the 1x state indicates one nodal plane perpendicular to the x-direction.

Table 1. Eigenenergy for the self-trapped hydrogen ε⁰ (meV).

Eigenstate	g	1x	1y	1z	2x
Eigenenergy ε0	232	327	373	385	412

Fig. 5.

Wave functions of hydrogen at the self-trapped position (upper panel) and those at coincidence state (lower panel) in bcc iron. (Online version in color.)

The wave function in the coincidence state is shown in the lower panel of Fig. 5. As mentioned in Section 2.2, two coincidence states in which the corresponding wave function has the same sign and that with opposite signs at the two T-sites for the SP configuration are obtained. For example, regarding the coincidence between the g-states, the former is expressed as g-g and the latter as (g-g)^*.

The tunneling matrix element J in each coincidence state, the energy in the coincidence state E^c, and the activation enthalpy E^a are listed in Table 2. It is found that J between the excited states (1x-1x, 1y-1y) is larger than that between the ground states (g-g). This can be attributed to the fact that the wave function spatially expands along a specific direction at an excited state, and thus, the overlap of the wave functions of the initial and final T sites of the migration increases.

Table 2. Tunneling matrix elements J, energies of the coincidence states E^c, and activation energies E^a (meV).

Coincidence state	g-g	1x-1x	1y-1y	1z-1z	2x-2x
Tunneling matrix element J	12	19	27	4	10
Coincidence State Energy Ec	264	362	396	468	450
Activation enthalpy Ea	32	35	12	94	38

The diffusion coefficient is calculated using Eq. (3). The Arrhenius plots of the contribution of each coincidence state shown by a broken line and the total diffusion coefficient shown by a thick line are shown in Fig. 6. The value obtained using the path integral centroid molecular dynamics (CMD) method¹⁾ is indicated by ▲, and that from experiments^{33,34,35,36,37,38)} is indicated by ○. In the CMD method, the interaction forces are calculated based on the interatomic potential, and the motion of hydrogen atom follows quantum mechanics. Figure 6 shows that the calculation in this study provided a slightly higher diffusion coefficient than the reference values. The difference may be due to the influences of lattice defects and impurity elements in the experiment, the accuracy of the interatomic potential in CMD, and the limitation of the model size in this study. A bend is observed in the total diffusion coefficient curve at about 400 K. This temperature is slightly lower than the transition temperature estimated via the CMD method,¹⁾ i.e., 500 K where the transition from quantum to classical diffusion occurs. It is found that the diffusion through the g-g coincidence state is dominant at low temperatures, whereas the contribution from the diffusion in the 1y-1y state is dominant at high temperatures. Here, the 2x-2x state showed a smaller contribution to the diffusion, although the values of J and E^a are close to those of the g-g state, as shown in Table 2. This is attributed to the fact that the eigenenergy ε⁰_2x is high, as shown in Table 1, and hence, additional energy is required for excitation.

Fig. 6.

Diffusion coefficients of hydrogen in bcc iron under no strain condition: ▲ and ○ indicate the CMD results¹⁾ and experimental data,^{33,34,35,36,37,38)} respectively.

3.2. Results under Volumetric Strain Conditions

The diffusion coefficients are calculated under 0.0, ±0.6, ±1.5, and ±3.0% volumetric strain conditions, as stated in the previous section. The volumetric strain dependence of the tunneling matrix element J and that of the activation enthalpy E^a are shown in Figs. 7 and 8, respectively. As shown in the figure, J increases slightly with compressive strain and decreases slightly with an increase in tensile strain. This is considered to be caused by the increase in the overlap of the wave functions by approaching two T-sites as a result of compressive strain. On the other hand, E^a tends to increase with tensile strain. The strain effect is stronger for the excited states. The jump frequency is an increasing function of J and a decreasing function of E^a. In other words, the trends of both J and E^a under volumetric strain contribute to the direction in which the jump frequency increases due to compressive strain and decreases due to tensile strain. Thus, the diffusion coefficient is also expected to increase with compressive strain and decrease with tensile strain.

Fig. 7.

Volumetric strain dependence of the tunneling matrix element J.

Fig. 8.

Volumetric strain dependency of the activation enthalpy E_a.

The Arrhenius plots of the diffusion coefficients obtained under each volumetric strain condition are shown in Fig. 9. Here, only the total diffusion coefficient and g, 1x, and 1y states, which show large contributions, are shown. As shown in the figure, the diffusion coefficient increases with compressive strain, and decreases with increasing tensile strain. This is consistent with the volumetric strain dependence of the diffusion coefficient predicted from the trends of J and E^a under the volumetric strain. The changes in the diffusion coefficient observed at 300 K under 3% compressive and tensile volumetric strain are about 23% and −14%, respectively.

Fig. 9.

Diffusion coefficients of hydrogen in bcc iron under volumetric strain. (Online version in color.)

The stress dependence of the diffusion coefficient of hydrogen in Pd was evaluated via the MD method by Nishimura.³⁹⁾ The diffusion coefficient was found to decrease with increasing compression. This demonstrates a trend opposite to that obtained in this study. While Pd has a relatively long jump distance for the hydrogen diffusion, bcc iron has the very short jump distance where the tunneling effect strongly appears. The diffusion is anticipated to be promoted in the compressive strain side of an edge dislocation in bcc iron, whereas it is promoted in the tensile strain side of the dislocation in Pd.⁴⁰⁾

4. Conclusion

In this study, we evaluated the diffusion coefficient of hydrogen in bcc iron using DFT and small-polaron theory. The effect of volumetric strain on hydrogen diffusion was also clarified.

We ascertained that in the no-strain condition, the ground state (g-state) and the low excited state (1y state) mainly contributed to the diffusion; the former dominated at low temperatures, while the latter dominated at high temperatures. A bend was observed in the Arrhenius plot of the diffusion coefficient at about 400 K, which is inferred to be caused by the change in the excited state, which is dominant in the diffusion. The contributions of the g and 1y states under volumetric strain, as well as under the no-strain condition, were high. The diffusion coefficient was found to increase with compressive volumetric strain and decrease with increasing tensile volumetric strain. The volumetric strain dependence of the diffusion coefficient was clarified by evaluating the volumetric strain dependence of the tunneling matrix element J and that of the activation enthalpy E^a.

References

• 1)   H.  Kimizuka,  H.  Mori and  S.  Ogata: Phys. Rev. B, 83 (2011), 094110. https://doi.org/10.1103/PhysRevB.83.094110
• 2)   Y.  Fukai: The Metal-Hydrogen System: Basic Bulk Properties, Springer-Verlag, Berlin, Heidelberg, (2005), 1. https://doi.org/10.1007/3-540-28883-X
• 3)   P.  Sofronis and  R. M.  McMeeking: J. Mech. Phys. Solids, 37 (1989), 317. https://doi.org/10.1016/0022-5096(89)90002-1
• 4)   H.  Kotake,  R.  Matsumoto,  S.  Taketomi and  N.  Miyazaki: Int. J. Press. Vessel. Pip., 85 (2008), 540. https://doi.org/10.1016/j.ijpvp.2008.02.002
• 5)   K.  Takayama,  R.  Matsumoto,  S.  Taketomi and  N.  Miyazaki: Int. J. Hydrog. Energy, 36 (2011), 1037. https://doi.org/10.1016/j.ijhydene.2010.10.046
• 6)   S.  Taketomi,  R.  Matsumoto and  N.  Miyazaki: Acta Mater., 56 (2008), 3761. https://doi.org/10.1016/j.actamat.2008.04.011
• 7)   S.  Nagase and  R.  Matsumoto: J. Soc. Mater. Sci., Jpn., 69 (2020), 119 (in Japanese). https://doi.org/10.2472/jsms.69.119
• 8)   S.  Nagase and  R.  Matsumoto: Mater. Trans., 61 (2020), 1265. https://doi.org/10.2320/matertrans.Z-M2020823
• 9)   T.  Holstein: Ann. Phys. (N. Y.), 8 (1959), 325. https://doi.org/10.1016/0003-4916(59)90002-8
• 10)   T.  Holstein: Ann. Phys. (N. Y.), 8 (1959), 343. https://doi.org/10.1016/0003-4916(59)90003-X
• 11)   C. P.  Flynn and  A. M.  Stoneham: Phys. Rev. B, 1 (1970), 3966. https://doi.org/10.1103/PhysRevB.1.3966
• 12)   H.  Sugimoto and  Y.  Fukai: Hyperfine Interact., 17 (1984), 139.
• 13)   O.  Yoshinari: J. Alloys Compd., 580 (2013), S36. https://doi.org/10.1016/j.jallcom.2012.12.108
• 14)   D.  Emin,  M. I.  Baskes and  W. D.  Wilson: Phys. Rev. Lett., 42 (1979), 791. https://doi.org/10.1103/PhysRevLett.42.791
• 15)   A.  Klamt and  H.  Teichler: Phys. Status Solidi, 134 (1986), 533. https://doi.org/10.1002/pssb.2221340211
• 16)   P. G.  Sundell and  G.  Wahnström: Phys. Rev. Lett., 92 (2004), 155901. https://doi.org/10.1103/PhysRevLett.92.155901
• 17)   P. G.  Sundell and  G.  Wahnström: Phys. Rev. B, 70 (2004), 224301. https://doi.org/10.1103/PhysRevB.70.224301
• 18)   D.  Di Stefano,  M.  Mrovec and  C.  Elsässer: Phys. Rev. B, 92 (2015), 224301. https://doi.org/10.1103/PhysRevB.92.224301
• 19)   J. P.  Hirth: Metall. Trans. A, 11 (1980), 861. https://doi.org/10.1007/BF02654700
• 20)   J.  Sanchez,  J.  Fullea,  C.  Andrade and  P. L.  De Andres: Phys. Rev. B, 78 (2008), 014113. https://doi.org/10.1103/PhysRevB.78.014113
• 21)   D. E.  Jiang and  E. A.  Carter: Phys. Rev. B, 70 (2004), 064102. https://doi.org/10.1103/PhysRevB.70.064102
• 22)   K. W.  Kehr: Hydrogen in Metals I: Basic Properties, eds. by G. Alefeld and J. Völkl, Springer-Verlag, Berlin, Heidelberg, (1978), 197. https://doi.org/10.1007/3540087052_47
• 23)   C. C.  Marston and  G. G.  Balint-Kurti: J. Chem. Phys., 91 (1989), 3571. https://doi.org/10.1063/1.456888
• 24)   A.  Klamt and  H.  Teichler: Phys. Status Solidi, 134 (1986), 103. https://doi.org/10.1002/pssb.2221340113
• 25)   P. G.  Shewmon: Diffusion in Solids, McGraw-Hill, New York, (1963), 41.
• 26)   G.  Kresse and  J.  Hafner: Phys. Rev. B, 48 (1993), 13115. https://doi.org/10.1103/PhysRevB.48.13115
• 27)   G.  Kresse and  J.  Furthmüller: Phys. Rev. B, 54 (1996), 11169. https://doi.org/10.1103/PhysRevB.54.11169
• 28)   G.  Kresse and  J.  Furthmüller: Comput. Mater. Sci., 6 (1996), 15. https://doi.org/10.1016/0927-0256(96)00008-0
• 29)   J. P.  Perdew,  K.  Burke and  M.  Ernzerhof: Phys. Rev. Lett., 77 (1996), 3865. https://doi.org/10.1103/PhysRevLett.77.3865
• 30)   P. E.  Blöchl: Phys. Rev. B, 50 (1994), https://doi.org/10.1103/PhysRevB.50.17953
• 31)   H. J.  Monkhorst and  J. D.  Pack: Phys. Rev. B, 13 (1976), 5188. https://doi.org/10.1103/PhysRevB.13.5188
• 32)   K.  Momma and  F.  Izumi: J. Appl. Crystallogr., 44 (2011), 1272. https://doi.org/10.1107/S0021889811038970
• 33)   W. L.  Bryan and  B. F.  Dodge: AIChE J., 9 (1963), 223. https://doi.org/10.1002/aic.690090217
• 34)   J. L.  Dillard: C. R. Hebd. Seances Acad. Sci., Ser. C (Rep. Acad. Sci.), 270 (1970), 669.
• 35)   T.  Heumann and  E.  Domke: Ber. Bunsenges. Phys. Chem. (Rep. Bunsen Soc. Phys. Chem.), 76 (1972), 825.
• 36)   M.  Nagano,  Y.  Hayashi,  N.  Ohtani,  M.  Isshiki and  K.  Igaki: Scr. Metall., 16 (1982), 973. https://doi.org/10.1016/0036-9748(82)90136-3
• 37)   Y.  Hayashi,  H.  Hagi and  A.  Tahara: Z. Phys. Chem. (J. Phys. Chem.), 164 (1989), 815. https://doi.org/10.1524/zpch.1989.164.Part_1.0815
• 38)   H.  Hagi: Mater. Trans., JIM, 35 (1994), 112. https://doi.org/10.2320/matertrans1989.35.112
• 39)   K.  Nishimura: J. Soc. Mater. Sci., Jpn., 63 (2014), 163 (in Japanese). https://doi.org/10.2472/jsms.63.163
• 40)   K.  Nishimura and  K.  Miyake: J. Soc. Mater. Sci., Jpn., 65 (2016), 148 (in Japanese). https://doi.org/10.2472/jsms.65.148