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

HAL
open science

# Energy variation and stress fields of spherical inclusions with eigenstrain in three-dimensional anisotropic bi-materials

Thiebaud Richeton

► To cite this version:

Thiebaud Richeton. Energy variation and stress fields of spherical inclusions with eigenstrain in three-dimensional anisotropic bi-materials. 2023. hal-04285488v2

HAL Id: hal-04285488

https://hal.univ-lorraine.fr/hal-04285488v2

Preprint submitted on 12 Feb 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Energy variation and stress fields of spherical inclusions with eigenstrain in three-dimensional anisotropic bi-materials

Thiebaud Richeton

► To cite this version:

Thiebaud Richeton. Energy variation and stress fields of spherical inclusions with eigenstrain in three-dimensional anisotropic bi-materials. 2023. hal-04285488v2

HAL Id: hal-04285488

https://hal.univ-lorraine.fr/hal-04285488v2

Preprint submitted on 12 Feb 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
CC-BY 4.0

# Energy variation and stress fields of spherical inclusions with eigenstrain in three-dimensional anisotropic bi-materials

T. Richeton<sup>a,*</sup>

<sup>a</sup>CNRS, Université de Lorraine, Arts et Métiers, LEM3, F-57000 Metz, France

## Abstract

In this study, expressions for the interior stress fields of a spherical inclusion with uniform eigenstrain embedded in an anisotropic bi-material featuring a planar interface are derived. These expressions involve surface integrals of the imaginary term of the first derivative of the Green tensor in an anisotropic bi-material which are well-suited for standard numerical integrations. Specific formulations are provided for cases where the inclusion belongs to either the same material or both materials. Additionally, expressions are presented for the equivalent Eshelby tensor. The interior stress fields and variations in elastic strain energy are computed using cubic elastic constants of Cu. Various inclusion positions relative to the interface, eigenstrain forms, and crystallographic orientations are considered. For instances involving dilational eigenstrain, the elastic strain energy variation with the inclusion's position may exhibit multiple extrema. The global minimum and maximum consistently occur when the inclusion spans both materials. In the context of symmetrical tilt boundaries, energy variations are perfectly symmetric, with a global minimum on the interface that decreases with the tilt angle. The

<footnote>*Corresponding author.</footnote>

<footnote>Email address: thiebaud.richeton@univ-lorraine.fr (T. Richeton)</footnote>
CC-BY 4.0

# Energy variation and stress fields of spherical inclusions with eigenstrain in three-dimensional anisotropic bi-materials

T. Richeton<sup>a,*</sup>

<sup>a</sup>CNRS, Université de Lorraine, Arts et Métiers, LEM3, F-57000 Metz, France

## Abstract

In this study, expressions for the interior stress fields of a spherical inclusion with uniform eigenstrain embedded in an anisotropic bi-material featuring a planar interface are derived. These expressions involve surface integrals of the imaginary term of the first derivative of the Green tensor in an anisotropic bi-material which are well-suited for standard numerical integrations. Specific formulations are provided for cases where the inclusion belongs to either the same material or both materials. Additionally, expressions are presented for the equivalent Eshelby tensor. The interior stress fields and variations in elastic strain energy are computed using cubic elastic constants of Cu. Various inclusion positions relative to the interface, eigenstrain forms, and crystallographic orientations are considered. For instances involving dilational eigenstrain, the elastic strain energy variation with the inclusion's position may exhibit multiple extrema. The global minimum and maximum consistently occur when the inclusion spans both materials. In the context of symmetrical tilt boundaries, energy variations are perfectly symmetric, with a global minimum on the interface that decreases with the tilt angle. The

<footnote>*Corresponding author.</footnote>

<footnote>Email address: thiebaud.richeton@univ-lorraine.fr (T. Richeton)</footnote>
significance of the observed energy variations for defects segregation is quan-
titatively assessed by comparing them with the interaction energy between
an eigenstrain and a grain boundary stress field.

Keywords: Inclusions, Interface, Anisotropic elasticity, Green function,
Image force

**1. Introduction**

The study of elastic fields generated by inclusions within an elastic ma-
trix stands as a cornerstone in micromechanics, significantly influencing the
development of advanced multi-components alloys (Mura, 1987; Zhou, 2013;
Barnett and Cai, 2018). An inclusion is associated with an eigenstrain, rep-
resenting the stress-free strain the inclusion would have if removed from the
encompassing matrix (Mura, 1987; Barnett and Cai, 2018). Inclusions are
generally considered in an infinite homogeneous medium as in the pioneering
work of Eshelby (1957). However, many works also focus on inclusion prob-
lems within finite spaces (e.g., (Eshelby, 1954, 1961; Kinoshita and Mura,
1984; Li et al., 2007)) and half-spaces (e.g., (Mindlin and Cheng, 1950; Chiu,
1978; Seo and Mura, 1979; Pan and Yang, 2001; Jiang and Pan, 2004)), neces-
sitating the introduction of ‘image’ terms, i.e. terms that should be added to
the elastic fields of the infinite medium in order to satisfy the boundary condi-
tions at outer surfaces. Despite this extensive exploration, limited attention
is given to inclusions in bi-materials, which also demand consideration of
‘image’ terms. Existing studies typically focus on dissimilar isotropic elastic
solids (Bacon, 1972; Gladwell, 1999; Selvadurai, 2000; Franciosi and Charles,
2016) and/or plane or anti-plane problems (Zhang and Chou, 1985; Ru, 2001;
Ru et al., 2001; Wang et al., 2007). In particular, to the best of the author’s
knowledge, the variations of stresses within a 3D inclusion in an anisotropic
bi-material and the variation of elastic strain energy with the inclusion’s
position have never been shown and analyzed. Yet, the 3D problem of an
inclusion embedded in an anisotropic bi-material holds potential significance
significance of the observed energy variations for defects segregation is quan-
titatively assessed by comparing them with the interaction energy between
an eigenstrain and a grain boundary stress field.

Keywords: Inclusions, Interface, Anisotropic elasticity, Green function,
Image force

**1. Introduction**

The study of elastic fields generated by inclusions within an elastic ma-
trix stands as a cornerstone in micromechanics, significantly influencing the
development of advanced multi-components alloys (Mura, 1987; Zhou, 2013;
Barnett and Cai, 2018). An inclusion is associated with an eigenstrain, rep-
resenting the stress-free strain the inclusion would have if removed from the
encompassing matrix (Mura, 1987; Barnett and Cai, 2018). Inclusions are
generally considered in an infinite homogeneous medium as in the pioneering
work of Eshelby (1957). However, many works also focus on inclusion prob-
lems within finite spaces (e.g., (Eshelby, 1954, 1961; Kinoshita and Mura,
1984; Li et al., 2007)) and half-spaces (e.g., (Mindlin and Cheng, 1950; Chiu,
1978; Seo and Mura, 1979; Pan and Yang, 2001; Jiang and Pan, 2004)), neces-
sitating the introduction of ‘image’ terms, i.e. terms that should be added to
the elastic fields of the infinite medium in order to satisfy the boundary condi-
tions at outer surfaces. Despite this extensive exploration, limited attention
is given to inclusions in bi-materials, which also demand consideration of
‘image’ terms. Existing studies typically focus on dissimilar isotropic elastic
solids (Bacon, 1972; Gladwell, 1999; Selvadurai, 2000; Franciosi and Charles,
2016) and/or plane or anti-plane problems (Zhang and Chou, 1985; Ru, 2001;
Ru et al., 2001; Wang et al., 2007). In particular, to the best of the author’s
knowledge, the variations of stresses within a 3D inclusion in an anisotropic
bi-material and the variation of elastic strain energy with the inclusion’s
position have never been shown and analyzed. Yet, the 3D problem of an
inclusion embedded in an anisotropic bi-material holds potential significance
in various applications, such as assessing the strength of laminated materials with reinforcements or better understanding the segregation at interfaces of defects like solute atoms, vacancies, precipitates (oxydes, carbides, etc.) or voids. In particular, solute atom segregation to grain boundaries markedly influences material properties, including grain boundary energy, mobility and cohesion (Dingreville and Berbenni, 2016). Solute atoms can be treated as spherical misfitting inclusions with purely dilatational eigenstrain (White and Coghlan, 1977; Cai et al., 2014; Dingreville and Berbenni, 2016). Point defects like vacancies and solute atoms (in substitution or in insertion) can also be modeled as a distribution of point-forces, i.e. like a force multipole that mimics the forces imposed on the atoms surrounding the defects (Siems, 1968; Bacon et al., 1980; Balluffi, 2016; Clouet et al., 2018). While this representation is actually only valid outside the defect, it becomes particularly relevant when considering the interaction energy with an external strain field. For instance, it is demonstrated that representing the defect with the first moment of the point force distribution, known as the elastic dipole, is equivalent to considering a uniform eigenstrain within a small inclusion but only in regions far from the inclusion (Clouet et al., 2018). The elastic dipole tensor can also be related to the displacement dipole tensor, which quantifies the strength of the point defect (Lazar, 2017). In the latter approach, the point defect is modeled with a three-dimensional Dirac $\delta$-singularity in the eigenstrain tensor such that only average stress values can be defined at the defect position due to the Dirac core (Lazar, 2017). Consequently, the representation of a defect with a dipole tensor cannot capture the full complexity of the interaction energy between a defect and an interface, as the interior stress fields of an inclusion within a bi-material are necessarily non-uniform (see, e.g., (Seo and Mura, 1979; Wang et al., 2007)), whereas a finite inclusion may also spans both materials.

The current study deals with the issue of an inclusion of finite size embedded in a bi-material considering linear heterogeneous anisotropic elastic-
in various applications, such as assessing the strength of laminated materials with reinforcements or better understanding the segregation at interfaces of defects like solute atoms, vacancies, precipitates (oxydes, carbides, etc.) or voids. In particular, solute atom segregation to grain boundaries markedly influences material properties, including grain boundary energy, mobility and cohesion (Dingreville and Berbenni, 2016). Solute atoms can be treated as spherical misfitting inclusions with purely dilatational eigenstrain (White and Coghlan, 1977; Cai et al., 2014; Dingreville and Berbenni, 2016). Point defects like vacancies and solute atoms (in substitution or in insertion) can also be modeled as a distribution of point-forces, i.e. like a force multipole that mimics the forces imposed on the atoms surrounding the defects (Siems, 1968; Bacon et al., 1980; Balluffi, 2016; Clouet et al., 2018). While this representation is actually only valid outside the defect, it becomes particularly relevant when considering the interaction energy with an external strain field. For instance, it is demonstrated that representing the defect with the first moment of the point force distribution, known as the elastic dipole, is equivalent to considering a uniform eigenstrain within a small inclusion but only in regions far from the inclusion (Clouet et al., 2018). The elastic dipole tensor can also be related to the displacement dipole tensor, which quantifies the strength of the point defect (Lazar, 2017). In the latter approach, the point defect is modeled with a three-dimensional Dirac $\delta$-singularity in the eigenstrain tensor such that only average stress values can be defined at the defect position due to the Dirac core (Lazar, 2017). Consequently, the representation of a defect with a dipole tensor cannot capture the full complexity of the interaction energy between a defect and an interface, as the interior stress fields of an inclusion within a bi-material are necessarily non-uniform (see, e.g., (Seo and Mura, 1979; Wang et al., 2007)), whereas a finite inclusion may also spans both materials.

The current study deals with the issue of an inclusion of finite size embedded in a bi-material considering linear heterogeneous anisotropic elastic-
ity (see Fig.1). As explained, the interface introduces modifications to the
elastic fields within and surrounding the inclusion compared to the infinite
medium. The gradient of the change in elastic strain energy associated with
the defect manifests as a force on the defect itself, known as an image force.
This configurational force, exhibiting either an attractive or repulsive na-
ture, depends on the defect’s position relative to the interface. It is known
that a softer material will attract the defect (whether it is a dislocation or
a point defect), while a harder one will repel it. This is rooted in the fact
that the elastic distortions induced by the defect yield lower contributions to
the elastic strain energy when situated within the softer material. However,
the scenario becomes more intricate when considering a grain boundary that
separates two crystals of the same material but with different orientations.
In addition, the details of how the elastic strain energy varies as the inclusion
is continuously displaced from one material to another remain unexplored.

The stress fields and elastic strain energy variations for an inclusion within
a bi-material can be effectively computed using the anisotropic elastic Green
functions in a bi-material with a planar interface, as developed by Pan and
Yuan (2000), and combining them with the Green's function formulation for
the inclusion stress fields (Mura, 1987). The expression of the anisotropic
elastic Green functions in a bi-material were derived via the Stroh formal-
ism (Stroh, 1958, 1962) and the 2D Fourier transform method. They encom-
pass the full-space Green's tensor and an additional complementary part.
The latter can be efficiently computed through a regular line-integral over
the range [0, π], amenable to standard numerical integration methods (Pan
and Yuan, 2000; Chu et al., 2011, 2012).

The paper is organized as follows. Section 2 presents the theoretical
framework. The expressions of the stress field due to a distribution of
eigenstrains in a general (i.e. homogeneous or heterogeneous) infinite elastic
medium are first introduced. Then, particular expressions of the stress fields
are derived considering a uniform eigenstrain within a spherical inclusion in
ity (see Fig.1). As explained, the interface introduces modifications to the
elastic fields within and surrounding the inclusion compared to the infinite
medium. The gradient of the change in elastic strain energy associated with
the defect manifests as a force on the defect itself, known as an image force.
This configurational force, exhibiting either an attractive or repulsive na-
ture, depends on the defect’s position relative to the interface. It is known
that a softer material will attract the defect (whether it is a dislocation or
a point defect), while a harder one will repel it. This is rooted in the fact
that the elastic distortions induced by the defect yield lower contributions to
the elastic strain energy when situated within the softer material. However,
the scenario becomes more intricate when considering a grain boundary that
separates two crystals of the same material but with different orientations.
In addition, the details of how the elastic strain energy varies as the inclusion
is continuously displaced from one material to another remain unexplored.

The stress fields and elastic strain energy variations for an inclusion within
a bi-material can be effectively computed using the anisotropic elastic Green
functions in a bi-material with a planar interface, as developed by Pan and
Yuan (2000), and combining them with the Green's function formulation for
the inclusion stress fields (Mura, 1987). The expression of the anisotropic
elastic Green functions in a bi-material were derived via the Stroh formal-
ism (Stroh, 1958, 1962) and the 2D Fourier transform method. They encom-
pass the full-space Green's tensor and an additional complementary part.
The latter can be efficiently computed through a regular line-integral over
the range [0, π], amenable to standard numerical integration methods (Pan
and Yuan, 2000; Chu et al., 2011, 2012).

The paper is organized as follows. Section 2 presents the theoretical
framework. The expressions of the stress field due to a distribution of
eigenstrains in a general (i.e. homogeneous or heterogeneous) infinite elastic
medium are first introduced. Then, particular expressions of the stress fields
are derived considering a uniform eigenstrain within a spherical inclusion in
a bi-material, whether it belongs to a single material or spans both mate-
rials (see Fig.1). The latter configuration might provide insights into the
impact of elastic heterogeneity when a large defect, such as a precipitate,
is within a grain boundary, or when a small defect, such as a solute atom,
is in the transition area between the bulk crystal and the grain boundary.
Furthermore, the expression for the elastic strain energy is provided, along
with simplifications emerging from the assumption of cubic elasticity and
dilatational eigenstrain. In Section 3, a convergence study is undertaken re-
garding the numerical integrations involved in the computation of the elastic
strain energy. Additionally, a comparative analysis with the known stresses
and energy variations in an isotropic half-space is performed, ensuring the
robustness of the numerical implementation. Moving on to Section 4, the
results in terms of interior stress fields profiles and variations in elastic strain
energy are shown assuming cubic elastic constants of Cu. Various inclusion
positions relative to the interface, diverse forms of eigenstrain, and distinct
crystallographic orientations are explored. Section 5 engages in a thorough
discussion of these results, with a specific focus on quantitatively assessing
the significance of the image effect. This assessment is made through com-
parisons with the interaction energy between an eigenstrain and an external
stress field, such as the one arising from a grain boundary. Summarizing
remarks follow in Section 6.

Throughout the paper, bold letters are used for vectors, tensors and ma-
trices. The superscript T denotes the transpose of a matrix and an overbar
the conjugate of a complex function. The Einstein convention over repeated
indices is used. i denotes the imaginary unit.

2. Theoretical framework

We consider an infinite general (i.e. homogeneous or heterogeneous) elastic medium of volume $V$ and a distribution of eigenstrains $\epsilon^*$. In the absence of body force, the displacement vector $\mathbf{u}$ follows the Navier equation
a bi-material, whether it belongs to a single material or spans both mate-
rials (see Fig.1). The latter configuration might provide insights into the
impact of elastic heterogeneity when a large defect, such as a precipitate,
is within a grain boundary, or when a small defect, such as a solute atom,
is in the transition area between the bulk crystal and the grain boundary.
Furthermore, the expression for the elastic strain energy is provided, along
with simplifications emerging from the assumption of cubic elasticity and
dilatational eigenstrain. In Section 3, a convergence study is undertaken re-
garding the numerical integrations involved in the computation of the elastic
strain energy. Additionally, a comparative analysis with the known stresses
and energy variations in an isotropic half-space is performed, ensuring the
robustness of the numerical implementation. Moving on to Section 4, the
results in terms of interior stress fields profiles and variations in elastic strain
energy are shown assuming cubic elastic constants of Cu. Various inclusion
positions relative to the interface, diverse forms of eigenstrain, and distinct
crystallographic orientations are explored. Section 5 engages in a thorough
discussion of these results, with a specific focus on quantitatively assessing
the significance of the image effect. This assessment is made through com-
parisons with the interaction energy between an eigenstrain and an external
stress field, such as the one arising from a grain boundary. Summarizing
remarks follow in Section 6.

Throughout the paper, bold letters are used for vectors, tensors and ma-
trices. The superscript T denotes the transpose of a matrix and an overbar
the conjugate of a complex function. The Einstein convention over repeated
indices is used. i denotes the imaginary unit.

2. Theoretical framework

We consider an infinite general (i.e. homogeneous or heterogeneous) elas-
tic medium of volume V and a distribution of eigenstrains ε*. In the absence
of body force, the displacement vector **u** follows the Navier equation
<figure><img src="image_3.png" /><figcaption>Figure 1: Infinite bi-material made of material (1) characterized by elasticity tensor $C^{(1)}$ in the region $x_3 > 0$ and of material (2) characterized by $C^{(2)}$ in the region $x_3 < 0$. On the left, the spherical inclusion with a uniform eigenstrain $\epsilon^*$ belongs only to material (2). On the right, the inclusion spans the two materials.</figcaption></figure>

$$C_{ijkl}(\mathbf{x}) \frac{\partial^2 u_k(\mathbf{x})}{\partial x_l \partial x_j} = C_{ijkl}(\mathbf{x}) \frac{\partial \epsilon_{kl}^*(\mathbf{x})}{\partial x_j} \quad (1)$$

where $C$ is the fourth-rank tensor of anisotropic elastic constants that possesses the symmetries $C_{ijkl} = C_{klij} = C_{jikl} = C_{ijlk}$. $\mathbf{x}$ denotes a particular point of $V$ and the indices refer to a rectangular Cartesian frame $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3)$. The elastic Green tensor $G$ of the Navier equation is defined by (Barnett, 1972; Mura, 1987; Lazar, 2016, 2017)

$$C_{ijkl}(\mathbf{x}) \frac{\partial^2 G_{km}}{\partial x_l \partial x_j}(\mathbf{x}, \mathbf{y}) + \delta_{im} \delta(\mathbf{x} - \mathbf{y}) = 0 \quad (2)$$

$\delta_{im}$ is the Kronecker delta tensor and $\delta(\mathbf{x} - \mathbf{y})$ is the three-dimensional Dirac delta function. The Green's function $G_{ij}(\mathbf{x}, \mathbf{y})$ represents the elastic displacement in the $\mathbf{e}_i$ direction at point $\mathbf{x}$ induced by a unit force in the $\mathbf{e}_j$ direction applied at point $\mathbf{y}$ (Mura, 1987; Barnett, 1972). For a distribution of eigenstrains, it satisfies the equation (Mura, 1987)

