Hill polarization tensors

The Eshelby problem reduces to computing one object, the Hill polarization tensor $\mathbb{P}(\boldsymbol{A},\mathbb{C})$. This page gives its closed forms.

The structure to keep in mind: $\mathbb{P}$ factors into a purely geometric part, which depends on the ellipsoid alone, and a purely material part, which depends on the reference moduli alone. The geometric part is a set of Newton-potential integrals; every shape — triaxial ellipsoid, spheroid, sphere, infinite cylinder — is one column of the same two tables.

This page follows the appendix Hill polarization tensors of the Echoes manual; expressions, conventions and bibliography are aligned on it. Extensions specific to MeanFieldHom are flagged as such.

Newton-potential integrals

Three integrals over the unit sphere, depending on $\boldsymbol{A}$ only, factor every analytical Hill formula:

\[\boldsymbol{I}^{\boldsymbol{A}} = \frac{\det\boldsymbol{A}}{4\pi} \int_{\|\underline{\xi}\|=1} \frac{\underline{\xi}\otimes\underline{\xi}} {\|\boldsymbol{A}\cdot\underline{\xi}\|^{3}}\,\mathrm{d}S_{\xi} = \frac{1}{4\pi} \int_{\|\underline{\zeta}\|=1} \frac{(\boldsymbol{A}^{-1}\!\cdot\underline{\zeta})\otimes (\boldsymbol{A}^{-1}\!\cdot\underline{\zeta})} {\|\boldsymbol{A}^{-1}\!\cdot\underline{\zeta}\|^{2}}\,\mathrm{d}S_{\zeta}\]

\[\mathbb{U}^{\boldsymbol{A}} = \frac{\det\boldsymbol{A}}{4\pi} \int_{\|\underline{\xi}\|=1} \frac{\underline{\xi}\otimes\underline{\xi}\otimes \underline{\xi}\otimes\underline{\xi}} {\|\boldsymbol{A}\cdot\underline{\xi}\|^{3}}\,\mathrm{d}S_{\xi}\]

\[\mathbb{V}^{\boldsymbol{A}} = \frac{\det\boldsymbol{A}}{4\pi} \int_{\|\underline{\xi}\|=1} \frac{\underline{\xi}\stackrel{s}{\otimes}\boldsymbol{1} \stackrel{s}{\otimes}\underline{\xi}} {\|\boldsymbol{A}\cdot\underline{\xi}\|^{3}}\,\mathrm{d}S_{\xi} = \frac{\boldsymbol{1}\stackrel{s}{\boxtimes}\boldsymbol{I}^{\boldsymbol{A}} + \boldsymbol{I}^{\boldsymbol{A}}\stackrel{s}{\boxtimes}\boldsymbol{1}}{2}\]

The two parametrizations of $\boldsymbol{I}^{\boldsymbol{A}}$ are related by the bijection of the unit sphere onto itself $\underline{\zeta}\mapsto\underline{\xi} = \boldsymbol{A}^{-1}\!\cdot\underline{\zeta}/ \|\boldsymbol{A}^{-1}\!\cdot\underline{\zeta}\|$, whose surface-element identity is

\[\mathrm{d}S_{\zeta} = \frac{\det\boldsymbol{A}}{\|\boldsymbol{A}\cdot\underline{\xi}\|^{3}}\, \mathrm{d}S_{\xi}.\]

An intrinsic proof is given in [7], as an alternative to the component reasoning of [8]. Both forms are useful: the $\underline{\xi}$ form makes the geometry explicit, the $\underline{\zeta}$ form is the one that degenerates cleanly in the cylinder limit below.

In MeanFieldHom these are tens_IA, tens_UA and tens_VA.

Principal coefficients

$\boldsymbol{A}$ and $\boldsymbol{I}^{\boldsymbol{A}}$ share their eigenvectors, so

\[\boldsymbol{I}^{\boldsymbol{A}} = \sum_{i=1}^{3} I_i^{\boldsymbol{A}}\, \underline{e}^{\boldsymbol{A}}_i\otimes\underline{e}^{\boldsymbol{A}}_i .\]

The coefficients $I_i^{\boldsymbol{A}}$, identified with Newton-potential integrals [9], [3], [10], and the secondary coefficients $I_{ij}^{\boldsymbol{A}}$ admit closed forms in every symmetry class:

Here $\mathcal{F} = \mathcal{F}(\theta,\kappa)$ and $\mathcal{E} = \mathcal{E}(\theta,\kappa)$ are the incomplete elliptic integrals of the first and second kind [11], of amplitude and parameter

\[\theta = \arcsin\sqrt{1-\frac{c^{2}}{a^{2}}}, \qquad \kappa = \sqrt{\frac{a^{2}-b^{2}}{a^{2}-c^{2}}}.\]

The secondary coefficients follow from the $I_i^{\boldsymbol{A}}$ by

\[I_{ij}^{\boldsymbol{A}} = \frac{I_j^{\boldsymbol{A}}-I_i^{\boldsymbol{A}}} {\rho_i^{2}-\rho_j^{2}} \quad (i\ne j), \qquad I_{ii}^{\boldsymbol{A}} = \frac{1}{3}\left(\frac{1}{\rho_i^{2}} - \sum_{j\ne i} I_{ij}^{\boldsymbol{A}}\right),\]

except where the denominator degenerates — each symmetry class then has its own regular expression:

with $I_{ij}^{\boldsymbol{A}} = I_{ji}^{\boldsymbol{A}}$. The circular cylinder $b=c$ gives $I_2^{\boldsymbol{A}}=I_3^{\boldsymbol{A}}=\tfrac{1}{2}$ and $I_{22}^{\boldsymbol{A}}=I_{33}^{\boldsymbol{A}}=I_{23}^{\boldsymbol{A}} =\tfrac{1}{4c^{2}}$.

Identities

Always satisfied, and useful as numerical checks [3]:

\[\sum_i I_i^{\boldsymbol{A}} = 1, \qquad 3\,I_{ii}^{\boldsymbol{A}} + \sum_{j\ne i} I_{ij}^{\boldsymbol{A}} = \frac{1}{\rho_i^{2}}, \qquad 3\,\rho_i^{2}\,I_{ii}^{\boldsymbol{A}} + \sum_{j\ne i}\rho_j^{2}\,I_{ij}^{\boldsymbol{A}} = 3\,I_i^{\boldsymbol{A}}.\]

Components of $\mathbb{U}^{\boldsymbol{A}}$ and $\mathbb{V}^{\boldsymbol{A}}$

Both are orthotropic along the ellipsoid axes. In the principal frame [7], [2]:

\[U^{\boldsymbol{A}}_{iiii} = \tfrac{3}{2}\bigl(I_i^{\boldsymbol{A}}-\rho_i^{2}I_{ii}^{\boldsymbol{A}}\bigr), \qquad U^{\boldsymbol{A}}_{iijj} = U^{\boldsymbol{A}}_{ijij} = U^{\boldsymbol{A}}_{ijji} = \tfrac{1}{2}\bigl(I_j^{\boldsymbol{A}}-\rho_i^{2}I_{ij}^{\boldsymbol{A}}\bigr) = \tfrac{1}{2}\bigl(I_i^{\boldsymbol{A}}-\rho_j^{2}I_{ij}^{\boldsymbol{A}}\bigr)\]

\[V^{\boldsymbol{A}}_{iiii} = I_i^{\boldsymbol{A}}, \qquad V^{\boldsymbol{A}}_{ijij} = V^{\boldsymbol{A}}_{ijji} = \tfrac{1}{4}\bigl(I_i^{\boldsymbol{A}}+I_j^{\boldsymbol{A}}\bigr) \qquad (i\ne j).\]

Sphere $\boldsymbol{A}=\boldsymbol{1}$:

\[\mathbb{U}^{\boldsymbol{1}} = \tfrac{1}{3}\mathbb{J} + \tfrac{2}{15}\mathbb{K}, \qquad \mathbb{V}^{\boldsymbol{1}} = \tfrac{1}{3}\mathbb{I}.\]

Infinite elliptic cylinder $a\to\infty$, axis $\underline{e}^{\boldsymbol{A}}_1$, transverse semi-axes $b\ge c$ ([8], §11.22) — substituting the cylinder column above:

\[\mathrm{Mat}\bigl(\mathbb{U}^{\mathrm{cyl}}\bigr) = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \frac{c(b+2c)}{2(b+c)^{2}} & \frac{bc}{2(b+c)^{2}} & 0 & 0 & 0\\ 0 & \frac{bc}{2(b+c)^{2}} & \frac{b(2b+c)}{2(b+c)^{2}} & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{bc}{(b+c)^{2}} & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \quad \mathrm{Mat}\bigl(\mathbb{V}^{\mathrm{cyl}}\bigr) = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \frac{c}{b+c} & 0 & 0 & 0 & 0\\ 0 & 0 & \frac{b}{b+c} & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{1}{2} & 0 & 0\\ 0 & 0 & 0 & 0 & \frac{b}{2(b+c)} & 0\\ 0 & 0 & 0 & 0 & 0 & \frac{c}{2(b+c)} \end{pmatrix}\]

in Kelvin–Mandel storage and in the frame $(\underline{e}^{\boldsymbol{A}}_i)$. The vanishing first row and column is the signature of the infinite cylinder: no polarization is transmitted along its axis. For the circular cylinder $b=c$ the non-zero entries reduce to $U^{\mathrm{cyl}}_{2222}=U^{\mathrm{cyl}}_{3333}=\tfrac{3}{8}$, $U^{\mathrm{cyl}}_{2233}=\tfrac{1}{8}$, $\bigl[\mathrm{Mat}(\mathbb{U}^{\mathrm{cyl}})\bigr]_{44}=\tfrac{1}{4}$, and $\mathrm{Mat}(\mathbb{V}^{\mathrm{cyl}})$ to the diagonal $\bigl(0,\tfrac{1}{2},\tfrac{1}{2},\tfrac{1}{2},\tfrac{1}{4},\tfrac{1}{4}\bigr)$.

Hill tensor in elasticity

Arbitrary anisotropy

\[\mathbb{P}(\boldsymbol{A},\mathbb{C}) = \frac{1}{4\pi} \int_{\|\underline{\zeta}\|=1} (\boldsymbol{A}^{-1}\!\cdot\underline{\zeta})\stackrel{s}{\otimes} \Bigl((\boldsymbol{A}^{-1}\!\cdot\underline{\zeta})\cdot\mathbb{C} \cdot(\boldsymbol{A}^{-1}\!\cdot\underline{\zeta})\Bigr)^{-1} \stackrel{s}{\otimes}(\boldsymbol{A}^{-1}\!\cdot\underline{\zeta}) \,\mathrm{d}S_{\zeta}\]

\[= \frac{\det\boldsymbol{A}}{4\pi} \int_{\|\underline{\xi}\|=1} \frac{\underline{\xi}\stackrel{s}{\otimes} \bigl(\underline{\xi}\cdot\mathbb{C}\cdot\underline{\xi}\bigr)^{-1} \stackrel{s}{\otimes}\underline{\xi}} {\|\boldsymbol{A}\cdot\underline{\xi}\|^{3}}\,\mathrm{d}S_{\xi}\]

[5], see also [8]. Inverting the acoustic tensor $\underline{\xi}\cdot\mathbb{C}\cdot\underline{\xi}$ pointwise is the source of all the computational work: in general no closed form exists and one resorts to numerical cubature [12], [13], [14]. MeanFieldHom offers three algorithm traits, mirroring the Echoes NUMINT / RESIDUES options:

• DECUHR — the surface integral is evaluated by the adaptive cubature for singular integrands of [15]. ForwardDiff-safe. Selected by :auto when its (weak-dependency) extension is loaded.
• NestedQuadGK — nested adaptive 1-D quadrature, type-generic and always available; the :auto choice otherwise, and the one used for Dual, Complex and symbolic coefficients whatever else is loaded.
• Residue — the inner $\varphi$ integral is reduced to a sum of residues by the Cauchy theorem, leaving a single 1-D quadrature [14]. The fastest of the three (~4 ms against ~11 ms and ~31 ms on a triaxial ellipsoid), but reachable on an explicit method = :residues only: it is Float64-only (the polynomial root finder it needs is not differentiable by ForwardDiff), and its acoustic polynomial degenerates when the reference is anisotropic in type while isotropic in value — returning NaN there instead of a number. That reference is what the differential and self-consistent schemes feed back at their first step, hence the choice of a cubature as the default.

Analytical paths exist in the literature for further anisotropy classes ([16], [17], [18], [19]) and are not all implemented yet.

Isotropic matrix — the shape/moduli factorization

With bulk modulus $k$, shear modulus $\mu$ and first Lamé parameter $\lambda = k-2\mu/3$, so that $\mathbb{C} = 3k\,\mathbb{J}+2\mu\,\mathbb{K} = 3\lambda\,\mathbb{I}+2\mu\,\mathbb{K}$, the general expression collapses to [5]

