Displacement field around a rigid sphere in a compressible elastic environment, corresponding higher-order Fax\'en relations, as well as higher-order displaceability and rotateability matrices

TL;DR

This paper presents an efficient method for calculating the displacement field around a rigid sphere in a compressible elastic medium, derives higher-order Faxén relations, and computes interaction matrices up to sixth order, with applications to hydrodynamics.

Contribution

It introduces a novel approach for displacement field calculation, explicitly derives higher-order Faxén relations, and computes displaceability and rotateability matrices up to sixth order.

Findings

Higher-order Faxén relations are explicitly derived for compressible media.

Displaceability and rotateability matrices are calculated up to sixth order.

The methods are applicable to both elastic media and low-Reynolds-number hydrodynamics.

Abstract

An efficient route to the displacement field around a rigid spherical inclusion in an infinitely extended homogeneous elastic medium is presented in a slightly alternative way when compared to some common textbook methods. Moreover, two Fax\'en relations of next-higher order beyond the stresslet are calculated explicitly for compressible media. They quantify higher-order moments involving the force distribution on rigid particles in a deformed elastic medium. Additionally, the displaceability and rotateability matrices are calculated up to (including) sixth order in inverse particle separation distance. These matrices describe the interactions mediated between the rigid embedded particles by the elastic environment. All methods and results can formally be transferred to the corresponding case of incompressible hydrodynamic low-Reynolds-number Stokes flow by considering the limit of an incompressible environment. The roles of compressibility of the embedding medium and of the here additionally derived higher-order contributions are highlighted by some selected example configurations.

Figures

Tables

Table 1. Table 1: Equation ( 17 ) for N = 0 , 1 , 2 , 3 𝑁 0 1 2 3 N=0,1,2,3 (from top to bottom).

	${\int }_{\partial V}d{S}^{\prime }$	$=$	$4\pi a^2$
	${\int }_{\partial V}d{S}^{\prime }{r}_j^{\prime }{r}_k^{\prime }$	$=$	$\frac{4\pi a^4}{3}{\delta }_{jk}$
	${\int }_{\partial V}d{S}^{\prime }{r}_j^{\prime }{r}_k^{\prime }{r}_l^{\prime }{r}_m^{\prime }$	$=$	$\frac{4\pi a^6}{15}\left[{\delta }_{jk}{\delta }_{lm}+{\delta }_{jl}{\delta }_{km}+{\delta }_{jm}{\delta }_{kl}\right]$
	${\int }_{\partial V}d{S}^{\prime }{r}_j^{\prime }{r}_k^{\prime }{r}_l^{\prime }{r}_m^{\prime }{r}_n^{\prime }{r}_o^{\prime }$	$=$	$\frac{4\pi a^8}{105}[{\delta }_{jk}({\delta }_{lm}{\delta }_{no}+{\delta }_{ln}{\delta }_{mo}+{\delta }_{lo}{\delta }_{mn})$
			$+{\delta }_{jl}\left({\delta }_{km}{\delta }_{no}+{\delta }_{kn}{\delta }_{mo}+{\delta }_{ko}{\delta }_{mn}\right)$
			$+{\delta }_{jm}\left({\delta }_{kl}{\delta }_{no}+{\delta }_{kn}{\delta }_{lo}+{\delta }_{ko}{\delta }_{ln}\right)$
			$+{\delta }_{jn}\left({\delta }_{kl}{\delta }_{mo}+{\delta }_{km}{\delta }_{lo}+{\delta }_{ko}{\delta }_{lm}\right)$
			$+{\delta }_{jo}({\delta }_{kl}{\delta }_{mn}+{\delta }_{km}{\delta }_{ln}+{\delta }_{kn}{\delta }_{lm})]$
	${\int }_{\partial V}d{S}^{\prime }{r}_j^{\prime }{r}_k^{\prime }{r}_l^{\prime }{r}_m^{\prime }{r}_n^{\prime }{r}_o^{\prime }{r}_p^{\prime }{r}_q^{\prime }$	$=$	$\frac{4\pi a^{10}}{945}\left[{\delta }_{jk}{\delta }_{lm}{\delta }_{no}{\delta }_{pq}+\text{all permutations}\right]$

Table 2. Table 2: Derivatives of the Green’s tensor, see Eq. ( 5 ).

	$16\pi \left(1-\nu \right)\mu {\nabla }_k{G}_{ij}$	$=$	$\frac{1}{r^3}\left[-\left(3-4\nu \right)r_k{\delta }_{ij}+r_j{\delta }_{ik}+r_i{\delta }_{jk}\right]-\frac{3}{r^5}{r}_i{r}_j{r}_k$
	$16\pi \left(1-\nu \right)\mu {\nabla }^2{G}_{ij}$	$=$	$\frac{2}{r^3}{\delta }_{ij}-\frac{6}{r^5}{r}_i{r}_j$

Table 3. Table 3: Here, Ψ k ​ l := ( ϑ ^ k ​ ϑ ^ l + φ ^ k ​ φ ^ l ) / 2 assign subscript Ψ 𝑘 𝑙 subscript ^ italic-ϑ 𝑘 subscript ^ italic-ϑ 𝑙 subscript ^ 𝜑 𝑘 subscript ^ 𝜑 𝑙 2 \Psi_{kl}:=(\hat{\vartheta}_{k}\hat{\vartheta}_{l}+\hat{\varphi}_{k}\hat{\varphi}_{l})/2 . For k = l 𝑘 𝑙 k=l , Ψ k ​ k = 1 subscript Ψ 𝑘 𝑘 1 \Psi_{kk}=1 and the results from Tab.  1 are recovered. In the first term on the right-hand side of the last equation, all permutations of δ i ​ j ​ δ m ​ n ​ δ o ​ p subscript 𝛿 𝑖 𝑗 subscript 𝛿 𝑚 𝑛 subscript 𝛿 𝑜 𝑝 \delta_{ij}\delta_{mn}\delta_{op} are added (i.e., δ i ​ j ​ δ m ​ o ​ δ n ​ p subscript 𝛿 𝑖 𝑗 subscript 𝛿 𝑚 𝑜 subscript 𝛿 𝑛 𝑝 \delta_{ij}\delta_{mo}\delta_{np} , δ i ​ j ​ δ m ​ p ​ δ n ​ o subscript 𝛿 𝑖 𝑗 subscript 𝛿 𝑚 𝑝 subscript 𝛿 𝑛 𝑜 \delta_{ij}\delta_{mp}\delta_{no} , and so on), which makes in total 15 terms. The second term on this right-hand side is a superposition of all permutations of δ i ​ j ​ δ k ​ m ​ δ l ​ n ​ δ o ​ p subscript 𝛿 𝑖 𝑗 subscript 𝛿 𝑘 𝑚 subscript 𝛿 𝑙 𝑛 subscript 𝛿 𝑜 𝑝 \delta_{ij}\delta_{km}\delta_{ln}\delta_{op} , except for those containing δ k ​ l subscript 𝛿 𝑘 𝑙 \delta_{kl} , thus in total 90 terms (altogether, the number of terms on the right-hand side is therefore 15 + 90 = 105 15 90 105 15+90=105 ).

	${\int }_{\partial V}dS\left(\widehat{𝐤}\right){\Psi }_{kl}$	$=$	$\frac{4\pi }{3}{\delta }_{kl}$
	${\int }_{\partial V}dS\left(\widehat{𝐤}\right){\Psi }_{kl}{\hat{k}}_m{\hat{k}}_n$	$=$	$\frac{2\pi }{15}\left[4{\delta }_{kl}{\delta }_{mn}-\left({\delta }_{km}{\delta }_{ln}+{\delta }_{kn}{\delta }_{lm}\right)\right]$
	${\int }_{\partial V}dS\left(\widehat{𝐤}\right){\hat{k}}_i{\hat{k}}_j{\Psi }_{kl}{\hat{k}}_m{\hat{k}}_n$	$=$	$\frac{2\pi }{105}\{6{\delta }_{kl}({\delta }_{ij}{\delta }_{mn}+{\delta }_{im}{\delta }_{jn}+{\delta }_{in}{\delta }_{jm})$
			$-[{\delta }_{ij}({\delta }_{km}{\delta }_{ln}+{\delta }_{kn}{\delta }_{lm})+{\delta }_{mn}({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})+{\delta }_{ik}({\delta }_{jm}{\delta }_{ln}+{\delta }_{jn}{\delta }_{lm})$
			$+{\delta }_{jk}({\delta }_{im}{\delta }_{ln}+{\delta }_{in}{\delta }_{lm})+{\delta }_{km}({\delta }_{il}{\delta }_{jn}+{\delta }_{in}{\delta }_{jl})+{\delta }_{kn}({\delta }_{il}{\delta }_{jm}+{\delta }_{im}{\delta }_{jl})]\}$
	${\int }_{\partial V}dS\left(\widehat{𝐤}\right){\hat{k}}_i{\hat{k}}_j{\Psi }_{kl}{\hat{k}}_m{\hat{k}}_n{\hat{k}}_o{\hat{k}}_p$	$=$	$\frac{2\pi }{945}\{8{\delta }_{kl}({\delta }_{ij}{\delta }_{mn}{\delta }_{op}+\text{all permutations})$
			$-[{\delta }_{ij}{\delta }_{km}{\delta }_{ln}{\delta }_{op}+\text{all permutations except for those containing }{\delta }_{kl}]\}$

