Ab Initio Investigation Of Phosphorus And Hydrogen Co-segregation And Embrittlement In {\alpha}-fe Twin Boundaries

TL;DR

This paper introduces a new statistical physics model to analyze how phosphorus and hydrogen co-segregate at grain boundaries in alpha-iron and how this affects embrittlement, using atomistic calculations for parameters.

Contribution

The study develops a general, efficient model for solute segregation at grain boundaries that incorporates multiple species, sites, and entropy, applied specifically to P-H co-segregation in alpha-Fe.

Findings

P-H interactions increase phosphorus segregation at grain boundaries.

Co-segregation of P and H leads to greater embrittlement than individual segregation.

The model accurately predicts segregation behavior using parameters from density-functional theory.

Abstract

We propose a new statistical physics model to study equilibrium solute segregation at grain boundaries and the resulting embrittlement effect. This low-temperature expansion model is general and efficient, and its parameters can be obtained from atomistic calculations. It is possible to take into account multiple species, multiple segregation sites with different segregation free energies, account for configurational entropy, grain radius and site competition between solutes. As an example, the model is then applied to the study of phosphorus and hydrogen co-segregation at $\Sigma3$ $109.5\deg[0\bar{1}1]\{111\}$ twin boundaries in $\alpha$-Fe, using energetic parameters from density-functional theory calculations. We show that P--H interactions may lead to increased P segregation at grain boundaries and cause additional embrittlement compared to the case where P and H are considered separately.

Figures

Equations

$$E_{zpe}^H=\frac{h}{2}\sum _i{\nu }_i=\frac{\hslash }{2}\sum _i\sqrt{\frac{k_i}{m_H}},$$

$$F_{vib}^H=k_{B}T\sum _i\ln \left(2\sinh \left(\frac{h{\nu }_i}{2k_{B}T}\right)\right),$$

$$E^b\left(\mathrm{PH}\right)$$

$$-E\left[54\mathrm{Fe}\right]-E\left[53\mathrm{Fe}+\mathrm{P}+\mathrm{H}\right],$$

$$E^b\left(\mathrm{PH}\right)$$

$$-2E\left[54\mathrm{Fe}\right]-E\left[\left(n-1\right)\mathrm{Fe}+\mathrm{P}+\mathrm{H},d\right],$$

$$\mathcal{A}$$

$$=-k_{B}T\ln \left(\sum _i{G}_i\exp \left(\frac{-E_i+{\sum }_{\alpha }{n}_i^{\alpha }{\mu }_{\alpha }}{k_{B}T}\right)\right),$$

$$\mathcal{A}={\mathcal{A}}_0-$$

$$k_{B}T\ln \left(1+\sum _{i\ne 0}{G}_i\exp \left(\frac{E_0-E_i+{\sum }_{\alpha }\delta n_i^{\alpha }{\mu }_{\alpha }}{k_{B}T}\right)\right),$$

$$\frac{\mathcal{A}}{N}=\frac{{\mathcal{A}}_0}{N}-k_{B}T\sum _{i\ne 0}{g}_i\exp \left(\frac{E_i^b+{\sum }_{\alpha }\delta n_i^{\alpha }{\mu }_{\alpha }}{k_{B}T}\right),$$

$$\left[\alpha \right]=-\frac{1}{N}\frac{\partial \mathcal{A}}{\partial {\mu }_{\alpha }}=\frac{n_0^{\alpha }}{N}+\sum _{i\ne 0}\delta n_i^{\alpha }{g}_i\exp \left(\frac{E_i^b}{k_{B}T}\right)\prod _{\alpha }{X}_{\alpha }^{\delta n_i^{\alpha }}.$$

$$\gamma =\frac{N_{\varphi }}{N_b}=\frac{4\pi r^2{V}_{at}}{\frac{4}{3}\pi r^3{S}_{at}},$$

$$g_i=\frac{\gamma }{1+\gamma }{g}_i^{\varphi }=\frac{1}{1+\gamma }{g}_i^b.$$

$$E_0^{\varphi }=N{E}_m+N_{\varphi }{\rho }_{\varphi }{\gamma }_{\varphi },$$

$$W_{sep}^0=\frac{E_0^{fs}-E_0^{gb}}{N_{gb}}=\rho \left(2{\gamma }_{fs}-{\gamma }_{gb}\right).$$

$$W_{sep}=W_{sep}^0-\left[X\right]EP,$$

$$GE{P}_b=\sum _{i\in GB}\left(E_{i,fs}^b-E_{i,gb}^b\right)g_i^{gb}\exp \left(\beta {\mathcal{A}}_{i,gb}^b\right),$$

$$GE{P}_f=\sum _{\alpha }\left[\alpha \right]\left({\mu }_{\alpha }^{fs}-{\mu }_{\alpha }^{gb}\right)$$

$$+\sum _{i\in GB}{g}_i^{\varphi }\left(2E_{i,fs}^b\exp \left(\beta {\mathcal{A}}_{i,fs}^b\right)-E_{i,gb}^b\exp \left(\beta {\mathcal{A}}_{i,gb}^b\right)\right)$$

$$+\sum _{i\in bulk}\frac{g_i^b{E}_{i,b}^b}{\gamma }\left((1-\gamma )\exp \left(\beta {\mathcal{A}}_{i,fs}^b\right)-\exp \left(\beta {\mathcal{A}}_{i,gb}^b\right)\right),$$

$$E^f\left(vac\right)=E\left[53\mathrm{Fe}+vac\right]-\frac{53}{54}E\left[54\mathrm{Fe}\right],$$

$$\left[H\right]=\frac{1}{1+\gamma }6Y+\frac{\gamma }{1+\gamma }Y{S}_H\simeq Y\left(6+\gamma S_H\right).$$

$${\theta }_H=Y{S}_H=\frac{\left[H\right]S_H}{6+\gamma S_H}=\frac{\left[H\right]}{\frac{6}{S_H(T)}+\frac{\sqrt{3}}{2r}}.$$

$$\ln \left[H\right]=\ln \left(\frac{6{\theta }_H}{S_H}\right)+\ln \left(1+\frac{\sqrt{3}{S}_H}{12r}\right),$$

$$G_1=N,$$

$$G_2=\frac{N\left(N-1\right)}{2!},$$

$$G_3=\frac{N\left(N-1\right)\left(N-2\right)}{3!},$$

$$G_i=\frac{N\left(N-1\right)\cdots \left(N-i+1\right)}{i!},$$

$$\mathcal{A}$$

$$={\mathcal{A}}_0-k_{B}TN\sum _{i\ne 0}\frac{(-1{)}^{i-1}}{i}{X}^i$$

$$={\mathcal{A}}_0-k_{B}TN\ln \left(1+X\right),$$

$${\left[m\right]}^{\varphi }$$

$${\left[\alpha \right]}^{\varphi }$$

$$\frac{{\left[\alpha \right]}^{\varphi }}{{\left[m\right]}^{\varphi }}=\exp \left(\frac{E_{\alpha }^b+{\mu }_{\alpha }-{\mu }_m}{k_{B}T}\right).$$

Full text

1 Introduction

Grain boundaries (GBs) are of prime importance in the study of materials because they form an interconnected network that spreads over the whole material and directly affects its mechanical properties. Hence, the study of GBs has drawn a lot of attention for more than 60 years [1, 2, 3]. Nowadays, a certain level of understanding of GB thermodynamics and kinetics is reached, such that microstructure optimization via GB engineering approaches are emerging [3, 4]. One of the most studied feature of GBs is their interaction with solutes from the bulk material, because segregation phenomena lead to large and very localized chemical inhomogeneities [5], and atomic bonding at GB planes are found very sensitive to the presence of solutes [6].

There are many difficulties in investigating GB segregation experimentally: ideally, we would need a technique that is able to probe very localized regions of the material while making sure that the chemical distribution is not altered by the measurement, and this technique should be able to give accurate chemical information–sometimes for very low concentrations–and crystallographic information, as well as being fast enough to produce sufficient statistics over a large number of GBs. Consequently, the study of GB segregation usually resorts to a combination of experimental techniques [7, 8] and results depend on the annealing temperature, bulk chemical composition, chemical environment and GB structure. The task becomes even more complicated when studying co-segregation effects, i.e. the simultaneous segregation of at least two chemical species which occurs all the time in real-life materials. There are several mechanisms at work: solute-GB interaction, solute-solute interaction at the GB and in the bulk, site competition, concentration ratio between solutes and kinetic properties of each solute [9]. For instance in Fe-P-C alloys, it was noted that intergranular P segregation decreases with increasing C content but the reason why this happens is not clear [10, 8]. If Cr is added to this alloy it seems that P segregation increases only when C concentration is above a certain level [10].

