Molecular Dynamics Study of Phosphorus Migration in Σ3(111) and Σ5(0-13) Grain Boundaries of α-Iron

Phosphorus atoms in steels accumulate at grain boundaries via thermal and irradiation effects and induce grain boundary embrittlement, which is experimentally confirmed by an increase in the ductile-brittle transition temperature. Quantitative prediction of phosphorus segregation at grain boundaries under various temperature and irradiation conditions is essential for preventing embrittlement. To develop a model of grain boundary phosphorus segregation in α-iron, we studied the migration of a phosphorus atom in two types of symmetrical tilt grain boundaries (Σ3[1-10](111) and Σ5[100](0-13) grain boundaries) using molecular dynamics simulations with an embedded atom method potential. The results revealed that, in the Σ3 grain boundary, phosphorus atoms migrate three-dimensionally mainly in the form of interstitial atoms, whereas in the Σ5 grain boundary, these atoms migrate one-dimensionally mainly via vacancy-atom exchanges. Moreover, de-trapping of phosphorus atoms and vacancies was investigated.


Introduction
Previous studies have shown that phosphorus (P) reduces the cohesive energy of grain boundaries (GBs) in steels [1,2] and leads to GB embrittlement. The segregation of P atoms is promoted by the thermal and irradiation effect [3][4][5]. The resulting embrittlement is referred to as temper embrittlement and irradiation embrittlement [6][7][8][9]. A previous study has reported that the ductile-brittle transition temperature (DBTT) of nuclear reactor pressure vessel (RPV) steels increases during the service period [8]. These steels are exposed to high-temperature and high neutron irradiation environments, and hence the embrittlement is attributed to partly GB P segregation. The GB embrittlement increases the damage risk of RPVs and must be considered during the aging management of the nuclear reactors. Therefore, quantitative prediction of the GB P segregation under various thermal and irradiation conditions is crucial for preventing this embrittlement [10,11].
We have been developing a rate theory model for numerical estimation of the GB P segregation based on the models reported in [12,13] . The model employs reaction-diffusion equations for three types of point defects (vacancies, self-interstitial atoms, and P atoms). P atoms migrate by interacting with vacancies or self-interstitial atoms. The migration of vacancies and self-interstitial atoms is also affected by the P atoms. Consequently, the partial diffusion coefficients in these equations are necessary for modeling the migration process and are evaluated using the kinetic Monte Carlo (kMC) simulation. These simulations incorporate the transition rate between the different states comprising a P atom and a vacancy or a P atom and a self-interstitial atom. The transition rate is evaluated using the barrier energy obtained from first-principles calculations. The diffusion coefficient of each defect is determined analytically or through kMC simulations [14,15]. Moreover, the rate theory model incorporates the model of trapping and de-trapping P atoms at GBs, and therefore, GB P segregation evaluated by the model can be compared with experimental measurement [16]. The trapping process is incorporated into the model based on the result of molecular dynamics (MD) simulations, in which a pair consisting of a P atom and a self-interstitial atom or a P atom and a vacancy dissociates before reaching a GB [17,18]. However, an intuitive, simple model is used for the de-trapping process, and hence a more detailed model is necessary. For a detailed consideration of the physical process of de-trapping, the migration of a P atom in a Σ3 [1][2][3][4][5][6][7][8][9][10](111) symmetrical tilt GB in α-iron (Fe) is examined using MD simulation [19]. However, owing to the abundance and diversity of GB types, obtaining a general understanding of the P migration in GBs and the P de-trapping process by considering only one GB type is difficult.
In this study, we analyzed a Σ5[100](0-13) GB and a Σ3[1-10](111) GB in α-Fe for MD simulation of P atom migration in these GBs and examined the effect of a vacancy on this migration. The simulation results obtained for the two types of GBs revealed some differences: P atoms migrate mainly via the interstitial mode in the Σ3 GB and mainly via the vacancy mode in the Σ5 GB. Moreover, the de-trapping phenomena of a P atom and a vacancy were investigated in these simulations.