Equations

$${\nabla }^2\mathbf{u}(\mathbf{r})+\frac{1}{1-2\nu }\nabla \nabla \cdot \mathbf{u}(\mathbf{r})+\frac{1}{\mu }{\mathbf{f}}_b(\mathbf{r})=0,$$

$${\nabla }^2\nabla \cdot \mathbf{u}(\mathbf{r})$$

$${\nabla }^4\mathbf{u}(\mathbf{r})$$

$$\nabla \times {\nabla }^2\mathbf{u}(\mathbf{r})$$

$$\underline{\mathbf{G}}(\mathbf{r})=\frac{1}{16\pi (1-\nu )\mu r}\left((3-4\nu )\underline{\hat{\mathbf{I}}}+\hat{\mathbf{r}}\hat{\mathbf{r}}\right),$$

$$\mathbf{u}(\mathbf{r})={\int }_{\partial V}\mathrm{d}{S}^{\prime }\underline{\mathbf{G}}(\mathbf{r}-{\mathbf{r}}^{\prime })\cdot \mathbf{f}({\mathbf{r}}^{\prime }),$$

$$G_{ij}(\mathbf{r}-{\mathbf{r}}^{\prime })=\sum _{n=0}^{\infty }\frac{(-1{)}^n}{n!}({\mathbf{r}}^{\prime }\cdot \nabla {)}^n{G}_{ij}(\mathbf{r}).$$

$$u_i(\mathbf{r})$$

$$F_j$$

$$A_{jk}$$

$$S_{jk}$$

$$M_{jkl}$$

$$N_{jklm}$$

$$T_l=-{ϵ}_{ljk}{A}_{jk},$$

$$\mathbf{T}={\int }_{\partial V}\mathrm{d}{S}^{\prime }{\mathbf{r}}^{\prime }\times \mathbf{f}({\mathbf{r}}^{\prime }).$$

$$\mathbf{u}(|\mathbf{r}|\to \infty )$$

$${\int }_{\partial V}\mathrm{d}{S}^{\prime }\prod _{k=1}^{2N}{r}_{i_k}^{\prime }$$

$$\mathbf{u}(\mathbf{r})=\left(1+\frac{a^2}{6}{\nabla }^2\right)\underline{\mathbf{G}}(\mathbf{r})\cdot \mathbf{F}.$$

$$\mathbf{u}(\mathbf{r})=-\frac{1}{2}(\mathbf{T}\times \nabla )\cdot \underline{\mathbf{G}}(\mathbf{r}).$$

$$\mathbf{u}(\mathbf{r})=\left(1+\frac{a^2}{6}{\nabla }^2\right)\underline{\mathbf{G}}(\mathbf{r})\cdot \mathbf{F}-\frac{1}{2}(\mathbf{T}\times \nabla )\cdot \underline{\mathbf{G}}(\mathbf{r}).$$

$$\mathbf{u}(\mathbf{r})$$

$$\mathbf{U}=\frac{5-6\nu }{24\pi (1-\nu )\mu a}\mathbf{F},\mathbf{\Omega }=\frac{1}{8\pi \mu a^3}\mathbf{T}.$$

$$\mathbf{U}+\mathbf{\Omega }\times \mathbf{r}={\mathbf{u}}_0(\mathbf{r})+{\int }_{\partial V}\mathrm{d}{S}^{\prime }\underline{\mathbf{G}}(\mathbf{r}-{\mathbf{r}}^{\prime })\cdot \mathbf{f}({\mathbf{r}}^{\prime }),$$

$$\mathbf{U}$$

$$\mathbf{\Omega }$$

$$\underline{\mathbf{S}}$$

$$u_{0,i}(\mathbf{r})$$

$$0$$

$${\int }_{\partial V}\mathrm{d}{S}^{\prime }{\int }_{\partial V}\mathrm{d}S{r}_k{r}_l{G}_{ij}(\mathbf{r}-{\mathbf{r}}^{\prime })f_j({\mathbf{r}}^{\prime })$$

$$0$$

$${\Gamma }_{jkl}$$

$$M_{jkl}^{\text{S}}$$

$$M_{jkl}^{\text{A}}$$

$$M_{jkl}^{\text{S}}$$

$$M_{jkl}^{\text{A}}$$

Full text

I Introduction

A basic analytical solution provided by low-Reynolds-number hydrodynamics concerns the flow past a uniformly translating or rotating rigid spherical particle in an incompressible, infinitely extended, and otherwise quiescent fluid. It is explicitly treated by many classical textbooks, see, for example, Refs. Happel and Brenner, 1981; Kim and Karrila, 1991; Dhont, 1996. The flow field that is created around the sphere under no-slip surface conditions is referred to as Stokes flow. A complementary problem concerns the situation, in which a sphere is immersed in a fluid, where it is subject to a given flow field. The sphere is then translated and rotated by the flow. In this case, to find the translational and angular velocity of the particle, the so-called Faxén relations may be employed Batchelor and Green (1972); Happel and Brenner (1981); Kim and Karrila (1991); Dhont (1996). These relations allow to determine the velocity, the angular velocity, and the stresslet that a rigid sphere acquires in an arbitrary imposed flow field. The solutions to both problems are important tools to calculate the hydrodynamic interactions between the spheres within semidilute many-sphere systems in the form of mobility matrices Kim and Karrila (1991); Dhont (1996); Mazur and Van Saarloos (1982). Refs. Happel and Brenner, 1981; Kim and Karrila, 1991; Dhont, 1996 themselves cite many related works.

