Forces on rigid inclusions in elastic media and resulting matrix-mediated interactions

TL;DR

This paper derives analytical expressions for the elastic interactions between spherical particles embedded in elastic media, supported by experiments with magnetic particles, aiming to extend hydrodynamic interaction methods to elasticity.

Contribution

It introduces a novel analytical approach to calculate elastic interactions between particles in elastic media, inspired by hydrodynamic interaction techniques.

Findings

Analytical formulas for particle interactions in elastic media

Experimental validation with magnetic particles

Potential applications in elastic composite materials

Abstract

To describe many-particle systems suspended in incompressible low-Reynolds-number fluids, effective hydrodynamic interactions can be introduced. Here, we consider particles embedded in elastic media. The effective elastic interactions between spherical particles are calculated analytically, inspired by the approach in the fluid case. Our experiments on interacting magnetic particles confirm the theory. In view of the huge success of the method in hydrodynamics, we similarly expect many future applications in the elastic case, e.g. for elastic composite materials.

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}).$$

$$\mathbf{\hskip 0.56917pt\underline{G}}(\mathbf{r})={}\frac{1}{{\color[rgb]{0,0,0}8\pi}\mu}\left[\frac{{\color[rgb]{0,0,0}1}}{r}\mathbf{\underline{\hat{I}}}+\frac{\mathbf{r}\mathbf{r}}{r^{3}}\right]$$

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

$$\mathbf{U}_{i}^{(0)}\,=\,\mathbf{u}_{i}^{(0)}({\color[rgb]{0,0,0}|\mathbf{r}-\mathbf{r}_{i}|}=a)\,=\,\frac{{\color[rgb]{0,0,0}1}}{{\color[rgb]{0,0,0}6\pi}\mu a}\mathbf{F}_{i}$$

$${\mathbf{U}}_j^{(1)}+{\mathbf{\Omega }}_j^{(1)}\times (\mathbf{r}-{\mathbf{r}}_j)={\mathbf{u}}_i^{(0)}(\mathbf{r})+{\int }_{\partial V_j}\underline{\mathbf{G}}(\mathbf{r}-{\mathbf{r}}^{\prime })\cdot \mathbf{f}({\mathbf{r}}^{\prime })\mathrm{d}{S}^{\prime }$$

$${\mathbf{U}}_j^{(1)}$$

$${\mathbf{\Omega }}_j^{(1)}$$

$$\mathbf{\underline{S}}_{j}^{(1)}\!=\frac{{\color[rgb]{0,0,0}10\pi}\mu a^{3}}{-{\color[rgb]{0,0,0}3}}\left(1+\frac{a^{2}}{10}\nabla^{2}\right)\!\left[\nabla\mathbf{u}_{i}^{(0)}(\mathbf{r})+\big{(}\nabla\mathbf{u}_{i}^{(0)}(\mathbf{r})\big{)}^{\!T}\right]\!\Big{|}_{\mathbf{r}_{j}}\!\!.$$

$${\mathbf{u}}_j^{(1)}(\mathbf{r})=-\left({\underline{\mathbf{S}}}_j^{(1)}\cdot \nabla \right)\cdot \underline{\mathbf{G}}(\mathbf{r}-{\mathbf{r}}_j).$$

$$\left(\begin{array}{c}{\mathbf{U}}_1\\ ⋮\\ {\mathbf{U}}_N\end{array}\right)=\left(\begin{array}{ccc}{\underline{\mathbf{M}}}_{11}& \dots & {\underline{\mathbf{M}}}_{1N}\\ ⋮& ⋮& ⋮\\ {\underline{\mathbf{M}}}_{N1}& \dots & {\underline{\mathbf{M}}}_{NN}\end{array}\right)\cdot \left(\begin{array}{c}{\mathbf{F}}_1\\ ⋮\\ {\mathbf{F}}_N\end{array}\right).$$

$${\underline{\mathbf{M}}}_{i=j}$$

$${\underline{\mathbf{M}}}_{i\ne j}$$

$$\mathbf{\underline{M}}_{i\not=j}^{(3)}\!=\!M_{0}\frac{{\color[rgb]{0,0,0}15}}{{\color[rgb]{0,0,0}8}}\!\sum\limits_{\begin{subarray}{c}k=1\\ k\not=i,j\end{subarray}}^{N}\!\bigg{(}\frac{a}{r_{ik}}\bigg{)}^{\!\!\!2}\bigg{(}\frac{a}{r_{jk}}\bigg{)}^{\!\!\!2}\Big{[}{\color[rgb]{0,0,0}1-3}(\mathbf{\hat{r}}_{ik}\cdot\mathbf{\hat{r}}_{jk})^{2}\Big{]}\mathbf{\hat{r}}_{ik}\mathbf{\hat{r}}_{jk}.$$

