Cite this Article
Li Xianglong, Wu Ping, Yang Ruijie, Yan Dan, Chen Sen, Zhang Shiping, Chen Ning. Boron diffusion in bcc-Fe studied by first-principles calculations. Chinese Physics B, 2016, 25(3): 036601  
Permissions
Boron diffusion in bcc-Fe studied by first-principles calculations
Li Xianglong1, Wu Ping1, †, , Yang Ruijie1, Yan Dan1, Chen Sen1, Zhang Shiping1, Chen Ning2
Department of Physics, School of Mathematics and Physics, University of Science and Technology Beijing, Beijing 100083, China
Department of Inorganic Non-metallic Materials, School of Materials Science and Engineering, University of Science and Technology Beijing, Beijing 100083, China

 

† Corresponding author. E-mail: pingwu@sas.ustb.edu.cn

Project supported by the National Natural Science Foundation of China (Grant No. 51276016) and the National Basic Research Program of China (Grant No. 2012CB720406).

Abstract
Abstract

The diffusion mechanism of boron in bcc-Fe has been studied by first-principles calculations. The diffusion coefficients of the interstitial mechanism, the B–monovacancy complex mechanism, and the B–divacancy complex mechanism have been calculated. The calculated diffusion coefficient of the interstitial mechanism is D0 = 1.05 × 10−7 exp (−0.75 eV/kT) m2 · s−1, while the diffusion coefficients of the B–monovacancy and the B–divacancy complex mechanisms are D1 = 1.22 × 10−6 f1 exp (−2.27 eV/kT) m2 · s−1 and D2 ≈ 8.36 × 10−6 exp (−4.81 eV/kT) m2 · s−1, respectively. The results indicate that the dominant diffusion mechanism in bcc-Fe is the interstitial mechanism through an octahedral interstitial site instead of the complex mechanism. The calculated diffusion coefficient is in accordance with the reported experiment results measured in Fe–3%Si–B alloy (bcc structure). Since the non-equilibrium segregation of boron is based on the diffusion of the complexes as suggested by the theory, our calculation reasonably explains why the non-equilibrium segregation of boron is not observed in bcc-Fe in experiments.

PACS: 66.30.J–;61.72.J–;61.72.S–;75.50.Bb;
Keyword:boron;diffusion coefficient;divacancy;first-principles calculation;
1. Introduction

It has been reported that the addition of trace amounts of boron can significantly improve the strength, hot ductility, and hardenability of ultra-low carbon microalloyed steel.[1,2] The beneficial effects may be related to the boron diffusion at high temperature and the boron segregation during the cooling process. Due to the experiment difficulties in detecting trace amounts of B (light element), the mechanism of the diffusion and segregation of boron still need to be investigated.

The non-equilibrium segregation of boron has drawn many researchers’ attention and a theory has been proposed.[3,4] In the theory, it is assumed that, in the grains, there is a certain concentration of boron–vacancy complexes. When the sample is quenched at high temperature, the supersaturated vacancies and complexes may diffuse to the grain boundaries (sink of vacancy). After the annihilation of the vacancies, the surplus boron atom will stay at the boundaries and the segregation occurs. According to Wu et al.,[5] the non-equilibrium segregation of boron is mainly observed in the steel cooled down from the fcc structure. Zhang et al.[6] did not detect non-equilibrium segregation of boron in bcc Fe–3%Si–B alloy by the particle-tracking autoradiograph (PTA) method. Therefore, it may be meaningful to study the diffusion mechanism of boron in the bcc structure.

Generally, a solute atom at a substitutional site may diffuse with a solute–vacancy mechanism, while that in an interstitial site may diffuse with an interstitial mechanism.[7] The experiments of Busby et al.[8] and Wang et al.[9] showed that boron tends to occupy a substitutional site, however, the calculations of Bialon et al.[10] and Fors et al.[11] showed that boron diffuses with an interstitial mechanism. Since the viewpoints on the diffusion mechanism are inconsistent with the solution type, a more detailed study on the diffusion mechanism of boron is needed.

Divacancies are more mobile than monovacancies, and can be more effective diffusion vehicles than monovacancies.[7] In the study of the self-diffusion of sodium (bcc structure) at high temperature, Ho[12] found that the divacancy mediated diffusion plays an important role in the diffusion of sodium. To clarify the diffusion mechanism of boron, both the contributions of the monovacancy and the divacancy should be taken into consideration.

