Magneto-Mechanical Coupling in Magneto-Active Elastomers

In the present work, the magneto-mechanical coupling in magneto-active elastomers is investigated from two different modeling perspectives: a micro-continuum and a particle–interaction approach. Since both strategies differ significantly in their basic assumptions and the resolution of the problem under investigation, they are introduced in a concise manner and their capabilities are illustrated by means of representative examples. To motivate the application of these strategies within a hybrid multiscale framework for magneto-active elastomers, their interchangeability is then examined in a systematic comparison of the model predictions with regard to the magneto-deformation of chain-like helical structures in an elastomer surrounding. The presented results show a remarkable agreement of both modeling approaches and help to provide an improved understanding of the interactions in magneto-active elastomers with chain-like microstructures.


Introduction
Field-controllable functional polymers represent a relatively new class of applied materials exhibiting a strong coupling of mechanical and additional-e.g., electric, magnetic, or thermal-external fields. The application of these external fields influences the interactions between different local material phases and causes an evolution of the microstructure. A prominent example of field-controllable functional polymers are magneto-active elastomers (MAEs) which feature mechanical moduli that can be enhanced under an applied magnetic field [1][2][3][4][5][6][7][8][9][10] as well as the ability for magnetically induced deformations [2,3,11] and actuation stresses. These properties make MAEs attractive for a variety of technical implementations: up to now, applications for actuators and sensors [12][13][14], energy harvesting [15][16][17], micro-robots [18], and -pumps [19] as well as prosthetic and orthotic devices with tunable stiffness [20] have been proposed.
In all of these applications, the effective magneto-mechanical material behavior is of special interest. To this end, an in-depth understanding of the materials' structure-property relationships is required. MAEs typically represent composite materials, in which micronsized magnetizable particles are embedded into a compliant polymer network [21][22][23][24][25][26]. The flexibility of polymer sub-chains between cross-links allows for a considerable degree of particle rearrangement under strong magnetic fields. This is especially true for sufficiently soft matrices, in which the particles are prone to organize themselves into elongated microstructures, that significantly influence the coupled magneto-mechanical properties of MAEs. The simultaneous consideration of inhomogeneous magnetic and mechanical fields and their interaction with the underlying, evolving microstructure is the particular challenge in the modeling of MAEs.
In order to provide a better understanding of the effective material behavior of MAEs, different theoretical approaches have been presented in the literature over the last few years. Following the schematic illustrations in Figure 1, they can be divided into particleinteraction as well as micro-and macro-continuum models. With few exceptions [27,28], particle-interaction models, see Figure 1a, are based on a calculation of the effective MAE behavior from overall energetic expressions that can be found using the dipole approximation, i.e., by not resolving the local magnetic fields [29][30][31][32]. Consequently, the approach allows for considering MAE specimens comprising a large number of particles with little computational effort. On the other hand, the method provides only an approximate description of the particle interactions within the material and, thus, is especially suited to characterize dilute systems with a low particle-volume fraction. As indicated in Figure 1b, micro-continuum models pursue a different approach: here, local magnetic and mechanical fields are resolved explicitly using a continuum formulation of the magneto-mechanical boundary value problem [33][34][35]. Since this yields a system of fully coupled, nonlinear partial differential equations which require computationally demanding numerical methods-such as a finite element (FE) analysis-for their solution, only comparably small MAE samples as well as representative microstructures in the surrounding of a material point can be considered. To this end, the approach allows for identifying mechanisms leading to magnetically induced deformations and enhanced mechanical moduli on the microstructural level but requires an appropriate computational homogenization procedure to predict the effective material behavior [36][37][38][39]. In macro-continuum approaches, the MAE is modeled as a homogeneous continuum in which microstructural information are captured via suitable coupling terms; see Figure 1c. Here, a phenomenological material model is typically fitted to data obtained from experiments [40,41] or more resolved, microscopic analyses [42][43][44]. Due to its phenomenological nature, this strategy cannot provide any understanding of microscopic mechanisms which drive the materials' effective response, but-as it allows a consideration of realistic MAE samples under complex loading conditions-yields relevant information on macroscopic shape effects [45].  Figure 1. Illustration of different modeling strategies for MAEs: (a) particle-interaction models with dipolar magnetic particles (red arrows) distributed in an elastomer matrix (light gray background), (b) micro-continuum models with fully resolved particles (dark gray circles) within the same matrix, and, (c) homogeneous macro-continuum models with no resolution of the underlying microstructure.
Considering experimental investigations on MAEs, their effective properties seem to be ambiguous. Especially for the magnetically induced deformation-a phenomenon which here is referred to as magneto-deformation to avoid confusion with the intrinsic magneto-striction of single-phase materials-samples with comparable shapes and particlevolume fractions are reported to show either a contraction [11,46] or an elongation [40,47,48] in the direction of the applied magnetic field. This shows that micro-and macrostructural effects require to be taken into account for an understanding of the material. Consequently, both structural levels have to be captured in the pursued modeling strategy. To allow this, a combination of all of the aforementioned approaches is advantageous: the microscale can be considered performing fast analyses with particle-interaction models while using the obtained results for a fitting of a macro-continuum model enables investigations of realistic MAE samples under complex loading conditions. However, due to their accuracy, microcontinuum models are still required to define the limits of the simplified particle-interaction approach-only by comparing both strategies for the range of targeted applications and microstructures, a reasonable characterization of MAEs can be ensured.
On the way towards the development of such a hybrid multi-scale approach, the main focus of the current contribution is a comprehensive comparison of the two micromodeling strategies. Up to now, these approaches mainly coexisted and have only been compared for simplified two-dimensional MAEs in a preceding article of the authors [49]. Here, this work is extended to the general three-dimensional case with an emphasis on chain-like microstructures. An analysis of ideal helical structures allows for investigating the magnetic response of chains with different angles and distances between the individual particles and, thus, facilitates an identification of distributions which are more prone to positive magneto-deformation, i.e., an elongation, than others. In view of the fact that structured MAEs are of significant interest within the research community, these results can help to provide an improved understanding of realistic chain-like structures.
The paper is organized as follows: in Section 2, the applied micro-modeling strategies, i.e., the micro-continuum and the particle-interaction models, are presented-relevant equations are specified and the capabilities of the individual approaches are summarized using representative examples. In the subsequent Section 3, the comparison of both approaches regarding their predictions for the magnetically induced deformation of helical chains is performed: after the problem under consideration is outlined, the results and their implications for the modeling of MAEs are discussed. Finally, the paper is concluded by a short summary in Section 4.