Methods
The atomic models of the Σ3[1-10](111)GB and Σ5[100](0-13)GB of α-Fe (hereafter referred to as Σ3 GB and Σ5 GB, respectively) are shown in Figure 1. These GBs are symmetrical tilt GBs with comparatively large GB energy and large GB volume, as shown in Table 1. The figure shows a magnified view of the GB region considered in the simulation model (size and number of atoms are 42 Å × 97 Å × 49 Å, 14,689 for Σ3 GB and 45 Å × 100 Å × 43 Å, 14,400 for Σ5 GB). The model for the Σ3 GB has been used in our previous studies [15,19]. The GB was set on the x-z plane at the center of the y-direction, and a vacuum region with a thickness of about 6 Å was set at both edges of the y-direction. The thickness of the vacuum region was chosen based on the cut-off length of the employed potential. In the figure, the light gray particle denotes the Fe atom occurring in the GB region that is detected by the common neighbor analysis [20]. The numbers in the figure show the segregation sites of a P atom used in this simulation. After structural relaxation of the atomic state, a P atom and a vacancy was placed at one of the segregation sites in each GB region. The time development of each system was analyzed using the MD code LAMMPS [21] with the widely used embedded atom method (EAM) potential [22], which is considered the best potential for dilute P-Fe systems [23]. The potential is built by adding the fitting data obtained by the first-principles calculations for more specific atomistic level configurations such as Fe bcc crystal and P-Fe systems with and without an interstitial atom and a vacancy. The atomic image was visualized by a tool called OVITO [24]. In the simulation, the atoms on both edges of x and z directions were fixed to suppress the thermal migration of GBs. The temperatures of 600 K, 700 K, and 800 K were used for observing P migration and 800 K to 1400 K for evaluating the migration barrier energy. The temperature of 600 K is close to the operating temperature of nuclear reactors, and the higher temperatures, at which P migration is more marked, were employed for examining the behavior and the mechanism of P migration in the GBs. Furthermore, a time step of 5 fs for the numerical integration was chosen to allow proper simulation of atomic migration for a long time at the simulation temperature.

GBs without P Atoms
First, we examine the behavior of Fe atoms in the GB region without P atoms and take this as the reference scenario. Figure 2 shows the displacement of the Fe atoms. The red particles in the figure are the atoms that migrated more than one lattice unit. As shown in Figure 2a, the Fe atoms in the Σ3 GB region migrate frequently (those near the GB region also migrate). However, the Fe atoms of the Σ5 GB undergo only slight migration ( Figure 2b). It is confirmed that the situation is the same at 800 K. It can be owing to that the atomic spacing differs between these two GBs: Atomic alignment in the former are loose, while those in the latter tight. We noticed that when a vacancy is added by removing a Fe atom of Site 1 to the Σ5 GB ( Figure 1b), a few Fe atoms around the vacancy migrate as shown in Figure 2c; thus, the vacancy-exchange mode of migration is possible in this GB.  Figure 3 shows the ratios of the number of migrating Fe atoms to the total number of Fe atoms (the GB and intragranular regions are both considered). Here Fe atoms returning to the position of less than one lattice unit from their initial site at a measurement time are not counted as migrating Fe atoms. The cases of Σ3 GB with and without a vacancy, shown in Figure 3a,b, respectively do not show remarkable differences. The increase in the ratio of the migrating Fe atoms mitigates over time for both cases. The ratios for 700 K and 800 K become almost the same after time passes, but are significantly larger than that obtained for 600 K. Only the diagram with a vacancy is shown for the Σ5 GB in Figure 3, because Fe atoms in this GB hardly migrate without a vacancy. As shown in Figure 3c, the ratio of migrating Fe atoms for Σ5 GB with a vacancy is smaller than that of Σ3 GB in general (note the smaller scale of the vertical axis in Figure 3c for the Σ5 GB compared with that in Figure 3a,b for the Σ3 GB). Moreover, for times greater than ∼0.1 µs the ratio of migrating Fe atoms for 700 K is larger than that obtained for 800 K. This suggests that the increase in the number of migrating Fe atoms at 800 K stopped because of de-trapping of the vacancy from the GB region at about 0.03 µs. The trajectory of the vacancy (Figure 4) confirms the occurrence of vacancy de-trapping.  These results suggest that in the Σ3 GB, the added vacancy is absorbed in the GB region through the redistribution of atomic positions in this region. In the Σ5 GB, the space corresponding to the added vacancy remains as a vacancy and migrates by exchanging the vacancy site with an Fe atom.