\[\boxed{\; \mathbb{P}\bigl(\boldsymbol{A},\,3\lambda\,\mathbb{I}+2\mu\,\mathbb{K}\bigr) = \frac{1}{\lambda+2\mu}\,\mathbb{U}^{\boldsymbol{A}} + \frac{1}{\mu}\,\bigl(\mathbb{V}^{\boldsymbol{A}}-\mathbb{U}^{\boldsymbol{A}}\bigr). \;}\]

This is the factorization announced at the top of the page: shape and orientation on one side ($\mathbb{U}^{\boldsymbol{A}}$, $\mathbb{V}^{\boldsymbol{A}}$), reference moduli on the other. Each shape column of the tables above therefore yields a closed-form $\mathbb{P}$ at once.

For a sphere, substituting $\mathbb{U}^{\boldsymbol{1}}$ and $\mathbb{V}^{\boldsymbol{1}}$ gives the classical Eshelby result

\[\mathbb{P}\bigl(\boldsymbol{1},\,3k\,\mathbb{J}+2\mu\,\mathbb{K}\bigr) = \frac{1}{3k+4\mu}\left(\mathbb{J} + \frac{3(k+2\mu)}{5\mu}\,\mathbb{K}\right).\]

For an infinite cylinder, substituting $\mathrm{Mat}(\mathbb{U}^{\mathrm{cyl}})$ and $\mathrm{Mat}(\mathbb{V}^{\mathrm{cyl}})$ gives a closed form with $P^{\mathrm{cyl}}_{1jkl}\equiv 0$, i.e. no polarization along the axis ([8], §11.22).

Implementation: src/Elasticity/hill_3d_iso.jl and src/Elasticity/hill_3d_cylinder_iso.jl, selected by method = :auto when $\mathbb{C}_0$ is a TensISO.

Transversely isotropic matrix coaxial with a spheroid

When the matrix is transversely isotropic and its symmetry axis is parallel to the spheroid axis, a fully analytical path exists [2]. The Hill tensor is transversely isotropic too, hence five Walpole coefficients (see Notation — there is no $P_4$ because $\mathbb{P}$ is major-symmetric):

\[\mathbb{P} = P_1\,\mathbb{W}_1 + P_2\,\mathbb{W}_2 + P_3\,(\mathbb{W}_3+\mathbb{W}_4) + P_5\,\mathbb{W}_5 + P_6\,\mathbb{W}_6 .\]

The five coefficients are closed-form combinations of acosh and complex square roots (equations 53–58 of [2]), depending on the aspect ratio $\omega$ (axial / transverse) and the five independent constants $(C_{1111}, C_{1122}, C_{1133}, C_{3333}, C_{2323})$.

The dispatcher routes a TensTI{4} matrix combined with a coaxial Ellipsoid{3, Spherical|Prolate|Oblate} to this path by default; coaxiality is detected by _ti_coaxial(C₀, ell). Non-coaxial spheroids and triaxial ellipsoids fall back to the anisotropic default, i.e. a cubature. Implementation: src/Elasticity/hill_3d_ti_coaxial.jl.

Anisotropic matrix, cylinder limit

Cylinder is a first-class inclusion type here (extension over Echoes). For an arbitrarily anisotropic matrix the Masson residue algorithm does not apply: it rests on the six complex roots of the acoustic polynomial along $\underline{\xi}_3$, and at the cylinder limit one root escapes to infinity, degenerating the polynomial.

The $\underline{\zeta}$ form of the Willis integral degenerates cleanly instead. The axial component of $\underline{\zeta}$ vanishes identically, so the surface integral collapses to a single quadrature over the transverse unit circle:

\[\mathbb{P}^{\mathrm{cyl}} = \frac{1}{2\pi}\int_{0}^{2\pi} \underline{\zeta}\stackrel{s}{\otimes} \bigl(\underline{\zeta}\cdot\mathbb{C}\cdot\underline{\zeta}\bigr)^{-1} \stackrel{s}{\otimes}\underline{\zeta} \,\mathrm{d}\varphi, \qquad \underline{\zeta}(\varphi) = \Bigl(0,\ \frac{\cos\varphi}{b},\ \frac{\sin\varphi}{c}\Bigr).\]