The pioneering work of Wu et al. [6] on $\Sigma 3[1\bar{1}0](111)$ GB in $\alpha$-Fe showed that ab initio methods were able to provide some insight on solute properties at GBs, as the authors were able to explain the strengthening by B additions and the embrittlement by P additions in terms of electronic structure at the GB and at free surfaces. These methods are fast-growing, and computation capabilities have increased a lot such that ab initio studies on GBs are now standard. For instance, it is now possible to study the segregation of various solutes (B,C,N,O,P) on various GBs ($\Sigma 3$ and $\Sigma 5$) at various concentrations (1.0, 0.5 and 0.25 mono-layers) in a systematic manner [11]. Moreover, there are some efforts to regroup the large collection of ab initio data computed over the years and extract some trends as a function of solute properties [12, 13, 14]. Most of this huge amount of data is restricted to the segregation of one type of solute, and most of the time one type of segregation site. Yet, there are some studies focusing on more complex segregation phenomena, e.g. strong co-segregation of Mn-C[15] or site exclusion between N and C[16]. These effects were shown using ab initio calculations but a comprehensive analysis using statistical physics tool is still missing.

Of course these ab initio methods have their drawbacks, mostly related to the complexity of properly taking into account temperature effects, and the limited number of atoms that can be included in a simulation box. Regarding the latter, a question was raised about the reliability of ab initio calculations for solutes with small solubility limits [17, 18]. To shorten the introduction, we briefly comment on this question in A.

A review of the various statistical models that have been developed to study GB segregation can be found in Refs. [19, 2]. The most widely used model is the simplest one, the Langmuir-McLean isotherm [20, 1], which was extended to multi-component systems by Guttmann [21, 9]. The limitation of these approaches is that they are derived from mean-field arguments. Even though there are some semi-empirical procedures to derive the mean-field interaction parameters [22], the mean-field model cannot always render the full complexity of the segregation and the connection with ab initio binding energies is not straightforward. There were some attempts to consider explicitly the short-range order in the first nearest-neighbor shell [23] but this model is not easy to generalize to more complicated case. A model for multi-site segregation was proposed [24], but again the extension to multi-species segregation is not straightforward. The embrittlement potency of a solute is defined as the solute segregation energy difference between the GB and the two free surfaces obtained after fracture [25]. As noted earlier [26], this thermodynamic definition is not valid for more complicated cases, namely multi-site segregation or co-segregation with competing embrittlement potencies and this data is now becoming available (e.g. with molecular statics simulations[27]). In a nutshell, it seems that nowadays we have access to a large number of atomic-scale data, but we are missing some efficient and general model to take full advantage of this data.

Let us now focus on the practical case of P and H co-segregation in $\alpha$-Fe. H [28] and P [29] are both known to alter fracture toughness of steels with very small amounts of impurities. Aucouturier studied a Fe-P-H steels and was not able to tell if P simply reduces the critical H concentration for hydrogen-induced cracking, or if GB embrittlement was due to P-H interactions at the GB [30]. According to Kameda et al. the effects of P and H are additive [31] meaning that P and H do not interact at all, or in other words, there is no co-segregation effect. According to Komazazki, P does not influence the H content of GBs but weakens the GB, which translates into increased H embrittlement [32]. This was somehow confirmed by McMahon who emphasizes the dynamic effect of H and its low segregation enthalpy to GBs [33]. Solubilities of H and P are very low in these materials which complicates the precise determination of P-H co-segregation effects. It seems that modern modeling tools can shed light on the interaction between P and H at the GB. $\Sigma 3$ $109.5\degree\left[0\bar{1}1\right]\left\{111\right\}$ twin boundaries are among the most common ones in lath martensite in low carbon steels [34] and they are quite easily modeled with ab initio methods because they are high-angle symmetric tilt GBs with a well-known structure [35]. Also, due to the dense coincidence site lattice, the range of elastic deformation around the GB is expected to be small such that segregation occurs mainly in the grain boundary plane and we can reduce the size of the simulation cell and the number of configurations to explore.

In this paper, we study the co-segregation of P and H at $\Sigma 3$ $109.5\degree\left[0\bar{1}1\right]\left\{111\right\}$ twin boundaries in $\alpha$-Fe. First, we introduce the modeling tools (Sec. 2), namely density functional theory to compute binding and segregation energies, and low-temperature expansions. The latter provide a new type of statistical model to study grain boundary segregation and has the advantage of being general and straightforward to use in conjunction with ab initio calculations. It can take into account competing segregations (fixed solute concentration distributed among bulk and various GBs) and also the effect of grain size. Then we present ab initio calculations of P-H interactions in the bulk, at the GB and at a free surface (Sec. 3). This data is finally used in the low-temperature expansion model to conclude about the effect of the simultaneous presence of H and P on segregation levels and embrittlement potency (Sec. 4). Finally, we discuss our modeling assumptions (Sec. 5).

2 Methods

2.1 Density-functional theory

Density functional theory (DFT) calculations were performed using the PWscf package from the Quantum ESPRESSO software suite [36, 37]. The exchange and correlation energy is computed with the PBE generalized gradient approximation [38]. Hydrogen $\left(1s^{1}\right)$, phosphorus $\left(\left[\mathrm{Ne}\right]3s^{2}3p^{3}\right)$ and iron $\left(\left[\mathrm{Ar}\right]4s^{2}3d^{6}\right)$ nuclei and core electrons are modeled by a projector augmented wave potential [39], as generated by Dal Corso [40, 41]. The kinetic energy cutoff for wave-functions is set to 870 eV, which ensures energy convergence below 0.01 eV/atom. Calculations are performed on 54, 76 and 144 atoms supercells for bulk, surface and GB calculations, respectively, with a $4\times 4\times 4$, $5\times 5\times 2$ and $5\times 5\times 2$ gamma centered k-point mesh. First-order Methfessel-Paxton smearing [42] with an energy set to 0.27 eV ensures accurate evaluation of forces for ionic relaxation. Binding energy calculations are relaxed with a conjugate gradient method, until the force on each ion is below 0.025 eV/Å. For GB calculations, supercell relaxation is performed along the direction perpendicular to the GB plane, with a pressure threshold equal to 0.1 kbar. For surface calculations, a void equal to the slab thickness was used. Using a void that is half or twice as thick gives identical surface energy values within less than 1 meV. All calculations are spin-polarized with collinear magnetization to reproduce iron ferromagnetic state. With these settings the lattice parameter for body-centered cubic Fe was found equal to $a=2.826$ Å, in good agreement with experimental measurements $a=2.87$ Å[43] and other DFT calculations $a=2.832$ Å[44]. We computed the elastic constants using three volume-conserving strains of a Fe primitive cell [45, 46] and found the following values $C_{11}=305$ GPa, $C_{12}=146$ GPa and $C_{44}=84$ GPa, in fair agreement with other DFT calculations ($C_{11}=278$ GPa, $C_{12}=148$ GPa and $C_{44}=98$ GPa [44]). Zero-point energies are included for all supercells containing hydrogen atoms, as this correction is not negligible for hydrogen atoms [47, 48]. We could not afford the rigorous calculation of the force constant matrix where all atoms in the simulation cell are displaced so we limit ourselves to displacing the hydrogen atom only, by 0.02 Å in six directions. This approximation is justified by the large mass difference between hydrogen and iron atoms [47, 49]. The zero-point energy correction at $T=0$ K is then expressed as:

[TABLE]

where $\hbar$ is the reduced Planck constant, $m_{H}$ the hydrogen atom mass and $k_{i}$ are the eigenvalues of the force-constant matrix. For any temperature, the H atom vibration contribution to the total energy of the system is [50, 51]:

[TABLE]

and Eq. 2 reduces to Eq. 1 at $T=0$ K. In the calculations performed in Sec. 4 Eq. 2 is used to compute the vibrational contribution to the free energy of the system.

