# On the effect of elastic anisotropy and polarizability on solute segregation at low-angle grain boundaries

Joé Petrazoller, Julien Guénolé, Stéphane Berbenni, Thiebaud Richeton

► To cite this version:

Joé Petrazoller, Julien Guénolé, Stéphane Berbenni, Thiebaud Richeton. On the effect of elastic anisotropy and polarizability on solute segregation at low-angle grain boundaries. 2024. hal-04793337

HAL Id: hal-04793337

https://hal.univ-lorraine.fr/hal-04793337v1

Preprint submitted on 20 Nov 2024

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
# On the effect of elastic anisotropy and polarizability on solute segregation at low-angle grain boundaries

Joé Petrazollerª, Julien Guénoléª, Stéphane Berbenniª, Thiebaud Richetonª,*,

ªUniversité de Lorraine, CNRS, Arts et Métiers, LEM3, F-57000 Metz, France

## Abstract

Solute segregation towards grain boundaries is investigated by modeling solute atoms as elastic dipoles interacting with the strain fields of symmetric tilt low-angle grain boundaries (LAGBs). Elastic dipoles are determined using molecular dynamics (MS) considering both the permanent second-rank tensor and the fourth-rank polarizability tensor, which is needed to capture the elastic dipole dependence on external strain. For cubic lattices, the latter tensors are related to size and modulus effects, respectively. The strain fields of LAGBs are evaluated either through MS or by considering arrays of edge dislocations within the frameworks of linear isotropic elasticity or heterogeneous anisotropic elasticity using the Stroh formalism. The interaction energies arising from the coupling between elastic dipoles and LAGB strain fields are compared to segregation energies computed on a site-by-site basis using MS. These comparisons are made for three LAGBs and two cubic systems (Cu and Ag) with solute atoms in substitution (Ag and Ni, respectively). The results underscore the critical role of anisotropic elasticity in accurately modeling solute segregation. Notably, variations in behavior between grain boundaries having a same tilt angle are only captured when anisotropic elasticity is considered. Furthermore, despite the inherent limitations in addressing non-linear effects at defect cores, the elastic dipole approximation proves to be an effective method for approximating segregation energy spectra in LAGBs obtained through atomistic simulations. Lastly, the estimation of overall solute concentration at grain boundaries highlights the prominent influence of the modulus effect.

**Keywords:** Solute segregation, Anisotropic elasticity, Elastic dipole, Polarizability, Grain boundaries, Dislocations, Molecular dynamics.

## 1. Introduction

Grain boundaries (GBs) play a key role on the mechanical properties of polycrystals [1, 2]. They can act as obstacles to dislocations motion, affecting the material plasticity. They can also interact with solute atoms, being the home of solute segregation or anti-segregation. Solute atoms segregation to GBs can affect intergranular fracture [3, 4], GB migration [5, 6] or can enhance GB disordering due to the premelting effect when the alloy composition and/or temperature approach the solidus line [7]. Intergranular solute segregation is often described through thermodynamic-based models based on the formalism of the Gibbs surface adsorption description [8, 9]. The

<footnote>*Corresponding author</footnote>

<footnote>Email address: thiebaud.richeton@univ-lorraine.fr (Thiebaud Richeton)</footnote>
Langmuir-McLean model [10, 11] was the first to be proposed, offering a description of the inter-
granular solute concentration for diluted binary system. The solute attraction was described using
one single constant segregation energy value, referring to the amount of energy the system is
lowered after solute segregation. This assumption was reexamined by White and Stein [12], con-
cluding that this single energy value was limited to describe the interaction of solute atoms with all
the sites of a GB, introducing afterwards a “per-site” segregation energy. White and Coghlan [13]
proposed a micromechanical model that considered both a size and a modulus effect to estimate
this per-site segregation energy. The size effect was based on the Eshelby’s inclusion result [14].
It was computed as the interaction energy between a solute atom modeled as a spherical inclusion
with a uniform dilatational eigenstrain and the stress field of a GB modeled as an infinite array of
edge dislocations in an infinite isotropic elastic medium. The modulus effect originates from the
influence of a solute atom on the modification of the elastic constants of the region in its vicinity.
It was also computed as an interaction energy with the GB stress field in the White and Coghlan’s
model [13].

Solute atoms can be treated as spherical misfitting inclusions with purely dilatational eigen-
strain [13, 15, 16]. As point defects, solute atoms can also be modeled as a distribution of point-
forces that mimics the forces imposed on the atoms surrounding the defects [17, 18, 19, 20, 21, 22,
23, 24, 25, 26]. The first moment of such a point force distribution is a second-rank tensor known
as the elastic dipole tensor. It was shown that the concept of elastic dipole is actually equivalent
to a uniform eigenstrain within a small inclusion, but only in regions far from the inclusion [25].
Moreover, under an external strain, such as the strain field of a grain boundary, the elastic dipole
is altered. The link between the applied strain and the elastic dipole change can be made by in-
troducing a polarizability tensor, usually denoted α [27, 28, 29, 30, 25]. The latter informs on the
change of the effective elastic constants of the alloy as a function of the solute concentration. Both
the elastic dipole tensor in the absence of external loading and the polarizability tensor can be
determined readily from molecular statics (MS) or ab initio calculations [22, 24, 25, 31]. The latter
are therefore convenient tools to bridge the scales between the atomic description of a point-defect
and the continuum micromechanics based modelling of its interaction with other defects like GBs.

Low-angle symmetric tilt GBs are interfacial defects that can be represented as arrays of edge dislocations, and that commonly exhibit misorientation angles below 15° [9]. In this study, low-angle symmetric tilt GB strain fields are evaluated by MS or by considering edge dislocation walls within the theory of linear elasticity assuming either isotropic elasticity [32] or heterogeneous anisotropic elasticity using the Stroh formalism [33, 34, 35, 36, 37, 38, 39]. The elastic dipole tensors, including the polarizability tensor, are computed by MS simulations for two systems: a solute of Ag in substitution in a matrix of Cu and a solute of Ni in substitution in a matrix of Ni. The same three LAGBs are considered for each system. Then, the interaction energy fields between elastic dipoles and the differently evaluated GB strain fields are determined. The latter are compared to segregation energy fields computed on a site-by-site basis using MS simulations [40].

The paper is organized as follows. Section 2 presents how the GBs are built and how their strain
fields and the segregation energy fields are computed by MS. The micromechanics based modelling
of LAGBs using dislocation walls is described in Section 3. Then, Section 4 briefly recalls the
theory behind the concepts of elastic dipole and polarizability tensors, and how these tensors can
be extracted from atomistic simulations. In Section 5.1, comparisons are first performed between
GB strain fields evaluated either from MS simulations, or from heterogeneous anisotropic/isotropic
elasticity theory. Section 5.2 pursues the discrete/continuum comparisons, but now between the
segregation energy fields computed by MS, and the continuum based interaction energy fields.
Finally, Section 6 proposes a discussion of these results, with a specific focus on the impact of
elastic anisotropy (vs. elastic isotropy) and modulus effect on segregation energies.

## 2. From discrete to continuum modelling methods

### 2.1. Grain boundary generation and characterization

Low angle symmetric tilt grain boundaries (LAGBs) with tilt axis along the [001] direction are considered in this study because they can be readily modeled in the continuum as an infinite array of non-dissociated edge dislocations. Moreover, by considering (010) or (110) GB planes before the rotation by $\theta$ around [001], symmetric tilt GBs can be built for any $\theta$ value, such that $(hkl)_{\text{grain1}} = (hkl)_{\text{grain2}}$ after rotation. The tilt angles $\theta$ can thus be chosen arbitrarily. In the following, GBs with tilt angles of 5° and 10° before relaxation are hence considered.

Using the Atomsk software [41], two fully periodic single crystalline atomistic configurations with identical orientations are created. The crystallographic directions are [100], [010] and [001] along the simulation box directions $x_1, x_2$ and $x_3$, respectively. One system is rotated by $+\frac{\theta}{2}$, while the other is rotated by $-\frac{\theta}{2}$ [42], the rotation being performed around the tilt axis [001] which is parallel to the $x_3$ axis (see Fig. 1). Each system is then cut in half and then merged together. With a periodic computational box, this procedure creates two identical GBs. Then, in-plane translations are performed within the GB plane in order to account for microscopic degrees of freedom in finding the configuration of minimum energy. The energy of each system is minimized using MS with the LAMMPS software (Version 2022-06-23 Update1) [43] using the FIRE algorithm [44] with a force criterion of $10^{-9}$ eV/atom. Once the system with the lowest energy is identified, the simulation box is relaxed in order to lower the residual stresses in the system. This relaxation is performed using the conjugate gradient algorithm, with a force criterion of $10^{-9}$ eV/atom. The configuration with the lowest energy and residual stresses is kept and is assumed to be the more thermodynamically stable. Note that the distance between the two GBs is chosen large enough to neglect their mutual interactions [45].

The mechanical stresses of the atomistic model correspond to the symmetric virial stresses. Per-atom stress tensors [46] are first computed directly from interatomic forces in units of pressure $\times$ volume and are converted in stresses by using the Voronoi volume associated to each atom. The elastic strains are obtained by the calculation of Green-Lagrange strain tensors from a reference FCC configuration, as implemented in the OVITO software [47]. Note that, as the reference configuration has to be a known crystallographic structure and is determined by common neighbor analysis (CNA), it can only be FCC, HCP, BCC or diamond. This implies that no elastic strain can be computed for atoms located at the dislocation core, e.g., the white atoms in Fig. 1.

