The Green tensor of Mindlin's anisotropic first strain gradient elasticity

We derive the Green tensor of Mindlin's anisotropic first strain gradient elasticity. The Green tensor is valid for arbitrary anisotropic materials, with up to 21 elastic constants and 171 gradient elastic constants in the general case of triclinic media. In contrast to its classical counterpart, the Green tensor is non-singular at the origin, and it converges to the classical tensor a few characteristic lengths away from the origin. Therefore, the Green tensor of Mindlin's first strain gradient elasticity can be regarded as a physical regularization of the classical anisotropic Green tensor. The isotropic Green tensor and other special cases are recovered as particular instances of the general anisotropic result. The Green tensor is implemented numerically and applied to the Kelvin problem with elastic constants determined from interatomic potentials. Results are compared to molecular statics calculations carried out with the same potentials.


Introduction
Green functions are objects of fundamental importance in field theories, since they represent the fundamental solution of linear inhomogeneous partial differential equations (PDEs) from which any particular solution can be obtained via convolution with the source term [Green (1828)]. Moreover, Green functions are the basis of important numerical methods for boundary value problems, such as the boundary element method [Becker (1992)], and they provide "flexible" boundary conditions for atomistic simulations [Trinkle (2008) ]. In the context of linear elasticity, the Green function is a tensor-valued function of rank two, also known as the Green tensor. When contracted with a concentrated force acting at the origin, the Green tensor yields the displacement field in an infinite elastic medium. [Lord Kelvin (1882)] first derived the closed-form expression of the classical Green tensor for isotropic materials. For anisotropic materials, [Lifshitz and Rosenzweig (1947)] and [Synge (1957)] were able to derive the Green tensor in terms of an integral expression over the equatorial circle of the unit sphere in Fourier space. [Barnett (1972)] extended this result to the first two derivatives of the Green tensor, and showed that the line-integral representation is well suited for numerical integration (see also [Bacon et al. (1979), Teodosiu (1982)]).
The Green tensor and its derivatives are singular at the origin, ultimately as a consequence of the lack of intrinsic length scales in the classical theory of elasticity. The unphysical singularities in the elastic fields derived from the Green tensor hinder their applicability in nano-mechanics, including the elastic theory of defects such as cracks, dislocations and inclusions [Mura (1987), Askes and Aifantis (2011)]. Generalized elastic field theories with intrinsic length scales have been proposed in the context of micro-continuum theories [Eringen (1999)], non-local theories [Eringen (2002)], and gradient theories [Kröner (1963), Mindlin (1964), , Mindlin (1972), Mindlin and Eshel (1968)]. In particular, Mindlin's anisotropic strain gradient elasticity has received renewed attention as a tool to solve engineering problems at the micro-and nano-scales for realistic materials [Polizzotto (2018) ]. Only recently, the structure of the gradient-elastic tensor has been rationalized for different material symmetry classes [Auffray et al. (2013)], and its atomistic representation and ensuing determination from interatomic potentials has become available [Admal et al (2016)].
The number of independent strain gradient elastic moduli ranges from 5 for isotropic materials, to 171 in the general case of triclinic materials. While simple expressions of the Green tensor exist for the isotropic case [Rogula (1973), ], and for simplified anisotropic theories , b], the Green tensor of the fully anisotropic theory of Mindlin's strain gradient elasticity has remained so far a rather elusive object. [Rogula (1973)] provided an expression for the Green tensor in gradient elasticity of arbitrary order, which involves a sum of terms associated with the roots of a certain characteristic polynomial. However, such representation renders its numerical implementation rather impractical, and it conceals the mathematical structure of the Green tensor in relationship to its classical counterpart.
The objective of this paper is to derive a simple representation of the Green tensor of Mindlin's anisotropic first strain gradient elasticity, whose integral kernel involves only matrix operations suitable for efficient numerical implementation. Following a brief summary of Mindlin's anisotropic first strain gradient elasticity in section 2, we derive the matrix representation of the Green tensor in section 3. It is shown that the Green tensor is non-singular at the origin, while its first gradient is finite but discontinuous at the origin. The classical tail of the Green tensor, as well as its classical limit for vanishing gradient parameters are easily recovered from the non-singular expression. In section 4 we demonstrate that the Green tensor generalizes other expressions found in the literature. In section 5 we consider the Kelvin problem and compare the prediction of the Green tensor to atomistic calculations.

