Density functional approach to elastic properties of three-dimensional dipole-spring models for magnetic gels

Magnetic gels are composite materials consisting of a polymer matrix and embedded magnetic particles. Those are mechanically coupled to each other, giving rise to the magnetostrictive effects as well as to a controllable overall elasticity responsive to external magnetic fields. Due to their inherent composite and thereby multiscale nature, a theoretical framework bridging different levels of description is indispensable for understanding the magnetomechanical properties of magnetic gels. In this study, we extend a recently developed density functional approach from two spatial dimensions to more realistic three-dimensional systems. Along these lines, we connect a mesoscopic characterization resolving the discrete structure of the magnetic particles to macroscopic continuum parameters of magnetic gels. In particular, we incorporate the long-range nature of the magnetic dipole–dipole interaction and consider the approximate incompressibility of the embedding media and relative rotations with respect to an external magnetic field breaking rotational symmetry. We then probe the shape of the model system in its reference state, confirming the dependence of magnetostrictive effects on the configuration of the magnetic particles and on the shape of the considered sample. Moreover, calculating the elastic and rotational coefficients on the basis of our mesoscopic approach, we examine how the macroscopic types of behavior are related to the mesoscopic properties. Implications for real systems of random particle configurations are also discussed.


I. INTRODUCTION
Ferrogels, magnetic gels as well as magnetorheological gels and elastomers, all referred to as magnetic gels, are soft elastic composite materials containing magnetic or magnetizable particles, both simply referred to as magnetic particles.Their mechanical properties are controllable by external magnetic fields [1][2][3][4] .The composite nature arises as the magnetic particles are mechanically coupled to a surrounding polymeric matrix [5][6][7][8][9] .Such a magnetomechanical coupling has even been enhanced by anchoring polymers directly on the surface of magnetic particles 2,5,6,9 .To understand the rheological properties of these materials, the dependence of their elastic moduli and magnetostrictive effects on external magnetic fields have been investigated in various settings [10][11][12][13][14][15][16][17] .In particular, induced changes in the configuration of the magnetic particles, especially the touching of adjacent particles and chain formation, have been repeatedly reported as prominent features in the response of magnetic gels to external magnetic fields 7,[18][19][20][21][22][23] .
Due to their inherent composite nature, a complete theoretical understanding of magnetic gels is still challenging 24 .At macroscopic scales, thermodynamic and a) Electronic mail: s.goh@fz-juelich.deb) Electronic mail: a.menzel@ovgu.dec) Electronic mail: rene.wittmann@hhu.de d Electronic mail: hlowen@hhu.dehydrodynamic theories have been developed [25][26][27] , in which both the magnetic and elastic components are modeled as homogeneous continua.Notably, the positive magnetostriction, namely elongation along the magnetic field direction, has been predicted using such continuum approaches [28][29][30] .However, potential effects stemming from the detailed configurations of the magnetic particles and the polymers as well as the mechanisms governing these effects are hardly resolved at this scale.Rather than that, phenomenological coefficients have to be determined via modeling or experiments.In a theoretical perspective, one may consider a model resolving all the magnetic particles and polymer molecules at a microscopic level.Indeed, numerical simulation studies have been performed at this scale, revealing the roles of the polymer network topology and the coupling between the orientation of magnetic particles and the surrounding polymers 31,32 , as well as the degree of cross-linking in the polymer matrix 33 .However, unifying all the ingredients of such models to derive macroscopic parameters remains a demanding task.
In this regard, mesoscopic approaches still address the configurations of the magnetic particles explicitly, while individual polymeric building blocks are not resolved.Indeed, the significance of detailed structures at mesoscopic length scales has been revealed, as the mesoscopic models predict, for instance, both positive and negative magnetostrictive effects depending on the specific configuration of the magnetic particles [34][35][36][37] .The rotational fluctuations of magnetic particles have also been addressed arXiv:2211.01219v1[cond-mat.soft] 2 Nov 2022 within a mesoscopic approach 38 .Specifically, the polymer matrix can be coarse-grained as an elastic continuum [39][40][41][42] .We note that the role of the magnetic particles can also be modeled using continuum fields 43,44 that describe the particle arrangements.
Alternatively, the elastic continuum can be discretized on the mesoscopic scale as a network of harmonic springs [45][46][47] .One advantage of this approach is that microscopic theories as well as simulation techniques developed in the framework of statistical mechanics are directly applicable.The interaction energies are explicitly defined in this case.As demonstrated in Ref.48  using a description of a uniaxial magnetic gel, and in Refs.24 and 49 using an approach for isotropic one-and two-dimensional systems, a bridging description between mesoscopic and macroscopic scales may unravel the role of the discrete mesoscopic structures in the materials for the macroscopic behavior.In this way, the gap between continuum theories and mesoscopic models is closed.
In the present study, we further explore the issue of scale-bridging and, in particular, the statistical mechanics of magnetic gels.As for the mesoscopic description, we employ a simple but tangible model consisting of magnetic dipolar particles and harmonic springs connecting them.Starting from the mesoscale model, we aim at calculating macroscopic elastic and rotational coefficients, the trend of which we then correlate with mesoscopic characteristics, i.e., the configuration of the magnetic particles.Specifically, we employ classical density functional theory (DFT) [50][51][52][53] , extending the concept of mapping the elastic interactions between the particles through the surrounding elastic medium onto pseudosprings 24,49 to three dimensions.The resulting free energy allows for a calculation of macroscopic elasticity parameters.
To this end, the following issues need to be addressed in advance.First of all, the dipolar magnetic interaction is strictly long-ranged, rendering the system thermodynamically ill-defined 54 .While the Ewald summation technique 55,56 can be adopted to numerically simulate systems with long-range interactions such as suspensions of magnetic particles in liquid crystalline matrices 57,58 , the shape dependence of the free-energy has to be clarified as in the studies of dipolar fluids 59,60 and of magneto-sensitive elastomers 61 .In addition, we note that the aforementioned magnetostrictive effects originate from the anisotropic nature of the magnetic dipoledipole interaction.These points require a careful choice of the DFT implementation.Second, thermal fluctuations of the magnetic particles need to be included.As we develop a statistical theory for the equilibrium free energy, i.e., DFT, this issue is resolved automatically.Third, when an external magnetic field explicitly breaks the rotational symmetry of the system, relative rotations with respect to the field direction should be considered in the underlying elasticity theory.Originally, this concept was introduced in the context of liquid-crystalline elastomers [62][63][64][65] , but has also been extended to uniax-ial magnetic gels 26,48 .Lastly, the role of incompressibility that may be inherent in many systems of magnetic gels should also be clarified in the description.Just as conventional gels, magnetic gels can swell/shrink by absorbing/releasing liquid.Otherwise, they are regarded as incompressible, for instance, due to the dispersed fluid.Such incompressibility may alter the mechanical properties of magnetic gels 15,66 , calling for a theory respecting volume conservation.
This paper is organized as follows.In Sec.II our model for magnetic gels is introduced.We then formulate elasticity theory for incompressible systems in Sec.III, including the components of relative rotations.The DFT for our model system and its implementation are described in Sec.IV.In Sec.V, we present the elastic and rotation coefficients as well as magnetostrictive effects obtained from the DFT calculation.Lastly, discussions are included in Sec.VI.