Recently, the first-principles calculations based on the density functional theory (DFT) have been widely applied to study the interactions and the complex of point defects in crystals.[13] Taiquan[14] calculated the interaction of a Ge atom with the vacancy in silicon single crystal and found the stable complexes. Liu et al.[15] investigated the diffusion behaviors of H isotopes in W, and calculated the diffusion factor of H, which is in agreement with the experimental value. Zhou et al.[16] studied the carbon trapping mechanism in copper, and found that a monovacancy is capable of trapping as many as four C atoms to form a CnV complex. Domain et al.[17] calculated the interactions between a C or N atom and a vacancy or other interstitials in α-Fe, and found that a vacancy can trap up to two C atoms. Ohnuma et al.[18] studied the interaction between solute atoms (C/N etc.) and the vacancy in bcc Fe, and obtained the binding energy and the stable configuration of the solute–vacancy complex.

In this work, the diffusion mechanism of boron in bcc-Fe is studied by first-principles calculations. The formation energies of boron at different sites and the binding energies between the point defects are calculated in detail. Then, the diffusion coefficients of the interstitial mechanism and the complex mechanism are calculated and the diffusion mechanism is analyzed.

2. Calculation details

All the calculations were performed by using the CASTEP code.[19] The electronic exchange–correlation terms were described with the generalized gradient approximation (GGA) proposed by Perdew et al. (PW91).[20] Ultrasoft pseudo-potentials for all the elements were used. The cutoff energy of the atomic wave functions was fixed at 450 eV for all calculations to achieve a good convergence. The calculations were performed using 3×3×3 supercells which are periodic in all three directions. The energy integration over a Brillouin zone was performed with 4×4×4 k-points grids according to Monkhorst–Pack.[21] The spin-polarized approximation was adopted and a ferromagnetic model was used.

2.1. Formation and binding energies of the defects

The formation energy of single point defect i at a substitutional (S) site is calculated as

Here, the point defect i is a boron atom or a monovacancy, N is the number of atoms in the perfect supercell, Ei,Fe is the energy of the supercell with a point defect i, and is the energy of the perfect supercell. When the point defect is a boron atom, is the energy of an isolated boron atom, and when the point defect is a monovacancy, is equal to 0.

The formation energy of an interstitial boron atom is calculated as

Here, EB,Fe is the energy of the supercell with an interstitial boron atom and is the energy of an isolated boron atom.

The binding energy of two defects Eb (i, j) is defined as

However, the defects i and j may interact with each other in the supercell in this calculation since the 3×3×3 supercell is not large enough. Therefore, the binding energy of two defects Eb (i, j) is calculated by

where Ei,j, Fe is the energy of the supercell with defects i and j.[18]

The configurations with the defects we are concerned with are shown in Fig. 1.

Fig. 1. The configurations used in this work. (a) Boron atom at the substitutional site (S), the tetrahedral interstitial site (T), and the octahedral interstitial site (O). (b) The configurations of a boron atom and a monovacancy. Around the monovacancy (V), the first four (1, 2, 3, and 4) nearest S, O, T sites of the boron atom are shown. (c) The configurations of a divacancy and an adjacent boron atom near the divacancy. The first five nearest neighboring divacancies (1N, 2N, 3N, 4N, 5N) with boron atoms at their first three nearest substitutional sites (1, 2, and 3) are shown.
2.2. Diffusion coefficients
2.2.1. Diffusion coefficient of the interstitial mechanism

The diffusion coefficient of boron diffusion through octahedral interstitial sites is

where d is the jumping distance, Z is the coordinate number, ν0 is the jumping frequency from an octahedral site to its first nearest neighboring octahedral site, k is the Boltzmann constant, and Em is the jumping barrier, which can be calculated as shown in Fig. 2.

Fig. 2. (a) The jumping atom was placed at the evenly spaced sites along the jumping path in sequence. Then, the supercell was fully optimized with the coordinates of the jumping atom fixed, and the energy of the supercell was calculated. (b) Finally, the jumping barrier was obtained from the results.
2.2.2. Diffusion coefficients of the complex mechanisms

For the complex mechanisms, the solute atom diffuses through a series of exchanges with neighbor vacancies, and the diffusion coefficient is related to the exchange frequencies, the correlation factor, and the concentration of vacancies. The correlation factor is a quantity introduced so that the successive jumps (exchanges) are correlated, and a larger correlation factor means a lower probability to jump backward. In this work, two mechanisms were taken into consideration, which are the boron–monovacancy complex mechanism and the boron–divacancy complex mechanism.

2.2.2.1. Diffusion coefficient of the boron–monovacancy complex mechanism

