The Newtonian potential inhomogeneity problem: non-uniform eigenstrains in cylinders of non-elliptical cross section

Understanding the fields that are set up in and around inhomogeneities is of great importance in order to predict the manner in which heterogeneous media behave when subjected to applied loads or other fields, e.g., magnetic, electric, thermal, etc. The classical inhomogeneity problem of an ellipsoid embedded in an unbounded host or matrix medium has long been studied but is perhaps most associated with the name of Eshelby due to his seminal work in 1957, where in the context of the linear elasticity problem, he showed that for imposed far fields that correspond to uniform strains, the strain field induced inside the ellipsoid is also uniform. In Eshelby’s language, this corresponds to requiring a uniform eigenstrain in order to account for the presence of the ellipsoidal inhomogeneity, and the so-called Eshelby tensor arises, which is also uniform for ellipsoids. Since then, the Eshelby tensor has been determined by many authors for inhomogeneities of various shapes, but almost always for the case of uniform eigenstrains. In many application areas in fact, the case of non-uniform eigenstrains is of more physical significance, particularly when the inhomogeneity is non-ellipsoidal. In this article, a method is introduced, which approximates the Eshelby tensor for a variety of shaped inhomogeneities in the case of more complex eigenstrains by employing local polynomial expansions of both the eigenstrain and the resulting Eshelby tensor, in the case of the potential problem in two dimensions.


