Elastic strain field due to an inclusion of a polyhedral shape with a non-uniform lattice misfit

TL;DR

This paper derives an analytical solution for the three-dimensional elastic strain field caused by a polyhedral inclusion with a spatially varying lattice misfit, applicable to quantum dots and similar structures.

Contribution

It introduces a closed-form analytical method to compute elastic strain for polyhedral inclusions with polynomial lattice mismatch, extending previous models to non-uniform misfits.

Findings

Strain tensor expressed as sum over faces, edges, vertices

Contributions involve polynomials and simple geometric functions

Method applicable to electrostatic and gravitational potentials

Abstract

An analytical solution in a closed form is obtained for the three-dimensional elastic strain distribution in an unlimited medium containing an inclusion with a coordinate-dependent lattice mismatch (an eigenstrain). Quantum dots consisting of a solid solution with a spatially varying composition are examples of such inclusions. It is assumed that both the inclusion and the surrounding medium (the matrix) are elastically isotropic and have the same Young modulus and Poisson ratio. The inclusion shape is supposed to be an arbitrary polyhedron, and the coordinate dependence of the lattice misfit, with respect to the matrix, is assumed to be a polynomial of any degree. It is shown that, both inside and outside the inclusion, the strain tensor is expressed as a sum of contributions of all faces, edges and vertices of the inclusion. Each of these contributions, as a function of the observation point's coordinates, is a product of some polynomial and a simple analytical function, which is the solid angle subtended by the face from the observation point (for a contribution of a face), or the potential of the uniformly charged edge (for a contribution of an edge), or the distance from the vertex to the observation point (for a contribution of a vertex). The method of constructing the relevant polynomial functions is suggested. We also found out that similar expressions describe an electrostatic or gravitational potential, as well as its first and second derivatives, of a polyhedral body with a charge/mass density that depends on coordinates polynomially.

Figures