Micro-Modeling Strategies
In order to clarify the common basis of the models presented in this section, Figure 2 can be used: it illustrates an MAE sample with a cylindrical shape which is frequently used in experiments [50][51][52]. In the vicinity of each material point, the microstructure is assumed to consist of stiff, magnetizable particles surrounded by a soft and non-magnetizable elastomer matrix. In Figure 2, a monodisperse random microstructure with spherical particles is shown-however, this distribution can vary depending on the used materials and the sample preparation process [11,53]. The experimentally observable magneto-deformation as well as changes of the samples' moduli are a result of mutual particle interactions caused by an external magnetic field H ∞ . Since, within this contribution, only magnetically soft materials are of interest, the local magnetization can be described via functions of the form which relate the magnetization M and magnetic field H at any material point r. Note that M and H are always aligned-their connecting function M represents any saturating function, e.g., of Langevin-, hyperbolic-tangent-or Frölich-Kennelly-type, which only depends on the norm of the magnetic field. For the mechanical properties of the matrix material, a purely elastic behavior is assumed. Its implementation within the individual approaches is presented in the subsequent subsections.

Micro-Continuum Model
The foundation of the micro-continuum model presented here is the work of de Groot and Suttorp [54]. Therein, long-and short-range contributions of atomic interactions are investigated for various macroscopic physical quantities-the resulting field equations feature additional coupling terms and are obtained by performing statistical averages to a level where the continuum hypothesis can be considered valid. For the interaction of stationary magnetic and mechanical fields, the governing equations will be briefly introduced, in the following.

Governing Equations
The magnetic part of the coupled problem requires Maxwells equations to be solved [55]. Assuming a material body which is free of any current densities and surface currents, they are given by and their associated jump conditions [56]. Here, B and H represent the magnetic induction as well as the magnetic field of a material point in the current configuration [57]-they are linked by the magnetization M via the equation with µ 0 being the permeability of free space. As shown in [56,57], the coupling of the mechanical fields into the magnetic equations can be pointed out by defining the reference fields using the deformation gradient F and its determinant J: Finally, a scalar potential approach with H = −∇ϕ is applied within this contribution to reduce the number of unknowns and equations to be solved [55].
For the mechanical part of the problem, the existence of magnetic field yields additional body force f pon = (∇B) T · M, couple c pon = M × B and power p pon = M ·Ḃ densities which have to be accounted for within the individual balance equations. Especially the additional couple density has severe consequences for the solution of the coupled magneto-mechanical boundary value problem: the mechanical stress σ is not necessarily symmetric. In order to ease the computation, the magnetic body force density can be expressed as the divergence of a ponderomotive stress via f pon = ∇ · σ pon with and I being the identity tensor [57]. This mathematical trick facilitates the introduction of the symmetric total stress σ tot = σ + σ pon and-if vanishing mechanical body force densities are assumed-leads to the following balance of linear momentum: Since thermal effects are not within the scope of this work, the balance of the internal energy yields no additional information and is omitted here for brevity, see [58] for further information. For the derivation of thermodynamically consistent constitutive models, only the Clausius-Duhem inequality-as a consequence of the entropy balance and the second law of thermodynamics-is still required. With C = F T · F being the right Cauchy-Green deformation tensor, the amended free energy extends the Helmholtz free energy Ψ by free field magnetic contributions [57]. Using this relation, the Clausius-Duhem inequality can be stated as:

Constitutive Models
In order to ensure constitutive models that are not only thermodynamically consistent but also objective, quantities of the reference configuration are chosen as independent variables. Here, the right Cauchy-Green deformation tensor C as well as the reference magnetic field H 0 are used. By performing an evaluation of the Clausius-Duhem inequality (9) according to the procedure of Coleman and Noll [59], the constitutive relations for the magnetic induction and the total stress can be found: Since both constituents in MAEs, the stiff magnetizable particles, and the non-magnetizable elastomer matrix show intrinsic magneto-mechanical coupling effects that are vanishingly small if compared to the compounds coupling behavior, the free energy term in Equation (8) can be split into purely magnetic and mechanical contributions: For the latter, an isotropic hyperelastic energy function is adequate-it allows for accounting for potentially large deformations of the matrix material and to capture the resulting rotations of the particles. The choices for such functions vary from simple neo-Hookean to elaborate Ogden models which allow for accurately describing a nonlinear elastic behavior over a wide range of deformations [60]. For the different results presented in this contribution, the following energetic contributions are applied: Therein, the quantities I k and I iso k represent the principle invariants of C and its isochoric counterpart C iso = J − 2 3 C, λ iso β are the eigenvalues of F iso = J − 1 3 F with algebraic multiplicity ν β and the sets µ p , α p , κ , {µ 1 , µ 2 , κ} as well as {µ, λ} are corresponding material parameters which are specified within the individual examples.
For the magnetic contribution of the Helmholtz free energy, Ψ mag = 0 holds within the matrix material, while the experimentally observed magnetization curves of carbonyl iron powder and nickel particles as described in [61,62] are applied to capture the nonlinear magnetization behavior of the particles. Thus, Ψ mag is given by Hyperbolic tangent: Ψ mag 2 Hyperbolic tangent:

Computational Homogenization Framework
In order to be able to determine physically meaningful macroscopic quantities from numerical simulations on the microscopic scale, an adequate scale transition has to be performed. Following the work of Hill [63] and its extension to the coupled magnetomechanical boundary value problem [36], the Hill-Mandel condition, i.e., the equivalence of the macroscopic and averaged microscopic energies, is applied here. If the averaging process of a microscopic quantity (·) is indicated by (·) and its macroscopic counterpart is denoted as( ·), it is given by Therein, P tot = JF −1 · σ tot is the first Piola-Kirchhof stress tensor. Within the FE simulations, the fulfillment of the condition is ensured by making use of representative volume elements (RVEs) for which periodic boundary conditions are applied [56].

Model Applications
In the following, the versatility of the presented micro-continuum approach will be outlined by means of three representative examples. In order to ensure that the strategy is eligible for the modeling of the strongly coupled magneto-mechanical behavior of MAEs, it is validated with experimental data for a simplified MAE specimen in the beginning. Afterwards, recent findings for the field-dependent behavior of MAEs with different ideal and random microstructures are presented. Since the micro-and macroscales need to be considered to capture all effects of realistic MAE samples, a first approach towards the development of a hybrid multiscale modeling strategy by fitting a macro-continuum model with data generated from micro-continuum simulations is briefly summarized as a last example.