The formal similarity between the incompressible low-Reynolds-number hydrodynamics problem and the corresponding problem in linear elasticity theory has been noted several times. In the hydrodynamic formulation Happel and Brenner (1981); Kim and Karrila (1991); Dhont (1996), there is only one material parameter, namely the shear viscosity $\eta$ of the incompressible fluid. In linear elasticity theory, an infinitely extended and isotropic compressible homogeneous medium is described by two material parameters Landau and Lifshitz (1986); Mura (1987). First, the shear modulus $\mu$ is the elastic analogon to the hydrodynamic shear viscosity. The second parameter determines the compressibility of the material, for instance, in the form of the Poisson ratio $\nu$. It is thus possible to obtain the hydrodynamic expressions from the linearly elastic equations by assuming the elastic material to be incompressible, which formally corresponds to taking the limit of $\nu\rightarrow 1/2$. Then, passing to the hydrodynamic case, the elastic displacement field is replaced by the hydrodynamic flow field and the shear modulus $\mu$ by the hydrodynamic viscosity $\eta$. The methods employed in the present work can in this way directly be transferred to hydrodynamic situations of incompressible fluid flows.

In the following, we address the two initially stated problems for compressible linear elasticity theory, thereby extending results in earlier works Puljiz et al. (2016); Puljiz and Menzel (2017). First, we present a straightforward, slightly alternative approach to the well-known displacement field created by a uniformly translated and rotated sphere, which we have not found in the mentioned textbooks of hydrodynamics (for the Stokes flow past a sphere) nor in corresponding literature for linear elasticity theory Mura (1987); Phan-Thien (1993). Here, the framework of the multipole expansion is used in combination with an ansatz for the force density (compare, e.g., to Refs. Kim and Karrila, 1991; Phan-Thien, 1993, which use the multipole expansion in the same context, but rather present the final expressions that satisfy the boundary conditions and the underlying equations). We use an ansatz for the surface force density on the sphere Dhont (1996) and show how this ansatz, when inserted into the multipole expansion, leads to the familiar results that satisfy the requirements on the solutions. As for the second problem, we calculate the third- and fourth-rank Faxén laws beyond the stresslet for a compressible medium, following the route outlined in Ref. Batchelor and Green, 1972 for the incompressible hydrodynamic case. We have not found these expressions in Refs. Kim and Karrila, 1991; Happel and Brenner, 1981; Dhont, 1996; Kim and Phan-Thien, 1995; Phan-Thien and Kim, 1994; Mura, 1987. As an application of these additional Faxén relations, the displaceability and rotateability matrices up to sixth order in inverse interparticle separation distance are calculated for a compressible or incompressible elastic environment. They quantify the interactions between embedded rigid spherical particles as mediated by the distorted environment, including three-body interactions. In the incompressible situation, the hydrodynamic mobility matrices as presented in Ref. Mazur and Van Saarloos, 1982 are formally recovered. The difference due to compressibility is demonstrated by exemplary induced interactions between paramagnetic spherical particles embedded in an elastic environment, which has already to lower order found application in the theoretical description of elastic magnetic gels Puljiz et al. (2016, 2018); Menzel (2017); Schopphoven et al. (2018); Puljiz and Menzel (2019).

In Sec. II the Green’s function and its multipole expansion are presented. The displacement field around a uniformly translated and rotated sphere is then calculated in Sec. III. In Sec. IV, the two Faxén laws beyond the stresslet are derived. The displacement field around a sphere that is both actively displaced and rotated by an external force and/or torque, respectively, as well as by a deformation of the surrounding medium is calculated. Then, we determine the displaceability and rotateability matrices in a compressible medium up to sixth order in inverse interparticle distance in Sec. V together with some simple examples that show the effect of compressibility and of the additionally derived higher-order contributions. We conclude in Sec. VI.

II Green’s function and multipole expansion

In linear elasticity theory, the displacement field $\mathbf{u}(\mathbf{r})$ of an unbounded, homogeneous, and isotropic elastic medium is in general described by the Navier–Cauchy equation Landau and Lifshitz (1986)

[TABLE]

with $\nu$ denoting the Poisson ratio, $\mu$ the elastic shear modulus, and $\mathbf{f}_{b}(\mathbf{r})$ a force volume density. From Eq. (1), the relations

[TABLE]

can be deduced for areas of constant or vanishing $\mathbf{f}_{b}(\mathbf{r})$ Puljiz and Menzel (2017). Equation (3) is also referred to as the biharmonic equation.

The displacement field induced by a constant point force $\mathbf{F}$ attacking at the origin, i.e., $\mathbf{f}_{b}(\mathbf{r})=\mathbf{F}\delta(\mathbf{r})$, can be written in the form $\mathbf{u}(\mathbf{r})=\mathbf{\underline{G}}(\mathbf{r})\cdot\mathbf{F}$. The corresponding Green’s function $\mathbf{\underline{G}}(\mathbf{r})$ is a symmetric second-rank tensor and reads Landau and Lifshitz (1986); Puljiz and Menzel (2017)

[TABLE]

with $\mathbf{\underline{\hat{I}}}$ denoting the identity matrix, $r=|\mathbf{r}|$, $\mathbf{\hat{r}}=\mathbf{r}/r$, and $\mathbf{\hat{r}}\mathbf{\hat{r}}$ representing a dyadic tensor. This solution satisfies Eqs. (1)–(4).

In the case of a spatially extended, displaced rigid particle in an otherwise unperturbed elastic medium, the displacement field can be written as

[TABLE]

where $\mathbf{r}^{\prime}$ are the points on the particle surface $\partial V$ and $\mathbf{f}(\mathbf{r}^{\prime})$ is the force surface density exerted by the particle onto the medium Phan-Thien (1993); Puljiz and Menzel (2017). Equation (6) follows from Eq. (5) via the superposition principle due to the linearity of the Navier–Cauchy equation. Thus, the displacement field is solely induced by the translating and/or rotating particle shell acting on the medium.

Switching to index notation and considering a particle located around the origin at $\mathbf{r}=\mathbf{0}$, the Taylor series of $\mathbf{\underline{G}}(\mathbf{r}-\mathbf{r}^{\prime})$ in $\mathbf{r}^{\prime}$ reads

[TABLE]

Inserting this expression into Eq. (6) yields

[TABLE]

with the constant expansion tensors

[TABLE]

Stopping at this order, we here assume that the position $\mathbf{r}$ at which the displacement field is evaluated is sufficiently far away from the surface of the particle, parametrized by $\mathbf{r}^{\prime}$.

In Eqs. (9)–(13), $\mathbf{F}$ is the total force acting on the rigid particle, resulting, for example, from an external potential. This force is transmitted to the surrounding medium and leads to deformation. In a balanced stationary state, the medium presses back by a counteracting elastic restoring force. Along the same lines, $\mathbf{\underline{A}}$ is an antisymmetric second-rank tensor that is related to the total torque $\mathbf{T}$ acting on the particle from outside by

[TABLE]

with $\epsilon_{jkl}$ denoting the Levi-Civita symbol, so that

[TABLE]

The symmetric second-rank tensor $\mathbf{\underline{S}}$ is called stresslet and would be related to possible deformations of the particle surface, if the particle were deformable. The same applies to the third-rank tensor $\mathbf{\underline{\underline{M}}}$, the fourth-rank tensor $\mathbf{\underline{\underline{\underline{N}}}}$, and higher-order terms that are not considered here.

III Displacement field around a uniformly translated and rotated rigid sphere