The systems chosen for our work are FCC copper and FCC silver using the interatomic potentials of William et al. [48] and Pan et al. [49], respectively (see Table 1). The latter was specifically developed to simulate the segregation of Ni atoms at Ag grain boundaries. All the considered GBs are characterized by a periodic arrangement of same defects (see Fig. 1). These defects are identified as being infinite straight edge dislocations, whose Burgers vector is either [010] or [110] (see Fig. 1). The distance between successive dislocations along the $x_1$ axis is denoted $h$. It is noteworthy that, due to the discrete nature of crystallographic planes, dislocations are not necessarily perfectly aligned within the GB, leading to a shift $\delta$ in the $x_2$ direction between two successive dislocations. These GB characteristics, along with the tilt angles $\theta$ before minimization and the GB energies are listed in Table 2.

### 2.2. Computation of the solute segregation energy by molecular statics (MS)

The distribution of segregation energy due to a solute in substitution can be computed through a series of MS calculations, where site-by-site, a solvant atom is substituted by a solute atom [40].
<table><caption>Table 1: Values of the lattice parameter *a* and cubic elastic constants for FCC Cu and FCC Ag, as computed with the potentials.</caption><thead><tr><th>System</th><th>Potential</th><th>a</th><th>C<sub>11</sub></th><th>C<sub>12</sub></th><th>C<sub>44</sub></th><th>Zener ratio</th></tr></thead><tbody><tr><td>Cu</td><td>Williams et al. [48]</td><td>3.615 Å</td><td>169.886 GPa</td><td>122.609 GPa</td><td>76.191 GPa</td><td>3.22</td></tr><tr><td>Ag</td><td>Pan et al. [49]</td><td>4.161 Å</td><td>123.438 GPa</td><td>94.161 GPa</td><td>47.103 GPa</td><td>3.21</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Figure 1: Stable atomistic structures of Cu low-angle grain boundaries. Green and white areas represent FCC and undefined crystallographic structures, respectively. The white ones are defects identified as edge dislocations. Burgers circuits are drawn on the structures to ease the interpretation. In the three cases shown, the dislocations are not perfectly aligned leading to a shift denoted *δ* between two successive dislocations.</figcaption></figure>

<table><caption>Table 2: Characteristics of the different studied GBs: *θ* is the tilt angle before minimization, *b* is the identified Burgers vector normal to the GB plane, *h* is the distance along *x*₁ between two successive dislocations, *δ* is the shift distance along *x*₂ between two successive dislocations and *E*<sup>GB</sup> is the GB energy. The tilt axis direction is [0 0 1] and is parallel to *x*₃ for all the GBs.</caption><thead><tr><th>System</th><th>θ</th><th><b>b</b></th><th>h (Å)</th><th>δ (Å)</th><th>E<sup>GB</sup> (mJ/m<sup>2</sup>)</th></tr></thead><tbody><tr><td>Cu</td><td>5°</td><td>[0 1 0]</td><td>40.7</td><td>1.8</td><td>431</td></tr><tr><td>Cu</td><td>5°</td><td><sup>1</sup>/<sub>2</sub>[1 1 0]</td><td>28.1</td><td>1.25</td><td>338</td></tr><tr><td>Cu</td><td>10°</td><td><sup>1</sup>/<sub>2</sub>[1 1 0]</td><td>15.1</td><td>0.6</td><td>537</td></tr><tr><td>Ag</td><td>5°</td><td>[0 1 0]</td><td>47.5</td><td>0</td><td>318</td></tr><tr><td>Ag</td><td>5°</td><td><sup>1</sup>/<sub>2</sub>[1 1 0]</td><td>32.3</td><td>1.3</td><td>239</td></tr><tr><td>Ag</td><td>10°</td><td><sup>1</sup>/<sub>2</sub>[1 1 0]</td><td>16.3</td><td>0</td><td>366</td></tr></tbody></table>
In order to perform accurate comparisons with the continuum approach (where the interaction energy between a solute and a GB strain field tends to zero at infinite distance from the GB, see further Eq. 31), the following relationship is used to compute the solute segregation energy field:

$$E^{\text{Seg}}(\mathbf{x}) = E^{\text{X}}(\mathbf{x}) - E^{\text{X}}(\infty), \qquad (1)$$

where $E^X(\mathbf{x})$ is the energy of the system containing the GB when a solvant atom $X$ is replaced by a solute atom at the position $\mathbf{x}$. $E^X(\infty)$ is a reference energy value corresponding to the energy of the system containing the GB when a solvant atom $X$ is replaced by a solute atom far from the GB. The energy relaxations are performed using the FIRE algorithm [44], with a force criterion of $10^{-9}$ eV/Å.

### 2.3. From discrete to continuous representation of the elastic fields

In the atomistic model, atoms are evidently not regularly spaced, with an average spacing of 2.5 Å for the systems studied in this work. For the continuum model, we observed that a 0.1 Å grid is required to well capture field variations, especially in the vicinity of GBs. Therefore, in order to be able to make comparisons at the exact same positions with the continuum approach, an interpolation of the atomistic (discrete) quantities on a continuum grid (field variables) should be performed. This new regular grid is created from the information of the atoms surrounding a given node of the regular grid through the use of a Gaussian kernel [50] as follows:

$$\begin{align}
\langle S \rangle &= \frac{\sum_i (S_i \times W_i)}{\sum_i W_i}, \nonumber \\
W_i &= \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(\frac{-d_i^2}{2\sigma^2}\right), \tag{2} \\
d_i &= \sqrt{(x_{1i} - x_1)^2 + (x_{2i} - x_2)^2}. \nonumber
\end{align}$$

In Eq. 2, $\langle S \rangle$ is the interpolated quantity, $W_i$ the weighting factors, $\sigma$ the standard deviation of the Gaussian law, $x_{1i}$ and $x_{2i}$ the coordinates of the atom i and $x_1$ and $x_2$ the coordinates of the considered node in the regular grid. In addition, an interpolation cutoff should be defined. In the following, this cutoff is set to $3\sigma$ while $\sigma$ is fixed at 1 Å for Cu systems and 1.2 Å for Ag systems. Hence, only the contribution of the closest atoms is taken into account in the interpolation (see Fig. 2).

Besides, for comparisons between atomistic and continuum models, one needs also to be able to define a common reference point. The latter is chosen as the point where the interpolated virial stress $\sigma_{12}$ along $x_2$ equals zero (see Fig. 3).

## 3. Continuum representation of low angle grain boundaries (LAGB)

### 3.1. Elastic fields in linear anisotropic elasticity

It is well-known that low-angle symmetrical tilt grain boundaries can be modelled by periodic arrays of edge dislocations [32]. An infinite wall of edge dislocations spaced along the $x_1$ axis is considered in Fig. 4. The Burgers vector lies along the $x_2$ axis, the dislocation line is along the $x_3$ axis, and the GB plane corresponds the $(x_1, x_3)$ plane. The spacing distance between dislocations is denoted $h$. This dislocation wall separates two anisotropic materials or crystals, which differs only by their stiffness tensor: $C^I$ and $C^{II}$ in the region $x_2 > 0$ and $x_2 < 0$, respectively. The
<figure><img src="image_3.png" /><figcaption>Figure 2: Interpolation of the volumetric strain field of a Cu symmetric tilt grain boundary ($\theta = 5^\circ$, $\mathbf{b} = [001]$) using a Gaussian law with $\sigma = 1$ Å, $\mu = 0$. The interpolation cutoff is set to $3\sigma$. The full width at half maximum of the Gaussian law is equal to $2.35$ Å.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Contour and plot profile of the $\sigma_{12}$ component for a Cu symmetric tilt grain boundary ($\theta = 5^\circ$, $\mathbf{b} = [010]$).</figcaption></figure>
interface is thus planar with its normal along the $x_2$ axis, to be consistent with Fig. 1. This interface is assumed perfectly bonded (i.e., satisfying the continuity of displacement and traction vectors).

<figure><img src="image_5.png" /><figcaption>Figure 4: Symmetric tilt grain boundary modeled as an infinite periodic array of edge dislocations.</figcaption></figure>

First of all, we recall how the elastic field of a straight dislocation in such an anisotropic bi-material can be computed using the Stroh's formalism in linear elasticity [34, 35, 36, 37, 38, 39]. For a dislocation located at point $(X_1, X_2)$, the displacement gradient components at point $(x_1, x_2)$ can be written as [51]:

$$ 
\begin{aligned}
u_{i,1} &= 2\mathrm{Re}(A_{ij}g_j(z_j)), \\
u_{i,2} &= 2\mathrm{Re}(A_{ij}p_jg_j(z_j)),
\end{aligned} 
\qquad (3) 
$$

where $z_j = x_1 + p_j x_2$ and $\mathbf{A}$ is a $3 \times 3$ matrix containing the eigenvectors related to the eigenvalues $p_j$ with positive part of the Stroh sextic equation:

$$ \det(C_{i1k1} + p(C_{i1k2} + C_{i2k1}) + p^2 C_{i2k2}) = 0. \qquad (4) $$

The function vector $\mathbf{g}$ is given by:

$$ g_j(z_j) = \begin{cases} g_j^0(z_j) + \overline{V_{ji}^{I,II}} \overline{g_i^0(z_j)} & \text{if } x_2 > 0, \\ W_{ji}^{I,II} g_i^0(z_j) & \text{if } x_2 < 0, \end{cases} \qquad (5) $$