In this paper, a positive binding energy reflects a configuration which is more stable than all solutes being isolated in their reference states. The reference states are chosen as the most stable P and H sites in bulk Fe, substitutional and interstitial tetrahedral, respectively. The binding energy $E^{b}$ between P and H in a bulk system is defined as:

[TABLE]

where $E\left[\alpha\right]$ is the energy of a supercell containing $\alpha$ atoms, computed using DFT. In order to keep the same references to compute binding energies at surfaces or GBs, the following definition is used:

[TABLE]

where $E\left[\alpha,\>d\right]$ is the energy of a supercell containing $\alpha$ atoms and an extended defect $d$ which in our case is either a surface or a GB. With this definition (Eq. 4), the binding energy at an interface contains the P-H interaction at the interface, as well as the P-interface and H-interface interactions.

2.2 Low-temperature expansion

The segregation of chemical species to interfaces is generally quantified using the McLean model [1, 19, 52]. Unfortunately, this model is limited to one segregated species at one specific segregation site. A generalization of this model to multi-component framework was derived by Guttmann [21, 9], based on a mean-field approach. The mean-field does not allow to take into account configurations with specific short-range order and segregation energies. In this paper, the concentrations of segregated elements are computed using a low-temperature expansion formalism (LTE) [53, 54, 55, 56]. This formalism is easier to generalize to several chemical species and a variety of segregation sites and configurations. It is also able to treat systems in both canonical (fixed concentration) and grand-canonical (fixed chemical potential) ensembles. Finally, it is in principle straightforward–even though practically challenging–to take into account the overlap and interaction between segregated species at large coverage. In B, we show that the LTE formalism is rigorously equivalent to McLean’s model in the limit of non-interacting particles having the same segregation energy.

For the system of interest, the grand potential $\mathcal{A}$ is related to the grand-canonical partition function $\left(Z\right)$ with the standard relation:

[TABLE]

where $k_{B}$ is the Boltzmann constant, $T$ is the absolute temperature and the sum runs over micro-states $i$ characterized by energy $E_{i}$ and containing $n_{i}^{\alpha}$ atoms of chemical species $\alpha$. $G_{i}$ corresponds to the degeneracy of micro-state $i$ and $\mu_{\alpha}$ is the chemical potential associated with species $\alpha$. Note that for substitutional solutes, $\mu_{\alpha}$ denotes the difference $\mu_{\alpha}-\mu_{m}$ because it is necessary to remove a matrix atom $m$ from the system in order to insert the solute while keeping the total number of sites constant. We now define a reference state with energy $E_{0}$ and containing $\{n_{0}^{\alpha}\}$ atoms for each species $\alpha$. In our case, the reference system is the ideal interface without any atoms other than matrix (Fe) atoms. The addition of solutes to the system will be treated as excitations with respect to this reference state.

[TABLE]

where $\mathcal{A}_{0}=E_{0}-\sum_{\alpha}n_{0}^{\alpha}\mu_{\alpha}$ and $\delta n_{i}^{\alpha}=n_{i}^{\alpha}-n_{0}^{\alpha}$. The logarithm function is then Taylor expanded to infinite order and the linked-cluster theorem shows that all terms that are non-linear in the number of sites $N$ cancel out [53]. The chemical potential is the energy required to take a given atom from a reservoir to a reference site in the system. Taking these reference sites identical to those in Eqs. 3 and 4, $-\left(E_{i}-E_{0}\right)=E_{i}^{b}$, the binding energy of atoms in configuration $i$. Also, because all remaining excited states are proportional to $N$, we can define $g_{i}=\lim_{N\rightarrow 0}\left[G_{i}/N\right]$ as the degeneracy per site.

[TABLE]

The concentration of atoms of species $\alpha$ is computed as the derivative of $\mathcal{A}/N$ with respect to $\mu_{\alpha}$:

[TABLE]

A similar equation must be written for each species $\alpha$, which gives a system of coupled polynomial equations whose unknowns are quantities $X_{\alpha}=\exp\left(\mu_{\alpha}/k_{B}T\right)$. Either these variables $X_{\alpha}$ are known (grand-canonical ensemble) and it is straightforward to compute the total solute concentrations, or nominal concentrations $\left[\alpha\right]$ are known and one needs to solve the system of equations. The quantity $g_{i}\exp\left(E_{i}^{b}/k_{B}T\right)\prod_{\alpha}X_{\alpha}^{\delta n_{i}^{\alpha}}$ represents the concentration per site of atoms in the specific configuration of micro-state $i$.

The sum over micro-states $i$ in Eqs. 7 and 8 is theoretically infinite but states with a low energy value (binding energy plus chemical potential contribution) are negligible. In other words, we only need to keep in the series states of low excitation energies, which are usually the states with a small number of solute atoms that are different from those of the reference states.

There are two equivalent ways to deal with Eqs. 7 and 8 in our case: either consider a single system containing bulk and interface sites, or two separates systems in equilibrium (i.e. with equal values of $\mathcal{A}$). In this paper, we will consider the first method, which then requires the proportionality factor between the number of bulk and interface sites. Degeneracy factors $g_{i}$ are given with respect to the total number of lattice sites in the system, which are in fact divided into $N_{b}$ bulk sites and $N_{\phi}$ interface sites: $N=N_{b}+N_{\phi}$. Assuming that the interface is a spherical grain boundary, we are able to assign specific values to $N_{b}$ and $N_{\phi}$ as a function of the grain radius $r$:

[TABLE]

where $S_{at}$ and $V_{at}$ are the atomic surface and atomic volume, respectively. In our case, Fe is body-centered cubic such that $V_{at}=a^{3}/2$ ($a$ is the lattice parameter) and the GB plane that we study has orientation $\left(111\right)$ hence $S_{at}=a^{2}\sqrt{3}$. In the end, $\gamma=a\sqrt{3}/2r$. This quantity is useful to define $g_{i}$ coefficients for configurations specific to bulk or interface environment because it is much easier to count segregation sites per interface site for instance. Thus, if we write $g_{i}^{\phi}$ and $g_{i}^{b}$ the multiplicity per interface and bulk site, respectively, the $g_{i}$ coefficients that appear in Eq. 8 are computed as:

[TABLE]

Moreover, it is straightforward to generalize to more complex microstructure with multiple GB types, for instance: $N=N_{b}+N_{\phi_{1}}+N_{\phi_{2}}$ and we compute the $\gamma$ factor for each interface. Also, in Eq. 9, the $\gamma$ factor is defined as the surface over volume ratio of a sphere, but one could include a more sophisticated interface shape.

2.3 Embrittlement potency

Without any solute segregation, the energy of a system containing an interface reads:

[TABLE]

where $E_{m}$ is the energy of one matrix atom, $N$ is the total number of sites in the system, $\gamma_{\phi}$ is the interface energy per unit surface and $\rho_{\phi}$ is the interface area per interface atom. Note that $\rho_{gb}=\rho_{fs}=\rho$ and $N_{fs}=2N_{gb}$ because the fracture of a GB creates two free surfaces. The ideal (reversible) work of separation per GB site ($W_{sep}^{0}$) quantifies the energy change as two grains are separated from each other,

[TABLE]

The embrittlement potency of a segregating element translates its ability to promote or reduce the probability of inter-granular fracture. It is commonly computed as the difference in segregation energy at the grain boundary and at the free surface, as originally proposed by Rice and Wang [25, 11]: $EP=-\left(E_{gb}^{b}-E_{fs}^{b}\right)$ and the minus sign is added because in our convention, a positive binding energy represents an attractive configuration (see Eq. 4). When a single solute segregates to a single type of interface sites without interacting with other solutes, the solute embrittlement potency quantifies the difference in the ideal work of separation due to the presence of the solute. If the fracture occurs fast enough so that solutes do not have time to redistribute between bulk and free surfaces:

[TABLE]

where $\left[X\right]$ is the solute concentration per GB site. Hence, a solute with a positive embrittlement potency (i.e. the solute is more stable at the free surface than at the GB) will promote inter-granular failure as it will reduce the absolute value of $W_{sep}$.