Model Validation by Means of a Simplified MAE Specimen
While the material behavior of the MAE's constituents is often well-known, macroscopic effects of inhomogeneous magnetic and mechanical fields make an experimental characterization of their compound behavior almost impossible [42,64]. To this end, a proper validation of the micro-continuum model can only be performed with well-defined, simplified MAE samples under controlled loading scenarios. In a previous study, this procedure was applied for systems with only few magnetizable particles [62]. It allows for a detailed analysis of the sample behavior with a limited number of influencing factors. Here, the results for the three-particle sample in the aforementioned work are summarized.
A schematic illustration of the system under consideration is shown in Figure 3a: three magnetizable nickel particles with diameters of approximately 200 µm were placed into the center of a Polydimethylsiloxane (PDMS) matrix with a quadratic cross section of 15 mm length and a thickness of 4 mm, see [62] for further information on the sample generation. During the experiment, a magnetic field B 0 with norm B 0 = 170 mT was applied. While its angle β was changed in steps of ∆β = 5°, the inter-particle distances d ij of the magnetizable inclusions were tracked using a microscope. The results are shown in Figure 3b. Within the FE simulations, the experimental setup was mimicked as closely as possible: only the shape of the magnetizable inclusions has been approximated by spheres-this results in a symmetry of the specimen with respect to the x 1 -x 2 -plane. For the FE simulations, the displacement of all outer specimen boundaries was fixed, while inhomogeneous Dirichlet boundary conditions were used for the magnetic scalar potential to apply the external magnetic field. Regarding the constitutive behavior, the energetic contributions introduced in Equations (11b) and (12b) have been used. For the magnetizable nickel particles, experimental investigations using a novel approach with a superconducting quantum interference device have shown that M s = 314.5 kA/m and ξ = 1.886 · 10 −2 m/kA are a reasonable choice to describe their magnetization curve, while µ p 1 = 76 GPa, µ p 2 = 0, and κ p = 164 GPa are applied to account for their elastic behavior. Since no full experimental characterization of the elastomer matrix was available, the numerical simulations have been performed for ν m = 0.49 and various matrix moduli-the optimal value E m opt = 6842 Pa was determined in a least-squares sense by comparing the experimental and numerical predictions for every distance and angle.
The comparison of the experimental and numerical results for the best fit in Figure 3b shows a very good qualitative and quantitative agreement: for all inter-particle distances d ij , the sample behavior is captured by the presented micro-continuum approach which demonstrates that it represents an adequate framework to describe microstructural interactions in MAEs. Variations of the particles initial positions that were also performed in [62], indicate a sensitivity of the simulations with regard to the sample geometry. To this end, the systematic error for d 23 and the small deviations for the distance of particles 1 and 2 can be eliminated by incorporating data from micro-tomography measurements [65] for the generation of more accurate simulation models.

Field-Dependent Behavior of MAE Microstructures
As a field-induced change of the macroscopic stiffness, the magneto-rheological (MR) effect is one of the most frequently investigated effects in MAEs [1,8,66,67]. In such experiments, MAE samples with different shapes and microstructures are typically exposed to a shear deformation for varying external magnetic fields: by comparing the modulus to the one of the magnetically unloaded sample, the MR effect is determined. On the microstructural level, the effect can be ascribed to attractive and repulsive forces between the interacting magnetizable particles which means that-apart from macroscopic shape effects due to, e.g., inhomogeneous magnetic fields-different particle distributions must entail differences in the materials MR effect as well.
In a recent study on the field-dependent behavior of MAEs, two microstructures have been investigated exemplarily regarding their influence on the MR effect-a lattice-like simple cubic and a more realistic random microstructure [68]. For the analysis, the experimental setup was mimicked within a finite element simulation: representative volume elements for the materials' microstructure were exposed to a macroscopic magnetic field H 0 while their macroscopic deformation was fixed. At the end of each load step, a small shear deformationγ was applied to determine the shear components of the macroscopic, i.e., compound, stiffness tensor. Using the nominal macroscopic mechanical stress tensorP and a loading of the following form: the analyzed macroscopic shear componentḠ and its relative change ∆Ḡ are defined by: To obtain the results depicted in Figure 4, the energy functions (11c) and (12a) have been employed. The stiff magnetizable particles are characterized by M s = 868 kA/m, ζ = 2.18 · 10 −2 m/kA as well as λ p = 121 GPa and µ p = 81 GPa. For the elastomer matrix, λ m = 1644 kPa and µ m = 34 kPa are applied. As can be seen from the results for the cubic microstructure in Figure 4a, a field-induced stiffening is found for all particle-volume fractions φ. All curves show the typical behavior of an initially quadratic course with a saturation for higher magnetic fields [56,69]. Qualitatively, the results are in good agreement with findings of other authors [70,71]. Quantitatively, an increase of the modulus by up to 40 % is observed with moderate particle-volume fractions-this effect can even be increased by considering a softer elastomer matrix as it is often used in realistic MAE samples.  Regarding the RVEs comprising 100 randomly placed inclusions with the same particle-volume fractions in Figure 4b, the situation is different. Although many experimental as well as numerical investigations report a field-induced increase of the materials modulus for such microstructures, a decreasing stiffness is found within this study. From the qualitative point of view, the course of the curves is similar to the one observed for the cubic microstructure. However, it is apparent that the MR effect is not systematic for the analyzed RVEs: the microstructure with φ = 10% shows a bigger effect than the one with φ = 15%. As will be pointed out in the subsequent example, recent findings of the authors show that-on the microstructural level-the effects in MAEs are driven by clusters of interacting particles which are very close to each other. Minimum distances, here 27.5% of the particle diameter, which act as a limit to impede intersecting particles during the random placement algorithm in the geometry generation can, as in this and other examples [72], entail such an unsystematic behavior, i.e., produce microstructures that are random but apparently not representative. Concerning the missing data points for φ = 5% as well as φ = 25% another drawback of the applied modeling strategy in this example can be seen: for increasing magnetic fields, the simple neo-Hookean material model (11c) is not able to resist the attractive forces between close inclusions, see [73] for further information. To this end, the numerical simulations start to fail.

Macroscopic Model Calibration Using a Decoupled Multiscale Framework
To facilitate a full characterization of MAEs on the microscopic and macroscopic scales, a multiscale framework is required. In the following, the procedure according to [42,68], where a macro-continuum model is calibrated from micro-continuum simulations for the simplified two-dimensional case, is briefly presented. The framework represents a decoupled multiscale scheme which-in contrast to, e.g., FE 2 -approaches [39,45]-can only account for microstructural evolution within the limits of the model calibration. Since isotropic MAEs are considered, it is sufficient to describe the macroscopic behavior by the averaged principal stretchesλ iso α and the three additional magneto-mechanical invariants ChoosingH 0 as the independent magnetic quantity, the amended free energy for the macroscopic model can be expressed asΩ =Ω(λ iso β ,J,Ī 4 ,Ī 5 ,Ī 6 ). Similar to the discussion in Section 2.1.2,Ω can be divided into purely magnetic and mechanical parts-however, since coupling effects cannot be neglected on the macroscopic scale, such an approach requiresΩ to take the form: Based on this, a possible model which is able to account for a strongly nonlinear behavior of the elastomer matrix as well as magnetic saturation effects of the particles is given by the following energetic contributions: For the fitting of the sets of macroscopic material parameters {μ p ,ᾱ p ,γ k ,δ k ,κ}, the systematic procedure described in [42,68] is applied: a total number of nine calibration and five validation load cases is performed with RVEs comprising 400 randomly distributed particles at different particle-volume fractions. To prevent an unsystematic RVE behavior, the minimum particle distance during the RVE geometry generation using a random sequential adsorption algorithm is set to 5% of the particle diameter. For the microscale simulations, the constitutive models according to Equations (11a) and (12a) were used-for the values of the individual parameters and the quality of the model calibration, the authors refer to [42,68].
To demonstrate the model capabilities, the magneto-deformation of a circular macroscopic sample with a particle-volume fraction φ = 30% is investigated. In order to mimic the situation of a homogeneous far field, the sample is embedded into a sufficiently sized free space with negligible mechanical properties and the magnetic properties of free space, i.e.,μ r = 1. Regarding the prediction of the magnetically induced strainε = ∆l/l 0 given in Figure 5a, an initially quadratic behavior which saturates for large external fieldsH ∞ is observable. As already shown in the previous example, this response is typical for MAEs and qualitatively coincides with experimental findings. A closer look on the embedded contour plot of the sample reveals a highly inhomogeneous macroscopic deformation field C: it results from the jump of the magnetic quantities on the MAE surface and illustrates the specific difficulties of experimental investigations on MAEs. Even macroscopic samples with a shape that allows for a homogenous magnetization cannot provide the basis for a shape independent material behavior. Magneto-deformation of a circular MAE sample with φ = 30% particle-volume fraction: (a) predicted magnetically induced deformation of the sample and contour plot of the macroscopic deformation field, and (b) comparison of the calibrated macro-and averaged micro-contnuum models regarding their predictions for the macroscopic mechanical stress at the center of the sample. Data taken from [74].
In order to show the quality of the calibrated model, a localization step is performed in Figure 5b: for the center point of the MAE sample, the predictions of the macroscopic model with regard to the mechanical stress are compared to the averaged results of a microscopic simulation with an RVE that is subjected to the same macroscopic fieldsH 0 andF. The high accuracy of the macro-continuum model over the whole range of applied magnetic fields not only documents that it implicitly captures microstructural information, but also shows that the proposed framework represents an appropriate multiscale scheme which is a computationally efficient alternative to conventional FE 2 -approaches.