$$u_i(\mathbf{x}) = - \int_V G_{ij}(\mathbf{x}, \mathbf{y}) C_{jlmn}(\mathbf{y}) \frac{\partial \epsilon_{mn}^*(\mathbf{y})}{\partial y_l} d\mathbf{y} \quad (3)$$

where the term $-C_{jlmn}(\mathbf{y}) \frac{\partial \epsilon_{mn}^*(\mathbf{y})}{\partial y_l}$ represents a fictitious body force in
<figure><img src="image_2.png" /><figcaption>Figure 1: Infinite bi-material made of material (1) characterized by elasticity tensor $C^{(1)}$ in the region $x_3 > 0$ and of material (2) characterized by $C^{(2)}$ in the region $x_3 < 0$. On the left, the spherical inclusion with a uniform eigenstrain $\epsilon^*$ belongs only to material (2). On the right, the inclusion spans the two materials.</figcaption></figure>

$$C_{ijkl}(\mathbf{x}) \frac{\partial^2 u_k(\mathbf{x})}{\partial x_l \partial x_j} = C_{ijkl}(\mathbf{x}) \frac{\partial \epsilon_{kl}^*(\mathbf{x})}{\partial x_j} \quad (1)$$

where $C$ is the fourth-rank tensor of anisotropic elastic constants that possesses the symmetries $C_{ijkl} = C_{klij} = C_{jikl} = C_{ijlk}$. $\mathbf{x}$ denotes a particular point of $V$ and the indices refer to a rectangular Cartesian frame $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3)$. The elastic Green tensor $G$ of the Navier equation is defined by (Barnett, 1972; Mura, 1987; Lazar, 2016, 2017)

$$C_{ijkl}(\mathbf{x}) \frac{\partial^2 G_{km}}{\partial x_l \partial x_j}(\mathbf{x}, \mathbf{y}) + \delta_{im} \delta(\mathbf{x} - \mathbf{y}) = 0 \quad (2)$$

$\delta_{im}$ is the Kronecker delta tensor and $\delta(\mathbf{x} - \mathbf{y})$ is the three-dimensional Dirac delta function. The Green's function $G_{ij}(\mathbf{x}, \mathbf{y})$ represents the elastic displacement in the $\mathbf{e}_i$ direction at point $\mathbf{x}$ induced by a unit force in the $\mathbf{e}_j$ direction applied at point $\mathbf{y}$ (Mura, 1987; Barnett, 1972). For a distribution of eigenstrains, it satisfies the equation (Mura, 1987)

$$u_i(\mathbf{x}) = - \int_V G_{ij}(\mathbf{x}, \mathbf{y}) C_{jlmn}(\mathbf{y}) \frac{\partial \epsilon_{mn}^*(\mathbf{y})}{\partial y_l} dy \quad (3)$$

where the term $-C_{jlmn}(\mathbf{y}) \frac{\partial \epsilon_{mn}^*(\mathbf{y})}{\partial y_l}$ represents a fictitious body force in
the $e_j$ direction. By integration by parts and assuming the boundary terms vanish, we have (Mura, 1987)

$$u_i(\mathbf{x}) = \int_V C_{jlmn} \varepsilon_{mn}^*(\mathbf{y}) \frac{\partial G_{ij}(\mathbf{x}, \mathbf{y})}{\partial y_l} d\mathbf{y} \quad (4)$$

Hence, from Hooke's law, the stress field due to a distribution of eigen-strains in a general infinite elastic medium may be expressed as

$$\sigma_{pq}(\mathbf{x}) = C_{pqim}(\mathbf{x}) \left( \int_V C_{jkln}(\mathbf{y}) \varepsilon_{ln}^*(\mathbf{y}) \frac{\partial^2 G_{ij}}{\partial y_k \partial x_m}(\mathbf{x}, \mathbf{y}) d\mathbf{y} - \varepsilon_{im}^*(\mathbf{x}) \right) \quad (5)$$

## 2.1. Homogeneous medium

For an infinite homogeneous body, we introduce $G^\infty$ as the full-space anisotropic elastic Green tensor. We have then $\frac{\partial G_{ij}^\infty}{\partial y_l}(\mathbf{x}, \mathbf{y}) = -\frac{\partial G_{ij}^\infty}{\partial x_l}(\mathbf{x}, \mathbf{y})$ (Mura, 1987), and thus the more classical expression of the stress field for a distribution of eigenstrains is retrieved

$$\sigma_{pq}(\mathbf{x}) = -C_{pqim} \left( \int_V C_{jkln} \varepsilon_{ln}^*(\mathbf{y}) \frac{\partial^2 G_{ij}^\infty}{\partial x_k \partial x_m}(\mathbf{x}, \mathbf{y}) d\mathbf{y} + \varepsilon_{im}^*(\mathbf{x}) \right) \quad (6)$$

## 2.2. Bi-material

We now consider an infinite bi-material made of material (1) characterized by elasticity tensor $C^{(1)}$ in the region $x_3 > 0$ and of material (2) characterized by $C^{(2)}$ in the region $x_3 < 0$. The two materials are thus separated by a planar interface whose normal is along $e_3$. The latter is assumed perfectly bonded. From Eq. 5, the stress field due to a distribution of eigenstrains can be written as

$$\sigma_{pq}(\mathbf{x}) = C_{pqim}(\mathbf{x}) \left( \int_V C_{jkln}(\mathbf{y}) \varepsilon_{ln}^*(\mathbf{y}) \frac{\partial^2 G_{ij}^{Bi}}{\partial y_k \partial x_m}(\mathbf{x}, \mathbf{y}) d\mathbf{y} - \varepsilon_{im}^*(\mathbf{x}) \right) \quad (7)$$
the $e_j$ direction. By integration by parts and assuming the boundary terms vanish, we have (Mura, 1987)

$$u_i(\mathbf{x}) = \int_V C_{jlmn} \varepsilon_{mn}^*(\mathbf{y}) \frac{\partial G_{ij}(\mathbf{x}, \mathbf{y})}{\partial y_l} d\mathbf{y} \quad (4)$$

Hence, from Hooke's law, the stress field due to a distribution of eigen-strains in a general infinite elastic medium may be expressed as

$$\sigma_{pq}(\mathbf{x}) = C_{pqim}(\mathbf{x}) \left( \int_V C_{jkln}(\mathbf{y}) \varepsilon_{ln}^*(\mathbf{y}) \frac{\partial^2 G_{ij}}{\partial y_k \partial x_m}(\mathbf{x}, \mathbf{y}) d\mathbf{y} - \varepsilon_{im}^*(\mathbf{x}) \right) \quad (5)$$

## 2.1. Homogeneous medium

For an infinite homogeneous body, we introduce $G^\infty$ as the full-space anisotropic elastic Green tensor. We have then $\frac{\partial G_{ij}^\infty}{\partial y_l}(\mathbf{x}, \mathbf{y}) = -\frac{\partial G_{ij}^\infty}{\partial x_l}(\mathbf{x}, \mathbf{y})$ (Mura, 1987), and thus the more classical expression of the stress field for a distribution of eigenstrains is retrieved

$$\sigma_{pq}(\mathbf{x}) = -C_{pqim} \left( \int_V C_{jkln} \varepsilon_{ln}^*(\mathbf{y}) \frac{\partial^2 G_{ij}^\infty}{\partial x_k \partial x_m}(\mathbf{x}, \mathbf{y}) d\mathbf{y} + \varepsilon_{im}^*(\mathbf{x}) \right) \quad (6)$$

## 2.2. Bi-material

We now consider an infinite bi-material made of material (1) characterized by elasticity tensor $C^{(1)}$ in the region $x_3 > 0$ and of material (2) characterized by $C^{(2)}$ in the region $x_3 < 0$. The two materials are thus separated by a planar interface whose normal is along $e_3$. The latter is assumed perfectly bonded. From Eq. 5, the stress field due to a distribution of eigenstrains can be written as

$$\sigma_{pq}(\mathbf{x}) = C_{pqim}(\mathbf{x}) \left( \int_V C_{jkln}(\mathbf{y}) \varepsilon_{ln}^*(\mathbf{y}) \frac{\partial^2 G_{ij}^{Bi}}{\partial y_k \partial x_m}(\mathbf{x}, \mathbf{y}) d\mathbf{y} - \varepsilon_{im}^*(\mathbf{x}) \right) \quad (7)$$
where $G^{Bi}$ is the 3D anisotropic elastic Green tensor which was derived by Pan and Yuan (2000). For such a bi-material, it is noteworthy that the derivative relationship is lost, $\frac{\partial G_{ij}^{Bi}}{\partial y_l}(\mathbf{x}, \mathbf{y}) \neq -\frac{\partial G_{ij}^{Bi}}{\partial x_l}(\mathbf{x}, \mathbf{y})$, as well as the symmetry relations, $G_{ij}^{Bi}(\mathbf{x}, \mathbf{y}) \neq G_{ji}^{Bi}(\mathbf{x}, \mathbf{y})$ (Chu et al., 2012).

In the present work, we write $G^{Bi}$ as the sum of a full-space term and an imaginary term $G^{Im}$ as follows

$$G_{ij}^{Bi}(\mathbf{x}, \mathbf{y}) = \begin{cases} G_{ij}^{\infty(1)}(\mathbf{x}, \mathbf{y}) + G_{ij}^{Im}(\mathbf{x}, \mathbf{y}) & \text{if } x_3 > 0 \\ G_{ij}^{\infty(2)}(\mathbf{x}, \mathbf{y}) + G_{ij}^{Im}(\mathbf{x}, \mathbf{y}) & \text{if } x_3 < 0 \end{cases} \quad (8)$$

The expressions of $G^{Im}$ depend on the signs of $x_3$ and $y_3$ (Pan and Yuan, 2000; Chu et al., 2011, 2012). They are explicitly given in the Appendix B.

### 2.2.1. Uniform eigenstrain within a spherical inclusion belonging to a same material

In the following, the problem is further restricted to the issue of a uniform eigenstrain $\varepsilon^*$ within a spherical inclusion of radius $R$ and volume $V^I$. If the whole inclusion completely belongs to a same material, then the stress field inside the spherical inclusion writes

$$\sigma_{pq}(\mathbf{x}) = C_{pqim}(\mathbf{x}) \left( C_{jkln}(\mathbf{x}) \varepsilon_{ln}^* \left( E_{ijkm}(\mathbf{x}) + \int_{V^I} \frac{\partial^2 G_{ij}^{Im}}{\partial y_k \partial x_m}(\mathbf{x}, \mathbf{y}) dy \right) - \varepsilon_{im}^* \right) \quad (9)$$

where the fourth-rank tensor $E$ is defined as (Kröner, 1990; Lazar, 2016)

$$E_{ijkm}(\mathbf{x}) = \begin{cases} E_{ijkm}^{(1)} = - \int_{V^I} \frac{\partial^2 G_{ij}^{\infty(1)}}{\partial x_k \partial x_m}(\mathbf{x}, \mathbf{y}) dy & \text{if } x_3 > 0 \\ E_{ijkm}^{(2)} = - \int_{V^I} \frac{\partial^2 G_{ij}^{\infty(2)}}{\partial x_k \partial x_m}(\mathbf{x}, \mathbf{y}) dy & \text{if } x_3 < 0 \end{cases} \quad (10)$$

$E$ is related to the Dirac $\delta(\mathbf{x}-\mathbf{y})$-term of the second derivative of the full-
where $G^{Bi}$ is the 3D anisotropic elastic Green tensor which was derived by Pan and Yuan (2000). For such a bi-material, it is noteworthy that the derivative relationship is lost, $\frac{\partial G_{ij}^{Bi}}{\partial y_l}(\mathbf{x}, \mathbf{y}) \neq -\frac{\partial G_{ij}^{Bi}}{\partial x_l}(\mathbf{x}, \mathbf{y})$, as well as the symmetry relations, $G_{ij}^{Bi}(\mathbf{x}, \mathbf{y}) \neq G_{ji}^{Bi}(\mathbf{x}, \mathbf{y})$ (Chu et al., 2012).

In the present work, we write $G^{Bi}$ as the sum of a full-space term and an imaginary term $G^{Im}$ as follows

$$G_{ij}^{Bi}(\mathbf{x}, \mathbf{y}) = \begin{cases} G_{ij}^{\infty(1)}(\mathbf{x}, \mathbf{y}) + G_{ij}^{Im}(\mathbf{x}, \mathbf{y}) & \text{if } x_3 > 0 \\ G_{ij}^{\infty(2)}(\mathbf{x}, \mathbf{y}) + G_{ij}^{Im}(\mathbf{x}, \mathbf{y}) & \text{if } x_3 < 0 \end{cases} \quad (8)$$

The expressions of $G^{Im}$ depend on the signs of $x_3$ and $y_3$ (Pan and Yuan, 2000; Chu et al., 2011, 2012). They are explicitly given in the Appendix B.

### 2.2.1. Uniform eigenstrain within a spherical inclusion belonging to a same material

In the following, the problem is further restricted to the issue of a uniform eigenstrain $\varepsilon^*$ within a spherical inclusion of radius $R$ and volume $V^I$. If the whole inclusion completely belongs to a same material, then the stress field inside the spherical inclusion writes

$$\sigma_{pq}(\mathbf{x}) = C_{pqim}(\mathbf{x}) \left( C_{jkln}(\mathbf{x}) \varepsilon_{ln}^* \left( E_{ijkm}(\mathbf{x}) + \int_{V^I} \frac{\partial^2 G_{ij}^{Im}}{\partial y_k \partial x_m}(\mathbf{x}, \mathbf{y}) dy \right) - \varepsilon_{im}^* \right) \quad (9)$$

where the fourth-rank tensor $E$ is defined as (Kröner, 1990; Lazar, 2016)

$$E_{ijkm}(\mathbf{x}) = \begin{cases} E_{ijkm}^{(1)} = -\int_{V^I} \frac{\partial^2 G_{ij}^{\infty(1)}}{\partial x_k \partial x_m}(\mathbf{x}, \mathbf{y}) dy & \text{if } x_3 > 0 \\ E_{ijkm}^{(2)} = -\int_{V^I} \frac{\partial^2 G_{ij}^{\infty(2)}}{\partial x_k \partial x_m}(\mathbf{x}, \mathbf{y}) dy & \text{if } x_3 < 0 \end{cases} \quad (10)$$

$E$ is related to the Dirac $\delta(\mathbf{x}-\mathbf{y})$-term of the second derivative of the full-
space Green tensor (Lazar, 2016). It can be computed almost instantaneously via an integral over the unit sphere (Kneer, 1965; Lazar, 2016) (see details in the Appendix A).

Using Gauss's theorem, the volume integral term is reduced to a surface integral

$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}(\mathbf{x}) \left( C_{jkln}(\mathbf{x}) \varepsilon_{ln}^* \left( E_{ijkm}(\mathbf{x}) + \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} \right) - \varepsilon_{im}^* \right) \quad (11) $$

$\partial V^I$ denotes the surface of the spherical inclusion and $\mathbf{n}$ its external unit normal. It must be underlined that the introduction of the tensor $\mathbf{E}$ and its computation from Eq. A.4 is essential from a numerical point of view. Otherwise, the numerical satisfaction with a good accuracy of Eq. 10 is very slow to achieve.

From Eq. 11, one may also define the Eshelby tensor for a spherical inclusion in a bi-material, $S^{Esh,Bi}(\mathbf{x})$, although being non-uniform. The latter should connect the total strain with the eigenstrain, $\varepsilon_{im}(\mathbf{x}) = S_{imln}^{Esh,Bi}(\mathbf{x})\varepsilon_{ln}^*$ and should possess the symmetries $S_{imln}^{Esh,Bi} = S_{mln}^{Esh,Bi} = S_{imnl}^{Esh,Bi}$. Hence, we might define

$$ S_{imln}^{Esh,Bi}(\mathbf{x}) = \frac{1}{2} C_{jkln}(\mathbf{x}) (E_{ijkm}(\mathbf{x}) + E_{mjk i}(\mathbf{x})) + \frac{1}{2} C_{jkln}(\mathbf{x}) \left( \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} + \int_{\partial V^I} \frac{\partial G_{mj}^{Im}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} \right) \quad (12) $$

The second term is zero for a homogenous material as it is only related to the imaginary parts of the Green tensor. Hence, the first term corresponds indeed to the expression of the Eshelby tensor for a spherical inclusion in an anisotropic homogeneous medium as can be found in (Lazar, 2016). In homogenization theory, the Eshelby tensor in a homogeneous
space Green tensor (Lazar, 2016). It can be computed almost instantaneously via an integral over the unit sphere (Kneer, 1965; Lazar, 2016) (see details in the Appendix A).

Using Gauss's theorem, the volume integral term is reduced to a surface integral

$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}(\mathbf{x}) \left( C_{jkln}(\mathbf{x})\varepsilon_{ln}^* \left( E_{ijkm}(\mathbf{x}) + \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y})n_k d\mathbf{y} \right) - \varepsilon_{im}^* \right) \quad (11) $$

$\partial V^I$ denotes the surface of the spherical inclusion and $\mathbf{n}$ its external unit normal. It must be underlined that the introduction of the tensor $\mathbf{E}$ and its computation from Eq. A.4 is essential from a numerical point of view. Otherwise, the numerical satisfaction with a good accuracy of Eq. 10 is very slow to achieve.

From Eq. 11, one may also define the Eshelby tensor for a spherical inclusion in a bi-material, $S^{Esh,Bi}(\mathbf{x})$, although being non-uniform. The latter should connect the total strain with the eigenstrain, $\varepsilon_{im}(\mathbf{x}) = S_{imln}^{Esh,Bi}(\mathbf{x})\varepsilon_{ln}^*$ and should possess the symmetries $S_{imln}^{Esh,Bi} = S_{miln}^{Esh,Bi} = S_{imnl}^{Esh,Bi}$. Hence, we might define

$$ S_{imln}^{Esh,Bi}(\mathbf{x}) = \frac{1}{2}C_{jkln}(\mathbf{x}) (E_{ijkm}(\mathbf{x}) + E_{mjk i}(\mathbf{x})) + \frac{1}{2}C_{jkln}(\mathbf{x}) \left( \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y})n_k d\mathbf{y} + \int_{\partial V^I} \frac{\partial G_{mj}^{Im}}{\partial x_i}(\mathbf{x}, \mathbf{y})n_k d\mathbf{y} \right) \quad (12) $$

The second term is zero for a homogenous material as it is only related to the imaginary parts of the Green tensor. Hence, the first term corresponds indeed to the expression of the Eshelby tensor for a spherical inclusion in an anisotropic homogeneous medium as can be found in (Lazar, 2016). In homogenization theory, the Eshelby tensor in a homogeneous
medium is more classically given as $S_{imln}^{Esh,\infty} = P_{imjk}^{\infty}C_{jkln}$ where $P_{imjk}^{\infty} = \frac{1}{2} \int_{V^I} (\Gamma_{imjk}^{\infty}(\mathbf{x}, \mathbf{y}) + \Gamma_{imkj}^{\infty}(\mathbf{x}, \mathbf{y})) \, dy$ with $\mathbf{P}^{\infty}$ the Hill's polarization tensor (Hill, 1965) and $\Gamma^{\infty}$ the modified Green function (Kröner, 1990). For a spherical inclusion, we have $P_{imjk}^{\infty} = \frac{1}{4}(E_{ijkm} + E_{mjki} + E_{ikjm} + E_{mkji})$. For the analogy, it is therefore interesting to re-write the Eshelby tensor for a spherical inclusion in a bi-material as

$$S_{imln}^{Esh,Bi}(\mathbf{x}) = P_{imjk}^{Bi}(\mathbf{x})C_{jkln}(\mathbf{x}) \qquad (13)$$

with

$$P_{imjk}^{Bi}(\mathbf{x}) = P_{imjk}^{\infty}(\mathbf{x}) + P_{imjk}^{Im}(\mathbf{x}) \qquad (14)$$

and