When various solutes interact and are located at various segregation sites, a generalized formula of the embrittlement potency is required. First of all, let us consider various fracture scenarios, as shown in Fig. 1. Figure 1a represents solutes segregated at the GB, at equilibrium where pairs and isolated solutes exist. Figures 1b-1g represent different scenarios of GB separation into two surfaces. In Figs. 1b-1d it is assumed that fracture is much faster than diffusion such that there is no solute equilibration at the surface. In Fig. 1b, all solutes are located on the same surface, while they are evenly distributed between both surfaces in Fig. 1c. If pairs and monomers do not interact (segregated solutes are sufficiently dilute), scenarios 1b and 1c will have the same energy, although scenario 1c cannot result in flat surfaces. Figure 1d shows a scenario similar to Fig. 1c except that solute pairs are possibly separated during the fracture, which might lead to a different free surface energy. In Figs. 1e-1g it is assumed that diffusion is fast compared with the fracture phenomenon, such that atomic species redistribute at free surfaces. In Fig. 1e, all solutes are located on the same free surface, while they are evenly separated between both surfaces in Fig. 1f. Both scenarios might not lead to the same energy because solute distribution depends on solute concentration. Contrary to previous cases, Fig. 1g shows an example where the system is at fixed chemical potential, such that solutes are added or removed in order to equilibrate with the free surface.

To evaluate the embrittlement potency, we will consider two fracture scenarios among those presented here (1b, 1f). Indeed, as our model assumes a dilute solid solution with flat interfaces, scenarios 1b and 1c will give the same energy difference with respect to the GB configuration, and we have no way to quantify the distribution of solutes after fracture in Fig. 1d, as the short-range order remaining at the free surfaces is unknown. Scenario 1e could in principle be computed but it is cumbersome and it is not expected to provide much insight compared with scenario 1f which is much easier to compute. Scenario 1g requires a very slow fracture to enable long-range diffusion during the fracture event, which, in most cases, does not seem realistic. Moreover, because segregation energies to free surfaces are usually much higher than segregation energies to GBs (and this is the case in the Fe-P-H system, cf. Table 2), such scenario would generate very large segregations to free surfaces, much larger than at GBs, thus a clear thermodynamic driving force for GB fracture, which translates into solute embrittlement.

To generalize the definition of the embrittlement potency, we use the results from the LTE calculations to compute the work of separation. We consider the average energy of the interface at thermodynamic equilibrium as being representative of the system. The complete derivation of the following expressions are given in C. We define the generalized embrittlement potency ($GEP$) as the deviation from the ideal work of separation for pure materials, a deviation that is caused by a given concentration of one or several chemical species. For scenario 1b, the concentration is conserved and fracture is fast with respect to the diffusion time scale of solutes, such that the free surface solute segregation is identical to the GB solute segregation:

[TABLE]

For scenario 1f, the fracture is slow compared to diffusion time scales, allowing solutes to reach local equilibrium at the free surfaces, with the constraint of constant local concentration:

[TABLE]

These expressions will be used in Sec. 4 to compute the general embrittlement potency in each case.

3 DFT computations

3.1 Bulk energies

First, we present our DFT results for bulk configurations. The purpose of these calculations is to check the correct convergence of our DFT calculations, and to quantify the interaction between P and H atoms in the bulk, which may affect their distribution and thus their segregation to GBs in fixed concentration samples. All calculations were performed at fixed volume.

The vacancy formation energy in a $3\times 3\times 3$ cubic supercell (54 atoms) of ferromagnetic body-centered cubic Fe was computed as:

[TABLE]

and we found $E^{f}\left(vac\right)=2.19$ eV, in good agreement with previous calculations (2.02-2.20 eV [57]). To ensure that our simulation cell was large enough we used the Aneto code which computes the elastic self-interaction of a point defect from elastic constants and residual stress in the supercell [58]. The energy of elastic interactions with supercell replicas amounts to 0.017 eV for the vacancy, which is of the order of the error bar resulting from our choices of DFT parameters. Hence, we consider that this 54-atom supercell is large enough for bulk calculations.

The energy of a bulk system containing one substitutional P atom was computed for 24- and 54-atom supercells. Both results agree within 0.02 eV and the energy due to elastic interaction with supercell replicas amounts to 0.007 eV. We also computed the energy of an interstitial octahedral P atom but the energy was 3.2 eV above that of substitutional P atom (taking into account elastic corrections), slightly higher than previous calculations (2.8 eV [59] and 3.1 eV [60]). Therefore we will only consider substitutional P atoms hereafter.

The H atom is more stable as an interstitial atom, and can occupy both octahedral and tetrahedral positions. The tetrahedral site was found more stable, 0.15 eV lower in energy than the octahedral site, in agreement with previous calculations (0.13 eV [47] and 0.12 eV [48]). Yet, as pointed out in Sec. 2.1, it is important to account for zero-point energy (ZPE) effects for light elements. This ZPE correction amounts to 0.14 eV and 0.24 eV for the octahedral and tetrahedral sites, respectively, again in agreement with previous results (0.12 eV and 0.24 eV, respectively[48]). After ZPE corrections, the energy difference between tetrahedral and octahedral sites reduces to 0.045 eV only, meaning that both sites must be taken into account. The elastic interaction computed with Aneto amounts to 0.017 eV and 0.002 eV, respectively for the octahedral and tetrahedral sites, again indicating that our 54-atom supercell is large enough. The octahedral site is found unstable (second-order saddle point) and will not be included in our calculations. Therefore, the tetrahedral site is chosen as the reference H configuration to compute binding energies (see Eqs. 3 and 4).

Figure 2 and Table 1 show the various P-H configurations and the associated binding energy. At small distance, P and H have a negative binding energy, resulting from a repulsive configuration. We also tested the tetrahedral site located in between octahedral sites 1 and 2, but this configuration spontaneously relaxes towards configuration 3. As the P-H distance increases, binding energies converge towards 0, which is the expected behavior with no long-range interactions. From this plot, we can estimate the range of P–H interactions to be about one Fe lattice parameter, with a strong repulsion at small P-H distances. The calculation of vibrational modes for H shows that configurations $\mathrm{PH_{1}}$ and $\mathrm{PH_{4}}$ correspond to second-order saddle-points, and configurations $\mathrm{PH_{2}}$ and $\mathrm{PH_{6}}$ correspond to first-order saddle-points. These configurations are therefore unstable and will not be included in the calculations.

3.2 Grain boundary and surface energies

This section presents the P and H segregation energies to $\Sigma 3$ $109.5\degree\left[0\bar{1}1\right]\left\{111\right\}$ GB and $\left\{111\right\}$ free surface, as well as P-H binding energies at these interfaces. The simulation cells are shown in Fig. 3. They contain 4 atoms per interface plane which is defined by two vectors of length 0.801 nm forming a 60°angle corresponding to $\left\langle 110\right\rangle$-type directions (see Fig. 5).

The simulation cells were relaxed only in the direction perpendicular to the interface plane. If we had a very large box in the direction normal to the interface plane, and we would relax the simulation cell in all three directions, the relaxation would only occur in the direction normal to the interface plane, because a relaxation in a direction parallel to the interface plane would modify the inter-atomic distance for all bulk atoms and thus increase their energy, and as they would be much more numerous than interface atoms, such relaxation would not be energetically favorable. In our calculations we are limited in the size of the simulation cell we study. Even though we have checked that the simulation cell was large enough to avoid elastic interactions between both interfaces (10 mJ/m2 variation in going from 12 to 24 atoms in between interfaces), the number of bulk atoms is not so large compared with the number of interface or near-interface atoms, such that some relaxation of the cell parallel to the interface plane might be energetically favorable. Hence, it seems more realistic to freeze the degrees of freedom parallel to the interface plane when relaxing the simulation cell. This may explain the difference in GB energy found with previous DFT calculations were the simulation cell was relaxed in all three directions (all values are given in J/m2): 1.46 [61] and 1.52 [35, 62] whereas we find 1.23. This value is consistent with molecular static simulations using semi-empirical interatomic potentials, which are performed on much larger cells than DFT calculations: 1.25 [63], 1.31 [27] and 1.55, 1.39, 1.30, 1.23, 1.22 using various potentials [35]. We find the $\{111\}$ free surface energy at 2.07 (previous DFT calculations give 2.69 [62]).

For the GB, we performed spin-polarized calculations with both collinear and non-collinear magnetic moments. Over all atoms, the maximum relative difference in magnetic moment was 0.6% which is negligible. Hence we performed all subsequent calculations with collinear magnetic moments, which is cheaper in terms of computational cost.