Introduction
For almost two centuries now, significant interest has been focused on the ability to predict fields arising inside and surrounding inhomogeneities embedded in otherwise uniform host materials. In particular, the isolated inhomogeneity problem where a single inhomogeneity resides inside a host medium that is considered to be of infinite extent, with conditions imposed in the far field, is a classical one. In 1826, Poisson studied the perturbed field due This work is associated with methodology. Therefore, all research data supporting this article are directly available within this publication.
D. Joyce · W. J. Parnell (B) School of Mathematics, University of Manchester, Oxford Road, Manchester M13 9PL, UK e-mail: william.parnell@manchester.ac.uk to an isolated ellipsoid in the context of the Newtonian potential problem [1]. It was shown that the induced electric (or magnetic) field inside the ellipsoid is uniform given a uniform electric polarisation (or magnetisation). In 1873, Maxwell derived explicit expressions for this field [2]. In 1931, Dive [3] and independently in 1932 Nikliborc [4] proved that given an arbitrary imposed uniform electric polarisation in the far field, the electric polarisation induced inside an ellipsoid is also uniform. In the modern language of inhomogeneity problems, this is considered to be proof of the weak Eshelby conjecture for the Newtonian potential problem. This language itself arises following the seminal paper of Eshelby in 1957, who considered the isolated inhomogeneity problem in the case of linear elastostatics and proved that the induced strain in an ellipsoid is uniform given the arbitrary uniform strains in the far field [5]. In 1961, Eshelby conjectured that the ellipsoid is the only shape that has this property [6]. This being true for a uniform arbitrary far-field strain is known as the weak Eshelby conjecture. On the other hand, this being true for a uniform specific far-field strain is known as the strong Eshelby conjecture. Given the similarity of the Newtonian potential and linear elastostatics problems, differing mainly by their tensorial order, they are often discussed simultaneously [7], and hence the use of this term in both contexts. The strong conjecture in the context of the potential problem is true in two dimensions [8,9], but it is not true in dimensions greater than two. This was illustrated in 2008 when Liu gave a non-ellipsoidal counterexample associated with a specific far-field loading [10].
For ellipsoidal inhomogeneity problems, the so-called Eshelby tensor arises naturally, as does its counterpart the Hill tensor. These can be shown to be uniform for homogeneous ellipsoidal inhomogeneities. These tensors arise specifically in the context of micromechanics and effective medium theory, with the aim of predicting the effective properties of inhomogeneous media [7,11]. Speaking in the context of thermal conductivity as we do in this article, components of the Eshelby tensor written with respect to Cartesian coordinates in the isotropic case take the form: where k 0 is the conductivity,D is the ellipsoidal inhomogeneity, andG = 1/(4π k 0 |x −ỹ|) is the associated Green's tensor. Hill's tensor P = k −1 0S . The associated two-dimensional problem is also frequently studied wheñ D is an ellipse, andG = −1/(2π k 0 ) log |x −ỹ|. This tensor arises naturally because it is uniform for ellipsoidal inhomogeneities, and therefore, assuming uniform interior potential gradient e 1 , the associated integral equation gives (as will be shown below) e 1 = e ∞ − (κ − 1)Se 1 , so that e 1 = [I + (κ − 1)S] −1 e ∞ , (1.2) where e ∞ is the uniform far-field potential gradient, and κ = k 1 /k 0 where k 1 is the conductivity of the inhomogeneity.
For non-ellipsoidal inhomogeneities or in studies of, e.g. interacting inhomogeneities, in general, the temperature gradient will not be uniform, and therefore, the Eshelby tensor associated with non-ellipsoidal inhomogeneities is not generally a natural quantity to work with, since it will not arise naturally in the associated integral equation. In spite of this, the non-ellipsoidal problem has been studied frequently in the literature usually in connection with Eshelby's conjecture, to prove that a tensor is not uniform for specific classes of inhomogeneities. Attention has been paid to polygonal and polyhedral inhomogeneities and the associated properties of Eshelby's tensor [12][13][14][15][16][17][18][19]. The case of the so-called supersphere was considered in [20], which built on the work in [21][22][23]. For the case of the Newtonian potential problem and planar elastostatics, some analytical expressions have been derived for twodimensional problems where inhomogeneities are either polygonal or their shape can be described by finite Laurent expansions [24,25]. Further properties of the Eshelby tensor were deduced in [26,27], including the relationship of the averaged Eshelby tensor for non-ellipsoidal inhomogeneities to their ellipsoidal counterparts.
The non-ellipsoidal problem is also considered since the so-called uniform eigenstrains can be considered in inclusion regions within a homogeneous host material. A non-zero eigenstrain field can be used to simulate the presence of an inhomogeneity inside a homogeneous host medium or can additionally be considered to incorporate some other physical effect, boundary presence, interacting inhomogeneities, or mismatch [28]. In order to simulate the presence of an isolated ellipsoidal inhomogeneity with uniform potential gradient imposed in the far field, the required eigenstrain is uniform [28]. Ru [29] developed a method based on conformal mappings to determine the strain field that arises in and around two-dimensional inhomogeneities of arbitrary cross section, due to uniform interior eigenstrain. Importantly, it should be noted that this is not equivalent to solving an inhomogeneity problem with uniform far-field loading because a uniform eigenstrain cannot account for the presence of an inhomogeneity of arbitrary cross section [28].
For the study of interior fields associated with non-ellipsoidal (non-elliptical in two dimensions) inhomogeneities, a more natural quantity to study is what we shall term here the generalised Eshelby tensor taking the form where e * (ỹ) is a component of the so-called eigenstrain. Note that the term eigenstrain is always used regardless of the application area. Rahman [30] derived explicit expressions for the polynomial eigenstrain problem associated with ellipsoidal inclusions in the elastostatics context. For arbitraryD in the context of the thermal problem, if one considers polynomial eigenstrains, then the generalised Eshelby tensors take the form where for generality the polynomial is expanded abouts. An interesting property for ellipsoids is that of polynomial conservation which states that Eshelby's tensor will be a polynomial of order N when the eigenstrain is of order N [28,30]. Regardless of the fact that a generalised Eshelby tensor of the form (1.3) or (1.4) is of more physical significance for non-ellipsoidal problems, such tensors do not appear to have been studied in the literature for general inho-mogeneitiesD. Mura ([28, Sect. 12]) states certain results associated with the polynomial conservation property which arise from results in the ellipsoidal potential theory context due to Ferrers [31] and Dyson [32]. For a more complete treatment of the inhomogeneous ellipsoid problem, see the more recent paper by Rahman [33]. Mura also uses these results in order to approximate the effect of non-uniformity of interior strain fields when two ellipsoidal inhomogeneities interact. However, even for isolated non-ellipsoidal cases, no results appear to have been derived except in the case when N = 0, the classical Eshelby tensor scenario. The case when N ≥ 1 for non-ellipsoidal inhomogeneities is therefore the focus of the present article, and our attention is restricted to the two-dimensional case, i.e. when the generalised Eshelby tensor takes the form withG = −1/(2π k 0 ) log |x −ỹ| and whereD is not restricted to being elliptical. Besides being of interest in their own right, the evaluation of such tensors would assist in developing approximate schemes for interacting inhomogeneity problems and in cases when eigenstrains are non-uniform. They may also be seen as approximations to cases when eigenstrains are non-polynomial, but Taylor series expansions are taken of the eigenstrains, expanded about locations insideD.
In this article, the focus is not in predicting internal fields within inhomogeneities, although this is clearly of interest and will be the focus of future work. Here polynomial approximations shall be determined for the generalised Eshelby tensors in the case of non-ellipsoidal inhomogeneities. In Sect. 2, we introduce the necessary mathematical background discussing both the inhomogeneity and inclusion problems and illustrating where Eshelby's tensor arises in these contexts. Polynomial and more complex eigenstrain fields are then introduced in order to define the concept of a generalised Eshelby tensor. The method of approximating this tensor by via polynomial expansions is then described in Sect. 3. Application of the method to elliptical and non-elliptical domains is then described in Sect. 4. We summarise and discuss future directions in Sect. 5.