II. MESOSCOPIC DIPOLE-SPRING MODEL
We consider a three-dimensional version of the previously studied dipole-spring system 46,67 as a mesoscopic model for magnetic gels.The model consists of N identical magnetic particles of diameter σ D and dipole moment m, which are connected by identical harmonic springs of spring constant k el and rest length a el .The position of the ith particle is denoted by r i (i = 1, . . ., N ).
The total Hamiltonian H tot of the system is introduced as the sum of the kinetic part H kin and the interaction part H int of the magnetic particles, the latter of which splits into three parts: Here, H m and H st , respectively, denote the energies of the magnetic dipole-dipole interaction and the steric repulsion, which are all-to-all pairwise additive.Therefore, with r ij = r j − r i , they take the form where u m and u st denote the two-body magnetic dipoledipole interaction and steric repulsion, respectively, as detailed below.
In stark contrast to Eq. ( 2), the elastic part H el does not simply take the form of a pairwise additive potential, namely no general two-body potential applying simultaneously to all pairs of particles can be introduced.Specifically, the elastic contribution is written in the form where i, j indicates that the sum only includes a predefined set of neighbors, which labels the particles such that they become distinguishable.Thus the potential energy cannot be written as a sum over pair potentials of indistinguishable particles.For the two-body potential u el a harmonic spring potential of spring constant k el is adopted, while a el is the rest length of the springs and r ij = |r ij |.Here, we assume a face-centered cubic (FCC) lattice structure with twelve nearest-neighbor particles, which is indicated by the angular bracket in Eq. ( 3).Therefore, in total 6N harmonic springs connect the nearest-neighboring pairs of magnetic particles in this specific model (except for boundary effects).
Next, for the two-body steric repulsion in Eq. ( 2), we assume a particle diameter σ D and adopt a hard-core potential in the form of As a dimensionless density we introduce the packing fraction η defined as the fraction of the volume occupied by the magnetic particles, i.e., η ≡ N (4π/3)(σ D /2) 3 /V where V is the volume of the system.Specifying the two-body magnetic dipole-dipole interaction energy, the two-body potential in Eq. ( 2) reads where µ 0 is the vacuum permeability.The magnetic moment m is determined by the applied magnetic field H = H ẑ, which is always directed along the z-direction, and the magnetization properties of each magnetic particle, see, e.g., Refs.68 and 69.When an external magnetic field is applied, we assume that m(H) H, i.e., m = mẑ, see Fig. 1(a) for illustration.In the absence of the applied field, magnetic particles may or may not retain their magnetization.Here we consider two simple cases.First, in Model I, we assume that the magnetic particles are ferromagnetic.There, the magnetic moment of each particle persists and is fixed with respect to the particle orientation, once the magnetic particles are magnetized.We then investigate elastic properties of the model system in the absence of external fields as depicted in Fig. 1(b).In this case, the dipole moment m rotates rigidly with the whole system.In Model II, we assume that the magnetic particles are paramagnetic.To retain the magnetization, the external field H needs to be persistently applied to the system in this case.In contrast to Model I, m is then always directed along H as shown in Fig. 1(c), even when the whole system rotates.Consequently, relative rotations between the magnetization direction and the rest of the system become relevant.We note that the magnetic dipole-dipole interaction breaks the isotropy of the system both in Model I and Model II.
As for the orientation of the system, we consider two cases in which the (0, 0, 1)-and (1, 1, 1)-orientations of the FCC lattice are directed along the z-axis 36 .When the lattices are elongated or contracted along the z-direction due to the anisotropic magnetic interaction, the resultant lattices of the (0, 0, 1)-and (1, 1, 1)-cases are tetragonal and rhombohedral, respectively.Lastly, we assume that our model system is incompressible, i.e., the volume of the whole system is fixed and persists even under deformations.Here we set V = ( √ 2/2)a 3 el N , at which the total Hamiltonian H tot is minimized for m = 0.