$$P_{imjk}^{Im}(\mathbf{x}) = \frac{1}{4} \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k dy + \frac{1}{4} \int_{\partial V^I} \frac{\partial G_{mj}^{Im}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k dy \\ + \frac{1}{4} \int_{\partial V^I} \frac{\partial G_{ik}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_j dy + \frac{1}{4} \int_{\partial V^I} \frac{\partial G_{mk}^{Im}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_j dy \qquad (15)$$

It is noteworthy that $\mathbf{P}^{\infty}(\mathbf{x})$ is uniform within a same material whereas $\mathbf{P}^{Im}(\mathbf{x})$ is not.

### 2.2.2. Uniform eigenstrain within a spherical inclusion spanning the two materials

If the inclusion spans the two materials, i.e. the center of the spherical inclusion is such that $y_3 - R \le 0$ and $y_3 + R \ge 0$, the elasticity tensor cannot be taken out of the integral. However, one can write $C(\mathbf{y}) = C^{(2)} + (C^{(1)} - C^{(2)}) H(y_3)$ where $H(y_3)$ is the Heaviside step function. For $x_3 < 0$, the following equation is thus obtained for the stress field inside the inclusion
medium is more classically given as $S_{imln}^{Esh,\infty} = P_{imjk}^{\infty}C_{jkln}$ where $P_{imjk}^{\infty} = \frac{1}{2} \int_{V^I} (\Gamma_{imjk}^{\infty}(\mathbf{x}, \mathbf{y}) + \Gamma_{imkj}^{\infty}(\mathbf{x}, \mathbf{y})) \, d\mathbf{y}$ with $\mathbf{P}^{\infty}$ the Hill's polarization tensor (Hill, 1965) and $\Gamma^{\infty}$ the modified Green function (Kröner, 1990). For a spherical inclusion, we have $P_{imjk}^{\infty} = \frac{1}{4}(E_{ijkm} + E_{mjki} + E_{ikjm} + E_{mkji})$. For the analogy, it is therefore interesting to re-write the Eshelby tensor for a spherical inclusion in a bi-material as

$$S_{imln}^{Esh,Bi}(\mathbf{x}) = P_{imjk}^{Bi}(\mathbf{x})C_{jkln}(\mathbf{x}) \qquad (13)$$

with

$$P_{imjk}^{Bi}(\mathbf{x}) = P_{imjk}^{\infty}(\mathbf{x}) + P_{imjk}^{Im}(\mathbf{x}) \qquad (14)$$

and

$$P_{imjk}^{Im}(\mathbf{x}) = \frac{1}{4} \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} + \frac{1}{4} \int_{\partial V^I} \frac{\partial G_{mj}^{Im}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} \\ + \frac{1}{4} \int_{\partial V^I} \frac{\partial G_{ik}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_j d\mathbf{y} + \frac{1}{4} \int_{\partial V^I} \frac{\partial G_{mk}^{Im}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_j d\mathbf{y} \qquad (15)$$

It is noteworthy that $\mathbf{P}^{\infty}(\mathbf{x})$ is uniform within a same material whereas $\mathbf{P}^{Im}(\mathbf{x})$ is not.

### 2.2.2. Uniform eigenstrain within a spherical inclusion spanning the two materials

If the inclusion spans the two materials, i.e. the center of the spherical inclusion is such that $y_3 - R \le 0$ and $y_3 + R \ge 0$, the elasticity tensor cannot be taken out of the integral. However, one can write $C(\mathbf{y}) = C^{(2)} + (C^{(1)} - C^{(2)}) H(y_3)$ where $H(y_3)$ is the Heaviside step function. For $x_3 < 0$, the following equation is thus obtained for the stress field inside the inclusion
$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}^{(2)} \varepsilon_{ln}^{*} C_{jkln}^{(2)} \int_{V^I} \frac{\partial^2 G_{ij}^{Bi}}{\partial y_k \partial x_m}(\mathbf{x}, \mathbf{y}) dy \\ + C_{pqim}^{(2)} \varepsilon_{ln}^{*} (C_{jkln}^{(1)} - C_{jkln}^{(2)}) \int_{V^{I(1)}} \frac{\partial^2 G_{ij}^{Bi}}{\partial y_k \partial x_m}(\mathbf{x}, \mathbf{y}) dy - C_{pqim}^{(2)} \varepsilon_{im}^{*} \quad (16) $$

where $V^{I(1)}$ is the volume of the inclusion located in material (1). Using Eq. 8 and Gauss's theorem, we have then

$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}^{(2)} C_{jkln}^{(2)} \varepsilon_{ln}^{*} \left( E_{ijkm}^{(2)} + \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k dy \right) \\ + C_{pqim}^{(2)} (C_{jkln}^{(1)} - C_{jkln}^{(2)}) \varepsilon_{ln}^{*} \int_{\partial V^{I(1)}} \frac{\partial G_{ij}^{Bi}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k dy - C_{pqim}^{(2)} \varepsilon_{im}^{*} \quad (17) $$

It must be underlined that the surface area $\partial V^{I(1)}$ represents the part of $\partial V^I$ that belongs to material (1) plus the intersection disk between the spherical inclusion and the plane $y_3 = 0$. In this case, the Eshelby tensor for a spherical inclusion in a bi-material, $S^{Esh,Bi}$, reads

$$ S_{imln}^{Esh,Bi} = \frac{1}{2} C_{jkln}^{(2)} (E_{ijkm}^{(2)} + E_{mjki}^{(2)}) \\ + \frac{1}{2} C_{jkln}^{(2)} \left( \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k dy + \int_{\partial V^I} \frac{\partial G_{mj}^{Im}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k dy \right) \\ + \frac{1}{2} (C_{jkln}^{(1)} - C_{jkln}^{(2)}) \left( \int_{\partial V^{I(1)}} \frac{\partial G_{ij}^{Bi}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k dy + \int_{\partial V^{I(1)}} \frac{\partial G_{mj}^{Bi}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k dy \right) \quad (18) $$

For $x_3 > 0$, the same procedure leads to
$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}^{(2)} \varepsilon_{ln}^{*} C_{jkln}^{(2)} \int_{V^I} \frac{\partial^2 G_{ij}^{Bi}}{\partial y_k \partial x_m}(\mathbf{x}, \mathbf{y}) dy \\ + C_{pqim}^{(2)} \varepsilon_{ln}^{*} (C_{jkln}^{(1)} - C_{jkln}^{(2)}) \int_{V^{I(1)}} \frac{\partial^2 G_{ij}^{Bi}}{\partial y_k \partial x_m}(\mathbf{x}, \mathbf{y}) dy - C_{pqim}^{(2)} \varepsilon_{im}^{*} \quad (16) $$

where $V^{I(1)}$ is the volume of the inclusion located in material (1). Using Eq. 8 and Gauss's theorem, we have then

$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}^{(2)} C_{jkln}^{(2)} \varepsilon_{ln}^{*} \left( E_{ijkm}^{(2)} + \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k dy \right) \\ + C_{pqim}^{(2)} (C_{jkln}^{(1)} - C_{jkln}^{(2)}) \varepsilon_{ln}^{*} \int_{\partial V^{I(1)}} \frac{\partial G_{ij}^{Bi}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k dy - C_{pqim}^{(2)} \varepsilon_{im}^{*} \quad (17) $$

It must be underlined that the surface area $\partial V^{I(1)}$ represents the part of $\partial V^I$ that belongs to material (1) plus the intersection disk between the spherical inclusion and the plane $y_3 = 0$. In this case, the Eshelby tensor for a spherical inclusion in a bi-material, $S^{Esh,Bi}$, reads

$$ S_{imln}^{Esh,Bi} = \frac{1}{2} C_{jkln}^{(2)} (E_{ijkm}^{(2)} + E_{mjki}^{(2)}) \\ + \frac{1}{2} C_{jkln}^{(2)} \left( \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k dy + \int_{\partial V^I} \frac{\partial G_{mj}^{Im}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k dy \right) \\ + \frac{1}{2} (C_{jkln}^{(1)} - C_{jkln}^{(2)}) \left( \int_{\partial V^{I(1)}} \frac{\partial G_{ij}^{Bi}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k dy + \int_{\partial V^{I(1)}} \frac{\partial G_{mj}^{Bi}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k dy \right) \quad (18) $$

For $x_3 > 0$, the same procedure leads to
$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}^{(1)} C_{jkln}^{(1)} \varepsilon_{ln}^{*} \left( E_{ijkm}^{(1)} + \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k \mathbf{d}y \right) + C_{pqim}^{(1)} (C_{jkln}^{(2)} - C_{jkln}^{(1)}) \varepsilon_{ln}^{*} \int_{\partial V^{I(2)}} \frac{\partial G_{ij}^{Bi}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k \mathbf{d}y - C_{pqim}^{(1)} \varepsilon_{im}^{*} \quad (19) $$

$$ S_{imln}^{Esh,Bi} = \frac{1}{2} C_{jkln}^{(1)} (E_{ijkm}^{(1)} + E_{mjki}^{(1)}) + \frac{1}{2} C_{jkln}^{(1)} \left( \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k \mathbf{d}y + \int_{\partial V^I} \frac{\partial G_{mj}^{Im}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k \mathbf{d}y \right) + \frac{1}{2} (C_{jkln}^{(2)} - C_{jkln}^{(1)}) \left( \int_{\partial V^{I(2)}} \frac{\partial G_{ij}^{Bi}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k \mathbf{d}y + \int_{\partial V^{I(2)}} \frac{\partial G_{mj}^{Bi}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k \mathbf{d}y \right) \quad (20) $$

### 2.2.3. Cubic elasticity and dilatational eigenstrain

In the case of cubic elasticity and for a dilatational eigenstrain of the form

$$ \boldsymbol{\varepsilon}^* = \begin{bmatrix} \delta^* & 0 & 0 \\ 0 & \delta^* & 0 \\ 0 & 0 & \delta^* \end{bmatrix} \quad (21) $$

it is noteworthy that $(C_{jkln}^{(2)} - C_{jkln}^{(1)}) \varepsilon_{ln}^{*} = 0 \forall j, \forall k$. Hence, Eqs. 17 and 19 reduce seemingly to Eq. 11, namely, for $x_3 < 0$

$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}^{(2)} \left( C_{jkln}^{(2)} \varepsilon_{ln}^{*} \left( E_{ijkm}^{(2)} + \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k \mathbf{d}y \right) - \varepsilon_{im}^{*} \right) \quad (22) $$

and for $x_3 > 0$
$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}^{(1)} C_{jkln}^{(1)} \varepsilon_{ln}^{*} \left( E_{ijkm}^{(1)} + \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} \right) + C_{pqim}^{(1)} (C_{jkln}^{(2)} - C_{jkln}^{(1)}) \varepsilon_{ln}^{*} \int_{\partial V^{I(2)}} \frac{\partial G_{ij}^{Bi}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} - C_{pqim}^{(1)} \varepsilon_{im}^{*} \quad (19) $$

$$ S_{imln}^{Esh,Bi} = \frac{1}{2} C_{jkln}^{(1)} (E_{ijkm}^{(1)} + E_{mjki}^{(1)}) + \frac{1}{2} C_{jkln}^{(1)} \left( \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} + \int_{\partial V^I} \frac{\partial G_{mj}^{Im}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} \right) + \frac{1}{2} (C_{jkln}^{(2)} - C_{jkln}^{(1)}) \left( \int_{\partial V^{I(2)}} \frac{\partial G_{ij}^{Bi}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} + \int_{\partial V^{I(2)}} \frac{\partial G_{mj}^{Bi}}{\partial x_i}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} \right) \quad (20) $$

### 2.2.3. Cubic elasticity and dilatational eigenstrain

In the case of cubic elasticity and for a dilatational eigenstrain of the form

$$ \boldsymbol{\varepsilon}^* = \begin{bmatrix} \delta^* & 0 & 0 \\ 0 & \delta^* & 0 \\ 0 & 0 & \delta^* \end{bmatrix} \quad (21) $$

it is noteworthy that $(C_{jkln}^{(2)} - C_{jkln}^{(1)}) \varepsilon_{ln}^{*} = 0 \forall j, \forall k$. Hence, Eqs. 17 and 19 reduce seemingly to Eq. 11, namely, for $x_3 < 0$

$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}^{(2)} \left( C_{jkln}^{(2)} \varepsilon_{ln}^{*} \left( E_{ijkm}^{(2)} + \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} \right) - \varepsilon_{im}^{*} \right) \quad (22) $$

and for $x_3 > 0$
$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}^{(1)} \left( C_{jkln}^{(1)} \varepsilon_{ln}^{*} \left( E_{ijkm}^{(1)} + \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k d\mathbf{y} \right) - \varepsilon_{im}^{*} \right) \quad (23) $$

However, at the difference of Eq. 11, the computation of $G^{Im}$ involves the use of full-space terms (see Eqs. B.14 and B.20) when the inclusion spans the two materials, which has consequences on the numerical integration speed (see further Section 3.1). Finally, it can be noticed that the stress components $\sigma_{pq}$ scale directly with $\delta^*$ in this case.

## 2.3. Elastic strain energy

The elastic strain energy is defined as

$$ W^e = \frac{1}{2} \int_V \sigma_{ij}(\mathbf{x}) \varepsilon_{ij}^e(\mathbf{x}) d\mathbf{x} \quad (24) $$

where $\varepsilon^e$ is the elastic strain tensor. When eigenstrains are prescribed in $V^I$ and in the absence of any external force and surface constraint, $W^e$ writes (Mura, 1987; Berbenni et al., 2008)

$$ W^e = -\frac{1}{2} \int_{V^I} \sigma_{ij}(\mathbf{x}) \varepsilon_{ij}^*(\mathbf{x}) d\mathbf{x} \quad (25) $$

For a uniform dilatational eigenstrain of the form of Eq. 21 ($\varepsilon_{ij}^* = \delta^* \delta_{ij}$), $W^e$ scales thus with $\delta^{*2}$ since $\sigma_{ij}$ scales also with $\delta^*$ (see Eqs. 22-23).

# 3. Numerical implementation

The numerical integrations over $[0, \pi]$ involved in the calculation of $E_{ijkm}$ and $\frac{\partial G_{ij}^{Im}}{\partial x_m}$ are performed with the function integral of the software MATLAB R2022a. They are used in conjunction with the eig function to solve the Stroh eigenrelation (Eq. B.7). The numerical integrations over $\partial V^I$, $\partial V^{I(1)}$, $\partial V^{I(2)}$ and $V^I$ are obtained from a Gauss-Legendre quadrature method Abramowitz
$$ \sigma_{pq}(\mathbf{x}) = C_{pqim}^{(1)} \left( C_{jkln}^{(1)} \varepsilon_{ln}^{*} \left( E_{ijkm}^{(1)} + \int_{\partial V^I} \frac{\partial G_{ij}^{Im}}{\partial x_m}(\mathbf{x}, \mathbf{y}) n_k \mathbf{dy} \right) - \varepsilon_{im}^{*} \right) \quad (23) $$

However, at the difference of Eq. 11, the computation of $G^{Im}$ involves the use of full-space terms (see Eqs. B.14 and B.20) when the inclusion spans the two materials, which has consequences on the numerical integration speed (see further Section 3.1). Finally, it can be noticed that the stress components $\sigma_{pq}$ scale directly with $\delta^*$ in this case.

### 2.3. Elastic strain energy

The elastic strain energy is defined as

$$ W^e = \frac{1}{2} \int_V \sigma_{ij}(\mathbf{x}) \varepsilon_{ij}^e(\mathbf{x}) \mathbf{d}\mathbf{x} \quad (24) $$

where $\varepsilon^e$ is the elastic strain tensor. When eigenstrains are prescribed in $V^I$ and in the absence of any external force and surface constraint, $W^e$ writes (Mura, 1987; Berbenni et al., 2008)

$$ W^e = -\frac{1}{2} \int_{V^I} \sigma_{ij}(\mathbf{x}) \varepsilon_{ij}^*(\mathbf{x}) \mathbf{d}\mathbf{x} \quad (25) $$

For a uniform dilatational eigenstrain of the form of Eq. 21 ($\varepsilon_{ij}^* = \delta^* \delta_{ij}$), $W^e$ scales thus with $\delta^{*2}$ since $\sigma_{ij}$ scales also with $\delta^*$ (see Eqs. 22-23).

## 3. Numerical implementation

The numerical integrations over $[0, \pi]$ involved in the calculation of $E_{ijkm}$ and $\frac{\partial G_{ij}^{Im}}{\partial x_m}$ are performed with the function integral of the software MATLAB R2022a. They are used in conjunction with the eig function to solve the Stroh eigenrelation (Eq. B.7). The numerical integrations over $\partial V^I$, $\partial V^{I(1)}$, $\partial V^{I(2)}$ and $V^I$ are obtained from a Gauss-Legendre quadrature method Abramowitz
<figure><img src="image_4.png" /><figcaption>Figure 2: Relative error made on the computation of the elastic strain energy (Eq. 25) as a function of $N_S$ for different values of $N_E$ for an inclusion centered at $y_3 = -1.1R$. The reference value is considered at $N_S = 52$ and $N_E = 10$.</figcaption></figure>

and Stegun (1968); von Winckel (2004). The order of the Legendre polynomials are denoted $N_S$ for integrals over surfaces to compute the stresses ($N_S \times N_S$ integration points) and $N_E$ for integrals over volumes to get the elastic strain energy ($N_E \times N_E \times N_E$ integration points).

### 3.1. Convergence study

For the convergence study, a uniform dilatational eigenstrain of the form of Eq. 21 with $\delta^* = 0.1$ is considered along with cubic elastic constants of Cu, $C_{11} = 170$ GPa, $C_{12} = 124$ GPa, $C_{44} = 75$ GPa, a spherical inclusion's radius $R = \frac{a_0\sqrt{2}}{4}$ with the lattice parameter $a_0 = 0.36149$ nm and crystallographic orientations corresponding to 'Orientation A' in Table 1.

Figure 2 shows the relative error made on the computation of the elastic strain energy (Eq. 25) as a function of $N_S$ for an inclusion centered at $y_3 = -1.1R$ and values of $N_E$ ranging from 6 to 10. In this case, the computations of the stresses involves Eq. 11 since the inclusion completely belongs to material (2). It is seen then that a rapid convergence is obtained. A very high accuracy is achieved for $N_S > 20$ with relative errors less than $5 * 10^{-6}$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Relative error made on the computation of the elastic strain energy (Eq. 25) as a function of $N_S$ for different values of $N_E$ for an inclusion centered at $y_3 = -1.1R$. The reference value is considered at $N_S = 52$ and $N_E = 10$.</figcaption></figure>

and Stegun (1968); von Winckel (2004). The order of the Legendre polynomials are denoted $N_S$ for integrals over surfaces to compute the stresses ($N_S \times N_S$ integration points) and $N_E$ for integrals over volumes to get the elastic strain energy ($N_E \times N_E \times N_E$ integration points).

### 3.1. Convergence study

For the convergence study, a uniform dilatational eigenstrain of the form of Eq. 21 with $\delta^* = 0.1$ is considered along with cubic elastic constants of Cu, $C_{11} = 170$ GPa, $C_{12} = 124$ GPa, $C_{44} = 75$ GPa, a spherical inclusion's radius $R = \frac{a_0\sqrt{2}}{4}$ with the lattice parameter $a_0 = 0.36149$ nm and crystallographic orientations corresponding to 'Orientation A' in Table 1.

Figure 2 shows the relative error made on the computation of the elastic strain energy (Eq. 25) as a function of $N_S$ for an inclusion centered at $y_3 = -1.1R$ and values of $N_E$ ranging from 6 to 10. In this case, the computations of the stresses involves Eq. 11 since the inclusion completely belongs to material (2). It is seen then that a rapid convergence is obtained. A very high accuracy is achieved for $N_S > 20$ with relative errors less than $5 * 10^{-6}$.
<figure><img src="image_5.png" /><figcaption>Figure 3: Relative error made on the computation of the elastic strain energy (Eq. 25) as a function of $N_S$ for different values of $N_E$ for an inclusion centered at $y_3 = -0.5R$ (left) and at $y_3 = 0$ (right). The reference values are considered at $N_S = 100$ and $N_E = 10$.</figcaption></figure>

Figure 3 shows then the relative error made on the computation of the elastic strain energy (Eq. 25) for an inclusion centered at $y_3 = -0.5R$ and at $y_3 = 0$. In this case, the computations of the stresses involves Eqs. 22 and 23 since the inclusion spans to the two materials. Numerically, the main difference with Eq. 11 lies in the computation of $G^{Im}$ which involves the use of full-space terms (see Eqs. B.14 and B.20) only when the inclusion spans the two materials. Hence, surface integrals of first derivative of full-space Green tensors must be computed over the part of $\partial V^I$ that belongs to material (1) or (2) and this increases considerably the time needed to achieve an acceptable numerical precision. The convergence is thus much more difficult. However, for $N_S > 60$, the relative error remains below $10^{-3}$ which is acceptable. In both Figures 2 and 3, it can be observed that increasing the value of $N_E$ above 6 has only a marginal effect. Therefore, it is pointless to increase $N_E$ if the stresses within the inclusion are not evaluated with enough accuracy.

## 3.2. Comparison with known solutions in an isotropic half-space

Seo and Mura (1979) have solved the problem of an ellipsoidal inclusion
<figure><img src="image_4.png" /><figcaption>Figure 3: Relative error made on the computation of the elastic strain energy (Eq. 25) as a function of $N_S$ for different values of $N_E$ for an inclusion centered at $y_3 = -0.5R$ (left) and at $y_3 = 0$ (right). The reference values are considered at $N_S = 100$ and $N_E = 10$.</figcaption></figure>

Figure 3 shows then the relative error made on the computation of the elastic strain energy (Eq. 25) for an inclusion centered at $y_3 = -0.5R$ and at $y_3 = 0$. In this case, the computations of the stresses involves Eqs. 22 and 23 since the inclusion spans to the two materials. Numerically, the main difference with Eq. 11 lies in the computation of $G^{Im}$ which involves the use of full-space terms (see Eqs. B.14 and B.20) only when the inclusion spans the two materials. Hence, surface integrals of first derivative of full-space Green tensors must be computed over the part of $\partial V^I$ that belongs to material (1) or (2) and this increases considerably the time needed to achieve an acceptable numerical precision. The convergence is thus much more difficult. However, for $N_S > 60$, the relative error remains below $10^{-3}$ which is acceptable. In both Figures 2 and 3, it can be observed that increasing the value of $N_E$ above 6 has only a marginal effect. Therefore, it is pointless to increase $N_E$ if the stresses within the inclusion are not evaluated with enough accuracy.

## 3.2. Comparison with known solutions in an isotropic half-space

Seo and Mura (1979) have solved the problem of an ellipsoidal inclusion
<figure><img src="image_6.png" /><figcaption>Figure 4: Comparison of energy variation towards a free surface for a spherical inclusion with a uniform dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. The solution of Seo and Mura (1979) is computed from Eq. 26 whereas the solution corresponding to the current work is obtained from Eqs. 25, 22, 27 and 28. $\mu = 23$ GPa, $\nu = 0.4218$ and $R = 0.1278$ nm.</figcaption></figure>

with a uniform dilatational eigenstrain of the form of Eq. 21 located in an isotropic half-space (see also (Mura, 1987)). With our notations, the elastic strain energy for a sphere reads in this case

$$W^e = 2\mu \frac{1 + \nu}{1 - \nu} \delta^{*2} V^I \left( 1 - \frac{1 + \nu}{6} \left( \frac{R}{y_3} \right)^3 \right) \quad (26)$$

where $\nu$ is the Poisson coefficient and $\mu$ the shear modulus. $y_3$ stands for the position of the inclusion center and $y_3 = 0$ corresponds to the free surface. This analytical solution can serve as a reference to check the correctness of the numerical implementation performed in this work. However, in isotropic media, the Stroh eigenvalues which are solutions of Eq. B.7 are not distinct (Bacon et al., 1980). This problem can be circumvented by considering the Stroh integral formalism (Bacon et al., 1980) in which such degeneracies never appear. An alternative approximate approach, that is nevertheless
<figure><img src="image_5.png" /><figcaption>Figure 4: Comparison of energy variation towards a free surface for a spherical inclusion with a uniform dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. The solution of Seo and Mura (1979) is computed from Eq. 26 whereas the solution corresponding to the current work is obtained from Eqs. 25, 22, 27 and 28. $\mu = 23$ GPa, $\nu = 0.4218$ and $R = 0.1278$ nm.</figcaption></figure>

with a uniform dilatational eigenstrain of the form of Eq. 21 located in an isotropic half-space (see also (Mura, 1987)). With our notations, the elastic strain energy for a sphere reads in this case

$$W^e = 2\mu \frac{1 + \nu}{1 - \nu} \delta^{*2} V^I \left( 1 - \frac{1 + \nu}{6} \left( \frac{R}{y_3} \right)^3 \right) \quad (26)$$

where $\nu$ is the Poisson coefficient and $\mu$ the shear modulus. $y_3$ stands for the position of the inclusion center and $y_3 = 0$ corresponds to the free surface. This analytical solution can serve as a reference to check the correctness of the numerical implementation performed in this work. However, in isotropic media, the Stroh eigenvalues which are solutions of Eq. B.7 are not distinct (Bacon et al., 1980). This problem can be circumvented by considering the Stroh integral formalism (Bacon et al., 1980) in which such degeneracies never appear. An alternative approximate approach, that is nevertheless
quite satisfactory from the numerical point of view, is to compute $C_{44}$ as

$$C_{44} = 1.0001 \frac{C_{11} - C_{12}}{2} \qquad (27)$$

In the same way, anisotropic Green's function for half-space do exist (Barnett and Lothe, 1975; Pan and Yuan, 2000). However, for the present purpose of checking the equations implementation, it was chosen to simply set

$$C^{(2)} = 10^{-10} C^{(1)} \qquad (28)$$

Hence, the comparison was made using Eqs. 25, 22, 27 and 28 on one side and Eq. 26 with $\mu = (C_{11} - C_{12})/2$ and $\nu = C_{12}/(2(C_{12} + \mu))$ on the other side. The same values as for the convergence study (Section 3.1) were considered for $C_{11}$, $C_{12}$, $\delta^*$ and $R$. Figure 4 displays the elastic strain energy variation with the distance to the free surface for the two solutions. A perfect numerical match is obtained. It is noteworthy that the energy variation scales as $(R/y_3)^3$ and that the energy always decreases when approaching the free surface meaning that the free surface attracts the inclusion.

Seo and Mura (1979) have also obtained the solutions for the stresses which are equivalent to the results of Mindlin and Cheng (1950) in the case of a sphere. Figure 5 shows the comparisons between the explicit solutions of Mindlin and Cheng (1950) and the current work for the stress distributions along $e_3$ within a spherical inclusion at $y_3 = R$ in an isotropic half-space. Again, a perfect numerical match is obtained. It is noteworthy that the $\sigma_{33}$ exhibits a non-monotonous behavior with a maximum value inside the inclusion.
quite satisfactory from the numerical point of view, is to compute $C_{44}$ as

$$C_{44} = 1.0001 \frac{C_{11} - C_{12}}{2} \qquad (27)$$

In the same way, anisotropic Green's function for half-space do exist (Barnett and Lothe, 1975; Pan and Yuan, 2000). However, for the present purpose of checking the equations implementation, it was chosen to simply set

$$C^{(2)} = 10^{-10} C^{(1)} \qquad (28)$$

Hence, the comparison was made using Eqs. 25, 22, 27 and 28 on one side and Eq. 26 with $\mu = (C_{11} - C_{12})/2$ and $\nu = C_{12}/(2(C_{12} + \mu))$ on the other side. The same values as for the convergence study (Section 3.1) were considered for $C_{11}$, $C_{12}$, $\delta^*$ and $R$. Figure 4 displays the elastic strain energy variation with the distance to the free surface for the two solutions. A perfect numerical match is obtained. It is noteworthy that the energy variation scales as $(R/y_3)^3$ and that the energy always decreases when approaching the free surface meaning that the free surface attracts the inclusion.

Seo and Mura (1979) have also obtained the solutions for the stresses which are equivalent to the results of Mindlin and Cheng (1950) in the case of a sphere. Figure 5 shows the comparisons between the explicit solutions of Mindlin and Cheng (1950) and the current work for the stress distributions along $e_3$ within a spherical inclusion at $y_3 = R$ in an isotropic half-space. Again, a perfect numerical match is obtained. It is noteworthy that the $\sigma_{33}$ exhibits a non-monotonous behavior with a maximum value inside the inclusion.
<figure><img src="image_7.png" /><figcaption>Figure 5: Comparison of stress distributions along *e*₃ (*x*₁ = *x*₂ = 0) within a spherical inclusion with a uniform dilatational eigenstrain ε*ᵢⱼ* = 0.1*δ*ᵢⱼ located at *y*₃ = *R*. The solutions of Mindlin and Cheng (1950) are given explicitly in their paper whereas the solutions corresponding to the current work are obtained from Eqs. 22, 27 and 28. μ = 23 GPa, ν = 0.4218 and *R* = 0.1278 nm.</figcaption></figure>

# 4. Results

## 4.1. Interior stress fields

### 4.1.1. Dilatational eigenstrain

Considering the same parameters as for the convergence study (Section 3.1), i.e. a spherical inclusion with a uniform dilatational eigenstrain of the form of Eq. 21 with *δ*\* = 0.1, cubic elastic constants of Cu and ‘Orientation A’ of Table 1, the stresses within the inclusion for different positions of the inclusion with respect to the interface were computed using Eqs. 22 and 23. The order of the Legendre polynomials was set to *N*<sub>S</sub> = 30 when the inclusion belongs to a same material and to *N*<sub>S</sub> = 100 when the inclusion spans the two materials. Figures 6, 7, 8 and 9 show the variations along *e*₃ of the stress components which are not negligible for the particular crystallographic orientations considered. ‘Orientation A’ is supposed to provide
<figure><img src="image_6.png" /><figcaption>Figure 5: Comparison of stress distributions along *e*₃ (*x*₁ = *x*₂ = 0) within a spherical inclusion with a uniform dilatational eigenstrain ε*ᵢⱼ* = 0.1δ*ᵢⱼ* located at *y*₃ = *R*. The solutions of Mindlin and Cheng (1950) are given explicitly in their paper whereas the solutions corresponding to the current work are obtained from Eqs. 22, 27 and 28. μ = 23 GPa, ν = 0.4218 and *R* = 0.1278 nm.</figcaption></figure>

**4. Results**

4.1. *Interior stress fields*

4.1.1. *Dilatational eigenstrain*

Considering the same parameters as for the convergence study (Sec-
tion 3.1), i.e. a spherical inclusion with a uniform dilatational eigenstrain
of the form of Eq. 21 with δ* = 0.1, cubic elastic constants of Cu and ‘Orien-
tation A’ of Table 1, the stresses within the inclusion for different positions
of the inclusion with respect to the interface were computed using Eqs. 22
and 23. The order of the Legendre polynomials was set to Ns = 30 when the
inclusion belongs to a same material and to Ns = 100 when the inclusion
spans the two materials. Figures 6, 7, 8 and 9 show the variations along
e3 of the stress components which are not negligible for the particular crys-
tallographic orientations considered. ‘Orientation A’ is supposed to provide
strong contrast between the two materials since the directions [111] and [001], respectively the strongest and the weakest direction for the directional Young modulus (Richeton, 2019), are facing each other along **e**_3.

<figure><img src="image_8.png" /><figcaption>Figure 6: Variation of $\sigma_{11}$ along **e**_3 ($x_1 = x_2 = 0$) for a dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium.</figcaption></figure>

The first observation is that the stresses are not uniform within the inclusion at the difference of the well-known Eshelby’s result of uniform stresses for ellipsoidal inclusions with uniform eigenstrain in homogeneous media (Eshelby, 1957). Besides, it is noteworthy that for a dilatational eigenstrain, due to the cubic symmetry, we have $(E_{jkln}^{(2)} - E_{jkln}^{(1)}) \epsilon_{ln}^* = 0 \ \forall j, \forall k$, in addition to $(C_{jkln}^{(2)} - C_{jkln}^{(1)}) \epsilon_{ln}^* = 0 \ \forall j, \forall k$ (see Section 2.2.3). Hence, the stress solutions in a homogenous medium (see Eqs. 22 or 23 with $G^{Im} = 0$) do not depend on the crystallographic orientation. They are the same considering material (2) or material (1) as the reference homogenous medium. Therefore the stress variations cannot be interpreted as smooth transitions from given values in material (2) to different values in material (1). The effect of the imaginary term $G^{Im}$ is much more complex. It is also important to underline that the stress variations are homothetic with respect to the inclusion radius $R$.
strong contrast between the two materials since the directions [111] and [001], respectively the strongest and the weakest direction for the directional Young modulus (Richeton, 2019), are facing each other along **e**_3.

<figure><img src="image_7.png" /><figcaption>Figure 6: Variation of $\sigma_{11}$ along **e**_3 ($x_1 = x_2 = 0$) for a dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium.</figcaption></figure>

The first observation is that the stresses are not uniform within the inclusion at the difference of the well-known Eshelby’s result of uniform stresses for ellipsoidal inclusions with uniform eigenstrain in homogeneous media (Eshelby, 1957). Besides, it is noteworthy that for a dilatational eigenstrain, due to the cubic symmetry, we have $(E_{jkln}^{(2)} - E_{jkln}^{(1)}) \epsilon_{ln}^* = 0 \ \forall j, \forall k$, in addition to $(C_{jkln}^{(2)} - C_{jkln}^{(1)}) \epsilon_{ln}^* = 0 \ \forall j, \forall k$ (see Section 2.2.3). Hence, the stress solutions in a homogenous medium (see Eqs. 22 or 23 with $G^{Im} = 0$) do not depend on the crystallographic orientation. They are the same considering material (2) or material (1) as the reference homogenous medium. Therefore the stress variations cannot be interpreted as smooth transitions from given values in material (2) to different values in material (1). The effect of the imaginary term $G^{Im}$ is much more complex. It is also important to underline that the stress variations are homothetic with respect to the inclusion radius $R$.
<figure><img src="image_9.png" /><figcaption>Figure 7: Variation of $\sigma_{22}$ along $e_3$ ($x_1 = x_2 = 0$) for a dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Variation of $\sigma_{33}$ along $e_3$ ($x_1 = x_2 = 0$) for a dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium.</figcaption></figure>

As it can be expected, when the inclusion belongs to an unique material, the further away from the interface, the more important is the deviation from the homogeneous reference solution. This deviation is asymmetric as it
<figure><img src="image_8.png" /><figcaption>Figure 7: Variation of $\sigma_{22}$ along $e_3$ ($x_1 = x_2 = 0$) for a dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Variation of $\sigma_{33}$ along $e_3$ ($x_1 = x_2 = 0$) for a dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium.</figcaption></figure>

As it can be expected, when the inclusion belongs to an unique material, the further away from the interface, the more important is the deviation from the homogeneous reference solution. This deviation is asymmetric as it
<figure><img src="image_11.png" /><figcaption>Figure 9: Variation of $\sigma_{23}$ along $e_3$ ($x_1 = x_2 = 0$) for a dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium.</figcaption></figure>

is much more pronounced on the side of the interface. Nevertheless, it should be noticed that the stress variations are not necessarily monotonous, see e.g., $\sigma_{11}$ and $\sigma_{22}$ at $y_3 = 1.1R$ and $\sigma_{23}$ at $y_3 = -1.1R$. The non-monotony of some interior stress components was already noticed in isotropic half-spaces (Seo and Mura, 1979) and in 2D dissimilar isotropic bi-materials (Wang et al., 2007). The stress variations can be significant close to the interface (i.e., for $y_3 = \pm 1.1R$) with stress differences greater than 600 MPa and relative variations greater than 5%. However, it seems that these stress differences rapidly decrease as one moves away from the interface as can be noticed from the relative small departures of the Eshelby's homogeneous solutions at $y_3 = \pm 2R$.

When the inclusion spans the two materials, it can be checked that the tractions continuity is well respected (see $\sigma_{33}$ and $\sigma_{23}$). It can be also observed that the stress jumps for $\sigma_{11}$ and $\sigma_{22}$ are more important at $y_3 = \pm 0.5R$ than at $y_3 = 0$. There is also clear departures from the reference homogeneous stress solutions, even at $y_3 = 0$, showing that there is no compensation of
<figure><img src="image_10.png" /><figcaption>Figure 9: Variation of $\sigma_{23}$ along $e_3$ ($x_1 = x_2 = 0$) for a dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium.</figcaption></figure>

is much more pronounced on the side of the interface. Nevertheless, it should be noticed that the stress variations are not necessarily monotonous, see e.g., $\sigma_{11}$ and $\sigma_{22}$ at $y_3 = 1.1R$ and $\sigma_{23}$ at $y_3 = -1.1R$. The non-monotony of some interior stress components was already noticed in isotropic half-spaces (Seo and Mura, 1979) and in 2D dissimilar isotropic bi-materials (Wang et al., 2007). The stress variations can be significant close to the interface (i.e., for $y_3 = \pm 1.1R$) with stress differences greater than 600 MPa and relative variations greater than 5%. However, it seems that these stress differences rapidly decrease as one moves away from the interface as can be noticed from the relative small departures of the Eshelby's homogeneous solutions at $y_3 = \pm 2R$.

When the inclusion spans the two materials, it can be checked that the tractions continuity is well respected (see $\sigma_{33}$ and $\sigma_{23}$). It can be also observed that the stress jumps for $\sigma_{11}$ and $\sigma_{22}$ are more important at $y_3 = \pm 0.5R$ than at $y_3 = 0$. There is also clear departures from the reference homogeneous stress solutions, even at $y_3 = 0$, showing that there is no compensation of
<figure><img src="image_12.png" /><figcaption>Figure 10: Variation of $\sigma_{11}$, $\sigma_{22}$ and $\sigma_{33}$ along $e_1$ ($x_2 = 0$, $x_3 = y_3$) for a dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium.</figcaption></figure>

the imaginary term contributions. Stress variations are never monotonous and are much more disturbed than when the inclusion belongs to a same material. Stress differences within the inclusion can also be very high, e.g., approaching 1 GPa for $\sigma_{22}$ at $y_3 = 0.5R$.

For illustration purposes, variations of stresses along $e_1$ through the center of the spherical inclusion are also shown in Figures 10 and 11 at $y_3 = \pm 2R$ and $y_3 = \pm 1.1R$. In this case too, the stresses are non homogeneous although the distance from the interface does not change. For the particular 'Orientacion A' considered, the stress variations are either symmetric for $\sigma_{11}$, $\sigma_{22}$, $\sigma_{33}$ and $\sigma_{23}$ or anti-symmetric for $\sigma_{31}$ and $\sigma_{12}$. As expected, the variations are much more pronounced at $y_3 = \pm 1.1R$ than at $y_3 = \pm 2R$. They are also globally weaker than along $e_3$, which makes sense since the distance from the interface is fixed.
<figure><img src="image_11.png" /><figcaption>Figure 10: Variation of $\sigma_{11}$, $\sigma_{22}$ and $\sigma_{33}$ along $e_1$ ($x_2 = 0$, $x_3 = y_3$) for a dilatational eigenstrain $\epsilon_{ij}^* = 0.1\delta_{ij}$. Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium.</figcaption></figure>

the imaginary term contributions. Stress variations are never monotonous and are much more disturbed than when the inclusion belongs to a same material. Stress differences within the inclusion can also be very high, e.g., approaching 1 GPa for $\sigma_{22}$ at $y_3 = 0.5R$.

For illustration purposes, variations of stresses along $e_1$ through the center of the spherical inclusion are also shown in Figures 10 and 11 at $y_3 = \pm 2R$ and $y_3 = \pm 1.1R$. In this case too, the stresses are non homogeneous although the distance from the interface does not change. For the particular 'Orientacion A' considered, the stress variations are either symmetric for $\sigma_{11}$, $\sigma_{22}$, $\sigma_{33}$ and $\sigma_{23}$ or anti-symmetric for $\sigma_{31}$ and $\sigma_{12}$. As expected, the variations are much more pronounced at $y_3 = \pm 1.1R$ than at $y_3 = \pm 2R$. They are also globally weaker than along $e_3$, which makes sense since the distance from the interface is fixed.
<figure><img src="image_13.png" /><figcaption>Figure 11: Variation of $\sigma_{23}$, $\sigma_{31}$ and $\sigma_{12}$ along $e_1$ ($x_2 = 0, x_3 = y_3$) for a dilatational eigenstrain (Eq. 21 with $\delta = 0.1$). Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium (i.e., $y_3 = \infty$).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12: Variation of $\sigma_{12}$ along $e_3$ ($x_1 = x_2 = 0$) for a shear eigenstrain (Eq. 29). Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed and dotted lines represent the solutions in the homogeneous medium.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Variation of $\sigma_{23}$, $\sigma_{31}$ and $\sigma_{12}$ along $e_1$ ($x_2 = 0, x_3 = y_3$) for a dilatational eigenstrain (Eq. 21 with $\delta = 0.1$). Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed line represents the solution in the homogeneous medium (i.e., $y_3 = \infty$).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Variation of $\sigma_{12}$ along $e_3$ ($x_1 = x_2 = 0$) for a shear eigenstrain (Eq. 29). Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed and dotted lines represent the solutions in the homogeneous medium.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 13: Variation of $\sigma_{31}$ along $e_3$ ($x_1 = x_2 = 0$) for a shear eigenstrain (Eq. 29). Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed and dotted lines represent the solutions in the homogeneous media.</figcaption></figure>

### 4.1.2. Shear eigenstrain

Still considering the cubic elastic constants of Cu and ‘Orientation A’ of Table 1, a spherical inclusion with a uniform shear eigenstrain of the form of Eq. 29 below is now assumed.

$$ \epsilon^* = \begin{bmatrix} 0 & 0.1 & 0 \\ 0.1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (29) $$

It is noteworthy that, in this case, the homogenous reference solutions are different between materials (1) and (2). Figures 12 and 13 show the variations along $e_3$ of $\sigma_{12}$ and $\sigma_{31}$ (i.e., the two components which are not negligible for the particular crystallographic orientations considered). Same general remarks as for the dilatational eigenstrain hold: non-homogeneous stresses and not necessarily monotonous variations are observed, as well as a rapid decrease of the departure from the homogeneous solution when one moves away from the interface. When the inclusion spans the two materials,
<figure><img src="image_14.png" /><figcaption>Figure 13: Variation of $\sigma_{31}$ along $e_3$ ($x_1 = x_2 = 0$) for a shear eigenstrain (Eq. 29). Different positions ($y_1 = y_2 = 0, y_3$) of the inclusion center are considered on both sides of the interface. The dashed and dotted lines represent the solutions in the homogeneous media.</figcaption></figure>

### 4.1.2. Shear eigenstrain

Still considering the cubic elastic constants of Cu and ‘Orientation A’ of Table 1, a spherical inclusion with a uniform shear eigenstrain of the form of Eq. 29 below is now assumed.

$$ \varepsilon^* = \begin{bmatrix} 0 & 0.1 & 0 \\ 0.1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (29) $$

It is noteworthy that, in this case, the homogeneous reference solutions are different between materials (1) and (2). Figures 12 and 13 show the variations along $e_3$ of $\sigma_{12}$ and $\sigma_{31}$ (i.e., the two components which are not negligible for the particular crystallographic orientations considered). Same general remarks as for the dilatational eigenstrain hold: non-homogeneous stresses and not necessarily monotonous variations are observed, as well as a rapid decrease of the departure from the homogeneous solution when one moves away from the interface. When the inclusion spans the two materials,
the continuity of $\sigma_{31}$ is well respected while the stress jumps displayed by $\sigma_{12}$ are very high with a jump greater than 2.5 GPa at $y_3 = 0.5R$. $\sigma_{31}$ exhibits very strong stress gradients which seem to be controlled by the difference between the homogenous reference solutions of materials (1) and (2).

## 4.2. Variation of elastic strain energy

The variation of the elastic strain energy (Eq. 25) with respect to the position of the inclusion center from one side to the other side of the interface are computed for the different crystallographic orientations considered in Table 1 and different uniform eigenstrains (dilatational (Eq. 21) and shear (Eq. 29)).

While in 'Orientation A', the two extreme directions for the directional Young modulus, [111] and [001], are facing each other along $e_3$, i.e., in a normal way to the interface, these two directions are facing each other along $e_1$, i.e., in a parallel way to the interface, in 'Orientation B'. 'Orientation A' and 'Orientation B' have the same misorientation angle of 56.6°. In 'Orientation C' the direction [111] is along $e_3$ in material (1) and along $e_1$ in material (2). 'Orientation C' yields a misorientation angle of 42.8°. In addition, crystallographic orientations corresponding to symmetrical pure tilt boundaries composed of an infinite array of edge dislocations are also considered (Hirth and Lothe, 1982). In the initial configuration, the Burgers vectors [111] is along $e_3$ while the dislocation line direction [112] is along $e_2$. Starting from this orientation, the stiffness tensors are rotated by $-\theta/2$ and $\theta/2$ around $e_2$ in materials (2) and (1), respectively, where $\theta$ is the tilt or misorientation angle (Hirth and Lothe, 1982). Two cases are computed, $\theta = 11.5°$ and $\theta = 19.1°$ which correspond approximately to $h/b = 5$ and $h/b = 3$, respectively, with $b$ the Burgers vector magnitude and $h$ the spacing between dislocations in the periodic wall.

For these applications, the order of the Legendre polynomials was set to $N_E = 6$ for integrals over volumes. For integral over surfaces, the order was set to $N_S = 20$ when the inclusion belongs to a same material. When the
the continuity of $\sigma_{31}$ is well respected while the stress jumps displayed by $\sigma_{12}$ are very high with a jump greater than 2.5 GPa at $y_3 = 0.5R$. $\sigma_{31}$ exhibits very strong stress gradients which seem to be controlled by the difference between the homogenous reference solutions of materials (1) and (2).

## 4.2. Variation of elastic strain energy

The variation of the elastic strain energy (Eq. 25) with respect to the position of the inclusion center from one side to the other side of the interface are computed for the different crystallographic orientations considered in Table 1 and different uniform eigenstrains (dilatational (Eq. 21) and shear (Eq. 29)).

While in 'Orientation A', the two extreme directions for the directional Young modulus, [111] and [001], are facing each other along $e_3$, i.e., in a normal way to the interface, these two directions are facing each other along $e_1$, i.e., in a parallel way to the interface, in 'Orientation B'. 'Orientation A' and 'Orientation B' have the same misorientation angle of $56.6^\circ$. In 'Orientation C' the direction [111] is along $e_3$ in material (1) and along $e_1$ in material (2). 'Orientation C' yields a misorientation angle of $42.8^\circ$. In addition, crystallographic orientations corresponding to symmetrical pure tilt boundaries composed of an infinite array of edge dislocations are also considered (Hirth and Lothe, 1982). In the initial configuration, the Burgers vectors [$\bar{1}11$] is along $e_3$ while the dislocation line direction [112] is along $e_2$. Starting from this orientation, the stiffness tensors are rotated by $-\theta/2$ and $\theta/2$ around $e_2$ in materials (2) and (1), respectively, where $\theta$ is the tilt or misorientation angle (Hirth and Lothe, 1982). Two cases are computed, $\theta = 11.5^\circ$ and $\theta = 19.1^\circ$ which correspond approximately to $h/b = 5$ and $h/b = 3$, respectively, with $b$ the Burgers vector magnitude and $h$ the spacing between dislocations in the periodic wall.

For these applications, the order of the Legendre polynomials was set to $N_E = 6$ for integrals over volumes. For integral over surfaces, the order was set to $N_S = 20$ when the inclusion belongs to a same material. When the
inclusion spans the two materials, $N_S = 60$ was used for ‘Orientations A, B and C’ (Table 1) and $N_S = 100$ for the two symmetrical tilt boundaries as the energy variation is much weaker in these cases and thus required a higher numerical precision.

<table><caption>Table 1: Different crystallographic orientations considered for the two materials.</caption><thead><tr><th></th><th>Material (1)</th><th>Material (2)</th><th>Misorientation angle</th></tr></thead><tbody><tr><td>Orientation A</td><td>$[111] \parallel e_3$<br>$[11\bar{2}] \parallel e_2$</td><td>$[001] \parallel e_3$<br>$[010] \parallel e_2$</td><td>56.6°</td></tr><tr><td>Orientation B</td><td>$[111] \parallel e_1$<br>$[1\bar{1}2] \parallel e_2$</td><td>$[100] \parallel e_1$<br>$[010] \parallel e_2$</td><td>56.6°</td></tr><tr><td>Orientation C</td><td>$[111] \parallel e_3$<br>$[1\bar{1}2] \parallel e_2$</td><td>$[1\bar{1}1] \parallel e_1$<br>$[11\bar{2}] \parallel e_2$</td><td>42.8°</td></tr></tbody></table>

### 4.2.1. Dilatational eigenstrain

For a uniform dilatational eigenstrain of the form of Eq. 21, $W^e$ scales with $\delta^{*2}$ as explained in Section 2.3 and for cubic elasticty, as for the stresses, the energy in a homogenous medium does not depend on the particular crystallographic orientation. Moreover, if the stress variations are homothetic with respect to the inclusion radius $R$, the elastic strain energy scales with the inclusion volume (cf. Eq. 25). For these reasons, Fig. 14 displays the relative energy variation for the five orientations considered with respect to the homogeneous reference energy. The quantitative aspect of the variation will be discussed later in Section 5. Since the energy should go back to the same value in each material, energy variations are necessarily non-homogeneous as shown by Fig. 14. There are however three to four local extrema depending on the orientation considered. The variations are really significant only within the range $-1.5R < y_3 < 1.5R$. In the range where the inclusion belongs to a single material (i.e., $y_3 < -R$ or $y_3 > R$), it can be conjectured that the departure from the homogeneous reference energy, i.e. the
inclusion spans the two materials, $N_S = 60$ was used for ‘Orientations A, B and C’ (Table 1) and $N_S = 100$ for the two symmetrical tilt boundaries as the energy variation is much weaker in these cases and thus required a higher numerical precision.

<table><caption>Table 1: Different crystallographic orientations considered for the two materials.</caption><thead><tr><th></th><th>Material (1)</th><th>Material (2)</th><th>Misorientation angle</th></tr></thead><tbody><tr><td>Orientation A</td><td>$[111] \parallel e_3$<br>$[1\overline{1}2] \parallel e_2$</td><td>$[001] \parallel e_3$<br>$[0\overline{1}0] \parallel e_2$</td><td>56.6°</td></tr><tr><td>Orientation B</td><td>$[1\overline{1}1] \parallel e_1$<br>$[1\overline{1}2] \parallel e_2$</td><td>$[100] \parallel e_1$<br>$[0\overline{1}0] \parallel e_2$</td><td>56.6°</td></tr><tr><td>Orientation C</td><td>$[1\overline{1}1] \parallel e_3$<br>$[1\overline{1}2] \parallel e_2$</td><td>$[1\overline{1}1] \parallel e_1$<br>$[1\overline{1}2] \parallel e_2$</td><td>42.8°</td></tr></tbody></table>

### 4.2.1. Dilatational eigenstrain

For a uniform dilatational eigenstrain of the form of Eq. 21, $W^e$ scales with $\delta^{*2}$ as explained in Section 2.3 and for cubic elasticty, as for the stresses, the energy in a homogenous medium does not depend on the particular crystallographic orientation. Moreover, if the stress variations are homothetic with respect to the inclusion radius $R$, the elastic strain energy scales with the inclusion volume (cf. Eq. 25). For these reasons, Fig. 14 displays the relative energy variation for the five orientations considered with respect to the homogeneous reference energy. The quantitative aspect of the variation will be discussed later in Section 5. Since the energy should go back to the same value in each material, energy variations are necessarily non-homogeneous as shown by Fig. 14. There are however three to four local extrema depending on the orientation considered. The variations are really significant only within the range $-1.5R < y_3 < 1.5R$. In the range where the inclusion belongs to a single material (i.e., $y_3 < -R$ or $y_3 > R$), it can be conjectured that the departure from the homogeneous reference energy, i.e. the
<figure><img src="image_16.png" /><figcaption>Figure 14: Relative energy variation with respect to the homogeneous reference energy for a dilatational eigenstrain (Eq. 21) when the position of the inclusion center ($y_3$) is displaced from one side to the other side of the interface. On the left, comparison between the three orientations of Table 1. On the right, comparison between the two symmetrical tilt boundaries ($\theta$ being the tilt angle).</figcaption></figure>

decrease of the imaginary term, probably follows the same scaling of $\left(\frac{R}{y_3}\right)^3$ (cf. Eq. 26) as the imaginary term in an isotropic half-space (Seo and Mura, 1979; Mura, 1987).

For the symmetrical tilt boundaries, the variations are perfectly symmet- ric with a global minimum at $y_3 = 0$ that decreases with the tilt angle. For the orientations of Table 1, a rapid transition is observed, around $y_3 = 0$, from the global minimum on one side of the interface to the global maximum on the other side of the interface. It is noteworthy that the global minimum and maximum always correspond to configurations where the inclusion spans the two materials. ‘Orientations A’ exhibits the lowest minimum, followed by ‘Orientations B’ and ‘Orientations C’, which agrees with the supposed elastic contrast for a dilatation displayed by these orientations.
<figure><img src="image_15.png" /><figcaption>Figure 14: Relative energy variation with respect to the homogeneous reference energy for a dilatational eigenstrain (Eq. 21) when the position of the inclusion center ($y_3$) is displaced from one side to the other side of the interface. On the left, comparison between the three orientations of Table 1. On the right, comparison between the two symmetrical tilt boundaries ($\theta$ being the tilt angle).</figcaption></figure>

decrease of the imaginary term, probably follows the same scaling of $\left(\frac{R}{y_3}\right)^3$ (cf. Eq. 26) as the imaginary term in an isotropic half-space (Seo and Mura, 1979; Mura, 1987).

For the symmetrical tilt boundaries, the variations are perfectly symmet- ric with a global minimum at $y_3 = 0$ that decreases with the tilt angle. For the orientations of Table 1, a rapid transition is observed, around $y_3 = 0$, from the global minimum on one side of the interface to the global maximum on the other side of the interface. It is noteworthy that the global minimum and maximum always correspond to configurations where the inclusion spans the two materials. ‘Orientations A’ exhibits the lowest minimum, followed by ‘Orientations B’ and ‘Orientations C’, which agrees with the supposed elastic contrast for a dilatation displayed by these orientations.
<figure><img src="image_17.png" /><figcaption>Figure 15: Energy variation for the shear eigenstrain of Eq. 29 when the position of the inclusion center ($y_3$) is displaced from one side to the other side of the interface. On the left, comparison between the three orientations of Table 1. On the right, comparison between the two symmetrical tilt boundaries ($\theta$ being the tilt angle).</figcaption></figure>

### 4.2.2. Shear eigenstrain

For the shear eigenstrain of Eq. 29, the energy in a homogenous medium depends on the particular crystallographic orientation. It is however identical in materials (1) and (2) when the latter refer to the orientations on both side of a symmetrical tilt boundary. Therefore, for ‘Orientations A, B and C’, the energy variation does not go through an extremum value. There is a smooth transition from the reference homogeneous energy in one material to the reference homogeneous energy in the other material over the range $-1.5R < y_3 < 1.5R$. For the symmetrical tilt boundaries, the energy variations are perfectly symmetric and go through a single minimum at $y_3 = 0$ that decreases with the tilt angle.

## 5. Discussion

It should be first pointed out that the developed computational framework is quite general. It allows to consider any kind of linear elasticity in
<figure><img src="image_16.png" /><figcaption>Figure 15: Energy variation for the shear eigenstrain of Eq. 29 when the position of the inclusion center ($y_3$) is displaced from one side to the other side of the interface. On the left, comparison between the three orientations of Table 1. On the right, comparison between the two symmetrical tilt boundaries ($\theta$ being the tilt angle).</figcaption></figure>

### 4.2.2. Shear eigenstrain

For the shear eigenstrain of Eq. 29, the energy in a homogenous medium depends on the particular crystallographic orientation. It is however identical in materials (1) and (2) when the latter refer to the orientations on both side of a symmetrical tilt boundary. Therefore, for ‘Orientations A, B and C’, the energy variation does not go through an extremum value. There is a smooth transition from the reference homogeneous energy in one material to the reference homogeneous energy in the other material over the range $-1.5R < y_3 < 1.5R$. For the symmetrical tilt boundaries, the energy variations are perfectly symmetric and go through a single minimum at $y_3 = 0$ that decreases with the tilt angle.

## 5. Discussion

It should be first pointed out that the developed computational framework is quite general. It allows to consider any kind of linear elasticity in
<figure><img src="image_18.png" /><figcaption>Figure 16: On the left, maximum energy variation due to image effect with respect to $\delta^*$ considering a dilatational eigenstrain $\epsilon_{ij}^* = \delta^*\delta_{ij}$ and the three orientations of Table 1, as well as the two symmetrical tilt boundaries defined by their tilt angle $\theta$. These variations are compared to the interaction energy (cf. Eq.31) between the same eigenstrain and an external stress field defined by the value of its trace, $\sigma_{kk}^{ext} = 1$ GPa. On the right is plotted the equivalent value of $\sigma_{kk}^{ext}$ such that the maximum energy variation due to the image effect equals the interaction energy between the dilatational eigenstrain and $\sigma_{kk}^{ext}$.</figcaption></figure>

each material. Hence, it can be applied to isotropic or anisotropic dissimilar materials (e.g., Cu-Ag) and is not restricted to crystalline materials. The shape of the eigenstrain tensor is also free. It can thus accounts for a thermal dilatation, a plastic strain or stands for an inhomogeneity through the equivalent inclusion method. In particular, the dilatational eigenstrain tensor considered in Eq. 21 was shown to be well-suited for modeling the effects of vacancies or solute atoms in substitution in cubic lattices (White and Coghlan, 1977; Cai et al., 2014; Dingreville and Berbenni, 2016; Clouet et al., 2018) although their size fall in the atomic scale with radius typically in the range of 0.1 nm to 0.2 nm. In the context of an inclusion representing a solute atom, interior stress variations have therefore no clear physical meaning but the latter are however needed to compute accurately the elastic strain energy through Eq. 25. Indeed, as seen from Eq. 24, the elastic strain energy depends physically almost exclusively on the exterior elastic fields. It
<figure><img src="image_17.png" /><figcaption>Figure 16: On the left, maximum energy variation due to image effect with respect to $\delta^*$ considering a dilatational eigenstrain $\epsilon_{ij}^* = \delta^*\delta_{ij}$ and the three orientations of Table 1, as well as the two symmetrical tilt boundaries defined by their tilt angle $\theta$. These variations are compared to the interaction energy (cf. Eq.31) between the same eigenstrain and an external stress field defined by the value of its trace, $\sigma_{kk}^{ext} = 1$ GPa. On the right is plotted the equivalent value of $\sigma_{kk}^{ext}$ such that the maximum energy variation due to the image effect equals the interaction energy between the dilatational eigenstrain and $\sigma_{kk}^{ext}$.</figcaption></figure>

each material. Hence, it can be applied to isotropic or anisotropic dissimilar materials (e.g., Cu-Ag) and is not restricted to crystalline materials. The shape of the eigenstrain tensor is also free. It can thus accounts for a thermal dilatation, a plastic strain or stands for an inhomogeneity through the equivalent inclusion method. In particular, the dilatational eigenstrain tensor considered in Eq. 21 was shown to be well-suited for modeling the effects of vacancies or solute atoms in substitution in cubic lattices (White and Coghlan, 1977; Cai et al., 2014; Dingreville and Berbenni, 2016; Clouet et al., 2018) although their size fall in the atomic scale with radius typically in the range of 0.1 nm to 0.2 nm. In the context of an inclusion representing a solute atom, interior stress variations have therefore no clear physical meaning but the latter are however needed to compute accurately the elastic strain energy through Eq. 25. Indeed, as seen from Eq. 24, the elastic strain energy depends physically almost exclusively on the exterior elastic fields. It
is thus assumed to be well predicted by the continuum linear elastic theory. Besides, examining the segregation of vacancies or solute atoms toward grain boundaries requires considering various other important factors. First, the interaction with the stress field resulting from the grain boundary should be taken into account as well. Two primary contributions are typically associated with this interaction (White and Coghlan, 1977; Clouet et al., 2018): a size effect linked to the dilatation considered in the eigenstrain and a modulus effect stemming from the fact that the presence of solutes alters the elastic constants of the material. Secondly, grain boundaries possess small, albeit finite, width (of the order of the nanometer), unlike the sharp interface examined in this study. Therefore, it is certainly worthwhile to consider a specific elastic stiffness tensor over the grain boundary width, as done in previous studies (Chen et al., 2019, 2020, 2021). For the symmetrical tilt boundaries, a smooth transition from $C^{(2)}$ to $C^{(1)}$ is nevertheless rather expected at the scale of a few nanometers. Additionally, elasticity may be nonlinear within certain regions of the grain boundary. Thus, the configuration where the inclusion spans the two materials does not reflect a direct reality but it might provide insights into the impact of elastic heterogeneity when the defect is within the grain boundary, especially as the defect size becomes comparable to the grain boundary size, as with large precipitates. In the case of solute atoms which are typically much smaller than the grain boundary width, the configuration where the inclusion spans the two materials may rather reflect a situation where the solute is in the transition area between the bulk crystal and the grain boundary.

Furthermore, it might be valuable to quantitatively compare the image effect felt by a dilatation center in heterogeneous anisotropic elasticity (Fig. 14) to the typical size effect associated with the interaction energy between a dilatation center and a grain boundary stress field. The latter can be determined from the general elastic interaction energy between an eigenstrain $\epsilon^*$ distributed in $V^1$ and an external stress field $\sigma^{ext}$ (such as the one arising
is thus assumed to be well predicted by the continuum linear elastic theory. Besides, examining the segregation of vacancies or solute atoms toward grain boundaries requires considering various other important factors. First, the interaction with the stress field resulting from the grain boundary should be taken into account as well. Two primary contributions are typically associated with this interaction (White and Coghlan, 1977; Clouet et al., 2018): a size effect linked to the dilatation considered in the eigenstrain and a modulus effect stemming from the fact that the presence of solutes alters the elastic constants of the material. Secondly, grain boundaries possess small, albeit finite, width (of the order of the nanometer), unlike the sharp interface examined in this study. Therefore, it is certainly worthwhile to consider a specific elastic stiffness tensor over the grain boundary width, as done in previous studies (Chen et al., 2019, 2020, 2021). For the symmetrical tilt boundaries, a smooth transition from $C^{(2)}$ to $C^{(1)}$ is nevertheless rather expected at the scale of a few nanometers. Additionally, elasticity may be nonlinear within certain regions of the grain boundary. Thus, the configuration where the inclusion spans the two materials does not reflect a direct reality but it might provide insights into the impact of elastic heterogeneity when the defect is within the grain boundary, especially as the defect size becomes comparable to the grain boundary size, as with large precipitates. In the case of solute atoms which are typically much smaller than the grain boundary width, the configuration where the inclusion spans the two materials may rather reflect a situation where the solute is in the transition area between the bulk crystal and the grain boundary.

Furthermore, it might be valuable to quantitatively compare the image effect felt by a dilatation center in heterogeneous anisotropic elasticity (Fig. 14) to the typical size effect associated with the interaction energy between a dilatation center and a grain boundary stress field. The latter can be determined from the general elastic interaction energy between an eigenstrain $\epsilon^*$ distributed in $V^1$ and an external stress field $\sigma^{ext}$ (such as the one arising
from a grain boundary), as given by (Mura, 1987)

$$W^{int} = - \int_{V^I} \sigma_{ij}^{ext}(\mathbf{x}) \varepsilon_{ij}^{*}(\mathbf{x}) d\mathbf{x} \qquad (30)$$

Considering that the volume $V^I$ is sufficiently small to disregard variations of the external stress field within it, we can write for a uniform dilatational eigenstrain of the form $\varepsilon_{ij}^* = \delta^* \delta_{ij}$

$$W^{int} = -V^I \delta^* \sigma_{kk}^{ext} \qquad (31)$$

Figure 16 displays on the left the maximum energy variation per unit volume due to the image effect with respect to $\delta^*$ considering the three orientations of Table 1 and the two symmetrical tilt boundaries. These variations are directly inferred from Fig. 14. They are then compared to the interaction energy per unit volume computed from Eq. 31 between the same dilatational eigenstrain and an external stress field with a trace value of 1 GPa. As anticipated, the image effect scales with $\delta^{*2}$, while the interaction energy with a stress field scales linearly with $\delta^*$. In the current scenario with $\sigma_{kk}^{ext} = 1$ GPa, the image effect becomes predominant for ‘Orientations A and B’ when $\delta^* > 0.1$. To quantitatively evaluate the significance of the image effect, Figure 16 then presents, on the right, the equivalent value of $\sigma_{kk}^{ext}$ such that the maximum energy variation due to the image effect equals the interaction energy between the dilatational eigenstrain and $\sigma_{kk}^{ext}$ according to Eq. 31. These $\sigma_{kk}^{ext}$ values scale therefore linearly with $\delta^*$. While these values remain very low for the symmetrical tilt boundaries, they can be comparable to typical values of grain boundary stress fields (Dingreville and Berbenni, 2016) for the other orientations considered (cf. Table 1). Furthermore, Figure 16 emphasizes that the larger the dilatation of the inclusion, the more likely the image effect is to be significant compared to the interaction energy with a stress field (Eq. 31) due to the different scaling. Another distinction should be pointed out. The image effect remains indeed the same at a
from a grain boundary), as given by (Mura, 1987)

$$W^{int} = - \int_{V^I} \sigma_{ij}^{ext}(\mathbf{x}) \varepsilon_{ij}^{*}(\mathbf{x}) d\mathbf{x} \qquad (30)$$

Considering that the volume $V^I$ is sufficiently small to disregard variations of the external stress field within it, we can write for a uniform dilatational eigenstrain of the form $\varepsilon_{ij}^* = \delta^* \delta_{ij}$

$$W^{int} = -V^I \delta^* \sigma_{kk}^{ext} \qquad (31)$$

Figure 16 displays on the left the maximum energy variation per unit volume due to the image effect with respect to $\delta^*$ considering the three orientations of Table 1 and the two symmetrical tilt boundaries. These variations are directly inferred from Fig. 14. They are then compared to the interaction energy per unit volume computed from Eq. 31 between the same dilatational eigenstrain and an external stress field with a trace value of 1 GPa. As anticipated, the image effect scales with $\delta^{*2}$, while the interaction energy with a stress field scales linearly with $\delta^*$. In the current scenario with $\sigma_{kk}^{ext} = 1$ GPa, the image effect becomes predominant for ‘Orientations A and B’ when $\delta^* > 0.1$. To quantitatively evaluate the significance of the image effect, Figure 16 then presents, on the right, the equivalent value of $\sigma_{kk}^{ext}$ such that the maximum energy variation due to the image effect equals the interaction energy between the dilatational eigenstrain and $\sigma_{kk}^{ext}$ according to Eq. 31. These $\sigma_{kk}^{ext}$ values scale therefore linearly with $\delta^*$. While these values remain very low for the symmetrical tilt boundaries, they can be comparable to typical values of grain boundary stress fields (Dingreville and Berbenni, 2016) for the other orientations considered (cf. Table 1). Furthermore, Figure 16 emphasizes that the larger the dilatation of the inclusion, the more likely the image effect is to be significant compared to the interaction energy with a stress field (Eq. 31) due to the different scaling. Another distinction should be pointed out. The image effect remains indeed the same at a
given distance from the interface. On the contrary, the interaction energy varies with the grain boundary stress field (Eq. 31), thus exhibiting strong variations along the grain boundary that range from negative to positive values (White and Coghlan, 1977; Dingreville and Berbenni, 2016).

In addition to the dilatation value, several other factors influence the significance of the image effect, such as the degree of elastic anisotropy and the extent of elastic heterogeneity. Elastic heterogeneity can be associated with misorientations between grains but also to the presence of different phases (heterophase interfaces) or to materials of different kinds as in composites. The present results were obtained for Cu, which has a relatively high Zener ratio with $A = \frac{2C_{44}}{C_{11} - C_{12}} = 3.26$. However, the full space of cubic orientations was not investigated in this study and other orientations might provide larger image effects. Besides, certain cubic materials can exhibit much higher values of A, like lithium metal with $A \approx 8$ (Aspinall et al., 2022), a crucial component in batteries, or some shape memory alloys like CuAlBe with $A \approx 13$ (Rios-Jara et al., 1991). Nevertheless, the image effect is expected to increase asymptotically with A since the maximum image effect should occur in presence of free or rigid surfaces as it was computed for dislocations (Chen et al., 2021).

The exploration of non-cubic elasticity could be intriguing also. First, in a homogeneous medium, the energy due to a dilatational eigenstrain will depend on the specific crystallographic orientation. Second, the dilatational eigenstrain tensor used to simulate the role of point defects as dilatation centers should not be considered isotropic, as in Eq. 21, but should be anisotropic with the same or lower symmetry as the host lattice (Nye, 1957; Balluffi, 2016). In other words, the solute atoms should be viewed as anisotropic centers of dilatation (Cochardt et al., 1955; Siems, 1968; Clouet et al., 2018).

Finally, it may be worth underlying that the image effect might impact interface motion. Indeed, when an interface, such as a grain boundary, starts moving in presence of a distribution of defects (solutes, precipitates or va-
given distance from the interface. On the contrary, the interaction energy varies with the grain boundary stress field (Eq. 31), thus exhibiting strong variations along the grain boundary that range from negative to positive values (White and Coghlan, 1977; Dingreville and Berbenni, 2016).

In addition to the dilatation value, several other factors influence the significance of the image effect, such as the degree of elastic anisotropy and the extent of elastic heterogeneity. Elastic heterogeneity can be associated with misorientations between grains but also to the presence of different phases (heterophase interfaces) or to materials of different kinds as in composites. The present results were obtained for Cu, which has a relatively high Zener ratio with $A = \frac{2C_{44}}{C_{11} - C_{12}} = 3.26$. However, the full space of cubic orientations was not investigated in this study and other orientations might provide larger image effects. Besides, certain cubic materials can exhibit much higher values of A, like lithium metal with $A \approx 8$ (Aspinall et al., 2022), a crucial component in batteries, or some shape memory alloys like CuAlBe with $A \approx 13$ (Rios-Jara et al., 1991). Nevertheless, the image effect is expected to increase asymptotically with A since the maximum image effect should occur in presence of free or rigid surfaces as it was computed for dislocations (Chen et al., 2021).

The exploration of non-cubic elasticity could be intriguing also. First, in a homogeneous medium, the energy due to a dilatational eigenstrain will depend on the specific crystallographic orientation. Second, the dilatational eigenstrain tensor used to simulate the role of point defects as dilatation centers should not be considered isotropic, as in Eq. 21, but should be anisotropic with the same or lower symmetry as the host lattice (Nye, 1957; Balluffi, 2016). In other words, the solute atoms should be viewed as anisotropic centers of dilatation (Cochardt et al., 1955; Siems, 1968; Clouet et al., 2018).

Finally, it may be worth underlying that the image effect might impact interface motion. Indeed, when an interface, such as a grain boundary, starts moving in presence of a distribution of defects (solutes, precipitates or va-
cancies), the elastic strain energy will be altered because it depends on the
relative distance of the defects to the interface in heterogeneous elasticity.
Hence, the elastic strain energy can increase or decrease and thus can ei-
ther promote or hinder interface motion. The significance of such an image
effect on interface motion was previously studied for distributions of disloca-
tions (Richeton et al., 2020), providing insights into the connection between
interface motion and crystallographic disorientations.

**6. Conclusions**

In this work, expressions of the interior stress fields of a spherical inclusion with a uniform eigenstrain embedded in an anisotropic bi-material are derived. They involve surface integrals of the imaginary term of the first derivative of the Green tensor in an anisotropic bi-material with a planar interface. The latter can be computed via the Stroh formalism from a regular line-integral over [0, π] as proposed by Pan and Yuan (2000). The entire expression is thus amenable to standard numerical integrations such as Gauss-Legendre quadrature methods. Specific expressions are provided based on whether the inclusion belongs to the same or the two materials. Additionally, expressions for the equivalent Eshelby tensor are presented. It is noteworthy that straightforward derivations from Eq. 5 allow for expressions of both external and internal stress fields for inclusions of any shape with non-uniform eigenstrain in anisotropic bi-materials. However, the exploration of efficient numerical integration methods for these cases remains an avenue for future investigation.

The stress fields within the spherical inclusion, along with the elastic
strain energy, are then computed, taking into account the cubic elastic con-
stants of Cu. Various inclusion positions relative to the interface, different
forms of eigenstrain (dilatational (Eq.21) and pure shear (Eq.29)), and dif-
ferent crystallographic orientations are considered. The following key con-
clusions emerge:
cancies), the elastic strain energy will be altered because it depends on the
relative distance of the defects to the interface in heterogeneous elasticity.
Hence, the elastic strain energy can increase or decrease and thus can ei-
ther promote or hinder interface motion. The significance of such an image
effect on interface motion was previously studied for distributions of disloca-
tions (Richeton et al., 2020), providing insights into the connection between
interface motion and crystallographic disorientations.

6. Conclusions

In this work, expressions of the interior stress fields of a spherical inclusion with a uniform eigenstrain embedded in an anisotropic bi-material are derived. They involve surface integrals of the imaginary term of the first derivative of the Green tensor in an anisotropic bi-material with a planar interface. The latter can be computed via the Stroh formalism from a regular line-integral over [0, π] as proposed by Pan and Yuan (2000). The entire expression is thus amenable to standard numerical integrations such as Gauss-Legendre quadrature methods. Specific expressions are provided based on whether the inclusion belongs to the same or the two materials. Additionally, expressions for the equivalent Eshelby tensor are presented. It is noteworthy that straightforward derivations from Eq. 5 allow for expressions of both external and internal stress fields for inclusions of any shape with non-uniform eigenstrain in anisotropic bi-materials. However, the exploration of efficient numerical integration methods for these cases remains an avenue for future investigation.

The stress fields within the spherical inclusion, along with the elastic
strain energy, are then computed, taking into account the cubic elastic con-
stants of Cu. Various inclusion positions relative to the interface, different
forms of eigenstrain (dilatational (Eq.21) and pure shear (Eq.29)), and dif-
ferent crystallographic orientations are considered. The following key con-
clusions emerge:
* In a bi-material, the stresses within the inclusion are non-homogeneous, and their variations are not necessarily monotonous even when the inclusion belongs to a same material. A rapid decrease in the imaginary term is observed as one moves away from the interface.

* With a dilatational eigenstrain, the variation of elastic strain energy can exhibit several extrema. The global minimum and maximum consistently correspond to configurations where the inclusion spans the two materials.

* For symmetrical tilt boundaries, energy variations are perfectly symmetric for both dilatational and shear eigenstrains, exhibiting a global minimum located on the interface that decreases with the tilt angle.

* The variations of energy are really significant only within the range $-1.5R < y_3 < 1.5R$ ($R$ being the radius of the inclusion and $y_3$ the position of the inclusion center relative to the interface). Once the inclusion belongs to a same material, the decrease of the imaginary term probably follows the same scaling of $(\frac{R}{y_3})^3$ (cf. Eq. 26) as the imaginary term in an isotropic half-space (Bacon, 1972; Seo and Mura, 1979; Mura, 1987).

* The significance of the image effect related to the described energy variations is quantitatively assessed from comparisons with the interaction energy between a dilatational eigenstrain ($\epsilon_{ij}^* = \delta^*\delta_{ij}$) and an external stress field (Eq. 31). The image effect remains almost insignificant for the symmetrical tilt boundaries studied but can be comparable to typical interaction energy values for the other orientations considered (cf. Table 1). Moreover, the bigger the dilatation of the inclusion, the higher the probability that the image effect is significant due to the different scaling (respectively, as $\delta^{*2}$ for the image effect and as $\delta^*$ for the interaction energy with a stress field).
* In a bi-material, the stresses within the inclusion are non-homogeneous, and their variations are not necessarily monotonous even when the inclusion belongs to a same material. A rapid decrease in the imaginary term is observed as one moves away from the interface.

* With a dilatational eigenstrain, the variation of elastic strain energy can exhibit several extrema. The global minimum and maximum consistently correspond to configurations where the inclusion spans the two materials.

* For symmetrical tilt boundaries, energy variations are perfectly symmetric for both dilatational and shear eigenstrains, exhibiting a global minimum located on the interface that decreases with the tilt angle.

* The variations of energy are really significant only within the range $-1.5R < y_3 < 1.5R$ ($R$ being the radius of the inclusion and $y_3$ the position of the inclusion center relative to the interface). Once the inclusion belongs to a same material, the decrease of the imaginary term probably follows the same scaling of $(\frac{R}{y_3})^3$ (cf. Eq. 26) as the imaginary term in an isotropic half-space (Bacon, 1972; Seo and Mura, 1979; Mura, 1987).

* The significance of the image effect related to the described energy variations is quantitatively assessed from comparisons with the interaction energy between a dilatational eigenstrain ($\epsilon_{ij}^* = \delta^*\delta_{ij}$) and an external stress field (Eq. 31). The image effect remains almost insignificant for the symmetrical tilt boundaries studied but can be comparable to typical interaction energy values for the other orientations considered (cf. Table 1). Moreover, the bigger the dilatation of the inclusion, the higher the probability that the image effect is significant due to the different scaling (respectively, as $\delta^{*2}$ for the image effect and as $\delta^*$ for the interaction energy with a stress field).
## 7. Acknowledgments

The author is grateful to the French National Research Agency (ANR) for financial support under contract ANR-22-CE92-0058-01 (SILA). The author acknowledges also the use of the Cassiopee Arts et Métiers Institute of Technology HPC Center made available for conducting the research reported in this paper.

## Appendix A. Full-space Green tensor

The full-space anisotropic Green tensor can be computed as (Barnett, 1972; Lazar, 2016)

$$G_{ij}^{\infty}(\mathbf{x}, \mathbf{y}) = \frac{1}{8\pi^2|\mathbf{x} - \mathbf{y}|} \int_{0}^{2\pi} (nCn)_{ij}^{-1} d\phi \quad (A.1)$$

where $\mathbf{n}$ is a unit-vector perpendicular to $\mathbf{x} - \mathbf{y}$ and $(nCn)_{ij} = n_k C_{ikjl} n_l$.

### Appendix A.1. First derivative of the full-space Green tensor

The first derivative of the full-space anisotropic Green tensor can be evaluated as (Barnett, 1972; Lazar, 2016)

$$\frac{\partial G_{ij}^{\infty}(\mathbf{x}, \mathbf{y})}{\partial x_k} = -\frac{1}{8\pi^2|\mathbf{x} - \mathbf{y}|^2} \int_{0}^{2\pi} \left( \tau_k (nCn)_{ij}^{-1} - n_k F_{ij} \right) d\phi \quad (A.2)$$

where $\tau = \frac{\mathbf{x} - \mathbf{y}}{|\mathbf{x} - \mathbf{y}|}$ is a unit-vector and $F_{ij} = (nCn)_{ip}^{-1} [(nC\tau)_{pq} + (\tau Cn)_{pq}] (nCn)_{qj}^{-1}$.

### Appendix A.2. Second derivative of the full-space Green tensor

The second derivative of the full-space anisotropic Green tensor can be obtained as the sum of a singular term and a $\frac{1}{|\mathbf{x} - \mathbf{y}|^3}$ term (Kröner, 1990; Lazar, 2016)

$$\frac{\partial^2 G_{ij}^{\infty}(\mathbf{x}, \mathbf{y})}{\partial x_k \partial x_m} = -\delta(\mathbf{x} - \mathbf{y}) E_{ijkm} + \frac{1}{|\mathbf{x} - \mathbf{y}|^3} H_{ijkm} \quad (A.3)$$
## 7. Acknowledgments

The author is grateful to the French National Research Agency (ANR) for financial support under contract ANR-22-CE92-0058-01 (SILA). The author acknowledges also the use of the Cassiopee Arts et Métiers Institute of Technology HPC Center made available for conducting the research reported in this paper.

# Appendix A. Full-space Green tensor

The full-space anisotropic Green tensor can be computed as (Barnett, 1972; Lazar, 2016)

$$G_{ij}^{\infty}(\mathbf{x}, \mathbf{y}) = \frac{1}{8\pi^2|\mathbf{x} - \mathbf{y}|} \int_{0}^{2\pi} (nCn)_{ij}^{-1} d\phi \quad (A.1)$$

where $\mathbf{n}$ is a unit-vector perpendicular to $\mathbf{x} - \mathbf{y}$ and $(nCn)_{ij} = n_k C_{ikjl} n_l$.

## Appendix A.1. First derivative of the full-space Green tensor

The first derivative of the full-space anisotropic Green tensor can be evaluated as (Barnett, 1972; Lazar, 2016)

$$\frac{\partial G_{ij}^{\infty}(\mathbf{x}, \mathbf{y})}{\partial x_k} = -\frac{1}{8\pi^2|\mathbf{x} - \mathbf{y}|^2} \int_{0}^{2\pi} \left( \tau_k (nCn)_{ij}^{-1} - n_k F_{ij} \right) d\phi \quad (A.2)$$

where $\tau = \frac{\mathbf{x} - \mathbf{y}}{|\mathbf{x} - \mathbf{y}|}$ is a unit-vector and $F_{ij} = (nCn)_{ip}^{-1} [(nC\tau)_{pq} + (\tau Cn)_{pq}] (nCn)_{qj}^{-1}$.

## Appendix A.2. Second derivative of the full-space Green tensor

The second derivative of the full-space anisotropic Green tensor can be obtained as the sum of a singular term and a $\frac{1}{|\mathbf{x} - \mathbf{y}|^3}$ term (Kröner, 1990; Lazar, 2016)

$$\frac{\partial^2 G_{ij}^{\infty}(\mathbf{x}, \mathbf{y})}{\partial x_k \partial x_m} = -\delta(\mathbf{x} - \mathbf{y}) E_{ijkm} + \frac{1}{|\mathbf{x} - \mathbf{y}|^3} H_{ijkm} \quad (A.3)$$
The tensor $\mathbf{E}$ is given as (Lazar, 2016)

$$E_{ijkm} = \frac{1}{4\pi} \int_{0}^{\pi} \int_{0}^{2\pi} n_m n_k (nCn)_{ij}^{-1} \sin\theta \, d\theta \, d\phi \quad (A.4)$$

where $\mathbf{n}$ is here the external unit normal of the unit sphere, i.e. $n = [\sin\theta \cos\phi, \sin\theta \sin\phi, \cos\theta]$. The tensor $\mathbf{H}$ is not needed for the present work. Its expression can be found in (Barnett, 1972; Lazar, 2016).

## Appendix B. Imaginary part of the anisotropic bi-material Green tensor

The purpose of this appendix is to explain briefly the general approach of Pan and Yuan (2000) to obtain the anisotropic elastic Green tensor in a bi-material and to provide the equations needed to perform the numerical calculations. All the details can be found in (Pan and Yuan, 2000).

A full-space bi-material is considered where the region $x_3 > 0$ is occupied by material (1) and the region $x_3 < 0$ by material (2). A point force is applied at $(y_1, y_2, y_3)$. To find the displacement vector $\mathbf{u}$ and thus the Green function tensor $\mathbf{G}^{Bi}$, a two-dimensional Fourier transform is applied

$$\tilde{u}_k(\xi_1, \xi_2, x_3) = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} u_k(x_1, x_2, x_3) e^{i(x_1\xi_1 + x_2\xi_2)} dx_1 dx_2 \quad (B.1)$$

In the Fourier transformed domain and in the absence of body force, the equilibrium equation then writes in each material

$$C_{i\alpha k\beta} \xi_{\alpha} \xi_{\beta} \tilde{u}_k + i(C_{i\alpha k3} + C_{i3k\alpha}) \xi_{\alpha} \tilde{u}_{k,3} - C_{i3k3} \tilde{u}_{k,33} = 0 \quad (B.2)$$

where $\alpha, \beta = 1, 2$. Writing $\xi = [\xi_1, \xi_2, 0] = \eta\mathbf{n}$, a general solution of Eq. B.2 can be expressed as

$$\tilde{u}_k(\xi_1, \xi_2, x_3) = a \exp(-ip\eta x_3) \quad (B.3)$$
The tensor $\mathbf{E}$ is given as (Lazar, 2016)

$$E_{ijkm} = \frac{1}{4\pi} \int_{0}^{\pi} \int_{0}^{2\pi} n_m n_k (nCn)_{ij}^{-1} \sin\theta \, d\theta \, d\phi \quad (A.4)$$

where $\mathbf{n}$ is here the external unit normal of the unit sphere, i.e. $n = [\sin\theta \cos\phi, \sin\theta \sin\phi, \cos\theta]$. The tensor $\mathbf{H}$ is not needed for the present work. Its expression can be found in (Barnett, 1972; Lazar, 2016).

## Appendix B. Imaginary part of the anisotropic bi-material Green tensor

The purpose of this appendix is to explain briefly the general approach of Pan and Yuan (2000) to obtain the anisotropic elastic Green tensor in a bi-material and to provide the equations needed to perform the numerical calculations. All the details can be found in (Pan and Yuan, 2000).

A full-space bi-material is considered where the region $x_3 > 0$ is occupied by material (1) and the region $x_3 < 0$ by material (2). A point force is applied at $(y_1, y_2, y_3)$. To find the displacement vector $\mathbf{u}$ and thus the Green function tensor $\mathbf{G}^{Bi}$, a two-dimensional Fourier transform is applied

$$\tilde{u}_k(\xi_1, \xi_2, x_3) = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} u_k(x_1, x_2, x_3) e^{i(x_1\xi_1 + x_2\xi_2)} dx_1 dx_2 \quad (B.1)$$

In the Fourier transformed domain and in the absence of body force, the equilibrium equation then writes in each material

$$C_{i\alpha k\beta} \xi_{\alpha} \xi_{\beta} \tilde{u}_k + i(C_{i\alpha k3} + C_{i3k\alpha}) \xi_{\alpha} \tilde{u}_{k,3} - C_{i3k3} \tilde{u}_{k,33} = 0 \quad (B.2)$$

where $\alpha, \beta = 1, 2$. Writing $\xi = [\xi_1, \xi_2, 0] = \eta\mathbf{n}$, a general solution of Eq. B.2 can be expressed as

$$\tilde{u}_k(\xi_1, \xi_2, x_3) = a \exp(-ip\eta x_3) \quad (B.3)$$
where $p$ and $a$ satisfy the Stroh eigenrelation (Stroh, 1958, 1962) for the oblique plane spanned by $n = [\cos\theta, \sin\theta, 0]^T$ and $m = [0, 0, 1]^T$

$$[Q + p(R + R^T) + p^2T]a = 0 \qquad (B.4)$$

with

$$Q_{ik} = C_{ijkl}n_j n_s, \quad R_{ik} = C_{ijkl}n_j m_s, \quad T_{ik} = C_{ijkl}m_j m_s \qquad (B.5)$$

Denoting $p_i$ and $a_i$ ($i = 1, 2, ..., 6$) the eigenvalues and the associated eigenvectors, respectively, we let

$$\mathrm{Im}(p_i) > 0, \quad p_{i+3} = \bar{p}_i, \quad a_{i+3} = \bar{a}_i \quad (i = 1, 2, 3) \qquad (B.6)$$
$$A = [a_1, a_2, a_3], \quad B = [b_1, b_2, b_3] \quad \text{with} \quad b_i = (R^T + p_i T)a_i$$

where there is no summation over the repeated index $i$ in Eq. B.6. From a practical perspective, the eigenvalues and the eigenvectors are actually computed from the following standard eigenrelation which is a rewriting of Eqs. B.4 and B.6 (Ting, 1992)

$$\begin{bmatrix} N_1 & N_2 \\ N_3 & N_1^T \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = p \begin{bmatrix} a \\ b \end{bmatrix} \qquad (B.7)$$

$$N_1 = -T^{-1}R^T, \quad N_2 = -T^{-1}, \quad N_3 = RT^{-1}R^T - Q \qquad (B.8)$$

Moreover, the eigenvectors are normalized according to the relation (Pan and Yuan, 2000)

$$b_i^T a_j + a_i^T b_j = \delta_{ij} \qquad (B.9)$$

Finally, a useful quantity that will appear in the analysis is the tensor $M$ defined by
where $p$ and $\mathbf{a}$ satisfy the Stroh eigenrelation (Stroh, 1958, 1962) for the oblique plane spanned by $\mathbf{n} = [\cos\theta, \sin\theta, 0]^T$ and $\mathbf{m} = [0, 0, 1]^T$

$$[\mathbf{Q} + p(\mathbf{R} + \mathbf{R}^T) + p^2\mathbf{T}] \mathbf{a} = \mathbf{0} \qquad (B.4)$$

with

$$Q_{ik} = C_{ijkl}n_j n_s, \quad R_{ik} = C_{ijkl}n_j m_s, \quad T_{ik} = C_{ijkl}m_j m_s \qquad (B.5)$$

Denoting $p_i$ and $\mathbf{a}_i$ ($i = 1, 2, ..., 6$) the eigenvalues and the associated eigenvectors, respectively, we let

$$\mathrm{Im}(p_i) > 0, \quad p_{i+3} = \bar{p}_i, \quad \mathbf{a}_{i+3} = \bar{\mathbf{a}}_i \quad (i = 1, 2, 3) \qquad (B.6) \\ \mathbf{A} = [\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3], \quad \mathbf{B} = [\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3] \quad \text{with} \quad \mathbf{b}_i = (\mathbf{R}^T + p_i \mathbf{T}) \mathbf{a}_i$$

where there is no summation over the repeated index $i$ in Eq. B.6. From a practical perspective, the eigenvalues and the eigenvectors are actually computed from the following standard eigenrelation which is a rewriting of Eqs. B.4 and B.6 (Ting, 1992)

$$\begin{bmatrix} \mathbf{N}_1 & \mathbf{N}_2 \\ \mathbf{N}_3 & \mathbf{N}_1^T \end{bmatrix} \begin{bmatrix} \mathbf{a} \\ \mathbf{b} \end{bmatrix} = p \begin{bmatrix} \mathbf{a} \\ \mathbf{b} \end{bmatrix} \qquad (B.7)$$

$$\mathbf{N}_1 = -\mathbf{T}^{-1}\mathbf{R}^T, \quad \mathbf{N}_2 = -\mathbf{T}^{-1}, \quad \mathbf{N}_3 = \mathbf{R}\mathbf{T}^{-1}\mathbf{R}^T - \mathbf{Q} \qquad (B.8)$$

Moreover, the eigenvectors are normalized according to the relation (Pan and Yuan, 2000)

$$\mathbf{b}_i^T \mathbf{a}_j + \mathbf{a}_i^T \mathbf{b}_j = \delta_{ij} \qquad (B.9)$$

Finally, a useful quantity that will appear in the analysis is the tensor $\mathbf{M}$ defined by
$$ M = -iBA^{-1} \qquad (B.10) $$

The general solution in the transformed domain is obtained as the superposition of the six eigensolutions. The unknown quantities are determined from continuity relations at the plane $x_3 = y_3$ where the point force is applied and from the continuity of the displacement and traction vectors at the interface $x_3 = 0$ (Pan and Yuan, 2000). Then, the solution in the physical domain, $G^{Bi}$, is derived by applying the inverse Fourier transform in the polar coordinate (Pan and Yuan, 2000). In the present work, $G^{Bi}$ is written as the sum of a full-space term and an imaginary term $G^{Im}$. The full-space homogeneous Green tensor is explicitly given by Eq. A.1 but it is noteworthy that it can also be expressed via the present Stroh formalism. The numerical computation is however a little bit longer due to the fact that an eigenproblem must be solved for each value of $\theta$ used in the line integral numerical approximation. This observation remains true for the computation of $G^{Bi}$.

The final expressions of $G^{Im}$ are finally given below. If $x_3 > 0$ and $y_3 \ge 0$, we have (Pan and Yuan, 2000; Chu et al., 2012)

$$ G^{Im}(x, y) = \frac{1}{2\pi^2} \int_0^\pi \overline{A_1} G^{(11)}_u A_1^T d\theta \qquad (B.11) $$

with

$$ (G_u^{(11)})_{ij} = \frac{(G_{11})_{ij}}{-\overline{p_i}^{(1)}x_3 + p_j^{(1)}y_3 - (x_1 - y_1)\cos\theta - (x_2 - y_2)\sin\theta} \qquad (B.12) $$

$$ G_{11} = -\overline{A_1}^{-1} (\overline{M_1} + M_2)^{-1} (M_1 - M_2) A_1 \qquad (B.13) $$

where $A_1$, $M_1$ and $p_i^{(1)}$ are obtained by solving the Stroh eigenrelation (Eq. B.7) in material (1), while $A_2$, $M_2$ and $p_i^{(2)}$ are obtained by solving the Stroh eigenrelation (Eq. B.7) in material (2).
$$ M = -iBA^{-1} \qquad (B.10) $$

The general solution in the transformed domain is obtained as the superposition of the six eigensolutions. The unknown quantities are determined from continuity relations at the plane $x_3 = y_3$ where the point force is applied and from the continuity of the displacement and traction vectors at the interface $x_3 = 0$ (Pan and Yuan, 2000). Then, the solution in the physical domain, $G^{Bi}$, is derived by applying the inverse Fourier transform in the polar coordinate (Pan and Yuan, 2000). In the present work, $G^{Bi}$ is written as the sum of a full-space term and an imaginary term $G^{Im}$. The full-space homogeneous Green tensor is explicitly given by Eq. A.1 but it is noteworthy that it can also be expressed via the present Stroh formalism. The numerical computation is however a little bit longer due to the fact that an eigenproblem must be solved for each value of $\theta$ used in the line integral numerical approximation. This observation remains true for the computation of $G^{Bi}$.

The final expressions of $G^{Im}$ are finally given below. If $x_3 > 0$ and $y_3 \ge 0$, we have (Pan and Yuan, 2000; Chu et al., 2012)

$$ G^{Im}(x, y) = \frac{1}{2\pi^2} \int_0^\pi \overline{A_1} G^{(11)}_u A_1^T d\theta \qquad (B.11) $$

with

$$ (G_u^{(11)})_{ij} = \frac{(G_{11})_{ij}}{-\overline{p_i}^{(1)}x_3 + p_j^{(1)}y_3 - (x_1 - y_1)\cos\theta - (x_2 - y_2)\sin\theta} \qquad (B.12) $$

$$ G_{11} = -\overline{A_1}^{-1} (\overline{M_1} + M_2)^{-1} (M_1 - M_2) A_1 \qquad (B.13) $$

where $A_1$, $M_1$ and $p_i^{(1)}$ are obtained by solving the Stroh eigenrelation (Eq. B.7) in material (1), while $A_2$, $M_2$ and $p_i^{(2)}$ are obtained by solving the Stroh eigenrelation (Eq. B.7) in material (2).
If $x_3 < 0$ and $y_3 > 0$, we have (Pan and Yuan, 2000; Chu et al., 2012)

$$G^{Im}(x, y) = -\frac{1}{2\pi^2} \int_{0}^{\pi} A_2 G_u^{(21)} A_1^T d\theta - G^{\infty 2}(x, y) \quad (B.14)$$

$$(G_u^{(21)})_{ij} = \frac{(G_{21})_{ij}}{-\overline{p_i}^{(2)}x_3 + p_j^{(1)}y_3 - (x_1 - y_1)\cos\theta - (x_2 - y_2)\sin\theta} \quad (B.15)$$

$$G_{21} = A_2^{-1} (\overline{M_1} + M_2)^{-1} (M_1 + \overline{M_1}) A_1 \quad (B.16)$$

If $x_3 < 0$ and $y_3 \le 0$, we have (Pan and Yuan, 2000; Chu et al., 2012)

$$G^{Im}(x, y) = -\frac{1}{2\pi^2} \int_{0}^{\pi} A_2 G_u^{(22)} \overline{A_2}^T d\theta \quad (B.17)$$

$$(G_u^{(22)})_{ij} = \frac{(G_{22})_{ij}}{-p_i^{(2)}x_3 + \overline{p_j}^{(1)}y_3 - (x_1 - y_1)\cos\theta - (x_2 - y_2)\sin\theta} \quad (B.18)$$

$$G_{22} = A_2^{-1} (\overline{M_1} + M_2)^{-1} (\overline{M_1} - \overline{M_2}) \overline{A_2} \quad (B.19)$$

If $x_3 > 0$ and $y_3 < 0$, we have (Pan and Yuan, 2000; Chu et al., 2012)

$$G^{Im}(x, y) = \frac{1}{2\pi^2} \int_{0}^{\pi} \overline{A_1} G_u^{(12)} \overline{A_2}^T d\theta - G^{\infty 1}(x, y) \quad (B.20)$$

$$(G_u^{(12)})_{ij} = \frac{(G_{12})_{ij}}{-\overline{p_i}^{(1)}x_3 + \overline{p_j}^{(2)}y_3 - (x_1 - y_1)\cos\theta - (x_2 - y_2)\sin\theta} \quad (B.21)$$

$$G_{12} = \overline{A_1}^{-1} (\overline{M_1} + M_2)^{-1} (\overline{M_2} + M_2) \overline{A_2} \quad (B.22)$$
If $x_3 < 0$ and $y_3 > 0$, we have (Pan and Yuan, 2000; Chu et al., 2012)

$$G^{Im}(x, y) = -\frac{1}{2\pi^2} \int_{0}^{\pi} A_2 G_u^{(21)} A_1^T d\theta - G^{\infty 2}(x, y) \quad (B.14)$$

$$(G_u^{(21)})_{ij} = \frac{(G_{21})_{ij}}{-\overline{p_i}^{(2)}x_3 + p_j^{(1)}y_3 - (x_1 - y_1)\cos\theta - (x_2 - y_2)\sin\theta} \quad (B.15)$$

$$G_{21} = A_2^{-1} (\overline{M_1} + M_2)^{-1} (M_1 + \overline{M_1}) A_1 \quad (B.16)$$

If $x_3 < 0$ and $y_3 \le 0$, we have (Pan and Yuan, 2000; Chu et al., 2012)

$$G^{Im}(x, y) = -\frac{1}{2\pi^2} \int_{0}^{\pi} A_2 G_u^{(22)} \overline{A_2}^T d\theta \quad (B.17)$$

$$(G_u^{(22)})_{ij} = \frac{(G_{22})_{ij}}{-p_i^{(2)}x_3 + \overline{p_j}^{(1)}y_3 - (x_1 - y_1)\cos\theta - (x_2 - y_2)\sin\theta} \quad (B.18)$$

$$G_{22} = A_2^{-1} (\overline{M_1} + M_2)^{-1} (\overline{M_1} - \overline{M_2}) \overline{A_2} \quad (B.19)$$

If $x_3 > 0$ and $y_3 < 0$, we have (Pan and Yuan, 2000; Chu et al., 2012)

$$G^{Im}(x, y) = \frac{1}{2\pi^2} \int_{0}^{\pi} \overline{A_1} G_u^{(12)} \overline{A_2}^T d\theta - G^{\infty 1}(x, y) \quad (B.20)$$

$$(G_u^{(12)})_{ij} = \frac{(G_{12})_{ij}}{-\overline{p_i}^{(1)}x_3 + \overline{p_j}^{(2)}y_3 - (x_1 - y_1)\cos\theta - (x_2 - y_2)\sin\theta} \quad (B.21)$$

$$G_{12} = \overline{A_1}^{-1} (\overline{M_1} + M_2)^{-1} (\overline{M_2} + M_2) \overline{A_2} \quad (B.22)$$
Appendix B.1. First derivatives of the imaginary part of the anisotropic bi-material Green tensor

The first derivatives of $G^{Im}$ are obtained by straightforward derivations.
Therefore, they are written explicitly only for the case $x_3 > 0$ and $y_3 \ge 0$.

$$ \frac{\partial G_{ij}^{Im}}{\partial x_k}(\mathbf{x}, \mathbf{y}) = \frac{1}{2\pi^2} \int_0^\pi \overline{(A_1)_{im}} \frac{\left(\frac{\partial G_u^{(11)}}{\partial x_k}\right)_{mn}}{\left(A_1^T\right)_{nj}} (A_1^T)_{nj} d\theta \quad (B.23) $$

$$ \frac{\partial \left( G_u^{(11)} \right)_{ij}}{\partial x_1} = \frac{\cos \theta (G_{11})_{ij}}{\left( -\overline{p_i}^{(1)} x_3 + p_j^{(1)} y_3 - (x_1 - y_1) \cos \theta - (x_2 - y_2) \sin \theta \right)^2} \quad (B.24) $$

$$ \frac{\partial \left( G_u^{(11)} \right)_{ij}}{\partial x_2} = \frac{\sin \theta (G_{11})_{ij}}{\left( -\overline{p_i}^{(1)} x_3 + p_j^{(1)} y_3 - (x_1 - y_1) \cos \theta - (x_2 - y_2) \sin \theta \right)^2} \quad (B.25) $$

$$ \frac{\partial \left( G_u^{(11)} \right)_{ij}}{\partial x_3} = \frac{\overline{p_i}^{(1)} (G_{11})_{ij}}{\left( -\overline{p_i}^{(1)} x_3 + p_j^{(1)} y_3 - (x_1 - y_1) \cos \theta - (x_2 - y_2) \sin \theta \right)^2} \quad (B.26) $$
Appendix B.1. First derivatives of the imaginary part of the anisotropic bi-material Green tensor

The first derivatives of $G^{Im}$ are obtained by straightforward derivations.
Therefore, they are written explicitly only for the case $x_3 > 0$ and $y_3 \ge 0$.

$$ \frac{\partial G_{ij}^{Im}}{\partial x_k}(\mathbf{x}, \mathbf{y}) = \frac{1}{2\pi^2} \int_0^\pi \overline{(A_1)_{im}} \frac{\left(\frac{\partial G_u^{(11)}}{\partial x_k}\right)_{mn}}{(A_1^T)_{nj}} (A_1^T)_{nj} d\theta \quad (B.23) $$

$$ \frac{\partial \left( G_u^{(11)} \right)_{ij}}{\partial x_1} = \frac{\cos \theta (G_{11})_{ij}}{\left( -\overline{p_i}^{(1)} x_3 + p_j^{(1)} y_3 - (x_1 - y_1) \cos \theta - (x_2 - y_2) \sin \theta \right)^2} \quad (B.24) $$

$$ \frac{\partial \left( G_u^{(11)} \right)_{ij}}{\partial x_2} = \frac{\sin \theta (G_{11})_{ij}}{\left( -\overline{p_i}^{(1)} x_3 + p_j^{(1)} y_3 - (x_1 - y_1) \cos \theta - (x_2 - y_2) \sin \theta \right)^2} \quad (B.25) $$

$$ \frac{\partial \left( G_u^{(11)} \right)_{ij}}{\partial x_3} = \frac{\overline{p_i}^{(1)} (G_{11})_{ij}}{\left( -\overline{p_i}^{(1)} x_3 + p_j^{(1)} y_3 - (x_1 - y_1) \cos \theta - (x_2 - y_2) \sin \theta \right)^2} \quad (B.26) $$
References

Abramowitz, M., Stegun, I.A., 1968. Handbook of mathematical functions with formulas, graphs, and mathematical tables. volume 55. US Government printing office.

Aspinall, J., Armstrong, D.E., Pasta, M., 2022. Ebsd-coupled indentation: nanoscale mechanics of lithium metal. Materials Today Energy 30, 101183.

Bacon, D., 1972. The mechanical interaction of a point defect with a surface layer. physica status solidi (b) 50, 607-617.

Bacon, D., Barnett, D., Scattergood, R.O., 1980. Anisotropic continuum theory of lattice defects. Progress in Materials Science 23, 51-262.

Balluffi, R.W., 2016. Introduction to elasticity theory for crystal defects. World Scientific Publishing Company.

Barnett, D., 1972. The precise evaluation of derivatives of the anisotropic elastic green's functions. physica status solidi (b) 49, 741-748.

Barnett, D., Cai, W., 2018. Properties of the eshelby tensor and existence of the equivalent ellipsoidal inclusion solution. Journal of the Mechanics and Physics of Solids 121, 71-80.

Barnett, D., Lothe, J., 1975. Line force loadings on anisotropic half-spaces and wedges. Physica Norvegica 8, 13-22.

Berbenni, S., Berveiller, M., Richeton, T., 2008. Intra-granular plastic slip heterogeneities: Discrete vs. Mean Field approaches. International Journal of Solids and Structures 45, 4147-4172.

Cai, W., Sills, R., Nix, W., 2014. Modeling a distribution of point defects as misfitting inclusions in stresses solids. Journal of the Mechanics and Physics of Solids 66, 154-171.
References

Abramowitz, M., Stegun, I.A., 1968. Handbook of mathematical functions with formulas, graphs, and mathematical tables. volume 55. US Government printing office.

Aspinall, J., Armstrong, D.E., Pasta, M., 2022. Ebsd-coupled indentation: nanoscale mechanics of lithium metal. Materials Today Energy 30, 101183.

Bacon, D., 1972. The mechanical interaction of a point defect with a surface layer. physica status solidi (b) 50, 607-617.

Bacon, D., Barnett, D., Scattergood, R.O., 1980. Anisotropic continuum theory of lattice defects. Progress in Materials Science 23, 51-262.

Balluffi, R.W., 2016. Introduction to elasticity theory for crystal defects. World Scientific Publishing Company.

Barnett, D., 1972. The precise evaluation of derivatives of the anisotropic elastic green's functions. physica status solidi (b) 49, 741-748.

Barnett, D., Cai, W., 2018. Properties of the eshelby tensor and existence of the equivalent ellipsoidal inclusion solution. Journal of the Mechanics and Physics of Solids 121, 71-80.

Barnett, D., Lothe, J., 1975. Line force loadings on anisotropic half-spaces and wedges. Physica Norvegica 8, 13-22.

Berbenni, S., Berveiller, M., Richeton, T., 2008. Intra-granular plastic slip heterogeneities: Discrete vs. Mean Field approaches. International Journal of Solids and Structures 45, 4147-4172.

Cai, W., Sills, R., Nix, W., 2014. Modeling a distribution of point defects as misfitting inclusions in stresses solids. Journal of the Mechanics and Physics of Solids 66, 154-171.
Chen, X., Richeton, T., C.Motz, S.Berbenni, 2019. 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, 141–156.

Chen, X., Richeton, T., C.Motz, S.Berbenni, 2020. Atomic force microscopy study of discrete dislocation pile-ups at grain boundaries in bi-crystalline micro-pillars. Crystals 10, 411.

Chen, X., Richeton, T., C.Motz, S.Berbenni, 2021. Surface effects on image stresses and dislocation pile-ups in anisotropic bi-crystals. International Journal of Plasticity 143, 102967.

Chiu, Y., 1978. On the stress field and surface deformation in a half space with a cuboidal zone in which initial strains are uniform. Journal of Applied Mechanics, Transactions ASME 44, 587–590.

Chu, H., Pan, E., Han, X., Wang, J., Beyerlein, I., 2012. Elastic fields of dislocation loops in three-dimensional anisotropic bimaterials. Journal of the Mechanics and Physics of Solids 60, 418–431.

Chu, H., Pan, E., Wang, J., Beyerlein, I., 2011. Three-dimensional elastic displacements induced by a dislocation of polygonal shape in anisotropic elastic crystals. International Journal of Solids and Structures 48, 1164–1170.

Clouet, E., Varenne, C., Jourdan, T., 2018. Elastic modeling of point-defects and their interaction. Computational Materials Science 147, 49–63.

Cochardt, A., Schoek, G., Wiedersich, H., 1955. Interaction between dislocations and interstitial atoms in body-centered cubic metals. Acta metallurgica 3, 533–537.

Dingreville, R., Berbenni, S., 2016. On the interaction of solutes with grain boundaries. Acta Materialia 104, 237–249.
Chen, X., Richeton, T., C.Motz, S.Berbenni, 2019. 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, 141–156.

Chen, X., Richeton, T., C.Motz, S.Berbenni, 2020. Atomic force microscopy study of discrete dislocation pile-ups at grain boundaries in bi-crystalline micro-pillars. Crystals 10, 411.

Chen, X., Richeton, T., C.Motz, S.Berbenni, 2021. Surface effects on image stresses and dislocation pile-ups in anisotropic bi-crystals. International Journal of Plasticity 143, 102967.

Chiu, Y., 1978. On the stress field and surface deformation in a half space with a cuboidal zone in which initial strains are uniform. Journal of Applied Mechanics, Transactions ASME 44, 587–590.

Chu, H., Pan, E., Han, X., Wang, J., Beyerlein, I., 2012. Elastic fields of dislocation loops in three-dimensional anisotropic bimaterials. Journal of the Mechanics and Physics of Solids 60, 418–431.

Chu, H., Pan, E., Wang, J., Beyerlein, I., 2011. Three-dimensional elastic displacements induced by a dislocation of polygonal shape in anisotropic elastic crystals. International Journal of Solids and Structures 48, 1164–1170.

Clouet, E., Varenne, C., Jourdan, T., 2018. Elastic modeling of point-defects and their interaction. Computational Materials Science 147, 49–63.

Cochardt, A., Schoek, G., Wiedersich, H., 1955. Interaction between dislocations and interstitial atoms in body-centered cubic metals. Acta metallurgica 3, 533–537.

Dingreville, R., Berbenni, S., 2016. On the interaction of solutes with grain boundaries. Acta Materialia 104, 237–249.
Eshelby, J., 1954. Distortion of a crystal by point imperfections. Journal of Applied Physics 25, 255-261.

Eshelby, J., 1957. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proceedings of the Royal Society London A 241, 376-396.

Eshelby, J., 1961. Elastic inclusions and inhomogenieties. Progr. Solid Mech. 2, 89.

Franciosi, P., Charles, Y., 2016. Mean green operators and eshelby tensors for hemispherical inclusions and hemisphere interactions in spheres. application to bi-material spherical inclusions in isotropic spaces. Mechanics Research Communications 75, 57-66.

Gladwell, G., 1999. On inclusions at a bi-material elastic interface. Journal of elasticity 54, 27-41.

Hill, R., 1965. Continuum micro-mechanics of elastoplastic polycrystals. Journal of the Mechanics and Physics of Solids 13, 89-101.

Hirth, J., Lothe, J., 1982. Theory of dislocations. 2nd ed ed., Wiley, New York.

Jiang, X., Pan, E., 2004. Exact solution for 2d polygonal inclusion problem in anisotropic magnetoelectroelastic full-, half-, and bimaterial-planes. International Journal of Solids and Structures 41, 4361-4382.

Kinoshita, N., Mura, T., 1984. Eigenstrain problems in a finite elastic body. SIAM Journal on Applied Mathematics 44, 524-535.

Kneer, G., 1965. Über die berechnung der elastizitätsmoduln vielkristalliner aggregate mit textur. Physica Status Solidi 9, 825-838.
Eshelby, J., 1954. Distortion of a crystal by point imperfections. Journal of Applied Physics 25, 255-261.

Eshelby, J., 1957. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proceedings of the Royal Society London A 241, 376-396.

Eshelby, J., 1961. Elastic inclusions and inhomogenieties. Progr. Solid Mech. 2, 89.

Franciosi, P., Charles, Y., 2016. Mean green operators and eshelby tensors for hemispherical inclusions and hemisphere interactions in spheres. application to bi-material spherical inclusions in isotropic spaces. Mechanics Research Communications 75, 57-66.

Gladwell, G., 1999. On inclusions at a bi-material elastic interface. Journal of elasticity 54, 27-41.

Hill, R., 1965. Continuum micro-mechanics of elastoplastic polycrystals. Journal of the Mechanics and Physics of Solids 13, 89-101.

Hirth, J., Lothe, J., 1982. Theory of dislocations. 2nd ed ed., Wiley, New York.

Jiang, X., Pan, E., 2004. Exact solution for 2d polygonal inclusion problem in anisotropic magnetoelectroelastic full-, half-, and bimaterial-planes. International Journal of Solids and Structures 41, 4361-4382.

Kinoshita, N., Mura, T., 1984. Eigenstrain problems in a finite elastic body. SIAM Journal on Applied Mathematics 44, 524-535.

Kneer, G., 1965. Über die berechnung der elastizitätsmoduln vielkristalliner aggregate mit textur. Physica Status Solidi 9, 825-838.
Kröner, E., 1990. Modified green functions in the theory of heterogeneous and/or anisotropic linearly elastic media, in: Micromechanics and inhomogeneity: The Toshio Mura 65th Anniversary. G.J. Weng et al., Springer, Berlin, pp. 197–211.

Lazar, M., 2016. Micromechanics and dislocation theory in anisotropic elasticity. Journal of Micromechanics and Molecular Physics 01, 1650011.

Lazar, M., 2017. Micromechanics and theory of point defects in anisotropic elasticity: Eshelby factor meets eshelby tensor. Journal of Micromechanics and Molecular Physics 02, 1750005.

Li, S., Sauer, R.A., Wang, G., 2007. The eshelby tensors in a finite spherical domain—part i: theoretical formulations. Journal of Applied Mechanics 74, 770–783.

Mindlin, R.D., Cheng, D.H., 1950. Thermoelastic stress in the semi-infinite solid. Journal of Applied Physics 21, 931–933.

Mura, T., 1987. Micromechanics of Defects in Solids. Martinus Nijhoff, Dordrecht.

Nye, J., 1957. Physical Properties of Crystals. London: Oxford University Press.

Pan, E., Yang, B., 2001. Elastostatic fields in anisotropic substrate due to a buried quatum dot. Journal of Applied Physics 90, 6190–6196.

Pan, E., Yuan, F., 2000. Three-dimensional green's functions in anisotropic bimaterials. International Journal of Solids and Structures 37, 5329–5351.

Richeton, T., 2019. Mechanical derivation of elastic constants of $\varepsilon$-martensite. Scripta Materialia 169, 14–18.
Kröner, E., 1990. Modified green functions in the theory of heterogeneous and/or anisotropic linearly elastic media, in: Micromechanics and inhomogeneity: The Toshio Mura 65th Anniversary. G.J. Weng et al., Springer, Berlin, pp. 197–211.

Lazar, M., 2016. Micromechanics and dislocation theory in anisotropic elasticity. Journal of Micromechanics and Molecular Physics 01, 1650011.

Lazar, M., 2017. Micromechanics and theory of point defects in anisotropic elasticity: Eshelby factor meets eshelby tensor. Journal of Micromechanics and Molecular Physics 02, 1750005.

Li, S., Sauer, R.A., Wang, G., 2007. The eshelby tensors in a finite spherical domain—part i: theoretical formulations. Journal of Applied Mechanics 74, 770–783.

Mindlin, R.D., Cheng, D.H., 1950. Thermoelastic stress in the semi-infinite solid. Journal of Applied Physics 21, 931–933.

Mura, T., 1987. Micromechanics of Defects in Solids. Martinus Nijhoff, Dordrecht.

Nye, J., 1957. Physical Properties of Crystals. London: Oxford University Press.

Pan, E., Yang, B., 2001. Elastostatic fields in anisotropic substrate due to a buried quatum dot. Journal of Applied Physics 90, 6190–6196.

Pan, E., Yuan, F., 2000. Three-dimensional green's functions in anisotropic bimaterials. International Journal of Solids and Structures 37, 5329–5351.

Richeton, T., 2019. Mechanical derivation of elastic constants of $\varepsilon$-martensite. Scripta Materialia 169, 14–18.
Richeton, T., Chen, X., Berbenni, S., 2020. Misorientation dependence of the grain boundary migration rate: role of elastic anisotropy. Philosophical Magazine 100, 2436-2457.

Rios-Jara, D., Planes, A., Mañosa, L., Ortin, J., Belkahla, S., Morin, M., Guénin, G., Macqueron, J., 1991. Martensitic transition entropy change and elastic constants of cu-al-be alloys. Le Journal de Physique IV 1, C4-283.

Ru, C., 2001. A two-dimensional eshelby problem for two bonded piezoelectric half-planes. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 457, 865-883.

Ru, C., Schiavone, P., Mioduchowski, A., 2001. Elastic fields in two jointed half-planes with an inclusion of arbitrary shape. Zeitschrift für angewandte Mathematik und Physik ZAMP 52, 18-32.

Selvadurai, A., 2000. An inclusion at a bi-material elastic interface. Journal of Engineering Mathematics 37, 155-170.

Seo, K., Mura, T., 1979. The elastic field in a half space due to ellipsoidal inclusions with uniform dilatational eigenstrains. Journal of Applied Mechanics 46, 568-572.

Siems, R., 1968. Mechanical interactions of point defects. physical status solidi (b) 30, 645-658.

Stroh, A., 1958. Dislocations and cracks in anisotropic elasticity. Philosophical Magazine 3, 625-646.

Stroh, A., 1962. Steady state problems in anisotropic elasticity. Journal of Mathematical Physics 41, 77-103.
Richeton, T., Chen, X., Berbenni, S., 2020. Misorientation dependence of the grain boundary migration rate: role of elastic anisotropy. Philosophical Magazine 100, 2436-2457.

Rios-Jara, D., Planes, A., Mañosa, L., Ortin, J., Belkahla, S., Morin, M., Guénin, G., Macqueron, J., 1991. Martensitic transition entropy change and elastic constants of cu-al-be alloys. Le Journal de Physique IV 1, C4-283.

Ru, C., 2001. A two-dimensional eshelby problem for two bonded piezoelectric half-planes. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 457, 865-883.

Ru, C., Schiavone, P., Mioduchowski, A., 2001. Elastic fields in two jointed half-planes with an inclusion of arbitrary shape. Zeitschrift für angewandte Mathematik und Physik ZAMP 52, 18-32.

Selvadurai, A., 2000. An inclusion at a bi-material elastic interface. Journal of Engineering Mathematics 37, 155-170.

Seo, K., Mura, T., 1979. The elastic field in a half space due to ellipsoidal inclusions with uniform dilatational eigenstrains. Journal of Applied Mechanics 46, 568-572.

Siems, R., 1968. Mechanical interactions of point defects. physical status solidi (b) 30, 645-658.

Stroh, A., 1958. Dislocations and cracks in anisotropic elasticity. Philosophical Magazine 3, 625-646.

Stroh, A., 1962. Steady state problems in anisotropic elasticity. Journal of Mathematical Physics 41, 77-103.
Ting, T., 1992. Image singularities of green's functions for anisotropic elastic half-spaces and bimaterials. The Quarterly Journal of Mechanics and Applied Mathematics 45, 119-139.

Wang, X., Sudak, L., Ru, C., 2007. Elastic fields in two imperfectly bonded half-planes with a thermal inclusion of arbitrary shape. Zeitschrift für angewandte Mathematik und Physik 58, 488-509.

White, C., Coghlan, W., 1977. The spectrum of binding energies approach to grain boundary segregation. Metallurgical Transactions A 8A, 1403-1412.

von Winckel, G., 2004. Legendre-gauss quadrature weights and nodes. URL: https://www.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes.

Zhang, H., Chou, Y., 1985. Antiplane eigenstrain problem of an elliptic inclusion in a two-phase anisotropic medium. Journal of Applied Mechanics, Transactions ASME 52, 87-90.

Zhou, K., 2013. A review of recent works on inclusions. Mechanics of Materials 60, 144-158.
Ting, T., 1992. Image singularities of green's functions for anisotropic elastic half-spaces and bimaterials. The Quarterly Journal of Mechanics and Applied Mathematics 45, 119-139.

Wang, X., Sudak, L., Ru, C., 2007. Elastic fields in two imperfectly bonded
half-planes with a thermal inclusion of arbitrary shape. Zeitschrift für
angewandte Mathematik und Physik 58, 488–509.

White, C., Coghlan, W., 1977. The spectrum of binding energies approach to
grain boundary segregation. Metallurgical Transactions A 8A, 1403–1412.

von Winckel, G., 2004. Legendre-gauss quadrature weights and nodes.
URL: https://www.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes.

Zhang, H., Chou, Y., 1985. Antiplane eigenstrain problem of an elliptic inclusion in a two-phase anisotropic medium. Journal of Applied Mechanics, Transactions ASME 52, 87-90.

Zhou, K., 2013. A review of recent works on inclusions. Mechanics of Materials 60, 144–158.