$$\mathbf{F}_{i}={}-\frac{3\mu_{0}m^{2}}{4\pi}\sum\limits_{\begin{subarray}{c}j=1\\ j\not=i\end{subarray}}^{N}\frac{5\mathbf{\hat{r}}_{ij}(\mathbf{\hat{m}}\cdot\mathbf{\hat{r}}_{ij})^{2}-\mathbf{\hat{r}}_{ij}-2\mathbf{\hat{m}}(\mathbf{\hat{m}}\cdot\mathbf{\hat{r}}_{ij})}{{\color[rgb]{0,0,0}r_{ij}^{\>4}}},$$

$${\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}).$$

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

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

$$\mathbf{U}_{i}^{(0)}\,=\,\mathbf{u}_{i}^{(0)}({\color[rgb]{0,0,0}|\mathbf{r}-\mathbf{r}_{i}|}=a)\,=\,\frac{5-6\nu}{24\pi(1-\nu)\mu a}\mathbf{F}_{i}$$

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

$${\mathbf{U}}_j^{(1)}$$

$${\mathbf{\Omega }}_j^{(1)}$$

$${\underline{\mathbf{S}}}_j^{(1)}$$

$${\mathbf{u}}_j^{(1)}(\mathbf{r})=-\left({\underline{\mathbf{S}}}_j^{(1)}\cdot \nabla \right)\cdot \underline{\mathbf{G}}(\mathbf{r}-{\mathbf{r}}_j).$$

$$\left(\begin{array}{c}{\mathbf{U}}_1\\ ⋮\\ {\mathbf{U}}_N\end{array}\right)=\left(\begin{array}{ccc}{\underline{\mathbf{M}}}_{11}& \dots & {\underline{\mathbf{M}}}_{1N}\\ ⋮& ⋮& ⋮\\ {\underline{\mathbf{M}}}_{N1}& \dots & {\underline{\mathbf{M}}}_{NN}\end{array}\right)\cdot \left(\begin{array}{c}{\mathbf{F}}_1\\ ⋮\\ {\mathbf{F}}_N\end{array}\right).$$

$${\underline{\mathbf{M}}}_{i=j}$$

$${\underline{\mathbf{M}}}_{i\ne j}$$

$${\underline{\mathbf{M}}}_{i\ne j}^{(3)}$$

Full text

Hydrodynamics determines our daily life. Examples are given by the flow of air into our lungs Zhang and Kleinstreuer (2002), drinking of beverages and digestive processes Meng et al. (2005); *ferrua2010modeling, technical applications such as microfluidic devices Squires and Quake (2005), or shape optimization of planes, vehicles, ships, and propellers Wald (2006); *campana2006shape; *muller2014aerodynamic. All these processes are described by the Navier-Stokes equations Navier (1822); *stokes1845theories or variants thereof. This set of equations typically poses significant challenges during solution due to a convective nonlinearity reflecting inertial effects. Basically, turbulence is driven by the inertial term. It often renders analytical solutions impossible.

The situation changes for small dimensions and velocities or high viscosity. Then, the relative strength of inertial effects, measured by the Reynolds number, is low. The nonlinearity can be neglected. A Green’s function in terms of the so-called Oseen matrix is then available, which formally solves the problem analytically Kim and Karrila (1991); Dhont (1996). In this way, semi-dilute colloidal suspensions, i.e. the dispersion of nano- to micrometer-sized particles in a fluid Dhont (1996); Felderhof (1977); *ermak1978brownian; *durlofsky1987dynamic; *zahn1997hydrodynamic; *meiners1999direct; *dhont2004thermodiffusion; *rex2008influence, or microswimmer suspensions Pooley et al. (2007); *baskaran2009statistical; *menzel2016dynamical; *lauga2009hydrodynamics; *drescher2010direct; *drescher2011fluid; *paxton2004catalytic are described effectively. The explicit role of the fluid is eliminated and replaced by effective hydrodynamic interactions between the suspended particles Kim and Karrila (1991); Dhont (1996).

Despite the success of this theoretical approach for colloidal suspensions, hardly any investigations consider a surrounding elastic solid instead of a suspending fluid. This is surprising, since, as we show below, the formalism can be adapted straightforwardly to linearly elastic matrices and is confirmed by our experiments. Our approach will, for instance, facilitate describing the response of elastic composite materials to external stimuli. Such materials consist of more or less rigid inclusions embedded in an elastic matrix. They are of growing technological interest and may serve, e.g., as soft actuators or sound attenuation devices An and Shaw (2003); *fuhrer2009crosslinking; *bose2012soft; *cheng2006observation; *still2011collective; *baumgartl2007tailoring; *baumgartl2007erratum.