We now consider a rigid spherical particle of radius $a$, centered at $\mathbf{r}=\mathbf{0}$, that is translated and rotated in the elastic medium by a force $\mathbf{F}$ and a torque $\mathbf{T}$, respectively. The medium shall stick to the particle surface under no-slip boundary conditions. Moreover, the displacement field $\mathbf{u}(\mathbf{r})$ induced by the sphere shall vanish at infinite distance from the particle. Thus, the boundary conditions read

[TABLE]

with $\mathbf{U}$ and $\boldsymbol{\Omega}$ the translation and rotation of the particle, respectively. Due to the linearity of the Navier–Cauchy equation, we may consider the effect of $\mathbf{F}$ and $\mathbf{T}$ separately and then, at the end, superimpose the solutions.

The integrals in this section are evaluated using Dassios (2012)

[TABLE]

with $N$ denoting the number of dyadics of $\mathbf{r}^{\prime}$ and all indices $i_{k}$ ($k=1,\dots,N$) being pairwise distinct. Only even polyadics of $\mathbf{r}^{\prime}$ survive the integration over the isotropic spherical surface. Therefore, all such integrals lead to combinations of Kronecker deltas, see also Tab. 1. Accordingly, the integrations in Eqs. (9)–(13) are performed and the results are then inserted back into the right-hand side of Eq. (8).

We start by considering only translations induced by a force $\mathbf{F}$ imposed on the particle. Due to the linearity of the problem, the force surface density $\mathbf{f}(\mathbf{r^{\prime}})$ ($\mathbf{r^{\prime}}\in\partial V$) exerted by the sphere onto the medium scales linearly with $\mathbf{F}$, as therefore does the resulting displacement field $\mathbf{u}(\mathbf{r})$. As a consequence, all contributing terms in Eq. (8) must scale linearly in $\mathbf{F}$. Obviously, the factors $\mathbf{f}(\mathbf{r^{\prime}})$ in the integrals of Eqs. (9)–(13) provide this scaling.

Since the solution of the Navier–Cauchy equation under the given boundary conditions is unique, we may continue by an appropriate ansatz for $\mathbf{f}(\mathbf{r^{\prime}})$. If the result satisfies both the Navier–Cauchy equation and the boundary conditions, we have found the correct solution.

For the translated sphere, displaced in response to a net force $\mathbf{F}$ applied to it, the ansatz corresponds to a uniform force surface density $\mathbf{f}(\mathbf{r^{\prime}})=\mathbf{F}/4\pi a^{2}$ ($\mathbf{r^{\prime}}\in\partial V$) as displayed in Fig. 1(a) Dhont (1996). Consequently, all contributions of uneven numbers of gradients in Eq. (8) vanish because of the integration of an uneven number of factors $\mathbf{r^{\prime}}$ over the surface of the sphere, see, e.g., Eqs. (10), (11), (13), and our comments above. Moreover, no even orders higher than the one given by $\mathbf{\underline{\underline{M}}}$ may contribute, because the Kronecker deltas resulting via Eq. (17) pairwise combine the gradients in Eq. (8) to Laplacian operators, which via Eq. (3) leads to vanishing terms. The integral in Eq. (12) can then readily be evaluated using Eq. (17) and Tab. 1. As a result, we obtain the familiar form Phan-Thien (1993); Puljiz and Menzel (2017)

[TABLE]

Inserting this expression into the Navier–Cauchy equation and testing the boundary conditions confirms this solution and the chosen ansatz.

We now proceed along the same lines for a torque $\mathbf{T}$ imposed on the particle. The force density for a uniform rotation leads to the ansatz $\mathbf{f}(\mathbf{r}^{\prime})\sim\mathbf{T}\times\mathbf{r}^{\prime}$ ($\mathbf{r}^{\prime}\in\partial V$), see also Fig. 1 (b). Inserting this $\mathbf{f}(\mathbf{r}^{\prime})$ into Eqs. (9)–(13), only Eq. (10) leads to a nonzero result due to antisymmetry. The expressions in Eq. (13) and higher-order terms vanish according to Eqs. (3) and/or (4) employed in Eq. (8) together with Eq. (17). The resulting displacement field reads

[TABLE]

Adding the two solutions, we obtain the overall displacement field

[TABLE]

Explicit calculation (see Tab. 2) leads to

[TABLE]

Following the boundary conditions in Eq. (16), we find

[TABLE]

IV Faxén relations

Now we consider the opposite situation, in which a given displacement field $\mathbf{u}_{0}(\mathbf{r})$ in the elastic medium translates and rotates a rigid spherical particle centered at the origin. Then, we can reformulate the second expression in Eq. (16) as

[TABLE]

which must hold true for all $\mathbf{r}\in\partial V$. The left-hand side corresponds to the translations and rotations of the points on the spherical surface shell, see Eq. (16). Due to stick boundary conditions, these must be equal to the displacement field of the surrounding elastic medium at all surface points $\mathbf{r}$, see the right-hand side of Eq. (23). However, the given displacement field $\mathbf{u}_{0}(\mathbf{r})$ will in general not be of a constant or antisymmetrically distributed form over the whole surface. Thus, it would in general imply a deformation of the spherical surface shell. This conflicts with the assumed rigidity of the particle. Therefore, the contributions in $\mathbf{u}_{0}(\mathbf{r})$ that would lead to such impossible deformations will be balanced by counteracting effects by the particle itself. The displacement fields resulting from these counteracting effects are represented by the additional integral term in Eq. (23). They are given by all terms in Eq. (8) that do not contribute to Eq. (20), i.e., by an infinite series. Multiplying Eq. (23) by polyadics of $\mathbf{r}$ and performing another integral over the surface of the sphere, the expansion tensors in Eqs. (9)–(13) appear from the integral term. Using the procedure described below, explicit expressions are calculated for the third-rank tensor $\mathbf{\underline{\underline{M}}}$ and the fourth-rank tensor $\mathbf{\underline{\underline{\underline{N}}}}$. Expressions for $\mathbf{U}$, $\boldsymbol{\Omega}$, and $\mathbf{\underline{S}}$ were already calculated in Ref. Puljiz and Menzel, 2017 and read

[TABLE]

respectively (${}^{\text{T}}$ denotes the transpose).

In order to calculate $\mathbf{\underline{\underline{M}}}$ as defined in Eq. (12), we insert the series expansion of $\mathbf{u}_{0}(\mathbf{r})$ around the origin (where the sphere is centered),

[TABLE]

into Eq. (23) and substitute the expression given in Eq. (24) into $\mathbf{U}$ in Eq. (23). Then integrating over $\int_{\partial V}\mathrm{d}S~{}r_{k}r_{l}$, we obtain

[TABLE]

Evaluation of the integral on the right-hand side yields (see Appendix A)

[TABLE]

noting that $M_{ikl}=M_{ilk}$, see Eq. (12). Reinserting Eq. (29) into Eq. (28), we obtain

[TABLE]

with

[TABLE]

Equation (30) is symmetric with respect to $k\leftrightarrow l$, but not with respect to $j\leftrightarrow k$ or $j\leftrightarrow l$. To solve Eq. (30) for $M_{jkl}$, it is useful to decompose it into a symmetric and a partially antisymmetric part. We define

[TABLE]

so that $M_{jkl}={}(M_{jkl}^{\text{S}}+M_{jkl}^{\text{A}})/2$.

Symmetrizing and antisymmetrizing Eq. (30) according to Eqs. (32) and (33), respectively, we obtain