Dipole Approach to Particle Interactions
To capture the leading effects of long-and short-ranged magnetomechanical coupling responsible for the effective behavior of MAEs in an applied magnetic field, several approximations are introduced. The major simplification consists of the dipole approximation. In this framework, the magnetic interactions among filler particles are reduced to mutual dipole-dipole interactions. Thus, instead of resolving the, in general inhomogeneous, magnetization field in the volume of each individual inclusion, the problem is simplified to a single relation between the particles' center positions. This reduces the computational effort considerably and allows for an approximate description of magnetic interactions among very large numbers of particles.

Effective Macroscopic Behavior-Minimum of Global Energy
An adequate description of field-induced changes in the shape and mechanical properties of polymer networks with embedded magnetizable microparticles represents quite a challenging problem, even in the framework of the dipole approach to particle interactions. It is relatively easy to consider the influence of particle distribution on static properties [71,75] and dynamic moduli [76,77], using different structures defined on infinite lattices. However, the initial shape of the sample is proved to have a comparable and even a predominant effect on the change in shape [4,29], for example in MAEs with an isotropic distribution of the particles. The close correlation between the local particle structure on the microscale and the global geometric shape of the sample on the macroscale is already apparent from the magnetostatic laws alone. This fundamental interrelation of short-ranged structure effect and long-ranged shape effect makes it necessary not only to be able to describe both effects independently of each other, but also to combine them fundamentally in a unified approach in order to achieve a comprehensive characterization of MAEs.
Such unified approach in the context of the dipole approximation has been originally developed for linearly magnetizable particles placed on the sites of different lattices [29]. In that study, the effective magneto-mechanical behavior has been predicted by minimizing the free energy functional: where U el is the elastic energy of a deformed MAE sample due to the entropic elasticity of polymer chains and U mag arises from the dipole-dipole interactions between magnetizable particles placed in an external magnetic field. This field is assumed to be homogeneous and hence direct interactions of the particles with the field may be skipped as they do not contribute to the minimization procedure. The elastic contribution U el represents the effective mechanical properties of the composite in the absence of any magnetic fields, and it is usually implemented via a neo-Hookean material model. In the linear magnetization regime, i.e., M = χH with the magnetic susceptibility χ, the magnetic energy can be represented in a rather concise form [29]: The strength of the applied magnetic field is denoted by H ∞ . In the above equation, two scalar parameters f micro and f macro are introduced, with the first one quantifying the short range structure effect in regular, i.e., lattice, structures, and the second one the long-range shape effect. Note that these terms appear in the denominator of Equation (23), together with χ and demagnetizing factor n of a single particle, due to a self-consistent treatment of magnetic interactions. Minimizing Equation (22) for different lattice structures (simple cubic, body-centered cubic, hexagonal close-packed), it was found that the magnetomechanical behavior very sensitively depends on the particular choice of lattice parameters upon introducing unrealistically strong long-range ordering among the particles.
Apparently, presuming a perfect ordering which is intrinsic in the lattice structures may lead to undesirable artifacts. This problem can be eliminated by complementing the unified approach with a completely different characterization of particle distributions. An adequate description of practically relevant microstructures has been achieved in the mean-field version of the unified approach [32] via introduction of a continually varying density field, both for stochastically isotropic and elongated microstructures. The methodology developed for this purpose allows in part analytical solutions, on the basis of which the behavior of linearly magnetized spheroidal MAE samples was estimated over large parameter spaces and a wide range of situations. This allowed to draw comprehensive phase diagrams of the deformation behavior of MAE samples with elongated, or columnar, microstructures. On the left side of Figure 6, we schematically illustrate the results. Additionally, a discontinuous shape change for very oblate samples has been predicted, see Figure 10 from [32]. For the deformation behavior of samples with stochastically isotropic particle distribution, a good agreement with more detailed 2D continuum-mechanical simulation calculations was recently demonstrated [49].  [32]. Whether an MAE sample contracts or elongates in direction of the applied magnetic field crucially depends on the sample form (oblate or prolate), the thickness of the columnar microstructure and number density of such structures in the sample, i.e., the overall amount of magnetizable particles. Right side: Systematic comparison between theory and experiment allowed for identifying the contribution from macroscopic sample shape and microscopic particle structure to the magnetically induced stress in confined MAE samples. Reproduced from Ref. [78] with permission from the Royal Society of Chemistry.
Earlier theoretical considerations in the framework of the dipolar approach were restricted to the assumption of affine deformation on local scales [29,32]. Non-affine displacement of magnetic particles due to mutual hindrance of neighboring particles and its effect on the magneto-induced deformation was analyzed in a general study [79]. Allowing the particles to rearrange in the presence of magnetic field, the enhancement of the magnetoinduced deformation in linearly magnetized spheroidal samples has been predicted. The unified approach has been further developed in order to calculate and analyze the interplay between short-ranged structure effect and long-ranged shape effect for the nonlinearly magnetized samples of realistic shapes (cylinder, cuboids). In such generalized formulation, the governing equations turn into tensorial forms [78]. In particular, it was shown that the average magnetization M in the axially symmetric samples, with the external field applied along the axis of symmetry, can be compactly expressed in scalar form as [78]: Here, M denotes the magnetization function, see Equation (1), f macro defines the demagnetizing factor of a sample along its axis of symmetry and f micro describes the actual, in general arbitrary or irregular, microstructure. This theoretical development allowed for quantifying the contribution of microstructure effect into the magnetically induced stress in confined MAE samples, as described in detail in [50,78]. Clear trends for the isotropic and structured samples have been established by fitting the theory predictions, see the right side of Figure 6, to the measured stress data in a specially designed experiment [78].
The presence of external magnetic field transforms initially isotropic MAE samples into transversely isotropic ones with an axis of symmetry defined by the direction of magnetic field. Uniaxial deformations applied along and perpendicular to the field direction lead to a strong anisotropy in the magnetically induced stress response, as predicted in a very recent work [80]. In this study, an attempt was made to derive an effective material model from the free energy functional based on the dipole approximation for magnetic interactions. Importantly, a strong magneto-mechanical coupling between the internal magnetic field and the sample shape is treated self-consistently from the very beginning, being a direct outcome of the minimization of free energy (Equation (22)). This is in contrast to the phenomenological models, in which the form of this coupling is postulated through introduction of multiple pseudo-invariants.
The unified approach is based on the minimization of free energy and is therefore restricted to description of static magneto-mechanical behavior. Moreover, the magnetoinduced changes in the microstructure can be captured only on average by considering an evolution of a single scalar parameter f micro . To be able to describe these changes in more detail, we also considered an alternative approach to the magnetic and elastic interactions, as will be described in the next subsection.

