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

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 dilatational 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


Introduction
The study of elastic fields generated by inclusions within an elastic matrix 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, representing 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 problems within finite spaces (e.g., (Eshelby, 1954(Eshelby, , 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)), necessitating 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 conditions 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 δ-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 nature, 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 formalism (Stroh, 1958(Stroh, , 1962) ) and the 2D Fourier transform method.They encompass 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., 2011Chu et al., , 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 materials (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 regarding 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 comparisons 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 matrices.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.

Theoretical framework
We consider an infinite general (i.e.homogeneous or heterogeneous) elastic medium of volume V and a distribution of eigenstrains ε * .In the absence of body force, the displacement vector u follows the Navier equation 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 ε * belongs only to material (2).On the right, the inclusion spans the two materials.
where C is the fourth-rank tensor of anisotropic elastic constants that possesses the symmetries x denotes a particular point of V and the indices refer to a rectangular Cartesian frame (e 1 , e 2 , e 3 ).The elastic Green tensor G of the Navier equation is defined by (Barnett, 1972;Mura, 1987;Lazar, 2016Lazar, , 2017) δ im is the Kronecker delta tensor and δ(x − y) is the three-dimensional Dirac delta function.The Green's function G ij (x, y) represents the elastic displacement in the e i direction at point x induced by a unit force in the e j direction applied at point y (Mura, 1987;Barnett, 1972).For a distribution of eigenstrains, it satisfies the equation (Mura, 1987) where the term −C jlmn (y) ∂ε * mn (y) ∂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) Hence, from Hooke's law, the stress field due to a distribution of eigenstrains in a general infinite elastic medium may be expressed as

Homogeneous medium
For an infinite homogeneous body, we introduce G ∞ as the full-space anisotropic elastic Green tensor.We have then , 1987), and thus the more classical expression of the stress field for a distribution of eigenstrains is retrieved

Bi-material
We know 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 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, y), as well as the symmetry relations, G Bi ij (x, y) ̸ = G Bi ji (x, 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 The expressions of G Im depend on the signs of x 3 and y 3 (Pan and Yuan, 2000;Chu et al., 2011Chu et al., , 2012)).They are explicitly given in the Appendix B.

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 ε * 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 where the fourth-rank tensor E is defined as (Kröner, 1990;Lazar, 2016) E is related to the Dirac δ(x−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 (11) ∂V I denotes the surface of the spherical inclusion and n its external unit normal.It must be underlined that the introduction of the tensor 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 (x), although being non-uniform.The latter should connect the total strain with the eigenstrain, ε im (x) = S Esh,Bi imln (x)ε * ln and should possess the symmetries S Esh,Bi imln = S Esh,Bi miln = S Esh,Bi imnl .Hence, we might define 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 Esh,∞ imln x, y) dy with P ∞ the Hill's polarization tensor (Hill, 1965) and Γ ∞ the modified Green function (Kröner, 1990).For a spherical inclusion, we have For the analogy, it is therefore interesting to re-write the Eshelby tensor for a spherical inclusion in a bi-material as whith and It is noteworthy that P ∞ (x) is uniform within a same material whereas P Im (x) is not.

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 ≤ 0 and y 3 + R ≥ 0, the elasticity tensor cannot be taken out of the integral.However, one can write ) 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 where V I( 1) is the volume of the inclusion located in material (1).Using Eq. 8 and Gauss's theorem, we have then It must be underlined that the surface area ∂V I(1) represents the part of ∂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 For x 3 > 0, the same procedure leads to pqim C (1) pqim C (2) jkln E (1)

Cubic elasticity and dilatational eigenstrain
In the case of cubic elasticity and for a dilatational eigenstrain of the form it is noteworthy that C (2) jkln ε * ln = 0 ∀j, ∀k.Hence, Eqs. 17 and 19 reduce seemingly to Eq. 11, namely, for x 3 < 0 and for x 3 > 0 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 σ pq scale directly with δ * in this case.

Elastic strain energy
The elastic strain energy is defined as where ε 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) For a uniform dilatational eigenstrain of the form of Eq.21 (ε * ij = δ * δ ij ), W e scales thus with δ * 2 since σ ij scales also with δ * (see .

Numerical implementation
The numerical integrations over [0, π]  Relative error 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 × N S integration points) and N E for integrals over volumes to get the elastic strain energy (N E × N E × N E integration points).

Convergence study
For the convergence study, a uniform dilatational eigenstrain of the form of Eq.21 with δ * = 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 with the lattice parameter a 0 = 0.361 49 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.1Rand 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 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.5Rand 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 ∂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.