[TABLE]

Thus, solving Eq. (30) for $M_{jkl}$, from $(M_{jkl}^{\text{S}}+M_{jkl}^{\text{A}})/2$ we obtain

[TABLE]

Setting and summing over $k=l$ in Eq. (36), we find the remaining unknown,

[TABLE]

This finally leads us to

[TABLE]

The corresponding incompressible-creeping-flow hydrodynamic result is obtained by setting $\nabla\cdot\mathbf{u}_{0}(\mathbf{r})=0$ as well as $\nu=1/2$, and by replacing the shear modulus $\mu$ by the viscosity $\eta$,

[TABLE]

where $\mathbf{u}_{0}(\mathbf{r})$ is interpreted as the flow field.

To obtain the fourth-rank tensor $\mathbf{\underline{\underline{\underline{N}}}}$, we substitute the expression from Eq. (25) into Eq. (23), use Eq. (4), insert Eq. (27), and integrate over $\int_{\partial V}\mathrm{d}S~{}r_{k}r_{l}r_{m}$, yielding (see Appendix B)

[TABLE]

with

[TABLE]

The formal solution of Eq. (40) reads

[TABLE]

and can be obtained by manipulating the whole equation in agreement with the definitions

[TABLE]

so that $N_{jklm}={}(N_{jklm}^{\text{S}}+N_{jklm}^{\text{A}})/2$, similarly to the procedure for $M_{jkl}$. $N_{nnlm}$ is obtained by setting and summing over $j=k$ in Eq. (42),

[TABLE]

Moreover, from the trace of Eq. (45),

[TABLE]

is found. Finally,

[TABLE]

with

[TABLE]

The corresponding result in incompressible low-Reynolds-number hydrodynamics is obtained for $\nabla\cdot\mathbf{u}_{0}(\mathbf{r})=0$, then setting $\nu=1/2$, and replacing $\mu$ by $\eta$. It reads

[TABLE]

with

[TABLE]

At this point, we add some brief remarks concerning our results. The $\mathbf{\underline{\underline{M}}}$-tensor in Eq. (38) for setting and summing over $k=l$ becomes $M_{jnn}=a^{2}F_{j}$, which is in agreement with the definition in Eq. (12). Therefore, $\mathbf{\underline{\underline{M}}}$ consists of two contributions: one resulting from the active displacement of the sphere by an external force $\mathbf{F}$ included by $M_{jnn}$; the other because of the rigidity of the sphere, which leads to a passive resistance to deformations that would be induced by the imposed matrix deformation $\mathbf{u}_{0}(\mathbf{r})$. It is therefore natural to split up the corresponding term in Eq. (8) accordingly, i.e.,

[TABLE]

Together with the other force term in Eq. (8), $F_{j}G_{ij}(\mathbf{r})$, and in the absence of $\mathbf{u}_{0}(\mathbf{r})$, the first term on the right-hand side of Eq. (51) brings us back to Eq. (18) and confirms the result in Sec. III from another angle (the same applies to the hydrodynamic case).

The $\mathbf{\underline{\underline{\underline{N}}}}$-tensor in Eq. (47) can be written in the more legible form

[TABLE]

Since $B_{jknn}=0$ (and also $B_{jknn}^{\text{h}}=0$ in the hydrodynamic case), the correct contraction $N_{jknn}={}a^{2}[A_{jk}+S_{jk}]$ is obtained in accordance with the definition of $\mathbf{\underline{\underline{\underline{N}}}}$ in Eq. (13). Inserting Eq. (52) into the corresponding term on the right-hand side of Eq. (8), we thus find the identity

[TABLE]

where the antisymmetric contribution from $A_{jk}$ vanishes due to Eqs. (4) and (14).

As a rule, expansion tensors of odd rank always contain the external force $\mathbf{F}$, whereas the tensors of even rank always contain the external torque $\mathbf{T}$. However, from rank four on, these contributions to the resulting displacement and rotation of the sphere vanish, if inserted into Eq. (8), because of Eqs. (2)–(4). For instance, in the rank-five tensor beyond $\mathbf{\underline{\underline{\underline{N}}}}$ in Eq. (8), the external force may appear in the form $\sim\mathbf{F}\mathbf{\hat{\underline{I}}}\,\mathbf{\hat{\underline{I}}}$ and permutations thereof. Upon contraction with $\nabla\nabla\nabla\nabla\mathbf{\underline{G}}(\mathbf{r})$ the next-higher-order contribution in Eq. (8) either contains $\nabla^{2}\nabla\cdot\mathbf{\underline{G}}(\mathbf{r})=\mathbf{0}$ or $\nabla^{4}\mathbf{\underline{G}}(\mathbf{r})=\mathbf{\underline{0}}$. Therefore, all expansion tensors of higher rank than $\mathbf{\underline{\underline{M}}}$ only contribute to the rigidity-induced resistance of the sphere to deformations.

Altogether, in differential form, the total displacement field around the sphere reads

[TABLE]

which is correct up to order $(\nabla^{n}\mathbf{u}_{0}(\mathbf{r})|_{\mathbf{r}=\mathbf{0}})r^{-m}$, with $n,m\in\mathbb{N}$ and $n+m\leq 6$. The first square bracket expresses the exact displacement field that is created by the rigid translation and rotation of the sphere due to an external force $\mathbf{F}$ and torque $\mathbf{T}$, respectively Puljiz and Menzel (2017). This here directly follows from the explicit analytical expressions for the tensors $\mathbf{\underline{\underline{M}}}$ and $\mathbf{\underline{\underline{\underline{N}}}}$ as derived above. The second square bracket contains the displacement field induced by the rigidity of the sphere in resistance to deformations from any other distortion of the embedding elastic matrix given by $\mathbf{u}_{0}(\mathbf{r})$. These conclusions directly apply to incompressible low-Reynolds-number hydrodynamics as well.

V Matrix-mediated interactions between spherical inclusions

We now consider a system consisting of $N$ identical rigid spheres of radius $a$ centered at positions $\mathbf{r}_{i}$, $i=1,2,\dots,N$, embedded under no-slip surface conditions in the linearly elastic matrix. If external forces and/or torques are applied to the particles, these forces and/or torques are transmitted to the elastic environment, there leading to deformations, which in turn results in matrix-mediated interactions between the particles Puljiz and Menzel (2017). These interactions can be calculated using the so-called method of reflections. In principle, the method proceeds by an expansion in the inverse interparticle distance. So far, the matrix-mediated interactions have been calculated up to (including) fourth order for a compressible linearly elastic matrix Puljiz and Menzel (2017). With the now-available $\mathbf{\underline{\underline{M}}}$- and $\mathbf{\underline{\underline{\underline{N}}}}$-tensors, we continue this iteration scheme up to (including) the sixth order in the inverse interparticle separation. We have described the principles of the underlying scheme in detail in Ref. Puljiz and Menzel, 2017. Again, for incompressible systems, an analogon exists for low-Reynolds-number hydrodynamics.

To proceed along these lines, the displacement field created by the force $\mathbf{F}_{i}$ and/or the torque $\mathbf{T}_{i}$ on the $i$th spherical inclusion in the absence of any other sphere is denoted as $\mathbf{u}_{i}^{(0)}(\mathbf{r})$. The resulting translation and rotation of the sphere is referred to as $\mathbf{U}_{i}^{(0)}$ and $\boldsymbol{\Omega}_{i}^{(0)}$, respectively, which according to Eq. (22) are given by