Background
The motivation for the consideration of generalised Eshelby tensors is twofold. Firstly in the context of the prediction of induced fields associated with the presence of an isolated inhomogeneity (with different properties to those of the surrounding medium) embedded in an otherwise homogeneous medium, with some imposed condition at infinity. Secondly for the prediction of induced fields due to an isolated inclusion region (with the same properties as those of the homogeneous host medium) but within which a so-called eigenstrain is imposed.
Note the language used here, as introduced by Mura [28] which distinguishes an inclusion (same properties as those of the surrounding medium) from an inhomogeneity (different properties from those of the surrounding medium).

Inhomogeneity problem
Consider an unbounded host medium of infinite extent into which a cylindrical inhomogeneity is embedded, which itself can be considered to be of infinite extent along its axis and place thex 1x2 plane coincident with the cross section of the cylinder. With translational invariance in thex 3 direction then, the problem is purely two-dimensional, as depicted in Fig. 1 where the cross section of the inhomogeneity is defined asD ∈ R 2 . Since it is often useful to consider a specific physical problem, certainly in terms of language and terminology, the potential problem here is described in the context of steady-state thermal conductivity.
The equation governing the steady-state temperature distribution T (x) wherex = (x 1 ,x 2 ) is the Cartesian position vector in the medium described above and depicted in Fig. 1 is where the assumption is that no heat sources are present. The free-space Green's function associated with the host phase satisfies Assuming continuity of temperature and normal flux across ∂D, the resulting temperature distribution may be straightforwardly derived in integral equation form as which holds for allx. Here T ∞ (x) is the solution to the equivalent problem satisfying (2.1) with no inhomogeneity present (or equivalently with C 1 i j = C 0 i j ). Upon taking derivatives of (2.3) with respect tox i and noting the property ∂G/∂ỹ i = −∂G/∂x i it is found that for allx, where the ith component of the temperature gradient has been defined as e i = ∂ T /∂x i . If attention is restricted to anisotropy where the principal axes of the host and inhomogeneity are aligned, all anisotropic potential problem inhomogeneity problems can be reduced to isotropic ones by rescaling the spatial variablex which modifies the inhomogeneity shape. The problem is then solved in the scaled domain and mapped back to physical space [7]. Interest here is for general shaped inhomogeneities, and therefore we consider the isotropic problem only. Therefore take where κ = k 1 /k 0 and we have introduced the associated two-dimensional isotropic Greens functionG satisfying It is important to note that different references place the Dirac delta function on either side of the governing equation for the Green's function. Of course, this merely changes the sign of the Green's function, but it is important in terms of lifting information directly from one reference to another. We follow the convention in [7,28]. Equation (2.5) is an integral equation for the potential gradient e and determining this insideD, which we shall denote as e 1 then tells us the field everywhere. If one assumes that e ∞ is a uniform vector, taking y ∈D and assume that e 1 is also uniform, we arrive at the expression This is simply (1.2), whereS is Eshelby's tensor with components and where deriving (2.7) has exploited the symmetry S i j = S ji . Under these assumptions of uniformity, Eq. (2.7) is then only consistent ifS is a uniform tensor, which, from Eshelby's strong conjecture in two dimensions, we know is only true in the case of elliptical inhomogeneitiesD. In two dimensions then, when the far-field temperature gradient if uniform, then one can use (2.7) in the case of isolated elliptical inhomogeneities to predict the uniform interior fields and Eshelby's tensor takes the form ([7, (4.36)]) where = a 2 /a 1 where a j is the semi-axis of the ellipse directed along the x j axis. However, for non-elliptical inhomogeneities one must return to (2.5) and solve this integral equation. This approach also applies even for elliptical inhomogeneities when the far-field temperature gradient is non-uniform! One approach would be to pose some form for the interior strain in some basis, e.g. a polynomial expansion, e.g. forx 1 ∈D, for some coefficients A j mn , j = 1, 2 to be determined. Inserting these expansions in the integral equation will lead to consistency conditions on the coefficients (a linear system) with coefficients that are integrals of the form Clearly such an expansion (2.10) could not be expected to converge everywhere but can certainly be of use. This method has been utilised to predict the fields interior to interacting ellipsoids in [28, Sect. 23] for example.