Explicit Particle Structures-Balance of Local Forces
Naturally, the energetic minimum of Equation (22), yielding the equilibrium state, corresponds to the balance of magnetic and elastic forces. However, to obtain analytic relations directly from the energy based approach, additional constraints were introduced. For example, in the case of affine elastic coupling to a uniaxial macroscopic deformation, it is possible to treat the 'collective' repositioning in terms of only one single deformation parameter, i.e., the relative length change of the sample in direction of applied field, ε = ∆L L . In order to account explicitly for the local restructuring under the action of an applied field, it is more convenient to calculate the mutual magnetic forces as driving forces for the particle rearrangements. Thus, as an alternative to the effective energy approach discussed before, the magneto-mechanical problem can be formulated equivalently in terms of magnetic and elastic forces acting on the individual filler particles. In the literature, different models based on mutual dipole-dipole forces in bead spring or continuum approaches have been introduced to simulate the behavior of MAEs [81][82][83][84][85]. In the following, a force-based formulation for a sample of filler particles will be specified in close correspondence to the preceding subsection.
In principal, the magnetic forces are obtained as the negative gradient field of the magnetic energy. Accordingly, the starting point here is the general equation for the magnetic energy of a material/sample of volume V s in some external field H ∞ [32,86]: In case of an MAE, representing a composite of non-magnetizable matrix, where M = 0, with N embedded magnetic/magnetizable inclusions, the magnetic energy becomes a sum of integrals over the corresponding inclusion volumes. The relation between magnetization M and local magnetic field H is given via Equation (1).
In the dipole approach, we assign each particle, labeled i ∈ N, an individual dipole moment m i located at its center position. This dipole moment represents the average magnetization of the inclusion: Here, v i denotes the volume of particle i and the average · i in the r.h.s. means that it has taken over this volume v i . Assuming a homogeneous magnetic field over v i , and thus also a homogeneous magnetization of the sphere, the dipole field exactly describes the magnetic field generated by the particle in its surrounding.
The dipole approximation is equivalent to the assumption that magnetic field H(r) can be considered as constant over the extent of each particle. Since the external field is usually considered homogeneous, or will not change notably on the range of micron-sized particles anyways, the error of the dipole approximation is essentially due to inhomogeneities of the demagnetization fields among neighboring particles. It has been shown [27,86,87] that the dipole approximation model starts to deviate from exact calculations only as particles come closer than ∼ 1.5 d p to each other, where d p denotes a particle diameter.
In the following, center positions of the particles are denoted by r i and we introduce the short notation X i = X(r i ). Then, applying the dipole approximation the volume integral in Equation (25) can be performed straightforwardly over the volume of each inclusion. The result is a sum over all N dipoles: Note, due to Equation (1), we immediately find m i H i , whereas H ∞ is not necessarily aligned with the dipole moments. In the dipole approximation, Equation (1) becomes: For a sample with given macroscopic shape and particle distribution, we obtain from Equation (27) the magnetic contribution to the total energy. Note that the application of Equation (27) requires the self-consistent solution of a set of 3N coupled equations of the form of Equation (28). As outlined in Section 2.2.1, the energetic description is useful to calculate the equilibrium state and also serves as the basis for an effective material model describing for example finite deformation processes.
It was shown explicitly [87] for two equally sized particles with a linear magnetization behavior, i.e., M(H) = χH, that the magnetic forces obtained as the gradient field of the magnetic energy in Equation (27) are identical to the well-known dipole-dipole interaction forces. In Appendix A, we demonstrate the validity of this identity for any number of differently sized particles assuming an arbitrary magnetization function M(H) in Equation (1). In the dipole approximation, the magnetic force f d k acting on some particle k in a sample of N magnetized inclusions is given as: Here, r ki = r i − r k denotes the vector connecting particles k and i and r ki = |r ki |. Once the magnetization, i.e., the magnetic dipoles m i , of the particles are specified, the magnetic forces acting on each of them are also known, and the result is consistent to the energy based approach.
The local field H i is the sum of the external field, the (self-)demagnetization field and the dipole fields of all the other particles: Consequently, to obtain the individual m i , we require the self-consistent solution of a 3N system of coupled nonlinear equations as formulated through Equations (28) and (30). This calculation is performed numerically via classical Newton-Raphson technique providing a fast and stable convergence after only few iteration steps.
In the energetic formulation, we usually implement a macroscopic mechanical model to describe the elastic properties of the sample. That is, we introduce effective moduli of the composite material in the absence of magnetic fields, i.e., for the filler reinforced matrix at given volume fraction of embedded hard particles. In contrast, to account for particle rearrangements on local scale we require a relation between the displacement ∆r i of some individual particle and the set of applied forces {f k , k ∈ N} that are acting on the inclusions. The embedding medium shall thereby be characterized by its 'pure' elastic modulus, containing a given distribution of surrounding hard particles. Using linear elasticity theory, Puljiz et.al. obtained an analytic expression considering in first approximation mutual 2and 3-body elastic interactions among spherical particles [30,88]. In our present work, we assume an isotropic and incompressible elastic medium. Accordingly, we set ν = 0.5 and 3µ = E, where ν denotes the Poisson ratio, µ the shear modulus, and E the elastic modulus. From [88], we then find for each particle in the sample a linear dependency of the form: Note that subscript 0 refers to the initial particle positions in the absence of any forces, i.e., in the absence of magnetic fields, and f j denotes the force acting on the spherical inclusion j. The so defined tensors D ij are dimensionless and read for i = j: and for i = j: Here, I represents the identity tensor. We denote the dyadic product of two vectors a and b in the form ab, whereas a · b implies the scalar product. The first term in Equation (33) describes the single particle contribution due to an elastic medium without any other inclusions. The last term in Equation (32)