b of Ref. Glas, [2001](#bib.bib12). That is, our method provides the same strain distribution as the closed-form analytical solution found by GlasGlas ([2001](#bib.bib12)) for the particular case of a truncated pyramid with a constant misfit strain.

Equations

$${\epsilon }_0(x,y,z)=\frac{a(x,y,z)-a_0}{a_0}.$$

$${\epsilon }_0(\mathbf{r})=\{\begin{array}{ll}f(\mathbf{r})& \text{if }\mathbf{r}\text{ is inside the inclusion},\\ 0& \text{if }\mathbf{r}\text{ is outside the inclusion},\end{array}$$

$$f(x,y,z)=\sum _{\begin{array}{c}m,n,p\\ m+n+p\le N\end{array}}{C}_{mnp}{x}^m{y}^n{z}^p.$$

$${\epsilon }_0(\mathbf{r})=f(\mathbf{r})\chi (\mathbf{r}),$$

$$\chi (\mathbf{r})=\{\begin{array}{ll}1& \text{if }\mathbf{r}\text{ is inside the inclusion},\\ 0& \text{otherwise}.\end{array}$$

$${\epsilon }_{\alpha \beta }(\mathbf{R})=?$$

$$\Delta \phi (\mathbf{r})=-4\pi {\epsilon }_0(\mathbf{r}),$$

$$\mathbf{u}(\mathbf{r})\propto \frac{\mathbf{r}}{r^3}.$$

$$\mathbf{u}(\mathbf{r})=-\Lambda \frac{\partial \phi (\mathbf{r})}{\partial \mathbf{r}},$$

$$\Lambda =\frac{1}{4\pi }\frac{1+\nu }{1-\nu }.$$

$${\epsilon }_{\alpha \beta }(\mathbf{r})=\frac{1}{2}\left(\frac{\partial u_{\alpha }}{\partial r_{\beta }}+\frac{\partial u_{\beta }}{\partial r_{\alpha }}\right)-{\delta }_{\alpha \beta }{\epsilon }_0(\mathbf{r}),$$

$${\epsilon }_{\alpha \beta }(\mathbf{r})=-\Lambda \frac{{\partial }^2\phi (\mathbf{r})}{\partial r_{\alpha }\partial r_{\beta }}-{\delta }_{\alpha \beta }{\epsilon }_0(\mathbf{r}).$$

$$\phi (\mathbf{R})=\iiint {\epsilon }_0(\mathbf{r})\frac{d^{3}r}{|\mathbf{r}-\mathbf{R}|}\equiv \underset{\text{inclusion}}{\iiint }f(\mathbf{r})\frac{d^{3}r}{|\mathbf{r}-\mathbf{R}|},$$

$$\frac{\partial \phi (\mathbf{R})}{\partial \mathbf{R}}=\iiint \frac{\partial {\epsilon }_0(\mathbf{r})}{\partial \mathbf{r}}\frac{d^{3}r}{|\mathbf{r}-\mathbf{R}|},$$

$$\frac{{\partial }^2\phi (\mathbf{R})}{\partial R_{\alpha }\partial R_{\beta }}=\iiint \frac{{\partial }^2{\epsilon }_0(\mathbf{r})}{\partial r_{\alpha }\partial r_{\beta }}\frac{d^{3}r}{|\mathbf{r}-\mathbf{R}|}.$$

$${\rho }_{mnp}(\mathbf{r};\mathbf{R})=(x-X{)}^m(y-Y{)}^n(z-Z{)}^p\chi (\mathbf{r}),$$

$${\epsilon }_0(\mathbf{r})=\sum _{\begin{array}{c}m,n,p\\ m+n+p\le N\end{array}}{\tilde{C}}_{mnp}(\mathbf{R}){\rho }_{mnp}(\mathbf{r};\mathbf{R}).$$

$$\begin{gathered} \sum _{m,n,p}{C}_{mnp}{x}^m{y}^n{z}^p \\ =\sum _{m,n,p}{\tilde{C}}_{mnp}(\mathbf{R})(x-X{)}^m(y-Y{)}^n(z-Z{)}^p, \end{gathered}$$

$$\begin{gathered} {\tilde{C}}_{mnp}(\mathbf{R})= \\ \sum _{\begin{array}{c}{m}^{\prime }⩾m\\ n^{\prime }⩾n\\ p^{\prime }⩾p\end{array}}{C}_{m^{\prime }{n}^{\prime }{p}^{\prime }}\left(\genfrac{}{}{0ex}{}{m^{\prime }}{m}\right)\left(\genfrac{}{}{0ex}{}{n^{\prime }}{n}\right)\left(\genfrac{}{}{0ex}{}{p^{\prime }}{p}\right)X^{m^{\prime }-m}{Y}^{n^{\prime }-n}{Z}^{p^{\prime }-p}, \end{gathered}$$

$$\phi (\mathbf{R})=\sum _{m,n,p}{\tilde{C}}_{mnp}(\mathbf{R}){\phi }_{mnp}(\mathbf{R}),$$

$$\frac{\partial \phi (\mathbf{R})}{\partial R_{\alpha }}=\sum _{m,n,p}{\tilde{C}}_{mnp}(\mathbf{R}){\phi }_{,\alpha }^{(mnp)}(\mathbf{R}),$$

$$\frac{{\partial }^2\phi (\mathbf{R})}{\partial R_{\alpha }\partial R_{\beta }}=\sum _{m,n,p}{\tilde{C}}_{mnp}(\mathbf{R}){\phi }_{,\alpha \beta }^{(mnp)}(\mathbf{R}),$$

$${\phi }_{mnp}(\mathbf{R})=\iiint (x-X{)}^m(y-Y{)}^n(z-Z{)}^p\chi (\mathbf{r})\frac{d^{3}r}{|\mathbf{r}-\mathbf{R}|},$$

$$\varphi^{(mnp)}_{,\alpha}(\mathbf{R})\\ =\iiint\frac{\partial\big{[}(x-X)^{m}(y-Y)^{n}(z-Z)^{p}\,\chi(\mathbf{r})\big{]}}{\partial r_{\alpha}}\;\frac{d^{3}r}{|\mathbf{r}-\mathbf{R}|}\,,$$

$$\varphi^{(mnp)}_{,\alpha\beta}(\mathbf{R})\\ =\iiint\frac{\partial^{2}\big{[}(x-X)^{m}(y-Y)^{n}(z-Z)^{p}\,\chi(\mathbf{r})\big{]}}{\partial r_{\alpha}\partial r_{\beta}}\;\frac{d^{3}r}{|\mathbf{r}-\mathbf{R}|}\,.$$

$${\phi }_{mnp}(\mathbf{R})=\underset{\text{inclusion}}{\iiint }\frac{(x-X{)}^m(y-Y{)}^n(z-Z{)}^p}{|\mathbf{r}-\mathbf{R}|}dV.$$

$${\Phi }_{mnp}^{(i)}(\mathbf{R})=\underset{\text{face}i}{\iint }\frac{(x-X{)}^m(y-Y{)}^n(z-Z{)}^p}{|\mathbf{r}-\mathbf{R}|}dS.$$

$${\Omega }_{mnp}^{(i)}(\mathbf{R})=\underset{\text{face}i}{\iint }(x-X{)}^m(y-Y{)}^n(z-Z{)}^p\frac{(\mathbf{R}-\mathbf{r})\cdot {\mathbf{n}}_i}{|\mathbf{r}-\mathbf{R}{|}^3}dS,$$

$$L_{mnp}^{(k)}(\mathbf{R})=\underset{{\text{A}}_k}{\overset{{\text{B}}_k}{\int }}\frac{(x-X{)}^m(y-Y{)}^n(z-Z{)}^p}{|\mathbf{r}-\mathbf{R}|}dl,$$

$${\phi }_{,\alpha }^{(mnp)}={\phi }_{(mnp),\alpha }-\sum _{\begin{array}{c}i\\ \text{(faces)}\end{array}}{n}_{i\alpha }{\Phi }_{mnp}^{(i)},$$

Full text

quantum dot, eigenstrain, inclusion, analytical solution, potential

I Introduction

Self-assembled quantum dots are inclusions of one semiconducting material within another material (a matrix). Due to different lattice constants of the two semiconductors, such structures possess built-in elastic strain, when grown coherently. The strain plays a key role both in the island self-assembly during heteroepitaxy,Stangl, Holý, and Bauer (2004) and in electronic properties of quantum dots due to the strain effect on the carrier dispersion law.Bir and Pikus (1974); Van de Walle (1989) The strain is especially important in type-II heterostructures with quantum dots. For example, the electron localization in Ge/Si quantum dots occurs in potential wells formed in Si matrix near the Ge inclusion due to a spatially inhomogeneous strain.Dvurechenskii, Nenashev, and Yakimov (2002); Zinovieva et al. (2013) The knowledge of the spatial distribution of the strain induced by quantum dots is, therefore, important for the analysis of their electronic structure. It is natural in this context that many theoretical works on electronic properties of epitaxial quantum dotsGrundmann, Stier, and Bimberg (1995); Pryor (1998); Stier, Grundmann, and Bimberg (1999); Dvurechenskii, Nenashev, and Yakimov (2002); Stoleru, Pal, and Towe (2002) begin with considerations of the strain distribution.

There exist a few inclusion shapes, for which the problem of finding the strain distribution has an analytical solution in a closed form. Examples are two-dimensional problems of inclusions having polygonal shapesRodin (1996); Faux, Downes, and O’Reilly (1997); Glas (2001); Jiang and Pan (2004); Zou, He, and Zheng (2011); Sun, Xu, and Pan (2012); Nenashev, Koshkarev, and Dvurechenskii (2013); Chen, Xu, and Pan (2014); Yue et al. (2015); Lee, Zou, and E. (2015) and some flat shapes with a curvilinear boundary;Faux, Downes, and O’Reilly (1997); Zou, He, and Zheng (2011) three-dimensional problems of ellipsoidalEshelby (1957, 1959, 1961); Kinoshita and Mura (1971); Asaro and Barnett (1975); Mura and Kinoshita (1978); Kinoshita and Mura (1986); Yu, Sanday, and Chang (1994); Rahman (2002) and polyhedralChiu (1977, 1978); Rodin (1996); Downes, Faux, and O’Reilly (1997); Pearson and Faux (2000); Nozaki and Taya (2000); Glas (2001); Li and Anderson (2001); Ovid’ko and Sheinerman (2005); Kuvshinov (2008); Nenashev and Dvurechenskii (2010) inclusions. An extensive list of bibliography on this issue can be found, for example, in Refs. Ovid’ko and Sheinerman, 2005; Kuvshinov, 2008; Gao and Liu, 2012; Wang, Yu, and Wang, 2016. The present work is devoted to a three-dimensional polygonal inclusion. This choice is motivated by the facts that quantum dots often have a faceted surface, and that a polyhedron is a convenient approximation to a body of an arbitrary shape.

The aim of the present study consists in finding an analytical expression (via elementary functions) for the three-dimensional distribution of the strain tensor $\varepsilon_{\alpha\beta}(x,y,z)$ within an inclusion (a quantum dot) and in the surrounding matrix, taking into account spatial inhomogeneity of the misfit strain in the inclusion. We will obtain, in this paper, an analytical solution for an arbitrary polyhedral-shaped inclusion with the lattice misfit $\varepsilon_{0}(x,y,z)$ described by any polynomial function of coordinates. The lattice misfit is defined as a relative deviation of the unstrained lattice constant $a$ (which is determined by the composition of the material in a given point inside the inclusion) from the lattice constant of the environment $a_{0}$:

[TABLE]

The quantity $\varepsilon_{0}$ is often referred to as eigenstrain. The content of epitaxial quantum dots generally represents a solid solution with a composition smoothly varying in space.Offermans et al. (2005); Stoffel et al. (2009); Picco et al. (2012) This provides the motivation for considering a coordinate-dependent lattice misfit $\varepsilon_{0}(x,y,z)$.

Our analysis is based on the following simplifying assumptions. The elastic strain is considered in the continuous-medium, isotropic approximation. The lattice misfit is also isotropic. The strain is small. The inclusion and the surrounding matrix possess the same elastic properties, i. e., have the same values of Young modulus and Poisson ratio. The matrix is spread to infinity in all directions.

Let us briefly survey what is currently known about analytical expressions for the strain induced by inclusions with a spatially varied lattice misfit. For ellipsoidal inclusions, such expressions have been known for several decades,Eshelby (1961); Asaro and Barnett (1975); Rahman (2002) and generalized for the case of anisotropic medium.Mura and Kinoshita (1978) For a two-dimensional problem, where the inclusion shape and composition depend on two coordinates only, analytical answers were also obtained in the cases of linear,Sun, Xu, and Pan (2012); Chen, Xu, and Pan (2014) quadraticYue et al. (2015) and, finally, arbitrary polynomial coordinate dependence of the lattice misfit,Lee, Zou, and E. (2015) taking into account also elastic anisotropy and piezoeffect.

The problem of the polyhedral shaped inclusion with a constant lattice misfit was solved analytically in particular cases (a cuboid, a pyramid),Chiu (1977, 1978); Downes, Faux, and O’Reilly (1997); Pearson and Faux (2000); Glas (2001); Li and Anderson (2001); Ovid’ko and Sheinerman (2005) as well as for the general polyhedron.Rodin (1996); Nozaki and Taya (2000); Kuvshinov (2008); Nenashev and Dvurechenskii (2010) KuvshinovKuvshinov (2008) showed that an analytical expression for the strain distribution should exist in the case of polynomial coordinate dependence of the lattice misfit, but he did not provide an explicit analytical formula. One can find, however, in the literature, the explicit solutions for a related problem—calculation of the Newtonian potential and its derivatives for a massive body of a polyhedral shape with an inhomogeneously distributed mass. This problem, being mathematically equivalent to the elasticity problem (see Section II), was solved however only in the simplest cases of linearPohánka (1998); Holstein (2003); Hamayun, Prutkin, and Tenzer (2009); D’Urso (2014); Conway (2015) and quadraticBhaskara Rao, Prakash, and Ramesh Babu (1990); Gallardo-Delgado, Pérez-Flores, and Gómez-Treviño (2003); Gokula and Sastry (2015) dependence of the mass density on coordinates. Thus, the task of the present paper has not been solved yet.

The present study is based on the approach developed in our previous work, Ref. Nenashev and Dvurechenskii, 2010. In this work, the strain due to a polyhedral inclusion with a constant lattice misfit is shown to be a sum of contributions of faces and edges of the polyhedron. A contribution of a face is proportional to the solid angle subtended by this face. (The role of solid angles in calculating the potential and the strain is mentioned also in Refs. Werner and Scheeres, 1996; Downes, Faux, and O’Reilly, 1997; Kuvshinov, 2008.) A contribution of an edge is proportional to the potential induced of this edge, as if it were uniformly charged (see also Ref. Werner and Scheeres, 1996). In the present paper, this result is generalized to the case of a non-uniform lattice misfit in the following way. At first, the coefficients of proportionality at face and edge contributions become polynomial functions of coordinates of the observation point. At second, the contributions of vertices of the polyhedron to the strain tensor are also introduced; for each vertex, its contribution is equal to the product of some polynomial function of coordinates and the distance between this vertex and the observation point. We prove that such modification of the results of Ref. Nenashev and Dvurechenskii, 2010 indeed provides a solution for the problem of the present work, and suggest an algorithm for finding the coefficients of all needed polynomials.

Before proceeding to the development of the solution, let us introduce some notations. The coordinate-dependent lattice misfit $\varepsilon_{0}(\mathbf{r})$ defined by Eq. (1) can be represented as

[TABLE]

with $f(\mathbf{r})$ being a polynomial function of coordinates $x,y,z$ of some degree $N$:

[TABLE]

Sometimes it is convenient to represent the function $\varepsilon_{0}(\mathbf{r})$ as

[TABLE]

where $\chi(\mathbf{r})$ is the characteristic function that determines whether the point $\mathbf{r}$ belongs to the inclusion:

[TABLE]

Our goal is to calculate the strain tensor $\varepsilon_{\alpha\beta}$ at an arbitrary point $\mathbf{R}=(X,Y,Z)$:

[TABLE]

Also, as we will see in Section II, an important function for our consideration is the electrostatic potential $\varphi(\mathbf{r})$ that is produced by the charge distribution defined by the function $\varepsilon_{0}(\mathbf{r})$. More formally, $\varphi(\mathbf{r})$ is the solution of Poisson equation

[TABLE]

which vanishes at infinity. We will also represent the analytical expressions for the potential $\varphi(\mathbf{r})$ and its first and second derivatives.

This paper is organized as follows. The connection between the strain $\varepsilon_{\alpha\beta}$ and the potential $\varphi$ is explained in Section II. In Sections III and IV we will express the potential $\varphi$ in terms of more simple quantities—four kinds of “primitives” introduced in Section IV. They can be considered as potentials induced by a charge inhomogeneously distributed over the inclusion, or over some its face, or over some its edge, or by a double layer inhomogeneously covering one of the inclusion faces. Then, in Sections V and VI we will show how to find analytical expressions for the “primitives”. All these results are summed up in Section VII, yielding the final analytical formula for the strain distribution $\varepsilon_{\alpha\beta}(\mathbf{R})$ (Subsection VII.1) and the algorithm for finding the polynomial coefficients contributing to this formula (Subsection VII.2). Analytical results for the potential $\varphi$ and its derivatives are also supplied. Section VIII provides some examples of applying the general answer to particular cases. Concluding remarks are presented in Section IX.

II Elastic-electrostatic analogy

The displacement field $\mathbf{u}(\mathbf{r})$ in an isotropic elastic medium caused by a point-like isotropic inclusion (put at $\mathbf{r}=0$) looks just like the electric field of a point charge:

[TABLE]

This simple fact provides a basis for the analogy between the problems of electrostatics and of elasticity theory. More precisely, if the elastic medium is infinite and isotropic, and its Young module and Poisson ratio are constant throughout the whole space, then the displacement induced by the eigenstrain distribution $\varepsilon_{0}(\mathbf{r})$ can be expressed as followsLandau and Lifshitz (1986); Davies (1998):

[TABLE]

where $\varphi(\mathbf{r})$ is the potential defined by Eq. (7), and the coefficient $\Lambda$ is related to the Poisson ratio $\nu$:

[TABLE]

Since the strain tensor $\varepsilon_{\alpha\beta}$ is defined via derivatives of the displacement,

[TABLE]

then it is a combination of second derivatives of $\varphi$:Davies (1998); Nenashev and Dvurechenskii (2010)

[TABLE]

Here $\delta_{\alpha\beta}$ is the Kronecker delta.

Obviously, the potential $\varphi$ has an integral representation:

[TABLE]

as well as its first and second derivatives:

[TABLE]

[TABLE]

(Triple integrals without specification of the domain are assumed to be over the whole space.) The latter equation provides also an integral representation for the strain tensor, by virtue of Eq. (12). It is important to note that integrals (13)–(15) were extensively studied in the context of geophysicalNagy (1966); Paul (1974); Okabe (1979); Pohánka (1988); Holstein and Ketteridge (1996); Pohánka (1998); Bhaskara Rao, Prakash, and Ramesh Babu (1990); Holstein (2002, 2003); Gallardo-Delgado, Pérez-Flores, and Gómez-Treviño (2003); García-Abdeslem (2005); Hamayun, Prutkin, and Tenzer (2009); D’Urso (2014); Gokula and Sastry (2015); D’Urso (2015) and astronomy/spacecraftWerner (1994); Werner and Scheeres (1996); Broucke (1999); Conway (2015) problems. But, as mentioned in Introduction, fully analytical expressions for polyhedral bodies were obtained only when the density $f(\mathbf{r})$ is a constant ($N$=0), a linear ($N$=1) or a quadratic ($N$=2) function of coordinates.

In the next Section, we transform these integrals into a form more convenient for further analysis.

III From polynomials to monomials

It is easier to deal with potentials induced by charge density distributions of the following form:

[TABLE]

where $m,n,p$ are non-negative integer numbers, and the characteristic function $\chi(\mathbf{r})$ is defined in Eq. (5). One can expand the “actual” density $\varepsilon_{0}(\mathbf{r})$ into a series of terms $\rho_{mnp}$:

[TABLE]

Each of the coefficients $\tilde{C}_{mnp}$ is a polynomial of $X,Y,Z$ of degree $N-m-n-p$. The explicit expression for polynomials $\tilde{C}_{mnp}(\mathbf{R})$ can be obtained by substituting Eqs. (2), (3) and (16) into Eq. (17), that leads to the identity

[TABLE]

whence

[TABLE]

where symbols $\binom{m^{\prime}}{m}$ are binomial coefficients.

Substituting the expansion (17) into Eqs. (13)–(15), one can get similar expansions for the potential $\varphi(\mathbf{R})$ and its derivatives:

[TABLE]

[TABLE]

[TABLE]

where symbols $\varphi_{mnp}$, $\varphi^{(mnp)}_{,\alpha}$, $\varphi^{(mnp)}_{,\alpha\beta}$ denote potentials and their derivatives induced by “monomial” charge distributions $\rho_{mnp}$ defined in Eq. (16):

[TABLE]

[TABLE]

[TABLE]

IV Four primitives $\varphi_{mnp}$, $\Phi^{(i)}_{mnp}$, $\Omega^{(i)}_{mnp}$, $L^{(k)}_{mnp}$

We will consider four types of “primitive constituents” of the potential and its derivatives. The first one, $\varphi_{mnp}$, is the potential provided at point $\mathbf{R}=(X,Y,Z)$ by the charged inclusion (Fig. 1a) with the charge density depending on coordinates $x,y,z$ as $(x-X)^{m}(y-Y)^{n}(z-Z)^{p}$:

[TABLE]

This expression is just a different form of Eq. (23).

The second primitive, $\Phi^{(i)}_{mnp}$, is the potential of one of the inclusion faces (number $i$) charged with surface charge density $(x-X)^{m}(y-Y)^{n}(z-Z)^{p}$ (Fig. 1b):

[TABLE]

The third primitive, $\Omega^{(i)}_{mnp}$, is the potential of one of $i$th face covered with a thin dipole layer with dipole moment density $(x-X)^{m}(y-Y)^{n}(z-Z)^{p}$:

[TABLE]

where $\mathbf{n}_{i}$ is the unit normal vector to $i$th face, directed outside the inclusion (see Fig. 1c).

And the fourth one, $L^{(k)}_{mnp}$, is the potential due to $k$th edge of the inclusion, being charged with linear charge density $(x-X)^{m}(y-Y)^{n}(z-Z)^{p}$:

[TABLE]

where $\text{A}_{k}$ and $\text{B}_{k}$ are two ends of the $k$th edge (Fig. 1d).

Significance of these primitives for our study consists in the fact that the quantities $\varphi_{mnp}$, $\varphi^{(mnp)}_{,\alpha}$, $\varphi^{(mnp)}_{,\alpha\beta}$ that appear in Eqs. (20)–(22) can be expressed through them. As shown in Appendix A, $\varphi^{(mnp)}_{,\alpha}$ has the representation

[TABLE]

where $\varphi_{(mnp),\alpha}$ is the following volume integral:

[TABLE]

One can conclude from comparison between Eqs. (23) and (31) that

[TABLE]

In order to write down a representation of $\varphi^{(mnp)}_{,\alpha\beta}$ in terms of the primitives, let us introduce, for each face $i$ and each edge $k$ adjacent to this face, a unit vector $\mathbf{b}_{ik}$ directed perpendicularly to $k$th edge, along $i$th face and out of this face (see Fig. 2). Also, for each edge $k$, we define a tensor $\lambda^{(k)}_{\alpha\beta}$ as follows:Nenashev and Dvurechenskii (2010)

[TABLE]

where $i_{1}$ and $i_{2}$ are numbers of the two faces, which intersect at $k$th edge, as illustrated in Fig. 2. Then, the quantity $\varphi^{(mnp)}_{,\alpha\beta}$ takes the following representation (for derivation, see Appendix A):

[TABLE]

where $\varphi_{(mnp),\alpha\beta}$ and $\Phi^{(i)}_{(mnp),\alpha}$ are defined as follows:

[TABLE]

[TABLE]

From the comparison between Eqs. (23) and (35) it is easy to express the quantities $\varphi_{(mnp),\alpha\beta}$ through primitives $\varphi_{m^{\prime}n^{\prime}p^{\prime}}$:

[TABLE]

and so on. Similarly, the comparison between Eqs. (27) and (36) provides the following relations:

[TABLE]

In the simplest case of $m=n=p=0$, Eq. (34) is simplified toWerner and Scheeres (1996); Nenashev and Dvurechenskii (2010)

[TABLE]

The primitives $\Omega^{(i)}_{000}$ and $L^{(k)}_{000}$ have simple physical and geometrical meanings and well-known analytical representations. Below we will drop the index “000” at them:

[TABLE]

The quantity $\Omega^{(i)}(\mathbf{R})$ is the potential at point $\mathbf{R}$ due to a flat polygon ($i$th face of the inclusion) covered by a uniform dipole layer of unit density, see Fig. 3. The absolute value of $\Omega^{(i)}(\mathbf{R})$ is equal to the solid angle, at which the polygon is seen from the point $\mathbf{R}$.Tamm (1979); Stratton (1941) And the sign of $\Omega^{(i)}(\mathbf{R})$ is positive if the normal vector $\mathbf{n}_{i}$ is directed towards the point $\mathbf{R}$, and negative otherwise. There are several recipes in the literature, how to calculate such solid angles analytically.Oosterom and Strackee (1983); Werner and Scheeres (1996); Nenashev and Dvurechenskii (2007)

It is important to note that the characteristic function of the inclusion $\chi(\mathbf{r})$ can be expressed in terms of solid angles:

[TABLE]

Indeed, if the point $\mathbf{r}$ is inside the inclusion, then each of the quantities $\Omega^{(i)}(\mathbf{r})$ is negative, so their sum is equal to $-4\pi$ (the total solid angle). When $\mathbf{r}$ is outside the inclusion, then some of quantities $\Omega^{(i)}(\mathbf{r})$ are positive, and some are negative; moreover, the sum of negative $\Omega$’s exactly compensates the sum of positive ones. From Eqs. (4) and (41) one can get a representation of the lattice misfit $\varepsilon_{0}$:

[TABLE]

The quantity $L^{(k)}(\mathbf{R})$ is the potential at point $\mathbf{R}$ of a thin, uniformly charged rod ($k$th edge) with unit linear charge density, see Fig. 3. Elementary integration provides the following expression for $L^{(k)}$:Werner and Scheeres (1996); Nenashev and Dvurechenskii (2010)

[TABLE]

where $\mathrm{RA}_{k}$, $\mathrm{RB}_{k}$ and $\mathrm{A}_{k}\mathrm{B}_{k}$ are distances between the point $\mathbf{R}$ and two end points $\mathrm{A}_{k}$ and $\mathrm{B}_{k}$ of $k$th edge.

The main results of this Section are equations (30), (32), (34), (37), (38), which (together with the results of the previous sections) allow one to express the strain tensor, the potential and derivatives of the potential via four kinds of primitives.

V Reducing primitives

$\varphi_{mnp}$, $\Phi^{(i)}_{mnp}$, $\Omega^{(i)}_{mnp}$ to solid angles and line integrals

In this Section, we will demonstrate how to simplify the functions $\varphi_{mnp}(\mathbf{R})$, $\Phi^{(i)}_{mnp}(\mathbf{R})$ and $\Omega^{(i)}_{mnp}(\mathbf{R})$ by reducing them to the solid angles $\Omega^{(i)}(\mathbf{R})$ and to some line integrals along the edges of the polyhedron. For this purpose, we will use two linear functions of $\mathbf{R}$: $h_{i}(\mathbf{R})$ and $B_{ik}(\mathbf{R})$, whose meaning is illustrated in Fig. 4. The function $h_{i}(\mathbf{R})$ is defined as follows:

[TABLE]

where $\mathbf{r}_{i}$ is a radius vector of some point on $i$th face. (The value of $h_{i}$ does not depend on the choice of the point $\mathbf{r}_{i}$, because the vector $\mathbf{n}_{i}$ is orthogonal to the face.) The absolute value of $h_{i}$ has the meaning of the distance from $\mathbf{R}$ to the plane of $i$th face.

The function $B_{ik}(\mathbf{R})$ has the following definition:

[TABLE]

where $\mathbf{r}_{k}$ is a radius vector of some point on $k$th edge. Since the vector $\mathbf{b}_{ik}$ is orthogonal to the edge, the choice of the point $\mathbf{r}_{k}$ does not matter. The absolute value of $B_{ik}$ is the distance between the projection of the point $\mathbf{R}$ to $i$th face and the line of $k$th edge, as illustrated in Fig. 4.

We start from the usual recipe of reducing volume integrals to surface integrals. Applying Gauss’s theorem to the vector

[TABLE]

one can reduce the volume integral $\varphi_{mnp}$ to integrals over inclusion faces. The details of the derivation can be found in Appendix B. The answer reads:

[TABLE]

Now let us consider the primitives $\Phi^{(i)}_{mnp}$ and $\Omega^{(i)}_{mnp}$ for one of the faces (say, $i$th face). It is convenient to choose temporarily a system of coordinates $\hat{x},\hat{y},\hat{z}$ related to $i$th face, so that the axis $\hat{z}$ is perpendicular to this face. For that, we choose three mutually orthogonal normal vectors $\hat{\mathbf{e}}_{x},\hat{\mathbf{e}}_{y},\hat{\mathbf{e}}_{z}$ such that $\hat{\mathbf{e}}_{z}=\mathbf{n}_{i}$, and define the coordinates $\hat{x},\hat{y},\hat{z}$ as follows:

[TABLE]

(see Fig. 4). Then, let us define quantities $\hat{\Phi}^{(i)}_{mnp}$, $\hat{\Omega}^{(i)}_{mnp}$ and $\hat{L}^{(k)}_{mnp}$ similar to the primitives $\Phi^{(i)}_{mnp}$, $\Omega^{(i)}_{mnp}$, $L^{(k)}_{mnp}$ but related to the “tilted” axes $\hat{x},\hat{y},\hat{z}$:

[TABLE]

where $\hat{r}=\left(\hat{x}^{2}+\hat{y}^{2}+\hat{z}^{2}\right)^{1/2}=|\mathbf{r}-\mathbf{R}|$. (The considered edges are those that surround $i$th face.) Obviously, the “untilted” quantities are linear combinations of “tilted” ones:

[TABLE]

with coefficients $T^{m^{\prime}n^{\prime}p^{\prime}}_{mnp}$ defined by the identity

[TABLE]

Note that, at $m=n=p=0$, the quantities $\Phi^{(i)}_{mnp}$, $\Omega^{(i)}_{mnp}$ and $L^{(k)}_{mnp}$ are invariant with respect to the choice of coordinates. Consequently,

[TABLE]

Then, it is easy to get rid of the index $p$ in $\hat{\Phi}^{(i)}_{mnp}$ and $\hat{\Omega}^{(i)}_{mnp}$. One can see from Eqs. (44) and (48b) that $i$th face lies in a plane of constant coordinate $\hat{z}=h_{i}$. One can therefore take the factor $\hat{z}^{p}$ outside the integral sign in Eqs. (49a), (49b), what gives

[TABLE]

Applying the planar variant of Gauss’s theorem to the vector $(G_{\hat{x}},G_{\hat{y}})$ defined as

[TABLE]

one can reduce the surface integral $\hat{\Phi}^{(i)}_{mn0}$ to the quantity $\hat{\Omega}^{(i)}_{mn0}$ and line integrals $\hat{L}^{(k)}_{mn0}$. The details of the calculation can be found in Appendix C, and the following result was obtained:

[TABLE]

Here and below in this section, the index $k$ runs over edges adjacent to $i$th face.

The next step is decreasing the index $n$ at $\hat{\Omega}^{(i)}_{mn0}$. If $n\geq 2$, one can apply the following relation derived in Appendix D,

[TABLE]

This recursive relation makes it possible to lower the second index at $\Omega$ down to $n=1$ or $n=0$.

In the case of $n=1$, the expression for $\hat{\Omega}^{(i)}_{m10}$ can be simplified by direct integration over $\hat{y}$ (see details in Appendix E), yielding

[TABLE]

The last case to consider is that of $n=0$. If $m\geq 2$, one can reduce $m$ by two as follows:

[TABLE]

This relation is derived in Appendix F. Applying Eq. (58) repeatedly, one arrives finally at $\hat{\Omega}^{(i)}_{100}$ or $\hat{\Omega}^{(i)}_{000}\equiv\Omega^{(i)}$ plus line integrals.

Finally, the value of $\hat{\Omega}^{(i)}_{100}$ can be found in full analogy with Eq. (57):

[TABLE]

Thus, we have seen in this section that any of quantities $\varphi_{mnp}$, $\Phi^{(i)}_{mnp}$, $\Omega^{(i)}_{mnp}$ can be reduced to line integrals over inclusion edges and/or to solid angles $\Omega^{(i)}$. A simple but not so trivial example of such reducing will be considered in Subsection V.1. Then, in Section VI, we will know how to evaluate the line integrals.

V.1 Example: evaluation of $\Omega^{(i)}_{001}$

Here we will apply the above-formulated considerations to finding an analytical expression for the quantity $\Omega^{(i)}_{001}$. The first step is choosing the “tilted” frame associated with $i$th face of the polyhedron, i. e., choosing three mutually orthogonal unit vectors $\hat{\mathbf{e}}_{x}$, $\hat{\mathbf{e}}_{y}$, $\hat{\mathbf{e}}_{z}$ along the axes of the “tilted” frame. The vector $\hat{\mathbf{e}}_{z}$ should be orthogonal to $i$th face:

[TABLE]

We choose the vector $\hat{\mathbf{e}}_{x}$ in the plane spanned onto $\mathbf{e}_{z}$ and $\hat{\mathbf{e}}_{z}$ (see Fig. 5):

[TABLE]

where $\theta$ is the angle between $\mathbf{e}_{z}$ and $\hat{\mathbf{e}}_{z}$:

[TABLE]

And the vector $\hat{\mathbf{e}}_{y}$ is to be orthogonal to $\hat{\mathbf{e}}_{x}$ and $\hat{\mathbf{e}}_{z}$.

Then, according to Eq. (50b), $\Omega^{(i)}_{001}$ is a linear combination of the “tilted” quantities $\hat{\Omega}^{(i)}_{100}$, $\hat{\Omega}^{(i)}_{010}$ and $\hat{\Omega}^{(i)}_{001}$, taken with their respective coefficients $T^{100}_{001}$, $T^{010}_{001}$ and $T^{001}_{001}$. In order to find these coefficients from Eq. (51), one should represent the value of $z-Z$ via the “tilted” coordinates $\hat{x},\hat{y},\hat{z}$:

[TABLE]

whence, $T^{100}_{001}=\sin\theta$, $T^{010}_{001}=0$, and $T^{001}_{001}=\cos\theta$. Therefore

[TABLE]

Then one can evaluate the quantity $\hat{\Omega}^{(i)}_{100}$ using Eq. (59), where scalar products $\mathbf{b}_{ik}\cdot\hat{\mathbf{e}}_{x}$ can be found from Eq. (61) taking into account that the vector $\mathbf{b}_{ik}$ is orthogonal to $\hat{\mathbf{e}}_{z}\equiv\mathbf{n}_{i}$:

[TABLE]

Consequently, Eq. (59) takes the following form:

[TABLE]

And the quantity $\hat{\Omega}^{(i)}_{001}$ is calculated according to Eq. (53):

[TABLE]

Finally, collecting Eqs. (64), (66), (67), one can get the following answer:

[TABLE]

This is an analytical representation of $\Omega^{(i)}_{001}$ in a closed form, because, as mentioned in Section IV, the primitives $L^{(k)}$ and $\Omega^{(i)}$ can be expressed through elementary functions.

VI Evaluation of line integrals $L^{(k)}_{mnp}$

Let $\mathrm{A}_{k}$ and $\mathrm{B}_{k}$ be two end points of $k$th edge (we denote their radius vectors as $\mathbf{r}_{A}^{(k)}$ and $\mathbf{r}_{B}^{(k)}$), and $\mathbf{l}_{k}$ be the unit vector directed from $\mathrm{A}_{k}$ to $\mathrm{B}_{k}$. To find analytical expressions for the line integrals $L^{(k)}_{mnp}$, it is convenient to reduce them to more simple integrals $\mathcal{L}^{(k)}_{t}$:

[TABLE]

where $t$ is a non-negative integer number, $\xi=\mathbf{l}_{k}\cdot(\mathbf{r}-\mathbf{R})$ is a coordinate along $k$th edge, $\rho_{k}$ is the distance from the point $\mathbf{R}$ to the line of $k$th edge:

[TABLE]

$\xi_{k1}$ and $\xi_{k2}$ are the values of $\xi$ at $\mathrm{A}_{k}$ and $\mathrm{B}_{k}$:

[TABLE]

In order to perform such reduction, let us express the position vector $\mathbf{r}$ of a point on the edge through the coordinate $\xi$:

[TABLE]

where the point $\mathbf{r}^{(k)}_{0}=\left(x^{(k)}_{0},y^{(k)}_{0},z^{(k)}_{0}\right)$ is the projection of the point $\mathbf{R}$ to $k$th edge:

[TABLE]

Using Eq. (72), one can represent the factor $(x-X)^{m}(y-Y)^{n}(z-Z)^{p}$ that appears in the definition of $L^{(k)}_{mnp}$ as follows:

[TABLE]

The right hand side of Eq. (74) is some polynomial on $\xi$:

[TABLE]

where each coefficient $c_{mnp,t}(\mathbf{R})$ is a polynomial of $X,Y,Z$ of degree $m+n+p-t$. Substituting this series decomposition into the line integral $L^{(k)}_{mnp}$, Eq. (29), and taking into account that $|\mathbf{r}-\mathbf{R}|=\sqrt{\rho_{k}^{2}+\xi^{2}}$ on the edge, one can easily see that

[TABLE]

A similar representation takes place for the “tilted” line integrals $\hat{L}^{(k)}_{mnp}$ defined by Eq. (49c):

[TABLE]

The coefficients $\hat{c}_{mnp,t}(\mathbf{R})$ are to be obtained from the identity

[TABLE]

where

[TABLE]

and the same for $y$- and $z$-components.

The last question is how to find analytical expressions for the integrals $\mathcal{L}^{(k)}_{0},\mathcal{L}^{(k)}_{1},\ldots$ As explained in Appendix G, integration of Eq. (69) by parts provides the following recursive formula for $t>1$:

[TABLE]

and a simple answer for $t=1$:

[TABLE]

Eqs. (80) and (81) provide a possibility to evaluate the quantities $\mathcal{L}^{(k)}_{t}$ for integer $t\geq 0$ in terms of $\mathcal{L}^{(k)}_{0}$ and coordinates of the ends $\mathrm{A}_{k}$ and $\mathrm{B}_{k}$ of the edge. In its turn, $\mathcal{L}^{(k)}_{0}$ is the same as $L^{(k)}$ and therefore has the analytical representation according to Eq. (43). So, the problem of evaluating of line integrals along inclusion edges is solved.

VII The main result

Here we summarize the results of Sections II–V and present the analytical expression for the strain distribution that solves the problem posed in the Introduction.

At first, let us recall the steps carried out in the above text. In Section II, the strain tensor $\varepsilon_{\alpha\beta}(\mathbf{R})$ was expressed in terms of second derivatives of the electrostatic potential $\varphi(\mathbf{R})$. In Section III, these second derivatives, as well as the first derivatives and the potential itself, were represented as combinations of the quantities $\varphi_{mnp}$, $\varphi^{(mnp)}_{,\alpha}$, $\varphi^{(mnp)}_{,\alpha\beta}$, which correspond to the charge distribution $(x-X)^{m}(y-Y)^{n}(z-Z)^{p}\,\chi(\mathbf{r})$. These quantities, in their turn, were reduced in Section IV to four kinds of “primitives”: volume integrals $\varphi_{mnp}$, surface integrals $\Phi^{(i)}_{mnp}$ and $\Omega^{(i)}_{mnp}$, and line integrals $L^{(k)}_{mnp}$. In Section V, the volume and surface integrals were expressed through the solid angles $\Omega^{(i)}$ and line integrals $\hat{L}^{(k)}_{mnp}$. Finally, in SectionVI all the line integrals $L^{(k)}_{mnp}$ and $\hat{L}^{(k)}_{mnp}$ were reduced to the simplest line integrals $L^{(k)}$, whose analytical representations are given by Eq. (43), and to the distances from the point $\mathbf{R}$ (where the strain tensor is to be found) to vertices of the polyhedron.

Each of these steps represents itself a linear transformation whose coefficients are either constants (such as components of unit vectors $\mathbf{n}_{i}$ in Eq. (30)) or polynomials of $\mathbf{R}$ (for example, the factor $h_{i}=(\mathbf{r}_{i}-\mathbf{R})\cdot\mathbf{n}_{i}$ in Eq. (47)).

As a result, the strain tensor $\varepsilon_{\alpha\beta}(\mathbf{R})$ appears to be a linear combination of three kinds of analytical functions:

• •

surface integrals (potentials of the inclusion faces covered by uniform double layers) $\Omega^{(i)}(\mathbf{R})\equiv\Omega^{(i)}_{000}(\mathbf{R})$ defined according to Eq. (28) and having a geometrical meaning of solid angles;

• •

line integrals (potentials of uniformly charged edges of the inclusion) $L^{(k)}(\mathbf{R})\equiv L^{(k)}_{000}(\mathbf{R})$ that have the analytical form according to Eq. (43);

• •

distances $|\mathbf{R}-\mathbf{r}_{s}|$ between the point $\mathbf{R}$ and vertices $\mathbf{r}_{s}$ of the inclusion (we will use the index $s$ to number the vertices).

The same is true for the potential $\varphi(\mathbf{R})$ and its first and second derivatives on $\mathbf{R}$. The coefficients at quantities $\Omega^{(i)}$, $L^{(k)}$ and $|\mathbf{R}-\mathbf{r}_{s}|$ come from substituting polynomial expressions one into the other, and therefore should be polynomials of $\mathbf{R}$.

The facts described in this Section constitute the main result of the present study. In Subsection VII.1 we represent them in a more formal way, writing down the analytical expressions that define the distributions of the strain, the potential, and derivatives of the potential, up to some polynomial coefficients. Then, in Subsection VII.2 we will show how to find these polynomial coefficients.

VII.1 The main result: Analytical expressions for the strain distribution, the potential, and its derivatives

Based on the above-presented considerations, one can formulate a universal expression for the strain distribution $\varepsilon_{\alpha\beta}(\mathbf{R})$ induced by a polyhedral inclusion with a lattice misfit $\varepsilon_{0}(\mathbf{R})$ that depends on coordinates polynomially:

[TABLE]

Here summations take place over all faces (index $i$), edges (index $k$), and vertices (index $s$) of the polyhedron. $\Omega^{(i)}(\mathbf{R})$ and $L^{(k)}(\mathbf{R})$ are analytical functions described in Section IV: $\Omega^{(i)}$ has a geometrical meaning of a solid angle, and $L^{(k)}$ is defined by Eq. (43). $\mathbf{r}_{s}$ is a position vector of $s$th vertex. $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}(\mathbf{R})$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}(\mathbf{R})$, $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}(\mathbf{R})$ are some polynomial functions of $\mathbf{R}$ that are defined by orientations of faces and edges and by coordinates of vertices. The degrees of polynomials $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$ and $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$ are equal to $N$, i. e., to the degree of the polynomial $f(\mathbf{r})$ that determines the distribution of lattice misfit, whereas degrees of $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$ are equal to $N-1$. (if $N=0$, that corresponds to a constant lattice misfit, then $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}=0$ at any $\mathbf{R}$.) A recipe for calculation of the polynomial coefficients of $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$ and $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$ will be described in Subsection VII.2. Of course, Eq. (82) is valid only if the assumptions listed in the Introduction are fulfilled.