In previous theoretical studies, the physics of one single rigid or deformable inclusion was addressed Eshelby (1957); *eshelby1959elastic; *walpole1991rotated; *walpole1991translated; *walpole2005green; Phan-Thien (1993), also under acoustic irradiation Oestreicher (1951); *norris2006impedance; *norris2008faxen. For more than a single inclusion, mainly the so-called load problem was analyzed theoretically for a pair of rigid inclusions: one prescribes displacements of two rigid inclusions in an elastic matrix, and then determines the forces necessary to achieve these given displacements Phan-Thien and Kim (1994); *kim1995faxen.

Here, we take the converse point of view, based on the cause-and-effect chain in our experiments: external forces are imposed onto the inclusions, or mutual forces between the inclusions are induced, for example to actuate the material or to tune its properties. In response to the forces, the inclusions are displaced. Since they cannot penetrate through the surrounding elastic matrix, they transmit the forces to the matrix and distort it. Such distortions lead to mutual long-ranged interactions between the inclusions, in analogy to hydrodynamic interactions in colloidal suspensions Kim and Karrila (1991); Dhont (1996); Tanaka and Araki (2000).

We present a basic derivation of analytical expressions for these interactions from the underlying elasticity equations. Then, we verify the theory by experiments on rigid paramagnetic particles embedded in soft elastic matrices. Mutual particle interactions are induced by an external magnetic field. As we demonstrate, theory and experiment are in good agreement, and also allow for microrheological measurements Ziemann et al. (1994); *bausch1999measurement; *waigh2005microrheology; *wilhelm2008out.

For simplicity, we assume a homogeneous, isotropic, infinitely extended elastic matrix, and low-amplitude deformations. Applying a bulk force density $\mathbf{f}_{b}(\mathbf{r})$ to the matrix, its equilibrated state satisfies the linear elastostatic Navier-Cauchy equations Cauchy (1828),

[TABLE]

This is the elastic analogue to the linearized Stokes equation in low-Reynolds-number hydrodynamics Dhont (1996). Instead of velocities, $\mathbf{u}(\mathbf{r})$ here denotes the displacement field, describing the reversible relocations of the volume elements from their initial positions during deformations. $\mu$ is the shear modulus of the matrix and $\nu$ its Poisson ratio, connected to its compressibility Landau and Lifshitz (1986). We consider an incompressible matrix, i.e. $\nabla\cdot\mathbf{u}(\mathbf{r})=0$ along with $\nu=0.5$. Yet, in contrast to the hydrodynamic case, also compressible elastic systems are readily addressed, and we present the corresponding expressions in the supplemental material sup .

Importantly, for a point force density $\mathbf{f}_{b}(\mathbf{r})=\mathbf{F}\delta(\mathbf{r})$ acting on the matrix, the resulting deformation field can be calculated analytically from Eq. (1) via Fourier transform as $\mathbf{u}(\mathbf{r})=\mathbf{\hskip 0.56917pt\underline{G}}(\mathbf{r})\cdot\mathbf{F}$. Here,

[TABLE]

is the corresponding Green’s function Landau and Lifshitz (1986), $\mathbf{\underline{\hat{I}}}$ the identity matrix, $r$=$|\mathbf{r}|$, and the underscore marks second-rank tensors and matrices. Still, it is practically impossible to explicitly solve Eq. (1) analytically in the presence of several rigid embedded particles of finite size. An iterative procedure resolves this problem, see Fig. 1.

We consider $N$ rigid spherical particles of radius $a$, with no-slip boundary conditions on their surfaces. First we only address the $i$th particle at position $\mathbf{r}_{i}$, subject to an external force $\mathbf{F}_{i}$. The embedded particle transmits this force to the surrounding matrix and induces a displacement field

[TABLE]

This field is the elastic analogue of hydrodynamic Stokes flow Kim and Karrila (1991); Dhont (1996), for elastic media. Inserting Eq. (2) reproduces a corresponding expression in Ref. Phan-Thien, 1993. Eq. (3) is confirmed as it satisfies Eq. (1), shows the correct limit \mathbf{u}_{i}^{(0)}(\mathbf{r})=\mathbf{\hskip 0.56917pt\underline{G}}(\mathbf{r}{\color[rgb]{0,0,0}-\mathbf{r}_{i}})\cdot\mathbf{F}_{i} for $|\mathbf{r}-\mathbf{r}_{i}|>a$ when $a\rightarrow 0$, and for {\color[rgb]{0,0,0}|\mathbf{r}-\mathbf{r}_{i}|}=a is constant on the particle surface. Thus, Eq. (3) for {\color[rgb]{0,0,0}|\mathbf{r}-\mathbf{r}_{i}|}=a reveals the rigid displacement