III. MACROSCOPIC DESCRIPTION
As incompressibility is assumed in our model system, we should address the maintained volume conservation when developing our macroscopic description.While our model system is initialized with the prescribed volume V at m = 0, an external magnetic field induces a magnetostrictive effect, which is not necessarily volume preserving.However, the imposed incompressibility constraint hinders the system to relax to a new volume upon magnetization.Such effects introduce a predeformation hidden behind the maintained volume, rendering our model system nonlinear elastic.Here, following the group theoretical approach proposed in Ref. 70, we consider nonlinear elastic responses of our model system and calculate elastic moduli accordingly.For self-containedness, we briefly summarize the formulation and introduce the second-order corrections to the deformation gradient tensors that are relevant for the computation of elastic moduli.

A. Nonlinear deformation gradient tensor
Under the incompressibility condition, the deformation gradients in three dimensions are elements of the special linear group SL(3, R), the Lie algebra of which is sl(3, R).Generally, the components of the deformation gradient tensor F are defined by F ij = ∂r i /∂r j , where r and r mark the positions of the material elements in the deformed and undeformed state, respectively.Then non-linear deformation gradient tensors F may be expressed via the exponential map where, λ i denote the SL(3, R) group generators and i are small coefficients indicating the magnitude of deformations generated by λ i .One should choose a set of generators, which is appropriate for the system considered.Accordingly, for our model system, we employ the generators of Here, the transformations associated with λ 1 stretch (compress) the system along the x-axis, combined with compressions (stretches) along the y-axis; the deformations generated by λ 2 stretch (compress) the system in the xy-plane combined with compressions (stretches) along the z-axis of twice the magnitude; λ 3 , λ 4 , and λ 5 generate shear deformations in the yz-, zy-, and xy-plane, respectively; λ 6 , λ 7 , and λ 8 generate rotations in the yz-, zy-, and xy-plane, respectively.We note that λ 5 can also be regarded as a generator of shear deformations in the xy-plane, but with orientations different from those generated by λ 1 .For the purpose of calculating elastic moduli, corrections up to the second order of i are relevant.Accordingly, we may truncate the expansion at the third order of { i } and use We note that, in general, generators do not commute, i.e., λ i • λ j = λ j • λ i .Within our approach, the free-energy density F (see Sec. IV for the definition based on density functional theory), equivalent to the deformation energy density in nonlinear elasticity, is regarded as a function of { i }.This choice naturally allows us to define the generalized elastic moduli as In the case of Model I, only the five generators λ i for i = 1, . . ., 5 are relevant, among which the shear deformations generated by λ 3 and λ 4 lead to identical contributions to F due to the symmetry of tetragonal and rhombohedral lattices.In addition to those, the relative rotations corresponding to λ 6 and λ 7 must be included for the description of Model II, whereas rotations in the xy-plane generated by λ 8 are still irrelevant.Again due to the symmetry, the rotations corresponding to λ 6 and λ 7 lead to identical contributions to F.