Inclusion problem
Consider now a different configuration to the previous section. An unbounded host medium of infinite extent contains within it an inclusion-a cylindrical region which is considered to be of infinite extent along its axis and place thẽ x 1x2 plane coincident with the cross section of the cylinder. The conductivity tensor of the inclusion region is the same as its exterior. So-called eigenstrains (independent ofx 3 ) are imposed in this region which induce a perturbed field. Note that unlike the previous section there is no imposed field at infinity. As with the previous section, with translational invariance in thex 3 direction then, the problem is purely two-dimensional. Following Mura [28] but in the potential context, the equation governing T subject to the eigenstrain field e * is where χ(x) = 1 whenx ∈D and is zero otherwise. Using the Green's function as defined in (2.2) the induced field can be written which holds for allx. As in the previous section, taking derivatives of (2.3) with respect tox i and noting the property ∂G/∂ỹ i = −∂G/∂x i it is found that for allx, (2.14) Following the same arguments as for the inhomogeneity problem if attention is restricted to anisotropy where the principal axes of the host and inhomogeneity are aligned, all anisotropic problems for generalD reduce to the isotropic case due to scaling arguments and therefore set C 0 i j = k 0 δ i j so that (2.14) becomes whereG is defined in (2.6). With eigenstrain components e * j in the form of polynomials, the right-hand side of (2.15) can clearly then be written as linear combinations of terms in the form (2.11).
One can employ eigenstrain to 'represent' the effects of an inhomogeneity [28]. First add some far-field forcing to the problem so that (2.15) becomes Now comparing (2.16) with (2.5), we see that these are equivalent if In the case of an elliptical inhomogeneity with uniform e ∞ i ,S i j is uniform and is given by (2.9), and we have where A = [I + (κ − 1)S] −1 is known as the concentration tensor [7]. For imposed fields that are of polynomial form a polynomial eigenstrain of the same order can be employed to represent an inhomogeneity [28]. For more complex shaped domains, however, a general eigenstrain would be required.