The diffusion coefficient of the boron–monovacancy complex mechanism D1 is[7]

Here, f1 is the correlation factor, a is the lattice parameter, ω2 is the solute–vacancy exchange frequency, as shown in Fig. 3, C1V is the concentration of monovacancy, and EB−V is the binding energy of the boron atom and the monovacancy.

Fig. 3. The four-frequency model of boron diffusion in bcc-Fe. Here, ω1 is the solvent-vacancy exchange frequency, ω2 is the solute–vacancy exchange frequency, ω3 is the dissociation frequency of the complex, and ω4 is the association frequency of the complex.

According to the four-frequency model of Manning,[22] the correlation factor f1 can be calculated as

with

where ω1 is the solvent–vacancy exchange frequency, ω3 is the dissociation frequency of the complex, while ω4 is the association frequency of the complex. Since the first nearest neighbor sites of the vacancy are taken into consideration, as shown in Fig. 3, the dissociation frequency and the association frequency are replaced with the following effective frequencies:[23]

The frequencies can be calculated as

Here, νi is the jumping frequency. For the solution atom, we supposed that ν1 = ν3 = ν4 and the jumping frequency of the solute atom can be estimated by

Here, m2 and m1 are the masses of the solute atom and the solution atom; and T2 and T1 are the melting points of the solute and the solution, respectively.[24] The Ei is the jumping barrier and can be calculated with the method shown in Fig. 2.

The concentration of the monovacancy C1V in Eq. (6) can be calculated as

Here, EV is the formation energy of the monovacancy.

2.2.2.2. Diffusion coefficient of the boron–divacancy complex mechanism

The diffusion mechanism of the divacancy in the bcc structure was proposed by Mehrer.[25] The migration of a divacancy involves three configurations (1N, 2N, 4N), which are shown in Fig. 4.

Fig. 4. The divacancy configurations involved in the migration and the exchange frequencies between divacancies: , .

The diffusion coefficient of the boron–divacancy complex mechanism can be calculated with the method proposed by Mehrer[25]

where f2 is the correlation factor of the boron–divacancy complex mechanism, is the concentration of the 2N divacancies, ω21 is the exchange frequency between the 2N and 1N divacancies, and ω24 is the exchange frequency between the 2N and 4N divacancies, as shown in Fig. 4.

The concentration of the divacancy can be calculated as

where EV−V is the binding energy of the 2N divacancy.

The ω21 and ω24 can be calculated by the method introduced in Eq. (11) with v21 = v24 = v2.

The correlation factor f2 can be calculated with the method proposed by Belova et al.[26]

where s = ω21/ω24 and s < 1.

3. Results and discussion
3.1. B and monovacancy

The formation energies of the point defects were calculated. A boron atom at a substitutional site (S), an octahedral site (O), and a tetrahedral site (T) was introduced into the supercell respectively. The configurations are shown in Fig. 1(a). The formation energy of the boron atom and the monovacancy was calculated using Eq. (1) and the results are shown in Table 1.

Table 1.

The formation energy of a monovacancy (V) and a boron atom at a substitutional site (S), an octahedral interstitial site (O), and a tetrahedral interstitial site (T).

.

The results show that the formation energy of the vacancy is higher than that obtained by Domain et al.[27] The formation energy of the boron atom at O site (EO) is higher than that at S site (ES), while the formation energy of the boron atom at T site (ET) is much higher, and these are in accordance with that obtained by Fors et al.[11] Hence, we conclude that the boron atom tends to occupy the S or O site instead of the T site.

The binding energies of a boron atom at different sites (S, O, and T) with its neighbor monovacancy (V) have also been calculated. The configurations are shown in Fig. 1(b), and the results are shown in Table 2.

Table 2.

The binding energy (in eV) of a boron atom (at S/O/T site) with a monovacancy at the first (1), the second (2), the third (3), and the fourth (4) nearest neighbor sites.

.

The binding energy of a boron atom at the S site with the monovacancy at the first and the second nearest neighbor sites (S–V1 and S–V2) are 0.25 eV and 0.21 eV, and they are higher than the others. The calculation of Bialon et al.[10] also showed that the S–V1 and S–V2 configurations are more stable.

In Table 2, the configurations without binding energy presented are unstable, which means that the configurations have changed after structure optimization. In the configurations of O–V1, O–V2, and T–V1, the interstitial boron atom enters the monovacancy site and changes into a substitutional atom. For the T–V2, T–V3, and T–V4 configurations, the boron atom enters the octahedral site and the configurations change into an O–V3 configuration, while the boron atom of the O–V4 configuration transfers to a substitution site by pushing an adjacent Fe atom into the monovacancy site.