Mindlin's anisotropic gradient elasticity
Let us consider an infinite elastic body in three-dimensional space and assume that the gradient of the displacement field u is additively decomposed into an elastic distortion tensor β and an inelastic 1 eigen-distortion tensor β * : In the linearized theory of Mindlin's form-II first strain gradient elasticity [Mindlin (1964), , Mindlin and Eshel (1968), Mindlin (1972)], the strain energy density of an homogeneous and centrosymmetric 2 material is given by The strain energy density (2) is a function of the infinitesimal elastic strain tensor and of its gradient e ij,m . The tensor C is the standard rank-four tensor of elastic constants. By virtue of the symmetries it possesses up to 21 independent constants with units of eV/Å 3 . The tensor D is the rank-six tensor of strain gradient elastic constants, with symmetries It has units of eV/Å. In the general case of triclinic materials the number of independent constants in the tensor D is equal to 171 [Auffray et al. (2013)].
The quantities conjugate to the elastic strain tensor and its gradient are the Cauchy stress tensor σ and the double stress tensor τ , respectively. These are defined as: In the presence of a body forces density b, the static Lagrangian density of the system becomes: where is the potential of the body force. The condition of static equilibrium is expressed by the Euler-Lagrange equation In terms of the Cauchy and double stress tensors, Eq. (10) takes the following form [Mindlin (1964)]: Using Eqs.
(1) (6) (7), Eq. (11) can be cast in the following equation for displacements: In Eq. (12), L M ik denotes the differential operator of Mindlin's anisotropic first strain gradient elasticity while is the forcing term. Note that the second term on the right hand side of Eq. (14) is an "effective" internal force due to the inelastic eigen-distortion, and arises in the presence of material defects, such as inclusions, cracks, and dislocations. This term is the gradient version of the internal force in Mura's eigen-strain theory [Mura (1987)].
The Green tensor of Mindlin's anisotropic first strain gradient elasticity 5 3 The Green tensor of Mindlin's first strain gradient elasticity In this section, we derive the three-dimensional Green tensor of the operator (13). To this end, we seek the solution to Eq. (12) in the form where the symbol * indicates convolution over the three-dimensional space, and G is the Green tensor of Mindlin's anisotropic differential operator L M . Substituting Eq. (15) into Eq. (12), one finds that G satisfies the following inhomogeneous PDE: In Eq. (16), δ ij is the Kronecker symbol, while δ is the three-dimensional Dirac δ-distribution.
Taking the Fourier transform 3 of Eq. (16), we obtain the following algebraic equation for the Green tensorĜ kj (k) in Fourier space where are symmetric matrices. If we further define the unit vector in Fourier space then (20) becomes: or equivalently, in matrix notation, 3 The Fourier transform and its inverse are defined as, respectively [Vladimirov (1971)]: For a real-valued function, the inverse Fourier transform is Stability of the differential operator L M requires that the matrix C(κ) + k 2 D(κ) be positive definite. Since this requirement must hold for all k and κ, then the matrices C(κ) and D(κ) must be individually positive definite. Under the assumption that C(κ) and D(κ) are symmetric positive definite (SPD) matrices, the solution of (25) in Fourier space clearly reads: The three-dimensional Green tensor in real space is obtained by inverse Fourier transform of Eq. (26). It reads: In Eq. (27), dV = k 2 dk dω indicates the volume element in Fourier space, and dω is an elementary solid angle on the unit sphere S. Our objective now is to obtain an alternative expression of the matrix inverse [C(κ) + k 2 D(κ)] −1 which allows us to carry out the the k-integral analytically. By doing so, the non-singular nature of the Green tensor at the origin is revealed. We start by observing that, by virtue of its SPD character, the matrix C(κ) admits the following eigen-decomposition where R(κ) is the orthogonal matrix of the eigenvectors of C(κ), while V 2 (κ) is the diagonal matrix of positive eigenvalues of C(κ). Moreover, the matrix is also SPD. Using (29), let us consider the following identity: where is a SPD matrix with units of length squared. With this decomposition, the Green tensor in Fourier space becomeŝ while in real space we obtain The Green tensor of Mindlin's anisotropic first strain gradient elasticity 7 In order to carry out the k-integral, we make use of the following matrix identity: 4 With this identity, the Green tensor takes the form Next, Eq. (36) is further simplified noting that the integration kernel is an even function of κ. Therefore, the integral over the unit sphere S is twice the integral over a hemisphere.
At the origin, any arbitrary hemisphere H can be chosen, and the Green tensor assumes the value This noteworthy result shows that the Green tensor is non-singular at the origin, in contrast to classical elasticity. Away from the origin, we can choose the hemisphere having the direction x as the zenith. This is a convenient choice because all points κ on such a where D 2 (κ) = diag λ 2 i (κ) is the diagonal matrix of the positive eigenvalues of Λ 2 (κ), and Q(κ) is the orthogonal matrix of its eigenvectors. With this observation, we immediately obtain With the help of the definite integral 3.767 in [Gradshteyn and Ryzhik (2007)], we obtain In the last step we have used the property that the matrix exponential is an isotropic tensor-valued function of its argument. hemisphere satisfy the condition κ · x ≥ 0. This hemisphere can be parameterized by the zenith angle θ and the azimuth angle φ, as shown in Fig. 1. In this reference system, the unit vector κ can be expressed as κ(θ, φ) = sin θ cos φê 1 + sin θ sin φê 2 + cos θê 3 , whereê 3 = x/x. Finally, letting q = cos θ, the elementary solid angle becomes and κ(q, φ) = 1 − q 2 cos φê 1 + 1 − q 2 sin φê 2 + qê 3 .
Therefore the Green tensor of the anisotropic Mindlin differential operator of first order finally becomes