In a fully analogous manner, one can write down a similar expression for the electrostatic potential $\varphi$ [a solution of Poisson’s equation (7)] produced by a charged polyhedron with a non-uniform charge distribution $\varepsilon(\mathbf{r})$, which is a polynomial function of coordinates:

[TABLE]

where polynomials $\mathcal{A}^{(i)}_{\varphi}(\mathbf{R})$ and $\mathcal{B}^{(k)}_{\varphi}(\mathbf{R})$ have the degree $N+2$, and polynomials $\mathcal{C}^{(s)}_{\varphi}(\mathbf{R})$ have the degree $N+1$.

The first and second derivatives of the potential $\varphi(\mathbf{R})$ also have similar representations:

[TABLE]

and

[TABLE]

where $\mathcal{A}^{(i)}_{\varphi\alpha}(\mathbf{R})$ and $\mathcal{B}^{(k)}_{\varphi\alpha}(\mathbf{R})$ are polynomials of degree $N+1$; $\mathcal{C}^{(s)}_{\varphi\alpha}(\mathbf{R})$, $\mathcal{A}^{(i)}_{\varphi\alpha\beta}(\mathbf{R})$ and $\mathcal{B}^{(k)}_{\varphi\alpha\beta}(\mathbf{R})$ are polynomials of degree $N$; and $\mathcal{C}^{(s)}_{\varphi\alpha\beta}(\mathbf{R})$ are polynomials of degree $N-1$. The coefficients of all these polynomials are to be calculated by the method of Subsection VII.2.