accounts for 3-body hydrodynamic interactions and the rest in Equations (32) and (33) refers to corresponding 2-body interactions.
Thus, the final equilibrium positions r i = r 0 i + ∆r i when applying some homogeneous H ∞ are found upon insertion of Equation (29) in Equation (31) assuring at the same time the self-consistency of the magnetic dipoles of all particles with respect to these new positions via Equations (30) and (28). The numerical, or computational, effort for the calculation increases with order O((3N) 3 ) and becomes quite cumbersome for large systems (N 1). The bottleneck in the computation arises from the 3-particle hydrodynamic interactions, i.e., the last term in Equation (32). However, due to its reference to initial positions, the contribution for each D ij only needs to be specified once for a given system. Up to N = 1000 particles, the solution of this algorithm is moderately accessible (O(10 2 seconds)). To test the accuracy of the present approximation approach, we will consider finite samples with few embedded particles in Section 3 and compare the results to the refined continuum description.
The above defined algorithm provides a static process from initial particle positions r 0 i (the prepared equilibrium in the absence of magnetic fields) to final positions r i (the equilibrium in the presence of magnetic fields). Alternatively, in Appendix B, we suggest a simulation model based on the force balancing description of the problem. This extension follows in a straightforward manner upon insertion of corresponding forces into the fundamental Newtons equation and thereby allows for including dynamical aspects of particle repositioning. Furthermore, in contrast to the micro-continuum approach and the static force balancing approach, the dynamical simulation model can naturally prevent the overlap or collapse of neighboring particles in very soft matrices upon introducing a hard sphere repulsion term (see Appendix B). In Figure 7, we display two snapshots of the dynamical simulation for particle rearrangement. The matrix is chosen to be so soft that neighboring particles would collapse into each other due to predominant magnetic forces over elastic ones. In the dynamical simulation model, the particles stick to form pairs/clusters; see the right side in Figure 7. Under such conditions, both aforementioned approaches would fail. As long as the elastic matrix is stiff enough to prevent such collapsing events, the final equilibrium positions of the particles when applying an external magnetic field are identical to the static force balancing approach. Example videos of rearrangement processes, as calculated via the suggested simulation model specified in Appendix B, are attached in the Supplementary Materials.

Comparison of the Approaches
In this section, the presented micro-modeling strategies are compared with regard to their predictions for the magneto-deformation of MAE samples with a specific particle distribution. Since chain-like microstructures are of significant interest in both experimental [8,90,91] and theoretical [34,79,92,93] investigations, idealized samples comprising helical chains with different angles between two neighboring particles are analyzed. Such structures represent the three-dimensional generalization of wavy chains, which, in the simplified two-dimensional case, have been found to be a possible microstructural explanation for the often observed positive magneto-deformation of MAEs with chain-like microstructures [94]. However, recent findings of the authors [56,79] show that two-dimensional analyses tend to overestimate the resulting effects in MAEs and experimentally observed effects cannot be explained with interactions in a simplified wavy microstructure. To this end, the subsequent investigation not only compares the two micro-modeling strategies but can also help to understand why some particle arrangements are more prone to positive magneto-deformation than others and serve as a basis for studies on more realistic, randomized chain-like structures.
Within this study, helical chains with angles ∆ϑ = {45°, 60°, 90°, 180°} between adjacent particles are embedded into an elastomer surrounding which has a size that is sufficient to apply a homogeneous far field B ∞ in the chain axis direction with B ∞ = 1 T. For every type of helical chains, a parameter α = r c b , with b being the particle distance in the chain axis direction, controls the chain radius r c , i.e., the distance of the particle centers to the chain axis. This allows for investigating helical structures with different distance ratios, see Figure 8 for an illustration of the chain geometry, and represents a continuation of the studies performed in [56]. For the numerical simulations, a modulus E m = 200 kPa is assumed for the quasiincompressible matrix (ν m = 0.49 within the FE simulations of the micro-continuum approach), while the magnetization function of the particles is characterized by M s = 841 kA/m and ζ = 2.18 × 10 −2 m/kA according to Equation (13a). In order to quantify the magneto-deformation ε, the relative change of the structure length l is compared to the length l 0 in the absence of magnetic fields, i.e., ε = ∆l l 0 . The diameter and the distance of the particles in the chain axis direction are fixed, d = 5 µm and b = 6 µm, and all types of helical structures are analyzed for a number of n p ∈ {3, . . . , 17} particles to allow an investigation of up to 2 full chain cycles within the elastomer surrounding.
In Figure 9, the results of the comparison are summarized for all helical chains. Within the individual subplots, data obtained with micro-continuum simulations are indicated by markers, whereas the results of the particle-interaction model are indicated by solid lines of the same color. Note that all results represent bare findings of the comparison of the two modeling approaches, i.e., they are obtained by just applying the same magnetic field H ∞ and material parameters but entail no further fitting or adjustment.
Starting with the widely investigated plane chains with ∆ϑ = 180°in Figure 9a, three main conclusions can be drawn. First of all, the comparison of the micro-continuum and particle-interaction models shows a remarkable agreement over almost the whole range of analyzed distance ratios α: only for very small values of α, i.e., for the case of almost straight chains, a difference between the two approaches is visible. Keeping in mind that, with b/d = 1.2, adjacent particles are very close to each other for small values of α; these discrepancies can be ascribed to errors caused by the dipole approximation and might also be related to the assumption of a linear elastic matrix material within the particle-interaction model. A second finding is that both approaches consistently predict a contraction of wavy chains for almost all distance ratios α. Only in a very small range 0.4 ≤ α ≤ 0.65 a positive magneto-deformation is observed, see the magnified inset shown in Figure 9a. Finally, it is apparent that the effect of magneto-deformation decreases with the number of particles n p considered within the models. While the particle interactions in short structures yield a strong effect, especially for almost straight chains, this interaction is reduced systematically if n p is increased. This is a clear evidence for boundary effects which dominate the behavior of short chains and are of minor importance if larger, elongated structures are considered. Altogether, the results for the analyzed geometries are in line with other studies [56] and point out that the experimentally observed behavior of MAEs with chain-like microstructures cannot be explained with oversimplified wavy chains.
Since helical chains with ∆ϑ = 180°only represent a plane problem in which the complex interactions of particles in MAEs cannot be captured, the analysis is expanded to more elaborate helices with ∆ϑ = {90°, 60°, 45°}, as can be seen in Figure 9b-d. Again, the consistency of the results obtained with the two different modeling approaches is striking: within all examples, they only differ for very small values of α. Comparable to the preceding example, all helices show a negative magneto-deformation if the distance ratio is small, i.e., the chains are almost straight. Interestingly, the transition from contraction to elongation is shifted towards higher values of α, if ∆ϑ decreases. Additionally, it can be seen that-independent of the number of particles within the chains-their higher complexity allows for finding positive magneto-deformation within a broad range of α. In contrast to the helices with ∆ϑ = 180°, there is no second transition from elongation to contraction for larger values of α. A closer look at the magnified insets in Figure 9b-d reveals that the maximum elongations are almost independent of ∆ϑ. Again, boundary effects lead to a much stronger magneto-deformation, especially if the analyzed geometry does not represent at least a full chain cycle. However, even if up to two full cycles are considered, an elongation of the structure is still found. In summary, the present study shows that more complex helical structures represent an adequate generalization of the well-investigated wavy chains: their complexity allows for a broader range of possible effects and supports the idea that the particle arrangement in chain-like structures can trigger a transition from contraction to elongation in an applied magnetic field. In contrast to previous works, the results of this study show that such a transition cannot be realized with wavy chains: to attain a positive magneto-deformation with simplified chain-like structures, the particles must be arranged in a helical shape with ∆ϑ = 180°. As the primary focus of the analysis is a systematic comparison of the micro-continuum and particle-interaction models presented in Section 2, the remarkable agreement of the results obtained with both approaches must again be emphasized. This finding is in line with the outcome of a former comparison of the approaches for the simplified two-dimensional case [49] and promotes the idea of a hybrid multiscale framework for the analysis of realistic MAE samples in which computationally expensive FE simulations can be replaced by the less costly particle-interaction approach.

