# Size-Dependent Solute Segregation at Symmetric Tilt Grain Boundaries in α-Fe: A Quasiparticle Approach Study

Helena Zapolsky, Antoine Vaugeois, Renaud Patte, Gilles Demange

► To cite this version:

Helena Zapolsky, Antoine Vaugeois, Renaud Patte, Gilles Demange. Size-Dependent Solute Segregation at Symmetric Tilt Grain Boundaries in α-Fe: A Quasiparticle Approach Study. Materials, 2021, 14 (15), pp.4197. 10.3390/ma14154197. hal-03472590

HAL Id: hal-03472590

https://hal.science/hal-03472590v1

Submitted on 26 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Size-Dependent Solute Segregation at Symmetric Tilt Grain Boundaries in α-Fe: A Quasiparticle Approach Study

Helena Zapolsky*,† ID, Antoine Vaugeois†, Renaud Patte † ID and Gilles Demange*,†

GPM, UMR CNRS 6634, Université de Rouen-Normandy, 76575 Saint Étienne du Rouvray, France;
antoine.vauggeois@univ-rouen.fr (A.V.); renaud.patte@univ-rouen.fr (R.P.)

* Correspondence: helena.zapolsky@univ-rouen.fr (H.Z.); gilles.demange@univ-rouen.fr (G.D.)

† These authors contributed equally to this work.

**Abstract:** In the present work, atomistic modeling based on the quasiparticle approach (QA) was performed to establish general trends in the segregation of solutes with different atomic size at symmetric (100) tilt grain boundaries (GBs) in α-Fe. Three types of solute atoms X₁, X₂ and X₃ were considered, with atomic radii smaller (X₁), similar (X₂) and larger (X₃) than iron atoms, respectively, corresponding to phosphorus (P), antimony (Sb) and tin (Sn). With this, we were able to evidence that segregation is dominated by atomic size and local hydrostatic stress. For low angle GBs, where the elastic field is produced by dislocation walls, X₁ atoms segregate preferentially at the limit between compressed and dilated areas. Contrariwise, the positions of X₂ atoms at GBs reflect the presence of tensile and compressive areal regions, corresponding to extremum values of the σ<sub>XX</sub> and σ<sub>YY</sub> components of the strain tensor. Regarding high angle GBs Σ5 (310) (θ = 36.95°) and Σ29 (730), it was found that all three types of solute atoms form Fe<sub>9</sub>X clusters within B structural units (SUs), albeit being deformed in the case of larger atoms (X₂ and X₃). In the specific case of Σ29 (730) where the GB structure can be described by a sequence of |BC.BC| SUs, it was also envisioned that the C SU can absorb up to four X₁ atoms vs. one X₂ or X₃ atom only. Moreover, a depleted zone was observed in the vicinity of high angle GBs for X₂ or X₃ atoms. The significance of this research is the development of a QA methodology capable of ascertaining the atomic position of solute atoms for a wide range of GBs, as a mean to highlight the impact of the solute atoms' size on their locations at and near GBs.

**Keywords:** grain boundary; quasiparticle approach (QA); segregation; phase-field crystal (PFC)

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Zapolsky, H.; Vaugeois, A.; Patte, R.; Demange, G. Size-Dependent Solute Segregation at Symmetric Tilt Grain Boundaries in α-Fe: A Quasiparticle Approach Study.
*Materials* **2021**, *14*, 4197. https://doi.org/10.3390/ma14154197

Academic Editors: Rebecca Janisch
and Sergiy V. Divinskiy

Received: 31 May 2021

Accepted: 22 July 2021

Published: 27 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The physical properties of polycrystalline materials (embrittlement, electric transport, and corrosion) are often driven by the segregation of impurities at grain boundaries (GBs) [1–5]. Manifold experimental studies have been carried out to investigate the interaction of impurities with GBs [6–10]. Nowadays, it is well established that sulfur (S) and phosphorous (P) atoms are embrittling elements in steel, while carbon (C) and boron (B) are known as flagship GB cohesion enhancers [11].

In addition to experiments, computer simulation of GB segregation have ushered in atomic level investigations of the segregation process. In this regard, significant progress has been made to decipher the interaction between solute atoms and GBs by means of molecular dynamic (MD) and ab initio calculations [12–17]. Ab-initio methods are frequently used to evaluate the most stable segregation sites. For instance, J. Wang et al. [18] looked into segregation at symmetric tilt grain boundaries in α-Fe. It was evidenced that the maximum fracture strength of a GB depends on the maximum carbon concentration that can be accommodated by these GBs. Y. Hu et al. [19] studied the segregation effects of six transition metal elements (Cr, Ni, Cu, Zr, Ta, and W) on the Σ3 (111) tilt boundary in bcc iron and concluded that the segregation of Zr, Ni and Cu elements decreases the
GB cohesive strength, as opposed to the strengthening effect of Cr, Ta, and W atoms. Ab initio calculations also entailed meaningful progresses regarding the comprehension of the phenomenon of hydrogen embrittlement of metals and alloys [20,21]. Despite paramount breakthroughs in the field of GB interactions with solute atoms enabled by ab initio calculations, this class of approaches remains computationally expensive, hence its restricted use to a few GBs at 0 K. Alternatively, MD modeling spearheaded a host of new results on GBs segregation. For example, N. R. Rhode et al. [22] examined a segregation of carbon atoms on a large number of GBs, including general low and high angle GBs. In particular, they underlined the influence of the local structure of the GB on the segregation energy. The main result of this work was to demonstrate that the atomic sites at and close to GBs show an asymmetric distribution of segregation energies displaying extreme values that extend over 10 Å from the grain boundary. Albeit significantly less computationally demanding than the ab initio calculations, the MD simulation is limited to the description of systems on relatively short time scales, so it struggles to reproduce the diffusion of solute atoms toward GB. Notwithstanding this shortcoming, MD modeling allowed to gather a comprehensive database for the segregation spectra of 250+ binary alloys [23], which was successfully exploited in the machine learning framework to predict the segregation energy of a solute atom in a GB site. However, new methods are now required to improve this database and provide insight on the full segregation phenomenon, starting from the solute diffusion to GB, up to equilibrium segregation at GBs.

In this regard, an orthogonal approach to MD is provided by the phase-field crystal (PFC) model [24–26]. This approach was first curtailed to bcc symmetry, and thereupon expanded to arbitrary crystallography [27] under the banner of XPFC [28]. Despite the great success of the XPFC method in describing GB structures [29], this approach is currently not suitable to account for atoms of different radii, which hinders the assessment of the solute atoms radius influence on GB segregation patterns. This shortcoming can be circumvented by an alternative model to XPFC, coined the quasiparticle approach (QA). This method was the outcome of the extension of the atomic density function approach (ADF) [30] to the continuum case [31]. With this, it opened a way to model the evolution of different aperiodic systems (such as GBs or glass) or displacive phase transformations. The QA has already been applied to model the symmetric tilt grain boundaries structure [32] and vacancies annihilation [33,34] at GBs in α-Fe, the self-assembly of atoms into complex structures [35], and fcc/bcc phase transformations [36].

In this work, we propose to use the QA to study the segregation of solute atoms at symmetric (100) tilt GBs in α-Fe. The main objective of this paper is to provide insight on how the local structure of GB can affect the segregation of various elements with different atomic radii. Specifically, the present work aims at disentangling the strain and chemical contributions of solute atoms segregation at GBs, upon modulating the atomic radii of solute atoms, while keeping the chemical interactions fixed. For this purpose, three types of solute atoms with different atomic radii and two classes of GBs were considered: low angle GBs (LAGBs), with $\theta = 7.15^{\circ}$ and $\theta = 9.53^{\circ}$, and two high angle GBs (HAGBs), $\Sigma5$ (310) ($\theta = 36.95^{\circ}$) and $\Sigma29$ (730) ($\theta = 46.40^{\circ}$). This choice was motivated by the fact that $\Sigma5$ (310) is a low fit index symmetric (special) GB in the coincident site lattice (CSL) theory [37] characterized by a low segregation tendency [4,38]. Contrariwise, $\Sigma29$ (730) is a high index (general) GB, which displays a high segregation tendency. Finally, LAGBs ($\theta < 15^{\circ}$) display a dislocation wall structure [39].

This paper is organized as follows. First, an overview of the quasiparticle approach with application to the binary system is presented, with special emphasis on the choice of model parameters. The QA is then applied to model the solute diffusion and segregation at symmetric (100) tilt GBs in α-Fe. Finally, we compare our results with the available ab initio calculations. We demonstrate that the proposed model gives a good description of the solute atoms distribution at the considered GBs in α-Fe.
## 2. Quasiparticle Approach (QA)