B. Irreducible representation for stiffness tensors
In the framework of linear elasticity theory, irreducible representations for stiffness tensors are determined by the underlying symmetry of the systems.We here consider the strains which are defined in linear elasticity as where u denotes the displacement field.For a tetragonal lattice [(0, 0, 1)-orientation], the stiffness matrix C in Mandel (or orthonormal) notation, where the stiffness matrix becomes a second-rank tensor 71,72 , takes the form We note that here the indices run from 0 to 5, not from 1 to 6.
We then turn to nonlinear elasticity.The infinitesimal group generators corresponding to Eq. ( 10) are { λi } for i = 0, . . ., 5, three of which are defined componentwise via [ λi−1 ] lm = δ il δ im for i, l, m = 1, 2, 3, and the others by λi = λ i for i = 3, 4, 5.However, such a choice of group generators is not compatible with the symmetry of systems under the incompressibility constraint.Therefore, we should introduce a transformation, which allows us to switch to the generators of Eq. (7).Since a set of infinitesimal group generators is a basis of a vector space, namely Lie algebra, we can find the form of generalized elastic constants corresponding to Eq. ( 7) via a linear transformation.Specifically, a unitary transformation connects, via = U • ˜ , the deformation vector ˜ = 1/ √ 2(˜ 0 , . . ., ˜ 5 ) in Mandel notation, corresponding to { λi }, to = ( 0 , . . ., 5 ), corresponding to the set of group generators consisting of {λ i } for i = 1, . . ., 5, and λ 0 ≡ 2/3I, where I is the 3 × 3 identity matrix, such that i i λ i = i ˜ i λi .Subsequently, the stiffness tensor C com for compressible systems can be computed from as is obvious from linear algebra.Within this representation, all deformations involving a volume change are associated with the generator λ 0 .Therefore, under the incompressibility constraint, the components in C com associated with λ 0 become irrelevant.Furthermore, predeformations also give rise to additional terms that are absent in linear elasticity, as demonstrated in Ref. 70.In our case, as only a predeformation in volume is involved, such nonlinear contributions are all diagonal and associated with the generalized pressure which only makes sense if a volume change is allowed.Finally, we conclude that the stiffness tensor of incompressible systems takes the form of whose components are defined by Eq. ( 9).Second, if the (1, 1, 1)-direction of the lattice is oriented along the z-axis, we have a rhombohedral lattice (RI Laue group), the stiffness tensor of which, again in Mandel notation, is given as 73,74 Then the stiffness tensor of the corresponding incompressible systems, within our notation, reads As we demonstrate in Sec.V, the macroscopic approach described here provides a precise and economic framework to investigate nonlinear elastic properties of incompressible anisotropic systems.In particular, our choice of generators given by Eq. ( 7) and the corresponding stiffness tensors given in Eqs. ( 14) and ( 16), respectively, determine all the necessary but only allowed deformations and elastic constants compatible with the underlying symmetry of the system and the imposed constraint.Sticking to the linear strain tensors as given by Eq. (10), instead of our nonlinear definition in Eq. ( 8), may involve errors in the second order, which are relevant for elastic constants.Indeed, the rotation coefficients C 66 and C 77 shown in Fig. 4(a) can become negative, if the volume conservation in the second order is not explicitly taken into account via Eqs.( 6) and (7).Alternatively, one may consider the method of Lagrange multipliers, which is, however, technically demanding, particularly in combination with the density functional calculation that we describe next.

IV. DENSITY FUNCTIONAL THEORY: BRIDGING SCALES
We now formulate a density functional theory (DFT) for the dipole-spring model by approximating the free energy functional F[ρ(r)] where ρ(r) denotes the onebody density field of the magnetic particles.Together with the ideal gas term where β ≡ (k B T ) −1 is the inverse temperature, the total free-energy functional subjected to minimization is given as where F ex [ρ(r)] denotes the excess functional describing the interparticle interactions (1).Following Ref. 75, we employ the Picard iteration algorithm with a mixing parameter α and where which is updated in each iteration step to ensure that the total (average) number of particles is kept fixed.Accordingly, for the verification of successful minimization, we use the relative chemical potential defined as In this way, F is minimized for a prescribed value of the vacancy concentration n vac .In principle, our model systems are defect-free, i.e., n vac = 0. To accelerate and enhance the robustness of the minimization processes, however, we consider lattices with vacancy concentration of n vac = 0.001 ± 10 −6 .
Regarding the geometry of the calculation box, we use the primitive unit cell in our calculations, consisting of only one particle, instead of the cubic unit cell of the FCC lattice, consisting of four particles, usually adopted in DFT studies of freezing.Accordingly, we use periodic boundary conditions in the directions of three primitive vectors.Both the primitive and reciprocal lattice vectors of undeformed and deformed systems are summarized in Appendix A. With this geometry, we are able to minimize our free-energy functional more precisely, (∆µ rel < 10 −15 in most cases), compared to the method using the cubic unit cell (∆µ rel ≈ 10 −8 for the tested cases).Now we turn to the excess functional F ex [ρ(r)], which is given as a sum of three functionals corresponding to the steric repulsion, magnetic dipole-dipole interaction, and harmonic spring potential.First, for the hard-core repulsion, we use the White-Bear II (WB-II) functional 76 with the Tarazona tensors 77 , which is one of the most precise versions among the fundamental measure theory for hard spheres 78 .Then, for the elastic and magnetic dipole-dipole interactions, we intend to adopt the simple mean-field functional in the form of where u(r) is an appropriate pair potential.However, the practical evaluation of the above functionals is not straightforward.In what follows, we describe how to construct the Fourier transforms of the elastic and magnetic energies, which allow us to perform DFT calculations in Fourier space.