Comparison with known solutions in an isotropic half-space
Seo and Mura (1979) have solved the problem of an ellipsoidal inclusion 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 where ν is the Poisson coefficient and µ 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 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 Hence, the comparison was made using Eqs. 25, 22, 27 and 28 on one side and Eq. 26 with µ = (C 11 − C 12 ) /2 and ν = C 12 / (2 (C 12 + µ)) on the other side.The same values as for the convergence study (Section 3.1) were considered for C 11 , C 12 , δ * 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 σ 33 exhibits a non-monotonous behavior with a maximum value inside the inclusion.

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 S = 30 when the inclusion belongs to a same material and to N S = 100 when the inclusion spans the two materials.Figures 6, 7, 8 and 9 show the variations along e 3 of the stress components which are not negligible for the particular crystallographic 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 .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 (2) jkln ε * ln = 0 ∀j, ∀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.    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  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., σ 11 and σ 22 at y 3 = 1.1R and σ 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 = ±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 = ±2R.
When the inclusion spans the two materials, it can be checked that the tractions continuity is well respected (see σ 33 and σ 23 ).It can be also observed that the stress jumps for σ 11 and σ 22 are more important at y 3 = ±0.5Rthan 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 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 σ 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 = ±2R and y 3 = ±1.1R.In this case too, the stresses are non homogeneous although the distance from the interface does not change.For the particular 'Orientation A' considered, the stress variations are either symmetric for σ 11 , σ 22 , σ 33 and σ 23 or anti-symmetric for σ 31 and σ 12 .As expected, the variations are much more pronounced at y 3 = ±1.1Rthan at y 3 = ±2R.They are also globally weaker than along e 3 , which makes sense since the distance from the interface is fixed.

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.
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 σ 12 and σ 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 σ 31 is well respected while the stress jumps displayed by σ 12 are very high with a jump greater than 2.5 GPa at y 3 = 0.5R.σ 31 exhibits very strong stress gradients which seem to be controlled by the difference between the homogenous reference solutions of materials ( 1) and (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 −θ/2 and θ/2 around e 2 in materials ( 2) and ( 1), respectively, where θ is the tilt or misorientation angle (Hirth and Lothe, 1982).Two cases are computed, θ = 11.5˚andθ = 19.1˚whichcorrespond 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.

Dilatational eigenstrain
For a uniform dilatational eigenstrain of the form of Eq.21, W e scales with δ * 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 Relative energy variation (%) igure 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 (θ being the tilt angle).
decrease of the imaginary term, probably follows the same scaling of R y 3 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 symmetric 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.igure 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 (θ being the tilt angle).

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.

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 16: On the left, maximum energy variation due to image effect with respect to δ * considering a dilatational eigenstrain ε * ij = δ * δ ij and the three orientations of Table 1, as well as the two symmetrical tilt boundaries defined by their tilt angle θ.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, σ ext kk = 1 GPa.On the right is plotted the equivalent value of σ ext kk such that the maximum energy variation due to the image effect equals the interaction energy between the dilatational eigenstrain and σ ext kk .
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 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(Chen et al., , 2020(Chen et al., , 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 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 ε * distributed in V I and an external stress field σ ext (such as the one arising from a grain boundary), as given by (Mura, 1987) 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 ε Figure 16 displays on the left the maximum energy variation per unit volume due to the image effect with respect to δ * 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 δ * 2 , while the interaction energy with a stress field scales linearly with δ * .In the current scenario with σ ext kk = 1 GPa, the image effect becomes predominant for 'Orientations A and B' when δ * > 0.1.To quantitatively evaluate the significance of the image effect, Figure 16 then presents, on the right, the equivalent value of σ ext kk such that the maximum energy variation due to the image effect equals the interaction energy between the dilatational eigenstrain and σ ext kk according to Eq. 31 .These σ ext kk values scale therefore linearly with δ * .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 = 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 ≈ 8 (Aspinall et al., 2022), a crucial component in batteries, or some shape memory alloys like CuAlBe with A ≈ 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 either promote or hinder interface motion.The significance of such an image effect on interface motion was previously studied for distributions of dislocations (Richeton et al., 2020), providing insights into the connection between interface motion and crystallographic disorientations.

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 constants of Cu.Various inclusion positions relative to the interface, different forms of eigenstrain (dilatational (Eq.21) and pure shear (Eq.29)), and different crystallographic orientations are considered.The following key conclusions 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 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 (ε * ij = δ * δ 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 δ * 2 for the image effect and as δ * for the interaction energy with a stress field).
where p and a satisfy the Stroh eigenrelation (Stroh, 1958(Stroh, , 1962) )  Denoting p i and a i (i = 1, 2, ..., 6) the eigenvalues and the associated eigenvectors, respectively, we let Im(p i ) > 0 , p i+3 = p i , a i+3 = a i (i = 1, 2, 3) 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) Moreover, the eigenvectors are normalized according to the relation (Pan and Yuan, 2000) b T i a j + a T i b j = δ ij (B.9) Finally, a useful quantity that will appear in the analysis is the tensor M defined by

Figure 1 :
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 ε * belongs only to material (2).On the right, the inclusion spans the two materials.
involved in the calculation of E ijkm and ∂G Im ij ∂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 ∂V I , ∂V I(1) , ∂V I(2) and V I are obtained from a Gauss-Legendre quadrature method Abramowitz

Figure 2 :
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.

Figure 3 :
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.

Figure 4 :
Figure 4: Comparison of energy variation towards a free surface for a spherical inclusion with a uniform dilatational eigenstrain ε * ij = 0.1δ 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.µ = 23 GPa, ν = 0.4218 and R = 0.1278 nm.

Figure 5 :
Figure 5: Comparison of stress distributions along e 3 (x 1 = x 2 = 0) within a spherical inclusion with a uniform dilatational eigenstrain ε * ij = 0.1δ ij located at y 3 = 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.

Figure 6 :
Figure 6: Variation of σ 11 along e 3 (x 1 = x 2 = 0) for a dilatational eigenstrain ε * ij = 0.1δ 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.

Figure 7 :
Figure 7: Variation of σ 22 along e 3 (x 1 = x 2 = 0) for a dilatational eigenstrain ε * ij = 0.1δ 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.

Figure 8 :
Figure 8: Variation of σ 33 along e 3 (x 1 = x 2 = 0) for a dilatational eigenstrain ε * ij = 0.1δ 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.

Figure 9 :
Figure 9: Variation of σ 23 along e 3 (x 1 = x 2 = 0) for a dilatational eigenstrain ε * ij = 0.1δ 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.

Figure 12 :Figure 13 :
Figure12: Variation of σ 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 media.

Table 1 :
Different crystallographic orientations considered for the two materials.