### 2.1. Basic Equation

In this section, the main equations of the quasiparticle approach (QA) are briefly introduced. This method can be seen as extension of the seminal atomic density function (ADF) approach [30] to the continuum case. In the ADF theory, the atomic configuration is described by a set of occupation probabilities $P_α(r, t)$, defined as the probability for a given lattice site $r$ to be occupied by the atom of type $α$, at a given time $t$. The temporal evolution of these variables is governed by Onsager-type diffusion equations, wherein it is proportional to the thermodynamic driving force [30]. In the ADF approach, the probability function is specified at each site of the underlying Ising lattice, which coincides with the simulation grid, thereby confining the application range of the model to isostructural phase transformations. This shortcoming was later circumvented in [31] upon choosing the simulation grid spacing several times smaller than the interatomic distance. With this, the so-called continuum atomic density function (CADF) theory allowed to account for the atomic movement in the continuous space, hence making the model applicable to structural phase transitions. In the CADF framework, atoms are no longer points, but rather spheres containing a certain number of the simulation grid nodes (see Figure 1). In this regard, a new interaction Hamiltonian should be defined to set the dynamics of the system. One such Hamiltonian was proposed in an upgrade of the CADF model referred to as the quasiparticle approach (QA) [35].

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic representation of the simulation grid (gray solid lines connecting grid nodes) for the ADF model (a) and the CADF theory and QA (b).</figcaption></figure>

