Towards a soft magnetoelastic twist actuator

Soft actuators allow to transform external stimuli to mechanical deformations. Because of their deformational response to external magnetic fields, magnetic gels and elastomers represent ideal candidates for such tasks. Mostly, linear magnetostrictive deformations, that is, elongations or contractions along straight axes are discussed in this context. In contrast to that, we here suggest to realize a twist actuator that responds by torsional deformations around the axis of the applied magnetic field. For this purpose, we theoretically investigate the overall mechanical response of a basic model system containing discrete magnetizable particles in a soft elastic matrix. Two different types of discrete particle arrangements are used as starting conditions in the nonmagnetized state. These contain globally twisted anisotropic particle arrangements on the one hand, and groups of discrete helical-like particle structures positioned side by side on the other hand. Besides the resulting twist upon magnetization, we also evaluate different other modes of deformation. Our analysis supports the construction of magnetically orientable and actuatable torsional mixing devices in fluidic applications or other types of soft actuators that initiate relative rotations between different components.

Soft actuators allow to transform external stimuli to mechanical deformations. Because of their deformational response to external magnetic fields, magnetic gels and elastomers represent ideal candidates for such tasks. Mostly, linear magnetostrictive deformations, that is, elongations or contractions along straight axes are discussed in this context. In contrast to that, we here suggest to realize a twist actuator that responds by torsional deformations around the axis of the applied magnetic field. For this purpose, we theoretically investigate the overall mechanical response of a basic model system containing discrete magnetizable particles in a soft elastic matrix. Two different types of discrete particle arrangements are used as starting conditions in the nonmagnetized state. These contain globally twisted anisotropic particle arrangements on the one hand, and groups of discrete helical-like particle structures positioned side by side on the other hand. Besides the resulting twist upon magnetization, we also evaluate different other modes of deformation. Our analysis supports the construction of magnetically orientable and actuatable torsional mixing devices in fluidic applications or other types of soft actuators that initiate relative rotations between different components.