It seems that the configuration becomes unstable when there is a monovacancy near the interstitial boron atom. Bialon et al.[10] also found that the configurations of O–V1 and O–V2 are unstable.

3.2. B and divacancy

The binding energies of the first five nearest neighboring divacancies were calculated and the configurations are shown in Fig. 1(c). The binding energies of the divacancies and the comparisons with the references are shown in Table 3.

Table 3.

The binding energies (in eV) of the divacancies.

.

The results show that the binding energies of the 1N, 2N, and 4N nearest neighboring divacancies are larger than the others, and the binding energy of the 2N divacancy is the largest, which is in accordance with the results in the references.

The interaction between the boron atom and the neighboring divacancy was studied. As the configuration becomes unstable when there is a vacancy near the interstitial boron atom, only the boron atom at a substitutional site was considered. The binding energies between the five types of divacancies (1N, 2N, 3N, 4N, 5N) and the substitutional boron atom at their first three nearest neighboring sites (1nn, 2nn, 3nn) were calculated. The configurations used in the calculations are shown in Fig. 1(c) and the binding energies are shown in Table 4.

Table 4.

The binding energy between the substitutional atom and its neighboring divacancy.

.

The results show that there are four complexes of the substitutional boron atom and the divacancy, whose binding energies are higher. They are the 1N, 2N, and 4N divacancies with the boron atom at the first nearest site (1N–1nn, 2N–1nn, 4N–1nn), and the 3N divacancy with the boron atom at the second nearest site (3N–2nn).

According to Mehrer,[25] there are three different configurations involved in the migration of a divacancy in the bcc structure, which are the 1N, 2N, and 4N divacancies.

In our results, the binding energies of the 1N/2N/4N divacancies and the 1nn boron atom are larger than most of the others, which means that they are likely to form complexes. Therefore, it is necessary to calculate the effect of the divacancy in the diffusion of the boron atom.

3.3. Diffusion coefficient of B

In this section, the diffusion coefficients of both the interstitial mechanism and the complex mechanism are studied.

3.3.1. Diffusion coefficient of the interstitial mechanism

The diffusion coefficient can be calculated with Eq. (5) and the parameters needed are listed in Table 5. Here, we supposed that ν0 is 8.0 THz.[11]

Table 5.

The parameters for the calculation of the diffusion coefficient.

.
3.3.2. Diffusion coefficient of the boron–monovacancy complex mechanism

As shown in Fig. 3, the diffusion coefficient of the boron–monovacancy mechanism has been calculated in the frame of the four-frequency mechanism with ω3 and ω4 replaced by the efficient frequencies.

The jumping barrier energy and the jumping frequency needed for the calculation of the diffusion coefficient are listed in Table 6. The jumping frequency of the solution atom has been set to be 5.7 THz,[11] while that of the solute atom can be calculated with Eq. (12).

Table 6.

The parameters needed for the calculation of the diffusion coefficient.

.

The correlation factor f1 can be calculated with the methods introduced before and the results are shown in Table 7.

Table 7.

The correlation factors of the boron–monovacancy mechanism (f1).

.

The correlation factors of the boron–monovacancy diffusion mechanism at different temperatures are shown in Table 7. It is easily found that ω2 is much larger than ω3 from Eq. (11) with the parameters shown in Table 5. Therefore, it is reasonable to get the tiny correlation factor from Eq. (7).

The diffusion coefficient of the boron–monovacancy complex mechanism can be obtained from Eq. (6) with the parameters given in Tables 6 and 7.

3.3.3. Diffusion coefficient of the boron–divacancy complex mechanism

The parameters for the calculation of the diffusion coefficient of the boron–divacancy complex are shown in Table 8.

Table 8.

The parameters needed for the calculation of the diffusion coefficient.

.

The correlation factor of the boron–divacancy mechanism has been calculated with Eq. (16) and the results are shown in Table 9. The correlation factors are approximately equal to 0.43 and change little in the temperature interval.

Table 9.

The correlation factors of the boron–divacancy mechanism (f2).

.

Thus, the diffusion coefficient of the boron–divacancy complex mechanism can be deduced from Eq. (14) with the parameters given in Tables 8 and 9.