[TABLE]

Corrections to the overall translation $\mathbf{U}_{i}$ and rotation $\boldsymbol{\Omega}_{i}$ of the $i$th sphere due to matrix-mediated interactions between the spheres are taken into account by adding them to $\mathbf{U}_{i}^{(0)}$ and $\boldsymbol{\Omega}_{i}^{(0)}$, respectively, in the form $\mathbf{U}_{i}=\mathbf{U}_{i}^{(0)}+\mathbf{U}_{i}^{(1)}+\dots$ and $\boldsymbol{\Omega}_{i}=\boldsymbol{\Omega}_{i}^{(0)}+\boldsymbol{\Omega}_{i}^{(1)}+\dots$. This superposition works on the basis of the linearity of the underlying equations. We can calculate these corrections from the displacement fields

[TABLE]

and from the displacement fields

[TABLE]

with $n=1,2,3,\dots$, induced by the other spheres. By the operator $\operatorname{Tr}{{}^{2\cdot 3}}$ we imply contraction of the second and third index of the third-rank tensor $\mathbf{\underline{\underline{M}}}_{k}^{(n)}$. We apply the Faxén relations, Eqs. (24) and (25), to the displacement fields generated by the other spheres, so that the corrections to the translation $\mathbf{U}_{i}^{(0)}$ and rotation $\boldsymbol{\Omega}_{i}^{(0)}$ of a single sphere read

[TABLE]

$n\geq 1$. The displacement field $\mathbf{u}_{i}^{(0)}(\mathbf{r})$ is in this context created directly by an external force $\mathbf{F}_{i}$ and/or an external torque $\mathbf{T}_{i}$. In Eq. (57), $\mathbf{r}_{kj}=\mathbf{r}_{k}-\mathbf{r}_{j}$,

[TABLE]

see Eq. (26), and analogously for $\mathbf{\underline{\underline{M}}}_{k}^{(n)}(\mathbf{r}_{kj})$ and $\mathbf{\underline{\underline{\underline{N}}}}_{k}^{(n)}(\mathbf{r}_{kj})$, as obtained from Eqs. (38) and (47), respectively. By the summations, many-particle interactions are taken into account. In this way, expressions for $\mathbf{U}_{i}$ and $\boldsymbol{\Omega}_{i}$ are obtained as functions of the inverse interparticle distances $r_{ik}^{-1}$, with $k=1,2,\dots,N$ and $k\neq i$. We here stop at order $r_{ik}^{-6}$. Thus, three-particle interactions are involved (four-particle interactions start to contribute only at order $r_{ik}^{-7}$). That is, one particle generates a displacement field by deforming the surrounding matrix, initiated by a force or torque imposed on this particle. A second particle resists the deformation implied by the resulting deformation field, which is expressed by the higher-order moments of its surface force density, see Eqs. (11)–(13). Accordingly, a higher-order displacement field is generated, see Eq. (57), to which a third particle (or again the first particle) is exposed. For rigidly translated and/or rotated spherical inclusions, we thus stop at the level of $\mathbf{U}_{i}^{(2)}$ and $\boldsymbol{\Omega}_{i}^{(2)}$. In the end, the effect of the coupled matrix-mediated interactions between the spheres can be written in the form of mathematical displaceability and rotateability matrices. Thus, the total translation and rotation of the $i$th particle to the given order are obtained as

[TABLE]

respectively, for $i=1,2,\dots,N$.

So far, the results up to (including) order $r_{ij}^{-4}$ had been calculated in Ref. Puljiz and Menzel, 2017. Our scope in the following subsections is to determine the contributions of the orders $r_{ij}^{-5}$ and $r_{ij}^{-6}$. The complete expressions for the displaceability and rotateability matrices up to (including) order $r_{ij}^{-6}$ are for completeness listed as well.

V.1 Displaceability matrices

First, we consider the situation, in which the $i$th particle is rigidly displaced by an external force $\mathbf{F}_{i}$. This leads to the particle displacement $\mathbf{U}_{i}^{(0)}$, see Eq. (55), and induces the displacement field $\mathbf{u}_{i}^{(0)}(\mathbf{r})$ in the elastic medium, see Eq. (56). The remaining contributions, to complete the order $r_{ik}^{-6}$ beyond the order $r_{ik}^{-4}$, see Ref. Puljiz and Menzel, 2017, follow from the fields

[TABLE]

see Eq. (57).

Here, the tensors are calculated via Eqs. (26), (38), and (47) by inserting into Eq. (60) and analogously into the corresponding expressions for $\mathbf{\underline{\underline{M}}}_{k}^{(1)}(\mathbf{r}_{kj})$ and $\mathbf{\underline{\underline{\underline{N}}}}_{k}^{(1)}(\mathbf{r}_{kj})$ the field $\mathbf{u}_{j}^{(0)}(\mathbf{r})$ obtained via Eq. (56). For $\mathbf{\underline{S}}_{k}^{(1)}(\mathbf{r}_{kj})$, we find

[TABLE]

Therefore, the sixth-order contribution to $\mathbf{U}_{i}^{(2)}$ arising from the stresslet generated on particle $k$ due to its resistance to deformation in a displacement field generated by a force $\mathbf{F}_{j}$ acting on particle $j$ is obtained from

[TABLE]

as

[TABLE]

with $M_{0}^{\text{t}}=(5-6\nu)/24\pi(1-\nu)\mu a$. (There is no fifth-order contribution involving $\mathbf{\underline{S}}_{k}^{(1)}$.)

Next, the $\mathbf{\underline{\underline{M}}}_{k}^{(1)}(\mathbf{r}_{kj})$-term is addressed. We find from

[TABLE]

in Eq. (38) to third order the contribution

[TABLE]

Here, for better readability, we have shifted the particle indices to the superscript, which are not summed over, and the coordinate indices to the subscript. Thus, we obtain to sixth order from

[TABLE]

the expression

[TABLE]

As for $\mathbf{\underline{\underline{\underline{N}}}}_{k}^{(1)}(\mathbf{r}_{kj})$, we find via Eq. (53) and the stresslet in Eq. (64) that from $-N_{jklm}^{k(1)}(\mathbf{r}^{kj})/6$ only

[TABLE]

adds to our considered order. This leads to a sixth-order contribution from

[TABLE]

of the form

[TABLE]

The additional rotation $\boldsymbol{\Omega}_{i}^{(2)}$ up to (including) sixth order is obtained solely from the stresslet term and, via Eqs. (59), (63), and (64), reads

[TABLE]

with $M_{0}^{\text{r}}=1/8\pi\mu$ (there is only this additional fifth-order contribution up to the considered level).

In total, together with our previously derived expressions to fourth order Puljiz and Menzel (2017), the sixth-order displaceability and rotateability matrices for translation–translation and rotation–translation couplings, respectively, read

[TABLE]

and

[TABLE]

V.2 Rotateability matrices

Likewise, an external torque $\mathbf{T}_{i}$ acting on the $i$th particle to zeroth order leads to a displacement field $\mathbf{u}_{i}^{(0)}(\mathbf{r})$ as given by Eq. (56) and the particle rotation $\boldsymbol{\Omega}_{i}^{(0)}$ as specified by Eq. (55). Overall, up to (including) order $r_{ij}^{-6}$, only third-order parts of the stresslets $\mathbf{\underline{S}}_{k}^{(1)}(\mathbf{r}_{kj})$ contribute to $\mathbf{U}_{i}^{(2)}$ and $\boldsymbol{\Omega}_{i}^{(2)}$, which here according to Eqs. (56) and (60) are of the form