A. Magnetic dipolar interaction
As discussed, there are two important properties inherent in the magnetic dipole-dipole interaction, Eq. ( 5), in three dimensions.It is long-range and anisotropic 79 , which has to be taken into account in DFT calculations.
When we switch m to m = 0, the systems elongate or contract, and so does the unit cell.Then the side lengths of the cubic unit cell are no longer the same, but satisfy the relation a x = a y = a z where a x , a y , and a z denote the side lengths in the x-, y-, and z-direction, respectively.Here, we define the aspect ratio as R asp ≡ a z /a x .We note that R asp characterizes the deformation of the internal lattice structure.Now, we address the long-range nature of the magnetic dipole-dipole interaction in three dimensions.The difficulty arises from the fact that the interaction energy, i.e., the integral of u m , diverges at both short and long distances.In our DFT calculation, this issue can be resolved rather easily.On the one hand, the steric repulsion hinders particles from approaching closer than their diameter and therefore prevents the divergence at short distances.On the other hand, as the DFT calculation is performed in Fourier space, the divergence at long distance can be handled directly as follows.While, for k = 0 Fourier modes, the Fourier transform of the dipole-dipole interaction can be obtained with the aid of the orthogonality of the spherical harmonics Y m l , the k = 0 mode, which dictates the long-range divergence, indeed depends on the shape of the whole material body (see Appendix B for more details).With such a shape dependent mode, which is related to the demagnetizing factor in continuum theory 61 , we are able to capture the long-ranged nature of the magnetic interaction.In general, we may consider a system with the initially spheroidal shape (at m = 0) of the shape parameter R sh ≡ R z /R x , where R x = R y and R z are the lengths of the semiaxes along the x-, y-and z-axis, respectively.In contrast to R asp , here R sh indicates the aspect ratio of the whole material.As we turn on the magnetic interaction applying a magnetic field, the initial aspect ratio of the whole system shape further changes to R asp R sh due to magnetostriction associated with a change in internal lattice structure.

B. Elastic energy
While the magnetic particles are strictly labeled due to fixation by the surrounding polymer matrix, namely the elastic interaction term given as Eq. ( 3), the conventional machinery of DFT calculation assumes the indistiguishability of particles, i.e., as if the potential u el was acting equally between all pairs of particles throughout the system.To nevertheless enable DFT calculations, a mapping of the harmonic spring potential onto a pseudospring potential u pseudo has been proposed in Ref. 24.There, only nearest-neighbor pairs of the resulting configuration are within the range of u pseudo and thus elastically coupled to each other, as in the original system based on the harmonic springs, see Eq. ( 3).In the present study, the mapping is extended to three dimensions.Notably, in two and three spatial dimensions, the success of applying the finite-ranged pseudospring potential is connected to the particle arrangement arising from a freezing transition 49 , which has been extensively investigated within density functional approaches [80][81][82][83][84] .Specifically, we consider the pseudospring potential In this expression, R c and u 0 denote the cut-off length and the offset for the pseudospring potential, respectively.The cutoff length R c is determined from corresponding Monte-Carlo simulations as the distance at which the pair correlation function g(r) is minimized, which turns out to be R c = 1.21a el .The value of u 0 is determined within the DFT calculations so as to match the vacancy concentration of the resultant lattice with the prescribed value of n vac = 0.001.The obtained u pseudo instead of u el is inserted into Eq.( 23) via u.We refer to our previous study 49 for the detailed description and verification of the mapping of the real onto the pseudospring potential.Moreover, due to the anisotropy of dipolar interactions as discussed in Sec.IV A, corresponding lattice structures may become anisotropic as well.Such anisotropy then should also be taken into account when we construct the pseudospring potential.In practice, we then cut the spring potential at the surface of the spheroid with the aspect ratio R asp , instead of at the surface of the sphere with the radius R c as in Eq. (24).In other words, the cutting is direction-dependent.The resultant Fourier components are presented in Appendix C explicitly.

V. MECHANICAL PROPERTIES
From now on, we measure lengths and energies in units of the rest length a el and the thermal energy k B T , respectively.Accordingly, the magnitude m of the magnetic moment and the spring constant k el are measured in units of m 0 ≡ k B T a 3 el /µ 0 and k B T /a 2 el , respectively.We consider systems with elastic constant k el = 100 and shape parameter R sh = 1, and investigate the effects of magnetization on the mechanical properties, by varying the magnitude of the magnetic moment m.The two models, described in Sec.II, give identical results as long as no rotations are considered, while only the paramagnetic Model II has a unique reference state with respect to rotations.
One can also probe steric effects by varying the volume packing fraction η.Naively speaking, while steric repulsion should affect the bulk modulus of a system, how and to what extent it would affect the mechanical properties under each specific deformation is still unclear.Moreover, there might also appear numerical artifacts due to several approximations employed.Therefore, leaving systematic investigations for further studies, we demonstrate that our method is valid for a reasonable range of η by employing two representative values of η = 0.1 and 0.3, which are relatively small when compared to the coexisting fluid (crystal) packing fractions 0.495 (0.544) 75 for the WB-II functional used in this study.We note that the pseudospring potential suffices to stabilize the FCC crystal for η = 0 within our model.Indeed, steric repulsion does not play a dominant role for these low packing fractions, as one may confirm from Fig. 2 in Sec.V A, as well as from Figs. 3 and 4 in Sec.V B.