This is a one-dimensional QuadGK integral, and it stays ForwardDiff-compatible. Calling hill_tensor(Cylinder(…), C₀; method = :residues) falls back to it silently. Implementation: src/Elasticity/hill_3d_cylinder_aniso.jl (CylinderQuadrature trait). The in-plane components coincide with the solution of the 2-D plane-strain problem — the cylinder is the 3-D realization of the 2-D ellipse.

2-D plane strain

Plane strain is handled directly (extension over Echoes), integrating over the unit circle $\underline{\xi}\in S^{1}$ with a $1/(2\pi)$ prefactor in place of $1/(4\pi)$. The isotropic case is analytical; the anisotropic one uses the Masson residue reduction on the line integral.

Hill tensor in conductivity

Arbitrary anisotropy — closed form

For a conductivity tensor $\boldsymbol{K}$ [5]:

\[\boldsymbol{P}(\boldsymbol{A},\boldsymbol{K}) = \frac{\det\boldsymbol{A}}{4\pi} \int_{\|\underline{\xi}\|=1} \frac{\underline{\xi}\otimes\underline{\xi}} {(\underline{\xi}\cdot\boldsymbol{K}\cdot\underline{\xi})\, \|\boldsymbol{A}\cdot\underline{\xi}\|^{3}}\,\mathrm{d}S_{\xi}.\]

Unlike the order-4 case, this integral has a closed form for any matrix anisotropy. Since $\boldsymbol{K}$ is symmetric positive definite it has a square root, $\boldsymbol{K}^{1/2} = \sum_i\sqrt{K_i}\, \underline{e}^{\boldsymbol{K}}_i\otimes\underline{e}^{\boldsymbol{K}}_i$, and the denominator can be absorbed into the change of variable $\underline{\zeta}\mapsto\boldsymbol{K}^{1/2}\!\cdot\boldsymbol{A}^{-1}\!\cdot \underline{\zeta}$, which turns the anisotropic problem into an isotropic one for a fictitious ellipsoid of shape tensor $\boldsymbol{A}\cdot\boldsymbol{K}^{-1/2}$:

\[\boxed{\; \boldsymbol{P}(\boldsymbol{A},\boldsymbol{K}) = \boldsymbol{K}^{-1/2}\cdot \boldsymbol{P}(\boldsymbol{A}\cdot\boldsymbol{K}^{-1/2},\boldsymbol{1})\cdot \boldsymbol{K}^{-1/2} = \boldsymbol{K}^{-1/2}\cdot \boldsymbol{I}^{\boldsymbol{A}\cdot\boldsymbol{K}^{-1/2}}\cdot \boldsymbol{K}^{-1/2}. \;}\]

This is the transformation derivation of [20]; an equivalent Green's-function derivation is given in [21]. Since $\boldsymbol{A}\cdot\boldsymbol{K}^{-1/2}$ need not be symmetric, the fictitious semi-axes and principal directions are obtained by diagonalizing $\boldsymbol{K}^{-1/2}\cdot\boldsymbol{A}^{\!T}\!\cdot\boldsymbol{A}\cdot \boldsymbol{K}^{-1/2}$.

Isotropic matrix — immediate

If $\boldsymbol{K} = K\,\boldsymbol{1}$ the prefactor comes straight out:

\[\boldsymbol{P}(\boldsymbol{A}, K\,\boldsymbol{1}) = \frac{\boldsymbol{I}^{\boldsymbol{A}}}{K}.\]

For a sphere, $\boldsymbol{I}^{\boldsymbol{1}} = \tfrac{1}{3}\boldsymbol{1}$ gives $\boldsymbol{P} = \tfrac{1}{3K}\boldsymbol{1}$ and $\boldsymbol{s} = \tfrac{1}{3}\boldsymbol{1}$ — independent of $K$. Implementation: src/Conductivity/hill_order2_3d.jl.

Dispatch

Entry point hill_tensor; shape tensor via shape_tensor; geometric auxiliaries via tens_IA, tens_UA, tens_VA.

Cylinder shape traits: CircularCylindrical when $b=c$ (transversely isotropic response, returned as TensTI{4} with axis $\underline{e}^{\boldsymbol{A}}_1$) and EllipticCylindrical when $b>c$ (orthotropic, returned as TensOrtho). Practical usage is covered in the manual page Cylindrical inclusions.