[TABLE]

of the $i$th particle in response to $\mathbf{F}_{i}$ in accord with the no-slip conditions at {\color[rgb]{0,0,0}|\mathbf{r}-\mathbf{r}_{i}|}=a.

To find the effective elastic interactions between particles $i$ and $j$ ($j\neq i$), we take the induced displacement field $\mathbf{u}_{i}^{(0)}(\mathbf{r})$ as given. We need to determine how particle $j$ reacts to the imposed field $\mathbf{u}_{i}^{(0)}(\mathbf{r})$. In general, particle $j$ can be rigidly translated by a displacement vector $\mathbf{U}_{j}^{(1)}$ and rigidly rotated by a rotation vector $\mathbf{\Omega}_{j}^{(1)}$. Taking into account the no-slip conditions on the surface $\partial V_{j}$ of the $j$th particle, the equality

[TABLE]

must hold for all $\mathbf{r}\in\partial V_{j}$. That is, the rigid displacement of each point on the surface shell of particle $j$ (l.h.s.) must equal the displacement field in the matrix at the same point (r.h.s.). The latter is given by the imposed displacement field, here $\mathbf{u}_{i}^{(0)}(\mathbf{r})$, plus the deformation that the particle surface itself induces in the matrix, i.e. the integral term. Also an externally imposed global displacement field could be included (on the r.h.s.). $\mathbf{f}(\mathbf{r}^{\prime})$ describes the surface force density exerted by the surface of particle $j$ onto the matrix.

Such an embedded particle will translate and rotate as dictated by the surrounding matrix. We obtain the expression for $\mathbf{U}_{j}^{(1)}$ by integrating Eq. (5) over $\partial V_{j}$. Similarly, for $\mathbf{\Omega}_{j}^{(1)}$, Eq. (5) is multiplied dyadically by $\mathbf{r}-\mathbf{r}_{j}$, and after integration over $\partial V_{j}$ the antisymmetric part is extracted. To perform the calculation, $\mathbf{u}_{i}^{(0)}(\mathbf{r})$ is Taylor expanded around $\mathbf{r}_{j}$. Moreover, we use that Eq. (1) for $\mathbf{r}\notin\partial V_{i}$ leads to $\nabla^{4}\mathbf{u}_{i}^{(0)}(\mathbf{r})=\mathbf{0}$ and $\nabla\times\nabla^{2}\mathbf{u}_{i}^{(0)}(\mathbf{r})=\mathbf{0}$. The last term in Eq. (5) vanishes at this stage as no total net external force or torque is applied to particle $j$ at the present step of iteration. In the end, we recover the elastic analogues of the hydrodynamic Batchelor and Green (1972); Kim and Karrila (1991); Dhont (1996) Faxén laws

[TABLE]

This is how particle $j$ is translated and rotated in the field $\mathbf{u}_{i}^{(0)}(\mathbf{r})$ induced by particle $i$. Yet, elastic retroaction occurs between the particles, as described in the following.

The force densities $\mathbf{f}(\mathbf{r}^{\prime})$ in Eq. (5) that the particles exert on their environment in general will not vanish identically. Since the particles are rigid, they resist any deformation that $\mathbf{u}_{i}^{(0)}$ would imply. Thus, they exert counteracting stresses onto the deformed matrix. The stresslet exerted by particle $j$ onto the matrix can be denoted as \mathbf{\underline{S}}_{j}=\int_{\partial V_{j}}\mathrm{d}S^{\prime}\{[\mathbf{f}(\mathbf{r}^{\prime})\mathbf{r}^{\prime}+(\mathbf{f}(\mathbf{r}^{\prime})\mathbf{r}^{\prime})^{T}]/2-{\color[rgb]{0,0,0}\mathbf{\underline{\hat{I}}}\,[\mathbf{f}(\mathbf{r}^{\prime})\cdot\mathbf{r}^{\prime}]/3}\}, where $[\bullet]^{T}$ marks the transpose. In our case, we can directly calculate from Eq. (5) the stresslet $\mathbf{\underline{S}}_{j}^{(1)}$ that particle $j$ exerts onto the matrix when it resists to the deformation described by $\mathbf{u}_{i}^{(0)}(\mathbf{r})$. To find the expression for $\mathbf{\underline{{S}}}_{j}^{(1)}$, one proceeds in the same way as described above for $\mathbf{\Omega}_{j}^{(1)}$ but eventually extracts the symmetric part. The latter contains the definition of $\mathbf{\underline{S}}_{j}^{(1)}$. We obtain