A. Magnetostriction
As a first step, we determine the reference equilibrium state of the undeformed system for a given magnetization m.Technically, we first determine the value of u 0 for which the vacancy concentration of system becomes equal to the prescribed value within the margin of tolerated error.Then, varying R asp while fixing u 0 , we find the value of the aspect ratio R asp at which the free energy functional is minimized.The resultant values of R asp are shown in Fig. 2.
The most prominent feature here is that the magnetostriction effects of the (0, 0, 1)-and (1, 1, 1)-orientations are opposite to each other.In line with the results reported in Ref. 36, our systems elongate when the dipole moments are directed along the (0, 0, 1)-orientation, while a contraction along the direction of the dipole moments is observed in the (1, 1, 1)-case, confirming that the internal configuration of magnetic particles is a decisive factor of the magnetostriction effect.In addition, we also note that the magnetostriction effect can be reversed, if large values of the shape parameter (R sh 2) are used in the case of the (0, 0, 1)-orientation (results not shown).Such shape-dependence is a trivial consequence of the long-range nature of the dipolar interaction.

B. Elastic coefficients
Now we determine the elastic constants C ij for i, j = 1, . . ., 5, defined in Eq. ( 9), from our DFT, explicitly deforming the primitive unit cell.Specifically, we numerically calculate the derivatives through finite differences and obtain the diagonal terms of the stiffness tensor from while the offdiagonal terms can be calculated as In most of the cases, we use 1 = 2 / √ 3 = 3 = • • • = 7 = 0.0001, except for the cases of the (1, 1, 1)orientation with η = 0.3 and m ≤ 2.5, in which the functional can be minimized up to the values of ∆µ rel between 10 −4 and 10 −8 at most.There, we use 1 = .0001 or 0.001 to obtain con-sistent results.Before proceeding to the results in Fig. 3, we recall from Sec. III that some coefficients vanish and others are equal.Specifically, we confirm in Figs.3(b) and (f) for rhombohedral lattices that C 11 = C 55 and C 13 = C 45 , respectively, in accordance with Eq. ( 16).
First, we take a closer look at the elastic constant C 55 , corresponding to shear deformations in the xy-plane, and C 11 , corresponding to stretches (compressions) along the x-axis combined with compressions (stretches) along the y-axis.C 11 can also be regarded as a shear modulus, but corresponding to shear deformations with orientations different from those for C 55 .In most cases, the dipolar interaction, which is repulsive in the plane perpendicular to the dipole moment, causes an increase of the elastic constants.Specifically, as shown in Fig. 3(b), values of both C 11 = C 55 increase as m increases in the (1, 1, 1)-case, while, in the (0, 0, 1)-case, only C 55 is an increasing function of m, as shown in Fig. 3(e).In sharp contrast, C 11 in the (0, 0, 1)-case is a decreasing function of m, as shown in Fig. 3(a).Furthermore, as m increases further, it drops towards zero, indicating instability of the tetragonal lattices.We notice here that hexagonal configurations can be obtained eventually by squeezing the tetragonal lattice in the xy-plane, if the whole lattice is projected on the xy-plane.In other words, as m increases, there might arise a growing tendency to match the lattice to the underlying symmetry of the magnetic dipole-dipole interaction, which prefers the hexagonal lattice over the tetragonal lattice in the plane perpendicular to the dipole moment.Therefore, we conclude that such a softening effect correlates with a rearrangement of the magnetic particles in the plane perpendicular to m.
Next we turn to the elastic constants of C 22 , corresponding to stretches (compressions) in the xy-plane combined with compressions (stretches) along the z-axis of twice the magnitude, and C 33 , corresponding to shear deformation in the yz-plane (or equivalently C 44 , corresponding to shear deformations in the xz-plane).All of them involve deformations in the z-direction.In both the (0, 0, 1)-and (1, 1, 1)-orientations, C 22 is an increasing function of m [Fig.3(c)], indicating hardening of the materials.Since there is no significant difference between the systems of η = 0.1 and 0.3, the phenomenon of hardening observed here has a purely elastic origin.Simultaneously, C 33 is always a decreasing function of m [Fig.3(d)].Moreover, at large m, the rhombohedral lattice becomes unstable as well, with the values of C 33 dropping towards zero.Such instabilities at large m and the decrease of C 33 in general may originate from the tendency towards pair formation 20,85 or similarly from the typical chain-like aggregates forming under strong dipolar interactions 48,86,87 .Indeed, we observe a shift of the energetic minimum in the landscape of two-body interaction energy from separated to touching configurations occurs between m = 2.5 and 3.0 in the (1, 1, 1)-case (not shown).This seems to confirm that the instability is the consequence of the formation of touching pairs.In the case of the (0, 0, 1)-orientation, the drop towards zero in C 11 occurs in advance of that in C 33 , compare Figs.3(a) and (d), indicating that rearrangement in the xy-plane is preferred over rearrangement in the z-direction.
Lastly, the values of C 13 and C 45 in the (1, 1, 1)-case are presented in Fig. 3  Rotation coefficients obtained from the (1, 1, 1)orientation.In (a), the coefficients C66 and C77, respectively, corresponding solely to the rotations in the yz-and xz-planes are presented, whereas the off-diagonal coefficients C36 and C47 as well as C16 and C57 are depicted in (b) and (c), respectively.
behavior, increasing from negative values for small m to positive ones for large m.We note, however, that these constants reflect a rather specific symmetry inherent in the lattice, and therefore, may not reflect the situation of real magnetic gels.