To reduce the number of segregation sites to study we analyzed the symmetry of our simulation cell, and found an order 3 rotation symmetry around the $[111]$ axis as well as a symmetry plane normal to the $[110]$ direction. These symmetries are clearly seen in Fig. 4 where black-diamond shape lines show the unit supercell for this GB, and our actual supercell is composed of 4 unit supercells, as shown in Fig. 5. Figure 4 additionally shows, for each point in the simulation cell, the distance to the closest Fe atom. This procedure aims at identifying potential interstitial segregation sites for H atoms. Red regions are the ones where there is the most space, so places where H atoms are expected to segregate. Hence we identified a collection of sites that are all equivalent when H is the only segregating species and infinitely dilute, represented by site 0 in Fig. 5. We tried to insert P on this interstitial position but the calculation would not converge and the site was unstable. We conclude that P segregation is substitutional at this GB, in agreement with previous work [64] but this might be true only for low temperatures [65]. For P–H pairs, when P is placed at a substitutional site at the interface, interstitial sites in the GB plane become nonequivalent (numbered from 1 to 6, a H atom located at site 0 relaxes toward site 3 when P is added). Other sites were tested (red/yellow regions which denote available space in Fig. 4) but relaxed towards one of these sites.

Table 2 shows the results of these DFT calculations in both GB and free surface for isolated P and H, P and H together but considered infinitely separated from each other, and P-H pairs with P being at a substitutional interface site and H being at an interstitial site numbered as in Fig. 5. The last three columns show the difference between the P–H pair being at short distance and P and H being infinitely separated from each other, and thus quantifies the magnitude of the P–H interaction at the interface. First of all, the binding energy of H to the GB (including ZPE correction) is much higher than the binding energy of P to the GB. Hence, at similar bulk concentration levels, P is expected to be much less segregated to the GB than H. A similar comment holds for a free surface, even though the difference between the two species is much smaller. However, the embrittlement potency ($EP=E_{fs}^{b}-E_{gb}^{b}$) is positive for both species and much smaller for H than for P, which means that these solutes thermodynamically favor free surfaces and thus GB fracture, but P has a much stronger effect than H. We find an embrittlement potency of 0.25 eV/at for H–in good agreement with previous work: 0.26 eV/at [66, 16], 0.33 eV/at [67] and 0.41 eV/at [61]–and 1.48 eV/at for P while previous DFT results show widely spread values: 1.63 eV/at [68], 0.79 eV/at [6], 0.17 eV/at [69].

Looking at co-segregation effects in detail, all P–H configurations are repulsive at the GB and only one configuration is attractive at the free surface. Thus, P–H interactions will tend to lower slightly the segregation levels at the GB and increase slightly the segregation levels at free surfaces, meaning that they provide an additional thermodynamic driving force for GB fracture. Looking at the last column of Table 2, the embrittlement potency difference is positive for 3 out of the 5 configurations tested, confirming the additional embrittlement of the GB caused by the short-range order between P and H. Note that configuration $\mathrm{PH}_{4}$ is in fact found unstable (first order saddle point) when computing vibrational modes for the H atom, and this configuration will not be included in our calculations.

In short, at equivalent bulk concentrations, P segregation to the GB is lower than that of H, but at equivalent segregation levels, P makes the GB more brittle than H. When P and H are included simultaneously in the system, segregation levels are expected to be mostly as if there was no interaction between P and H but some configurations cause more GB embrittlement than isolated P and H. The next step consists in computing how much of these P–H pairs form as a function of temperature and nominal solute concentrations and thus how the embrittlement potency of the system is affected. To this end, we will use the LTE model presented in Sec. 2.

4 Effect of co-segregation on embrittlement

4.1 Effect of grain size on segregation

One of the nice features of the LTE model is that it is straightforward to take into account the effect of grain size, which translates into a contribution to configurational entropy: the larger the grain, the lower the number of interface sites at constant total number of sites, hence the lower the configurational entropy associated with interface sites. In this study, the ratio between the number of interface sites and the total number of sites is computed in the ideal case of spherical grains (cf. Eq. 9).

Figure 6 shows the effect of grain radius for isolated species (no P–H pair formation is allowed). As shown in Table 2, the H binding energy to the GB is about 4 times larger than the P binding energy to the GB. Hence H segregation is expected to be much more important than P segregation. This is indeed found from the two first sub-figures, where solute concentrations in the GB plane and in the bulk were computed at T=300 K for a 10 appm nominal solute concentration. At identical temperature and nominal solute concentrations, GB segregation is higher for H than for P, and for both species it decreases with increasing grain radius. This behavior is also expected because higher grain size favor bulk sites over GB sites in terms of configurational entropy.

A less-intuitive conclusion drawn from these calculations is that GB coverage increases and then seems to saturate with increasing grain size. GB coverage is the number of solutes at GBs per GB sites, while GB concentration is the number of solutes at GBs per total number of sites in the system. As grain radius increases, the average solute energy in the system–its chemical potential–increases because more and more solutes are added to the bulk where they are less stable than at GBs. At some point, most solutes are in the bulk such that increasing the grain radius further does not really affect the average solute energy in the system, hence the saturation of GB coverage. Figure 6 shows that GB coverage and bulk concentration have the same variations, which comes from the fact that they are both proportional to the exponential of the solute chemical potential. This variation of GB coverage with grain size is clearly seen from LTE expressions, considering only one species (H) without geometrical frustration effects, to keep the discussion simple. Also we assume that $1/(1+\gamma)\simeq 1$, which is true except for very small grains (a few nm). Defining $S_{H}=3\exp\left(F_{gb}^{b}\left(H^{gb}\right)/k_{B}T\right)$ and $Y=\exp\left(\mu_{\mathrm{H}}/k_{B}T\right)$, the total H concentration is the sum of bulk H atoms and H atoms segregated at the GB:

[TABLE]

The GB coverage by H atoms is expressed as:

[TABLE]

Hence at given H nominal composition and temperature, higher grain radius $r$ results in high GB coverage until a saturation value is reached, $\theta_{H}\simeq[H]S_{H}/6$.

The second sub-figure in Fig. 6 (H segregation at T=300K) is particularly interesting because GB coverage saturates at 1, which means that all interstitial sites in the GB plane are occupied by H atoms. When concentrations reach high levels on a specific type of sites, one must be cautious with LTEs, because geometrical frustration effects must be added. They consist in adding excitation states where solutes are non-interacting. These states have a negligible probability at low solute concentrations but become necessary when concentration rise above $\simeq 1\%$. Generally speaking it is a difficult task to compute these frustration effects. Yet, we can perform the exact calculation when solutes are considered non-interacting, as it is done in B and C. In this Section, all calculations were performed taking into account geometrical frustration in this ideal case for solutes at interfaces, treated in a mean-field way. This way of dealing with interface concentration issues is not exact but allows to keep calculations not too complicated. Note that only geometrical frustrations are treated in a mean-field way, segregation being treated exactly and independently for each type of segregation site. D presents the derivation of the mean-field geometrical frustration and an example of the equations solved in this study.

4.2 Effect of temperature on co-segregation

Figure 7 shows how GB coverage evolves as a function of temperature in two cases: either P and H solutes are treated separately, i.e. they do not interact (solid lines); or P–H interactions–summarized in Tables 1 and 2–are accounted for (symbols). Hence, the difference between symbols and solid lines represents the effect of co-segregation.

The decrease of GB coverage with temperature is trivial as segregation is proportional to the exponential of the segregation energy over $k_{B}T$. It is interesting to note the asymmetry in H and P behavior: P–H interactions increase P segregation at GBs by a factor lower than 2, and this effect is temperature dependent, with a maximum at $T\simeq 650K$ (see inset). On the other hand, H segregation is not affected by 100 appm of P, whatever the temperature, and additional calculations not shown here confirm that H segregation is unaffected by the presence of P in the system. This can be understood from Fig. 6: at a grain radius of $10\mu$m, H concentration at GBs is 1-2 orders of magnitude higher than that of P. Hence, most P atoms can pair with H atoms while leaving most H atoms isolated, i.e. without changing the H chemical potential. The reverse asymetric behavior is expected at free surfaces because P segregation energy is much higher than H segregation energy.