[TABLE]

This stresslet leads to additional distortions of the matrix, see Fig. 1, described by a displacement field $\mathbf{u}^{(1)}_{j}(\mathbf{r})$ that overlays $\mathbf{u}^{(0)}_{i}(\mathbf{r})$. We find $\mathbf{u}^{(1)}_{j}(\mathbf{r})$ from the general expression $\mathbf{u}_{j}(\mathbf{r})=\int_{\partial V_{j}}\mathrm{d}S^{\prime}\mathbf{\hskip 0.56917pt\underline{G}}(\mathbf{r}-\mathbf{r}^{\prime})\cdot\mathbf{f}(\mathbf{r}^{\prime})$ by Taylor expanding the Green’s function in $\mathbf{r}^{\prime}$ around $\mathbf{r}^{\prime}=\mathbf{r}_{j}$. The definition of $\mathbf{\underline{S}}_{j}$ shows up as the symmetric part of the second-order term of the series, similarly to the hydrodynamic case Batchelor and Green (1972); Kim and Karrila (1991), leading to

[TABLE]

This expression completes our first step of iteration. In the second step, it is particle $i$ that is exposed to the field $\mathbf{u}_{j}^{(1)}(\mathbf{r})$. Correspondingly, we find its reaction from Eqs. (6)–(9) by replacing $(\mathbf{u}_{i}^{(0)},\mathbf{U}_{j}^{(1)},\mathbf{\Omega}_{j}^{(1)},\mathbf{\underline{S}}_{j}^{(1)},\mathbf{u}_{j}^{(1)},\mathbf{r}_{j})$ with $(\mathbf{u}_{j}^{(1)},\mathbf{U}_{i}^{(2)},\mathbf{\Omega}_{i}^{(2)},\mathbf{\underline{S}}_{i}^{(2)},\mathbf{u}_{i}^{(2)},\mathbf{r}_{i})$. Particle $i$ now feels the consequences of its self-generated field $\mathbf{u}_{i}^{(0)}(\mathbf{r})$ reflected by particle $j$ in the form of $\mathbf{u}_{j}^{(1)}(\mathbf{r})$. Therefore, the procedure was termed method of reflections in hydrodynamics Kim and Karrila (1991); Dhont (1996). The displacement $\mathbf{U}_{i}^{(2)}$ in Fig. 1 results from this reflection. We have not found in the hydrodynamic derivation Dhont (1996) the above reasoning of explicitly imposing on the matrix environment the rigidity-induced stress.

In principle, this refinement of the deformation field via back-and-forth reflections between the two particles can be continued, leading to increasingly-higher-order corrections in $a/r_{ij}$, where $r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$. For our example systems below, these iterations converge quickly, see Fig. 3(c), so that it is sufficient to consider contributions up to (including) order $r_{ij}^{-4}$.

Due to the linearity of Eq. (1), we can sum up the particle displacements obtained from the different steps of iteration. Moreover, we can consider external forces $\mathbf{F}_{i}$ on all particles and calculate the resulting net displacements $\mathbf{U}_{i}$ due to the mutual elastic interactions ($i=1,...,N$). These contributions superimpose. In analogy to the hydrodynamic Dhont (1996) mobility matrix we express the result by an elastic displaceability matrix $\mathbf{\underline{M}}$:

[TABLE]

Limiting ourselves to contributions up to (including) order $r_{ij}^{-4}$, we find

[TABLE]

where M_{0}={}{\color[rgb]{0,0,0}1}/{\color[rgb]{0,0,0}6\pi}\mu a and $\mathbf{\hat{r}}_{ij}=\mathbf{r}_{ij}/r_{ij}$ ($i,j=1,2,...,N$).

In Eq. (11), $\mathbf{\underline{\hat{I}}}$ represents the immediate displacement of particle $i$ due to the force $\mathbf{F}_{i}$ ($\mathbf{U}^{(0)}_{i}$ in Fig. 1). The second term $\sim\!r_{ik}^{-4}$ describes the rigidity-induced reflection from another particle ($\mathbf{U}^{(2)}_{i}$ in Fig. 1). It counteracts $\mathbf{U}^{(0)}_{i}$, or vanishes for $\mathbf{F}_{i}\perp\mathbf{\hat{r}}_{ik}$.