C. Rotation coefficients
Finally, we investigate the rotation coefficients, which are relevant only in Model II.Alike the elastic constants, the rotation coefficients are calculated from Eqs. ( 25) and (26).As the (0, 0, 1)-orientation turns out to be unstable with respect to rotations in xz-and yz-planes, we only analyze the results for the (1, 1, 1)-orientation.
First, the coupling of the model systems to the applied magnetic field is captured by the rotation coefficients C 66 and C 77 , corresponding to rotations in the yz-and xzplane, respectively.As shown in Fig. 4(a), the values of C 66 and C 77 increase as m increases, indicating an enhanced resistance to the rotations.
As shown in Fig. 4(b), the mixed coefficients of C 36 and C 47 , corresponding to mixed shear deformations and rotations in the yz-and xz-plane, respectively, first exhibit an increase as a function of m for small values of m.Then, the increasing trend is reversed for large m.We note that, in Ref. 48, where chain-like aggregates are assumed, only a decreasing tendency in the form of −m 2 has been predicted for D 2 , which is equivalent to C 36 and C 47 in the present study.Presumably, as already mentioned for magnetostrictive effects in Sec.V A, different behaviors may be due to the internal configuration of the magnetic particles.We also note that the values of C 66 and C 77 are approximately 10 3 times smaller than those of C 33 and C 44 .In Ref. 48, the rotation coefficient D 1 (equivalent to C 66 and C 77 in the present study) is even larger than ∆c 5 (C 33 and C 44 in the present study).Again, this may be caused by the different internal structure, which is chain-like in Ref. 48.
Lastly, the additional mixed coefficients C 16 and C 57 increase monotonically, as shown in Fig. 4(c), seems to be a simple consequence of enhancement of both hardening in the xy-plane (C 11 and C 55 ) and resistance to rotations in the xz-or yz-plane (C 66 and C 77 ).

VI. CONCLUSION
So far we have constructed and evaluated a DFT for three-dimensional dipole-spring models, which bridges from the discretized mesoscopic model to a macroscopic elasticity theory of magnetic gels.Based on the scalebridging description, we have determined the elastic and rotational material coefficients.They depend on the mesoscopic configuration of the magnetic particles.Notably, we have observed softening responses to magneti-zation both in the external field direction and in the plane perpendicular to the external field, which indicates a tendency towards an instability.We have proposed that such behaviors imply changes in overall symmetry, accompanied by rearrangement of magnetic particles.Such rearrangements might be decomposed into the formation of a hexagonal-like arrangement in the plane perpendicular to the magnetic field and pair formation along the magnetic field direction.To verify our conclusion, the decreasing behavior of C 11 needs to be tested experimentally.Notably, in a previous study, where random configurations for magnetic particles are assumed [see Fig. 14(b) of Ref. 46], a decrease of the shear modulus has been observed, suggesting that the idea of rearrangement may also be valid for real magnetic gels with disordered configurations.
Conversely, one could equally well think about synthesizing a sample with the regular arrangement adopted in this study.In particular, the prescribed FCC-based connectivity shows certain characteristics as explained above.For instance, 6 among 12 nearest neighboring particles are located in the same plane perpendicular to the magnetization in the case of (1, 1, 1)-orientation, and thereby, the repulsive interaction in the plane seems to dominate the response of the magnetic particles.This leads to the contraction along the magnetization direction.We note that there have been attempts to synthesize thin ferrogel films 88 .Since in planar configurations, magnetic particles form hexagonal arrangements in the plane perpendicular to the external magnetic field 89 , it would be possible to obtain ferrogel films with a hexagonal configuration in such a way.Then, by stacking two-dimensional layers, a magnetic gel with a threedimensional hexagonal structure might be fabricated experimentally.Our results of the (1, 1, 1)-case may then provide an insight into such systems.
At the same time, regarding future work on our theory, one important direction is to address systems with random configurations.An important additional ingredient to model the heterogeneity inherent in real samples is polydispersity of the magnetic particles 90 .The idea of the replica DFT 91,92 might be used to address directly disordered configurations.Lastly, dynamical density functional theory [93][94][95] should provide a route to investigate the dynamics of the systems.
The reciprocal lattice vectors for the (0, 0, 1)-orientation.As for the infinitesimal parameters { i}, see Eqs. ( 6) and (7) in which the deformation gradient tensor as well as the generators are defined.