Generalised Eshelby tensors
Let us call (2.11) the generalised Eshelby tensor of order (m, n). Upon defining we havẽ (2.20) Determining explicit results for these tensors for generalD is clearly difficult, if not impossible, unless we restrict attention to specific classes ofD. Instead, then consider an approach that is based on polynomial expansions. The general methodology shall be developed in the next section before it is applied to specific cases. As we shall discuss in the Sect. 3, we will seek polynomial expansions as approximations toJ mn for a givenD. As should be expected, it transpires that the choice of location about which the polynomial is expanded is extremely influential as to the domain of validity of the polynomial expansion. Suppose therefore that the polynomial expansion is centred onx =r. The methodology to follow relies upon orthogonality and thereforeJ mn (x;s) shall first be written as a linear combination ofJ mn (x;s), i.e. we writẽ Further, in an effort to yield a polynomial expansion that is convergent over a larger domain, one can expand about distinct pointsr j , j = 1, 2, . . . , p and then split the domain into non-intersecting support regions sayD j (such thatD =D 1 ∪D 2 ∪ · · · ∪D p ) over which the expansion is locally valid. For example, with reference to Fig. 2 write (2.23) as Fig. 1 Configuration of the problem: a cylindrical inhomogeneity or inclusion regionD is embedded perfectly inside a host medium of infinite extent. Thex 1x2 plane is aligned with the cross section of the cylinder, which is denotedD Fig. 2 The regionD is split into two non-intersecting regionsD 1 andD 2 , boundary of which here is depicted by the dashed line.
A polynomial expansion approximation for the function J mn is sought locally in each region via an expansion aboutx =r j , where j = 1, 2. An expression for the function J mn is then obtained by combining these local expansions via support functions χ(x) where

More complex eigenstrains
Note that in principal more complex eigenstrain fields can be considered in (2.15); particularly in the context of the local expansion approach (2.24). One would Taylor expand the eigenstrain aboutr j , j = 1, 2, . . . , p insideD as above, which leads to generalised Eshelby tensors (polynomial eigenstrains). For example, if we have then inD j , Taylor expand f (ỹ) aboveỹ =r to givẽ so that upon truncating  Fig. 3 Depicting lengthscales of the inhomogeneityD. The local expansion aboutr 1 is scaled on L 1 for convenience. Once this local expansion is determined, in order to contribute to the general expansion scheme, we reintroduce L 1 . Finally, once the global expansion has been determined, this is scaled on some lengthscale, e.g. the radius L of the smallest circle that can be inscribed intoD

Evaluation of the generalised Eshelby tensor for polynomial eigenstrains
Consider the situation where we need to determineJ (x;r j ) as in (2.24), with reference to Fig. 2. As shall be explained below we will first re-scale these local coordinates and then seek polynomial expansions in the scaled coordinates x − r j for this local expansion. In order to use this expression in (2.24) one needs to re-introduce the scaling so that an expansion is determined in each local subregionD j terms ofx −r j . Once a global expansion is determined it is then sensible to non-dimensionalise on some associated lengthscale of the inhomogeneity region L, e.g. smallest circle that can be inscribed in the region whilst still touching the boundary at two points, as depicted in Fig. 3.