[TABLE]

They yield the additional contributions to the translation and rotation, respectively, of

[TABLE]

according to Eqs. (58) and (63), and of

[TABLE]

according to Eqs. (59) and (63) (there is no sixth-order contribution to the translation resulting from an imposed torque). We stress that the sixth-order contribution to the rotation is nonzero also for incompressible systems, i.e., for $\nu\rightarrow 1/2$.

In summary, the translation–rotation and rotation–rotation coupling matrices, together with our previously-derived expressions to fourth order Puljiz and Menzel (2017), are given by

[TABLE]

and

[TABLE]

respectively.

For all derived displaceability and rotateability matrices, the expressions of corresponding low-Reynolds-number hydrodynamic mobility matrices for incompressible fluid flows are obtained by setting $\nu=1/2$ in the limit of an incompressible environment and replacing $\mu$ by the hydrodynamic viscosity $\eta$. Previously, for such incompressible fluid surroundings, the expressions for the hydrodynamic mobility matrices had been derived in a different way Mazur and Van Saarloos (1982). These expressions are recovered from ours for $\nu=1/2$.

V.3 Examples: Magnetic particles

We now present examples for two- and three-particle systems. Specifically, the difference between particle interactions through an elastic compressible or an incompressible environment is illustrated, together with the effect of the higher-order terms. For simplicity, we assume magnetic dipole–dipole interactions between $N$ identical particles embedded in an elastic medium, reminiscent of the situation in ferrogels and magnetorheological elastomers Jarkova et al. (2003); Attaran et al. (2017); Weeber et al. (2018); Metsch et al. (2018). Corresponding magnetic dipole forces are given by Jackson (1999)

[TABLE]

$i=1,\dots,N$. Here, $\mu_{0}$ is the magnetic vacuum permeability and $\mathbf{m}$ represents the magnetic dipole moment that we assume to be identical for all particles ($m=|\mathbf{m}|$, $\mathbf{\hat{m}}=\mathbf{m}/m$). This can be achieved, for instance, by a strong homogeneous saturating external magnetic field $\mathbf{B}$ magnetizing the spherical particles Puljiz et al. (2016), with the magnetic moments oriented along $\mathbf{B}$. If particle displacements occur, Eq. (90) should be evaluated with the modified particle distances considered Puljiz et al. (2016). An iterative numerical loop is employed for this purpose to calculate the final displacements Puljiz et al. (2016), which in the regime of linear elasticity converges quickly. In the following, we set $\mu_{0}m^{2}/\mu a^{6}=5\times 10^{2}$. Moreover, all initial interparticle center-to-center distances are set to $3a$.

Assuming uniaxial magnetic anisotropy of the particles, magnetic torques are calculated from the idealized Stoner–Wohlfarth model Stoner and Wohlfarth (1948); Puljiz and Menzel (2017); Roeder et al. (2015),

[TABLE]

In this expression, $\mathbf{\hat{B}}=\mathbf{B}/|\mathbf{B}|$ and $\mathbf{\hat{n}}_{i}$ represents the nonpolar axis of magnetic anisotropy of the $i$th particle, rigidly anchored to the particle frame. We set the anisotropy parameter $K/\mu=1/3$.

In Figs. 2–5, the $xy$-plane contains the particle centers. The external magnetic field is rotated counterclockwise within the plane (around the $\mathbf{\hat{z}}$-axis). Particle displacements and rotations due to the induced forces and torques are shown for $\nu=0.5$, $0.3$, and [math] (solid, dashed, and dotted lines, respectively). It can be seen in all examples that the translational displacements get significantly larger for decreasing $\nu$, i.e., increasing compressibility of the elastic environment. In Figs. 2 and 3, both forces and torques are induced. The differences in the rotations depicted for the different values of $\nu$, both in Figs. 2 (c) and in 3 (c), are very small and here hardly visible. This is because the Poisson ratio does not enter the $\mathbf{\underline{M}}_{ij}^{\text{rt}}$- and $\mathbf{\underline{M}}_{ij}^{\text{rr}}$-matrices below order $r_{ij}^{-5}$ and $r_{ij}^{-6}$, see Eqs. (80) and (89), respectively.

In Fig. 4, the particles do not feature an axis of magnetic anisotropy. Therefore, only forces are induced by the magnetic field. Thus the rotations in Fig. 4 (c) stem solely from the $\mathbf{\underline{M}}_{ij}^{\text{rt}}$-couplings, see Eqs. (78)–(80). The pronounced $\nu$-dependence partly originates from the magnetic forces, which were calculated as described above. Since for the lower considered values of $\nu$ the displacements are larger, also the changes of the magnetic forces can become larger in magnitude according to the increased changes in the mutual particle distances. This, in turn, affects the magnitudes of the induced rotations via the $\mathbf{\underline{M}}_{ij}^{\text{rt}}$-matrices. Additionally, the higher-order contributions from the $\mathbf{\underline{M}}_{ij}^{\text{rt(5)}}$-matrices, see Eq. (80), which include three-particle interactions and which depend on $\nu$, are now more exposed. The particular shapes of the rotation curves for $\nu=0.5$ are in fact to a larger extent affected by the $\mathbf{\underline{M}}_{ij}^{\text{rt(5)}}$-matrices. For smaller $\nu$, the lower-order $\mathbf{\underline{M}}_{i\neq j}^{\text{rt}}$-matrices, see Eq. (79), here lead to more regular rotation curves. For comparison, the same configuration as in Fig. 4 is repeated in Fig. 5, there with permanent axes of magnetic anisotropy of the particles included. Thus, magnetic torques are additionally induced. The rotations in Fig. 5 (c) are then again dominated by the magnetic torques via the $\mathbf{\underline{M}}_{i=j}^{\text{rr}}$-matrices, see Eq. (87).

Summarizing, we observe that compressibility can play a significant role for the coupling of the translations and rotations of rigid particles mediated by an embedding elastic environment. This has at least two sources. First, a compressible environment can allow for larger displacements of already a single particle on its own. The resulting larger displacements can lead to stronger interactions between the particles. Second, the coupling via the environment itself is affected directly, see the $\nu$-dependence of the expressions for the coupling matrices in Secs. V.1 and V.2. We note that the Poisson ratio enters the rotation–translation and rotation–rotation couplings only at fifth and sixth order, respectively. Therefore, in many cases the corresponding effects associated with compressibility of the elastic environment will be masked if torques are imposed directly on each embedded particle.

V.4 Higher-order effects

Already during our derivation of the additional fifth- ($r_{ij}^{-5}$) and sixth- ($r_{ij}^{-6}$) order contributions to the coupling matrices in Secs. V.1 and V.2, we have seen that qualitatively new features arise from these higher orders. Particularly, the rotation–translation, translation–rotation, and rotation–rotation couplings now become explicitly dependent on the compressibility of the embedding elastic environment. These effects are absent up to (including) the fourth order ($r_{ij}^{-4}$) considered before Puljiz and Menzel (2017).