VII.2 The main result: Algorithm of finding the polynomial coefficients

Equations (82)–(85) determine the sought-for functions $\varepsilon_{\alpha\beta}(\mathbf{R})$, $\varphi(\mathbf{R})$, $\partial\varphi(\mathbf{R})/\partial R_{\alpha}$ and $\partial^{2}\varphi(\mathbf{R})/\partial R_{\alpha}\partial R_{\beta}$ up to the coefficients of the corresponding polynomials $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$. Here we represent the algorithm (a sequence of steps) for calculating these coefficients. The algorithm is just a summary of equations obtained in previous sections.

On the input of the algorithm there are: the polynomial function $f(x,y,z)$, which defines the distribution of eigenstrain inside the inclusion, and the geometry of the inclusion—unit vectors $\mathbf{n}_{i}$, $\mathbf{b}_{ik}$, $\mathbf{l}_{k}$ and coordinates of vertices.

On the output there are polynomials $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}(\mathbf{R})$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}(\mathbf{R})$, $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}(\mathbf{R})$ that determine the distribution of elastic strain $\varepsilon_{\alpha\beta}(\mathbf{R})$ according to Eq. (82). For finding the polynomials related to the potential and its derivatives, one should modify Steps 1–3, as explained in the end of this Subsection.

The algorithm consists of eleven steps listed below. Each step contains only arithmetic operations on polynomials. Therefore, it is easy to implement the algorithm in a program code.

Step 1. Express strain components $\varepsilon_{\alpha\beta}$ via second derivatives $\partial^{2}\varphi/\partial R_{\alpha}\partial R_{\beta}$ using Eq. (12); express the function $\varepsilon_{0}$ that appears in the right-hand side of Eq. (12) via solid angles $\Omega^{(i)}$ using Eq. (42).

Step 2. Find the coefficients $\tilde{C}_{mnp}(\mathbf{R})$ from Eq. (19), and represent the second derivatives $\partial^{2}\varphi/\partial R_{\alpha}\partial R_{\beta}$ as combinations of $\varphi^{(mnp)}_{,\alpha\beta}$ using Eq. (22).

Step 3. Express each of the quantities $\varphi^{(mnp)}_{,\alpha\beta}$ through integrals $\varphi_{m^{\prime}n^{\prime}p^{\prime}}$, $\Phi^{(i)}_{m^{\prime}n^{\prime}p^{\prime}}$, $\Omega^{(i)}_{m^{\prime}n^{\prime}p^{\prime}}$ and $L^{(k)}_{m^{\prime}n^{\prime}p^{\prime}}$ with different $m^{\prime},n^{\prime},p^{\prime}$ according to Eqs. (34), (37), (38).

Step 4. Reduce each of the volume integrals $\varphi_{mnp}$ to a combination of surface integrals $\Phi^{(i)}_{mnp}$ using Eq. (47). (After this step the strain tensor, as a function of coordinates, is represented as a linear combination of functions $\Phi^{(i)}_{mnp}$, $\Omega^{(i)}_{mnp}$, $L^{(k)}_{mnp}$ with polynomial coefficients.)

Steps 5–9 should be performed for each face $i$ of the inclusion:

Step 5. Choose a “tilted” coordinate system $\hat{x},\hat{y},\hat{z}$ with axis $\hat{z}$ along the normal $\mathbf{n}_{i}$ to $i$th face. Then, transform the integrals $\Phi^{(i)}_{mnp}$ and $\Omega^{(i)}_{mnp}$ into their “tilted” versions $\hat{\Phi}^{(i)}_{mnp}$ and $\hat{\Omega}^{(i)}_{mnp}$ using Eqs. (50a) and (50b). The coefficients of transformation $T^{m^{\prime}n^{\prime}p^{\prime}}_{mnp}$ are to be found from the identity (51).

Step 6. Decrease the index $p$ at $\hat{\Phi}^{(i)}_{mnp}$ and $\hat{\Omega}^{(i)}_{mnp}$ down to zero by means of Eq. (53).

Step 7. Express the integrals $\hat{\Phi}^{(i)}_{mn0}$ via $\hat{\Omega}^{(i)}_{mn0}$ and $\hat{L}^{(k)}_{mn0}$ using Eq. (55). (After this step, the strain distribution is expressed in terms of surface integrals $\hat{\Omega}^{(i)}_{mn0}$ and line integrals $L^{(k)}_{mnp}$, $\hat{L}^{(k)}_{mn0}$.)

Step 8. For each $\hat{\Omega}^{(i)}_{mn0}$, decrease the index $n$ recursively down to $n=1$ or $n=0$ by means of Eq. (56). Then, if $n=1$, express the quantity $\hat{\Omega}^{(i)}_{m10}$ via $\hat{L}^{(k)}_{m00}$ using Eq. (57).

Step 9. For each $\hat{\Omega}^{(i)}_{m00}$, decrease the index $m$ recursively down to $m=1$ or $m=0$ by means of Eq. (58). Then, if $m=1$, express the quantity $\hat{\Omega}^{(i)}_{100}$ via $L^{(k)}$ using Eq. (59). (After this step, the strain distribution appears to be a combination of solid angles $\Omega^{(i)}$ and line integrals $L^{(k)}_{mnp}$, $\hat{L}^{(k)}_{mn0}$.)

Step 10. For each edge $k$ of the inclusion, decompose each line integral $L^{(k)}_{mnp}$ into a series $c_{mnp,0}\mathcal{L}^{(k)}_{0}+c_{mnp,1}\mathcal{L}^{(k)}_{1}+\ldots$ according to Eq. (76), where coefficients $c_{mnp,0}(\mathbf{R}),c_{mnp,1}(\mathbf{R}),\ldots$ are to be found from Eq. (75). Do the same with “tilted” line integrals $\hat{L}^{(k)}_{mn0}$ by means of Eqs. (77), (78).

Step 11. For each edge $k$ and each line integral $\mathcal{L}^{(k)}_{t}$, decrease the index $t$ recursively down to $t=1$ or $t=0$ by means of Eq. (80). Then, if $t=1$, evaluate the integral $\mathcal{L}^{(k)}_{1}$ using Eq. (81).

After performing Steps 1–11, one obtains the strain distribution $\varepsilon_{\alpha\beta}(\mathbf{R})$ in the form of a combination of solid angles $\Omega^{(i)}(\mathbf{R})$, logarithmic expressions $L^{(k)}(\mathbf{R})$, and distances from $\mathbf{R}$ to inclusion vertices, i. e., in the form of Eq. (82) with the known polynomial coefficients $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}(\mathbf{R})$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}(\mathbf{R})$, $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}(\mathbf{R})$. Hence, these steps actually constitute an algorithm for calculating the polynomials $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$.

In Appendix J, one can find an example of implementation of this algorithm in the form of an informal “pseudocode”.

In a slightly modified form, this algorithm can be also used to calculate the potential $\varphi(\mathbf{R})$, its derivatives and its second derivatives. For that, it is enough to omit Step 1 and make obvious modifications of Steps 2 and 3. Namely, for finding an expression for the potential $\varphi(\mathbf{R})$, one can use Eq. (20) instead of Eq. (22), and omit Step 3. For finding the expressions for the first derivatives $\partial\varphi(\mathbf{R})/\partial R_{\alpha}$, one can use Eq. (21) instead of Eq. (22) at Step 2, and Eqs. (30), (32) instead of Eqs. (34), (37), (38) at Step 3.

VIII Examples

This Section has two aims: (i) to illustrate the results of Section VII on concrete examples, and (ii) to provide the evidence that our results are correct.

VIII.1 Constant lattice misfit

For the simplest example, let us suppose that the lattice parameter is the same in any point of the inclusion. Therefore

[TABLE]

and the lattice misfit $\varepsilon_{0}$ is equal to $f$ inside the inclusion, and to 0 outside. Applying the algorithm of Section VII to this case, one can see that the answer for the strain distribution in the form (82) appears already after Step 3:

[TABLE]

where values of $\Lambda$ and $\lambda^{(k)}_{\alpha\beta}$ are defined by Eqs. (10) and (33). A similar result was reported in Ref. Nenashev and Dvurechenskii, 2010.

A comparison between Eq. (82) and Eq. (87) provides the values of the coefficients $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$ for the case of a constant lattice misfit inside the inclusion:

[TABLE]

For better understanding of the way leading to this answer, one can visualize the algorithm in the form of a tree, as shown in the left part of Fig. 6. The sought-for function $\varepsilon_{\alpha\beta}$ is placed at the top of this graph. At Step 1, this function acquires the following representation:

[TABLE]

that is symbolized in the graph by arrows connecting $\varepsilon_{\alpha\beta}$ to $\partial^{2}\varphi/\partial R_{\alpha}\partial R_{\beta}$ and to $\Omega^{(i)}$. A value near an arrow means the factor, with which the expression on the head of the arrow contributes to the expression on the tail. Then, Step 2 represents $\partial^{2}\varphi/\partial R_{\alpha}\partial R_{\beta}$ as follows:

[TABLE]

This connection is shown in the graph by the arrow with mark “$f$”. Finally, Step 3 provides the expression for $\varphi^{(000)}_{,\alpha\beta}$:

[TABLE]

which is depicted by two lower arrows in the graph.

The formulas for the coefficients $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$ can be easily constructed on the basis of this graph. Each appearance of $\Omega^{(i)}$ or $L^{(k)}$ in the graph provides the contribution to $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$ or $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$, correspondingly. The value of the contribution is equal to the product of the quantities on arrows along the way from the root of the tree to $\Omega^{(i)}$ or $L^{(k)}$. For example, the way to $L^{(k)}$ consists of three arrows with values $-\Lambda$, $f$ and $\lambda^{(k)}_{\alpha\beta}$ on them; hence, their product $-\Lambda f\lambda^{(k)}_{\alpha\beta}$ appears in Eq. (88b).