Polynomial approximations
For general f (i.e. for a general shaped inhomogeneity D), it is clearly not possible to obtain an analytical form for the expression (3.6); generally, it must be evaluated numerically for any point x ∈ D. Given its importance in a number of inhomogeneity and inclusion problems as indicated above, it would therefore be useful if approximate schemes can be developed in order to generate expressions for the right-hand side. Let us suppose that we may expand J in powers of its argument, e.g. in the form J mn x; r j ≈ J mn x; r j = p+q=P p,q: p+q=0 for non-negative integers p, q and some positive integer P, where the notation J mn has been used to indicate that this is an approximation to J mn . We will see in fact that this approach is useful for many practical cases, especially with an appropriate choice of r j . Note that the upper limit P → ∞ in general. For this local expansion to be useful in (2.24), we need to re-introduce the dimensional variables, so that (3.7) becomes J mn x;r j = where the coefficients a mn j , b mn j are determined explicitly in terms of ρ and D mn pq by employing the plane polar form (3.2) I for x − r j and then expand the powers of the resulting trigonometric functions in their multiple angle forms (this is straightforwardly done algorithmically in a symbolic package such as mathematica) so that for example a mn 0 (ρ) = D mn 00 + and analogously for a mn j , b mn j . Next note that the D mn pq are not independent. We can see this by applying the Laplacian operator to (3.7) and exploiting the form of (3.6) with ∇ 2 G + δ(y − x) = 0 to find Comparing coefficients of x a 1 x b 2 provides equations linking the coefficients D δξ pq , so that for a fixed pair m, n we have for p, q ≥ 0 ( p + 2)( p + 1)D mn ( p+2)q + (q + 2)(q + 1)D mn p(q+2) = δ pm δ qn . (3.12) Next introduce the operators C k ( * ) and S k ( * ) as those that multiply * by cos(kθ) and sin(kθ), respectively and then integrate over θ ∈ [0, 2π). Equate (3.3) and (3.9) and apply the operators C m and S m in turn to each side of the resulting equations to yield (3.13) where δ m0 is the Kronecker delta. Without loss of generality, set ρ = 1, to find that Note again here that by construction, ζ = f j (φ) prescribes the boundary of the fibre with f j (φ) ≥ 1 so that ζ ≥ 1. By employing some exact results for certain integrals, we now show that we are able to write (3.14)-(3.16) in terms of a single integral in φ involving the boundary function f j (φ). To proceed, define (3.19) and this substitution allows us to simplify the expression with various terms reducing to zero by employing double angle formulae for the trigonometric functions. We then obtain  (3.23) where k ≥ 1 in the latter. The ζ integration can then be straightforwardly carried through using integration by parts, to determine that noting that we have introduced the notation F mn 0 for this integral form of a mn 0 so that we can refer to it later on. Proceeding analogously, defining it can be straightforwardly shown that, for k ≥ 1, We now employ the equations relating a m , b m and D pq and equate these to the integral forms just developed above, so that, e.g. · · · = · · · = · · · , (3.31)

32)
· · · = · · · = · · · (3.33) These equations coupled with the conditions (3.12) appropriately truncated gives a closed system from which the coefficients D mn pq are determined, i.e. this can be written where D is a vector of the unknowns D mn pq and F is a vector of the right-hand sides associated with (3.29)-(3.33) and (3.12). It should be noted that the matrix A is identical regardless of the shape of the inhomogeneity and the expansion point r. The latter affect only the right-hand side F. Therefore one can immediately write the explicit form (3.35) and the only computation in order to form the polynomial approximation are the integrals in F. In the next section several examples are considered with various eigenstrains and a number of inhomogeneities of distinct shape, in order to illustrate the efficacy of this scheme.

Method implementation
The aim then is to determineS mn i j (2.20) via polynomial representations of the integral defined in (2.19) for a given domainD, integers m, n and position vectors. We shall always choose min f (θ ) = 1 so that in terms of the original dimensional coordinates L = 1. This means in effect that there is no difference between variables that have a 'tilde' and those that do not above. We shall talk in terms of variables without the tilde for convenience.

Elliptical inclusions: polynomial conservation
First consider the case when s = 0 and the inhomogeneity is elliptical, with where a 1 and a 2 are the semi-axes of the ellipse in the x 1 and x 2 directions, respectively. Note that for convenience, we choose min(a 1 , a 2 ) = 1. As stated above, for elliptical inhomogeneities, polynomial eigenstrains of order m + n will give rise to J mn integrals that are polynomials of order m + n + 2. Our first check is therefore to ensure that the method above reproduces this so-called property of polynomial conservation. Let us compute J 30 which corresponds to a cubic eigenstrain. The method yields J 30 as a fifth-order polynomial and gives rise to the three components of the Eshelby tensor S 30 20 , S 30 11 and S 30 02 . With all coefficients taken to five significant figures, for a 1 = 2, a 2 = 1, these polynomials take the form S 30 20 = −0.098765x 1 + 0.13374x 3 1 + 0.098765x 1 x 2 2 , and surface plots of these components are provided in Fig. 4.