I. INTRODUCTION
Torsional actuators respond by a twist-type deformation to external stimuli. Most studies are concerned with linear actuators that contract or elongate along a certain axis upon actuation. However, there are several important prospective applications of twist actuators, for example microfluidic mixing, microscopic surgery tools, and prosthetics [1]. Depending on the application, a certain degree of softness of the actuator in combination with a certain degree of biocompatibility may be beneficial or even mandatory, particularly when it comes to medical applications. This is one of the reasons why socalled magnetic gels and elastomers (also commonly referred to as magnetorheological elastomers or ferrogels) [2][3][4][5][6][7][8][9][10][11] were introduced as important candidates for the construction of soft actuators [2,[12][13][14][15][16][17][18][19][20][21][22][23][24]. These materials usually consist of magnetic or magnetizable colloidal particles embedded in an elastic, typically polymeric matrix. Such magnetic gels have the advantage that their distortions can be induced by external magnetic fields and the resulting deformation is typically reversible [25].
To now generate magnetoelastic twist actuators in the form of magnetic gels or elastomers, see Fig. 1, we suggest to build on the following previously explored insights. When the materials are fabricated in the presence of strong homogeneous external magnetic fields, chainlike structures of the inserted particles may form before the surrounding polymeric matrix is permanently established through corresponding chemical processes. Once the elastic matrix has reached its elastic solid state, these particle structures remain locked in the material, as can * lfischer@thphy.uni-duesseldorf.de † menzel@thphy.uni-duesseldorf.de be seen in many experimental realizations [13,[26][27][28][29][30][31][32][33].
One possible route to generate torsional actuators may be to additionally twist these chain-like aggregates, before the particle positions are fixed in the material by the final chemical crosslinking and establishing of the elastic polymeric matrix. This leads to self-supported torsional actuators. Such a concept is different from materials that are clamped at one end, contain anisotropic nonchiral structures, and are twisted by external magnetic fields that exert torques on the contained anisotropic aggregates [34,35]. Naturally, the situation that we consider is also different from studying how magnetic fields modify the stiffness of magnetic gels and elastomers when distorted by externally imposed torsional deformations [36][37][38][39][40][41][42].
To realize soft torsional actuators, in our case, on the one hand, one may think of a globally, collectively twisted state of the whole set of embedded chain-like aggregates in the initial, cured state of the materials. On the other hand, one may consider each individual chain-like aggregate to show an initially twisted structure.
We start by considering globally twisted particle arrangements as initial states. To generate corresponding samples, a procedure of the following kind could be realistic. The approach is inspired by a protocol of synthesizing monodomain nematic liquid-crystalline elastomers [43][44][45], consisting of liquid-crystal molecules that are chemically attached to or part of crosslinked polymeric networks [45][46][47]. Its scheme follows a two-step crosslinking process [43][44][45], employing two crosslinkers of different speed of chemical reaction. The action of the first crosslinker generates a weakly crosslinked elastomeric sample that is stiff enough to already be uniaxially stretched. Maintaining this stretched state, in which the liquid-crystal molecules are on average uniaxially oriented in response to the imposed strain, the second crosslinker reacts and locks in this configuration. Along these lines, monodomain nematic samples, featuring an average uniaxial molecular liquid-crystalline alignment, are obtained. Such materials show pronounced nonlinear stress-strain properties when stretched perpendicular to the direction of nematic alignment [43,44,48,49].
In our case of magnetic gels and elastomers, the two-step crosslinking process may be performed accordingly. First, under the presence of strong homogeneous external magnetic fields, uniaxially ordered chainlike aggregates of the magnetized particles form. They get locked into the sample by the generated surrounding elastic environment resulting from the quick action of the first crosslinker [26-28, 31-33, 50, 51]. Finite gaps between the particles as considered below may result from previous coating of the particles, or by using surface-functionalized particles themselves as crosslinkers [19,[52][53][54][55]. In a next step, this pre-crosslinked system is twisted around the anisotropy axis. This leads to a global twist of the contained chain-like particle aggregates. The sample is maintained in this state while the second, slower crosslinker is reacting chemically and establishing the final elastic matrix. In this way, the twisted structure gets permanently locked in.
Another, possibly more academic procedure to generate example systems for investigations of the effects that we here predict might be to deposit the particles in a controlled way, maybe even by hand, at prescribed positions while generating the elastic environment layer by layer [56,57]. Even macroscopic spherical particles could be used for such proofs of concept [58]. In this case, besides implementing globally twisted structures, one could also arrange the magnetizable particles in individual helices, positioned in an aligned way side by side. Maybe, in the future, such a deposition process can be automatized, as has recently been achieved for the production of magnetic microhelices [59,60].
In the present work, we use such twisted discrete particle configurations as an input to calculate resulting magnetically induced overall deformations of corresponding elastic composite systems. Our theoretical approach is analytical, based on linear elasticity theory, and then evaluated numerically. To achieve such an analytical approach, we concentrate on elastic systems of overall spherical shape. The degree of initial twist is varied and the consequences of such variations are analyzed, both by numerical evaluations and by simplified analytical considerations. Both the globally twisted structures as well as several individual helical-like aggregates arranged side by side are addressed.
We give a brief overview on our theoretical approach in Sec. II, together with a motivation of our chosen parameter values. After that, in Sec. III, the torsional actuation of systems containing globally twisted particle configurations are addressed. In Sec. IV, we consider particle arrangements of helical aggregates positioned side by side. To further facilitate the understanding, we compare the resulting twisting deformation to a minimal analytical consideration in Sec. V. We conclude in Sec. VI.
FIG. 1. Illustration of the general idea and setup. The considered soft magnetoelastic composite system is spherical in overall shape. Upon application of a homogeneous external magnetic field B, it shows a reversible torsional twist deformation as indicated by the curved arrows on the right-hand side.

II. THEORETICAL FRAMEWORK
To perform the following evaluations, we build on our methods developed in Ref. 61. We assume that the elastic material used for the magnetorheological elastomer and containing the magnetizable particles is spatially isotropic as well as homogeneous. Moreover, we confine ourselves to small deformations (up to a couple of percent) so that we can use linear elasticity theory. This allows us to superimpose the deformations resulting from each internal force center.
Consequently, we describe this material via only two elastic coefficients, namely the shear modulus µ and the Poisson ratio ν. They quantify the stiffness and compressibility of the material, respectively. A Poisson ratio ν of 1/2, representing an upper bound [62], describes incompressible materials. However, the Poisson ratio can reach negative values as well, down to −1 [62]. In these cases, the corresponding material is called auxetic, implying that when stretched along one axis it will show expansion to the lateral directions instead of contraction.
Generally, the response of the elastic material to an applied force density f (r) inside it is then quantified by the so-called Navier-Cauchy equations [63], where u(r) denotes the displacement field at position r. In our case, the elastic material forms a free-standing elastic sphere of radius R. Fortunately, an analytical solution for Eq. (1) in this case is available in terms of the corresponding Green's function. f (r) then specifies the effect of a point-like force center inside the elastic sphere. We were able to transfer this solution to the case of a free-standing sphere of free surface [61], starting from previous work that considered the sphere embedded in an elastic background material [64]. This analytical solution for the elastic part of the problem was afterwards implemented numerically. Next, to include the magnetic effects of magnetorheological gels and elastomers, we distributed magnetic in-clusions at prescribed positions inside the elastic material, see Sec. III and Sec. IV. We always assume the magnetic inclusions to be sufficiently far apart from each other so that we can describe their magnetic signature as magnetic dipoles. As a further simplification, we assume that the magnetic dipolar moment m = mm, where m = |m|, is identical for all inclusions. In experiments, such a situation could be realized by applying a strong external magnetic field that magnetizes all (identical) inclusions to saturation.
In this case, the magnetic dipole-dipole forces are given by [65] where F i is the force exerted by all other inclusions on the ith inclusion. Moreover, µ 0 denotes the magnetic vacuum permeability,r i marks the position of the ith inclusion, the difference vector of positions is given byr ij =r i − r j =r ijrij withr ij = |r ij | (i, j = 1, ..., N ), and N sets the number of magnetized inclusions. The resulting force density inserted into Eq. (1) based on Eq. (2) is where δ(r) represents the Dirac delta function and we thus assume point-like magnetic force centers.
After rescaling lengths by R and forces by µR 2 , the strength of the magnetic forces relative to the elastic restoring forces is characterized by a dimensionless force coefficient 3µ 0 m 2 /4πµR 6 . Its value is set to 5.4×10 −8 for all that follows, as inspired by realistic experimental parameters [61]. The inclusions are assumed to be of spherical shape as well, with their radius set to a = 0.02R.
To include the effect of the induced elastic distortions on the positions of the magnetized inclusions and thus on the resulting magnetic forces and vice versa, an iterative scheme had been developed, see Ref. 61. Finally, to characterize the induced overall deformations and capabilities of actuation, we evaluate the resulting displacement field on 49152 surface points of the elastic sphere. These points are approximately evenly distributed with positions generated by the HEALPix package (http://healpix.sourceforge.net) [66].
For the problem at hand, we choose the z-axis to always coincide with the magnetization direction of the magnetic inclusions, i.e.m =ẑ. Moreover, we express the displacement of each surface point using spherical coordinates as Thus, the components u ⊥ , u θ , and u ϕ describe displacements inwards or outwards of the elastic surface, tangential deformations along the polar direction and tangential deformations along the azimuthal direction, respectively. Below, the latter coefficient u ϕ will become particularly important to quantify the overall twisting deformation.
To associate the resulting displacement field with different modes of overall deformation, we perform spherical harmonic expansions of u ⊥ , u θ , and u ϕ . We use the same definitions for spherical harmonics, especially concerning the Condon-Shortley phase, as in Ref. 65. The most relevant spherical harmonics for our analysis are given by Y 00 = 1/4π, Y 10 = 3/4π cos θ, and Y 20 = 5/16π 3 cos 2 θ − 1 .
As announced above, we then focus on the resulting overall torsional deformations for two types of spatial arrangements of the magnetizable inclusions: globally twisted and side-by-side aligned helical structures, see Secs. III and IV, respectively. The degree of initial structural twist in the nonmagnetized state is quantified by a parameter γ, see below for its definition. In both cases, we confine the initial positions of the inclusions by requiring a minimal distance of 3a = 0.06R to the elastic spherical surface.

III. GLOBALLY TWISTED STRUCTURES
To numerically generate the globally twisted structures, we start from layers of hexagonally arranged magnetic inclusions [51,[67][68][69]. These layers are all oriented parallel to the xy-plane and spaced equally from each other in their normal direction by a distance d layer = 0.11R. The center layer is located in the plane z = 0. In the initial, nonmagnetized situation, the hexagonal particle arrangements within each layer are in a state rotated by an angle of γz/d layer relative to the arrangement in the plane z = 0. This corresponds to a globally twisted configuration of the inclusions when compared to straight chain-like aggregates aligned parallel to the Illustration of two layers of hexagonally arranged magnetizable inclusions inside the elastic material. d layer sets the vertical distance between two layers, d chain the in-plane particle distance. We set d chain = 0.25R > d layer = 0.11R, which implies vertically aligned chain-like aggregates. Here, for illustration, d layer is exaggerated. The upper arrangement shows a rotation by an angle γ relative to the lower arrangement, where we chose γ = π/6 for reasons of visibility. To emphasize the twist from layer to layer, we plot the positions corresponding to the lower layer in the upper layer as gray spheres, together with a dotted arrow that shows their vertical identification. Having applied a rotation by γ to the structure from the lower layer (gray), the positions marked by dark spheres result. We indicate this in-plane rotational displacement by blue in-plane arrows. In the teal triangle, we illustrate the definition of the angle γ.
z-axis. Here, we consider small angles γ 0.159π to preserve the chain-like structure, see Fig. 2. The lattice constant within each plane, which equals the lateral distance between the chains, is set to d chain = 0.25R. Overall, this leads to 623 magnetizable inclusions in 55 chains. An illustration of an initial structure is presented in Fig. 3, where we have, however, increased d chain to 0.5R for clarity. In the numerical evaluation, we consider the range 0 γ 0.159π in steps of approximately 0.0016π. We distinguish four possible values of the Poisson ratio: ν = 0.5 (incompressible), ν = 0.3, ν = 0, and ν = −0.5 (auxetic).
As a first step, we focus on the following spherical harmonic expansion parameters for the resulting overall surface distortions: u ⊥ 00 , u ⊥ 20 , and u ϕ 10 . The coefficient u ⊥ 00 quantifies overall changes in volume of the composite material. Positive values correspond to an increase in volume, while negative values correspond to a decrease in volume. Next, the coefficient u ⊥ 20 describes a relative elongation (u ⊥ 20 > 0) or contraction (u ⊥ 20 < 0) along the direction of magnetization, here along the z-axis. Most important for our investigation in the present context is the parameter u ϕ 10 . This coefficient is set by the lowest mode of a twist-type deformation around the z-axis. For a counter-clockwise rotation of the upper hemisphere against the lower hemisphere, it becomes u ϕ 10 > 0. For a 3. Illustration of an example for the initial structure of the magnetizable inclusions, indicated as small spheres, inside the larger elastic sphere. This structure is generated from hexagonally arranged parallel chain-like aggregates of particles, where each horizontal layer of particles is rotated relative to the next particle layer underneath by an angle γ, see Fig. 2.
In this illustration, we chose γ ≈ 0.019π. Moreover, for better visibility, we here set d chain = 0.5R. Instead, for our actual numerical evaluation, we used a value of d chain = 0.25R.
reversed mutual sense of rotation, one obtains u ϕ 10 < 0. The three coefficients u ⊥ 00 , u ⊥ 20 , and u ϕ 10 are shown in Fig. 4 when the aforementioned particle structures are magnetized. We have not included the data for negative values of γ because the curves in Figs. 4(a) and 4(b) are mirror symmetric with respect to the line γ = 0, while the curve in Fig. 4(c) features a point symmetry with respect to the origin.
As a first result, we infer from Fig. 4(a) that the overall volume is constant (u ⊥ 00 ≈ 0) for ν = 0.5, as expected for an incompressible material. With decreasing Poisson ratio, we find that the elastic sphere shrinks more and more upon magnetization. Naturally, this volume decrease is maximal for γ = 0, i.e. straight chains of magnetizable inclusions. In this case, the induced attraction between the particles along each chain is strongest. When increasing γ, the volume decrease becomes smaller and oscillates for higher values of γ.
Similarly, we infer from Fig. 4(b) that the overall contraction along the magnetization direction relative to a lateral expansion, as quantified by u ⊥ 20 , is strongest for γ = 0 for the same reason as above. This effect is most pronounced for incompressible materials because the contraction along the field implies lateral expansions for reasons of volume conservation. In contrast to that, the auxetic nature for ν = −0.5 counteracts the lateral ex- To quantify perpendicular surface displacements, we plot the two coefficients (a) u ⊥ 00 and (b) u ⊥ 20 , indicating overall volume changes and overall elongation along the magnetization direction relative to lateral contraction, respectively. To quantify the lowest mode of an overall twist deformation around the magnetization axis, we plot the coefficient (c) u ϕ 10 . In all three cases, we display the behavior with increasing angle γ, characterizing the global twist of the initial nonmagnetized structure of inclusions (see Fig. 2 for the definition of γ). Moreover, we show graphs for the four different values of the Poisson ratio, namely ν = 0.5, ν = 0.3, ν = 0, and ν = −0.5. pansion for γ = 0. The oscillations for increasing values of γ can be found in this coefficient as well.
When we focus on the behavior of u ϕ 10 in Fig. 4(c), we observe that it is almost independent of the Poisson ratio. This is expected because a pure twist-type deformation leaves the total volume unchanged. A small effect of the Poisson ratio is still present and can most likely be attributed to nonlinear effects revealed by our iterative scheme, i.e. to the effects of the resulting displacements of the magnetic inclusions which are larger for more compressible materials. Furthermore, we do not observe any torsional deformation for γ = 0 because our initial config-uration is not twisted in this case. Increasing γ from zero, we see that the corresponding values of u ϕ 10 first become more and more negative. The sign here represents the sense of the induced torsional deformation of the composite which is opposing the sense of initial twist of the initial structure. We reach a maximum magnitude of this twist deformation at γ ≈ 0.019π. For larger values of γ, the magnitude of u ϕ 10 again decreases. This effect results from the increasing distance between the inclusions with increasing γ implying a decreasing magnetic interaction. At even larger values of γ, u ϕ 10 oscillates around zero. We return to this feature in Sec. V.
In practice, one would typically be interested in the situation of maximum observed effect. We therefore concentrate on the system for γ ≈ 0.019π. First, we checked how the magnitude of the induced torsion around the zaxis varies with the height z above or below the horizontal center plane (the xy-plane). For this purpose, we calculated the average azimuthal angular displacement of the horizontal plane parallel to the xy-plane at height z as where . . . z denotes an average over all surface points at which u ϕ was evaluated at a given height z. We found that this quantity is approximately proportional to z. Furthermore, we find that it is nearly independent of the Poisson ratio, in agreement with the behavior of u ϕ 10 in Fig. 4(c).
Next, in Fig. 5, we provide additional information on the importance of different modes involved in the overall surface displacement, obtained by our expansion of the perpendicular and tangential components of the surface displacement field into spherical harmonics. Again, we concentrate on the value of γ ≈ 0.019π and we use the same four values of the Poisson ratio as in Fig. 4. We select the expansion coefficients a lm of ten representative spherical harmonic modes for each component of the displacement field according to the following scheme. First, for each mode the value of a lm of highest magnitude is identified from the four values for the different Poisson ratios ν. These a lm are then ordered according to their absolute values, and we find the labels l, m for the ten largest ones. Due to the high degree of symmetry of our configurations, the most dominant modes are those of m = 0. However, we observe nonvanishing modes that depend on ϕ as well, characterized by m = 0. This leads to complex expansion coefficients. Since u ⊥ , u θ , and u ϕ are real, we can find values for negative m via the relation a l(−m) = (−1) m a * lm , where the star denotes complex conjugation. Consequently, for real a lm the corresponding spherical harmonics result together with a l(−m) in a cos (mϕ) mode, while purely imaginary a lm result in a −sin (mϕ) mode. The real and the imaginary part of a lm are shown separately in the plots. Figure 5 confirms that those coefficients that we have been concentrating on so far indeed dominate the spec- trum. For u ⊥ , see Fig. 5(a), these correspond to an overall volume change (l = m = 0), especially for auxetic materials and except for ν = 0.5, and to an overall contraction along the magnetization direction relative to a lateral expansion (l = 2, m = 0), with small higherorder corrections. All coefficients odd in l for m = 0 are approximately zero here. We observe some very small contributions related to the six-fold symmetry about the z-axis in the modes of l = 15, m = 6 and l = 24, m = 6. Turning to u θ in Fig. 5(b), significantly smaller absolute values of the expansion coefficients are obtained. Here, as for u ϕ in Fig. 5(c), the coefficients even in l van-ish approximately for m = 0, in contrast to the case for u ⊥ . The most important contribution to u θ in the mode l = 1, m = 0 corresponds to an overall surface displacement towards the equator on both the upper hemisphere and the lower hemisphere upon magnetization. In the incompressible case, this effect is most pronounced as we then have the strongest expansion of the sphere in the lateral directions. Again, higher-order contributions emerge which strengthen the aforementioned effect in the vicinity of the equatorial plane.
Considering u ϕ in Fig. 5(c) reveals the most important mode in the present context, associated with the twist deformation through magnetization. As noted already above, the mode l = 1, m = 0 is associated with a rotation around the magnetization direction of the upper hemisphere relative to the lower hemisphere. This mode dominates the overall behavior by its absolute value [only exceeded by the mode corresponding to overall volume expansion for the auxetic case ν = −0.5 in Fig. 5(a)]. Near the equatorial plane, higher-order modes in combination still support the effect of the upper hemisphere being rotated relatively to the lower hemisphere.

IV. HELICAL STRUCTURES
As a next step, we address helical structures of the magnetizable particles embedded in the same elastic spheres as before, arranged side by side. In contrast to the globally twisted structure of parallel chain-like aggregates investigated in Sec. III, we now consider each chainlike element by itself to feature an initial helical shape. To set up our numerical systems, we again start from hexagonal arrangements of aligned chain-like aggregates as before, this time for d chain = 0.5R, i.e. for double the distance to each other. As above, the vertical distance of the horizontal layers of particles is set to d layer = 0.11R. However, instead of initiating each layer rigidly rotated relatively to its upper and lower neighboring one, we now rigidly displace each layer laterally by adding a vector This lateral shift introduces an additional parameter, namely r helix . Here, we show results for structures corresponding to two different values r helix = 0.05R and r helix = 0.1R, see Figs. 6 and 7, respectively. In both cases, we fit 95 magnetizable inclusions into our elastic sphere, avoiding inclusions that would need to be deleted for particular values of γ. Importantly, the overall structure in each case is no longer six-fold rotationally symmetric about the z-axis nor globally screw-symmetric within the spherical boundaries. In the center layer for z = 0, all helices start with a particle deflection in the xdirection, r helix (0) = r helixx , according to Eq. . Furthermore, we here chose the radius of each helix to be r helix = 0.05R. In the depicted case, we set γ = π/8. depicted in Figs. 6 and 7. Using our numerical approach, we evaluate the full range 0 ≤ γ ≤ 2π in steps of π/360. As in Sec. III, we first address the expansion coefficients u ⊥ 00 , u ⊥ 20 , and u ϕ 10 for the overall displacements upon magnetization as functions of γ. The curves in Figs. 8(a), 8(b), 9(a), and 9(b) show a mirror symmetry with respect to the vertical line γ = π, while those in Figs. 8(c) and 9(c) feature a point symmetry with respect to the point (γ, u ϕ 10 ) = (π, 0). This is expected because Y 00 and Y 20 are even in z, while Y 10 is odd. Obviously, the results for helices of an initial twist π < γ < 2π can be mapped onto those for a corresponding initial twist of 2π − γ. Illustratively, this corresponds to helices that only differ by their sense of twist. The resulting displacements are of much smaller magnitude when compared to the results for the globally twisted arrangements in Fig. 4, which can already be expected from the lower total number of inclusions for the helical structures (95 here versus 623 inclusions in Fig. 4).
We start by considering the configurations of r helix = 0.05R. In Fig. 8(a), we again find that the elastic sphere, except for ν = 0.5, shrinks as a whole, specifically for the smallest and largest values of γ. For these values, the chains are straightest and therefore show the maximal internal longitudinal attractive forces. Furthermore, the absolute magnitude of overall contraction strongly increases with decreasing Poisson ratio, i.e. for more com-   pressible spheres. Next, we address in Fig. 8(b) the elongation along the magnetization relative to the lateral contraction. Qualitatively, we infer a similar behavior as in Fig. 8(a). Here, we observe a further, much smaller minimum for γ = π because we have effectively generated two chains of distance 2r helix = 0.1R out of each helix. Apart from that, auxetic materials show stronger relative contractions along the magnetization axis.
Concerning the magnitude of the twist actuation quantified by Fig. 8(c), we again find a pronounced minimum, here around γ ≈ 0.24π. In line with the point symmetry of the curve mentioned above, the corresponding maximum is located at γ ≈ 1.76π. As in Sec. III, u ϕ 10 as a measure for the twist actuation is approximately independent of the Poisson ratio. This behavior will also be discussed in Sec. V. Figure 9 shows corresponding results for r helix = 0.1R. The qualitative picture is similar to Fig. 8, with the same symmetries of the curves. We notice that the aforementioned minimum at γ = π is more pronounced and can be found in u ⊥ 00 [ Fig. 9(a)] as well. Concerning the coefficient u ϕ 10 quantifying the twist actuation, we see that the minimum is shifted to smaller values of γ, namely to γ ≈ 0.13π, and is increased in magnitude by a factor 0.0 of approximately 2.23. Moreover, some oscillations together with positive values of u ϕ 10 occur at higher values of γ < π. Again, we will return to this topic in Sec. V.
We continue with a discussion on the coefficients obtained from an expansion into spherical harmonics at that value of γ representing the minima in the curves of Figs. 8(c) and 9(c). Corresponding values are displayed in Figs. 10 and 11 for the configurations of r helix = 0.05R and r helix = 0.1R, respectively. Again, we plot ten relevant modes for each of the three components of the surface displacement field identified according to the same scheme as in Sec. III. Figure 10(a) shows that u ⊥ 00 and u ⊥ 20 dominate the overall behavior (for ν = 0.5 we correctly find u ⊥ 00 ≈ 0). Some higher-order contributions to u ⊥ are observed, however, of a relative magnitude of less than 15 % of the dominant mode, given by either u ⊥ 00 or u ⊥ 20 . The configurations are less symmetric than those in Sec. III, and we observe a stronger influence of the modes of m = 0, particularly for m = 1, which characterizes the lowest-order nontrivial dependence on ϕ.  Fig. 10, but for a setup with r helix = 0.1R and for γ ≈ 0.13π. The latter value identifies the minimum on the curve of Fig. 9(c).
Next, Fig. 10(b) identifies u θ 10 as a dominating mode of u θ for ν ≥ 0. The same was observed in Fig. 5(b). In general, higher-order modes enter as well, especially for auxetic materials. As before, the maximal magnitude of the modes described by u θ is smaller than the magnitude of the dominating mode for u ⊥ .
The modes relevant to torsional deformations of the elastic material are addressed in Fig. 10(c). We observe again the most dominant mode to be the lowest one, i.e. u ϕ 10 . However, we also find another mode to be almost equally as strong, namely u ϕ 41 . This is most likely an effect related to the specific helical structure that was used in our investigation. Nevertheless, both modes are of smaller magnitude when compared to the modes of u θ and even smaller when compared to u ⊥ . Thus, the twisting actuation for this structure is only of secondary importance when compared, for instance, to the global volume change or relative elongation along the magnetization direction.
In addressing the results for the structures of r helix = 0.1R in Fig. 11, we mainly focus on the differences when compared to the situation in Fig. 10. Due to the larger magnitude of r helix , the asymmetry of the configurations with respect to rotations around the z-axis by π/3 is still more pronounced and we thus observe even stronger modes for m = 0. This trend concerns all three components of the surface displacement field in Figs. 11(a), 11(b), and 11(c). Differences between Figs. 10 and 11, especially in the modes for m = 0, can to some extent be traced back to the different value of γ of the investigated structure, according to the different locations of the minima in Figs. 8(c) and 9(c).
Particularly when focusing on the torsional deformation addressed in Fig. 11(c), we observe that the mode u ϕ 10 identifying a global twist deformation is not even the strongest one here. Instead, the strongest mode is u ϕ 31 . This mode is symmetric for z → −z, implying that it cannot describe an overall twist deformation corresponding to a relative rotation of the top hemisphere with respect to the bottom hemisphere. However, the lowest mode of twist actuation u ϕ 10 has a much higher absolute magnitude when compared to the structures of r helix = 0.05R in Fig. 10(c). Apparently, the radius of the helical elements can have a pronounced effect, partly of antagonistic consequences. If such systems are transferred to actual applications, it is therefore important to adjust the radius of the helical elements to the desired behavior.
Overall, we observe a significantly more pronounced influence of higher-order modes and particularly modes depending on ϕ for the displacement fields in Figs. 10 and 11 when compared to the results in Fig. 5. Importantly, the ratio of the magnitudes of u ϕ 10 to the magnitudes of u ⊥ 00 (except for ν = 0.5) or u ⊥ 20 is much smaller. Thus, the twist actuation is not as pure for the investigated structures composed of helical elements and we conclude that the globally twisted structures of Sec. III are in general more promising candidates to construct a magnetoelastic twist actuator.

V. MINIMAL ANALYTICAL MODEL
Having presented our numerical results for the functions u ϕ 10 (γ) in Secs. III and IV, shown in Fig. 4(c) as well as in Figs. 8(c) and 9(c), respectively, we here discuss how we can understand the behavior qualitatively in simpler terms. To this end, we propose a minimal analytic model based on the dipole-dipole force between the inclusions, see Eq. (2). If we only concentrate on the magnetic interactions between two nearest neighbors on a single chain, the geometry can be parameterized as depicted in Fig. 12.
Obviously, the situation in reality is more complex as magnetic dipole interactions are long-ranged, leading to magnetic interactions between all particles. Furthermore, due to the magnetically induced deformations the particle positions are affected as well, which changes in turn the magnetic interactions, see Sec. II. Nevertheless, considering pairwise nearest-neighbor interactions along one chain will allow for a basic qualitative description,  12. In a simplified discussion, we consider the interactions between the magnetized nearest-neighboring particles i and j on an initially twisted chain-like aggregate. Their dipole moments, aligned with the center axisẑ, are depicted by small arrows. We denote the vector from the position of j to the position of i byrij. In (a), their distance along the z-axis is given by d layer and their lateral distance is denoted as ∆. ζ quantifies the angle between the z-axis andrij at the site of particle i. In (b), we show a bottom-view of the configuration. We introduce two right-angled triangles to relate the lateral distance ∆ between the particles to the radial distance ρ of the particles from the center axis around which the initial twist of the structure was set. The angle γ was defined previously for both the globally twisted structures and the helically twisted structural elements arranged side by side, see Fig. 2 and Eq. (7), respectively. see below.
Since we are interested in the magnetically induced overall twist deformation, we here focus on the magnetic force components perpendicular to the magnetization direction, i.e. in the xy-plane. These in-plane force components are the source of torsional deformations around the z-axis. Instead, the z-components of the magnetic forces are associated with axial contractions. For initially twisted particle configurations and not too large values of γ, the in-plane force components represent restoring forces that aim to straighten the chains. Defining ζ as the angle betweenm and the connecting vectorr ij between two nearest-neighboring particles i and j, see Fig. 12(a), the magnitude F xy of the in-plane force component on particle i, exerted by particle j, see Eq. (2), is given by Here, we have insertedm ·r ij = cos ζ,r ij = d layer / cos ζ, and sin ζ for the component ofr ij in the xy-plane. Next, we maximize F xy (ζ) with respect to ζ to find out which configuration of particles i and j leads to a maximized restoring force, which supports a maximized twist actuation. The maximum is found for If we now restrict the solutions to the range 0 < ζ < π/2, we find the unique solution When we compare to our previous results, we can use the relations deduced from Fig. 12(a) tan ζ = ∆ d layer (11) and Fig. 12 where we have introduced ρ as the distance of the inclusions i and j from the axis of the initial twist of the corresponding structure. To relate the result of this analytical consideration to our numerical evaluation, we find from Eqs. (11) and (12) γ max = 2 arcsin d layer 2ρ tan ζ max (13) where γ max implies a maximized in-plane torsional force component, based on this simplified analytical consideration. For the systems addressed in Sec. IV, to compare these analytical and the numerical results, we can simply set ρ = r helix . For the globally twisted configurations in Sec. III, the situation is more complex as there is not a single value of ρ that is equal for all chain-like aggregates, but the value of ρ depends on which chain we consider.
To illustrate this more complex dependence for the globally twisted structures on the angle γ, quantified by u ϕ 10 (γ), we have generated additional globally twisted configurations while removing from the systems considered in Sec. III those chain-like elements that have a value of ρ smaller than a certain threshold. Illustratively, this corresponds to only considering those chains that are located outside a coaxial cylinder of diameter 2ρ. In Fig. 13, we present results for cut-off values of ρ of R/2, 2R/3, and √ 13 d chain ≈ 0.901R, where the latter value marks the outermost chains. For comparison, we have added in Fig. 13 the results for the configurations of Sec. III as well. For this evaluation, we restrict ourselves to incompressible elastic materials (ν = 0.5) for clarity.
The main result of Fig. 13 is that as we increase the lower threshold value of ρ, the global minimum is shifted towards lower values of γ. For all chains considered, see Sec. III, the value of γ corresponding to a maximized twist deformation is γ ≈ 0.019π. Introducing a cut-off for ρ of R/2, this value is reduced to γ ≈ 0.018π. Moving on to a cut-off for ρ of 2R/3, it is further reduced to γ ≈ 0.016π. When keeping only the outermost chains, we obtain γ ≈ 0.014π for the location of the maximized twist deformation. Moreover, we observe a decrease in magnitude of the minimum of u ϕ 10 . This contains, however, a trivial effect as we decrease the number of inclusions for  Fig. 4(c), but for configurations for which we only consider those chains that have a minimal distance ρ from the axis of twist of the initial nonmagnetized structure. We show a comparison between the results of Fig. 4(c), here labeled as "all chains", and corresponding configurations that only include those chains for which ρ > R/2 and ρ > 2R/3. Furthermore, we show results for only keeping the outermost chains, i.e. chains of ρ = √ 13 d chain ≈ 0.901R. In all cases, we only display the results for incompressible systems, i.e. for ν = 0.5, for clarity. Particularly, we note how the position of the global minimum is slightly shifted towards smaller values of γ for configurations of larger average values of ρ for the chains. The vertical gray dashed line marks the value γmax ≈ 0.015π as obtained from Eqs. (10) and (13) for the outermost chains.
increasing cut-off values for ρ. More precisely, we find 623, 324, 168, and 60 inclusions for the four different systems addressed in Fig. 13.
When we now compare our numerical results to the minimal analytical model according to Eqs. (10) and (13), we consider the configurations of only keeping the outermost chains. In this case, inserting ρ ≈ 0.901R into Eq. (13), we obtain a value of γ max ≈ 0.015π, see the vertical dashed line in Fig. 13. This is only slightly bigger than the numerical value of γ ≈ 0.014π. It shows a fair agreement, considering for instance the assumptions of including only nearest-neighbor particle interactions and rigid particle positions in the analytical model. Next, we compare the numerical and analytical results for the structures composed of helical elements as considered in Sec. IV. Setting ρ = r helix , we find from the analytical consideration γ max ≈ 0.28π and γ max ≈ 0.14π for r helix = 0.05R and r helix = 0.1R, respectively. The results of our numerical investigation for u ϕ 10 were γ ≈ 0.24π and γ ≈ 0.13π, respectively, see Sec. IV. While showing fair agreement concerning the involved approximations, our analytical model again shows a tendency of overestimating the numerical results, see above.
Within our minimal analytical model, we may equally well estimate analytically the lowest value of γ > 0 for which u ϕ 10 becomes zero. Again, we require a fixed value of ρ. From Eq. (8), we find that F xy = 0 for a value ζ 0 > 0 of corresponding to Inserting the value of ρ for the outermost chains in Fig. 13 implies γ 0 ≈ 0.078π, while the numerical result for u ϕ 10 suggests γ ≈ 0.057π. As before, we observe that the analytically determined value of γ exceeds the one determined numerically. For the systems containing the helical structural elements, our analytical estimate does not imply any value of 0 < γ < π at which F xy = 0 because d layer > r helix in both cases. This is in line with our numerical results for r helix = 0.05R, for which u ϕ 10 < 0 for all 0 < γ < π. However, our numerical investigation reveals a value of γ ≈ 0.62π at which u ϕ 10 = 0 for r helix = 0.1R.
As had become obvious above and from Fig. 13, comparing the simple analytical model approach to the numerical results for the complete globally twisted structures is less direct. The different chain-like aggregates in the system are located at different radial distances ρ from the center axis. These varying distances need to be taken into account.
To find a reasonable measure, we start from the magnetic forces F i according to Eq. (2) on each particle i. We denote byφ i the local azimuthal unit vector in the spherical coordinate system at the position of particle i. To identify those force components that supposedly directly support the overall twist deformation, we project F i ontô ϕ i on the upper hemisphere and onto −φ i on the lower hemisphere. Particles i located on the center axis and on the equatorial plane are not taken into account. Finally, the sum over all force components obtained in this way is denoted as Σ xy . It is plotted as the dotted line in Fig. 14, together with the results for u ϕ 10 as displayed in Fig. 4(c). Since the initial positions of the particles are used for this basic analytical evaluation, the curves for Σ xy are independent of the value of the Poisson ratio ν.
Comparing these graphs, we notice that the force component Σ xy as well as the curves for u ϕ 10 have a pronounced minimum at approximately the same value of γ ≈ 0.019π. Moreover, for the lowest value of γ > 0 at which Σ xy = 0 and u ϕ 10 = 0, we find γ ≈ 0.068π and γ ≈ 0.072π, respectively.
In summary, we can estimate certain characteristic points on the curves of u ϕ 10 (γ) by simple analytical model considerations. Often, it is sufficient to focus on the interactions between neighboring dipoles only. For the globally twisted structures, see Sec. III, the different distances of the chain-like aggregates from the center axis of the elastic sphere need to be taken into account for more quantitative evaluations. Nonaffine elastic deformations have not been included in the simple analytical model. Furthermore, our simplified analytical approach does not account for the change in magnetized particle positions during deformations included in our numerical description.

VI. CONCLUSIONS
To conclude, we have suggested a way to construct soft torsional actuators using magnetic gels and elastomers. For this purpose, we have addressed two different structural arrangements of the magnetizable inclusions in the elastic material: globally twisted structures and side-byside arrangements of helical elements. Both are generated from initially hexagonally arranged parallel chainlike elements. For both configurations, we have explicitly calculated the resulting magnetostrictive distortion of the overall system upon magnetization. In this context, for reasons of analytical accessibility, we have here concentrated on systems of overall spherical shape. Particularly, we have focused on the degree of induced twist actuation, which we quantified using a spherical harmonic mode of expansion of the surface displacement field.
Among the systems that we investigated, we found the globally twisted structures to show a significantly larger twist actuation when compared to the systems containing helical elements arranged side by side. For the studied globally twisted structures, the overall deformational response is indeed dominated by a twist-type distortion. Instead, the overall twist response in the case of the embedded helical elements arranged side by side was less pure. Thus, it appears that the considered globally twisted structures are better suited to construct a soft torsional actuator. In the near future, these might also be the ones requiring less additional effort for actual fabrication.
Furthermore, we have quantified which degree of initial structural twist in the nonmagnetized state leads to a maximized twist actuation. Such an optimized value arises from two antagonistic tendencies. If the internal structure is not twisted at all, then an overall torsional deformation cannot be induced. However, if the initial structure is twisted too much, then the interactions between the inclusions upon magnetization even become repulsive for a too large lateral separation of the contained magnetizable particles. We find an optimized value in between. In fact, these properties can be understood already on a qualitative basis by addressing the magnetic interactions between two neighboring particles on one initially deformed chain-like structural aggregate.
Overall, we hope that our study will inspire experimental realizations in the future. Corresponding devices may find possible applications, for instance, as microfluidic mixing actuators. Not only can twist deformations and thus torsional flows around such an element be induced upon request from outside by alternating magnetic fields. But also, due to the existence of an overall structural anisotropy identified, for example, by the axis of global twist, can the mixing element simultaneously be oriented by the direction of the external field. Moreover, as long as dynamic effects like leaking electrical currents do not play an important role, our results equally apply to the construction of corresponding devices from electrorheological gels and elastomers [14,70,71], using external electric fields for actuation.