In a similar manner, one can calculate the potential of a uniformly charged polyhedral body $\varphi(\mathbf{R})$, and its first and second derivatives. Since $f$ is a constant ($N=0$), only primitives with $m=n=p=0$ will appear in the course of the algorithm. For this reason, one can omit Steps 5, 6, 8–11 (as well as Step 1, which is needed only for calculating the strain). The result for the potential is

[TABLE]

where $i_{1}$ and $i_{2}$ are numbers of the two faces that intersect at $k$th edge; the linear functions $h_{i}(\mathbf{R})$ and $B_{ik}(\mathbf{R})$ are defined by Eqs. (44) and (45), correspondingly. Expressions (92) are to be inserted into the general formula (83).

For the first derivatives of the potential, the answer is:

[TABLE]

and there is the answer for the second derivatives of the potential:

[TABLE]

In Fig. 6 one can see also graphical representations for the potential $\varphi$ (middle part), its derivative $\partial\varphi/\partial R_{\alpha}$ (right part), and its second derivative $\partial^{2}\varphi/\partial R_{\alpha}\partial R_{\beta}$ (the part of the left graph, beginning from $\partial^{2}\varphi/\partial R_{\alpha}\partial R_{\beta}$).

Our expressions (92)–(94) for the potential of a homogeneous polyhedron and its first and second derivatives are the same as the results by Werner and Scheeres.Werner and Scheeres (1996) One can consider the coincidence of our results for $f=\text{const}$ with the published ones as a simplest test on the correctness of our approach. According to the present Subsection, this test is successfully passed.

Also we have performed a comparison with analytical calculations made by GlasGlas (2001) for the case of a truncated pyramid with a constant misfit strain. The results of calculations shown in Fig. 3b of Ref. Glas, 2001 are fully reproduced by our method. For more details, see Appendix H.

VIII.2 Constant, vertically directed gradient of the lattice misfit

Now let us examine the case of linear dependence of the eigenstrain $\varepsilon_{0}$ on the coordinates inside the inclusion. For simplicity, we suppose in this Subsection that the eigenstrain depends on coordinate $z$ only:

[TABLE]

where $a$ (the gradient of the eigenstrain) and $b$ are constants. For such a function $f(\mathbf{r})$, the algorithm of Section VII is shown in a graphical form in Fig. 7. The meaning of arrows and expressions in the nodes of the graph is explained above, see the caption to Fig. 6. In the part of the tree below $\Omega^{(i)}_{001}$, we use the results of Subsection V.1. For the arrow connecting $\mathcal{L}^{(k)}_{1}$ to $|\mathbf{R}-\mathbf{r}_{s}|$, the factor is equal to $+1$ if the vector $\mathbf{l}_{k}$ is directed from the opposite vertex of $k$th edge to $s$th vertex, and $-1$ otherwise.

The expressions for the polynomials $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$, $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$ can be derived from the graph in Fig. 7 by the method described in Subsection VIII.1. The only new feature is the existence of a node $|\mathbf{R}-\mathbf{r}_{s}|$ that provides a contribution to the polynomial $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$, which is equal to the product of the quantities near the arrows on the way from $\varepsilon_{\alpha\beta}$ to $|\mathbf{R}-\mathbf{r}_{s}|$. As a result, one can get the following answers for $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$ and $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$:

[TABLE]

where index $k$ runs over the edges that enter $s$th vertex. The sign in the factor $(\pm l_{kz})$ is chosen to be plus if the vector $\mathbf{l}_{k}$ is directed from the opposite vertex of $k$th edge to $s$th vertex, and minus otherwise. Constants $\Lambda$, $\lambda^{(k)}_{\alpha\beta}$, and linear functions $h_{i}(\mathbf{R})$, $B_{ik}(\mathbf{R})$, $z^{(k)}_{0}(\mathbf{R})$ are defined by Eqs. (10), (33), (44), (45), (73), correspondingly.

In the next Subsection, these equations will be generalized to the case of an arbitrarily directed gradient of the lattice misfit.

VIII.3 Constant gradient of the lattice misfit: general case

Consider the general case of a constant gradient of the lattice misfit inside the polyhedral inclusion:

[TABLE]

where a vector $\mathbf{a}$ and a scalar $b$ are constants. It is easy to generalize formulas (96a)–(96c) to this case. For this purpose, it is enough to represent them in a coordinate-independent form, i. e., to replace the product $aZ$ by $\mathbf{a}\cdot\mathbf{R}$, the product $an_{iz}$ by $\mathbf{a}\cdot\mathbf{n}_{i}$, etc. This gives rise to the following result:

[TABLE]

In order to get the strain distribution $\varepsilon_{\alpha\beta}(\mathbf{R})$, one should insert these expressions into the universal equation (82).

By the same method, one can get the analytical formulas for the potential $\varphi$ (as well as its first and second derivatives) of a non-uniformly charged polyhedral body with a charge density $\mathbf{a}\cdot\mathbf{r}+b$. The polynomials $\mathcal{A}^{(i)}_{\varphi}$, $\mathcal{B}^{(k)}_{\varphi}$ and $\mathcal{C}^{(s)}_{\varphi}$ related to the potential are expressed as follows:

[TABLE]

Eqs. (99a)–(99c) are to be inserted into Eq. (83) in order to get the potential $\varphi(\mathbf{R})$.

Analytical expressions for the gravitational potential of a polyhedral body having a linearly varying density distribution were reported by many authors.Pohánka (1998); Holstein (2003); Hamayun, Prutkin, and Tenzer (2009); D’Urso (2014); Conway (2015) However these expressions were not represented in terms of solid angles. For this reason, it seems to be extremely difficult to compare them directly with our answer [a combination of Eqs. (83) and (99)] at the level of formulas. To ensure that the results of this Subsection are correct, we performed numerical tests, considering a pyramid with a square base as an example of an inclusion (see Fig. 8). The following statements were checked numerically:

• •

$\varphi(x,y,z)$ is a continuous function, and its first derivatives are continuous;

• •

$\varphi(x,y,z)$ obeys Poisson’s equation $\Delta\varphi=-4\pi\varepsilon_{0}$, where $\varepsilon_{0}(\mathbf{r})=\mathbf{a}\cdot\mathbf{r}+b$ inside the pyramid, and $\varepsilon_{0}=0$ outside the pyramid;

• •

the functions $\varphi(x,y,z)$ and $\varepsilon_{\alpha\beta}(x,y,z)$ are connected to each other by Eq. (12).

The potential $\varphi$ and the strain $\varepsilon_{\alpha\beta}$ were calculated by Eqs. (83), (99) and Eqs. (82), (98), correspondingly. The pyramid height $h$, the length $l$ of the edge of the pyramid base, the gradient $\mathbf{a}$ of the lattice misfit, and the free term $b$ were varied independently. An example of the calculated potential $\varphi$ and $xx$-component of the strain is shown in Fig. 8 for the case of $h=5$ nm, $l=10$ nm, $a_{x}=a_{y}=0$, and $a_{z}=-b/2h$. (Such a choice of $\mathbf{a}$ means that the lattice misfit linearly decreases from the base to the apex of the pyramid, so that at the apex it is twice smaller than at the base.) For this case, as well as for many other combinations of parameters $h$, $l$, $\mathbf{a}$ and $b$, the numerical test has been passed successfully.

In order to provide more evidences of validity of analytical results, we have compared the distribution of $\varepsilon_{xx}$ calculated analytically and shown in Fig. 8 with analogous distributions calculated numerically using the finite element method. This comparison is presented in Appendix I. One can see that for large enough size of the “box”, in which the inclusion is incorporated, the numerical results perfectly agree with the analytical ones.

VIII.4 Sinusoidal profile of the lattice misfit

As an example of a more complicated distribution of the lattice misfit $\varepsilon_{0}(x,y,z)$, let us consider the sinusoidal profile inside the inclusion:

[TABLE]

so that $f(\mathbf{r})=A\sin(kx)$, where $A$ is a constant. In order to apply the results of the present paper, one has to approximate the function $f(\mathbf{r})$ by a polynomial. As such an approximation, we choose first five terms of the Fourier series,

[TABLE]

that provides accuracy not worse than one percent of $A$ within one period, i. e., for $|kx|\leq\pi$.

In Fig. 9, an example of the potential $\varphi(x,y,z)$, its derivatives $\partial\varphi/\partial x$ and $\partial^{2}\varphi/\partial x^{2}$, and the $xx$-component of the strain tensor is shown for the same pyramidal inclusion as in Fig. 8. The height $h$ and the lateral size $l$ of the pyramid are 5 nm and 10 nm, correspondingly. The parameter $k$ is chosen to be $k=2\pi/l$, so that the pyramid size $l$ corresponds to the period of the sinusoid.

The values of the strain tensor, the potential, and its first and second derivatives were calculated by Eqs. (82)–(85). The coefficients of polynomials $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$, which contribute to these equations, were calculated numerically, directly following the algorithm of Subsection VII.2. Then, the same numerical test as for a linearly varying eigenstrain (see the previous Subsection) was performed and successfully passed. In addition, it was checked that the results for the potential and its first and second derivatives are consistent with each other.

Such a numerical test with a degree $N$ of the polynomial $f(\mathbf{r})$, larger than 1, is important for the following reason. Equations (56), (58) and (80) are used at Steps 8, 9, 11 of the algorithm only when index $m$, $n$ or $t$ becomes larger than 1, that occurs only if $N\geq 2$. Hence, to check the algorithm, it is necessary to consider at least quadratic ($N=2$) coordinate dependence of the eigenstrain. In the example of this Subsection, $N$ is equal to 9, that provides a good test of the algorithm as a whole.

IX Concluding remarks

Two problems, closely related to each other, are considered in this paper. The first one is the search for an analytical solution (in a closed form) for the spatial distribution of the elastic strain induced by a lattice-mismatched inclusion inserted in an infinite matrix, in the case of a polyhedral shape of the inclusion. The second problem is finding an analytical solution for an electrostatic of gravitational potential created by a polyhedral body, and for the derivatives of the potential. It is well known that the strain tensor is expressed via the second derivatives of the potential, see Eq. (12), assuming that the inclusion and the matrix are elastically isotropic and have the same elastic constants. Under these assumptions, we obtained the universal analytical expression describing the strain distribution, Eq. (82), for a polynomial coordinate dependence of the lattice misfit within the inclusion. Similar analytical formulas, Eqs. (83)–(85), were derived for the potential and its first and second derivatives in the case of a polynomial distribution of the charge/mass density over a polyhedral-shaped body.

The results of the present work allow one to understand the structure of the known analytical solutions for different particular cases, e. g. strain distributions due to homogeneous polyhedral inclusions, potentials of polyhedral bodies with a constant density or a constant gradient of the density. Typically, the solutions found in the literature are combinations of logarithmic terms and terms containing inverse trigonometric functions, mostly arctangents. Our expressions (83)–(85) clarify that logarithmic terms arise from potentials of uniformly charged edges $L^{(k)}$, and arctangent terms—from solid angles $\Omega^{(i)}$.

Finally, it is worth highlighting several ways of generalizing the results of the present study. The strain distribution for an inclusion in a half-space can be found by the “mirror image” method.Davies (2003) The effects of anisotropy can be taken into account either by series expansion of the elastic Green’s function,Faux and Pearson (2000); Faux and Christmas (2005) or by the method of stretching the coordinates.Kuvshinov (2008); Nenashev, Koshkarev, and Dvurechenskii (2013) In the case of different elastic constants of the matrix and the inclusion (an inhomogeneous inclusion), one can apply the perturbation method based on the concept of effective inclusion.Chu et al. (2011); Ma and Korsunsky (2014)

Acknowledgements.