where:

$$ g_i^0(z_j) = \frac{\mathbf{q}^0}{z_j - s_i}, \qquad (6) $$

$$ \mathbf{q}^0 = -\frac{1}{2\pi} i (\mathbf{B}^T \mathbf{b}), \qquad (7) $$

$$ B_{ij} = (C_{i2k1} + p_j C_{i2k2}) A_{kj}, \qquad (8) $$

and $s_j = X_1 + p_j X_2$ is related to the dislocation position. The tensors $\mathbf{V}^{I,II}$ and $\mathbf{W}^{I,II}$ are defined as [36, 38, 39]:

$$ 
\begin{aligned}
\mathbf{V}^{I,II} &= (\mathbf{B}^{II} \mathbf{A}^{II-1} \overline{\mathbf{A}^I} - \overline{\mathbf{B}^I})^{-1} (\mathbf{B}^I - \mathbf{B}^{II} \mathbf{A}^{II-1} \mathbf{A}^I) \text{ and} \\
\mathbf{W}^{I,II} &= (\overline{\mathbf{B}^I} \mathbf{A}^{II-1} \mathbf{A}^I - \mathbf{B}^{II})^{-1} (\overline{\mathbf{B}^I} \mathbf{A}^{II-1} \mathbf{A}^I - \mathbf{B}^I),
\end{aligned} 
\qquad (9) $$
where $A^I$ and $B^I$ refer to the resolution of the Stroh eigenrelation in material I, while $A^{II}$ and $B^{II}$ refer to material II.

Then, for a periodic array of dislocations and since linear elasticity holds, the elastic contributions of all the dislocations located at $s_j^n = nh + p_j X_2$ are simply added, $n$ being an integer going from $-\infty$ to $+\infty$. Hence, both displacement gradient components become:

$$ 
\begin{aligned}
u_{i,1} &= 2\mathrm{Re} \left( A_{ij} \sum_{n=-\infty}^{\infty} g_j^n(z_j) \right), \\
u_{i,2} &= 2\mathrm{Re} \left( A_{ij} \sum_{n=-\infty}^{\infty} p_j g_j^n(z_j) \right),
\end{aligned} 
\qquad (10) 
$$

where:

$$ g_j^n(z_j) = \begin{cases} g_j^{0n}(z_j) + \overline{V_{ji}^{I,II}} g_i^{0n}(z_j) & \text{if } x_2 > 0, \\ W_{ji}^{I,II} g_i^{0n}(z_j) & \text{if } x_2 < 0, \end{cases} \qquad (11) $$

and:

$$ g_i^{0n}(z_j) = \frac{q_i^0}{z_j - s_i^n}. \qquad (12) $$

This sum can be easily performed with the help of the following mathematical formula [32]:

$$ \sum_{-\infty}^{\infty} \frac{1}{n + a} = \pi \cot \pi a. \qquad (13) $$

(14)

Accordingly, the sum $\sum_{n=-\infty}^{\infty} g_i^{0n}(z_j)$ is computed as:

$$ 
\begin{aligned}
\sum_{-\infty}^{\infty} g_i^{0n}(z_j) &= q_i^0 \sum_{-\infty}^{\infty} \frac{1}{x_1 + p_j x_2 - nh - p_j X_2} \\
&= -\frac{q_i^0}{h} \sum_{-\infty}^{\infty} \frac{1}{n + \frac{p_j X_2 - p_j x_2 - x_1}{h}} \\
&= -\frac{q_i^0}{h} \pi \cot \frac{\pi (p_j X_2 - p_j x_2 - x_1)}{h}.
\end{aligned} 
\qquad (15) $$

It is important to underline that the latter expression cannot be factorized by $p_j$ since $p_j$ can refer to material I and/or material II. For instance, if we want to compute the elastic fields in material II for an array of dislocations located in material I, we should consider the expression $p_j^I X_2 - p_j^{II} x_2 - x_1$ where $p_j^I$ contains the eigenvalues of material I whereas $p_j^{II}$ contains the eigenvalues of material II. Finally, the strain, the rotation and Cauchy stress tensors are obtained from the classical relations in small strain setting:
$$ \varepsilon_{ij} = \frac{1}{2} (u_{i,j} + u_{j,i}), \qquad (16) $$

$$ \omega_{ij} = \frac{1}{2} (u_{i,j} - u_{j,i}), \qquad (17) $$

$$ \sigma_{ij} = C_{ijkl}\varepsilon_{kl}. \qquad (18) $$

In homogeneous isotropic elasticity, the expressions of the stress fields of such a dislocation wall are given by [32] :

$$ \begin{aligned} \sigma_{11} &= \sigma_0 \sin 2\pi\lambda (\cosh 2\pi\xi - \cos 2\pi\lambda - 2\pi\xi \sinh 2\pi\xi) \\ \sigma_{22} &= \sigma_0 \sin 2\pi\lambda (\cosh 2\pi\xi - \cos 2\pi\lambda + 2\pi\xi \sinh 2\pi\xi) \\ \sigma_{33} &= 2\nu\sigma_0 \sin 2\pi\lambda (\cosh 2\pi\xi - \cos 2\pi\lambda) \\ \sigma_{12} &= -2\pi\xi\sigma_0 (\cosh 2\pi\xi \cos 2\pi\lambda - 1) \\ \sigma_{31} &= 0 \\ \sigma_{23} &= 0 \end{aligned} \qquad (19) $$

where $\xi = x_2/h$, $\lambda = x_1/h$ and $\sigma_0 = \frac{\mu b}{2h(1-\nu)(\cosh 2\pi\xi - \cos 2\pi\lambda)^2}$. $\mu$ is the shear modulus and $\nu$ is the Poisson's ratio. Strains are obtained using the Hooke's law in isotropic elasticity.

## 3.2. Computation details

The computations of the elastic fields of the infinite dislocation walls are performed using the MATLAB (R2020a) software for both the isotropic and the anisotropic cases. For the isotropic case, the values of the shear modulus and the Poisson's ratio are obtained from a Voigt-Reuss-Hill average [52, 53] of the elastic constants associated to the used interatomic potential [48, 49] (cf. Table 1) in order to carry out relevant comparisons between continuum and discrete approaches. The values of these isotropic constants are given in Table 3.

<table><caption>Table 3: Shear modulus $\mu$ and Poisson's ratio $\nu$ obtained from the Voigt-Reuss-Hill average [52, 53] of the elastic constants associated to the interatomic potentials used in the atomistic simulations.</caption><thead><tr><th>System</th><th>μ (GPa)</th><th>ν</th></tr></thead><tbody><tr><td>Cu</td><td>47.750</td><td>0.345</td></tr><tr><td>Ag</td><td>29.539</td><td>0.370</td></tr></tbody></table>

For the anisotropic case, the used elastic constants are directly those of the atomic potentials (see Table 1). The following procedure is used to set the crystallographic orientations of both grains. First, the infinite dislocation wall is considered in an initial frame matching the crystallographic directions of the dislocations, i.e., with the Burgers vectors normal to the GB plane and the dislocation lines parallel to the tilt axis. The (elastic) rotation fields are first computed in homogeneous anisotropic elasticity following the equations of Section 3.1 with $C^I = C^{II}$. Checking that the obtained rotations can be consistently approximated by uniform values in each crystal (see Fig. 5), values far from the dislocations wall at $x_2 = 5h$ are considered to rotate $C^I$ and $C^{II}$. Then, the rotation fields are computed once again considering the new tensors $C^I$ and $C^{II}$. The procedure is repeated until convergence assuming a relative tolerance in rotations less than $10^{-5}$.
Finally, all the elastic fields are computed considering the last updated values of $C^I$ and $C^{II}$ after convergence.

Besides, in some cases showed in Table 2, a shift $\delta$ along $x_2$ between two successive dislocations is observed in the MS simulations (cf. Fig. 1 and Table 2). This shift is taken into account in the continuum based approach to compute the elastic fields by first evaluating the rotations of $C^I$ and $C^{II}$ with zero shift and then by replacing the infinite dislocations wall with separation distance $h$ located at $x_2 = 0$ by two infinite dislocations walls with separation distance $2h$ located at $x_2 = -\frac{\delta}{2}$ and $x_2 = +\frac{\delta}{2}$.

<figure><img src="image_6.png" /><figcaption>Figure 5: Rotation profiles perpendicular to the three Cu GBs obtained by the continuum approach in heterogeneous anisotropic elasticity.</figcaption></figure>

## 4. Elastic dipole and interaction energy

Point defects like solute atoms or vacancies are very often modeled as a distribution of point-forces, i.e. like a force multipole that mimics the forces imposed on the atoms surrounding the defects [17, 18, 19, 20, 21, 23, 24, 25, 26]. The first moment of the point force distribution is known as the elastic dipole $P$. When applying an external strain (e.g. the strain field of a GB), the elastic dipole tensor depends on the applied strain as follows [29, 25, 31]:

$$P_{ij}(\varepsilon^{ext}) = P_{ij}^{0} + \alpha_{ijkl}\varepsilon_{kl}^{ext}, \qquad (20)$$

where $P_{ij}^0$ is the permanent elastic dipole tensor in the absence of any external applied strain field [25, 24]. The polarizability tensor $\alpha$ characterizes how the elastic dipole tensor is altered by an applied external strain field. It describes how a given concentration of point defects modifies the elastic response of the alloy (see further Eqs. 26 and 28).