GBs with a P Atom
Next, we examine the migration of a P atom in the GB regions. According to firstprinciples calculations [26], the trapping energy of a P atom to GBs depends on segregation sites. In Figure 1, Site 0 is an interstitial site and all the other segregation sites are substitutional. The P trapping energy of each site evaluated from the first-principles calculations is shown in Table 2 for each GB associated with the energy determined from the molecular static (MS) simulations with the EAM potential. Despite certain discrepancies between the results from the first-principles calculation and the MS simulation, both methods reveal the same tendency; namely, the trapping energies of Site 0 and Site 2 are larger than that of Site 3, and the energy of Site 1 is the smallest. Thus, in this study, we focus on the cases where a P atom is placed at Site 0 and at Site 2, which are the most stable segregation sites.  Figure 5 shows the P migration behavior as the magnitude of displacement of a P atom from the initial position. According to the figure, the P atom initially at Site 2 of the case of the Σ3 GB migrate frequently at a high temperature of 800 K, and barely migrates at a low temperature of 600 K (Figure 5b). When the P atom is inserted at the interstitial site (i.e., Site 0) or a vacancy is inserted at the nearest neighbor site of the substitutional P atom at Site 2, the atom started migrating even at the low temperature; however, the atom stops immediately (Figure 5a,c) in the Σ3 GB. We observed that the stopped P atom remains at the position of a substitute site in both cases. In the Σ5 GB without a vacancy, even at the high temperature, the migration of the P atom at Site 0 is less frequent than that in the Σ3 GB (Figure 5d). No migration of the atom initially placed at Site 2 is observed at either temperature (Figure 5e). However, when a vacancy is initially created at a nearest neighbor site of this atom (Figure 5f), the atom migrates at 800 K and the migration is more frequent than the no vacancy case shown in Figure 5d. At 600 K, the P atom oscillates around the initial position by exchanging the site with a vacancy. The results suggest that P atoms in Σ3 GB migrate mainly via the interstitial mode, while in the Σ5 GB they do via the vacancy mode. Figure 6 shows the trajectory of the P atom at 800 K for the case shown in Figure 5a,f. According to the figure, the P atom migrates through three-dimensional motion in the Σ3 GB region and via one-dimensional motion (linearly along the z-direction, i.e., normal to the plane of Figure 2) in the Σ5 GB region. In addition, Figure 6a shows that the P atom migrates outside of the GB region. The difference in P migration behavior between these two GBs can be attributed to the difference in the migration behavior of Fe atoms and of a vacancy described in the previous section.
We also analyzed the effect of a P atom on the migration of the GBs. For this analysis, some particular MD simulations are conducted, in which the atoms on both edges of the x and z directions are not stationary so that the GBs can migrate thermally. The results are shown in Figure 7. According to the figure, in the case without a P atom, the Σ3 GB migrated in the upper direction whereas the Σ5 GB migrates only slightly. However, Σ3 GB region stays in the initial position when a P atom is placed in the GB region, indicating that P atoms suppress the thermal migration of this GB. It is because Fe atoms around the P atom prefer energetically the atomic arrangement of the GB region rather than that of the intragranular region. The results also suggest that the behavior of thermal migration is different even in GBs with similar large GB energy and GB volume.  graphs (a,d) show the front view that is the x-y plane, the centers (b,e) are the side view that is the y-z plane, and the right is the top view that is x-z plane. The large red particle is the P atom at the last step. The color of Fe atoms is the same meaning as Figure 1. The green line shows the trajectory of the P atom.