Conclusions
Within this work, two different modeling strategies for magneto-active elastomers are presented and compared with regard to their predictions for the magneto-deformation of chain-like helical structures in an elastomer surrounding. After both modeling strategies are introduced and their applications are illustrated by means of representative examples, the problem of particle interactions in elongated microstructures is analyzed. The investigation of helical chains with different angles ∆ϑ between adjacent particles represents a consistent continuation of former studies and allows for identifying effects of local particle rearrangements in MAEs on their macroscopic behavior. In line with a former comparison of the presented modeling approaches, the results show a remarkable agreement for all structures under investigation. Only in situations where basic assumptions of the particleinteraction model lose their validity, small deviations are found. Regarding the influence of different particle arrangements on the magneto-deformation of the structure, complex helical chains with ∆ϑ = 180°are found to allow for contraction as well as elongation depending on the chain geometry.
Certainly, the obtained results cannot fully explain the magneto-mechanical coupling behavior in MAEs with chain-like microstructures, since it also-and in many situations predominantly-depends on macroscopic shape effects as well as the elastic properties of the matrix material, among others. However, the findings of this study show that, at least on the microstructural level, the chain geometry can have a significant influence on the coupling behavior. To this end, the results should be a basis for further studies on MAEs with chain-like microstructures and, moreover, emphasize again that the complex interactions in those materials cannot be explained by just considering planar wavy chains.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Here, for the sake of generality, we account for the possibility of a position-dependent external field H ∞ . To simplify the derivation, we will make use of the form of Equation (1), resp. Equation (28), where the magnetization, i.e., the dipole moment, is aligned with the local magnetic field H i .
We start from Equation (27) and denote the position of particle k as r k . Hence, the magnetic force on particle k is given via the negative gradient with respect to r k : Note that the integral term in Equation (27) depends explicitly only on the local magnetic field H i up to which one has to integrate. Hence, we used the chain rule identity . Already at this step, one notes that the derivation is independent of the magnetization function M(H). From Equation (A1), follows: immediately follows.
Although we have m i H i , one must be aware that H ∞ i is not necessarily aligned with the dipole moments. Therefore, we change to tensor notation and rearrange Equation (A2) further to: Note the expressions of the form ∇ k H represent tensors of rank 2. We use the operator ∇ in standard notation as ∇ = e α ∂ ∂r α , where e α forms the basis vectors for example in Cartesian coordinates (α = x, y, z). In case of a constant external field, the last term in Equation (A3) vanishes, otherwise it only contributes if i = k and accordingly it represents the well-known force on a fixed magnetic dipole in varying external field (non-vanishing gradient of the field): Splitting the total magnetic force in Equation (A3) accordingly, we write: Here, f d k contains the first two terms of Equation (A3) representing the dipole-dipole interaction force. In analogy to Equation (30), we find in the dipole approximation the following relation: (A6) Parameter n i denotes the demagnetization factor of particle i. For a spherical inclusion, we have n i = 1/3 , ∀i, but, for generality, it is left as a parameter here. Note that in the case of ellipsoidal inclusions, the presented dipole approach is fully equivalent and works with the same accuracy as for spheres, only n i turns into the tensorial form to describe the demagnetization factor of an ellipsoid. The derivative of Equation (A6) with respect to the position of particle k gives: (A7) Here, the last factor on the r.h.s. represents a tensor of rank 3. Using Equation (A7) and (A6) in Equation (A3), the result for the dipole-dipole contribution f d k reads: (A8) Note that we used the notation " : " for the double dot product, also known as double inner product. Due to the symmetry i ↔ j and noticing that the application of the gradient ∇ k is only non-zero if either i = k or j = k, we further obtain: (A9) With r ki = r i − r k , and accordingly ∇ k = −∇ ki (denoting the derivative with respect to r ki ), the application of the nabla operator yields:  m i m k : 3r 2 ki r ki I + ∑ α e α r ki e α + Ir ki − 15r ki r ki r ki r 7 Finally, carrying out the double dot product, we find the well-known interaction force for mutually interacting dipoles: Thus, independent of the magnetization behavior, particle sizes, and particle demagnetization factors, the dipole-dipole interaction in Equation (A11) determines the force acting on an individual particle due to the presence of the surrounding inclusions in the first order dipole approximation as formulated in Equation (26) and (28). The form of Equation (A11) is rewritten in Equation (29) in the main text.