Squares, pentagons and hexagons
Now consider domains of more complex form than elliptical and for now take s = 0. Gielis [35] utilised what is referred to as the superellipse shape function to reproduce many shapes commonly occurring in nature. The function is defined as Here γ can be considered a parameter associated with rotational symmetry, so for a given γ the shape produced will have γ -fold rotational symmetry. Further, a and b are associated with the aspect ratio of the shape, and δ 1 , δ 2 and δ 3 help to further control the form of the shape. For instance, to reproduce an asymmetric shape, δ 1 , δ 2 and δ 3 should be chosen to have distinct values and to produce 'star' shapes, δ 1 should be chosen less than δ 2 and δ 3 . Note that choosing δ 1 = δ 2 = δ 3 = 2 and γ = 4 returns the shape function governing an ellipse, with a = b further reducing the equation to that of a circle or radius a. The function (4.1) is incredibly versatile, and can be used to produce close approximations to many boundary shapes, including polygons, with the essential rule of thumb being to choose a = b = 1, γ to be the degree of polygon sought, δ 2 and δ 3 to be large and equal, and δ 1 to be larger still-except for the square case, where all three δ j parameters are chosen to be equally large. So far no "re-expansion" has been required in the sense that the natural "origin" of the domains considered was the origin itself. Let us now consider the case of the so-called limacon, whose boundary is governed by the polar function  should remain along the x 1 axis. In fact, the most obvious way of maximising the radius of convergence would be to move the basis to r = (1, 0), as a circle of radius 2 based at this point will not only fit inside the limacon, but also cover the vast majority of the area of this domain. Figure 8b illustrates how the error improves when this re-expansion point is chosen, with the regions previously featuring the most extreme errors now featuring errors in the region of 10 −7 .

The 'hourglass' and hypogenis shapes
Let us now consider some more complex shapes of non-polygonal form. By choosing γ = 4, a = 1, b = 10, δ 1 = 4.2, δ 2 = 17 and δ 3 = 1.5 in the superellipse expression (4.1), we obtained what we shall term the 'hourglass' due to its resemblance to the traditional time-keeping device, as can be seen in Fig. 9a and take s = 0.
As one would expect, basing an expansion of the polynomial approximation at the origin yields large errors in the extreme north and south regions of the domain, as illustrated in Fig. 10a, focusing on the x 2 ≥ 0 region and illustrated in the case of J 12 . The symmetry of the domain is such that it makes sense to split it into two domains: x 2 ≥ 0 and x 2 < 0, and use expansions about two points on the x 2 axis equidistant from the origin in each. Once again, focusing on x 2 ≥ 0, Fig. 10b shows how the error improves upon centering the expansion about r = (0, 1.13).
An alternative to consider is the hypogenis, as depicted in Fig. 9b, boundary of which is described by with the factor of 3/2 giving rise to f (θ ) ≥ 1. In a similar vein to the 'hourglass', here it is advantageous to split the domain into the sub-domains x 1 ≥ 0 and x 1 < 0 and use expansions about two points on the x 1 axis equidistant from the origin for each. Figure 11 focuses on the x 1 ≥ 0 region and shows that an expansion about r = (0.82, 0) improves upon an expansion about the origin.  = (1, 0), but the inclusion shape remains centred at the origin. Here the domain to be examined is once again a square. In the case of, for instance, J 20 (x; (1, 0)), the expansion formula (2.23) will need to be utilised for r = 0 in order to express J 20 (x; (1, 0)) in terms of J mn (x; 0) integrals-terms which themselves shall need to be approximated by polynomial expansions about the origin. Figure 12 illustrates how increasing the order of these expansions improves the overall accuracy of the approximation of J 20 (x; (1, 0)).

Concluding remarks
A new scheme to determine polynomial approximations to integrals over non-elliptical two-dimensional domains has been developed. We call these integrals generalised Eshelby tensors for the Newtonian potential problem since they correspond to taking non-uniform eigenstrains over the domain, and when the eigenstrain is uniform, this corresponds to the classical Eshelby tensor. Excellent performance of the method is illustrated by comparison with the numerical evaluation of the integral, noting that this has to be evaluated for every point in the domain of interest, whereas once determined, the polynomial approximation can be rapidly evaluated at any point. Extensions to the two-dimensional in-plane elasticity case are underway and also to the full three-dimensional scalar and tensor scenarios. Further investigation is also underway that will explain theoretically the domains of convergence of the polynomial expansions.