Appendix A: Reciprocal lattices
For the (0, 0, 1)-, and (1, 1, 1)-orientations of the FCC lattice, the primitive vectors read and respectively.Here, a = √ 2a el denotes the side length of the cubic unit cell.In practice, the DFT calculations are performed with the reciprocal lattice vectors in Fourier space.For the (0, 0, 1)-orientation, the reciprocal vectors read while for the (1, 1, 1)-orientation, we obtain Under deformation, the reciprocal vectors are transformed accordingly.We expand the reciprocal vectors of deformed lattices with respect to { i } to compute the corresponding reciprocal lattice vectors in the form The correction terms ∆b 1 , ∆b 2 , and ∆b 3 for the (0, 0, 1)-and (1, 1, 1)-cases are given in Tables I and II, respectively, which are sufficient for the pure deformations that do not involve mixed terms, i.e., i j for i = j.When more than two different types of deformations are applied, Eq. ( 6) still provides a correct formulation.However, such mixed terms are irrelevant for our incompressible systems because second-order corrections only enter via the diagonal terms in the stiffness tensor, as we describe in Sec.III B (see Ref. 70 for details).Therefore, for the calculation of off-diagonal components in the stiffness tensors, we simply add the second-order corrections from two different types of pure deformations.
where j l and Y m l are spherical Bessel functions and spherical harmonics, respectively, and the superscript asterisk * denotes complex conjugate.Since the dipole-dipole interaction energy [Eq.(5)] is proportional to Y 0 2 for m = mẑ, i.e., where γ = R asp R sh .Apparently, the k = 0 Fourier mode depends on the shape of systems, namely the aspect ratio γ.

Deformed systems
When k = 0, the Fourier transformation is shape independent.For the k = 0 mode, however, the Fourier transform of the deformed system is in general different from the undeformed one, due to the dependence on the sample shape.To calculate the correction, we first clarify how a deformation F modifies the integration via Ω(r,θ,φ) where the prime indicates that the region of integration has been changed according to the deformation.Then, we recover the original shape of the system by changing the variables via r = F • r where the center dot • denotes matrix multiplication, and subsequently, rewriting the integration as while the boundaries of integration region stemming from the hard-core repulsion must be modified accordingly.We note that the differential d 3 r = dr r 2 sin θ remains unchanged because |det F| = 1.Also note that F ≡ dr /dr.Then the above integration can be performed up to the second order of { i } with straightforward algebra, which has been performed using Mathematica 96 .Here, with U 0 (θ) = − sin θ − 3 sin 3θ 16 ln (cos 2 θ + γ 2 sin 2 θ), (B7) we write the integrand in Eq. (B6) as where ∆U is the correction due to deformation.First, for stretches (compressions) along the x-axis, combined with compressions (stretches) along the y-axis associated with λ 1 , ∆U reads ∆U (θ, 1 ) = − ≡ ∆U 1 (θ, 1 ).(B9) At the same time, the correction stemming from the deformation associates with λ 2 are already reflected in Eq. (B4), as we have calculated the values of u m (k = 0) for arbitrary aspect ratios.Because of the uniaxial symmetry of the magnetic dipolar interaction, the correction due to the shear deformations in the xy-plane takes the same form as Eq.(B9), namely, ∆U (θ, 5 ) = ∆U 1 (θ, 5 ).Next, the correction due to the shear deformations in the yz-plane is given as ∆U (θ, 3 ) =

FIG. 1 .
FIG. 1.(a) Under an external magnetic field H, the magnetic dipole moment m (black solid arrow) is aligned along the external field direction H (red dashed arrow).Here the whole system (solid ellipse) has the same orientation (green dashed line at the center of the ellipse) as the initial magnetization, which is the reference state of our model system.Regarding rotations, we consider two models, (b) Model I for ferromagnetic particles and (c) Model II for paramagnetic particles, see text for details.

3 FIG. 2 .
FIG.2.The aspect ratio Rasp of the system is presented as a function of m.Converse magnetostriction effects are observed, depending on the orientation of the lattice.

C 57 FIG. 4 .
FIG. 4.Rotation coefficients obtained from the (1, 1, 1)orientation.In (a), the coefficients C66 and C77, respectively, corresponding solely to the rotations in the yz-and xz-planes are presented, whereas the off-diagonal coefficients C36 and C47 as well as C16 and C57 are depicted in (b) and (c), respectively.