The Evaluation of Migration Barrier Energy
Based on a previously reported method [27], which is described briefly below, we evaluate the barrier energy for P diffusion in both GBs for the cases without and with a vacancy. The barrier energy is evaluated from the temperature dependence of the diffusion coefficient D. Generally, D can be evaluated from the mean square displacement, msd as D = msd/6t where t is the time for recording the P atom migration process. The accuracy of evaluation, however, is limited by the length of the MD simulation. Thus, to achieve a certain level of accuracy, one trajectory of a P atom in a simulation is equally divided into many strips indexed by i, and D i = msd i /6t i is evaluated for each strip. We obtainD as the average of D i . Here, if the t i is too long, then the sample number of D i is not enough; if it is too short, then the value of D i includes too much statistical uncertainty. Thus, as shown in Figure 8a, the value ofD depends on the strip length. However, we observe thatD becomes nearly constant when the relation between t i and msd i is linear. Therefore, we evaluate the diffusion coefficient in the range of values corresponding to this linear relationship. The barrier energy for the obtainedD is evaluated for temperatures ranging from 800 K to 1400 K in the Arrhenius plot as shown in Figure 8b,c. The evaluated diffusion barrier energy is also shown in Table 3. According to the table, for the Σ3 GB, the P barrier energy started from a substitutional site (Site 2) with a vacancy is the almost same as that started from the interstitial site (Site 0). The result can indicate that the inserted vacancy disappears through redistribution of atomic positions in the Σ3 GB region, and P atoms migrate by the interstitial mode, as described in the previous section. For the Σ5 GB, the barrier energy associated with vacancy insertion at a nearest neighbor sites of a substitutional P atom is smaller than that of the interstitial P atom. Thus, unlike Σ3 GB, the result can indicate that the added vacancy migrates in the form of a vacancy and P atoms migrate via the vacancy mode in the GB. To recapitulate, the migration barriers derived here reflect the different modes of P migration in the different GBs.

Concluding Remarks
We investigated the P atom migration in two kinds of symmetrical tilt GBs, i.e., Σ3[1-10](111) GB and Σ5[100](0-13) GB, of α-Fe using MD simulations with an EAM potential. The effect of vacancy addition on P migration was also examined. In addition, the migration behavior of Fe atoms in the region of these two GBs without P atoms was investigated. The simulation results obtained are summarized as follows: 1.
In the Σ3 GB, Fe atoms migrated at 600 K, 700 K, and 800 K, and the addition of a vacancy has no effect on the migration behavior at these temperatures. Moreover, in the Σ5 GB, no Fe atom migration was observed even at 800 K; migration occurs when a vacancy was added, but was less frequent than that in the case of Σ3 GB. 2.
The de-trapping of a vacancy from the GB region of the Σ5 GB was observed at 800 K. 3.
The P atom migrated frequently regardless of the initial site, that is substitutional or interstitial, in the Σ3 GB at 800 K. Moreover, in the Σ5 GB, the P migration was observed when a vacancy was placed at one of the nearest neighbor sites of the substitutional P atom at 800 K. In addition, the P atom migration was a one-dimensional motion in the Σ5 GB and a three-dimensional motion in the Σ3 GB.

4.
In the Σ3 GB, the de-trapping of the P atom from the GB region was observed at 800 K.

5.
The Σ3 GB undergoes significant thermal migration, which is suppressed by a P atom, while the Σ5 GB undergoes only slight thermal migration. 6.
The diffusion barrier energy of a P atom started from the interstitial site was compared with that of a P atom starting from a substitutional site with a vacancy nearby; the evaluated value of both cases was almost the same in the case of Σ3 GB, while in the case of Σ5 GB, the value of the former was larger than that of the latter.
From these results, we can consider that vacancies inserted in the Σ3 GB disappear by redistributing the atomic positions in the GB region and P atoms migrate mainly by the interstitial mode. However, the form of the vacancy in the Σ5 GB is maintained and the vacancy migrates as a vacancy, and thus P atoms migrate mainly by the vacancy mode. In addition, since the de-trapping of the P atom and of a vacancy was observed, in the rate theory model of GB P segregation, the de-trapping process of P atoms at GBs should be considered, and it is necessary to consider GBs not only as the vacancy sink but also as the vacancy source. Data Availability Statement: This article has no additional data. All data included in this study are available upon request by contact with the corresponding author.