The first two gradients of the Green tensor
The first two gradients of the Green tensor are computed directly by differentiation of (36). The first gradient reads The Green tensor of Mindlin's anisotropic first strain gradient elasticity 9 In components this is: Note that, because of the presence of the sign function, the gradient of the Green tensor is finite but discontinuous at the origin. From a computational perspective, it is more convenient to express this result in reference system of Fig. 1. Doing so we find the alternative representation The second gradient of the Green tensor reads Letting n(φ) = κ(π/2, φ) be a unit vector on the equatorial plane κ · x = 0, we finally obtain Note that the second gradient of the Green tensor is singular at the origin.

The classical limit
It is now shown that Green tensor (36) converges to the classical Green tensor G 0 [Lifshitz and Rosenzweig (1947), Synge (1957)] when the field point x is sufficiently far from the origin compared to the characteristic length scales, that is when where λ i is an eigenvalue of Λ, and i = 1, 2, 3. This important property guarantees that the non-singular Green tensor (41) regularizes the classical anisotropic Green tensor in the far field. Moreover, as a special case satisfying condition (47), the classical Green tensor G 0 is also recovered in the limit of vanishing tensor of strain gradient coefficients D. The classical Green tensor G 0 is readily recovered if we consider the limit 5 wherex = x/x and I is the identity tensor. In fact, the substitution of (48) into (36) yields Here we used again the notation n(φ) = κ(π/2, φ) to indicate a unit vector on the equatorial plane κ · x = 0. Note that the span of integration can be reduced to the range 0 ≤ φ ≤ π using the symmetry C −1 (n) = C −1 (−n).

Special cases
In this section we show that the Green tensor (36) generalizes other results obtained in the literature.

The weakly non-local Green tensor G NL
Lazar and Po b] have considered a simplified strain gradient elasticity theory under the assumption a framework which was named Mindlin's strain gradient elasticity with weak non-locality because of its relation to non-local theories , Lazar and Agiasofitou (2011)]. The Green tensor (36) recovers our previous result as a special case. In fact, under the previous assumption, we have 5 Using the eigen-decomposition (34): and Therefore the Green tensor becomes which is the expression given in b].

The Green tensor of anisotropic gradient elasticity of Helmholtz type G H
An even simpler theory, named Mindlin's gradient elasticity of Helmholtz type, has been proposed by ]. The theory is characterized by only one gradient length scale parameter ℓ, which renders the tensor L diagonal: The non-singular Green tensor of this theory is obtained by substituting (53) in (52), thus yielding which coincides with the expression given in ].

The isotropic Green tensor G I
The isotropic tensor C has components where λ and µ are the Lamé constants. On the other hand, the isotropic tensor D reads where a 1 , a 2 , a 3 , a 4 , a 5 are the gradient parameters in isotropic Mindlin's first strain gradient elasticity theory [Mindlin (1964)] (see also , b]). Therefore, the matrices C(κ) and D(κ) become, respectively D ik (κ) = 2(a 1 + a 2 + a 3 + a 4 + a 5 )κ i κ k The two characteristic lengths ℓ 1 and ℓ 2 introduced above are defined as The Green tensor of Mindlin's anisotropic first strain gradient elasticity 13 Owing to the special structure 6 of C(κ) and D(κ), the following results are easily obtained: The matrix Λ −1 admits the eigenvalue 1/ℓ 1 , corresponding to the eigenvector v 1 = κ. The degenerate eigenvalue 1/ℓ 2 has multiplicity two, corresponding to two arbitrary eigenvectorsv 2 andv 3 perpendicular to κ. Choosing such eigenvectors to be mutually orthogonal, the matrix Λ −1 admits the eigen de- is an orthogonal matrix whose columns are the eigenvectors of Λ −1 , and is the diagonal matrix of its eigenvalues. This special form of Q yields the identity Using these results in (36), we obtain 6 Consider a matrix A with structure If a > b > 0, then the matrix is SPD, and a unique SPD square root of A ij exists with form Moreover, the inverse of A ij reads Because they form an orthonormal basis, the three eigenvectors satisfy the identityv 1 ⊗v 1 +v 2 ⊗v 2 +v 3 ⊗v 3 = I, hence we have The integral over the unit sphere is carried out using the relation where the scalar function A(x, ℓ) is The scalar function A(x, ℓ) can be regarded as a regularized distance function in the sense that A(x, ℓ) tends to x when x/ℓ ≫ 1, while it smoothly approaches to 2ℓ for small x, as shown in Fig. 2. By sake of (71), the Green tensor finally becomes: This result can also be obtained by direct inverse Fourier transform of (26), as shown in Appendix A. A more detailed analysis of the isotropic Green tensor (73) can be found in ].