This work is funded by RFBR (grants #16-29-14031 and #16-02-00397, development of analytical expressions) and by State Programme (Grant No. 0306-2016-0015, creating the computer program for calculation of the strain).

Appendix A Integrals containing the derivatives of characteristic function $\chi(\mathbf{r})$

This Appendix is devoted to reducing the integrals $\varphi^{(mnp)}_{,\alpha}$ and $\varphi^{(mnp)}_{,\alpha\beta}$ defined by Eqs. (24) and (25) to the primitive introduced in Section IV.

Let us define a function $f(\mathbf{r})$ as

[TABLE]

then

[TABLE]

and

[TABLE]

First we will focus on the integral $\varphi^{(mnp)}_{,\alpha}$. Let us open the brackets in the expression $\partial(f\chi)/\partial r_{\alpha}$:

[TABLE]

Substituting this expansion into Eq. (103), one can see that $\varphi^{(mnp)}_{,\alpha}$ is equal to the sum of two integrals, the first of which being the quantity $\varphi_{(mnp),\alpha}$ defined by Eq. (31). Therefore

[TABLE]

In the remaining integral, the derivative $\partial\chi/\partial r_{\alpha}$ differs from zero only on the surface of the polyhedron; thus, it is in fact a surface integral. To be more precise, let us assume that the polyhedron is convex (generalization to non-convex polyhedra does not lead to any difficulties), and define the characteristic function $\chi(\mathbf{r})$ as follows:

[TABLE]

where $\theta$ denotes the Heaviside step function, $\mathbf{n}_{i}$ is the outward normal to $i$th face, and $\mathbf{r}_{i}$ is a position vector of some (arbitrarily chosen) point on $i$th face. It is easy to check that Eqs. (5) and (107) define the same function $\chi(\mathbf{r})$. Differentiation of each factor in Eq. (107) is straightforward:

[TABLE]

where $\delta$ denotes the Dirac delta function. Using Eq. (108), one can represent the derivative $\partial\chi/\partial r_{\alpha}$ as follows:

[TABLE]

where functions $\psi_{i}$ are defined as

[TABLE]

One can easily see that the function $\psi_{i}$ vanishes everywhere except $i$th face of the polyhedron’s surface, where it has a $\delta$-like singularity. Hence, any volume integral containing the function $\psi_{i}$ is in fact a surface integral over $i$th face. In particular,

[TABLE]

Using Eqs. (109) and (111), one can reduce the integral in the right-hand side of Eq. (106) to surface integrals $\Phi^{(i)}_{mnp}$:

[TABLE]

Substitution of this representation to Eq. (106) gives rise to Eq. (30).

Hence, we have proved Eq. (30), i. e., have found out a representation of the integral $\varphi^{(mnp)}_{,\alpha}$ in terms of primitives.

Now let us consider the integral $\varphi^{(mnp)}_{,\alpha\beta}$. The second derivative in its definition, Eq. (104), can be expanded as follows:

[TABLE]

After substituting this expansion to Eq. (104), one gets four integrals corresponding to four terms of the right-hand side of Eq. (113). The first integral is equal to the quantity $\varphi_{(mnp),\alpha\beta}$ introduced in Eq. (35):

[TABLE]

The second integral, which contains the derivative $\partial\chi/\partial r_{\alpha}$, can be treated as above, via Eqs. (109) and (111):

[TABLE]

where $\Phi^{(i)}_{(mnp),\beta}$ is the quantity defined by Eq. (36). The third integral is the same as the second one, up to the index permutation $\alpha\leftrightarrow\beta$. Therefore the expression for $\varphi^{(mnp)}_{,\alpha\beta}$ takes the following form:

[TABLE]

In order to go further, one needs a method of dealing with integrals containing the second derivative $\partial^{2}\chi/\partial r_{\alpha}\partial r_{\beta}$. Using Eq. (109), this second derivative can be represented as

[TABLE]

Then, it is more convenient to define the function $\psi_{i}(\mathbf{r})$ in an alternative way:

[TABLE]

where index $k$ runs over the edges that surround $i$th face, and unit vectors $\mathbf{b}_{ik}$ are directed out of $i$th face perpendicularly to their corresponding edges, as shown in Fig. 2. One can easily figure out that both equations (110) and (118) define the same function, which has delta-function-like behaviour on $i$th face and is equal to zero outside this face. Consequently, the derivative $\partial\psi_{i}/\partial r_{\beta}$ takes the form

[TABLE]

where

[TABLE]

and

[TABLE]

(Notably, the function $\eta_{k}(\mathbf{r})$ is one and the same for both faces adjacent to $k$th edge.) The term containing $\omega_{i}$ arises in Eq. (119) from differentiating the delta-functional factor in Eq. (118), and the terms with $\eta_{k}$—from differentiating Heaviside $\theta$-functions in a manner similar to Eq. (108). With Eqs. (117) and (119) one can represent second derivatives of $\chi$ as

[TABLE]

where index $i$ runs over all faces, and index $k$—over edges adjacent to face $i$. Changing the order of summations in the last sum, and introducing the tensor $\lambda^{(k)}_{\alpha\beta}$ in accordance with Eq. (33), one can rewrite this representation in a more convenient form:

[TABLE]

where index $k$ runs over all edges of the polyhedron. Hence, the integral in the right-hand side of Eq. (116) is a linear combination of similar integrals containing $\omega_{i}$ or $\eta_{k}$ instead of $\partial^{2}\chi/\partial r_{\alpha}\partial r_{\beta}$.

Let us consider the integral with $\omega_{i}$. For convenience, we temporarily assume that $i$th face is perpendicular to the coordinate axis $z$. Then, the function $\omega_{i}(\mathbf{r})$ can be expressed as

[TABLE]

where $z_{i}$ is the value of $z$ at $i$th face, and the function $\tilde{\chi}_{i}(x,y)$ defines the polygon of $i$th face: $\tilde{\chi}_{i}(x,y)=1$ if the point $(x,y,z_{i})$ lies inside the face, otherwise $\tilde{\chi}_{i}(x,y)=0$. Therefore one can represent the integral with $\omega_{i}$ as follows:

[TABLE]

The inner integral can be taken by parts,

[TABLE]

and substitution of Eq. (126) into Eq. (125) leads to the following result:

[TABLE]

where the quantity $\Phi^{(i)}_{(mnp),z}$ is defined in Eq. (36). Though Eq. (127) was derived in a special coordinate frame, one can easily make it independent of the choice of the frame. For that, it is enough to rewrite the term $\Phi^{(i)}_{(mnp),z}$ as $n_{i\gamma}\Phi^{(i)}_{(mnp),\gamma}$, i. e., in the invariant manner:

[TABLE]

In order to calculate a similar integral with $\eta_{k}$, we choose (temporarily) a coordinate frame with axis $z$ along the vector $\mathbf{n}_{i}$, and axis $x$ along the vector $\mathbf{b}_{ik}$, so that axis $y$ is directed along $k$th edge. In this frame, the expression (121) can be rewritten as

[TABLE]

where $x_{k}$ and $z_{k}$ are coordinates of $k$th edge, and the function $\tilde{\chi}_{k}(y)$ is equal to 1 if the point $(x_{k},y,z_{k})$ lies on the edge, otherwise $\tilde{\chi}_{k}(y)=0$. Then it is obvious that

[TABLE]

according to the definition of $L^{(k)}_{mnp}$, Eq. (29). Of course, this result is invariant with respect to choice of the coordinate frame.

Finally, the combination of equations (123), (128), (130) provides the following result for the integral that appears in Eq. (116):

[TABLE]

Substituting this expression into Eq. (116), one obtains Eq. (34), which represents the quantity $\varphi^{(mnp)}_{,\alpha\beta}$ in terms of primitives.

Hence, Eq. (34) is proven, as well as Eq. (30).

Appendix B Reducing volume integrals to surface integrals

The aim of this Appendix consists in deriving Eq. (47). For simplicity, let us consider the case of

[TABLE]

that does not lead to any loss of generality. (In order to account for nonzero $\mathbf{R}$, one can make the substitution $\mathbf{r}\rightarrow\mathbf{r}-\mathbf{R}$.) Let us define a vector field $\mathbf{F}(\mathbf{r})$ as

[TABLE]

and apply Gauss’s theorem:

[TABLE]

The integrand in the left-hand side of Eq. (134) is

[TABLE]

Therefore the left-hand side of Eq. (134) is equal to $(m+n+p+2)\,\varphi_{mnp}$. The integrand in the right-hand side of Eq. (134) is

[TABLE]

At any point of $i$th face, the factor $\mathbf{r}\cdot\mathbf{n}_{i}$ has the same value

[TABLE]

where the quantity $h_{i}$ is defined in Eq. (44). Therefore this factor can be carried out of the integral, and the integral over $i$th face in Eq. (134) becomes equal to $h_{i}\,\Phi^{(i)}_{mnp}$. Hence, Eq. (134) takes the following form:

[TABLE]

that proves validity of Eq. (47).

Appendix C Converting primitives $\hat{\Phi}^{(i)}_{mn0}$ to $\hat{\Omega}^{(i)}_{mn0}$

Here we will derive Eq. (55). In order to simplify notations, we omit the hats at $\hat{x},\hat{y},\hat{z},\hat{r}$. Consider the two-dimensional vector field $(G_{x},G_{y})$ whose components are defined as

[TABLE]

and write down a two-dimensional version of Gauss’s theorem in the plane of $i$th face:

[TABLE]

where summation is over the edges surrounding face $i$; $b_{ikx}$ and $b_{iky}$ are components of the vector $\mathbf{b}_{ik}$ in the “tilted” frame.

The integrand of the surface integral in Eq. (140) can be rewritten as

[TABLE]

and, according to the definitions (49a) and (49b), the left-hand side of Eq. (140) obtains the following form:

[TABLE]

On the other hand, the integrand of the line integral in Eq. (140) is

[TABLE]

The factor $\mathbf{r}\cdot\mathbf{b}_{ik}$ is constant along the edge $k$, and has the value

[TABLE]

due to Eq. (45). Thus,

[TABLE]

Collecting together Eqs. (140), (142), and (145), one can get the relation

[TABLE]

which is equivalent to Eq. (55).

Appendix D Decreasing the index $n$ at $\hat{\Omega}^{(i)}_{mn0}$

In this Appendix we will derive the recursive relation (56) that allows a step-by-step decrease of the index $n$ at $\hat{\Omega}^{(i)}_{mn0}$. As in the previous Appendix, we omit the hats at $\hat{x},\hat{y},\hat{z},\hat{r}$.

Let us make simple transformations of the quantity $h_{i}\hat{\Phi}^{(i)}_{mn0}$:

[TABLE]

(Integrals are over the $i$th face.) On the other hand, one can use Eq. (55) for reducing the value of $h_{i}\hat{\Phi}^{(i)}_{mn0}$ to quantities $\hat{\Omega}^{(i)}_{mn0}$ and $\hat{L}^{(k)}_{mn0}$:

[TABLE]

where summation is over edges adjacent to face $i$. Since left-hand sides of Eqs. (147) and (148) are the same, then right-hand sides are equal to each other:

[TABLE]

Expressing $\hat{\Omega}^{(i)}_{m,n+2,0}$ from this equation, and taking into account that $\hat{\Omega}^{(i)}_{mn2}=h_{i}^{2}\hat{\Omega}^{(i)}_{mn0}$ due to Eq. (53), one can get the following result:

[TABLE]

Equation (56) can be obtained from Eq. (150) by means of the substitution $n\to n-2$.

Appendix E Reducing surface integrals $\hat{\Omega}^{(i)}_{m10}$ to line integrals

In this Appendix, we will derive Eq. (57) that simplifies the surface integral $\hat{\Omega}^{(i)}_{m10}$. As above, we will omit here the hats at $\hat{x},\hat{y},\hat{z},\hat{r}$.

Let us perform integration of $\hat{\Omega}^{(i)}_{m10}$ over $y$:

[TABLE]

where symbols $x_{1},x_{2},y_{1},y_{2}$ are explained in Fig. 10; $r(x,y)=(x^{2}+y^{2}+h_{i}^{2})^{1/2}$; $h_{i}$ is $z$-coordinate of $i$th face. The first integral in the right-hand side of Eq. (151) is taken over the upper part $y_{2}(x)$ of the periphery of $i$th face (between points A and B, see Fig. 10), and the second integral is over the lower part $y_{1}(x)$. It is convenient to express $dx$ via the line element $dl$ along the edge:

[TABLE]

By doing so, one can combine both integrals together, yielding

[TABLE]

i. e.

[TABLE]

Therefore we have obtained Eq. (57).

Appendix F Decreasing the index $m$ at $\hat{\Omega}^{(i)}_{m00}$

This Appendix is aimed to deduce the recursive rule (58) for reducing the index $m$ at $\hat{\Omega}^{(i)}_{m00}$. As above, we will omit here the hats at $\hat{x},\hat{y},\hat{z},\hat{r}$.

For this purpose, we rewrite the surface integral $\hat{\Omega}^{(i)}_{m20}$ in a manner similar to Eq. (151):

[TABLE]

and make the integration by parts, taking into account that $y\,dy=d(r^{2})/2$ and $d(y/r^{3})=(1/r^{3}-3y^{2}/r^{5})dy$:

[TABLE]

whence

[TABLE]

Substituting this into Eq. (155), one can obtain

[TABLE]

The first term in the right hand side is equal to $-h_{i}\hat{\Phi}^{(i)}_{m00}$ (by definition of $\hat{\Phi}^{(i)}_{m00}$); the second and the third terms can be combined together, literally repeating the steps that lead from Eq. (151) to Eq. (154). This results in the following representation of $\hat{\Omega}^{(i)}_{m20}$:

[TABLE]

where summation is over edges of $i$th face. The term $(-h_{i}\hat{\Phi}^{(i)}_{m00})$ can be reduced to $\Omega$’s by means of Eq. (147):

[TABLE]

Substituting this into Eq. (159), and expressing $\hat{\Omega}^{(i)}_{m02}$ as $h_{i}^{2}\hat{\Omega}^{(i)}_{m00}$ according to Eq. (53), one can obtain the following relation:

[TABLE]

from which one can get Eq. (58) by means of the substitution $m\rightarrow m-2$.

Appendix G Recursive formula for line integrals

Here we will obtain expressions (80), (81) for the integrals $\mathcal{L}^{(k)}_{1},\mathcal{L}^{(k)}_{2},\ldots$ This can be done using integration by parts, taking into account that $\xi\,d\xi=\frac{1}{2}\,d(\rho_{k}^{2}+\xi^{2})$:

[TABLE]

According to the definition (69), one can represent the last term as a combination of $\mathcal{L}^{(k)}_{t-2}$ and $\mathcal{L}^{(k)}_{t}$, yielding

[TABLE]

whence

[TABLE]

Square roots $\sqrt{\rho_{k}^{2}+\xi^{2}}$ at $\xi=\xi_{k1},\xi_{k2}$ have the meaning of distances from the point $\mathbf{R}$ to the ends of $k$th edge:

[TABLE]

Therefore, Eq. (164) takes the following form:

[TABLE]

Substituting $\xi_{k1}$ and $\xi_{k2}$ from Eq. (71) into Eq. (167), and resolving the latter equation with respect to $\mathcal{L}^{(k)}_{t}$, one can get the recursive relation (80).

At $t=1$, Eq. (167) is simplified to

[TABLE]

Appendix H Comparison with a closed-form solution found in the literature

In order to check the validity of our analytical results, it is natural to compare them with similar results obtained by other authors. But it seems that there are no analytical calculations in the literature, related to polyhedral-shaped inclusions with a smoothly varied lattice misfit. For this reason, we choose as a “benchmark” the analytical results by GlasGlas (2001) for the strain distribution due to an inclusion having the shape of a truncated pyramid with a constant misfit strain. Such an analytical calculations can be found in Fig. 3b of Ref. Glas, 2001.

In Fig. 11, we reproduce these calculations by our method—namely, by means of Eqs. (82) and (98). For the details of the inclusion geometry, see the caption to the figure. One can see that isolines in Fig. 11 are identical to ones in Fig. 3b of Ref. Glas, 2001. That is, our method provides the same strain distribution as the closed-form analytical solution found by GlasGlas (2001) for the particular case of a truncated pyramid with a constant misfit strain.

Appendix I Comparison with numerical calculations

The other useful test is comparison with numerical calculations. One difficulty in the way of such a comparison is that numerical methods deal with an elastic body of final size, whereas our analytical approach is valid for an inclusion in an infinitely large matrix. Therefore, in order to ensure that the analytical method and is compatible with numerical calculations, one should perform numerical tests for different sizes of the matrix and check whether the strain distribution converges to the analytical solution when the matrix size grows.

An example of such a comparison is shown in Fig. 12. As a test sample, a pyramid-shaped inclusion with a linearly varying misfit strain is chosen. The parameters of the inclusion (see details in the caption to the figure) are taken the same as ones in Fig. 8. Therefore, the analytical solution in Fig. 12c is just the same as shown in Fig. 8. In order to obtain the numerical solutions, we used the COMSOL Multiphysics software. Simulations were performed by the finite-difference method for two different sizes of the “box”, which contains the pyramidal inclusion: $20\times 20\times 15$ nm (Fig. 12a) and $30\times 30\times 25$ nm (Fig. 12b), whereas the inclusion size is $10\times 10\times 5$ nm. One can see that, indeed, with increasing the size of the “box”, the strain distribution converges to the analytical solution.

It is worth to note that numerical results shown in Fig. 12a,b demanded about one hour of machine time on PC and about 20 gigabytes of memory, whereas the analytical calculation (Fig. 12c) lasts less than one minute and demands less then one megabyte of memory.

Appendix J Pseudocode representation of the calculation algorithm

Calculation of the strain tensor $\varepsilon_{\alpha\beta}(\mathbf{r})$ for a set of points $\mathbf{r}$ is a two-step procedure. The first step is the calculation of polynomials $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$, $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$ for all faces, edges and vertices of the polyhedral inclusion, and for all values of tensor indices $\alpha\beta$. The second step is applying Eq. (82) for each point $\mathbf{r}$, that gives the value of $\varepsilon_{\alpha\beta}(\mathbf{r})$.

Below we represent these steps in a somewhat informal style, as Algorithm 1 and Algorithm 2, respectively. Algorithm 1 should be called six times—namely, for each pair $xx,yy,zz,xy,xz$ and $yz$ of the indices $\alpha\beta$. Then, Algorithm 2 is called as many times, as many points $\mathbf{r}$ there are, at which the strain tensor $\varepsilon_{\alpha\beta}$ is to be calculated.

Algorithm 2 does not demand any further comments, since it simply calculates $\varepsilon_{\alpha\beta}$ via Eq. (82). But Algorithm 1 surely demands explanations. The idea of Algorithm 1 consists in the fact that, at each stage of the transformations described in Sections II–VI, the function $\varepsilon_{\alpha\beta}(\mathbf{R})$ is represented as a sum of contributions of such functions as $\varphi_{mnp}(\mathbf{R})$, $\Omega^{(i)}_{mnp}(\mathbf{R})$, etc., each of them being multiplied by some polynomial. Let us denote the polynomial corresponding to the function $\varphi_{mnp}(\mathbf{R})$ as $[\varphi_{mnp}]$, etc. Therefore,

[TABLE]

Then, let there be some identity connecting a function $A(\mathbf{R})$ with other functions $B(\mathbf{R}),\,C(\mathbf{R})$:

[TABLE]

where $b(\mathbf{R}),\,c(\mathbf{R})$ are some polynomials. Then, one can rewrite the contribution $[A]\;A$ of the function $A$ into $\varepsilon_{\alpha\beta}$ as

[TABLE]

i.  e. transform it into contributions of $B$ and $C$. Such a transformation can be implemented in the algorithm as a transformation of polynomial $[A]$ into polynomials $[B],\;[C]$ as follows:

1:$[B]=[A]\;b$

2:$[C]=[A]\;c$

3:$[A]=0$

Implementing different identities from Sections II–VI in a due order, one can transform all the polynomials like $[\varphi_{mnp}]$, $[\Omega^{(i)}_{mnp}]$, etc. into the following ones: $[\Omega^{(i)}_{000}]\equiv\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$, $[L^{(k)}_{000}]\equiv\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$, and $[|\mathbf{R}-\mathbf{r}_{s}|]\equiv\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$ for different $i$, $k$, $s$. This is how one can find out the polynomials $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$. This is the idea, on which Algorithm 1 is based.

It is important to note that Algorithm 1 is by no means optimized. Its performance and memory requirements can be easily improved, but we leave it here in its non-optimized version, in order to keep its form as close to equations of Sections II–VI as possible. The numbers of used equations are indicated in comments.

Algorithm 1: finding out the polynomials $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$, $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$

Input:

tensor indices $\alpha,\beta$ determining which strain component $\varepsilon_{\alpha\beta}$ we are going to calculate;

polynomial $C$, which defines the coordinate dependence of the eigenstrain within the inclusion;

degree $N$ of the polynomial $C$;

the geometry of the polyhedral inclusion: coordinates of vertices; the set of unit vectors $\mathbf{n}_{i},\mathbf{b}_{ik},\mathbf{l}_{k}$ that describes orientations of faces and edges; for each face $i$ (edge $k$), one arbitrarily chosen point $\mathbf{r}_{i}$ ($\mathbf{r}_{k}$) on it;

Poisson ratio $\nu$.

Output: polynomials $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$ for each face $i$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$ for each edge $k$, $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$ for each vertex $s$.

Variables:

integer numbers: $i,i_{1},i_{2},k,m,n,p,m^{\prime},n^{\prime},p^{\prime},t,\gamma$;

real numbers, vectors and arrays: $\Lambda,\lambda^{(k)}_{\alpha\beta},\hat{\mathbf{e}}_{x},\hat{\mathbf{e}}_{y},\hat{\mathbf{e}}_{z},\mathbf{r}_{\mathrm{A}},\mathbf{r}_{\mathrm{B}},T^{m^{\prime}n^{\prime}p^{\prime}}_{mnp}$;

polynomials (of degree $N$) and arrays of polynomials: $[\varphi_{m,n,p}]$, $[\Phi^{(i)}_{m,n,p}]$, $[\Omega^{(i)}_{m,n,p}]$, $[L^{(k)}_{m,n,p}]$, $[\hat{\Phi}^{(i)}_{m,n,p}]$, $[\hat{\Omega}^{(i)}_{m,n,p}]$, $[\hat{L}^{(k)}_{m,n,p}]$, $[\varphi^{(mnp)}_{\alpha\beta}]$, $[\varphi_{(mnp),\alpha\beta}]$, $[\Phi^{(i)}_{(mnp),\gamma}]$, $[\mathcal{L}^{(k)}_{t}]$, $h_{i}$, $B_{ik}$, $\mathbf{r}_{0}^{(k)}$, $\rho^{2}_{k}$, $\tilde{C}_{mnp}$, $P_{\mathrm{A}}$, $P_{\mathrm{B}}$.

(Note 1: polynomials can be stored as real arrays $(N+1)\times(N+1)\times(N+1)$ of their coefficients.)

(Note 2: this algorithm can be modified for calculating the potential or its first derivatives, but in this case some polynomials will have degrees $N+2$ and $N+1$, respectively.)

1:$\triangleright$ Initialization:

2:set polynomials $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$, $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$ to zero

3:for all $m\geq 0,\;n\geq 0,\;p\geq 0\;\;(m+n+p\leq N)$ do

4:     $[\varphi_{m,n,p}]=0$

5:     for all faces $i$ do

6:         $[\Phi^{(i)}_{m,n,p}]=0$

7:     end for

8:end for

9:for $t=0$ to $N$ do

10:     for all edges $k$ do

11:         $[\mathcal{L}^{(k)}_{t}]=0$

12:     end for

13:end for

14:for all faces $i$ do

15:     define polynomial $h_{i}$ as follows: $h_{i}(\mathbf{R})=(\mathbf{r}_{i}-\mathbf{R})\cdot\mathbf{n}_{i}$ $\triangleright$ Eq. (44)

16:     for all edges $k$ adjacent to face $i$ do

17:         define polynomial $B_{ik}$ as follows: $B_{ik}(\mathbf{R})=(\mathbf{r}_{k}-\mathbf{R})\cdot\mathbf{b}_{ik}$ $\triangleright$ Eq. (45)

18:     end for

19:end for

20:for all edges $k$ do

21:     define polynomials \mathbf{r}_{0}^{(k)}=\big{\{}x_{0}^{(k)},y_{0}^{(k)},z_{0}^{(k)}\big{\}} and $\rho^{2}_{k}$ as follows:

22:     \mathbf{r}_{0}^{(k)}(\mathbf{R})=\mathbf{r}_{k}-\mathbf{l}_{k}\big{[}\mathbf{l}_{k}\cdot(\mathbf{r}_{k}-\mathbf{R})\big{]}, $\triangleright$ Eq. (73)

23:     \rho^{2}_{k}(\mathbf{R})=|\mathbf{r}_{k}-\mathbf{R}|^{2}-\big{[}\mathbf{l}_{k}\cdot(\mathbf{r}_{k}-\mathbf{R})\big{]}^{2}. $\triangleright$ Eq. (70)

24:end for

25:$\Lambda=\frac{1}{4\pi}\frac{1+\nu}{1-\nu}$ $\triangleright$ Eq. (10)

26:

27:$\triangleright$ Steps 1, 2:

28:for all $m\geq 0,\;n\geq 0,\;p\geq 0\;\;(m+n+p\leq N)$ do

29:     calculate polynomial $\tilde{C}_{mnp}$ according to Eq. (19)

30:     $[\varphi^{(mnp)}_{\alpha\beta}]=-\Lambda\tilde{C}_{mnp}$ $\triangleright$ Eqs. (12), (22)

31:end for

32:for all faces $i$ do

33:     add $\frac{\delta_{\alpha\beta}}{4\pi}C$ to $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$ $\triangleright$ Eqs. (12), (42), (3)

34:end for

35:$\triangleright$ Step 3:

36:for all $m\geq 0,\;n\geq 0,\;p\geq 0\;\;(m+n+p\leq N)$ do

37:     $[\varphi_{(mnp),\alpha\beta}]=[\varphi^{(mnp)}_{\alpha\beta}]$ $\triangleright$ begin of Eq. (34)

38:     for all faces $i$ do

39:         for all $\gamma\in\{x,y,z\}$ do

40:              $[\Phi^{(i)}_{(mnp),\gamma}]=(-n_{i\alpha}\delta_{\beta\gamma}-n_{i\beta}\delta_{\alpha\gamma}+n_{i\alpha}n_{i\beta}n_{i\gamma})\;[\varphi^{(mnp)}_{\alpha\beta}]$

41:         end for

42:         $[\Omega^{(i)}_{m,n,p}]=n_{i\alpha}n_{i\beta}\;[\varphi^{(mnp)}_{\alpha\beta}]$

43:     end for

44:     for all edges $k$ do

45:         find two faces $i_{1}$, $i_{2}$ adjacent to edge $k$

46:         $\lambda^{(k)}_{\alpha\beta}=n_{i_{1}\alpha}b_{i_{1}k\beta}+n_{i_{2}\alpha}b_{i_{2}k\beta}$ $\triangleright$ Eq. (33)

47:         $[L^{(k)}_{m,n,p}]=\lambda^{(k)}_{\alpha\beta}\;[\varphi^{(mnp)}_{\alpha\beta}]$

48:     end for

49:     $[\varphi^{(mnp)}_{\alpha\beta}]=0$ $\triangleright$ end of Eq. (34)

50:     if $(\alpha\beta)=(xx)$ and $m\geq 2$ then $\triangleright$ begin of Eq. (37)

51:         add $m(m-1)\;[\varphi_{(mnp),\alpha\beta}]$ to $[\varphi_{m-2,n,p}]$

52:     else if $(\alpha\beta)=(yy)$ and $n\geq 2$ then

53:         add $n(n-1)\;[\varphi_{(mnp),\alpha\beta}]$ to $[\varphi_{m,n-2,p}]$

54:     else if $(\alpha\beta)=(zz)$ and $p\geq 2$ then

55:         add $p(p-1)\;[\varphi_{(mnp),\alpha\beta}]$ to $[\varphi_{m,n,p-2}]$

56:     else if $(\alpha\beta)=(xy)$ or $(yx)$ and $m\geq 1$ and $n\geq 1$ then

57:         add $mn\;[\varphi_{(mnp),\alpha\beta}]$ to $[\varphi_{m-1,n-1,p}]$

58:     else if $(\alpha\beta)=(xz)$ or $(zx)$ and $m\geq 1$ and $p\geq 1$ then

59:         add $mp\;[\varphi_{(mnp),\alpha\beta}]$ to $[\varphi_{m-1,n,p-1}]$

60:     else if $(\alpha\beta)=(yz)$ or $(zy)$ and $n\geq 1$ and $p\geq 1$ then

61:         add $np\;[\varphi_{(mnp),\alpha\beta}]$ to $[\varphi_{m,n-1,p-1}]$

62:     end if

63:     $[\varphi_{(mnp),\alpha\beta}]=0$ $\triangleright$ end of Eq. (37)

64:     for all faces $i$ do

65:         if $m\geq 1$ then $\triangleright$ begin of Eq. (38)

66:              add $m\;[\Phi^{(i)}_{(mnp),x}]$ to $[\Phi^{(i)}_{m-1,n,p}]$

67:         end if

68:         if $n\geq 1$ then

69:              add $n\;[\Phi^{(i)}_{(mnp),y}]$ to $[\Phi^{(i)}_{m,n-1,p}]$

70:         end if

71:         if $p\geq 1$ then

72:              add $p\;[\Phi^{(i)}_{(mnp),z}]$ to $[\Phi^{(i)}_{m,n,p-1}]$

73:         end if

74:         $[\Phi^{(i)}_{(mnp),x}]=0$; $\;\;[\Phi^{(i)}_{(mnp),y}]=0$; $\;\;[\Phi^{(i)}_{(mnp),z}]=0$ $\triangleright$ end of Eq. (38)

75:     end for

76:end for

77:$\triangleright$ Step 4:

78:for all $m\geq 0,\;n\geq 0,\;p\geq 0\;\;(m+n+p\leq N)$ do

79:     for all faces $i$ do

80:         add $\frac{1}{m+n+p+2}h_{i}\;[\varphi_{m,n,p}]$ to $[\Phi^{(i)}_{m,n,p}]$ $\triangleright$ begin of Eq. (47)

81:     end for

82:     $[\varphi_{m,n,p}]=0$ $\triangleright$ end of Eq. (47)

83:end for

84:

85:for all faces $i$ do $\triangleright$ Now we are going to the tilted frame

86:     

87:     for all $m\geq 0,\;n\geq 0,\;p\geq 0\;\;(m+n+p\leq N)$ do

88:         $[\hat{\Phi}^{(i)}_{m,n,p}]=0$ $\triangleright$ Initialize polynomials related to the tilted frame

89:         $[\hat{\Omega}^{(i)}_{m,n,p}]=0$

90:         for all edges $k$ adjacent to face $i$ do

91:              $[\hat{L}^{(k)}_{m,n,p}]=0$

92:         end for

93:     end for

94:     $\triangleright$ Step 5:

95:     choose an orthonormal set of vectors $\hat{\mathbf{e}}_{x},\hat{\mathbf{e}}_{y},\hat{\mathbf{e}}_{z}$ so that $\hat{\mathbf{e}}_{z}=\mathbf{n}_{i}$

96:     for all $m\geq 0,\;n\geq 0,\;p\geq 0\;\;(m+n+p\leq N)$ do

97:         for all $m^{\prime}\geq 0,\;n^{\prime}\geq 0,\;p^{\prime}\geq 0$ such that $m^{\prime}+n^{\prime}+p^{\prime}=m+n+p$ do

98:              calculate coefficient $T^{m^{\prime}n^{\prime}p^{\prime}}_{mnp}$ from Eq. (51)

99:              add $T^{m^{\prime}n^{\prime}p^{\prime}}_{mnp}\;[\Phi^{(i)}_{m,n,p}]$ to $[\hat{\Phi}^{(i)}_{m^{\prime},n^{\prime},p^{\prime}}]$ $\triangleright$ Eq. (50a)

100:              add $T^{m^{\prime}n^{\prime}p^{\prime}}_{mnp}\;[\Omega^{(i)}_{m,n,p}]$ to $[\hat{\Omega}^{(i)}_{m^{\prime},n^{\prime},p^{\prime}}]$ $\triangleright$ Eq. (50b)

101:         end for

102:         $[\Phi^{(i)}_{m,n,p}]=0$; $\quad[\Omega^{(i)}_{m,n,p}]=0$ $\triangleright$ Eq. (50)

103:     end for

104:     $\triangleright$ Step 6:

105:     for all $m\geq 0,\;n\geq 0,\;p\geq 1\;\;(m+n+p\leq N)$ do

106:         add $(h_{i})^{p}\;[\hat{\Phi}^{(i)}_{m,n,p}]$ to $[\hat{\Phi}^{(i)}_{m,n,0}]$; $\quad[\hat{\Phi}^{(i)}_{m,n,p}]=0$ $\triangleright$ Eq. (53)

107:         add $(h_{i})^{p}\;[\hat{\Omega}^{(i)}_{m,n,p}]$ to $[\hat{\Omega}^{(i)}_{m,n,0}]$; $\quad[\hat{\Omega}^{(i)}_{m,n,p}]=0$ $\triangleright$ Eq. (53)

108:     end for

109:     $\triangleright$ Step 7:

110:     for all $m\geq 0,\;n\geq 0\;\;(m+n\leq N)$ do

111:         add $\frac{h_{i}}{m+n+1}\;[\hat{\Phi}^{(i)}_{m,n,0}]$ to $[\hat{\Omega}^{(i)}_{m,n,0}]$ $\triangleright$ begin of Eq. (55)

112:         for all edges $k$ adjacent to face $i$ do

113:              add $\frac{B_{ik}}{m+n+1}\;[\hat{\Phi}^{(i)}_{m,n,0}]$ to $[\hat{L}^{(k)}_{m,n,0}]$

114:         end for

115:         $[\hat{\Phi}^{(i)}_{m,n,0}]=0$ $\triangleright$ end of Eq. (55)

116:     end for

117:     $\triangleright$ Step 8:

118:     if $N\geq 2$ then

119:         for $m=0$ to $N-2$ do

120:              for $n=N-m$ down to 2 do

121:                  add $-[\hat{\Omega}^{(i)}_{m,n,0}]$ to $[\hat{\Omega}^{(i)}_{m+2,n-2,0}]$ $\triangleright$ begin of Eq. (56)

122:                  add $-\frac{m+n}{m+n-1}(h_{i})^{2}\;[\hat{\Omega}^{(i)}_{m,n,0}]$ to $[\hat{\Omega}^{(i)}_{m,n-2,0}]$

123:                  for all edges $k$ adjacent to face $i$ do

124:                       add $-\frac{h_{i}B_{ik}}{m+n-1}\;[\hat{\Omega}^{(i)}_{m,n,0}]$ to $[\hat{L}^{(k)}_{m,n-2,0}]$

125:                  end for

126:                  $[\hat{\Omega}^{(i)}_{m,n,0}]=0$ $\triangleright$ end of Eq. (56)

127:              end for

128:         end for

129:     end if

130:     if $N\geq 1$ then

131:         for $m=0$ to $N-1$ do

132:              for all edges $k$ adjacent to face $i$ do

133:                  add $h_{i}(\mathbf{b}_{ik}\cdot\hat{\mathbf{e}}_{y})\;[\hat{\Omega}^{(i)}_{m,1,0}]$ to $[\hat{L}^{(k)}_{m,0,0}]$ $\triangleright$ begin of Eq. (57)

134:              end for

135:              $[\hat{\Omega}^{(i)}_{m,1,0}]=0$ $\triangleright$ end of Eq. (57)

136:         end for

137:     end if

138:     $\triangleright$ Step 9:

139:     if $N\geq 2$ then

140:         for $n=N$ down to 2 do

141:              add $-(h_{i})^{2}\;[\hat{\Omega}^{(i)}_{m,0,0}]$ to $[\hat{\Omega}^{(i)}_{m-2,0,0}]$ $\triangleright$ begin of Eq. (58)

142:              for all edges $k$ adjacent to face $i$ do

143:                  add $-h_{i}(\mathbf{b}_{ik}\cdot\hat{\mathbf{e}}_{y})\;[\hat{\Omega}^{(i)}_{m,0,0}]$ to $[\hat{L}^{(k)}_{m-2,1,0}]$

144:              end for

145:              $[\hat{\Omega}^{(i)}_{m,0,0}]=0$ $\triangleright$ end of Eq. (58)

146:         end for

147:     end if

148:     if $N\geq 1$ then

149:         for all edges $k$ adjacent to face $i$ do

150:              add $h_{i}(\mathbf{b}_{ik}\cdot\hat{\mathbf{e}}_{x})\;[\hat{\Omega}^{(i)}_{1,0,0}]$ to $[\hat{L}^{(k)}_{0,0,0}]$ $\triangleright$ begin of Eq. (59)

151:         end for

152:         $[\hat{\Omega}^{(i)}_{1,0,0}]=0$ $\triangleright$ end of Eq. (59)

153:     end if

154:     add $[\hat{\Omega}^{(i)}_{0,0,0}]$ to $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$; $\quad[\hat{\Omega}^{(i)}_{0,0,0}]=0$ $\triangleright$ (since $\hat{\Omega}^{(i)}_{000}\equiv\Omega^{(i)}$)

155:     

156:     $\triangleright$ Step 10 in the tilted frame:

157:     for all edges $k$ adjacent to face $i$ do

158:         for all $m\geq 0,\;n\geq 0\;\;(m+n\leq N)$ do

159:              $p=0$

160:              for $t=0$ to $m+n+p$ do

161:                  calculate polynomial $\hat{c}_{mnp,t}$ according to Eq. (78)

162:                  add $\hat{c}_{mnp,t}\;[\hat{L}^{(k)}_{m,n,p}]$ to $[\mathcal{L}^{(k)}_{t}]$ $\triangleright$ begin of Eq. (77)

163:              end for

164:              $[\hat{L}^{(k)}_{m,n,p}]=0$ $\triangleright$ end of Eq. (77)

165:         end for

166:     end for

167:     

168:end **for$\triangleright$ Return back from the tilted frame **

169:

170:$\triangleright$ Step 10 in the non-tilted frame:

171:for all edges $k$ do

172:     for all $m\geq 0,\;n\geq 0,\;p\geq 0\;\;(m+n+p\leq N)$ do

173:         for $t=0$ to $m+n+p$ do

174:              calculate polynomial $c_{mnp,t}$ according to Eq. (75)

175:              add $c_{mnp,t}\;[L^{(k)}_{m,n,p}]$ to $[\mathcal{L}^{(k)}_{t}]$ $\triangleright$ begin of Eq. (76)

176:         end for

177:         $[L^{(k)}_{m,n,p}]=0$ $\triangleright$ end of Eq. (76)

178:     end for

179:end for

180:$\triangleright$ Step 11:

181:for all edges $k$ do

182:     Let A and B be numbers of the two vertices connected by edge $k$, and unit vector $\mathbf{l}_{k}$ be directed from point A to point B. Let $\mathbf{r}_{\mathrm{A}}$ and $\mathbf{r}_{\mathrm{B}}$ be position vectors of these points.

183:     Define two polynomials $P_{\mathrm{A}}$ and $P_{\mathrm{B}}$ as follows:

184:     $P_{\mathrm{A}}(\mathbf{R})=\mathbf{l}_{k}\cdot(\mathbf{r}_{\mathrm{A}}-\mathbf{R})$,  $P_{\mathrm{B}}(\mathbf{R})=\mathbf{l}_{k}\cdot(\mathbf{r}_{\mathrm{B}}-\mathbf{R})$.

185:     if $N\geq 2$ then

186:         for $t=N$ down to 2 do

187:              add $-\frac{t-1}{t}\rho^{2}_{k}\;[\mathcal{L}^{(k)}_{t}]$ to $[\mathcal{L}^{(k)}_{t-2}]$ $\triangleright$ begin of Eq. (80)

188:              add $\frac{(P_{\mathrm{B}})^{t-1}}{t}\;[\mathcal{L}^{(k)}_{t}]$ to $\mathcal{C}^{(\mathrm{B})}_{\varepsilon\alpha\beta}$

189:              add $-\frac{(P_{\mathrm{A}})^{t-1}}{t}\;[\mathcal{L}^{(k)}_{t}]$ to $\mathcal{C}^{(\mathrm{A})}_{\varepsilon\alpha\beta}$

190:              $[\mathcal{L}^{(k)}_{t}]=0$ $\triangleright$ end of Eq. (80)

191:         end for

192:     end if

193:     if $N\geq 1$ then

194:         add $[\mathcal{L}^{(k)}_{1}]$ to $\mathcal{C}^{(\mathrm{B})}_{\varepsilon\alpha\beta}$ $\triangleright$ begin of Eq. (81)

195:         add $-[\mathcal{L}^{(k)}_{1}]$ to $\mathcal{C}^{(\mathrm{A})}_{\varepsilon\alpha\beta}$

196:         $[\mathcal{L}^{(k)}_{1}]=0$ $\triangleright$ end of Eq. (81)

197:     end if

198:     add $[\mathcal{L}^{(k)}_{0}]$ to $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$; $\quad[\mathcal{L}^{(k)}_{0}]=0$ $\triangleright$ (since $\mathcal{L}^{(k)}_{0}\equiv L^{(k)}$)

199:end for

Algorithm 2: calculating the strain tensor $\varepsilon_{\alpha\beta}(\mathbf{r})$ at one point $\mathbf{r}$

Input:

coordinates of the point $\mathbf{r}$ where the strain tensor is to be calculated;

coordinates of vertices of the polyhedral inclusion, connections between vertices, edges and faces of the polyhedron;

the set of polynomials $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}$, $\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}$, $\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}$ for all faces, edges and vertices of the polyhedron (as calculated by Algorithm 1).