Binding energies being all repulsive for P–H pairs at GBs (Table 2), the increase of P segregation in the presence of H is due to P–H pairs that are stabilized by configurational entropy only. The inset shows that the co-segregation effect–quantified by the ratio of P segregated concentration with and without H in the system–presents a maximum at T$\simeq$650 K. As temperature increases, the configurational entropy contribution to the free energy of the system becomes more and more significant, and H concentration at GBs is more or less constant, close to a full coverage of the GB plane. At some temperature, the H coverage at GBs decreases, which reduces the probability of forming P–H pairs, and therefore the stabilization of P atoms by configurational entropy also decreases.

To summarize, co-segregation effects on segregation are mainly independent of P concentration and triggered by a large H concentration at GBs, which increases P segregation. It is therefore interesting to take a look at H segregation at GBs which depends on H nominal concentration, temperature and grain radius. Figure 8 shows the H GB coverage (H concentration per GB site at the interface plane) as a function of grain radius and H nominal concentration for two temperatures.

At low temperature, constant GB coverage contours appear as straight lines in Fig. 8, either vertical or with a slope equal to -1. This may be explained from the LTE expressions. Taking the natural logarithm of Eq. 18 we get:

[TABLE]

from which we identify two limiting behaviors: if $\sqrt{3}S_{H}\ll 12r$ (large grain radius and/or high temperature), Eq. 19 simplifies to $\ln\left[H\right]\simeq\ln\left(6\theta_{H}/S_{H}\right)$, which does not depend on the grain radius and should indeed give vertical constant coverage contour in Fig. 8; on the contrary, if $\sqrt{3}S_{H}\gg 12r$ (low grain radius and/or low temperature) we get a linear relation with a slope equal to -1 between the logarithms of grain radius and nominal H concentration $\ln\left[H\right]\simeq\ln\left(\sqrt{3}\theta_{H}/2\right)-\ln\left(r\right)$. Hence, Eq. 19 qualitatively explains the trends observed in Fig. 8 and shows that the comparison between the H segregation energy over $k_{B}T$ and the grain radius is the key quantity to predict H coverage and thus the possible P–H co-segregation effects.

4.3 Combined effects of grain radius, temperature and H nominal concentration on co-segregation

As in the inset of Fig. 7, the magnitude of co-segregation effects is expressed as the ratio between the computed P concentration at the GB and the P concentration obtained when H is removed from the calculation (hence no P–H interactions, denoted by subscript 0). The additional segregation of P to GB due to H is a function of temperature, grain radius and H nominal concentration. In this section we want to understand the interplay between these parameters. To this end, Fig. 9 presents the $[P]^{gb}/[P]^{gb}_{0}$ ratio as a function of temperature for various combinations of grain radius and H concentration. The lower plot shows the corresponding GB coverage of H, which is a function of the three same parameters. The curves that we obtain have the same shape as the one in Fig. 7, with a maximum located more or less at the point where H concentration at the GB becomes lower than unity.

Focusing first on the lower plot of Fig. 9, the curves show three values of H coverage at room temperature. This behavior is explained in Eq. 19 by the comparison between $S_{H}$ and $r$. By comparing some of the curves, we see that increasing the grain radius by one order of magnitude while decreasing the nominal concentration by the same amount does not change the low temperature coverage, but modifies the temperature evolution of $\theta_{H}$.

Looking now at the top plot of Fig. 9, the curves also seem to form three groups of curves, the ones with a maximum $[P]^{gb}/[P]^{gb}_{0}$ ratio between 1.4 and 1.7, the ones with a maximum ration between 1.1 and 1.2, and finally the one with almost constant ratio, and these "groups" correspond to the ones identified for H coverage. Generally speaking, increasing the grain radius increases the co-segregation effects (higher maximum of the curves) and translates the maximum ratio towards lower temperatures; increasing the H nominal concentration increases the co-segregation effect and translates the maximum ratio towards higher temperatures. Hence, these two variables do not have the same effect, because they do not affect the temperature evolution of the H coverage in the same way. The comparison between the top and bottom plots of Fig. 9 clearly demonstrates that the magnitude of co-segregation effects is directly linked with the fraction of interface sites occupied by H atoms.

4.4 Generalized Embrittlement Potency: rapid fracture at constant concentration

We now discuss the general embrittlement potency (GEP) obtained for scenario 1b where segregated solutes are assumed to be immobile during the fracture phenomena, such that they are simply transferred from the GB to the free surface, and the GEP is expressed by Eq. 14. Figure 10 shows the ratio between the GEP and the one that would be obtained without P–H interactions, so as to characterize co-segregation effects on embrittlement potency. This ratio is shown as a function of solute nominal concentrations for various temperatures and grain sizes. From previous paragraphs, we know that P–H interactions lead to larger concentrations of segregated P for H coverage close to unity, hence low temperature and/or large grain radius, such that the qualitative trends shown in Fig. 10 are expected. In addition to larger quantities of segregated P, ab initio calculations showed that P–H pairs mostly lead to increased embrittlement potency (cf. Table 2), adding up with increased segregation levels and leading to larger GEP, i.e. more brittle interfaces. Nevertheless, note that the GEP potency increase due to co-segregation remains quite low, always lower than 10% and occurs only at fairly high nominal concentrations of P (at least 100 appm). This is at variance from the effect on segregation levels which was mainly independent from the concentration of segregated phosphorus. In the P–H system, the embrittlement potency is essentially an average of the embrittlement potency of each species, weighted by their concentration. Hence, if H is much more concentrated than P at the interface, the contribution of P atoms on the embrittlement potency will be essentially negligible. If H concentration at the interface is low, then there is no co-segregation effect.

4.5 Generalized Embrittlement Potency: slow fracture at constant concentration

The general embrittlement potency obtained for scenario 1f characterizes what happens during a slow fracture, when segregated solutes may diffuse locally and equilibrate at the free surface during the fracture. The corresponding GEP is expressed by Eq. 15. As shown in Fig. 11, the effect of co-segregation on embrittlement is lower in this scenario than in the one describe in the previous paragraph, as in the range of parameters investigated, the GEP do not exceed a 5% increase due to co-segregation, even though it occurs at lower nominal concentrations. The plots are also more symmetric than the ones obtained for fracture scenario 1b, because, as H would increase P segregation at GBs, P increases H segregation at free surfaces, and this feature is now accounted for since solutes are equilibrated at free surfaces. Hence for each species there is some threshold concentration value above which co-segregation translates into increased GEP, either because the GB plane is covered with H atoms, or because the free surface plane is covered with P atoms. When both species are simultaneously in concentrations close to these threshold values, the increased GEP is also observed.

5 Discussion about current model limitations

The aim of this paper was mainly to introduce the new LTE statistical model for solute segregation at interfaces. The application to P–H co-segregation effects in $\alpha$-Fe was limited by the amount of computational time available for DFT calculations.

The first way to improve the model is to perform additional DFT calculations, either on other sites around the GB because some out-of-plane sites may be found very stable [64, 26]. Also, a similar type of work should be performed on other GBs because a real material contains various types of GBs [34]. Moreover, highly coherent GB such as first-order twins considered in this study may exhibit huge property changes if the GB plane is misaligned by only a few degrees [4]. These types of complex structures–for instance non-symmetric GBs–are now becoming accessible, especially using molecular dynamics approaches based on semi-empirical interatomic potentials [63, 70, 71]. Note that having a statistical model such as the one proposed in this paper that is able to treat each segregation site individually but also collections of segregation sites (i.e. short-range order) is particularly useful for grain boundaries in ordered structures, e.g. oxides [18].

Another point which is crucial–especially for comparison with experiments–is to take into account temperature effects. For instance, DFT calculations at 0 K indicate substitutional segregation of phosphorus while there is experimental evidence for interstitial segregation at higher temperature [65]. The discrepancy may be due to the temperature dependence of the segregation energy (entropy contribution) and magnetism evolution as the Curie temperature is approached. More generally speaking about finite temperature effects, segregation entropy is important to consider for relevant modeling/experiment comparison[52, 18]. In our LTE model, we can easily use a free energy for each configuration. As a matter of fact, we already included configurational entropy and a first-order approximation of vibrational entropy for H atoms in the present work. Other forms of entropy (e.g. electronic, magnetic and vibrational) can sometimes be computed with state-of-the-art DFT tools, but these remain computationally expensive, especially for complex GB structures. Temperature also affects the lattice parameter, magnetic ordering and GB structure which are all non trivial contributions. Because DFT calculations cannot consider all these contributions altogether, a closely bound experimental/numerical work is probably the best way to obtain meaningful data on solute segregation at GBs.