The diffusion coefficients of the interstitial mechanism and the boron–vacancy mechanism are shown in Fig. 5. It is clear that the diffusion coefficient of the complex mechanism is far less than that of the interstitial mechanism. This may be caused by the large vacancy formation energy and the small binding energy between the boron atom and the vacancy or divacancy. Therefore, the dominating diffusion mechanism of boron in bcc-Fe is the interstitial mechanism.

Fig. 5. The calculated diffusion coefficients of different diffusion mechanisms and the results from experiments.

Our calculations together with the calculation of Bialon et al.[25] show that the interstitial solute atom may enter a substitutional site when there is a vacancy near it. However, as mentioned above, the high formation energy of a vacancy makes the concentration of the vacancies too low to have an effect on the interstitial diffusion mechanism. Therefore, the dominant diffusion mechanism in boron is still the interstitial diffusion mechanism.

The experiment results of Busby et al.[8] and Wang et al.[9] are shown in Fig. 5. Our results show good agreement with the experiment results of Wang et al.,[9] and are in the interval of the results obtained by Busby et al.[8] Here, it should be noted that the magnetic structure of iron changes at its Curie point (1043 K). As shown in Fig. 5, the diffusion coefficients from the experiments do not change a lot around 1043 K, thus we did not take the magnetic structure change into consideration. Therefore, we may conclude that B diffuses through the interstitial mechanism in bcc-Fe and the diffusion coefficient is D = 8.04 × 10−8 × exp(−0.75 eV/kT) m2·s−1.

Our results show that the dominant diffusion mechanism is the interstitial mechanism. According to the theory of non-equilibrium segregation, the non-equilibrium segregation cannot form without the diffusion of the boron–vacancy complex. Zhang et al.[6] have studied the segregation of boron in bcc-Fe by the PTA method, no observable non-equilibrium segregation was detected in bcc-Fe, and this agrees with our conclusion.

4. Conclusion

The boron diffusion coefficients of the interstitial mechanism and the complex mechanism in bcc-Fe have been obtained by first-principles calculations. In particular, for the complex mechanism, both the monovacancy and the divacancy mechanisms were considered. The results are summarized as follows.

Reference
1He XChu YKe J1983Acta Metallurgica Sinica19A459
2Naderi MKetabchi MAbbasi MBleck W 2011 Procedia Engineering 10 460
3Anthony T RHanneman R E 1968 Scripta Metallurgica 2 611
4Aust K THanneman R ENiessen PWestbrook J H 1968 Acta Metallurgica 16 291
5Wu PHe X L1999Acta Metallurgica Sinica351009
6Sanhong ZXinlai HYouyi CJun K1993Acta Metallurgica Sinica29A241
7Mehrer H 2007 Diffusion in Solids: Fundamentals, Methods, Materials, Diffusion-controlled Processes 94
8Busby P EWells C1954Journal of Metals6972
9Wang WZhang SHe X 1995 Acta Metallurgica et Materialia 43 1693
10Bialon AHammerschmidt TDrautz R 2013 Phys. Rev. B 87 104109
11Fors D H RWahnstrom G 2008 Phys. Rev. B 77 132102
12Ho P S 1972 Philosophical Magazine 26 1429
13Becquart C SDomain C 2012 Current Opinion in Solid State and Materials Science 16 115
14Tai Q W 2011 Acta Phys. Sin. 61 063101 (in Chinese)
15Liu Y LLu WGao A YGui L JZhang Y 2012 Chin. Phys. B 21 126103
16Zhou H BJin S 2013 Chin. Phys. B 22 076104
17Domain CBecquart CFoct J 2004 Phys. Rev. B 69 144112
18Ohnuma TSoneda NIwasawa M 2009 Acta Materialia 57 5947
19Segall M DD P JLindanProbert M JPickard C JHasnip P JClark S JPayne M C 2002 J. Phys.: Condens. Matter 14 2717
20Perdew J PJackson K APederson M RSingh D JFiolhais C 1992 Phys. Rev. B 46 6671
21Monkhorst H JPack J D 1976 Phys. Rev. B 13 5188
22Manning J R 1964 Phys. Rev. 136 A1758
23Choudhury SBarnard LTucker J DAllen T RWirth B DAsta MMorgan D 2011 J. Nucl. Mater 411 1
24Qiong WShu S LYue MSheng K G 2012 Chin. Phys. B 21 109102
25Mehrer H 1973 J. Phys. F: Metal Phys. 3 543
26Belova I VGentle D SMurch G E 2002 Philosophical Magazine Letters 82 37
27Domain CBecquart C 2001 Phys. Rev. B 65 024103
28Johnson R 1964 Phys. Rev. 134 A1329