Output: the value of the strain tensor $\varepsilon_{\alpha\beta}$ at point $\mathbf{r}$.

Variables:

integer numbers: $i,k,s,\alpha,\beta$;

real numbers: $A,B,C,\Omega,L,R,l_{A},l_{B},l_{AB}$.

1:for $\alpha\in\{x,y,z\}$ do

2:     for $\beta\in\{x,y,z\}$ do

3:         $\varepsilon_{\alpha\beta}=0$

4:     end for

5:end for

6:

7:for all faces $i$ do $\triangleright$ the first sum in rhs of Eq. (82)

8:     Calculate the signed solid angle $\Omega$ subtended by face $i$ at point $\mathbf{r}$.

9:     (The sign of $\Omega$ is positive (negative) when the outer (the inner) side of face $i$ is seen from point $\mathbf{r}$.)

10:     $\triangleright$ Recipes for calculating solid angles can be found in Refs. Nenashev and Dvurechenskii, 2007; Werner and Scheeres, 1996.

11:     for all pairs $\alpha\beta\in\{xx,yy,zz,xy,xz,yz\}$ do

12:         $A=\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}(\mathbf{r})$ $\triangleright$ here we calculate the value of the polynomial function $\mathcal{A}^{(i)}_{\varepsilon\alpha\beta}(x,y,z)$ at point $\mathbf{r}$