In Eq. (12), {\color[rgb]{0,0,0}\mathbf{\underline{\hat{I}}}+}\mathbf{\hat{r}}_{ij}\mathbf{\hat{r}}_{ij} expresses the consequence of the force $\mathbf{F}_{j}$ acting on particle $j$: particle $i$ is relocated in the induced displacement field ($\mathbf{U}^{(1)}_{j}$ in Fig. 1, for $i\leftrightarrow j$). The term $\sim\!r_{ij}^{-3}$ corrects this displacement field due to the finite size of particle $j$, in analogy to the Rodne-Prager formula in the fluid case Dhont (1996). $\mathbf{\underline{M}}_{i\not=j}^{(3)}$ describes additional three-body interactions $\sim\!r^{-4}_{ij}$ calculated in full analogy to the above procedure for the two-body interaction,

[TABLE]

That is, the deformation field induced by a force on a first particle $j$ spreads to a second particle $k$, from where it is reflected towards the third particle $i$. The angular dependence allows for configuration-dependent attractive, repulsive, or bypass contributions, see Fig. 2.

Eqs. (10)–(13) represent the central theoretical result. Up to (including) order $r_{ij}^{-4}$ it is exact, higher-body interactions for $N>3$ do not enter (see Ref. sup, for $N=4$). To confirm and illustrate the merit of the theory, we performed experiments on small groups of paramagnetic particles embedded in a soft elastic gel matrix. Applying an external magnetic field induced mutual magnetic forces between the particles. Rotating the magnetic field tuned these forces. The resulting relative displacements of the particles were tracked by optical microscopy.

We used paramagnetic Nickel (Ni) particles obtained from Alfa Aesar ($-100+325$ mesh, purity 99.8%). The magnetic hysteresis curves (measured by a vibrating sample magnetometer, Lake Shore 7407) showed a low remanence of $\sim$ 7.5 kA/m, a low coercive field of $\sim$ 2.4 mT, and a volume magnetization of 291$\pm$17 kA/m under an external magnetic field of $\sim$ 216 mT. We carefully selected Ni particles of similar sizes (deviation less than 2% within each group) and a roundness $\gtrsim 0.91$ (measured by image analysis ima ). These particles were embedded in the middle plane of a soft elastic polydimethylsiloxane-based Huang et al. (2016) gel, see Fig. 3(a). First, a bottom gel layer with a thickness of 3.3 mm and a diameter of 24 mm was prepared in a plastic mold. Second, after sufficient stiffening ($\sim$ 0.5 h), the Ni particles were carefully deposited on its top around the center. Third, a top gel layer with the same composition and size as the bottom layer was added. To ensure good connection between the two layers, at least 7 days of cross-linking were allowed.

Using a 32-magnet Halbach array to generate a homogeneous magnetic field Huang et al. (2016), we applied $\sim$ 216 mT to the embedded Ni particles, which is close to saturation. Starting from the initial direction, the magnetic field was rotated clockwise for $180^{\circ}$ in 18 steps within the plane containing the Ni particles. Their center-of-mass positions were tracked by a CCD camera (MATRIX VISION mvBlueCOUGAR-S) with the zoom macro lens (Navitar Zoom 7000) mounted above the samples and subsequent image analysis ima .

We measured the changes in particle distance $\Delta r_{ij}$ ($i\neq j$) for a two- and three-particle system, see Figs. 3 and 4, respectively, when rotating the external magnetic field.

Forces $\mathbf{F}_{i}$ on the particles result from mutual magnetic interactions. Due to substantial particle separations, we approximate the induced magnetic moments as point dipoles Klapp (2005); *biller2015mesoscopic. Thus, we find Jackson (1999)

[TABLE]

with $\mu_{0}$ the vacuum permeability and $\mathbf{m}=m\mathbf{\hat{m}}$ the induced magnetic moments, considered identical for all particles in the close-to-saturating homogeneous external magnetic field. Using as input parameters the experimentally determined particle positions, sizes, and magnetization, we extracted the elastic shear modulus and calculated all changes $\Delta r_{ij}$ from Eqs. (10)–(14). The magnetic forces $\mathbf{F}_{i}$ after displacement are determined iteratively. Perfect agreement between theory and experiment in Figs. 3 and 4 supports the significance of the theoretical approach and highlights its potential for microrheological measurement of the shear modulus.

In summary, we considered rigid spherical particles displaced against a surrounding elastic matrix by externally induced forces. We derived analytical expressions to calculate the resulting particle displacements. Mutual interactions due to induced matrix deformations are effectively included. This renders the procedure a promising tool to describe the behavior of elastic composite materials Ilg (2013); *odenbach2016microstructure; *menzel2016hydrodynamic. Our experiments on paramagnetic particles in a soft elastic gel matrix and subject to tunable magnetic interactions confirm the potential of the theory.