One limitation of LTE models is their inability to treat phase transitions in the case of solute precipitation at the GB. There is a workaround though, but two separate LTE calculations must be performed–one with a reference free interface and the other with a reference precipitate phase at the GB–and then both systems must be put in equilibrium to solve concentrations and/or chemical potentials[72]. The main difficulty in performing LTE calculations is the computation of the so-called "counter-terms", i.e. the geometrical multiplicity of states that are dissociated[55]. These states are usually negligible for dilute concentrations (typically $<1at\%$ for a given type of site). For more concentrated cases these terms should be added, but it is always possible to check whether or not the next term in the expansion is negligible. Hence, one can safely use LTE expressions in more and more concentrated states without knowing the full (infinite) partition function of the system.

Finally, we were only focused on equilibrium segregation, while real materials contain multiple segregating elements, each having their own kinetic properties [73]. Thus out-of-equilibrium segregations may arise, for instance during quenching [74, 75], under irradiation [76, 77] or because of GB migration [78]. Note that in all these cases the elastic field around the GB modifies the diffusion properties of solutes [79] and thus the segregation.

This rather long list of model limitations does not mean that DFT calculations are useless to study GB segregation, but rather points towards some critical points to keep in mind when using such models and comparing them with experiments.

6 Conclusion

The main contribution of this paper is the development of the low-temperature expansion model to study interface segregation. The strength of this model is its general formulation which makes it suitable to address various segregation problems, for instance taking into account multi-species and multi-sites per species segregation, as well the effect of grain size on segregation.

To demonstrate the applicability of the LTE formulation, we applied it to the study of P–H co-segregation effects in $\alpha$-Fe twin boundaries. We performed ab initio calculations to obtain the input data necessary for the LTE model, mainly binding and segregation energies.

We showed the importance of grain size on GB coverage, especially for small grains. The classic McLean segregation isotherm [1] is basically restricted to the large grain size limit, where the solute chemical potential is given by the bulk solution energy. P and H solutes show asymmetric co-segregation behavior at the GB plane. H segregation levels are insensitive to the presence of P, while H atoms can slightly increase P segregation to GBs due to configurational entropy. For this to happen, the GB plane must be almost saturated with H atoms, to maximize the probability of forming P–H pairs. Hence, the H coverage at GBs is the key quantity, and depends on the grain radius, temperature and H nominal concentration. The situation is reversed at free surfaces where H may be stabilized by P atoms.

A general formula for embrittlement potency is proposed, and we show that co-segregation can lead to increased embrittlement potency, up to 10% depending on solute concentrations, temperature, grain size and the fracture scenario considered (comparison between fracture and local diffusion time scale).

Larger co-segregation effects would be obtained if P and H showed attractive interaction in the GB plane. This might be the case in other GBs, or for other solutes at a twin boundary.

Because H easily saturate GBs, future work on segregation at $\alpha$-Fe twin boundaries should focus on concentration effects, introducing other clusters (H–H interactions for instance) and out-of-plane configurations. For other types of GBs, the strain field around the GB may favor these out-of-plane configurations. Finally, let us stress that the ideal work of interfacial separation is predicted to have an important but probably not exclusive role in controlling embrittlement [25]: for instance, McMahon emphasizes the dynamic effect of H as an explanation to GB embrittlement [33]. Hence, equilibrium segregation studies provide necessary but not sufficient insight into embrittlement phenomena.

Acknowledgments

The authors gratefully acknowledge the various computing centers where calculations were performed for this study: the HPCaVe center at UPMC-Université, the CENTAURE cluster at École Nationale Supérieure des Mines de Saint-Étienne, the P2CHPD center at Université de Lyon and the CALMIP center at Université de Toulouse (UMS CNRS 3667). Part of this study was funded by the project GIBBS-ANR-15-CE30-0016, and the authors acknowledge financial support from the NEEDS-SEDI project and D. Connétable for his contribution in getting access to computing resources.

Appendix A Discussion about the reliability of DFT calculations for solutes with low solubility limits

Ab initio methods have their drawbacks, mostly related to the complexity of properly taking into account temperature effects, and the limited number of atoms that can be included in a simulation box. Regarding the latter, let us comment on some question that was raised about the reliability of ab initio calculations for solutes with small solubility limits [17, 18]. We do not question the discrepancy between experimentally measured and computed segregation energy that was identified by the authors, but the "low solubility limit" argument is not convincing, mainly because it is not an intrinsic parameter of a solute in a matrix: the solubility limit depends on the phases in the material, not only on the thermodynamically stable ones but on the ones that are actually formed. So even though a solute with a solubility around 1 appm is severely supersaturated in a 100-atom simulation cell, it cannot precipitate because at best it will form a lattice with more than 1 nm between each atom, which is usually not the stable phase required to impose this solubility limit. The convergence of the energy of the system with the cell size demonstrates that there is no interaction between a solute and its replicas such that in the simulation cell, the solute cannot be out-of-equilibrium because solubility limit is not really defined. In our opinion, the discrepancy is rather due to complicated temperature effects and experimental difficulties in measuring segregation energies, all of which being emphasized for the solutes with "small solubility limits". But the ab initio calculation itself is correct at $T=0$ K.

Appendix B Equivalence between LTE and McLean formalisms

In this Appendix, we demonstrate the equivalence between the low-temperature expansion formalism that is used in this paper, and the standard McLean model which is commonly used to compute the amount of segregated species [20, 1, 19, 52]. In the Langmuir-McLean model, there is only one type of solute, segregating at one type of site (single segregation energy) and segregated atoms do not interact with one another. Hence, the solute-GB binding energy of state $i$ containing $n_{i}^{\alpha}=i$ solutes $\alpha$ is $E_{i}^{b}=iE^{b}_{\alpha}$, where $E^{b}_{\alpha}$ is the segregation energy of species $\alpha$ at the GB.

From Eq. 6, we need to compute the degeneracy for each state to be able to compute the partition function of the system. Fortunately, for the specific case we are interested in, it is possible to give a general formula for each $G_{i}$ factor and write the infinite series of the partition function. The change in $G_{i}$ values as a function of the number of solutes in the system is mainly related to exclusion effects, because two solutes cannot occupy the same site. These exclusion effects are taken into account via "counter-terms" in the LTE formalism, which corresponds to $G_{i}$ factors for non-interacting of dissociated configurations. Also note that all solutes are not distinguishable hence the division by the number of permutations for $i$ solutes.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The total energy of each state is given by $i\left(E_{i}^{b}+\mu_{\alpha}\right)$. Hence, defining $X=\exp\left(\left(E^{b}_{\alpha}+\mu_{\alpha}-\mu_{m}\right)/k_{B}T\right)$, Eq. 6 becomes:

[TABLE]

and the second equality makes use of the linked-cluster theorem [53].

Matrix $m$ and solute $\alpha$ concentrations at the interface $\phi$ are obtained from the derivatives of Eq. 24 with respect to the corresponding chemical potentials:

[TABLE]

such that $\left[m\right]^{\phi}+\left[\alpha\right]^{\phi}=1$ and :

[TABLE]

The exact same work can be performed in the bulk rather than at the interface, the only difference being that there is no segregation energy in the bulk: $E^{b}_{\alpha}=0$. Thus, we recover the Langmuir-McLean segregation isotherm:

[TABLE]

Appendix C Generalized embrittlement potency

C.1 Average energy of a system

In this Appendix we derive the expressions presented in Sec. 2.3. To generalize the definition of the embrittlement potency, we use the results from the LTE calculations to compute the work of separation. We consider the average energy of the interface at thermodynamic equilibrium as being representative of the system. Let $\beta=1/k_{B}T$ and $\mathcal{A}_{i}=E_{i}-\sum_{\alpha}n_{i}^{\alpha}\mu_{\alpha}$. Then the partition function of the system reads $Z=\sum_{i}G_{i}\exp(-\beta\mathcal{A}_{i})$, which is nothing but the sum of the probability of each micro-state. The average energy of the system (including chemical potential contributions in the grand-canonical ensemble) is:

[TABLE]

Let us define the binding grand potential of a micro-state, defined with respect to reference state $i=0$: $\mathcal{A}_{i}^{b}=-\left(\mathcal{A}_{i}-\mathcal{A}_{0}\right)=E_{i}^{b}+\sum_{\alpha}\delta n_{i}^{\alpha}\mu_{\alpha}$. Using the expression from Eq. 7 and these notations:

[TABLE]

leading to:

[TABLE]

which is conveniently rewritten as:

[TABLE]

where $\left[\alpha\right]$ is the total concentration of species $\alpha$ in the system, given by Eq. 8. The transition from Eq. 29 to Eq. 30 is exact as long as all entropy contributions in $\mathcal{A}_{i}^{b}$ are temperature-independent, as the configurational entropy for instance. In our calculations, we consider temperature-dependent vibrational entropy for H atoms (cf. Eq. 2) but for simplicity we still use Eq. 30 to compute the general embrittlement potency.

C.2 Comparison with the ideal embrittlement potency in the dilute case

With this expression of the average energy of the system, let us recover the basic expression of the embrittlement potency in the case where there is only one segregating species, one segregation site, and the interface concentration of solute is very low.

[TABLE]

If all the atoms segregated at the GB stay at free surfaces without short-range order re-organization (fast fracture, as in scenario 1b or 1c), then the probability of each configuration–roughly $\exp\left(\beta\left(E_{i}^{b}+\mu_{\alpha}\right)\right)$ is identical at the GB and at the free surface. Also the chemical potential is not really defined once the system is fractured because the system is not at equilibrium. Yet, if we assume that there are more solute atoms in the bulk than in the GBs or free surfaces-–because grains are large enough so that the number of bulk sites is much larger than the number of interface sites—then we can assume that the solute chemical potential is mainly given by bulk solute atoms such that $\mu_{\alpha}^{gb}\simeq\mu_{\alpha}^{fs}$. For scenario 1b, we will always make this assumption since the chemical potential is ill-defined after fracture. The ideal work of separation per GB site is:

[TABLE]

with $X_{gb}=\exp\left(\beta\left(E_{gb}^{b}+\mu_{\alpha}\right)\right)$ is the concentration of solute atoms at the interface, sometimes called the interface coverage (this expression is valid for low coverage only) and the embrittlement potency for solute $\alpha$ is defined as in Sec. 2.3: $EP_{\alpha}=E_{fs}^{b}-E_{gb}^{b}$. If the solute is more stable at the GB than at the free surface ($E_{gb}^{b}>E_{fs}^{b}$), the embrittlement potency will be negative and solute segregation increases $W_{sep}$ meaning it strengthens the GB.

C.3 Comparison with the ideal embrittlement potency including geometrical frustration effects

Now let us do a comprehensive calculation for this specific case of one type of solute segregating at one type of interface site, taking full account of concentration effect. The expressions for $G_{i}$ coefficients are the ones shown in B (Eq. 23), with $N=N_{gb}$ and we want to compute the ideal work of separation per unit GB site, using $g_{i}=\lim_{N_{gb}\rightarrow 0}\left[G_{i}/N_{gb}\right]$. The energy of state $i$ where $n_{i}^{\alpha}=i$ solute atoms have been added is $E_{i}^{b}=iE^{b}$ because solutes are non-interacting.

[TABLE]

Knowing that

[TABLE]

we find the following expression, which is similar to Eq. 32 for low $X_{gb}$ values:

[TABLE]

C.4 Generalized embrittlement potency expressions

We now derive the expressions presented in Sec. 2.3 to define the generalized embrittlement potency for scenarios b, f and g in Fig. 1. We define the generalized embrittlement potency ($GEP$) as the deviation from the ideal work of separation for pure materials, a deviation that is caused by a given concentration of one or several chemical species:

[TABLE]

The GEP is expressed as a per GB site quantity and the minus sign is to conserve the convention that a positive $GEP$ translates into GB embrittlement.

For scenario 1b, concentration is conserved because the fracture event is much faster than the solute diffusion time scale. Also, because the final state is not at equilibrium, we assume that the chemical potential is mainly fixed by bulk atoms and thus unchanged during fracture. The contribution of bulk atoms and solutes to the overall energy of the system is also identical before and after fracture because solutes do not have time to diffuse. The probability of finding a segregated solute at a free surface is identical to the probability of finding a segregated solute at the GB, only the segregation energy changes.

[TABLE]

where the sum runs over excited states $i$ involving GB sites only.

Scenario 1f is a bit more complicated because upon fracture the system initially containing $N_{b}$ bulk sites, $N_{gb}$ GB sites and $\{N_{\alpha}\}$ solutes is split into two independent sub-systems, each containing $(N_{b}-N_{gb})/2$ bulk sites, $N_{gb}$ free surface site and $\{N_{\alpha}/2\}$ solutes. The concentration in each sub-system after fracture is:

[TABLE]

Still there is a difference when we solve the LTE equations because the $N_{\phi}/N_{b}$ ratio is not the same in a sub-system containing a free surface and in the initial system containing the GB:

[TABLE]

with $\gamma^{gb}$ being defined in Eq. 9. When using Eq. 10 to switch from per interface site multiplicities to per total number of sites multiplicities, the value $\gamma^{fs}$ defined in Eq. 39 will be used to compute the grand potential of the system containing a free surface. The GEP for scenario 1f requires two separate LTE calculations (where the same total species concentrations are imposed) and is written as:

[TABLE]

Appendix D Mean-field derivation of geometrical frustration at the interface

In this Appendix we derive a mean-field version of the geometrical frustration at the interface. This contribution must be taken into account in LTE calculations when one type site or configuration has a concentration per site above a few percent. Due to the effect of grain size on segregation (cf. Sec. 4.1) and the high segregation energies of solutes to interfaces (even more so for free surfaces), interface coverage close to unity often arises. With LTEs, one can take into account chemical short-range order exactly; while it is straightforward to do so for a set of well-identified configurations, it is much more complicated to do the combinatorics for "dissociated" configurations (meaning all configurations that are not well-identified). Therefore, in this study, we simplify the calculation of the geometrical multiplicity of these dissociated configurations using a mean-field approach. Short-range order for well-defined configurations is still taken into account exactly. Let us introduce some handy notations to write the derivation in a more concise manner.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $F^{b}$ and $F^{b}_{\phi}$ are the free energies (including vibrational entropy contributions for H) for bulk and interface configurations, with values taken from Tables 1 and 2. With these handy notations, let us introduce all these configurations–without introducing any geometrical frustration effect–in the expression of the grand potential of the system (Eq. 6), the $\gamma$ factors being introduced from Eq. 10 to take into account the effect of grain radius.

[TABLE]

for which we define the total concentration of segregated solutes at the interface:

[TABLE]

This is where the mean-field approximation comes in: we assume that each segregated solute–whatever its chemical species and/or short-range order–occupies a number of sites at the interface which are therefore no longer available for other species. The exact calculation would require to distinguish between solute type and local configuration. The advantage of this mean-field approximation is that we can use the expressions that have already been used in the previous appendices to take into account geometrical frustration for a single type of solute, the concentration of the unique type of solute becoming the mean quantity of solutes segregated at the interface, $Z_{\phi}$. Hence, Eq. 47 becomes:

[TABLE]

As previously, we use the linked-cluster theorem to get rid of all terms that are non-linear in the number of sites in the system $N$,

[TABLE]

and the total P concentration is obtained as follows:

[TABLE]

The geometrical frustration appears in the $1/(1+Z_{\phi})$ division, and we verify that it is equal to 1 at low interface coverage (i.e. $Z_{\phi}\ll 1$). Note that bulk configurations are always dilute in the system under study, such that a similar treatment is not required. Equation 51 can be recast into a second-order polynomial of variable $X$ (first unknown) where each coefficient is a polynomial function of variable $Y$ (second unknown).

[TABLE]

A similar equation is obtained for the total H concentration:

[TABLE]

and Eqs. 52 and 53 introduce $\tilde{\gamma}_{\alpha}$ for $\alpha=$H or P:

[TABLE]

Equations 52 and 53 correspond to the system of coupled polynomial equations that we solve to obtain P and H chemical potentials at fixed nominal concentrations, and determine whether the presence of one solute alters the segregation behavior of the other.