Appendix B
The starting point to set up a simulation model for particle rearrangement processes in MAEs is Newton's equation: The left-hand side represents the sum of all forces acting on some particle i at some time t. We denote m i as the particle mass andr i as the instantaneous acceleration. For monodisperse systems, clearly m i = m, ∀i. In the following, we discuss the meaning and the implementation of each term in the l.h.s. of Equation (A12) for the case of micron-sized magnetizable particles embedded in an incompressible linear elastic medium.
The magnetic force f d i acting on any particle i is again described via the dipole model Equation (29) where we assure self-consistency of the magnetic moments through Equations (28) and (30).
Contribution f el i shall provide the 'restoring' force due to the elastic matrix when the particles change their positions ({r 0 j } → {r j }). To implement f el i , we make use of the approximate results from the linear elasticity model as obtained by Puljiz et al. [30,88]. For an incompressible, isotropic matrix material, Equation (31) provides the superimposed relations between the set of locally applied forces {f j } and the set of displacements {∆r j } including hydrodynamic 2-and 3-body interactions. These linear relations must be apparently unique, i.e., the displacements ∆r j can only become zero for all j simultaneously if, and only if, no forces are applied at all, f j = 0 ∀j, as otherwise the elastostatic model is ill-defined. Thus, likewise, as a set of forces provides the corresponding set of displacements, also a given set of displacements {∆r j , j ∈ [1, N]} provides the corresponding set of forces that must be applied to each inclusion. Accordingly, since Equation (31) determines the 'new' elastically equilibrated positions with respect to the force free situation {r 0 j , j ∈ [1, N]}, we know that the elastic restore forces upon repositioning are determined as: Here, we defined the elasticity prefactor k := πEd p . The dimensionless tensors K ij are the entries i, j of that matrix of tensors found upon inversion of the matrix of tensors [D] ij = D ij . i.e., besides the minus sign in Equation (A13), it represents the inversion of Equation (31).
Since under general conditions the elastic restore forces f el i can not prevent neighboring particles from overlapping, or even merge into each other due to the singularity of the dipole-dipole force as r ij → 0 [73], the hard core repulsion f hc i prevents such artifactual events from happening. We will choose a centrosymmetric, short-ranged repulsive poten-tial of the type of a generalized Lennard-Jones potential. Such generalized potentials have been introduced long ago [95] and are known as Mie potentials: Originally, such potentials have been introduced to characterize inter-atomic, or inter-and intra-molecular, interactions, where r denotes the distance between the center positions of atoms, resp. molecules. Equation (A14) contains the well-known Lennard-Jones potential as special case ν = 12 and µ = 6. Typically, to describe purely repulsive behavior, i.e., athermal interactions, the Lennard-Jones potential is cut and shifted at its equilibrium, i.e., minimum, position r * [96]. Jover et.al. [97] considered a generalized cut-and-shifted Mie potential: λr −λa , (A15) and zero otherwise. Here, ∈ defines the energy scale. The authors in Ref. [97] suggest λ r = 50 and λ a = 49 to be used in molecular dynamics simulation when modeling hard sphere potentials as it proved to yield reasonable results and far more adequately resembles the ideal hard core potential. The form of U hc 50,49 (r) and the cut-and-shifted Lennard-Jones potential (U hc 12,6 (r)) are plotted in Figure A1 against the ideal hard core potential drawn as a gray dashed-dotted line. Unfortunately, this form U hc 50,49 (r) has two disadvantages when applying in our approach. Firstly, it contains an uneven power function, λ a = 49, which requires the computation of the square root of r 2 = r · r. Secondly, U hc 50,49 (r) raises very steeply in the onset of repulsive interaction (r r * ) requiring very small iteration steps to prevent 'breakdowns' or 'bursts' of the algorithm as two particles are approaching each other. Thus, for our purposes, the U hc 50,49 (r) would drastically slow down the computation performance and we suggest an alternative form U hc 180,2 (r) which copes with our requirements and at the same time it likewise resembles the ideal hard core repulsion as can be seen from Figure A1.   Figure A1. Plotting the generalized cut-and-shifted potentials from Equation (A15) for selected λ r and λ a . As a guide to the eye, we additionally sketch the ideal hard core potential as a grey dashed-dotted line. The classical Lennard-Jones potential (blue) can only very poorly resemble the ideal form-especially, since neighboring microspheres should only start to 'feel' each other (beside the hydrodynamic interaction included in f el in Equation (A13)), as they really come to physical contact at r ≈ d p , making alternative forms like U hc 50,49 (r) or U hc 180,2 (r) much more reasonable.
The corresponding force f hc i for each inclusion i is then constructed straightforwardly: The contribution f th i refers to thermal or Brownian motion of the particles and can be disregarded due to the comparatively large size of the inclusions d p ∼ 5 µm and the condition that they stick to an elastic matrix (i.e., no free motion). The energy ∆U el to deflect a particle from its equilibrium position by some ∆r can be estimated in leading order as ∆U el ∼ k∆r 2 . Equating this elastic energy with the typical order of thermal energy (∼k B T), we conclude that, due to thermal fluctuations, the corresponding deflections are of order ∆r ∼ (10 −4 − 10 −5 )d p considering typical values for the elastic modulus E ∼ (1-500) kPa. Thus, we neglect f th i in Equation (A12). The last term on the left side of Equation (A12), f fr i , refers to viscous friction when a particle is dragged through the embedding medium upon local repositioning. Primarily, this term defines the dynamical evolution of the rearrangement processes and accounts for the dissipation of kinetic energies. Its existence assures the particles to come to rest at some 'new' equilibrium positions after a typical time scale associated with the viscous properties of the system. Analogue to the linear elasticity model, we also consider the viscous properties of the matrix as linear, i.e., Newtonian fluid. Accordingly, in the leading order, we describe the viscous friction of a particle moving through the matrix via Stokes' law: Here, η denotes the viscosity parameter and we introduced the friction constant b := 3πηd p . Typically, for polymer melts, we have η ∼ 10 kPas. Since the chains are additionally crosslinked, we can safely assume that such η represents a lower limit for our system. We note that, already for η ∼ 10 kPas, the dynamics are highly overdamped. For example, the time scale ∆t until a particle reaches the limiting overdamped velocityṙ ∞ can be estimated as ∆t ∼ m b ∼ ρ Fe πd 2 p 18πη ∼ 10 −12 s, where we exemplary utilize the mass density of iron ρ Fe . Considering the magnetic dipole-dipole forces as typical 'drag' forces acting on the particles, we can evaluate the typical distance ∆r a particle moves within such time scale to reachṙ ∞ as ∆r 10 −7 d p . Alternatively, we can take recourse from well-known relations for an overdamped harmonic oscillator upon considering the particles embedded in the linear elastic medium as pinned to a Hookean spring with constant k. Then, the condition for overdamped oscillations reads b 2 > 4km which requires η 2 > 2 27 Ed 2 p ρ Fe . Apparently, for E ∼ 100 kPa and η ∼ 10 kPas, this holds by several orders of magnitude, 10 8 10 3 . Therefore, we may describe the motion of embedded particles as strongly overdamped and safely neglect the inertia term, i.e., the right-hand side, in Equation (A12). Finally, the model discussed here with its corresponding parameters, i.e., micron-sized inclusions in linear elastic medium, yields an equation of motion for each particle as dependent on the actual distribution of the other particles. Since in good approximation we can consider the motion as strongly overdamped, it represents a first order differential equation: We introduce the dimensionless time scale T := k b t = E 3η t and the dimensionless length scale R = 1 d p r. Thus, parameter η, which so far has only been specified by its rough order of magnitude η 10 kPas, is adsorbed in the timescale T and is no longer a relevant parameter for the simulation. It does not affect the rearrangement process when referring to dimensionless times T. With respect to the real time scale, i.e., t, the value of η rather universally tunes the speed of the rearrangement process. The remaining relevant forces, i.e., magnetic, elastic and hard-core repulsion, are summed in F i . To approximate the time evolution of the system, we apply a one-step discretization with fixed step-width ∆T = T n − T n−1 , ∀n and find the following iteration scheme: In the Supplementary Materials, we attached three example videos to visualize rearrangement processes as computed by the simulation approach suggested here.