Upon dynamic extension, a prospective application concerns macroscopic rheology Denn and Morris (2014), or nano- and microrheology Ziemann et al. (1994); *bausch1999measurement; *waigh2005microrheology; *wilhelm2008out where the matrix properties are tested by external agitation of embedded probe particles. Also biological and medical questions are addressable in this way, for instance cytoskeletal properties Ziemann et al. (1994); *bausch1999measurement; *waigh2005microrheology; *wilhelm2008out; Mizuno et al. (2007). An extension of the theory to include imposed torques on the particles, e.g., due to magnetic anisotropy, is straightforward and will be presented in the near future.

Supplemental material

As stressed in the main text, the derivation of the displaceability matrix can likewise be performed for compressible systems. Following the same steps of derivation as in the main text, we present below the corresponding expressions for completeness. Apart from that, we add further experimental results and comparison with the theory for a four-particle system, in complete analogy to our presentation for the three-particle system in the main text.

.1 Expressions for a compressible elastic matrix

For clarity and to facilitate the comparison with the hydrodynamic fluid case, we have presented in the main text the expressions for an incompressible elastic system. That is, the system tends to locally preserve the volume of all its volume elements during any type of elastic deformation. However, and in contrast to the hydrodynamic fluid case Batchelor and Green (1972); Kim and Karrila (1991); Dhont (1996), for elastic matrices it is straightforward to allow for compressibility in the derivation. This extended derivation proceeds in direct analogy to the one presented in the main text.

We again assume a homogeneous and isotropic elastic matrix of infinite extension. Once more, we start from the linear elastostatic Navier-Cauchy equations Cauchy (1828),

[TABLE]

As in the main text, $\mathbf{u}(\mathbf{r})$ denotes the displacement field, $\mu$ the shear modulus of the matrix Landau and Lifshitz (1986), $\nu$ the Poisson ratio Landau and Lifshitz (1986), and $\mathbf{f}_{b}(\mathbf{r})$ the bulk force density. Now, we do not restrict our analysis to incompressible materials that locally adhere to $\nabla\cdot\mathbf{u}(\mathbf{r})=0$, and we do not assign a specific value to $\nu$.

The resulting Green’s function for a point force density $\mathbf{f}_{b}(\mathbf{r})=\mathbf{F}\delta(\mathbf{r})$ then reads Landau and Lifshitz (1986)

[TABLE]

Using this expression, if an external force $\mathbf{F}_{i}$ is acting on a rigid spherical particle $i$ of radius $a$ embedded in the matrix with no-slip boundary conditions on its surface, a displacement field

[TABLE]

is induced. Eq. (17) has the same form as in the main text, but $\mathbf{\hskip 0.56917pt\underline{G}}$ is different, see Eq. (16). Again, the validity of Eq. (17) is confirmed as it satisfies Eq. (15), shows the correct limit \mathbf{u}_{i}^{(0)}(\mathbf{r})=\mathbf{\hskip 0.56917pt\underline{G}}(\mathbf{r}{\color[rgb]{0,0,0}-\mathbf{r}_{i}})\cdot\mathbf{F}_{i} for $|\mathbf{r}-\mathbf{r}_{i}|>a$ when $a\rightarrow 0$, and for {\color[rgb]{0,0,0}|\mathbf{r}-\mathbf{r}_{i}|}=a is constant on the particle surface. For {\color[rgb]{0,0,0}|\mathbf{r}-\mathbf{r}_{i}|}=a, it reveals the rigid displacement

[TABLE]

of the $i$th particle in response to $\mathbf{F}_{i}$.

The no-slip condition under our assumptions applies on the surface $\partial V_{j}$ of a particle $j$ also for compressible matrices. Thus Eq. (5) in the main text preserves its shape, i.e.

[TABLE]

where $\mathbf{U}_{j}^{(1)}$ denotes the translation of particle $j$, $\mathbf{\Omega}_{j}^{(1)}$ is its rotation, the displacement field $\mathbf{u}_{i}^{(0)}(\mathbf{r})$ is induced by particle $i$, and $\mathbf{f}(\mathbf{r}^{\prime})$ denotes the surface force density that particle $j$ exerts on the surrounding matrix. The derivation of expressions for $\mathbf{U}_{j}^{(1)}$ and $\mathbf{\Omega}_{j}^{(1)}$ in the form of the Faxén laws follows the same strategy as described in the main text and leads to

[TABLE]