13:         $\varepsilon_{\alpha\beta}=\varepsilon_{\alpha\beta}+A\Omega$

14:     end for

15:end for

16:

17:for all edges $k$ do $\triangleright$ the second sum in rhs of Eq. (82)

18:     Let $\mathbf{r}_{A}$ and $\mathbf{r}_{B}$ be radius vectors of the two ends of edge $k$.

19:     $l_{A}=|\mathbf{r}-\mathbf{r}_{A}|$; $\;\;l_{B}=|\mathbf{r}-\mathbf{r}_{B}|$; $\;\;l_{AB}=|\mathbf{r}_{A}-\mathbf{r}_{B}|$

20:     $L=\log\frac{l_{A}+l_{B}+l_{AB}}{l_{A}+l_{B}-l_{AB}}$ $\triangleright$ Eq. (43)

21:     for all pairs $\alpha\beta\in\{xx,yy,zz,xy,xz,yz\}$ do

22:         $B=\mathcal{B}^{(k)}_{\varepsilon\alpha\beta}(\mathbf{r})$

23:         $\varepsilon_{\alpha\beta}=\varepsilon_{\alpha\beta}+BL$

24:     end for

25:end for

26:

27:for all vertices $s$ do $\triangleright$ the third sum in rhs of Eq. (82)

28:     Let $\mathbf{r}_{s}$ be the radius vector of vertex $s$.

29:     $R=|\mathbf{r}-\mathbf{r}_{s}|$;

30:     for all pairs $\alpha\beta\in\{xx,yy,zz,xy,xz,yz\}$ do

31:         $C=\mathcal{C}^{(s)}_{\varepsilon\alpha\beta}(\mathbf{r})$

32:         $\varepsilon_{\alpha\beta}=\varepsilon_{\alpha\beta}+CR$

33:     end for

34:end for

35:

36:$\varepsilon_{yx}=\varepsilon_{xy}$; $\quad\varepsilon_{zx}=\varepsilon_{xz}$; $\quad\varepsilon_{zy}=\varepsilon_{yz}$