A comparison with Molecular Statics: The Kelvin problem
In this section, we compare the Green tensor obtained from Mindlin's strain gradient elastic theory to that obtained from an atomistic system. This study was carried out using Minimol ] which is a KIMcompliant molecular dynamics (MD) and molecular statics (MS) program. The Open Knowledgebase of Interatomic Models (KIM) is a project focused on creating standards for atomistic simulations including an application programming interface (API) for information exchange between atomistic simulation codes and interatomic potentials ),Tadmor et al. (2013]. We choose face-centered-cubic Aluminum and Copper for this comparison, and consider the following two interatomic potentials: the modified-embeddedatom-method (MEAM) by [Lee (2001)], and the embedded-atom-potential by ], which are archived in the OpenKIM repository. Elastic and gradient-elastic constants for these potentials were computed using the method described in [Admal et al (2016)], and they are available on the KIM repository [Lee (2001), ]. For convenience, the values of the independent elastic and gradient-elastic constants are reported in table 1. These components are used to populate the elastic tensors C and D The Green tensor of Mindlin's anisotropic first strain gradient elasticity 15 [Admal et al (2016), Auffray et al. (2013)]. The Voigt structure of the resulting tensors C and D is shown in Fig. 3.
The atomistic system is constructed by stacking together 15 × 15 × 15 unit cells resulting in 13500 atoms. A force of 0.0116 eV/Å in the x 1 direction is imposed on the central atom of the system, and displacement boundary conditions are imposed on five layers of atoms close to the boundary using the classical solution given in Eq. (49). The padding atoms thickness is 0.15 times the size of the box. A MS simulation is carried out using the abovementioned boundary conditions resulting in a deformed crystal. The resulting displacement field normalized with respect to the force on the central atom yields the atomistic Green tensor component fields.
Simulation results are shown in Fig. 4, where we compare the Green tensor components G 11 (x 1 , 0, 0) and G 22 (x 1 , 0, 0). Despite the fact that these potentials were never fitted to gradient-elastic constants, it can be observed that the analytical predictions are in good agreement with MS calculations, with a maximum error at the origin in the order of 5-30%, depending on the potential used. It should be noted that, compared to the EAM potential, the MEAM potential better compares to the analytical results, possibly as a result of artifacts in gradient-elastic constants evaluated by EAM potentials [Admal et al (2016)].

Summary and Conclusions
In this paper we have derived an expression for the Green tensor of Mindlin's anisotropic strain gradient elasticity, which possesses up to 21 elastic constants and 171 gradient elastic constants in the general case of triclinic media. The Green tensor is found in terms of a matrix kernel integrated over the unit   sphere in Fourier space. Such representation is similar to that of the classical anisotropic Green tensor, which requires integration over the equatorial plane of the unit sphere. In contrast to its classical counterpart, however, the Green tensor of Mindlin's anisotropic strain gradient elasticity is non-singular at the origin, while its gradient is finite but discontinuous at the origin. It is shown that the non-singular Green tensor converges to the classical tensor a few characteristic lengths away from the origin. Therefore, the Green tensor of Mindlin's first strain gradient elasticity can be regarded as a physical regularization of the classical anisotropic Green tensor. Moreover, existing expressions of the Green tensor found in the literature are recovered as special cases. Because the Green tensor regularizes its classical counterpart without unphysical singularities, it offers a more realistic description of near-core elastic fields of defects in micro-mechanics, and it provides more accurate boundary conditions for atomistic and ab-initio energy-minimization calculations. As an illustrative example, we have computed the displacement field induced by a concentrated force acting at the origin (Kelvin problem), and compared the analytical pre- dictions to atomistic calculations when the elastic and gradient-elastic moduli are consistently derived from the interatomic potentials. Despite the fact that these potentials were not fitted to gradient-elastic constants, it is shown that the analytical predictions are in good agreement with MS calculations, with a maximum error at the origin in the order of 5-30%, depending on the potential used.

Declarations
Availability of data and materials. Elastic and gradient-elastic material constants used to obtain the results in section 5 are freely available as part of the Open Knowledgebase of Interatomic Models (KIM).