Also the stresslet $\mathbf{\underline{S}}_{j}$ exerted by particle $j$ onto the matrix is derived in analogy to what is described in the main text. In general, for compressible systems, this stresslet is given by the expression $\mathbf{\underline{S}}_{j}=\int_{\partial V_{j}}\mathrm{d}S^{\prime}[\mathbf{f}(\mathbf{r}^{\prime})\mathbf{r}^{\prime}+(\mathbf{f}(\mathbf{r}^{\prime})\mathbf{r}^{\prime})^{T}]/2$. This expression slightly differs from the one introduced below Eq. (7) in the main text for incompressible systems. There, a trace-free definition was used to exclude compressions and dilations of the matrix, which needs to be the case for volume-conserving systems. It can be seen from the main text that the difference in definitions plays no actual role for our derivation. The reason is Eq. (9), where the extra term $\sim\mathbf{\underline{\hat{I}}}$ in the incompressible case only leads to a contribution $\sim\nabla\cdot\mathbf{\hskip 0.56917pt\underline{G}}$. Yet, $\nabla\cdot\mathbf{\hskip 0.56917pt\underline{G}}$ vanishes in the incompressible case. Therefore, following the same strategy as described in the main text, we obtain

[TABLE]

Likewise, the displacement field $\mathbf{u}_{j}^{(1)}(\mathbf{r})$ resulting from the rigidity of particle $j$ and its resistance to deformation, expressed by the stresslet $\mathbf{\underline{S}}_{j}^{(1)}$, is calculated as described in the main text. Eq. (22) here contains a term $\sim 1/(1-2\nu)$, which would diverge for $\nu\rightarrow 0.5$. However, it gets canceled by a counter-factor $\sim(1-2\nu)$ in the calculation. More precisely, upon inserting Eq. (17) into Eq. (22), the expression $\nabla\cdot\mathbf{\hskip 0.56917pt\underline{G}}$ appears; straightforward calculation of $\nabla\cdot\mathbf{\hskip 0.56917pt\underline{G}}$ via Eq. (16) leads to a factor $\sim(1-2\nu)$. In the end, $\mathbf{u}_{j}^{(1)}(\mathbf{r})$ has the same form as Eq. (9) in the main text,

[TABLE]

In the next step, again, the reaction of particle $i$ in response to the field $\mathbf{u}_{j}^{(1)}(\mathbf{r})$ is obtained from Eqs. (20)–(23) by replacing $(\mathbf{u}_{i}^{(0)},\mathbf{U}_{j}^{(1)},\mathbf{\Omega}_{j}^{(1)},\mathbf{\underline{S}}_{j}^{(1)},\mathbf{u}_{j}^{(1)},\mathbf{r}_{j})$ with $(\mathbf{u}_{j}^{(1)},\mathbf{U}_{i}^{(2)},\mathbf{\Omega}_{i}^{(2)},\mathbf{\underline{S}}_{i}^{(2)},\mathbf{u}_{i}^{(2)},\mathbf{r}_{i})$.

Summing up the contributions from the different steps of iteration and considering all $N$ particles simultaneously leads to an expression in the form of an elastic displaceability matrix $\mathbf{\underline{M}}$ as given in the main text:

[TABLE]

Limiting ourselves to contributions up to (including) order $r_{ij}^{-4}$, we find for a compressible system the more general expressions

[TABLE]

where $M_{0}={}(5-6\nu)/24\pi(1-\nu)\mu a$ and $\mathbf{\hat{r}}_{ij}=\mathbf{r}_{ij}/r_{ij}$ ($i,j=1,2,...,N$). Here, the three-body interactions contribute as given by $\mathbf{\underline{M}}_{i\not=j}^{(3)}$ in the form

[TABLE]

The corresponding expressions for incompressible systems in the main text readily follow from Eqs. (24)–(27) by setting the Poisson ratio $\nu=0.5$. Here, we derived and listed the more general expressions for compressible elastic matrices.

.2 Four-particle system

In addition to the two- and three-particle samples described in the main text, we also generated and analyzed four-particle systems. Their preparation, experimental analysis, and the corresponding comparison with the theory are in complete analogy to the three-particle system described in the main text. We recall that our theoretical description in the main text up to the investigated order (including $r_{ij}^{-4}$) is exact for arbitrary particle numbers. No higher-body interactions appear to this order. Therefore, Eqs. (10)–(14) in the main text also apply to systems of particle numbers $N>3$ up to (including) order $r_{ij}^{-4}$, i.e. if the particle separations are not significantly reduced.

Thus, our four-particle results predominantly provide a supplement to the results presented in the main text. Our experimental and theoretical results for the four-particle system are depicted in Fig. 5. One could continue to further increasing particle numbers in the same way.