### 4.1. Permanent elastic dipole determination from atomistic simulations

Different techniques can be used to identify the permanent elastic dipole tensors from atomistic simulations [25, 24]. Using MS simulations, the elastic dipole tensor can be determined from the residual stress tensor, as proposed by Clouet et al. [25, 24]. In this approach, a bulk FCC crystal in a fully periodic cubic simulation box of volume $V$ with no residual stress is first created. Then, a point defect is introduced in the simulation box. It is noteworthy that the dimensions of the box should be large enough to neglect the interactions between solutes due to periodic boundary
conditions. The permanent elastic dipole tensor is then obtained from the homogeneous stress tensor averaged over the simulation box, $\langle \sigma \rangle$, and weighted by $V$ [25, 24] as follows:

$$P_{ij}^{0} = -V\langle\sigma_{ij}\rangle. \qquad (21)$$

In the last equation, $\langle \sigma \rangle$ corresponds to the average of the stress tensor over the periodic box, induced by the solute atom placed inside the box in the absence of applied strain. Besides, $P^0$ can be related to an eigenstrain tensor, which is another widespread approach to model point defects. Indeed, solute atoms can also be considered as misfitting spherical inclusions with purely dilatational eigenstrain [13, 15, 16]. In the limit of an infinitesimal inclusion and far from the inclusion, Clouet et al. [25] have shown that the permanent elastic dipole tensor $P^0$ is related to the average eigenstrain tensor $\bar{\epsilon}^*$ over the inclusion's volume $\Omega$ by (see also [54]):

$$P_{ij}^{0} = \Omega C_{ijkl} \bar{\epsilon}_{kl}^{*}. \qquad (22)$$

Thus, the equivalent average eigenstrain tensor $\bar{\epsilon}^*$ can be computed as:

$$\bar{\epsilon}_{ij}^{*} = -\frac{V}{\Omega} S_{ijkl} \langle \sigma_{kl} \rangle. \qquad (23)$$

Considering an Ag atom in substitution in a Cu matrix and a Ni atom in substitution in an Ag matrix, denoted respectively “Cu/Ag” and “Ag/Ni”, the described method provides the following values for the components of the permanent elastic dipole tensor (in eV) and the eigenstrain tensors:

$$P^{0,Cu/Ag} = \begin{pmatrix} 5.51 & 0 & 0 \\ 0 & 5.51 & 0 \\ 0 & 0 & 5.51 \end{pmatrix}, \qquad P^{0,Ag/Ni} = \begin{pmatrix} -5.14 & 0 & 0 \\ 0 & -5.14 & 0 \\ 0 & 0 & -5.14 \end{pmatrix}, \qquad (24)$$

$$\bar{\epsilon}^{*Cu/Ag} = \begin{pmatrix} 0.10 & 0 & 0 \\ 0 & 0.10 & 0 \\ 0 & 0 & 0.10 \end{pmatrix}, \qquad \bar{\epsilon}^{*Ag/Ni} = \begin{pmatrix} -0.15 & 0 & 0 \\ 0 & -0.15 & 0 \\ 0 & 0 & -0.15 \end{pmatrix}. \qquad (25)$$

Here, both the permanent elastic dipole and the eigenstrain tensors are isotropic (only diagonal terms with a same value), as expected for cubic lattices [23]. This means that these solutes in substitution induce size effects but no shape effect. Evidently, these size effects are not sensitive to the crystallographic orientations of the host lattice.

## 4.2. Polarizability tensor determination from atomistic simulations

The procedure to identify the polarizability tensor $\alpha$ from MS simulations is now summarized. Note that it is similar to the one used to evaluate the elastic constants of the material modeled by an interatomic potential, but with a system containing a solute atom. A bulk crystal in a fully periodic simulation box containing a single point defect is strained. Six different symmetric applied strain tensors with only one non-zero component each time are considered, along with different strain values (e.g., from $\epsilon_{11}^{ext} = 0.01\%$ to $\epsilon_{11}^{ext} = 0.1\%$). It's important to ensure that calculations remain within the linear elastic regime and that forces and energy are minimized for each simulations. The homogeneous stress tensor averaged over the simulation box $\langle \sigma \rangle$ is measured each time and the components of $\alpha$ are obtained from linear regressions according to the modification of Eq. 21 in the presence of applied strain [25]:

$$P_{ij} = P_{ij}^{0} + \alpha_{ijkl} \epsilon_{kl}^{ext} = V (C_{ijkl} \epsilon_{kl}^{ext} - \langle \sigma_{ij} \rangle). \qquad (26)$$
<figure><img src="image_7.png" /><figcaption>Figure 6: Variation of $P_{11}$ with ten different $\varepsilon_{22}$ values applied in a Cu matrix with a substitutional Ag solute. The $\alpha_{1122}$ coefficient of the polarizability tensor is the slope value.</figcaption></figure>

For example, determining the $\alpha_{1122}$ component implies to consider $P_{11}$ as a linear function of $\varepsilon_{22}^{ext}$ strain (see Fig. 6). Besides, a simple rearrangement of Eq. 26 yields:

$$ \langle \sigma_{ij} \rangle = \left( C_{ijkl} - \frac{1}{V} \alpha_{ijkl} \right) \varepsilon_{kl}^{ext} - \frac{1}{V} P_{ij}^{0} = C_{ijkl}^{eff} \varepsilon_{kl}^{ext} - \frac{1}{V} P_{ij}^{0}, \quad (27) $$

and thus:

$$ \alpha_{ijkl} = V \left( C_{ijkl} - C_{ijkl}^{eff} \right). \quad (28) $$

$C^{eff}$ denotes the effective elastic stiffness tensor of the alloy as a function of the solute concentration [25, 31].

It is important to underline that it is not trivial to apply pure shear strain in MS simulations. Indeed, only a simple shear deformation is straightforward. In order to achieve pure shear deformation conditions, positive and negative strains are applied for each component and the final polarisability tensor is computed as the average of the two intermediate tensors, which annihilates the undesired rotations inherent to simple shear deformation. Using Voigt's notation, the polarizability tensors in the crystal's frame for both studied systems are (in eV):

$$ \alpha^{Cu/Ag} = \begin{pmatrix} -12.06 & -13.98 & -13.98 & 0 & 0 & 0 \\ -13.98 & -12.06 & -13.98 & 0 & 0 & 0 \\ -13.98 & -13.98 & -12.06 & 0 & 0 & 0 \\ 0 & 0 & 0 & -5.62 & 0 & 0 \\ 0 & 0 & 0 & 0 & -5.62 & 0 \\ 0 & 0 & 0 & 0 & 0 & -5.62 \end{pmatrix} \quad \text{and} \quad (29) $$

$$ \alpha^{Ag/Ni} = \begin{pmatrix} 35.95 & 19.21 & 19.21 & 0 & 0 & 0 \\ 19.21 & 35.95 & 19.21 & 0 & 0 & 0 \\ 19.21 & 19.21 & 35.95 & 0 & 0 & 0 \\ 0 & 0 & 0 & 15.17 & 0 & 0 \\ 0 & 0 & 0 & 0 & 15.17 & 0 \\ 0 & 0 & 0 & 0 & 0 & 15.17 \end{pmatrix}. \quad (30) $$
Unlike the permanent dipole tensor, these tensors are not isotropic and depend on the crystallographic orientation of the host lattice. Hence, they must be rotated in the same coordinates frame as the elastic fields of the GB to compute the interaction energy (see section 4.3).

### 4.3. Interaction energy between an elastic dipole defect and a GB strain field

Following Schober [29] and Clouet et al. [25], the elastic interaction energy between a point-defect located at position $\mathbf{x}$ and the strain field $\boldsymbol{\varepsilon}$ of a GB is at the first order:

$$E^{\text{int}}(\mathbf{x}) = -P_{ij}^{0}\varepsilon_{ij}(\mathbf{x}) - \frac{1}{2}\alpha_{ijkl}\varepsilon_{ij}(\mathbf{x})\varepsilon_{kl}(\mathbf{x}). \qquad (31)$$

Considering the isotropic permanent elastic dipole tensor for cubic lattices, the first term $-P_{ij}^{0}\varepsilon_{ij}(\mathbf{x})$ simplifies to $-\frac{1}{3}P_{kk}\varepsilon_{kk}(\mathbf{x})$. Hence, it can be identified as a size effect since it is only influenced by the volumetric strain field $\varepsilon_{kk}(\mathbf{x})$. As originally defined, the second term related to the polarizability tensor leads to a modulus effect. Considering the cubic symmetry of the polarizability tensor $\boldsymbol{\alpha}$ and using the Walpole's decomposition for cubic crystals [55, 56], the following expression of $E^{\text{int}}(\mathbf{x})$ can be readily derived:

$$E^{\text{int}}(\mathbf{x}) = -\frac{1}{3}P_{kk}\varepsilon_{kk}(\mathbf{x}) - \frac{1}{2}\alpha_K\varepsilon_{kk}^2(\mathbf{x}) - \alpha_{\mu'}(\varepsilon_{11}^{\prime 2}(\mathbf{x}) + \varepsilon_{22}^{\prime 2}(\mathbf{x}) + \varepsilon_{33}^{\prime 2}(\mathbf{x})) \\ -2\alpha_{\mu''}(\varepsilon_{12}^{\prime 2}(\mathbf{x}) + \varepsilon_{31}^{\prime 2}(\mathbf{x}) + \varepsilon_{23}^{\prime 2}(\mathbf{x})), \qquad (32)$$

where the components of $\boldsymbol{\alpha}$ and $\boldsymbol{\varepsilon}$ are expressed in the crystal's frame, $\boldsymbol{\varepsilon}'$ is the deviatoric part of $\boldsymbol{\varepsilon}$ while $\alpha_K$, $\alpha_{\mu'}$ and $\alpha_{\mu''}$ are defined similarly as their equivalent of the cubic stiffness tensor $\mathbf{C}$, i.e. the bulk modulus $K$ and shear moduli $\mu'$ and $\mu''$, reflecting the $(001)[110]$ and $(001)[100]$ shear resistances, respectively [56]. Using Voigt's notation, these polarizability moduli hence write:

$$\alpha_K = \frac{\alpha_{11} + 2\alpha_{12}}{3}, \quad \alpha_{\mu'} = \frac{\alpha_{11} - \alpha_{12}}{2}, \quad \alpha_{\mu''} = \alpha_{44}. \qquad (33)$$

For an isotropic polarizability tensor, $\alpha_{\mu'} = \alpha_{\mu''} = \alpha_{\mu}$ and the elastic interaction energy becomes:

$$E^{\text{int}}(\mathbf{x}) = -\frac{1}{3}P_{kk}\varepsilon_{kk}(\mathbf{x}) - \frac{1}{2}\left(\alpha_K\varepsilon_{kk}^2(\mathbf{x}) + 2\alpha_\mu\varepsilon_{ij}^{\prime}(\mathbf{x})\varepsilon_{ij}^{\prime}(\mathbf{x})\right). \qquad (34)$$

Considering in addition isotropic elasticity, the interaction energy can then be written as:

$$E^{\text{int}}(\mathbf{x}) = -\frac{1}{3}\Omega\delta^*\sigma_{kk}(\mathbf{x}) - \frac{1}{2}\left(\frac{\alpha_K}{9K^2}\sigma_{kk}^2(\mathbf{x}) + \frac{\alpha_\mu}{2\mu^2}\sigma_{ij}^{\prime}(\mathbf{x})\sigma_{ij}^{\prime}(\mathbf{x})\right), \qquad (35)$$

where $\delta^* = \bar{\varepsilon}_{kk}^*$ (see Eq. 23), $\boldsymbol{\sigma}$ is the GB stress field associated to $\boldsymbol{\varepsilon}$, and $\boldsymbol{\sigma}'$ its deviatoric part. $K$ and $\mu$ are the isotropic bulk and shear moduli of the host material, respectively. It is noteworthy that Eq. 35 is actually equivalent to the binding energy (i.e., the negative of the interaction energy) described by White and Coghlan [13], except of an inverse sign convention for the isotropic polarizability constants. Following [31], the values of $\alpha_\mu$ are obtained considering a Voigt average:

$$\alpha_{\mu} = \frac{2}{5}\alpha_{\mu'} + \frac{3}{5}\alpha_{\mu''}. \qquad (36)$$

The determined values of $\alpha_K$, $\alpha_{\mu'}$, $\alpha_{\mu''}$ and $\alpha_\mu$ for both systems studied are reported in Table 4.
<figure><img src="image_8.png" /><figcaption>Figure 7: Comparisons of volumetric strain fields for the three grain boundaries presented in Section 2.1 in Cu (a) and Ag (b). Three different ways of computation are considered: from interpolation of OVITO fields after relaxation by MS and from continuum approaches based on walls of edge dislocations, either in linear heterogeneous anisotropic elasticity or in isotropic elasticity. The box size is *h* along *x*₁ and 2*h* along *x*₂. Thus, it must be noticed that the box size changes according to the GB and the system.</figcaption></figure>
<table><caption>Table 4: Values of the bulk and shear polarizability moduli for both systems studied (see text).</caption><thead><tr><th>System</th><th>α<sub>K</sub> (eV)</th><th>α<sub>μ</sub>' (eV)</th><th>α<sub>μ</sub>'' (eV)</th><th>α<sub>μ</sub> (eV)</th></tr></thead><tbody><tr><td>Cu/Ag</td><td>-13.34</td><td>0.96</td><td>-5.62</td><td>-2.99</td></tr><tr><td>Ag/Ni</td><td>24.79</td><td>8.37</td><td>15.17</td><td>12.45</td></tr></tbody></table>

# 5. Results

## 5.1. Comparisons of LAGB strain fields

First of all, the volumetric strain fields $\varepsilon_{kk}$ of the three GBs presented in Section 2.1 are compared in Fig. 7 according to three different ways of computation: using the finite Green-Lagrange strain tensors from atomistic structure (see Section 2.1) and from small strain continuum approaches considering walls of edge dislocations, either in linear heterogeneous anisotropic elasticity (see Section 3.1) or in isotropic elasticity [32]. The fields extracted from atomistic simulations are interpolated using the method described in Section 2.3 on a grid with a spatial resolution of 0.1 Å. Comparisons are made for both Cu and Ag systems (see Table 1). It is noteworthy that linear elastic fields are singular at the dislocation positions, which leads to an overestimation of strain values within the dislocation cores. On the contrary, atomistic strain values are zero within the dislocation cores since no elastic strain information can be obtained at these points, as explained in Section 2.1 (see Figure 2). Therefore, comparisons should be made at a minimum distance of 4 Å away from the dislocation positions.

Figure 7 shows the 2D contour plots of $\varepsilon_{kk}$ obtained from the three methods of computation. Contour shapes obtained in heterogeneous anisotropic elasticity are very much alike to those interpolated from atomistics. In particular, the change of shape of $\varepsilon_{kk}$ fields between the two GBs having the same tilt angle value ($\theta = 5^{\circ}$) is well reproduced by the continuum approach in anisotropic elasticity. On the contrary, isotropic elasticity cannot capture such differences between GBs having a same tilt angle. Although $b$ and $h$ values are different between the two GBs with $\theta = 5^{\circ}$, the contour shapes of $\varepsilon_{kk}$ computed in isotropic elasticity are almost homothetic with same rounded shapes (Fig. 7). Minor discrepancies arise however because of different values of the shift distance $\delta$ (see Table 2).

Furthermore, profiles are also plotted perpendicular and parallel to GBs in Fig. 8. From the profiles perpendicular to GBs, it is clear that the strains computed in heterogeneous anisotropic elasticity better fit the strains from atomistics than the strains computed in isotropic elasticity. In particular, the perpendicular profiles obtained by MS show positive humps of $\varepsilon_{kk}$ in the transition area from minimum to zero values for the $\theta = 5^{\circ}$ and $b = [0 1 0]$ GB. The latter are quite well reproduced by the continuum-based approach in anisotropic elasticity. By comparison, the approach in isotropic elasticity exhibits similar smooth transitions from minimum to zero values for both the $\theta = 5^{\circ}$ and $b = [0 1 0]$ GB and the $\theta = 5^{\circ}$ and $b = \frac{1}{2}[1 1 0]$ GB. Nevertheless, profiles parallel to GBs show that isotropic elasticity is as good as anisotropic elasticity in the exact middle of GBs and sufficient as well to capture atomistic strains along the GB outside the dislocation cores.

## 5.2. Comparisons between segregation energy and interaction energy fields

Using the method previously described in Section 2.2, the solute segregation energy fields $E^{seg}(\boldsymbol{x})$ are computed for the 6 configurations considered (3 GBS and 2 systems: Cu/Ag and Ag/Ni). Figure 9 displays the fields obtained for the Cu/Ag system. Segregation is favored when the energy is negative (blue areas) while anti-segregation is promoted when the energy is positive (red areas). Atomistic segregation energy fields are interpolated using the procedure presented
<figure><img src="image_9.png" /><figcaption>Figure 8: Volumetric strain profiles perpendicular to GB at 5 Å from a dislocation (see inset) (a) and along the middle of the GB (see inset) (b).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Segregation energy fields (Eq.1) for the 3 grain boundaries presented in Section 2.1 in a Cu/Ag system.</figcaption></figure>

in Section 2.3 in order to make comparisons with approaches based on linear elasticity, i.e. the computation of the interaction energy $E^{int}(\mathbf{x})$ (see Eq.31). The comparisons between $E^{seg}(\mathbf{x})$ and $E^{int}(\mathbf{x})$, where the strains are computed according to the three different methods (see Section 5.1) are shown in Figs. 10 and 11. For the continuum-based approach in isotropic elasticity, an isotropic polarizability tensor is also considered (see Eqs. 34 and 36). The permanent elastic dipole and the polarizability tensors used in the computation of the interaction energy are always those obtained directly by molecular statics (see Eqs. 24, 29, 30 and Table 4 for the approach in isotropic elasticity).

The segregation and anti-segregation domains associated to $E^{seg}$ fields, respectively in blue and yellow in Fig. 10, are rather well-reproduced for all the GBs of both systems by using the $E^{int}$ fields, when strains are taken from atomistics or from the continuum approach based on heterogeneous anisotropic elasticity. As expected, some discrepancies can however be noticed close to the dislocation singularities (i.e. at less than 5 Å from the dislocation center) between full atomistic calculations and the methods based on linear elasticity. Besides, the segregation and anti-segregation domains are much less well-described by the continuum approach based on isotropic elasticity. These conclusions are even more obvious when analyzing the profiles perpendicular to the two $\theta = 5^\circ$ GBs in Fig. 11. For the $\theta = 5^\circ$ and $\mathbf{b} = [010]$ GB in the Cu/Ag system, the $E^{seg}$ profile shows a strong anti-segregation zone (i.e., positive energies) which is surrounded by segregation zones (i.e. negative energies). This feature is very well-captured by the $E^{int}$ fields when strains are taken from atomistics or from the continuum approach based on heterogeneous anisotropic elasticity. On the contrary, the approach based on isotropic elasticity displays a positive energy all along the profile. The reverse effect is observed for the Ag/Ni system, a strong segregation zone surrounded by anti-segregation zones, which is again well-captured by the linear elasticity approaches when strains are computed by atomistics or in anisotropic elasticity whereas a negative energy is observed all along the profile in isotropic elasticity. For the $\theta = 5^\circ$ and $\mathbf{b} = \frac{1}{2}[110]$ GBs, the $E^{seg}$ profiles show smooth transitions from maximal (for the Cu/Ag system) or minimal energy (for the Ag/Ni system) to zero energy zones far from the GB with no change of sign. These smooth transitions are again very well-reproduced by the linear elasticity approaches when strains
<figure><img src="image_11.png" /><figcaption>Figure 10: Comparisons between the segregation energy fields $E^{seg}$ (Eq.1) after interpolation and the interaction energy fields $E^{int}$ (Eq.31) for Cu/Ag system (a) and for Ag/Ni system (b). For the latter, three different methods are considered to compute the strains: from interpolation of OVITO fields after relaxation by MS, and from continuum approaches based on walls of edge dislocations, either in linear heterogeneous anisotropic elasticity or in isotropic elasticity. The permanent elastic dipole and the polarizability tensors are always those obtained directly by MS simulations unless for the shear polarizability $\alpha_\mu$ used in isotropic elasticity for which a Voigt average is performed (see Eq. 36).</figcaption></figure>
are computed by atomistics or in anisotropic elasticity. On the contrary, the approach based on isotropic elasticity exhibit much sharper transitions. Finally, it is noteworthy that the match between the $E^{seg}$ profiles and the $E^{int}$ profiles is only slightly better when strains are taken from atomistics as compared to the case when they are directly computed in linear anisotropic elasticity.

<figure><img src="image_12.png" /><figcaption>Figure 11: Segregation energy and interaction energy profiles perpendicular to the grain boundaries at 5 Å from a dislocation (see inset) for the 3 grain boundaries presented in Section 2.1 and the two systems: Cu/Ag (a) and Ag/Ni (b).</figcaption></figure>

### 5.3. Spectrum of segregation energies

Figure 12 displays the segregation and interaction energy spectra for the 3 GBs of the Ag/Ni system. The histograms based on segregation energies (Fig. 12, first column) represent the fraction of atomic positions that have a segregation energy $E^{seg}$ within a specific range. Such energy spectra are directly computed from segregation maps similar to Fig. 9 for the Ag/Ni system, except that the dimensions of the box are fixed at one $h$ period along $x_1$ and at two times the greatest value of $h$ for all the systems and GBs considered, i.e. $2 \times 47.5$ Å (see Table 2), along $x_2$. As a result, relevant comparisons in terms of atomic positions fraction between the different GBs can be performed. The histograms based on interaction energies (Fig. 12, second column) represent the fraction of nodes of the interpolation grid (step size of 0.1 Å) that have an interaction energy $E^{int}$ within a specific range. All histograms have a bin size of 0.02 eV. For more clarity, the values close to zero (from -0.06 eV to 0.06 eV) have been removed from the plots, since only the highest and lowest energy values are expected to affect anti-segregation and segregation, respectively.

All the different approaches reproduce the strong asymmetry of the 3 GBs which exhibit much more areas favorable to segregation than to anti-segregation. The continuum approaches based either on isotropic or anisotropic elasticity overestimate the fraction of nodes with very high energy values due to the consideration of singularities at the dislocation positions. On the contrary, the
<figure><img src="image_13.png" /><figcaption>Figure 12: Segregation and interaction energy spectra for the 3 GBs of the Ag/Ni system. The width of the box centered on the GB that is used to evaluate the fraction of atomic positions or nodes is fixed at 95 Å.</figcaption></figure>

approach based on the interaction energy with strains computed by atomistics underestimate the
fraction of nodes with very high energy values because no value of strains are available in the
dislocation cores.

**6. Discussion**

6.1. *Impact of modulus effect on segregation energies*

So far, both size and modulus effects have been considered together to evaluate the interaction
energy fields E<sup>int</sup>(x) from Eq. 31. Both effects can be easily analyzed separately due to the
additive nature of their contributions. Figure 13 shows contour plots of the size and modulus
effect contributions to E<sup>int</sup>(x) for a same GB in both studied systems. For cubic lattices, the
size effect, -1/3 P<sub>kk</sub>ε<sub>kk</sub>(x), favor both segregation and anti-segregation zones as it scales with the
volumetric strain field arising from arrays of edge dislocations. On the contrary, the expression
of the interaction energy for cubic lattices given by Eq. 32 shows that the modulus effect will in
general favor only segregation or anti-segregation. This is definitely the case when α<sub>K</sub>, α'<sub>μ</sub> and
α''<sub>μ</sub> have all the same sign, as for the Ag/Ni system (see Table 4), since the latter coefficients are
multiplied by squared strain components. Hence, for the Ag/Ni system, the modulus effect only
promotes segregation. Nevertheless, it might occur that α<sub>K</sub>, α'<sub>μ</sub> and α''<sub>μ</sub> have different signs, as
for the Cu/Ag system (see Table 4), in which case the situation might become more intricate. For
the Cu/Ag system, α'<sub>μ</sub> which is positive, is however much smaller in magnitude than α<sub>K</sub> and α''<sub>μ</sub>
which are negative. Hence, the modulus effect should promote only segregation for the Cu/Ag
system as it is actually observed in Fig. 13. It should also be noticed that the values of α<sub>K</sub>, α'<sub>μ</sub>
and α''<sub>μ</sub> are much higher in magnitude for the Ag/Ni system than for the Cu/Ag system which is
consistent with the contours of Fig. 13 showing a negligible modulus effect for the Cu/Ag system
and a small but significant effect for the Ag/Ni system. It is noteworthy also that most of the
<figure><img src="image_14.png" /><figcaption>Figure 13: Contour plots of the size and modulus effect contributions on the interaction energy fields with strains obtained by MS for two different systems: (a) Cu/Ag and (b) Ag/Ni.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: Segregation energy and interaction energy profiles along $x_2$ direction at 5 Å from a dislocation for two different GBs in a Ag/Ni system, with and without the consideration of the modulus effect.</figcaption></figure>
polarizability coefficients are much lower than the value of 20 eV tested by White and Coghlan in their original paper [13]. Accordingly, in the present study, the size effect appears as the main effect controlling the solute segregation for both systems. Besides, Fig. 13 indicates that the contribution of the modulus effect seems to increase with the tilt angle of the grain boundary.

In addition, Fig. 14 shows profiles of $E^{\text{int}}(\mathbf{x})$ perpendicular to the GB, with and without the consideration of the modulus effect, and evidences that the consideration of modulus effect leads to a better description of the segregation energy $E^{\text{Seg}}(\mathbf{x})$ (Eq. 1). Besides, the asymmetry between segregation and anti-segregation zones observed in the spectrum of segregation energies in Fig. 12 can be partly explained by the modulus effect. Indeed, Fig. 15 compares the spectrum of segregation energies with and without the modulus effect for the three methods of computing the interaction energy. Considering isotropic elasticity and size effect only, no asymmetry occurs. A small asymmetry is visible with the size effect alone when considering the computation of GB strain fields in anisotropic elasticity. However, the addition of the modulus effect creates or reinforces considerably the asymmetry of the spectrum and leads to results in better agreement with the spectrum of segregation energies obtained by MS (Fig. 15).

<figure><img src="image_16.png" /><figcaption>Figure 15: Histogram of the fraction of nodes having specific energy with and without consideration of the modulus effect for two GBs: (a) $\theta = 5^\circ$, $\mathbf{b} = \frac{1}{2}[110]$ and (b) $\theta = 10^\circ$, $\mathbf{b} = \frac{1}{2}[110]$ in the Ag-Ni system.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16: Variation with the tilt angle $\theta$ of the overall solute concentration at GB, $X^{GB}$, for two different kinds of LAGBs identified by their edge dislocation Burgers vector. The system considered is Ag/Ni. $X^{GB}$ is computed from Eqs. 37 and 38 with $X^0 = 0.01\%$ (small bulk concentration) and $T = 300K$ and considering fully continuum approaches either in anisotropic elasticity or in isotropic elasticity to evaluate $E^{Seg}(\mathbf{x})$. The GB zone of influence is defined as all the positions with $E^{int}$ values smaller than -0.01 eV or greater than 0.01 eV.</figcaption></figure>
## 6.2. Estimate of solute concentration at GBs

Following the extension of McLean's treatment [10] by White and Stein [13, 12], the solute concentration at a given position $\boldsymbol{x}$ can be expressed as:

$$X(\boldsymbol{x}) = \frac{X^0 \exp\left(-\frac{E^{\text{int}}(\boldsymbol{x})}{k_B T}\right)}{1 - X^0 + X^0 \exp\left(-\frac{E^{\text{int}}(\boldsymbol{x})}{k_B T}\right)}, \quad (37)$$

where $X^0$ is the fraction of solute in the bulk lattice under zero pre-existing stress, $k_B$ is the Boltzmann constant and $T$ is the considered temperature. Then, disregarding the interactions between solutes, the overall solute concentration at a grain boundary, $X^{GB}$, can be estimated from an average of $X(\boldsymbol{x})$ over the GB zone of influence [13, 12]:

$$X^{GB} = \langle X(\boldsymbol{x}) \rangle_{GB}. \quad (38)$$

Hereafter, the GB zone of influence is defined as all the positions with $E^{\text{int}}$ values smaller than -0.01 eV or greater than 0.01 eV. This means that points with values between these two bounds are not considered as being part of the GB zone of influence since they are not likely to be sites of segregation or anti-segregation. The major interest to have a continuum-based model of segregation is the possibility to study the effects of different parameters, such as the temperature [13] or the tilt angle $\theta$. Hence, Fig. 16 shows the variation of $X^{GB}$ with $\theta$ for two different kinds of LAGBs of the Ag/Ni system, considering fully continuum approaches to evaluate $E^{\text{int}}(\boldsymbol{x})$, either in heterogeneous anisotropic elasticity or in isotropic elasticity. As in the spectra energy plots of Section 5.3, the box size along $x_2$ is fixed at $2 \times 47.5$ Å.

The first point to be noticed is that the continuum approach in isotropic elasticity predicts virtually exactly the same overall concentration for both GBs, as opposed to the approach in anisotropic elasticity that exhibits clear distinct behaviors between both GBs (Fig. 16). The deviation of solute concentration between both GBs becomes more and more pronounced as the tilt angle is increased. Also, the discrepancy between $X^{GB}$ predictions in anisotropic and in isotropic elasticity increases with the tilt angle. At $\theta = 15^\circ$, the predicted value of $X^{GB}$ from isotropic elasticity is more than twice the value obtained in anisotropic elasticity for one of the GBs (Fig. 16). More importantly, the impact of the modulus effect on $X^{GB}$ appears to be very significant. This is evidenced in Fig. 16 from comparisons with results considering only the size effect in the evaluation of $E^{\text{int}}(\boldsymbol{x})$ in Eq. 37. As for the effect of anisotropic elasticity, the impact of the modulus effect increases with the tilt angle $\theta$.

## 6.3. Limitations of the proposed approach

Evidently, the present approach in linear elasticity has not for objective to handle non-linear effects at defects cores. In the present study, both the dislocations and the elastic dipoles are singular at their cores because of the use of traditional elastic Green's functions. Moreover, the expression of the elastic interaction energy of Eq. 31 is obtained from a series expansion at the first order of the GB displacement field [25, 54], i.e. under the assumption that the GB displacement field varies slowly close to the solute. Obviously, this assumption becomes less and less valid when approaching the GB dislocations cores, even when considering the GB strain or stress fields obtained by molecular statics. The use of non-singular treatments of dislocations (e.g., [57]) and/or regularized Green's function methods could help to overcome such issues and model the interactions
of solutes with dislocation cores [26]. Despite these aspects, results of Section 5.2 show that the interaction energies with GB strains obtained by MS or computed in anisotropic elasticity match very well the segregation energy for distances as small as 4 Å from the GB dislocations cores. This demonstrates the relevance of the elastic dipole approximation to model the interactions between defects outside the close vicinity of their cores. Besides, it should be underlined that a term is missing in the evaluation of the interaction energy between a solute and a GB in heterogeneous linear elasticity. Indeed, even if the absence of GB strain fields, the solute interacts with the GB because of difference of stiffness tensors on both sides which induce the addition of image terms to the elastic fields of a solute in the infinite medium in order to satisfy the boundary conditions at the interface [51, 58]. It is noteworthy that, if the elastic strain energy does not depend on the crystallographic orientation when modeling a solute as a pure dilatation center, the interaction energy with an interface separating two differently oriented grains of a same material is nevertheless not non-zero because of the presence of the additional image term [58]. However, the latter effect was shown to be negligible in the case of low angle symmetrical tilt boundaries as considered in the present work [58].

Moreover, several studies have reported that elastic dipole and polarizability tensors obtained by interatomic potentials can differ significantly from tensors evaluated from *ab initio* methods [24, 31]. However, the purpose of the present study is not to make quantitative predictive calculations but to show to what extent an approach based on linear elasticity is relevant to model solute segregation at GBs. Elastic dipole and polarizability tensors calculated from *ab initio* methods like the density functional theory (DFT) or inferred from experiments can be readily used in the developed continuum framework to make up-scaling predictions.

## 7. Conclusions

In this study, three symmetric tilt low-angle grain boundaries and two cubic systems (Cu and Ag) with solute atoms in substitution (Ag and Ni, respectively) are considered to investigate the continuum modelling of solute segregation at grain boundaries. The elastic dipole approximation is used to model the interaction energies between solutes and GBs. The strain fields are evaluated either through molecular statics (MS) simulations, or by considering arrays of edge dislocations within the frameworks of linear isotropic elasticity or heterogeneous anisotropic elasticity using the Stroh formalism. These interaction energies are compared to segregation energies computed on a site-by-site basis using MS. The following key conclusions emerge:

* The importance of considering anisotropic elasticity to model solute segregation in highly anisotropic materials like Cu or Ag has been evidenced. The interaction energy model based on isotropic elasticity is unable to well describe the segregation and anti-segregation domains around grain boundaries for both systems studied. As opposed, the same model in heterogeneous anisotropic elasticity leads to much better matches with segregation energy fields computed by MS. Moreover, variations in behavior between grain boundaries having a same tilt angle can only be captured when anisotropic elasticity is considered.

* Despite the inherent limitations in addressing non-linear effects at defect cores, the elastic dipole approximation proves to be an effective method for approximating segregation energy spectra in symmetric tilt low-angle grain boundaries obtained through atomistic simulations. The interaction energies with GB strains obtained by MS or computed from linear anisotropic elasticity match very well the segregation energy for distances as small as 4 Å from the GB dislocations cores.
* Once standardized, the procedures for obtaining the permanent second-rank tensor and the fourth-rank polarizability tensor of the elastic dipole via MS simulations are remarkably efficient. The components of these tensors can serve as valuable tools for comparing various atomic potentials and, ultimately, might be compared with experimental data. Such comparisons may help evaluate the effectiveness of different potentials in modeling solute segregation without performing time-consuming atomistic simulations.

* The polarizability tensor is associated with a modulus effect, which is shown to have a significant impact on estimating the overall solute concentration at grain boundaries. The expression of the interaction energy for cubic lattices given by Eq. 32 reveals that the modulus effect generally promotes either segregation or anti-segregation, as the coefficient $\alpha_K$, $\alpha_{\mu'}$ and $\alpha_{\mu''}$ are multiplied by squared strain components. This observation partially accounts for the asymmetry observed in the segregation energy spectra.

## 8. Data availability

The codes to compute elastic fields of straight dislocation walls in isotropic and anisotropic elasticity are available at: https://github.com/AniPlas/DislocationsWall.

The code to compute the elastic dipole of a single solute atom in substitution is available at: https://github.com/JoePetrazoller/Elastic-dipole.

The codes to compute the polarizability tensor of a single solute atom in substitution are available at: https://github.com/JoePetrazoller/Polarizability.

The code to perform the interpolation of Virial stress fields is available at: https://github.com/JoePetrazoller/Fields-Interpolation.

The dataset of all the grain boundaries generated by the softwares Atomsk and LAMMPS, as well as the associated segregation energy maps computed by molecular statics are available at: https://doi.org/10.57745/ACQWYZ.

## 9. Author contributions: CRediT

JP: Conceptualization, Software, Data curation, Formal analysis, Investigation, Visualization, Writing - original draft

JG: Conceptualization, Supervision, Funding acquisition, Validation, Methodology, Writing - review & editing

SB: Conceptualization, Supervision, Funding acquisition, Validation, Writing - review & editing

TR: Conceptualization, Software, Supervision, Funding acquisition, Validation, Methodology, Investigation, Project administration, Writing - review & editing

## 10. Acknowledgments

The authors are grateful to the French National Research Agency (ANR) for financial support under contracts ANR-22-CE92-0058-01 (SILA) and ANR-21-CE08-0001 (ATOUUM). The authors acknowledge the use of the Cassiopee Arts et Métiers Institute of Technology HPC Center made available for conducting the research reported in this paper.
References

[1] A. Sutton, R. W. Balluffi, Interfaces in Crystalline Materials, corr. ed. Edition, no. 51 in Monographs on the Physics and Chemistry of Materials, Clarendon Press Oxford university press, Oxford Oxford New York, 1996.

[2] P. M. Anderson, P. M. Anderson, J. P. Hirth, J. Lothe, Theory of Dislocations, third edition Edition, Cambridge University Press, New York, NY, USA, 2017.

[3] R. Latanision, H. Opperhauser, The intergranular embrittlement of nickel by hydrogen: the effect of grain boundary segregation, Metallurgical and Materials Transactions B 5 (1974) 483-492.

[4] A. Kumar, B. L. Eyre, Grain boundary segregation and intergranular fracture in molybdenum, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 370 (1743) (1980) 431-458.

[5] G. Schoeck, Moving dislocations and solute atoms, Physical Review 102 (6) (1956) 1458.

[6] N. Ma, S. Dregia, Y. Wang, Solute segregation transition and drag force on grain boundaries, Acta materialia 51 (13) (2003) 3687-3700.

[7] R. Koju, Y. Mishin, Atomistic study of grain-boundary segregation and grain-boundary diffusion in al-mg alloys, Acta Materialia 201 (2020) 596-603.

[8] R. A. Alberty, On the derivation of the Gibbs adsorption equation, Langmuir 11 (9) (1995) 3598-3600.

[9] L. Priester, Joints de grains et plasticité cristalline, Mécanique et ingénierie des matériaux, Hermès science publications-Lavoisier, 2011.

[10] D. McLean, A. Maradudin, Grain boundaries in metals, American Institute of Physics, 1958.

[11] I. Langmuir, The adsorption of gases on plane surfaces of glass, mica and platinum., Journal of the American Chemical Society 40 (9) (1918) 1361-1403.

[12] C. White, D. Stein, Sulfur segregation to grain boundaries in ni 3 al and ni 3 (ai, ti) alloys, Metallurgical Transactions A 9 (1978) 13-22.

[13] C. L. White, W. A. Coghlan, The spectrum of binding energies approach to grain boundary segregation, Metallurgical Transactions A 8 (9) (1977) 1403-1412.

[14] J. D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. R. Soc. Lond. A (241) (1957) 376-396.

[15] W. Cai, R. Sills, D. Barnett, W. Nix, Modeling a distribution of point defects as misfitting inclusions in stressed solids, Journal of the Mechanics and Physics of Solids 66 (2014) 154-171.

[16] R. Dingreville, S. Berbenni, On the interaction of solutes with grain boundaries, Acta Materialia 104 (2016) 237-249.

[17] A. Love, A Treatise on the Mathematical Theory of Elasticity, the University Press, 1959.
[18] E. Kröner, Allgemeine Kontinuumstheorie der Versetzungen und Eigenspannungen, Archive for Rational Mechanics and Analysis 4 (1) (1959) 273–334.

[19] R. Siems, Mechanical Interactions of Point Defects, Physica Status Solidi (b) 30 (2) (1968) 645–658.

[20] G. Leibfried, N. Breuer, Point defects in metals I: introduction to the theory, Springer, 1978.

[21] D. Bacon, D. Barnett, R. O. Scattergood, Anisotropic continuum theory of lattice defects, Progress in Materials Science 23 (1980) 51–262.

[22] A. Vattré, T. Jourdan, H. Ding, M.-C. Marinica, M. Demkowicz, Non-random walk diffusion enhances the sink strength of semicoherent interfaces, Nature communications 7 (1) (2016) 10424.

[23] R. W. Balluffi, Introduction to Elasticity Theory for Crystal Defects, 2nd Edition, World Scientific, 2016.

[24] C. Varvenne, E. Clouet, Elastic dipoles of point defects from atomistic simulations, Physical Review B 96 (22) (2017) 224103.

[25] E. Clouet, C. Varvenne, T. Jourdan, Elastic modeling of point-defects and their interaction, Computational Materials Science 147 (2018) 49–63.

[26] G. Gengor, O. Celebi, A. Mohammed, H. Sehitoglu, Continuum strain of point defects, Journal of the Mechanics and Physics of Solids 188 (2024) 105653.

[27] E. Kröner, Dia-and para-elasticity, Theory of crystal defects édité par Kroupa F. (Academia, 1964) p (1964) 215–230.

[28] K. H. Robrock, W. Schilling, Diaelastic modulus change of aluminium after low temperature electron irradiation, Journal of Physics F: Metal Physics 6 (3) (1976) 303–314.

[29] H. R. Schober, Polarizabilities of point defects in metals, Journal of Nuclear Materials 126 (1984) 220–225.

[30] G. Ackland, Theoretical study of the effect of point defects on the elastic constants of copper, Journal of Nuclear Materials 152 (1) (1988) 53–63.

[31] D. Da Fonseca, F. Onimus, F. Mompiou, M.-C. Marinica, E. De Sonis, E. Clouet, T. Jourdan, Numerical investigation of dislocation climb under stress and irradiation, Acta Materialia 242 (2023) 118431.

[32] J. P. Hirth, J. Lothe, Theory of dislocations, 2nd Edition, Wiley, New York, 1982.

[33] A. Stroh, Dislocations and cracks in anisotropic elasticity, Philosophical magazine 3 (30) (1958) 625–646.

[34] M. Tucker, Plane boundaries and straight dislocations in elastically anisotropic materials: The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics: Vol 19, No 162, Philosophical Magazine 19 (1969) 1141–1159.
[35] D. M. Barnett, J. Lothe, An image force theorem for dislocations in anisotropic bicrystals, Journal of Physics F: Metal Physics 4 (10) (1974) 1618.

[36] Suo, Singularities, interfaces and cracks in dissimilar anisotropic media, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 427 (1873) (1990) 331–358.

[37] T. Ting, D. Barnett, Image force on line dislocations in anisotropic elastic half-spaces with a fixed boundary, International Journal of Solids and Structures 30 (3) (1993) 313–323.

[38] S. T. Choi, Y. Y. Earmme, Elastic study on singularities interacting with interfaces using alternating technique Part I. Anisotropic trimaterial, International Journal of Solids and Structures (2002).

[39] X. Chen, T. Richeton, C. Motz, S. Berbenni, Elastic fields due to dislocations in anisotropic bi-and tri-materials: applications to discrete dislocation pile-ups at grain boundaries, International Journal of Solids and Structures 164 (2019) 141–156.

[40] R. Pei, Z. Xie, S. Yi, S. Korte-Kerzel, J. Guénolé, T. Al-Samman, Atomistic insights into the inhomogeneous nature of solute segregation to grain boundaries in magnesium, Scripta Materialia 230 (2023) 115432.

[41] P. Hirel, Atomsk: A tool for manipulating and converting atomic data files, Computer Physics Communications 197 (2015) 212–219.

[42] D. Hua, Q. Zhou, Y. Shi, S. Li, K. Hua, H. Wang, S. Li, W. Liu, Revealing the deformation mechanisms of (110) symmetric tilt grain boundaries in CoCrNi medium-entropy alloy, International Journal of Plasticity 171 (2023) 103832.

[43] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. In 'T Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, S. J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications 271 (2022) 108171.

[44] J. Guénolé, W. G. Nöhring, A. Vaid, F. Houllé, Z. Xie, A. Prakash, E. Bitzek, Assessment and optimization of the fast inertial relaxation engine (FIRE) for energy minimization in atomistic simulations and its implementation in LAMMPS, Computational Materials Science 175 (2020) 109584, arXiv:1908.02038 [cond-mat-physics:physics].

[45] M. A. Tschopp, S. P. Coleman, D. L. McDowell, Symmetric and asymmetric tilt grain boundary structure and energy in Cu and Al (and transferability to other fcc metals), Integrating Materials and Manufacturing Innovation 4 (1) (2015) 176–189.

[46] D. M. Heyes, Pressure tensor of partial-charge and point-dipole lattices with bulk and surface geometries, Physical Review B 49 (2) (1994) 755–764.

[47] A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO—the Open Visualization Tool, Modelling and Simulation in Materials Science and Engineering 18 (1) (2010) 015012.
[48] P. L. Williams, Y. Mishin, J. C. Hamilton, An embedded-atom potential for the Cu–Ag system, Modelling and Simulation in Materials Science and Engineering 14 (5) (2006) 817–833.

[49] Z. Pan, V. Borovikov, M. I. Mendelev, F. Sansoz, Development of a semi-empirical potential for simulation of Ni solute segregation into grain boundaries in Ag, Modelling and Simulation in Materials Science and Engineering 26 (7) (2018) 075004.

[50] A. Torres-Sánchez, J. M. Vanegas, M. Arroyo, Examining the Mechanical Equilibrium of Microscopic Stresses in Molecular Simulations, Physical Review Letters 114 (25) (2015) 258102.

[51] T. Richeton, X. Chen, S. Berbenni, Misorientation dependence of the grain boundary migration rate: role of elastic anisotropy, Philosophical Magazine 100 (19) (2020) 2436–2457.

[52] R. Hill, The Elastic Behaviour of a Crystalline Aggregate, Proceedings of the Physical Society. Section A 65 (5) (1952) 349–354.

[53] R. deWit, Theory of disclinations: II. Continuous and discrete disclinations in anisotropic elasticity, Journal of Research of the National Bureau of Standards Section A: Physics and Chemistry 77A (1) (1973) 49.

[54] R. Santos-Güemes, G. Álvarez, J. Segurado, An fft based adaptive grid framework to represent non-singular dislocations, Mechanics of Materials 194 (2024) 105004.

[55] L. Walpole, Elastic behavior of composite materials: theoretical foundations, Advances in applied mechanics 21 (1981) 169–242.

[56] R. deWit, Elastic constants and thermal expansion averages of a nontextured polycrystal, Journal of mechanics of materials and structures 3 (2) (2008) 195–212.

[57] W. Cai, A. Arsenlis, C. R. Weinberger, V. V. Bulatov, A non-singular continuum theory of dislocations, Journal of the Mechanics and Physics of Solids 54 (3) (2006) 561–587.

[58] T. Richeton, Energy variation and stress fields of spherical inclusions with eigenstrain in three-dimensional anisotropic bi-materials, International Journal of Solids and Structures 293 (2024) 112748.