To also obtain basic quantitative numbers on the typical role of the additional contributions of orders $r_{ij}^{-5}$ and $r_{ij}^{-6}$ derived in Secs. V.1 and V.2, we consider two minimal example situations. First, one rigid sphere is subjected to a constant force $\mathbf{F}_{1}$ that pushes it towards another rigid, nearby, “passive” sphere, see Fig. 6. Second, one rigid sphere is rotated by a constant torque $\mathbf{T}_{1}$ in the presence of another rigid, nearby, “passive” sphere, see Fig. 7. In practice, such a situation could be realized, for instance, by a configuration of two nearby rigid spherical particles of equal size. However, the first one of them is magnetic, the second one is not. Exposing the magnetic sphere to a magnetic field gradient, it will be subject to a net force Jackson (1999), while already a homogeneous field misaligned with an axis of internal magnetic anisotropy will impose a torque, see Eq. (91). Here, we set $|\mathbf{F}_{1}|=\mu a^{2}$ and $|\mathbf{T}_{1}|=\mu a^{3}$. We reduce the distance $r_{12}$ between the centers of the two spheres and determine the deviation of the solution up to (including) order $r_{12}^{-6}$ derived in Secs. V.1 and V.2 from the previously considered Puljiz and Menzel (2017) solution up to (including) order $r_{12}^{-4}$.

In both situations, the presence of the second, “passive” sphere affects the translation or rotation of the first one because its rigidity opposes to deformations in the environment generated by the first sphere. Accordingly, part of the displacement field induced by the first sphere is “reflected back” towards the latter. Concerning the results depicted in Fig. 6 for the case of an imposed force, this effect is included by the $\mathbf{\underline{M}}_{ij}^{\text{tt}(4,6)}$-matrices, see Eq. (77). There are already fourth-order contributions ($r_{12}^{-4}$) of this kind, with further terms arising to sixth order ($r_{12}^{-6}$). Obviously, when decreasing the distances $r_{12}$ towards $2a$, the sixth-order contributions become more pronounced. Moreover, for the considered values of the Poisson ratio $\nu=0$, $0.3$, and $0.5$, the strongest relative deviation between the sixth- and fourth-order solutions is found for $\nu=0.5$, see Fig. 6 (c).

Turning to the situation of an imposed torque on the first particle, see Fig. 7, no reflections of the kind described above occur up to (including) the fourth order, to which the rotation $\boldsymbol{\Omega}_{1}$ of sphere 1 is independent of $r_{12}$, see Eq. (87). [Since there is no magnetic torque acting on sphere 2, the $\mathbf{\underline{M}}_{i\neq j}^{\text{rr}}$-matrices given by Eq. (88) do not contribute to $\boldsymbol{\Omega}_{1}$.] Nevertheless, the presence of the rigid second sphere is felt by the first sphere to sixth order ($r_{12}^{-6}$), see Eq. (89). Moreover, a slightly different behavior is found for the different considered values of $\nu$, see Fig. 7 (b) and (c). The observed attenuation of the $z$-component of $\boldsymbol{\Omega}_{1}$ is most pronounced for $\nu=0.5$, and it grows in magnitude when decreasing the interparticle center-to-center distance $r_{12}$ towards $2a$.

Altogether, we observe in our examples that the higher-order contributions derived above become increasingly relevant when the Poisson ratio is growing from $\nu=0$ towards the limit of incompressibility of $\nu=0.5$. Moreover, the additional contributions notably grow when decreasing the center-to-center interparticle distance below approximately $3a$. For larger interparticle distances, fourth-order solutions for the particle displacements and rotations appear to be already of a reasonable quantitative precision Puljiz et al. (2016, 2018). At some point, when the particles come too close to each other, depending on the particular material and system investigated, induced heterogeneities around the particles Huang et al. (2016) or nonlinear elastic effects Puljiz et al. (2018) may become important.

VI Conclusions

In summary, we have presented a modified route of calculating both the classical solution of the displacement field and the incompressible Stokes flow around a translated and rotated rigid sphere of no-slip surface conditions in a surrounding possibly compressible linearly elastic solid or incompressible low-Reynolds-number fluid, respectively. (The latter result formally follows from the former by considering the limit of an incompressible surrounding medium.) For this purpose, the displacement field in an elastic environment has been expanded around the spherical particle. Then an explicit ansatz has been inserted for the surface force density as also presented, e.g., in Ref. Dhont, 1996. Afterwards, we showed by explicit calculation that the ansatz renders the expansion finite.

Next, the two Faxén laws beyond the stresslet describing the reaction of the sphere to an imposed deformation or flow field in its environment were derived, both for linear elasticity theory (within a possibly compressible elastic solid) and for low-Reynolds-number hydrodynamics (within an incompressible fluid), respectively.

Based on these results, we calculated the displaceability and rotateability matrices up to the sixth order in inverse particle separation distance, characterizing the particle interactions mediated by the surrounding compressible elastic environment. Differences when compared to the incompressible situation were demonstrated via basic example situations, together with a brief illustration of the effects of the additional higher-order contributions derived in the present work.

Our results will be important in the future, for instance, to quantify to higher-order precision the properties and behavior of complex elastic composite materials and their response to external stimuli Filipcsei et al. (2007); Menzel (2015); Odenbach (2016); Puljiz et al. (2018). Since, frequently, in soft complex elastic composite or multi-component materials a pronounced nonlinear stress-strain behavior is strongly dominated by interactions resulting from the additional components Menzel et al. (2009); Cremer et al. (2016) and not necessarily by the nonlinear elasticity of the purely elastic component, our results derived in the framework of linear elasticity theory may to some extent even serve to characterize more in detail the nonlinear stress-strain behavior of such systems.

Acknowledgements.

The authors thank the Deutsche Forschungsgemeinschaft (DFG) for support of this work through the priority program SPP 1681, Grant No. ME 3751/3.

Appendix A

The integral on the right-hand side of Eq. (28),

[TABLE]

can be solved by Fourier forth and back transform. The Fourier transform of the Green’s function reads Puljiz and Menzel (2017)

[TABLE]

where we introduced the wave vector $\mathbf{k}$, $k=|\mathbf{k}|$, and $\mathbf{\hat{k}}=\mathbf{k}/k$. First, only the $\mathrm{d}S$-integral is considered,

[TABLE]

Using the identity

[TABLE]

where

[TABLE]

and

[TABLE]

we obtain

[TABLE]

After solving the $\mathrm{d}k$-integral, only the part even in $\mathbf{\hat{k}}$ can survive the subsequent $\mathrm{d}S(\mathbf{\hat{k}})$-integral. Thus, the limits of integration may here be rewritten as $\int_{0}^{\infty}\mathrm{d}k=\frac{1}{2}\int_{\mathbb{R}}\mathrm{d}k$. With

[TABLE]

for $-1<\mathbf{\hat{k}}\cdot\mathbf{r}^{\prime}/a<1$, which can be obtained by Mathematica Wolfram Research Inc. (2016) as well as by partial integration and further proceeding as outlined in the appendices of Refs. Dhont, 1996 and Puljiz and Menzel, 2017, respectively, Eq. (A.7) becomes

[TABLE]

Evaluation with the help of Tabs. 1 and 3 directly leads to Eq. (29).

Appendix B

The same procedure as in Appendix A is applied to evaluate the integral

[TABLE]

now using the identity

[TABLE]

This leads to

[TABLE]

As in Appendix A, we may equivalently use $\int_{0}^{\infty}\mathrm{d}k=\frac{1}{2}\int_{\mathbb{R}}\mathrm{d}k$ here. Integrating by parts, we find

[TABLE]

The remaining integrals are the same as in Eqs. (A.9) and (A.8), respectively. Altogether, we obtain

[TABLE]

finally resulting in Eq. (40) after evaluation of the integrals using Tabs. 1 and 3.