The salient feature of this second extension of the ADF theory relies on the treatment of lattice points belonging to atomic spheres as non-traditional dynamic variables called fratons. In the QA, the configurational degrees of freedom are occupation numbers $c(r)$. The function $c(r)$ is equal to 1 if the lattice site in position $r$ lies within the atomic sphere, and is 0 otherwise. In this formalism, the occupation numbers of the fratons are dynamic variables of the system rather than the coordinates of atomic spheres in the conventional description of the configurational phase space. Therefore, a $m$-components system can be characterized by the values of $m$ stochastic numbers $c_α(r)$ in each lattice site $r$, where $α = 1, 2, ..m$ labels the fratons related to the corresponding atom of atomic species $α$. Occupation variables $c_α$ are then averaged over the time-dependent Gibbs ensemble into the occupation probability $ρ_α(r, t) = ⟨c_α(r, t)⟩$, where the $⟨·⟩$ symbol denotes the Gibbs ensemble average at temperature $T$ and time $t$. With this definition, the function $ρ_α(r, t)$ is the probability that a lattice point in $r$ is located anywhere inside the atomic sphere of any atom of the kind $α$ at time $t$. Therefore, the atomic configuration of the system can be fully described by the density function $ρ_α(r, t)$. Moreover, the temporal evolution of the system is given by the microscopic diffusion equation [31]:
$$ \frac{\partial \rho_{\alpha}}{\partial t}(\mathbf{r}, t) = \sum_{\beta=1}^{m} \sum_{\mathbf{r}' \in I} L_{\alpha\beta}(\mathbf{r} - \mathbf{r}') \frac{\delta F}{\delta \rho_{\beta}(\mathbf{r}', t)} \quad (1) $$

where the summation is carried out over the $N_0$ grid sites of the Ising lattice I. In Equation (1), $L_{\alpha\beta}(\mathbf{r} - \mathbf{r}')$ is the matrix of kinetic coefficients, where indices $\alpha$ and $\beta$ label fratons describing the different kinds of atoms ($\alpha = 1, 2, m$), and $F$ is the non-equilibrium Helmholtz free energy functional. Kinetic Equation (1) approximates the evolution rate of the density functions $\rho_{\alpha}(\mathbf{r}, t)$ by the first non-vanishing term of its expansion with respect to the thermodynamic driving force in the small driving force limit. This microscopic diffusion equation is significantly nonlinear with respect to the density field $\rho_{\alpha}(\mathbf{r}, t)$, but it is linear with respect to the driving force. To guarantee a conservation of the total number of fratons of the kind $\alpha$, the kinetic coefficients matrix should satisfy the following condition for all $\alpha, \beta = 1, \dots, m$:

$$ \sum_{\mathbf{r} \in I} L_{\alpha\beta}(\mathbf{r}) = 0. \quad (2) $$

As for the free energy $F$ of the system, it is defined under mean-field approximation by the following:

$$ F = \sum_{\alpha=1}^{m} \sum_{\substack{\beta=1 \\ \beta \ge \alpha}}^{m} \left[ \frac{1}{2} \sum_{\mathbf{r}, \mathbf{r}' \in I} W_{\alpha\beta}(\mathbf{r}-\mathbf{r}') \rho_{\alpha}(\mathbf{r}, t) \rho_{\beta}(\mathbf{r}', t) \right] \\ + k_B T \sum_{\mathbf{r} \in I} \left[ \sum_{\alpha=1}^{m} \rho_{\alpha} \ln(\rho_{\alpha}) + \left(1 - \sum_{\alpha=1}^{m} \rho_{\alpha}\right) \ln\left(1 - \sum_{\alpha=1}^{m} \rho_{\alpha}\right) \right], \quad (3) $$

where $k_B$ is the Boltzmann constant and $T$ is the temperature. In this expression, the first term on the right of the equality corresponds to the internal energy, while the term on the right stands for the configurational entropy. $W_{\alpha\beta}$ is the pairwise interaction potential between fratons of type $\alpha$ and $\beta$ separated by a distance $|\mathbf{r} - \mathbf{r}'|$. It is noteworthy that the mean-field approximation [38] posited in Equation (3) is asymptotically accurate at low and high temperatures, while its precision peaks when the interaction radius largely exceeds the distance between interacting particles [38]. Thus, the smaller the grid spacing, the more accurate the description of the continuous atomic movements. However, refining the simulation lattice feeds through an increment of the computational cost, so a trade-off between a fine simulation and computational efficiency should be reached.

In the QA, the model fraton–fraton pair potential $W_{\alpha\beta}$ embodies the so-called short range (SR) and long range (LR) interactions, respectively written $\theta_{\alpha}$ and $W_{\alpha\beta}^{LR}$:

$$ W_{\alpha\beta}(\mathbf{r} - \mathbf{r}') = \theta_{\alpha}(\mathbf{r} - \mathbf{r}')\delta_{\alpha\beta} + \lambda_{\alpha\beta}W_{\alpha\beta}^{LR}(\mathbf{r} - \mathbf{r}'). \quad (4) $$

Here, $\delta_{\alpha\beta}$ is the Kronecker delta function. The parameter $\lambda_{\alpha\alpha}$ ($\alpha = \beta$) is the relative amplitude between LR and SR interactions. In this work, we set $\lambda_{\alpha\alpha} \equiv \lambda$ for all components $\alpha = 1, \dots, m$. No SR cross interactions ($\alpha \neq \beta$) are considered, so $W_{\alpha\beta} \equiv \lambda_{\alpha\beta}W_{\alpha\beta}^{LR}$. The SR contribution $\theta_{\alpha}$ allows the spontaneous condensation of fratons into atomic spheres, and prevents the overlap or atoms. To reproduce this specific behavior of fratons at short distances, the step function depicted in Figure 2 was used for SR interactions [35]:

$$ \theta_{\alpha}(\mathbf{r}) = \begin{cases} -1 & \text{if } r \le R_{\alpha} \\ \zeta & \text{if } R_{\alpha} < r \le R_{\alpha} + \Delta R_{\alpha} \\ 0 & \text{otherwise,} \end{cases} \quad (5) $$
<figure><img src="image_5.png" /><figcaption>Figure 2. SR interaction $\theta_{\alpha}(r)$ step function profile.</figcaption></figure>

Here, $R_{\alpha}$ sets the width of the attractive part of the SR potential, which determines the radius of atomic spheres for component $\alpha$. Then, $\Delta R_{\alpha}$ and $\zeta$ are respectively the width and height of the SR potential barrier for each component $\alpha$. The introduction of a repulsive contribution in the SR potential not only prevents the overlap of atomic spheres, but also contributes to adjust the elastic properties of the system. For convenience, the interaction potential $W_{\alpha\beta}$ is implemented in reciprocal space by means of the Fourier transforms $\hat{\theta}_{\alpha}(\mathbf{k})$ and $\hat{W}_{\alpha\beta}^{LR}(\mathbf{k})$ of the SR and LR interactions, where $\mathbf{k}$ is the k-vector defined by $\mathbf{k} = (k_x, k_y, k_z) = (\frac{2\pi}{N_x}h, \frac{2\pi}{N_y}k, \frac{2\pi}{N_z}l)$ with $(h, k, l)$ being dimensionless coordinates, and $N_x, N_y$ and $N_z$ being the number of grid nodes on each edge of the simulation domain. Equation (4) thus becomes the following:

$$ \hat{W}_{\alpha\beta}(\mathbf{k}) = \hat{\theta}_{\alpha}(\mathbf{k})\delta_{\alpha\beta} + \lambda_{\alpha\beta}\hat{W}_{\alpha\beta}^{LR}(\mathbf{k}), \quad (6) $$

where the Fourier transformation $\hat{\theta}_{\alpha}(\mathbf{k})$ of the SR contribution reads as follows:

$$ \hat{\theta}_{\alpha}(\mathbf{k}) = \frac{4\pi}{k^3} \left[ -\sin(kR_{\alpha}) + kR_{\alpha}\cos(kR_{\alpha}) + \zeta \left\{ \sin(k(R_{\alpha} + \Delta R_{\alpha})) \right. \right. \\ \left. - k(R_{\alpha} + \Delta R_{\alpha})\cos(k(R_{\alpha} + \Delta R_{\alpha})) \right. \\ \left. - \sin(kR_{\alpha}) + kR_{\alpha}\cos(kR_{\alpha}) \right\} \right]. \quad (7) $$

In this work, the segregation of one solute species X at a GB in bcc iron is addressed. Henceforth, the first component ($\alpha = 1$) will conventionally pertain to Fe, while the second component ($\alpha = 2$) will refer to the solute species X. The corresponding occupation probabilities will then be written $\rho_1(\mathbf{r}, t) \equiv \rho_{\text{Fe}}(\mathbf{r}, t)$ and $\rho_2(\mathbf{r}, t) \equiv \rho_{\text{X}}(\mathbf{r}, t)$ in the remaining part of the paper.

The long range interaction potential $\hat{W}_{\alpha\beta}^{LR}(\mathbf{k})$ defines the crystal structure, the elastic properties, and the chemical interactions between atoms in the system. In this work, only spherically symmetric potentials were used as a means to allow the formation of crystallographic structures with arbitrary orientation. The potential ansatz introduced in [31] was selected to describe the (bcc) crystallographic structure and the elastic properties of $\alpha$-Fe. As was shown in [32] for a single-component system, the LR contribution can be conveniently fitted on the structure factor $S(k)$ of a given system close to the melting point. In this work, the $S(k)$ function for the bcc-iron calculated in [40] was fitted by the following function:

$$ \hat{W}_{11}^{LR}(k) \equiv \hat{W}_{\text{Fe}}^{LR}(k) = 1 - \frac{k^4}{(k^2 - k_1^2)^2 + k_2^4}. \quad (8) $$
Here, $k_1$ and $k_2$ are two model parameters that notably set the characteristic wave-length $k_0 = \frac{\sqrt{k_1^4+k_2^4}}{k_1}$ of the bcc crystal lattice, for the latter value minimizes the function in Equation (8). In the bcc structure, each atom is surrounded by 8 nearest neighbors, and the first neighbor distance between atoms is given by $d = a_0\sqrt{3}/2$, where $a_0$ is the lattice parameter of α-iron. Recalling that the maximum allowed radius for an atom embedded in a crystal structure is half the nearest neighbor distance $d$, the following condition for the step functions $\theta_{Fe,X}$ defined in Equation (5) holds:

$$R_{Fe,X} + \Delta R_{Fe,X} \le a_0 \sqrt{3}/4. \qquad (9)$$

The reciprocal lattice of the bcc crystal is the fcc structure with the reciprocal lattice parameter $4\pi/a_0$. It ensues that the first structural reflection is located at $k_0 = 2\sqrt{2\pi}/a_0$. As was discussed in [35,41], a multi-minima potential is required to model more complex structures. From a general perspective, the number of minima of such LR potential should be equal to the number of non-equivalent structural reflections in the first Brillouin zone of the system.

Under the dilute solution hypothesis, we assume in this work that $\hat{W}_{22}^{LR}(k) \equiv \hat{W}_{X}^{LR}(k) = 0$. Moreover, a fully repulsive interaction was considered for the long-range potential $\hat{W}_{12}^{LR}(k)$ between solute and Fe atoms. This contribution was implemented by a simple Gaussian function centered in $k = 0$ as a means to foster the phase separation of the different chemical species:

$$\hat{W}_{12}^{LR}(k) \equiv \hat{W}_{Fe-X}^{LR}(k) = -\exp\left(-\frac{k^2}{2\sigma^2}\right). \qquad (10)$$

The fitting parameter $\sigma$ tailors the range of the Fe-X repulsion, depending on the solute species X. As for the cross interaction potential, we recall that $\hat{W}_{12} = \lambda_{12}\hat{W}_{12}^{LR}$ in Equation (6). Using the same notations as in Equations (8)–(10) (index 1 for Fe, index 2 for X), this relation reads $\hat{W}_{Fe-X} = \lambda_{Fe-X}\hat{W}_{Fe-X}^{LR}$. With this, the parameter $\lambda_{Fe-X}$ weights the relative influence of the iron structural contribution $\hat{W}_{Fe}$ with respect to the chemical repulsion $\hat{W}_{Fe-X}^{LR}$ between iron and solute X atoms.

It should be also pointed out that each grid node can be occupied or not by a fraton. Then, in order to describe the atomic configuration in a binary system, the distribution of Fe, X and vacancies (V) should be considered. However, according to the conservation condition $\rho_{Fe} + \rho_X + \rho_V = 1$, only two fraton density functions should be defined. With indices 1 and 2 referring to Fe and X atoms, kinetic Equation (1) for the density probability functions $\rho_{Fe,X}$ was solved in Fourier space:

$$\begin{aligned} \frac{\partial \hat{\rho}_{Fe}}{\partial t}(\mathbf{k}, t) &= \hat{L}_{Fe}(k) \left[ \hat{W}_{Fe}(k) \hat{\rho}_{Fe}(\mathbf{k}, t) + \hat{W}_{Fe-X}(k) \hat{\rho}_{X}(\mathbf{k}, t) + k_B T \left\{ \ln \left( \frac{\rho_{Fe}}{1 - \rho_{Fe} - \rho_X} \right) \right\}_k \right] \\ &\quad + \hat{L}_{Fe-X}(k) \left[ \hat{W}_{X}(k) \hat{\rho}_{X}(\mathbf{k}, t) + \hat{W}_{Fe-X}(k) \hat{\rho}_{Fe}(\mathbf{k}, t) + k_B T \left\{ \ln \left( \frac{\rho_X}{1 - \rho_{Fe} - \rho_X} \right) \right\}_k \right] \\ \frac{\partial \hat{\rho}_{X}}{\partial t}(\mathbf{k}, t) &= \hat{L}_{X}(k) \left[ \hat{W}_{X}(k) \hat{\rho}_{X}(\mathbf{k}, t) + \hat{W}_{Fe-X}(k) \hat{\rho}_{Fe}(\mathbf{k}, t) + k_B T \left\{ \ln \left( \frac{\rho_X}{1 - \rho_{Fe} - \rho_X} \right) \right\}_k \right] \\ &\quad + \hat{L}_{Fe-X}(k) \left[ \hat{W}_{Fe}(k) \hat{\rho}_{Fe}(\mathbf{k}, t) + \hat{W}_{Fe-X}(k) \hat{\rho}_{X}(\mathbf{k}, t) + k_B T \left\{ \ln \left( \frac{\rho_{Fe}}{1 - \rho_{Fe} - \rho_X} \right) \right\}_k \right], \end{aligned} \qquad (11)$$

where $\hat{\rho}_{Fe,X}(\mathbf{k}, t) = \sum_r \rho_{Fe,X}(\mathbf{r}, t) \exp(i\mathbf{k} \cdot \mathbf{r})$ is the Fourier transform of the density function $\rho_{Fe,X}$, and $\{\cdot\}_k$ is the discrete Fourier transform operator. Moreover, $\hat{L}_{Fe,X,Fe-X}(k) = -L_{Fe,X,Fe-X}^{Ons}k^2$, where the coefficients $L_{Fe,X,Fe-X}^{Ons}$ are Onsager diffusion coefficients. In the present model, the interaction between vacancies and atoms as well the vacancy–vacancy interactions were neglected. Under this simplification, the Onsager diffusion coefficient matrix reads as follows:

$$L^{\text{Ons}} = \begin{pmatrix} L_{\text{Fe}}^{\text{Ons}} & L_{\text{Fe-X}}^{\text{Ons}} \\ L_{\text{Fe-X}}^{\text{Ons}} & L_{X}^{\text{Ons}} \end{pmatrix}. \qquad (12)$$
Different values of the coefficients of $L^{\text{Ons}}$ matrix were assessed in this work (not shown). It was observed that the equilibrium state obtained by solving Equation (11) remains quite insensitive to the exact values chosen for these coefficients, provided that the matrix is positive definite. With this, Onsager coefficients were set to $L_{\text{Fe}}^{\text{Ons}} = L_{\text{X}}^{\text{Ons}} = 1$ and $L_{\text{Fe-X}}^{\text{Ons}} = 0.5$ in all simulations.

Equation (11) was solved in reduced units. The average density of probabilities $\bar{\rho}_{\text{Fe}}$ and $\bar{\rho}_{\text{X}}$ of matrix (Fe) and solute (X) atoms is defined as $4\pi R_{\text{Fe,X}}^3 N_{\text{Fe,X}} / (3V)$, where $V = (\Delta x)^3 N_x N_y N_z$ is the total volume of the system, and $N_{\text{Fe}}$ ($N_{\text{X}}$) is the total number of Fe (X) atoms at ground state. The input parameters related to the energy $\zeta$, $\lambda$ and $\lambda_{\text{Fe-X}}$ are expressed in $k_B T_m$ units, where $T_m$ is the melting temperature of the pure iron system with composition $\bar{\rho}_{\text{Fe}}$. Lengths are expressed in units of the lattice parameter $\Delta x$, so the grid spacing $\Delta x$ of the underlying Ising lattice is defined as a fraction of the lattice parameter $a_0$ of Fe.

## 2.2. Model Parameters

In this work, the effect of the solute atom radius on the segregation at GBs in $\alpha$-iron ($a_0 = 2.87$ Å) was investigated. With this in mind, solute atoms $X_1$, $X_2$ and $X_3$ with three different ionic radius tantamount to that of phosphorus (P), antimony (Sb) and tin (Sn) were considered. The corresponding values of $R_{\text{Fe,X}_{1,2,3}}$ and $\Delta R_{\text{Fe,X}_{1,2,3}}$ complying with condition (9) are provided in Table 1. The resulting SR interaction profiles defined in Fourier space in Equation (7) could then be set, as displayed in Figure 3a.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Pairwise interaction potentials used in this work in Fourier space: (a) short range interaction potential $\hat{\theta}_{\alpha}(k)$ for Fe, $X_1$, $X_2$ and $X_3$ atoms, (b) long range interaction potential $\hat{W}_{\text{Fe}}^{\text{LR}}$ (blue) and $\hat{W}_{\text{Fe-X}}^{\text{LR}}$ (red) for $X=\text{X}_1$, $X_2$ and $X_3$ using parameters of Table 1.</figcaption></figure>

<table><caption>**Table 1.** Parameters used in this work for QA simulations.</caption><thead><tr><th></th><th>a<sub>0</sub>(Δx)</th><th>R(Δx)</th><th>ΔR</th><th>ζ</th><th>λ</th><th>k<sub>0</sub></th><th>k<sub>1</sub></th><th>k<sub>2</sub></th><th>λ</th><th>λ<sub>Fe-X</sub></th><th>σ</th><th>k<sub>B</sub>T</th><th>ρ̄</th><th>L<sup>Ons</sup><sub>Fe</sub></th><th>L<sup>Ons</sup><sub>Fe-X</sub></th></tr></thead><tbody><tr><td>Fe</td><td>16</td><td>6.15</td><td>1.17R</td><td>5.0</td><td>0.5</td><td>0.555</td><td>0.242</td><td>0.348</td><td>0.5</td><td>0.1</td><td>0.1</td><td>0.0275</td><td>0.104</td><td>1.0</td><td>0.5</td></tr><tr><td>X<sub>1</sub></td><td>-</td><td>4.2</td><td>1.17R</td><td>3.0</td><td>0.5</td><td>-</td><td>-</td><td>-</td><td>0.5</td><td>0.1</td><td>0.1</td><td>0.0275</td><td>0.02</td><td>1.0</td><td>0.5</td></tr><tr><td>X<sub>2</sub></td><td>-</td><td>5.7</td><td>1.17R</td><td>3.0</td><td>0.5</td><td>-</td><td>-</td><td>-</td><td>0.5</td><td>0.1</td><td>0.1</td><td>0.0275</td><td>0.02</td><td>1.0</td><td>0.5</td></tr><tr><td>X<sub>3</sub></td><td>-</td><td>6.64</td><td>1.17R</td><td>3.0</td><td>0.5</td><td>-</td><td>-</td><td>-</td><td>0.5</td><td>0.1</td><td>0.1</td><td>0.0275</td><td>0.02</td><td>1.0</td><td>0.5</td></tr></tbody></table>

To fix the model parameters introduced in Equation (10) for LR interactions between iron atoms, the elastic constants of $\alpha$-Fe were reckoned, using the following procedure. For
small deformations where materials exhibit a linear elastic behavior, the free energy can be expanded in Taylor series with respect to the deformation:

$$F(\{\epsilon_k\}) = F_0 + \frac{V_0}{2} \sum_{m,n=1}^{6} C_{mn} \epsilon_m \epsilon_n, \qquad (13)$$

where $\epsilon_k$ is the component of the rank 2 strain tensor $\bar{\epsilon}$, $F_0$ and $V_0$ are the free energy and volume of the unconstrained system respectively, and $C_{mn}$ are the elastic constants in Voigt notation for a cubic system. The latter can be evaluated through the second derivative of the free energy:

$$C_{mn} = \frac{\partial^2 F}{\partial \epsilon_m \partial \epsilon_n}. \qquad (14)$$

To characterize the elastic properties of cubic system, only three independent elastic constants, namely $C_{11}$, $C_{12}$ and $C_{44}$ must be defined. This can be achieved by applying the three characteristic deformations with specific transformation of coordinates to the system:

* Hydrostatic: $(x, y, z) \rightarrow ((1 - \epsilon)x, (1 - \epsilon)y, (1 - \epsilon)z)$

* Orthorhombic: $(x, y, z) \rightarrow ((1 + \epsilon)x, (1 - \epsilon)y, z)$

* Monoclinic: $(x, y, z) \rightarrow (x + \epsilon y, y, z)$

where the coefficient $\epsilon$ is the amplitude of the deformation. It is therefore possible to express the free energy of the system associated with these three deformations as a function of the strain, according to Equation (13):

$$ \begin{aligned} F_{\text{hydro}} &= F_0 + \frac{9}{2} V_0 B \epsilon^2 \\ F_{\text{ortho}} &= F_0 + 2 V_0 C \epsilon^2 \\ F_{\text{mono}} &= F_0 + \frac{1}{2} V_0 C_{44} \epsilon^2, \end{aligned} \qquad (15) $$

where $B = (C_{11} + 2C_{12})/3$ is the bulk modulus, and $C' = (C_{11} - C_{12})/2$. Another important parameter called the Zener anisotropy ratio $A_Z$ can be estimated from the elastic constants as follows:

$$A_Z = \frac{2C_{44}}{C_{11} - C_{12}} = \frac{C_{44}}{C'}. \qquad (16)$$

The case $A_Z = 1$ corresponds to an isotropic material, whereas $A_Z \neq 1$ indicates that the crystal is elastically anisotropic. To evaluate the free energy under three different deformations, a simulation box of 128³ was used with the next set of model parameters: $a_0 = 16\Delta x$, $R = 6.15\Delta x$ and $\Delta R = 0.7\Delta x$, while different values of $\zeta$ in the SR interactions (Equation (7)) and $k_1$ and $k_2$ in the LR interactions (Equation (8)) for $\alpha$-Fe were analyzed. For $\zeta = 5$, $k_1 = 0.435k_0$ and $k_2 = 0.626k_0$, the following set of elastic constants was found: $C_{11} = 163$ GPa, $C_{12} = 67$ GPa and $C_{44} = 95$ GPa, upon preliminarily fitting the numerical value of $C'$ on its ab initio counterpart [42]. The corresponding values of $B$ and $A_Z$ are 99 GPa and 2.0, respectively. A perfect match with [42] is achieved for $C_{44}$ ($C_{44} = 96$ GPa), while the obtained value for $A_Z$ falls within the range of numerical (ab initio) and experimental values spanned by the literature ($1.5 \le A_Z \le 2.7$ [42–46]). However, the present bulk modulus $B$ is too low with respect to the values reported in the literature ($168 \le B \le 189$ [42–46]). This discrepancy for $B$ is likely to reflect the small number of free parameters stepping in the expression of LR interactions. One should keep in mind that this underestimation of the bulk modulus of bcc iron in the QA makes the presently modeled iron too soft, compared to experiments. In turn, this may slightly influence the tendency of solute atoms to segregate at GBs depending on their ionic radius. We surmised in the present work that a first qualitative connection between the radius of solute atom and their segregation tendency at GBs could be made. This discrepancy should nonetheless be addressed in a future study, using a more sophisticated interaction potential.
The LR interaction potentials for α-Fe hereby fitted is depicted in Figure 3b (blue), along with the cross interaction potential $\hat{W}_{\text{Fe-X}}^{\text{LR}}(k)$ (red). In the present work, the same values for $\lambda_{\text{Fe-X}}$ and $\sigma$ were used for X₁, X₂ and X₃ in Equation (10). In this manner, the different solute species only depart from one another through their atomic radius, hence inducing different atomic misfits in the bcc iron lattice. In doing so, we also tacitly ascribe a specific segregation tendency at GBs to atomic the size effect, for purposes of disentangling the chemical and elastic mutual influence on the segregation at GBs. All parameters are compiled in Table 1.

In this work, the specific case of the segregation of one chemical species X at $\langle 100 \rangle$ symmetric tilt grain boundaries with the (010) interface plane in α-Fe was studied using the following procedure. First, the GB structure with specific misorientation angle $\theta$ was obtained by crystallizing a liquid layer placed in between two bcc crystal grains rotated by $\theta/2$ around the $\langle 100 \rangle$ axis. In order to find the minimum energy of the system during the crystallization stage, Equation (1) was integrated until equilibrium was reached. Then solute atoms (X) were introduced in substitutional position with a density of presence of 10%. To satisfy the periodic boundary conditions, two GB were introduced in the simulation box.

The space scale was chosen as the grid spacing $\Delta x = 0.018$ nm, as imposed by the number of grid lattices (16) spanning one bcc iron lattice parameter $a_0 = 16\Delta x = 2.87$ Å. Simulations were performed in three dimensions on a $600 \times 600 \times 64$ grid lattice ($N_x = N_y = 600, N_z = 64$) equipped with periodic boundary conditions. For the chosen length scale ($\Delta x \approx 0.018$ nm), the latter corresponds to a volume of (11 nm)$^3$. Kinetic Equation (11) was solved by the spectral-Eyre scheme [41] with the reduced time step $\Delta t = 0.005$, on the supercomputer CRIANN of Normandy.

# 3. Results

## 3.1. Segregation Pattern at GB

### 3.1.1. Low Angle GB

In the present work, we have assessed how the size of solute atoms influences the segregation phenomenon, starting with LAGBs. For that purpose, the segregation of X₁ and X₂ atoms was investigated at two LAGBs. As a reminder, the radius of X₁ and X₂ solute atoms were chosen close to the ionic radius of phosphorus (P) and antimony (Sb), respectively (radius $R_{X_1} < R_{X_2}$). Tilt symmetrical LAGBs with misorientation angle $\theta < 15^\circ$ can be described by a wall of edge dislocations [39]. These dislocations alter the elastic field around GB and hereby significantly influence solute atoms diffusion and segregation. In Figure 4, the GB segregation of X₁ (Figure 4a) and X₂ (Figure 4b) atoms for the low misorientation angles $\theta = 7.15^\circ$ and $\theta = 9.53^\circ$ respectively, are displayed after projection on the (100) plane.

Therein, Fe atoms in two successive (100) planes ($n$ and $n + 1$) are colored in white and black respectively. X₁ solute atoms are observed in (100) Fe planes, as well as (100) Fe interplanes. X₁ atoms lying in $n, n + \frac{1}{2}$ ($n/n + 1$ interplane) and $n + 1$ (100) Fe planes are indicated by red, orange and yellow spheres, respectively. As for X₂ solute atoms, they only occupy (100) Fe planes. Accordingly, X₂ atoms lying in $n$ and $n + 1$ (100) Fe planes are indicated by red and yellow spheres. In addition, $\langle 0 \pm 10 \rangle$ edge dislocations are spotted by red $\perp$ marks. On the one hand, the vast majority of X₁ atoms segregate in interstitial position between two (100) planes of Fe atoms as demonstrated in Figure 4a (orange spheres in $n + \frac{1}{2}$ plane). This observation is consistent with previous ab initio calculations on the segregation of P atoms at symmetric tilt grain boundaries in α-iron [17,47]. Only X₁ atoms nesting at the dislocation core are situated in (100) Fe planes (red and yellow spheres). On the other hand, X₂ atoms only segregate in (100) Fe planes (red and yellow spheres in Figure 4b).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Solute atoms segregation at 〈100〉 symmetric tilt LAGBs consisting of a wall of 〈100〉 edge dislocations (red ⊢ marks), as provided by QA simulations (*t* = 400). (a) X₁ (small) atoms segregation for a misorientation angle θ = 7.15°. (b) X₂ (larger) atoms for θ = 9.53°. Visualization via OVITO [48] after projection on the (100) plane. Black and white spheres correspond to Fe atoms in two successive (100) planes (*n* and *n* + 1), while red, orange and yellow spheres correspond to solute (X) atoms in three successive (100) planes (*n*, *n* + ½ and *n* + 1). The size of atoms follows the ordering of solute atoms radius (*R*(X₁) < *R*(X₂) < *R*(X₃)).</figcaption></figure>

In both cases, solute atoms segregate around the GB dislocations in interstitial positions, which leads to the formation of so-called Cottrell atmospheres [49] (see Figure 5a for X₁ atoms and (c) for X₂ atoms). The latter is widely used to describe carbon atoms segregation at GB in iron [50], as it roots dislocation pinning and dynamic strain aging in steels [51]. This peculiar distribution of solute atoms reflects the stress generated by edge dislocations at GB. This is compressive above the 〈010〉 slip plane (blue shades for atoms on the left of dislocations), and tensile below this plane (red shades on the right of dislocations). With this, solute atoms (orange in Figure 5a,c, green dots in Figure 5b,d) are mainly distributed close to the dislocation core, or at its border, between the dilated and compressed region. Therefore, the periodic pattern for solute atoms segregation stems from the periodic positioning of edge dislocations forming the wall of dislocations at LAGBs.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Formation of Cottrell atmospheres at (100) edge dislocations (red ⊥ marks) after segregation of solute atoms at (100) symmetric tilt LAGBs. (a,b) X₁ (small) atoms segregation for a misorientation angle θ = 7.15°. (c,d) X₂ (larger) atoms for θ = 9.53°. (a,c) solute atoms distribution (orange) and Fe atoms (transparent). (b,d) Volume per atom variation ΔV/V (Voronoi analysis). Red—dilatation, blue—compression. X atoms are spotted by green dots in (b,d).</figcaption></figure>

Complementary information can be deduced from the strain field in the (100) and (010) directions at the GB as provided by OVITO (finite strain theory). The strain field is respectively shown in Figure 6a,c for X₁ atoms, and Figure 6b,d for X₂ atoms. Albeit X₁ and X₂ atoms (green dots) displaying a rather similar segregation pattern in the close vicinity of the dislocation core (⊥ mark) where the dilatation is strong (see Figure 5b,d), they respond differently to the different components of the strain field tensor. First, X₂
atoms segregate where $\sigma_{XX} < 0$ is minimal (blue atoms emphasized by black dashed line in Figure 6b, and $\sigma_{YY} > 0$ is maximal (red atoms emphasized by black dashed line in Figure 6d). Contrariwise, $X_1$ segregates preferentially at the limit between compressed and dilated areas in the x and y directions (black dashed line in Figure 6a,c). These departing behaviors between $X_1$ and $X_2$ atoms is latched to the nature of the interstitial position they are prone to occupy. $X_1$ atoms are small enough to be interspersed in <111> bcc directions. In this regard, this distribution of small solute atoms ($X_1$) is reminiscent of crowdions formed by P atoms in irradiated α-Fe notably [52]. As for $X_2$ atoms, their radius is too big to be located in the same directions as $X_1$ atoms. They rather occupy octahedral interstitial sites. The pivotal factor for their preferential location is thus the tension of the bcc lattice in one direction, so the volume of octahedral sites in the corresponding direction is increased.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Strain field in the vicinity of one <100> edge dislocations after segregation of $X_1$ (a,c) and $X_2$ (b,d) atoms at a low angle GB ($\theta = 7.15^\circ$ for $X_1$ and $\theta = 9.53^\circ$ for $X_2$, as provided by QA simulations ($t = 400$). (a,b) Strain field $\sigma_{XX}$ in the <100> direction (red: tension in the x direction, blue: compression in the x direction). (c,d) Strain field $\sigma_{YY}$ in the <010> direction for $X_2$ (red—tension in the y direction, blue—compression in the y direction). X atoms are spotted by green dots for the sake of clarity. Black dashed lines: guide for the eye for the preferential area of segregation.</figcaption></figure>

### 3.1.2. High Angle GB

The effect of atomic radius on the position of solute atoms at high angle <100> symmetric tilt GB is now studied. Two representative misorientation angles were selected: one special GB, $\Sigma5$ (310) ($\theta = 36.95^\circ$) in Figures 7 and 8, and one general GB, $\Sigma29$ (730) ($\theta = 46.40^\circ$) in Figures 9 and 10. The structure of HAGBs is usually described using the structural unit (SU) representation [53,54]. In these figures, SUs are highlighted by solid black (plane $n$) and gray (plane $n+1$) lines. The virgin GB structure of $\Sigma5$ (310) is characterized by the structural units $|B.B|$ [55], while the structure of $\Sigma29$ (730) consists of a sequence of $|BC.BC|$ SUs.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Segregation of solute X₁ (a–c) X₂ (d–f) and X₃ (g–i) atoms at Σ5 (310) (θ = 36.87°) as provided by QA simulations (t = 400). B structural units are materialized by black (n plane) and gray (n + 1 plane) lines. (a,d,g) Fe and X atoms positions. (b,e,h) Fe₉X clusters. (c,f,i) Variation of X solute atom position within B unit, or equivalently Fe₉X cluster. The depleted (over-compressed) and segregated (less compressed) areas are delineated by dashed black lines, and referred to as zones 1 and 2 in (a,d,g).</figcaption></figure>

First, the segregation of X₁, X₂ and X₃ atoms (atomic radius R<sub>X1</sub> < R<sub>X2</sub> < R<sub>X3</sub>, comparable to P, Sb and Sn ionic radii, respectively) at Σ5 (310) is addressed, using the same color coding as that for LAGBs. At the GB, each B SU hosts one solute atom only, which lies on the same (100) α-Fe plane as the first nearest Fe neighbors, disregarding its atomic radius. Now, B units corresponds to the projection on the (100) plane of a three dimensional capped trigonal prisms structure. With this, Fe₉X clusters are formed upon hosting solute X atoms, as shown in Figure 7b,e,h. This typical structure for Σ5 (310) GB was notably witnessed for phosphorus (same atomic radius as X₁) and boron segregation in α-iron in previous MD studies [47]. However, while the structure of Fe₉(X₁) clusters closely resembles that of Fe₉P and Fe₉B (small atomic radius), a significant deviation can be observed for the Fe₉(X₂) and Fe₉(X₃) cluster structure as displayed in Figure 7c,f,i. Compared to X₁ atoms, the equilibrium position of X₂ and X₃ atoms is shifted to a more off-centered position within their hosting capped trigonal prism. The deformation of the Fe₉X cluster is maximal for X₃ atoms (largest atomic radius) in Figure 7i. This migration of X₂ and X₃ solute atom within their respective Fe₉X clusters is accompanied by a deformation of the capped trigonal prism. This is highlighted by the map of displacement vectors of Fe atoms in Figure 7c,f,i, where red arrows stand for Fe atoms displacement (amplified by a factor of 10 in the figure), between the equilibrium structure of virgin GB, and the GB structure after segregation. For small solute atoms (X₁), the capped trigonal prism is
roughly unaffected by segregation, once again in agreement with [47] for atoms with an atomic radius close to P. For bigger atoms, however (X₂ and X₃), noticeable displacements of Fe atom can be observed in the (001) plane. The general deformation trend is the dilatation of capped trigonal prism for X₂ atoms (intermediate atomic radius), versus shear displacements of Fe atoms for X₃ atoms (largest atomic radius). This observation may suggest that above a certain threshold for solute atoms radius, the SUs accommodate differently the introduction of solute atoms of different sizes.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Influence of the volume per atom (Voronoi analysis) on the segregation of X₁ (a) X₂ (b) and X₃ (c) solute atoms at Σ5 (310) (θ = 36.87°). X₁ atoms are spotted by green dots. The depleted (over-compressed) and segregated (less compressed) areas are delineated by dashed black lines, and referred to as zones 1 and 2 in (b,c).</figcaption></figure>
Interesting features can also be observed in the vicinity of GBs. For small solute atoms ($X_1$), a significant segregation can be observed up to five (010) atomic planes on both sides of the GB (see Figure 7c). This observation is consistent with previous molecular static simulations, where a substantially lower segregation energy for C atoms was obtained for distances up to 5 Å away from Σ5 (310), and 10 Å away from Σ29 (730) [22]. A similar pattern can be envisioned for $X_2$ atoms located far from the GB (zone 2 in Figures 7d and 8b). However, a narrow depleted zone is present at the immediate border of the GB (zone 1 in Figures 7d and 8b). The same observation applies for the largest solute atoms ($X_3$) in Figures 7g and 8c, where the thickness of the depleted region (zone 1) is even increased.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Segregation of $X_1$ atoms at Σ29 (730) ($\theta = 46.40^{\circ}$) as provided by QA simulations ($t = 400$). (a) Atoms positions. (b) Strain field of Fe atoms in the 〈100〉 direction. $X_1$ atoms are spotted by green dots. B and C structural units are materialized by black ($n$ plane) and gray ($n + 1$ plane) lines.</figcaption></figure>

A first interpretation, relying on the volume per atom, can be proposed, based on Figure 8. Generally speaking, a marked dilatation of the host Fe structure at the GB (atoms in red) goes along with an extended compression area bordering the GB (atoms in blue belonging to zones 1 and 2). The dilatation zone reflects the presence of the capped trigonal prism (B units), which are prone to host solute atoms disregarding their size. Now, this dilatation at the GB is balanced by the compression of the bcc Fe lattice in the abutting region. However, the connection between the position of solute atoms and the elastic field in the compressed region depends on the atomic radius of solute atoms in the following manner: small ($X_1$) atoms are present in the compressive region, up to a 5–10 Å distance from the GB (green dots in Figure 8a). Intermediate size atoms ($X_2$) are absent in the significantly compressed area (zone 1), but present in the moderately compressed area (zone 2). Finally, large atoms ($X_3$) are absent almost everywhere in the compressed regions, but some can be spotted in further regions (zone 2 in Figure 8c). Based on this observation, we suggest that small ($X_1$) atoms might be less sensible to the compression of the matrix due to their size, hence their recurrent segregation in this area, in spite of the compression.
On the contrary, X₂ and X₃ atoms remain in octahedral interstitial sites, and may be too big to segregate in regions where the host structure is sufficiently compressed.

In more details, the elastic field at SUs is not rigorously identical in the three cases. Indeed, the dilatation at B SUs is more intense for X₂ atoms (~10% in Figure 8b) than for X₁ and X₃ atoms (~7% in Figure 8a,c). This observation is consistent with the nature and the amplitude of the deformation of the capped trigonal prism depicted in Figure 7c,f,i, depending on the atomic radius of the solute atoms: a negligible deformation for X₁ atoms and a shear deformation for X₃ atoms, versus a planar dilatation for X₂ atoms. One consequence of this might be the over compression of the bcc structure in region 1 for X₂ segregation in Figure 8b, as entailed by the over dilatation at the GB, and a reduced compression in region 2. We surmise that the over compression in region 1 alternatively precludes the segregation of X₂ atoms in this area, and yet fosters their segregation in region 2 (see Figures 7b and 8b).

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Segregation of X₂ atoms at Σ29 (730) (θ = 46.40°) as provided by QA simulations (t = 400). (a) Atoms positions. (b) Strain field of Fe atoms in the 〈100〉 direction. (c) Strain field of Fe atoms in the 〈010〉 direction. X₂ atoms are spotted by green dots. B and C structural units are materialized by black (n plane) and gray (n + 1 plane) lines.</figcaption></figure>
One last remark touches upon the periodic distribution of solute atoms in the 〈100〉 direction near the Σ5 (310) GB. This periodicity stems from the structural periodicity of the sequence of B units at the GB. If this conclusion is obvious for solute atoms belonging to the Fe₉X clusters, we believe that it remains valid in the vicinity of the GB as well. In this case, the structural periodicity at the GB feeds through to the periodicity of the strain deformation tensor close to the GB.

In order to examine the connection between the structural periodicity of the GB and the positions of solute atoms, the segregation of X₁ and X₂ atoms at Σ29 (730) (θ = 46.40°) was also prospected in Figures 9 and 10. In this case, the majority of the conclusions formulated for Σ5 (310) remain valid for this GB (segregation at and near GBs, preferential interstitial positions for solute atoms, influence of volumetric strain and atomic radius, etc.). Now, notwithstanding an identical position of solute atoms within B units for both Σ5 (310) and Σ29 (730) GB, a different segregation tendency for X₁ and X₂ atoms within C units can be observed for Σ29 (730). In this case, four X₁ atoms are located within C units, versus only one X₂ atom. Here again, this phenomenon is linked to the smaller atomic radius of X₁ atoms, which allows the interspersed segregation in less dilated, or even compressed areas of the host structure.

One step further, the Σ29 (730) GB presents a more complicated segregation pattern than Σ5 (310). This likely stems from the presence of two SUs, which generates a more complex stress field in the vicinity of the GB, compared to a GB consisting of B units only. Indeed, it clearly appears in Figures 9b and 10b,c that the periodic segregation of solute atoms reflects the periodic modulation of the strain tensor. These variations of the strain field consists of the alternation of compressed (tensed) areas in the x (y) direction, with compressed (tensed) areas in the y (x) direction. The spatial periodicity of this variation of the field is precisely the length of one BC SU in the x direction. In detail, X₁ atoms are preferentially distributed in regions where σ<sub>XX</sub> (Figure 9b) and σ<sub>YY</sub> (not shown) switch signs, while X₂ atoms circumscribe areas where σ<sub>XX</sub> ≫ 0 (in red in Figure 10b) and σ<sub>YY</sub> ≪ 0 (in blue in Figure 10b). This is consistent with the observations of LAGBs.

## 4. Conclusions

In this study, we applied a new atomistic model based on the quasiparticle approach to explore the relationship between the local grain boundary structure and the size of the segregation solute atom in the α-iron crystal. Three types of solute atoms were considered, with three different atomic radii: X₁ atoms with a radius much smaller than Fe atom (R = 0.68R(Fe)) corresponding to phosphorus (P), X₂ atoms with an atomic radius close to Fe (R = 0.93R(Fe)) corresponding to antimony (Sb), and X₃ with a larger atomic size than Fe (R = 1.08R(Fe)) corresponding to tin (Sn). Two cases were investigated: segregation at LAGBs (θ < 15°) on the one hand, and at HAGBs Σ5 (310) (θ = 36.95°) and Σ29 (730) on the other hand. It was evidenced that all three sorts of atoms segregate at LAGBs in interstitial positions and generate Cottrell atmospheres around GB dislocations. X₁ and X₂ atoms respond differently to the different components of the strain field tensor. It was indeed shown that X₁ (small) solute atoms segregate preferentially at the limit between compressed and dilated areas in the x and y directions, whereas X₂ atoms are rather located where σ<sub>XX</sub> is minimal and σ<sub>YY</sub> is maximal (or the opposite). In the case of Σ5 (310) (θ = 36.95°) HAGBs, the three types of solute atoms form Fe₉X clusters in B units, with a capped trigonal prism structure. Upon increasing the size of solute atoms, a dilatation of this capped prism was observed. In detail, X₂ (larger) solute atoms induce an homogeneous dilatation of the hosting prism, while X₃ (largest) atoms entail a shear displacements of Fe atoms. One last noteworthy point regarding HAGBs touches upon the presence of a depleted zone at the immediate border of the GB. The width of this area rises with the atomic radius of solute species. This peculiar distribution of solute atoms was shown to nicely reflect the periodic amplitude variations of the elastic field around the GB. In the case of Σ29 (730), which hosts a series of |BC.BC| SUs, a similar segregation trend within B SUs as in Σ5 (310) was observed. However, the position of X₁ and X₂ atoms
within C SUs is different, insofar as four X₁ atoms are encountered therein, as opposed to a single X₂ (or X₃) atom. In details, X₁ atoms are preferentially distributed in regions where σ<sub>XX</sub> and σ<sub>YY</sub> switch signs, while X₂ atoms are located in the areas where σ<sub>XX</sub> ≪ 0 and σ<sub>YY</sub> ≫ 0 (or the opposite).

**Author Contributions:** A.V., G.D., R.P. and H.Z. developed the QA model. A.V. and G.D. performed the simulations. All authors participated in the redaction of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the French National Research Agency (ANR), contract C-TRAM ANR18-CE92-0021.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data sets generated in the current study are available from the corresponding authors upon reasonable request.

**Acknowledgments:** Part of this work was performed using computing resources of CRIANN (Nornandy, France) where simulations were performed as Project No. 2012008.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. McLean, D.; Maradudin, A. Grain boundaries in metals. *Phys. Today* **1958**, *11*, 35. [CrossRef]

2. Balluffi, R.; Sutton, A. Why should we be interested in the atomic structure of interfaces? In *Materials Science Forum*; Trans Tech Publications Ltd.: Freienbach, Switzerland, 1996; Volume 207, pp. 1–12.

3. Watanabe, T.; Tsurekawa, S. The control of brittleness and development of desirable mechanical properties in polycrystalline systems by grain boundary engineering. *Acta Mater.* **1999**, *47*, 4171–4185. [CrossRef]

4. Lejček, P.; Hofmann, S.; Paidar, V. Solute segregation and classification of [100] tilt grain boundaries in α-iron: Consequences for grain boundary engineering. *Acta Mater.* **2003**, *51*, 3951–3963. [CrossRef]

5. Randle, V. Twinning-related grain boundary engineering. *Acta Mater.* **2004**, *52*, 4067–4081. [CrossRef]

6. Krakauer, B.; Seidman, D. Subnanometer scale study of segregation at grain boundaries in an Fe (Si) alloy. *Acta Mater.* **1998**, *46*, 6145–6161. [CrossRef]

7. Raabe, D.; Herbig, M.; Sandlöbes, S.; Li, Y.; Tytko, D.; Kuzmina, M.; Ponge, D.; Choi, P.P. Grain boundary segregation engineering in metallic alloys: A pathway to the design of interfaces. *Curr. Opin. Solid State Mater. Sci.* **2014**, *18*, 253–261. [CrossRef]

8. Xing, W.; Kalindi, A.R.; Amram, D.; Schuh, C.A. Solute interaction effects on grain boundary segregation in ternary alloys. *Acta Mater.* **2018**, *161*, 285–294. [CrossRef]

9. Divinski, S.V.; Edelhoff, H.; Prokofjev, S. Diffusion and segregation of silver in copper Σ 5 (310) grain boundary. *Phys. Rev. B* **2012**, *85*, 144104. [CrossRef]

10. Herbig, M.; Raabe, D.; Li, Y.; Choi, P.; Zaefferer, S.; Goto, S. Atomic-scale quantification of grain boundary segregation in nanocrystalline material. *Phys. Rev. Lett.* **2014**, *112*, 126103. [CrossRef] [PubMed]

11. Yamaguchi, M. First-principles study on the grain boundary embrittlement of metals by solute segregation: Part I. iron (Fe)-solute (B, C, P, and S) systems. *Metall. Mater. Trans. A* **2011**, *42*, 319–329. [CrossRef]

12. Scheiber, D.; Pippan, R.; Puschnig, P.; Romaner, L. Ab initio calculations of grain boundaries in bcc metals. *Model. Simul. Mater. Sci. Eng.* **2016**, *24*, 035013. [CrossRef]

13. Yamaguchi, M.; Ebihara, K.I.; Itakura, M.; Kadoyoshi, T.; Suzudo, T.; Kaburaki, H. First-principles study on the grain boundary embrittlement of metals by solute segregation: Part II. Metal (Fe, Al, Cu)-hydrogen (H) systems. *Metall. Mater. Trans. A* **2011**, *42*, 330–339. [CrossRef]

14. Tahir, A.; Janisch, R.; Hartmaier, A. Ab initio calculation of traction separation laws for a grain boundary in molybdenum with segregated C impurities. *Model. Simul. Mater. Sci. Eng.* **2013**, *21*, 075005. [CrossRef]

15. Razumovskiy, V.I.; Ruban, A.V.; Razumovskii, I.; Lozovoi, A.; Butrim, V.; Vekilov, Y.K. The effect of alloying elements on grain boundary and bulk cohesion in aluminum alloys: An ab initio study. *Scr. Mater.* **2011**, *65*, 926–929. [CrossRef]

16. Razumovskiy, V.I.; Lozovoi, A.; Razumovskii, I. First-principles-aided design of a new Ni-base superalloy: Influence of transition metal alloying elements on grain boundary and bulk cohesion. *Acta Mater.* **2015**, *82*, 369–377. [CrossRef]

17. Yamaguchi, M.; Nishiyama, Y.; Kaburaki, H. Decohesion of iron grain boundaries by sulfur or phosphorous segregation: First-principles calculations. *Phys. Rev. B* **2007**, *76*, 035418. [CrossRef]

18. Wang, J.; Janisch, R.; Madsen, G.K.; Drautz, R. First-principles study of carbon segregation in bcc iron symmetrical tilt grain boundaries. *Acta Mater.* **2016**, *115*, 259–268. [CrossRef]
19. Hu, Y.J.; Wang, Y.; Wang, W.Y.; Darling, K.A.; Kecskes, L.J.; Liu, Z.K. Solute effects on the Σ3 111 [11-0] tilt grain boundary in BCC Fe: Grain boundary segregation, stability, and embrittlement. *Comput. Mater. Sci.* **2020**, *171*, 109271. [CrossRef]

20. Daw, M.S.; Baskes, M.I. Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals. *Phys. Rev. Lett.* **1983**, *50*, 1285. [CrossRef]

21. Barrera, O.; Bombac, D.; Chen, Y.; Daff, T.; Galindo-Nava, E.; Gong, P.; Haley, D.; Horton, R.; Katzarov, I.; Kermode, J.; others. Understanding and mitigating hydrogen embrittlement of steels: A review of experimental, modeling and design progress from atomistic to continuum. *J. Mater. Sci.* **2018**, *53*, 6251–6290. [CrossRef]

22. Rhodes, N.; Tschopp, M.; Solanki, K. Quantifying the energetics and length scales of carbon segregation to α-Fe symmetric tilt grain boundaries using atomistic simulations. *Model. Simul. Mater. Sci. Eng.* **2013**, *21*, 035009. [CrossRef]

23. Wagih, M.; Larsen, P.M.; Schuh, C.A. Learning grain boundary segregation energy spectra in polycrystals. *Nat. Commun.* **2020**, *11*, 1–9. [CrossRef]

24. Elder, K.; Katakowski, M.; Haataja, M.; Grant, M. Modeling elasticity in crystal growth. *Phys. Rev. Lett.* **2002**, *88*, 245701. [CrossRef] [PubMed]

25. Elder, K.; Grant, M. Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals. *Phys. Rev. E* **2004**, *70*, 051605. [CrossRef]

26. Elder, K.R.; Provatas, N.; Berry, J.; Stefanovic, P.; Grant, M. Phase-field crystal modeling and classical density functional theory of freezing. *Phys. Rev. B* **2007**, *75*, 064107. [CrossRef]

27. Greenwood, M.; Ofori-Opoku, N.; Rottler, J.; Provatas, N. Modeling structural transformations in binary alloys with phase field crystals. *Phys. Rev. B* **2011**, *84*, 064104. [CrossRef]

28. Smith, N.; Provatas, N. Generalization of the binary structural phase field crystal model. *Phys. Rev. Mater.* **2017**, *1*, 053407. [CrossRef]

29. Mellenthin, J.; Karma, A.; Plapp, M. Phase-field crystal study of grain-boundary premelting. *Phys. Rev. B* **2008**, *78*, 184110. [CrossRef]

30. Khachaturian, A. Ordering in substitutional and interstitial solid solutions. *Prog. Mater. Sci.* **1978**, *22*, 1–150. [CrossRef]

31. Jin, Y.M.; Khachaturyan, A.G. Atomic density function theory and modeling of microstructure evolution at the atomic scale. *J. Appl. Phys.* **2006**, *100*, 013519. [CrossRef]

32. Kapikranian, O.; Zapolsky, H.; Domain, C.; Patte, R.; Pareige, C.; Radiguet, B.; Pareige, P. Atomic structure of grain boundaries in iron modeled using the atomic density function. *Phys. Rev. B* **2014**, *89*, 014111. [CrossRef]

33. Kapikranian, O.; Zapolsky, H.; Patte, R.; Pareige, C.; Radiguet, B.; Pareige, P. Point defect absorption by grain boundaries in α-iron by atomic density function modeling. *Phys. Rev. B* **2015**, *92*, 224106. [CrossRef]

34. Vaugeois, A. Modélisation de l'influence de la Structure des Joints de Grains sur les Phénomènes de Ségrégation. Ph.D. Thesis, Normandie Université, Caen, France, 2017.

35. Lavrskyi, M.; Zapolsky, H.; Khachaturyan, A.G. Quasiparticle approach to diffusional atomic scale self-assembly of complex structures: From disorder to complex crystals and double-helix polymers. *NPJ Comput. Mater.* **2016**, *2*, 1–9. [CrossRef]

36. Demange, G.; Lavrskyi, M.; Chen, K.; Chen, X.; Wang, Z.; Patte, R.; Zapolsky, H. Fcc-> bcc phase transition kinetics in an immiscible binary system: Atomistic evidence of the twinning mechanism of transformation. *arXiv* **2021**, arXiv:2103.12384.

37. Bollmann, W. *Crystal Defects and Crystalline Interfaces*; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012.

38. Sutton, A.; Balluffi, R. *Interfaces in Crystalline Solids, Clarendon*; Oxford University Press: Oxford, UK, 1995.

39. Read, W.T.; Shockley, W. Dislocation models of crystal grain boundaries. *Phys. Rev.* **1950**, *78*, 275. [CrossRef]

40. Jaatinen, A.; Achim, C.; Elder, K.; Ala-Nissila, T. Thermodynamics of bcc metals in phase-field-crystal models. *Phys. Rev. E* **2009**, *80*, 031602. [CrossRef]

41. Demange, G.; Chamaillard, M.; Zapolsky, H.; Lavrskyi, M.; Vaugeois, A.; Luneville, L.; Simeone, D.; Patte, R. Generalization of the fourier-spectral eyre scheme for the phase-field equations: Application to self-assembly dynamics in materials. *Comput. Mater. Sci.* **2018**, *144*, 11–22. [CrossRef]

42. Müller, M.; Erhart, P.; Albe, K. Analytic bond-order potential for bcc and fcc iron—Comparison with established embedded-atom method potentials. *J. Phys. Condens. Matter* **2007**, *19*, 326220. [CrossRef]

43. Simonelli, G.; Pasianot, R.; Savino, E. Embedded-atom-method interatomic potentials for bcc-iron. *MRS Online Proc. Libr. Arch.* **1992**, *291*. [CrossRef]

44. Ackland, G.; Bacon, D.; Calder, A.; Harry, T. Computer simulation of point defect properties in dilute Fe-Cu alloy using a many-body interatomic potential. *Philos. Mag.* **A** **1997**, *75*, 713–732. [CrossRef]

45. Mendelev, M.; Han, S.; Srolovitz, D.; Ackland, G.; Sun, D.; Asta, M. Development of new interatomic potentials appropriate for crystalline and liquid iron. *Philos. Mag.* **2003**, *83*, 3977–3994. [CrossRef]

46. Dudarev, S.; Derlet, P. A 'magnetic' interatomic potential for molecular dynamics simulations. *J. Phys. Condens. Matter* **2005**, *17*, 7097. [CrossRef]

47. Hashimoto, M.; Ishida, Y.; Yamamoto, R.; Doyama, M. Atomistic studies of grain boundary segregation in Fe-P and Fe-B alloys—I. Atomic structure and stress distribution. *Acta Metall.* **1984**, *32*, 1–11. [CrossRef]

48. Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO—the Open Visualization Tool. *Model. Simul. Mater. Sci. Eng.* **2009**, *18*, 015012. [CrossRef]

49. Cottrell, A.H.; Bilby, B.A. Dislocation theory of yielding and strain ageing of iron. *Proc. Phys. Soc. Sect. A* **1949**, *62*, 49. [CrossRef]
50. Wilde, J.; Cerezo, A.; Smith, G. Three-dimensional atomic-scale mapping of a Cottrell atmosphere around a dislocation in iron. *Scr. Mater.* **2000**, *43*, 39–48. [CrossRef]

51. Portavoce, A.; Tréglia, G. Theoretical investigation of Cottrell atmosphere in silicon. *Acta Mater.* **2014**, *65*, 1–9. [CrossRef]

52. Nguyen-Manh, D.; Horsfield, A.; Dudarev, S. Self-interstitial atom defects in bcc transition metals: Group-specific trends. *Phys. Rev. B* **2006**, *73*, 020101. [CrossRef]

53. Bishop, G.H.; Chalmers, B. A coincidence—Ledge—Dislocation description of grain boundaries. *Scr. Metall.* **1968**, *2*, 133–139. [CrossRef]

54. Sutton, A.P.; Vitek, V. On the structure of tilt grain boundaries in cubic metals I. Symmetrical tilt boundaries. *Philos. Trans. R. Soc. Lond. Ser. Math. Phys. Sci.* **1983**, *309*, 1–36.

55. Tschopp, M.A.; Solanki, K.; Gao, F.; Sun, X.; Khaleel, M.A.; Horstemeyer, M. Probing grain boundary sink strength at the nanoscale: Energetics and length scales of vacancy and interstitial absorption by grain boundaries in α-Fe. *Phys. Rev. B* **2012**, *85*, 064108. [CrossRef]