[orcid= 0000-0003-2888-5878] [1]

[cor1]Corresponding author

# Magnetostriction in the magneto-sensitive elastomers with inhomogeneously
magnetized particles: pairwise interaction approximation

D. Yaremchuk<EMAIL_ADDRESS>Institute for Condensed Matter Physics of
the National Academy of Sciences of Ukraine, 1, Svientsitskii Str., 79011,
Lviv, Ukraine  D. Ivaneyko  Leibniz-Institut für Polymerforschung Dresden e.
V., Hohe Straße 6, 01069 Dresden, Germany  J. Ilnytskyi  Institute of Applied
Physics, Technische Universität Dresden, 01062 Dresden, Germany

###### Abstract

We analyze the magnetostriction effect occurring in the magneto-sensitive
elastomers (MSEs) containing inhomogeneously magnetized particles. As it was
shown before, the expression for the interaction potential between two
magnetic spheres, that accounts for their mutual inhomogeneous magnetization,
can be obtained from the Laplace equation. We use this potential in the
approximation formula form to construct magnetic energy of the sample in terms
of the pairwise interactions of the particles. We show that this form of
magnetic energy leads to the same demagnetizing factor as predicted by the
continuum mechanics, confirming that only dipole-dipole magnetic interactions
are important on a large scale. As the next step, we examine the role played
by the particles arrangement on the magnetostriction effect. We consider
different spatial distributions of the magnetic particles: a uniform one, as
well as several lattice-type distributions (SC, BCC, HCP and FCC
arrangements). We show that the particles arrangement affects significantly
the magnetostriction effect if the separation between them became comparable
with the particles’ dimensions. We also show that, typically, this
contribution to the magnetostriction effect is of the opposite sign to the one
related with the initial elastomer shape. Finally, we calculate the
magnetostriction effect using the same interaction potential but expressed in
a form of a series expansion, qualitatively confirming the above findings.

###### keywords:

magnetostriction inhomogeneous bulk magnetization magneto-sensitive elastomer

##  1 Introduction

The phenomenon when the magneto-responsive material changes its shape upon
application of an external magnetic field, is termed “the magnetostriction
effect” [1, 2, 3, 4, 5]. Studies of this phenomenon acquired new momentum in
advent of the interest in composite materials [6, 7, 8], started a few decades
ago. A class of magneto-sensitive elastomers (MSEs) includes composite
materials containing the elastic and magnetic subsystems and, therefore, can
be reversibly controlled via an external magnetic field. It turns out that
uniaxial magnetostriction for MSEs is three orders of magnitude larger than in
pure metals and comparable or higher than “giant magnetostriction” in rare-
earth alloys [9, 10, 11]. This and others unusual mechanical properties of
MSEs made them candidates for application in different devices, e.g., sensors,
dampers and actuators [12, 13, 14, 15]. In that practical context, prediction
of the mechanical properties for these novel materials based on the details of
their physical and structural properties is of a great importance.

The question of: how exactly the inner structure of MSE influences its
macroscopic behavior is a challenging theoretical task, which requires
thorough description of its both elastic and magnetic subsystems. The
properties of the elastic component viewed from a microscopic perspective are
discussed in Ref. [16, 17, 18] and literature therein. More recently, novel
technics, e.g, graph theory [19] and machine learning [20] were used to study
their structure-property relationship. On the other hand, the effects of the
structural order in the magnetic subsystem on the macro-scale response of the
material should also be studied. Relevantly, the investigations of sample
properties affected by the magnetic filler distribution [3, 21, 22, 23, 24]
are performed in recent years. The most common choice for the form of the
magnetic energy is the dipole-dipole interaction [3, 23, 7, 25]. However, it
is being argued in the literature [4, 26] that at the inter- particle
separations comparable with their dimensions, the inhomogeneous bulk
magnetization effects become important. Additionally, it was shown [27, 26],
that inhomogeneously magnetized particles models and the microscopic approach
based in dipole-dipole interaction predict different angular dependencies for
the forces between pair of particles. As a result, one may expect somehow
different predictions of how exactly the magnetic particles distribution
inside an elastomer matrix influence its macroscopic response on the magnetic
field.

Yet another aspect, as was shown in a similar problem of inhomogeneously
polarizable spheres [28], is a non-additivity of the interaction potential and
importance of the three- particle interactions. As a result, the dipole-dipole
interaction model is applicable to the MSEs in the limit of essentially
diluted magnetic filler. At the same time, a general theory which takes into
account the inhomogeneous magnetization of many particles may be impractically
complex.

In the current paper we generalize previous magnetostriction studies of
different spatial distributions of magnetic particles inside an MSE [23]
towards the case of the inhomogeneously magnetizable spheres. The equilibrium
strain of the spheroidal sample in the external magnetic field is considered
as a measure of magnetostriction effect. Dimensionless coefficient
(magnetostriction factor), that defines the sign of the strain, can be split
into two terms: one is related to the initial shape of the elastomer, whereas
the other one accounts for the initial distribution of the magnetic particles
[23, 3]. The magnetostriction factor is evaluated for our model by using the
expression for a pairwise potential energy obtained earlier [29, 4, 26] by
solving the Laplace equation. In our study, we use two forms of this potential
energy: the explicit expression for the first $n=9$ terms of the series
expansion [29], and the numerical approximation based on the first $n=100$
terms [26]. The obtained results are compared with the predictions of the
dipole-dipole interaction model [23] and in each case the magnetostriction
factor is calculated assuming the affine deformation of an MSE sample under
external magnetic field. This comparison is aimed on finding the range of
applicability for the dipole-dipole interaction model for the magnetostriction
effect. We are also interested in the question of which results, provided by
the dipole-dipole model, hold after the inhomogeneous magnetization effect is
taken into account in the pairwise interaction approximation.

Outline of the paper is as follows. In the Sec. 2 we derive the magnetic
energy of the elastomer within the pairwise interaction approximation (PIA).
In the Sec. 3 we show, that for macroscopic samples the shape dependent
magnetostriction is described well by the dipole interactions, while other
terms in the potential energy can be neglected. In the Sec. 4 the distribution
dependent magnetostriction effect is considered for the uniform and the
lattice-like arrangement of the particles. The effects of inhomogeneous
magnetization are analyzed by comparing the magnetostriction factors with
dipole-dipole model prediction. In the Sec. 5 we analyze briefly the
equilibrium strain for various values of external magnetic field and sample
parameters such as its initial shape and the volume fraction of magnetic
filler. The study is eventuated by Conclusions.

##  2 Magnetic energy density and equilibrium strain

Magnetic energy of two equally sized linearly magnetizable spheres can be
obtained by solving the Laplace equation [29, 4, 26]. The solution has a form
of infinite series [27, 4]

$U^{\mathrm{int}}_{\mathrm{LE}}=2U_{0}\sum_{n=3}^{\infty}\bigg{(}\frac{a}{r}\bigg{)}^{n}(a_{n}\cos^{2}{\theta}+b_{n}),$
(1)

where $a$ is the radius of particles and $a/r$ is dimensionless inverse
distance between their centers. The angle between the direction of the
external magnetic field $\bm{H}_{0}$ and the radius-vector $\bm{r}$ is denoted
by $\theta$. The coefficients $a_{n}$, $b_{n}$ are constants, they are
calculated analytically in Ref. [29] up to $n=9$. The first four of them are:
$a_{3}=3\beta_{1}$, $b_{3}=-\beta_{1}$, $a_{4,5}=b_{4,5}=0$,
$a_{6}=3\beta_{1}^{2}$, $b_{6}=\beta_{1}^{2}$.

These series is found to converge rather slowly and, therefore, it was
suggested to keep the terms up to $n=100$ [4, 26] and to approximate it with
the following expression

$U^{\mathrm{int}}=2U_{0}\sum_{k=3}^{7}\\!\beta_{1}^{p_{k}-1}\bigg{[}\frac{a^{k}A_{k}}{(r-aB_{k})^{k}}+\frac{a^{k}C_{k}}{(r-aD_{k})^{k}}\cos^{2}{\theta}\bigg{]},$
(2)

where the exact values of the coefficients $A_{k}$, $B_{k}$, $C_{k}$ and
$D_{k}$ are given in the Table 1 (from [26]).

The coefficient $U_{0}$ in Eqs. (1, 2) is proportional to the magnetic energy
of an isolated sphere under external field

$U_{0}=-\frac{3}{2}\mu\mu_{0}\beta_{1}H^{2}_{0}V_{p},$ (3)

where $\mu=\mu_{\mathrm{medium}}/\mu_{0}$ is the relative permeability of a
medium and $\mu_{0}$ is the permeability of a vacuum. In the case of MSE, the
medium is typically non-magnetic [7, 8] therefore, one can set $\mu=1$. The
factor $\beta_{n}=\frac{n\chi}{n\chi+2n+1}$ is the function of magnetic
susceptibility $\chi$. The factor $V_{p}$ is the volume of a particle
$V_{p}=\frac{4\pi}{3}a^{3}$.

Table 1: Coefficients entering the expression for the potential energy [26] k | $A_{k}$ | $B_{k}$ | $C_{k}$ | $D_{k}$ | $p_{k}$
---|---|---|---|---|---
3 | -1 | 0 | $3$ | $0$ | $2$
4 | 0 | 0 | $3.42\times 10^{-2}$ | $1.2976$ | $3$
5 | 0.111 | -0.689 | $2.83\times 10^{-6}$ | $1.8947$ | $11$
6 | 0.509 | 0.589 | $1.8\times 10^{-13}$ | $1.9898$ | $13$
7 | -0.424 | 0.592 | $0$ | $0$ | $20$

Now, we would like to find the energy density of the sample. To do so, let us
consider $N_{p}$ particles interacting via the pairwise potential (2). The
corresponding magnetic energy of $N_{p}$ particles reads

$U_{N}=N_{p}U_{0}+\frac{1}{2}\underset{i\neq
j}{\sum_{j,i}}U_{ij}^{(\mathrm{int})}.$ (4)

We should note here, that the same equation (1) holds for electrically
polarizable particles. And for that case, the interaction potential is found
to be non-additive [28]. Thus, the results obtained by means of PIA may not
represent accurately the case of very dense systems. However, it is
instructive to estimate the upper limit for the particles density where
dipolar interaction contribution is prevailing. Following our previous study
[27], we consider the sample of ellipsoidal shape which is characterized by a
constant demagnetizing field [30, 31]. Moreover, we consider the arrangement
of particles, such that, the configuration surrounding the $i$-th particle is
independent on the actual number $i$. When both conditions are satisfied, the
expression that remains after summation over $j$, is independent of $i$ and
can be evaluated easily. This results in the expression

$U_{N}=N_{p}U_{0}\bigg{(}1+3\beta_{1}\phi\frac{1}{4\pi c}\sum_{j\neq
i}F(r_{ij})\bigg{)}$ (5)

where $\phi=cV_{p}$ is the volume fraction of the particles and $c=N_{p}/V$ is
their number density. Here we introduced the short-hand functions $F(r_{ij})$
and $F_{k}(r_{ij})$:

$F(r_{ij})=\sum_{k=3}^{7}F_{k}(r_{ij})$ (6)

and

$F_{k}(r_{ij})=\beta_{1}^{p_{k}-2}\bigg{[}\frac{a^{k-3}A_{k}}{(r_{ij}-aB_{k})^{k}}+\frac{a^{k-3}C_{k}}{(r_{ij}-aD_{k})^{k}}\cos^{2}{\theta_{ij}}\bigg{]}.$
(7)

From Eqs. (5)-(7), one can obtain the energy density $u_{N}=U_{N}/V$ as
follows

$u_{N}=u_{0}\bigg{(}1+3\beta_{1}\phi F\bigg{)},$ (8)

where $u_{0}=-\frac{3}{2}\mu_{0}\beta_{1}\phi H_{0}^{2}$. The dimensionless
factor

$F=\frac{1}{4\pi c}\sum\limits_{j\neq i}F(r_{ij})$ (9)

is determined by both the distribution of the particles and the initial shape
of the sample. The similar expression for $u_{N}$ was obtained by us earlier
[27] using the infinite series Eq. (1) which involves evaluation of the factor
$F_{\mathrm{LE}}=\sum_{k}f^{(k)}$. The explicit forms of the first few terms,
$f^{(3)}$, $f^{(6)}$, $f^{(8)}$ and $f^{(9)}$, are known from the previous
studies [29] and will be used in the current study.

_Equilibrium strain._ The total energy density of the MSE under magnetic field
is as follows

$u_{\mathrm{tot}}=\frac{E\varepsilon^{2}}{2}+u_{0}\bigg{(}1+3\beta_{1}\phi
F(\varepsilon)\bigg{)},$ (10)

where the first term represents elastic contribution, and the second one
corresponds to the magnetic energy density of the sample. The
$\varepsilon=(l-l_{0})/l$ is the dimensionless strain of the sample along the
direction of applied magnetic field. The equilibrium strain
$\varepsilon_{\mathrm{eq}}$ is the result of minimization of the energy
density, $\frac{\partial}{\partial\varepsilon}u_{\mathrm{tot}}=0$, which
yields

$\varepsilon_{\mathrm{eq}}=-\frac{3\beta_{1}\phi
u_{0}}{E}\bigg{[}\frac{\partial}{\partial\varepsilon}F(\varepsilon)\bigg{]}_{\varepsilon_{\mathrm{eq}}}.$
(11)

Because $u_{0}\leq 0$, see Eq. (3), the prefactor in front of square brackets
is always positive. In the zero approximation for
$\frac{\partial}{\partial\varepsilon}F(r_{ij}(\varepsilon))$ with respect to
$\varepsilon$, the above equation has the following solution

$\varepsilon_{\mathrm{eq}}^{(1)}=\bigg{.}\varepsilon_{\mathrm{eq}}\bigg{|}_{\varepsilon=0}.$
(12)

It can be seen that the function

$\bigg{[}\frac{\partial}{\partial\varepsilon}F(\varepsilon)\bigg{]}_{\varepsilon=0}=K$
(13)

determines the sign of the strain. The case
$\varepsilon_{\mathrm{eq}}^{(1)}<0$ corresponds to the contraction, whereas
$\varepsilon^{(1)}_{\mathrm{eq}}>0$ corresponds to the elongation of the
sample along the direction of the external field $\bm{H}_{0}$.

Let us compare the obtained expression (11) for the equilibrium strain with
the analogous result obtained in Ref. [3] using dipole-dipole interactions
only. It was assumed that $\bm{m}||\bm{H}_{0}$ and the magnetic energy of a
sample is derived in a self-consistent manner. As a result, the following
expression was obtained in Ref. [3]

$\varepsilon_{\mathrm{eq}}=\frac{1}{2E}\mu_{0}\bigg{[}\frac{\phi^{2}H_{0}^{2}}{\big{(}\chi^{-1}+1/3-\phi
f^{(3)}\big{)}^{2}}\frac{\partial}{\partial\varepsilon}f^{(3)}\bigg{]}_{\varepsilon_{\mathrm{eq}}}.$
(14)

Rewriting the result from Ref. [3] in terms of our notations

$u_{0}=-\frac{3}{2}\mu_{0}\beta_{1}\phi
H^{2}_{0};\quad\beta_{1}=\frac{\chi}{\chi+3},$ (15)

one obtains

$\varepsilon_{\mathrm{eq}}=-\frac{3\beta_{1}\phi
u_{0}}{E}\bigg{[}\frac{1}{\big{(}1-3\beta_{1}\phi
f^{(3)}\big{)}^{2}}\frac{\partial}{\partial\varepsilon}f^{(3)}\bigg{]}_{\varepsilon_{\mathrm{eq}}}.$
(16)

For isotropic particle distributions, including the cases of certain lattice
arrangements, the factor $f^{(3)}$ is bounded from above by $1/3$ [3], and
$\beta_{1}\approx 1$ for carbonyl iron particles [4]. Thus, for small volume
fraction of particles $\phi\ll 1$, one can obtain the following approximation

$\varepsilon_{\mathrm{eq}}\approx-\frac{3\beta_{1}\phi
u_{0}}{E}\bigg{[}\frac{\partial}{\partial\varepsilon}f^{(3)}\bigg{]}_{\varepsilon_{\mathrm{eq}}}.$
(17)

At the same conditions, inhomogeneous magnetization can be ignored and one
obtains the following approximation

$\frac{1}{4\pi c}\sum\limits_{j\neq i}F(r_{ij}(\varepsilon))\approx f^{(3)}.$
(18)

Now, comparing the Eqs. (11, 18) with Eq. (17), one can conclude that both
models predict the same result for the equilibrium strain in the limit of
small volume fraction $\phi\leavevmode\nobreak\ \ll\leavevmode\nobreak\ 1$.

##  3 The influence of the initial shape of the sample on the magnetostriction
effect

It was shown previously [3], that the summation $\frac{1}{4\pi
c}\sum\limits_{j\neq i}$ over the particles positions inside the spheroidal
sample can be split into two terms. First term involves the summation within
the proximity region of the chosen particle (termed as a micro-sphere), where
the result is highly dependent on the particles distribution in space. Second
term is the one outside the micro-sphere, where exact particles positions are
not important, and the integration over this spatial region can be undertaken
instead of summation. Minimal radius of such a micro-sphere was estimated (for
the uniform distribution) to be $R_{\mathrm{ms}}\sim 10\langle r_{n.n}\rangle$
[3], where $\langle r_{n.n}\rangle$ is an average distance between the nearest
neighbors. This formalism is convenient, since it allows a separate treatment
of different effects. For example, the dimensionless factor $F(\varepsilon)$
can be split into two terms

$\displaystyle F(\varepsilon)\\!=\\!\frac{1}{4\pi
c}\sum_{r_{ij}<R_{\mathrm{ms}}}F(r_{ij}(\varepsilon))\\!+\\!\frac{1}{4\pi}\int_{V_{\mathrm{ms}}}^{V_{\mathrm{s.sample}}}dVF(r(\varepsilon))=$
$\displaystyle=F_{\mathrm{micro}}(\varepsilon)+F_{\mathrm{macro}}(\varepsilon),$
(19)

where integration is performed over the entire volume of a sample save the
micro-sphere. The first term provides the contribution to the magnetic energy
from the particles distribution, whereas the second one accounts for the role
of the sample shape. The effect on the magnetostriction can be treated
similarly. The function $K$ from the Eq. (13) can be split as follows [23]

$\\!K\\!\\!=\\!\\!\bigg{[}\frac{\partial}{\partial\varepsilon}F_{\mathrm{micro}}(\varepsilon)\bigg{]}_{\varepsilon=0}\\!\\!+\\!\\!\bigg{[}\frac{\partial}{\partial\varepsilon}F_{\mathrm{macro}}(\varepsilon)\bigg{]}_{\varepsilon=0}\\!\\!\\!=\\!\\!K_{\mathrm{micro}}+K_{\mathrm{macro}},$
(20)

where $K_{\mathrm{micro}}$ and $K_{\mathrm{macro}}$ represent the same
respective contributions as in Eq. (19).

_Shape factor._ Let us consider now the contribution related to the initial
shape in more detail. To do so, we repeat here the procedure outlined in Ref.
[3]. Assuming that the micro-sphere is situated in the center of the
spheroidal sample, the integration in Eq. (19) is greatly simplified

$\frac{1}{4\pi}\int_{V_{\mathrm{ms}}}^{V_{\mathrm{s.sample}}}dV=\frac{1}{4\pi}\int_{0}^{2\pi}d\varphi\int_{0}^{\pi}d\theta\sin{\theta}\int_{R_{\mathrm{ms}}}^{r_{\theta}}drr^{2},$
(21)

where
$r_{\theta}=A\big{[}(1-\gamma^{2}_{\varepsilon})\cos^{2}{\theta}+\gamma^{2}_{\varepsilon}\big{]}^{-\frac{1}{2}}$
is the upper bound of integration over $r$. Parameter
$\gamma_{\varepsilon}=A/B$ is the aspect ratio of a spheroid with the semi-
axes $A$ and $B$ which explicitly depends on the sample strain $\varepsilon$.
According to Table 1, the first term in Eq. (7) is

$F_{3}(r_{ij}(\varepsilon))=\frac{1}{r_{ij}^{3}(\varepsilon)}\bigg{[}3\cos^{2}{\theta_{ij}}(\varepsilon)-1\bigg{]}.$
(22)

It is a dipolar term, and after performing the integration one obtains

$F^{(3)}_{\mathrm{macro}}(\varepsilon)=\frac{1}{3}-N(\gamma_{\varepsilon}),$
(23)

where $N(\gamma_{\varepsilon})$ is a well known demagnetizing factor [23, 30]:

$N(\gamma_{\varepsilon})=\frac{\gamma_{\varepsilon}}{(\gamma_{\varepsilon}^{2}-1)^{3/2}}\ln{\big{(}\gamma_{\varepsilon}+\sqrt{\gamma_{\varepsilon}^{2}-1}\big{)}}-\frac{1}{\gamma_{\varepsilon}^{2}-1}.$
(24)

The demagnetizing factor, $N(\gamma_{\varepsilon})$, is related to the
external surface of the entire sample, whereas $1/3$ term in Eq. (23) is its
analogue, related to the micro-sphere. Let us note that the later term will
cancel out with its counterpart originated from summation over the micro-
sphere. It is clear, that $F^{(3)}_{\mathrm{macro}}(\varepsilon)$ depends on
the aspect ratio $\gamma_{\varepsilon}$ only and is dimensionless. Contrary to
that, the other terms, $F^{(k>3)}_{\mathrm{macro}}(\varepsilon)$, contain the
radius of the particle, $a$, in their respective expressions, which is
microscopic lengths scale. In general, for some $k>3$ one can write

$F_{\mathrm{macro}}^{(k>3)}=\beta_{1}^{p_{k}-2}\\!\\!\int_{0}^{1}\\!\\!\\!dx\int_{R_{\mathrm{ms}}}^{r_{x}}\\!\\!\\!drr^{2}\bigg{[}\frac{a^{k-3}A_{k}}{(r\\!-\\!aB_{k})^{k}}+\frac{a^{k-3}C_{k}}{(r\\!-\\!aD_{k})^{k}}x^{2}\bigg{]},$
(25)

where
$r_{x}=A\big{[}(1-\gamma^{2}_{\varepsilon})x^{2}+\gamma^{2}_{\varepsilon}\big{]}^{-\frac{1}{2}}$
is proportional to the $A$ semi-axis of a spheroid. This integration leads to
the negligible contributions to $F_{\mathrm{macro}}(\varepsilon)$ (see Eq.
(A.1)). As the result, for macroscopic sample one can use Eq. (23)

$F_{\mathrm{macro}}(\varepsilon)=\frac{1}{3}-N(\gamma_{\varepsilon}).$ (26)

In other words, due to the short-range nature of the interactions in
$F^{(k>3)}$, the shape factor is determined solely by a dipolar potential.

_Shape-dependent magnetostriction._ Now, let us turn our attention to the
$K_{\mathrm{macro}}$. First of all, we need to specify the dependency of the
aspect ratio $\gamma_{\varepsilon}$ on the strain $\varepsilon$. If one
assumes the linear affine deformation with respect to a chosen axis, then the
requirement of the volume conservation yields
$\gamma_{\varepsilon}=\gamma_{0}(1+\varepsilon)^{3/2}$. Thus, the derivative
over strain is
$\big{[}\frac{\partial}{\partial\varepsilon}F_{\mathrm{macro}}(\varepsilon)\big{]}_{\varepsilon=0}=\frac{3}{2}\gamma_{0}\frac{\partial}{\partial\gamma_{0}}F_{\mathrm{macro}}(\gamma_{0})$.
One can check (see Eq. (B.9)), that when calculating $K_{\mathrm{macro}}$, the
Eq. (26) can be used. Following the Eq. (24), one can obtain
$K_{\mathrm{macro}}(\gamma_{0})=\frac{\partial}{\partial\varepsilon}F_{\mathrm{macro}}(\varepsilon=0)$
in the following form

$K_{\mathrm{macro}}(\gamma_{0})=\frac{3}{2}\frac{2\gamma_{0}^{3}+\gamma_{0}}{(\gamma_{0}^{2}\\!-\\!1)^{5/2}}\ln{\big{[}\gamma_{0}+\sqrt{\gamma_{0}^{2}\\!-1\\!}\big{]}}-\frac{3}{2}\frac{3\gamma_{0}^{2}}{(\gamma_{0}^{2}\\!-\\!1)^{2}}.$
(27)

This expression was obtained in Ref. [23] previously for the case of dipolar
interactions. As we have shown, it holds for a more general case as well. It
can be seen that the shape dependent magnetostriction factor,
$K_{\mathrm{macro}}$, is always positive, what means the expansion of the
sample along the magnetic field.

In the current chapter we demonstrated the agreement of our results with the
ones of the continuum mechanics theory [32, 31, 30], namely, that the shape
dependent magnetostriction is determined solely by the dipolar interactions.

##  4 The influence of the spatial distribution of magnetic particles on the
magnetostriction effect

In this section we will discuss the contribution to the magnetostriction which
originates from spatial distribution of magnetic particles, given by the
$K_{\mathrm{micro}}$ term in Eq. (20). We assume the case of strong coupling
between magnetic particles and a host polymeric matrix, as in case of magnetic
particles embedded into a polymer via covalent bonds. In this case,
application of the external magnetic field results in an affine deformation of
a sample. The case of weak coupling is also possible when magnetic particles
are just physically embedded inside a soft elastic matrix [33] and move
independently of the matrix. In our case, the effect of different spatial
distributions on the micro-sphere magnetostriction factor,
$K^{\mathrm{affine}}_{\mathrm{micro}}$, is exemplified by the uniform and four
lattice-like particle arrangements.

###  4.1 Approach based on the approximation of the solution for the Laplace
equation

To this end we evaluate the function $F_{k}(r_{ij})$ and the first few terms
$f^{(k)}$, all originated from the interaction energy Eq. (1). Let us consider
the former case in more detail here. The linear affine deformation along the
$Ox$ axis [23], and the volume conservation requirement would lead to the
following relations: $r_{x}=r^{(0)}_{x}(1+\varepsilon)$,
$r_{y}=r^{(0)}_{y}\frac{1}{\sqrt{1+\varepsilon}}$ and
$r_{z}=r^{(0)}_{z}\frac{1}{\sqrt{1+\varepsilon}}$ between the initial
$(r^{(0)}_{x},r^{(0)}_{y},r^{(0)}_{z})$ and after-the- deformation
$(r_{x},r_{y},r_{z})$ positions of each particle. The
$F_{\mathrm{micro}}(\varepsilon)$ is, according to Eq. (19), a sum over the
micro-sphere as follows

$F_{\mathrm{micro}}(\varepsilon)=\frac{1}{4\pi
c}\sum_{r_{ij}<R_{\mathrm{ms}}}\sum_{k=3}^{7}F_{k}(r_{ij}(\varepsilon)).$ (28)

Using the explicit form for the $F(r_{ij}(\varepsilon))$ from Eqs. (6,7), and
rewriting $\cos^{2}{\theta_{ij}}$ as the fraction
$\frac{(r_{x})^{2}_{ij}}{r^{2}_{ij}}$, one can obtain the following

$\displaystyle\bigg{[}\frac{\partial}{\partial\varepsilon}F_{\mathrm{micro}}(\varepsilon)\bigg{]}_{\varepsilon=0}=\frac{1}{\phi}\\!\\!\sum_{r_{ij}<R_{\mathrm{ms}}}\sum_{k=3}^{7}\bigg{[}g_{4,r^{(0)}_{ij}}^{(k)}\cos^{4}{\theta_{ij}}+$
$\displaystyle+g_{2,r^{(0)}_{ij}}^{(k)}\cos^{2}{\theta_{ij}}+g_{0,r^{(0)}_{ij}}^{(k)}\bigg{]},$
(29)

where we used the relation between volume fraction of particles, $\phi$, and
their number density, $c$, in the form of $\frac{3}{4\pi ca^{3}}=1/\phi$. The
functions $g^{(k)}_{m,r^{(0)}_{ij}}$ are, in general, dependent on the
dimensionless inverse distance $a/r^{(0)}_{ij}$. With the help of some
ancillary notation,

$g_{s}(P;R)=\big{(}\frac{a}{r^{(0)}_{ij}}\big{)}^{k}\frac{\beta_{1}^{p_{k}-2}P}{(1-\frac{a}{r^{(0)}_{ij}}R)^{k}}\big{(}\frac{k}{s}\frac{1}{1-\frac{a}{r^{(0)}_{ij}}R}+1\big{)},$
(30)

these functions can be written in relatively compact form as follows

$\displaystyle\\!\\!\\!g^{(k)}_{4,r^{(0)}_{ij}}\\!=$ $\displaystyle-
g_{2}(C_{k};D_{k}),\quad g^{(k)}_{0,r^{(0)}_{ij}}=g_{6}(A_{k};B_{k}),$ (31)
$\displaystyle
g^{(k)}_{2,r^{(0)}_{ij}}\\!\\!=\\!g_{6}(C_{k};D_{k})-g_{2}(A_{k};B_{k}).$

The constants $A_{k},\dots,D_{k}$ can be found in the Table. (1).

One should note, however, that in this affine deformation formalism we
implicitly take into account the particles that are located outside the micro-
sphere after the deformation, or were inside it before the deformation occurs.
We should exclude the effect of the boundary crossing, denoted by
$C_{\mathrm{b.c}}$, from the final result. More details can be found in [23]
as well as in Eq. (B.10). As the result, we obtain the following expression

$K_{\mathrm{micro}}^{\mathrm{affine}}=\bigg{[}\frac{\partial}{\partial\varepsilon}F_{\mathrm{micro}}(\varepsilon)\bigg{]}_{\varepsilon=0}-\frac{2}{5}.$
(32)

Let us consider the contribution to $K_{\mathrm{micro}}^{\mathrm{affine}}$
from the $k=3$ term in Eq. (29). In this case, the coefficients
$g_{m,r^{(0)}_{ij}}^{(3)}$ are greatly simplified and one can obtain the
following

$\\!\\!\\!K_{\mathrm{micro}}^{\mathrm{affine}}(3)=\\!\frac{1}{\phi}\\!\\!\\!\sum_{r^{(0)}_{ij}<R_{\mathrm{ms}}}\\!\\!\\!\\!\frac{-15\cos^{4}{\theta_{ij}}+12\cos^{2}{\theta_{ij}}-1}{2(r^{(0)}_{ij})^{3}}-\frac{2}{5},$
(33)

which coincides with the expression for $K_{\mathrm{micro}}^{\mathrm{affine}}$
obtained for the dipole-dipole interaction in Ref. [23]. Other terms, which
reflect the effect of the inhomogeneous magnetization of the particles can be
expressed as follows

$\displaystyle
K^{\mathrm{affine}}_{\mathrm{micro}}(k>3)=\frac{1}{\phi}\sum_{k=4}^{7}\bigg{[}g_{4,r^{(0)}_{ij}}^{(k)}\cos^{4}{\theta_{ij}}+g_{2,r^{(0)}_{ij}}^{(k)}\cos^{2}{\theta_{ij}}+$
$\displaystyle+g_{0,r^{(0)}_{ij}}^{(k)}\bigg{]}.$ (34)

As a next step, let us consider the sign of the contribution
$K_{\mathrm{micro}}^{\mathrm{affine}}$ to the magnetostriction factor for
different particle distributions.

_Uniform arrangement of magnetic particles._ Micro-sphere is centered around
the randomly chosen $i$-th particle and the pair distances are
$R_{\mathrm{ms}}\geq r^{(0)}_{ij}\geq r^{(0)}_{\mathrm{1s}}$, where
$r^{(0)}_{\mathrm{1s}}$ defines the first coordination sphere. Because of the
finite size of the particles and related to it excluded volume effects, the
genuinely uniform (UNF) distribution can be hardly achieved here. However, we
can assume such distribution on average, namely by substitution the sum as
follows: $\sum_{j,i}U_{ij}^{(\mathrm{int})}=N_{p}\sum_{j}\langle
U_{ij}^{(\mathrm{int})}\rangle$. Conceptually similar smearing off for the
particle positions is reported in Refs. [5, 24]. Thus, we replace the
summation in Eq. (32) by the integration $\frac{1}{c}\sum\rightarrow\int dV$.
The positions of the particles inside integration domain are not correlated
with the absolute value of the radius-vector under uniform distribution
assumption, and the integration over the distance and the solid angle are
independent. Applying the expressions $\frac{1}{4\pi}\int
d\Omega\cos^{2}{\theta}=1/3$ and $\frac{1}{4\pi}\int
d\Omega\cos^{4}{\theta}=1/5$ to the Eq. (33), one can see that
$K_{\mathrm{micro}}^{\mathrm{affine}}$ turns into a constant

$K_{\mathrm{micro}}^{\mathrm{affine}}(3)|_{\mathrm{UNF}}=-\frac{2}{5},$ (35)

which is the known result, obtained earlier in Ref. [23]. To calculate the
other terms $K_{\mathrm{micro}}^{\mathrm{affine}}(k>3)|_{\mathrm{UNF}}$ in Eq.
(34), we use the Eqs. (31,30) and, as the result, obtain the following
expression

$\\!\\!\\!K^{\mathrm{affine}}_{\mathrm{micro}}(k\\!>\\!3)\big{|}_{\mathrm{UNF}}\\!=\\!\sum_{k=4}^{7}\\!\\!\frac{2\beta_{1}^{p_{k}-2}}{15}\frac{C_{k}(y+D_{k})^{3}}{y^{k}}\bigg{|}^{R_{\mathrm{ms}}/a-D_{k}}_{r^{(0)}_{\mathrm{1s}}/a-D_{k}}\\!\\!,$
(36)

where the constants $C_{k}$, $D_{k}$ and $p_{k}$ can be found in the Table 1,
and the change of the variable $y=r^{(0)}_{ij}/a-D_{k}$ was done during the
integration. The result in Eq. (36) does not depend on the constants $A_{k}$
and $B_{k}$ because the corresponding terms cancel out when
$\frac{1}{4\pi}\int d\Omega\cos^{2}{\theta}=1/3$.

Figure 1: Dependencies of the function
$K^{\mathrm{affine}}_{\mathrm{micro}}(k>3)\big{|}_{\mathrm{UNF}}$ from Eq.
(36) on various integration boundaries. In all cases the function is negative.
If lower bound is small, $r^{(0)}_{\mathrm{1s}}=2a$, the function converges to
the value $\sim-0.54$ already for the micro-sphere radius of
$R_{\mathrm{ms}}\sim 3a$. Also, one can observe that the magnitude of the
function diminishes considerably even for small increment of the loser bound
of the integration $r^{(0)}_{\mathrm{1s}}$.

The dependence of this function on the lower and upper integration bounds is
plotted in Fig. (1). It is clear that the function is negative for all the
cases. It converges rather quickly to a finite value when the radius of the
micro-sphere $R_{\mathrm{ms}}/a$ gets larger. For example, it converges to the
value $\sim-0.54$ for the micro-sphere of the radius $R_{\mathrm{ms}}\sim 3a$
when the lower bound of the integral is equal to $r^{(0)}_{\mathrm{1s}}=2a$.
As the lower integration bound increases, the function converges to the lower
value as well, and for $r^{(0)}_{\mathrm{1s}}=2.5a$ the contribution from the
$k>3$ terms is significantly smaller than the dipole $k=3$ contribution, as
given by Eq. (35). And, being the constant, the dipole term only shifts all
graphs down, and do not change their dependence on distances.

_Lattice-like arrangement of magnetic particles._ The lattice-like
arrangements of particles presents great theoretical and experimental
interest. The presence of the well-defined spatial patterns of the particles
positions in this case helps to link them with the properties of the MSE. Such
arrangements can be achieved experimentally. In particular, Zhang et al. [25]
reported the process of preparation the sample with arbitrary lattice-like
arrangement of the particles that are 400-1000 µm in diameter. Quick creation
and transfer into different substrates of $2$D patterns with smaller metallic
particles was reported as well [34], though without mentioning of MSE
preparation.

In order to analyze the influence of the different distributions of the
particles on the magnetostriction effect, it may be convenient to rearrange
the sum, like the one in Eq. (32), into the series of sums performed over the
coordination spheres ($1s,2s,\ldots$). For instance, the summation of the
$\cos^{q}{\theta_{ij}}$ term can be rearranged in the following way

$\sum_{r^{(0)}_{ij}<R_{\mathrm{ms}}}\cos^{q}{\theta_{ij}}=\sum_{r^{(0)}_{ij}\in
1s}\cos^{q}{\theta_{ij}}+\sum_{r^{(0)}_{ij}\in
2s}\cos^{q}{\theta_{ij}}+\cdots$ (37)

where the summation indices, $r^{(0)}_{ij}\in ns$, denote summation over the
$n$-th coordination sphere. In this way one can separate the angular and
radial dependent parts, because in each $n$-th coordination sphere the radius
vectors $r^{(0)}_{ij}$ are constant, $r^{(0)}_{ns}$. Additionally, introducing
the $N_{\mathrm{ns}}$ as the number of particles on the $n$-th coordination
sphere, one can define the average of some angular dependent function
$\langle\cos^{q}{\theta_{ij}}\rangle_{\mathrm{ns}}$ as follows

$\langle\cos^{q}{\theta_{ij}}\rangle_{\mathrm{ns}}=\frac{1}{N_{\mathrm{ns}}}\sum_{r^{(0)}_{ij}\in
ns}\cos^{q}{\theta_{ij}},$ (38)

which is simply the arithmetic mean of such function on the $n$-th
coordination sphere. In general, the average will depend on the particular
spatial distribution of particles.

Figure 2: Dependency of the function
$K^{\mathrm{affine}}_{\mathrm{micro}}\big{|}_{\mathrm{distr.}}$ given by Eq.
(42), on the radius of the first coordination sphere
$r^{(0)}_{\mathrm{1s}}\in[2.02a,2.52a]$. It is calculated for simple cubic
(SC), body-centered cubic (BCC), hexagonal close-packed (HCP) and face-
centered cubic (FCC) lattices. Its asymptotic values, given by the dipolar
contribution including the border crossing term, $2/5$, [23], are shown as
dashed lines. The inset provides a magnified view for the BCC lattice results
at $r^{(0)}_{\mathrm{1s}}\approx 2a$. For this case, the contribution to the
magnetostriction factor from a micro-sphere changes its sign at around
$r^{(0)}_{\mathrm{1s}}\approx 2.1a$.

In the case of the dipolar interparticle interaction $(k\leavevmode\nobreak\
=\leavevmode\nobreak\ 3)$, the magnetostriction factor
$K_{\mathrm{micro}}^{\mathrm{affine}}(k=3)$ was already calculated in Ref.
[23] for the cases of the simple cubic (SC), body-centered cubic (BCC) and
hexagonal close-packed (HCP) lattices inside the micro-sphere of relatively
large size $R_{\mathrm{ms}}$. We complemented these findings by the case of
the FCC lattice, where the magnetic field, $\bm{H}_{0}$, was chosen to be
collinear to one of the edges of a cube. The micro-sphere of the radius
$R_{\mathrm{ms}}=100r_{1s}$ was considered.

The functions $g^{(k)}_{m,r^{(0)}_{ij}}$, entering Eqs. (29,34) and defined
via Eqs. (31,30), alongside with the functions

$\displaystyle
M_{k}^{\mathrm{distr.}}(r^{(0)}_{ns})=\big{\langle}g^{(k)}_{4,r^{(0)}_{ns}}\cos^{4}{\theta_{ij}}+g^{(k)}_{2,r^{(0)}_{ns}}\cos^{2}{\theta_{ij}}+$
$\displaystyle+g^{(k)}_{0,r^{(0)}_{ns}}\big{\rangle}_{\mathrm{ns}}^{\mathrm{distr.}},$
(39)

can be used to define the contribution from the beyond dipolar interactions at
each coordination sphere as follows

$M_{k>3}^{\mathrm{distr.}}(r^{(0)}_{\mathrm{ns}})=\sum_{k=4}^{7}M^{\mathrm{distr.}}_{k}(r^{(0)}_{\mathrm{ns}}).$
(40)

This allows us to write the respective contributions to the magnetostriction
factor

$\displaystyle
K^{\mathrm{affine}}_{\mathrm{micro}}(k>3)\big{|}_{\mathrm{distr.}}=\frac{1}{\phi}\big{[}N_{\mathrm{1s}}M_{k>3}^{\mathrm{distr.}}(r^{(0)}_{\mathrm{1s}})+$
$\displaystyle+N_{\mathrm{2s}}M_{k>3}^{\mathrm{distr.}}(r^{(0)}_{\mathrm{2s}})+\cdots\big{]}.$
(41)

Using Table 2 one can calculate $K^{\mathrm{affine}}_{\mathrm{micro}}(k>3)$
given by the Eq. (41) for different lattice types.

The results discussed above allow us to consider the complete expression for
the magnetostriction factor inside the micro-sphere

$K_{\mathrm{micro}}^{\mathrm{affine}}=K^{\mathrm{affine}}_{\mathrm{micro}}(k=3)+K^{\mathrm{affine}}_{\mathrm{micro}}(k>3).$
(42)

The corresponding plots are shown in the Fig. (2) for the interval
$r^{(0)}_{\mathrm{1s}}=2.02a-2.52a$. It can be seen, that in all cases the
magnitude of this function decreases with the increase of the radius of the
first coordination sphere $r^{(0)}_{\mathrm{1s}}$. This reflects a short-range
nature of the interactions leading to the inhomogeneous magnetization effects.
Similarly, it will decrease as the volume fraction of the particles, $\phi$,
decreases. For the SC, FCC, and HCP lattices, this factor is always negative,
indicating the contribution from a micro-sphere towards the sample contraction
along the field. Contrary to that, $K_{\mathrm{micro}}^{\mathrm{affine}}$
changes sign from positive to negative at around $r^{(0)}_{\mathrm{1s}}\approx
2.1a$ for the BCC lattice case. Beside this feature, in the case of BCC
lattice, the contributions from the first and from the following coordination
spheres are visibly different (see inset in Fig. (2)). We explain this by the
proximity of the first, $r^{(0)}_{\mathrm{1s}}$, and of the second,
$r^{(0)}_{\mathrm{2s}}$, coordination spheres in the BCC lattice:
$r^{(0)}_{\mathrm{2s}}=\sqrt{4/3}r^{(0)}_{\mathrm{1s}}$. Because their
respective contributions to the magnetostriction factor,
$K_{\mathrm{micro}}^{\mathrm{affine}}$, are of the opposite sign, they compete
strongly in a narrow range of distances.

The dipolar contributions are represented by their respective asymptotic
values, $-1.51$, $-0.06$, $-0.21$ and $-0.10$, for the case of SC, BCC, HCP
and FCC lattices, shown as dashed lines in Fig. (2). One can see that they
approximate the function $K_{\mathrm{micro}}^{\mathrm{affine}}$ adequately
when the nearest neighbors distance reaches certain characteristic value. This
is discussed in more detail latter, when comparing magnitude of dipolar and
higher contributions to the micro-sphere magnetization factor
$K_{\mathrm{micro}}^{\mathrm{affine}}$.

Figure 3: Dependence of the higher than dipolar magnetostriction factor,
$K^{\mathrm{affine}}_{\mathrm{micro}}(k>3)$, given by the Eq. (41), on the
radius of the first coordination sphere in the interval
$r^{(0)}_{1s}\in[2a;2.02a]$.

At the interparticle separation equal to the first coordination sphere radius,
$r^{(0)}_{\mathrm{1s}}=2a$, the function
$K^{\mathrm{affine}}_{\mathrm{micro}}$, except for the case of the BCC
lattice, acquires gigantic negative values (see the Fig. (3)). Despite the
fact that the forces between inhomogeneously magnetized particles can be much
larger than the dipolar ones [26], we expect that this prediction for
magnetostriction factor is unphysical. This, in part, is due to the fact, that
in our model we cut off the interaction potential at the contact distances and
do not introduce the short-range elastic repulsion [35]. Accounting for such
repulsion terms may correct the model prediction at small interparticle
distances. Moreover, at this close separations between particles, one expects
a deviation from the affine deformation, as well as the increased significance
of the three- and higher particle interactions.

Figure 4: Dependencies of the functions $M(r^{(0)}_{\mathrm{1s}};\theta)$ on
the angle $\theta$ for various values of the first coordination sphere radius
$r^{(0)}_{\mathrm{1s}}/a$. As distance between particles become smaller, the
differences of the magnitudes become more apparent. The effect is more
pronounced for the small angles $\theta\sim 0^{\circ}$, which corresponds to
the chain-like conformation of the particles.

As was discussed above, the behavior of the micro-sphere magnetostriction
factor, $K^{\mathrm{affine}}_{\mathrm{micro}}$, exhibits a strong dependence
on the lattice symmetry. This is especially true at the distances comparable
to the first coordination sphere, and the contributions to
$K^{\mathrm{affine}}_{\mathrm{micro}}(k>3)$ provided by the interactions on
this lengths scale, could be sufficient for the description of their
qualitative behavior. Concerning the dipolar term in Eq. (42), its properties
also can be inferred from the first coordination sphere contribution [23]. As
can be seen in Eq. (41), the contribution to the magnetostriction factor from
the nearest neighbors, ($r^{(0)}_{ij}=r^{(0)}_{\mathrm{1s}}$), is determined
by the functions $M^{\mathrm{distr.}}_{k}(r^{(0)}_{\mathrm{1s}})$ given by the
Eq. (39). Therefore, let us introduce functions
$M(r^{(0)}_{\mathrm{1s}};\theta)$ defined via relations

$M(r^{(0)}_{\mathrm{1s}};\theta)=\sum_{k=3}^{7}M_{k}(r^{(0)}_{\mathrm{1s}};\theta),$
(43) $\langle
M_{k}(r^{(0)}_{\mathrm{1s}};\theta)\rangle^{\mathrm{distr.}}_{\mathrm{1s}}=\frac{3}{4}(r^{(0)}_{\mathrm{1s}}/a)^{3}M^{\mathrm{distr.}}_{k}(r^{(0)}_{\mathrm{1s}}).$
(44)

The angle $\theta$ is defined in Eq. (1) above. The prefactor
$(r^{(0)}_{\mathrm{1s}}/a)^{3}$ is used in order to account for the volume
fraction of the particles, $\phi\propto(a/r^{(0)}_{\mathrm{1s}})^{3}$,
entering the Eq. (41). This definition is convenient for the purpose of the
direct comparison with the earlier results. Namely, using the information from
the Table 1 and Eq. (39), the dipolar contribution can be written explicitly
as follows

$M_{3}(r^{(0)}_{\mathrm{1s}};\theta)=\frac{3}{8}\big{(}-15\cos^{4}{\theta}+12\cos^{2}{\theta}-1\big{)}.$
(45)

This function is studied in Ref. [23], and now we can compare it to the other
contributions, arising from the inhomogeneous bulk magnetization of the
particles. The dipolar term, $M_{3}(r^{(0)}_{\mathrm{1s}};\theta)$, which is
independent on the radius of the first coordination sphere,
$r^{(0)}_{\mathrm{1s}}$, is showed as a dashed line in the Fig. (4). The
function $M(r^{(0)}_{\mathrm{1s}};\theta)$ start to deviate from the dipolar
one when the distance $r_{\mathrm{1s}}$ is less than five particle radii. This
effect is more pronounced for the chain-like ($\theta\sim 0^{\circ}$)
arrangement of the particles. Qualitatively, the chain-like ($\theta\sim
0^{\circ}$) and plane-like ($\theta\sim 90^{\circ}$) conformations correspond
to the negative values of the function $M(r^{(0)}_{\mathrm{1s}};\theta)$. For
the intermediate particle conformation the interval exists where the function
is positive. In the case of the dipolar term
$M_{3}(r^{(0)}_{\mathrm{1s}};\theta)$, this region can be identified as
$\theta\sim 40^{\circ}-70^{\circ}$. For the combined function
$M(r^{(0)}_{\mathrm{1s}};\theta)$ the interval where it is positive can
slightly vary with increasing radius of the first coordination sphere,
exhibiting shift towards the larger values of the angle $\theta$.

###  4.2 Approach based on the accounting for explicit terms in the solution
of the Laplace equation

As was briefly discussed in Sec. 2, the magnetic energy density of the sample
can be written in the following form

$u_{N}=u_{0}\big{(}1+3\beta_{1}\phi F_{\mathrm{LE}}\big{)},$ (46)

where $F_{\mathrm{LE}}$ was obtained from the Eq. (1) and has the form of
infinite series $F_{\mathrm{LE}}=\sum_{k}f^{(k)}$ of dimensionless terms
$f^{(k)}$ [27]. The explicit expressions for the first four non-zero terms
$f^{(k)}$, entering the expression for $F_{\mathrm{LE}}$, are known [29]. This
allows us to derive the micro-sphere magnetostriction factor
$K^{\mathrm{affine}}_{\mathrm{micro}}$ in an alternative way to the one based
on Eq. (8) and considered in this paper above.

In a way similar to Eq. (19), the dimensionless terms $f^{(k)}$ each can be
split into two parts, namely $f^{(k)}_{\mathrm{micro}}$, and
$f^{(k)}_{\mathrm{macro}}$. The “macro” contribution is discussed in more
detail in Ref. [27]. In this section we consider the “micro” terms. The first
four of them are as follows [29, 27]

$f_{\mathrm{micro}}^{(3)}=\frac{1}{4\pi
c}\sum_{r_{ij}<R_{\mathrm{ms}}}\frac{3(r_{ij})_{x}^{2}-r^{2}_{ij}}{r_{ij}^{5}},$
(47) $f_{\mathrm{micro}}^{(6)}=\frac{\beta_{1}a^{3}}{4\pi
c}\sum_{r_{ij}<R_{\mathrm{ms}}}\frac{3(r_{ij})_{x}^{2}+r^{2}_{ij}}{r_{ij}^{8}},$
(48) $f_{\mathrm{micro}}^{(8)}=\frac{3\beta_{2}a^{5}}{4\pi
c}\sum_{r_{ij}<R_{\mathrm{ms}}}\frac{2(r_{ij})_{x}^{2}+r^{2}_{ij}}{r_{ij}^{10}},$
(49) $f_{\mathrm{micro}}^{(9)}=\frac{\beta^{2}_{1}a^{6}}{4\pi
c}\sum_{r_{ij}<R_{\mathrm{ms}}}\frac{9(r_{ij})_{x}^{2}-r^{2}_{ij}}{r_{ij}^{11}}.$
(50)

Assuming the affine deformation, as it was the case in the previous sections,
one can examine the respective magnetostriction factor
$K_{\mathrm{micro}}^{\mathrm{affine}}(k)$ for each of these terms. As
expected, the case of $k=3$ corresponds to a dipolar contribution

$K_{\mathrm{micro}}^{\mathrm{affine}}(3)=\bigg{[}\frac{\partial}{\partial\varepsilon}f_{\mathrm{micro}}^{(3)}(\varepsilon)\bigg{]}_{\varepsilon=0}-C^{(3)}_{\mathrm{b.c}}=$
(51) $\\!\\!\\!=\frac{1}{4\pi
c}\\!\\!\\!\sum_{r^{(0)}_{ij}<R_{\mathrm{ms}}}\\!\\!\\!\frac{1}{2(r^{(0)}_{ij})^{3}}\bigg{[}-45\cos^{4}{\theta_{ij}}+36\cos^{2}{\theta_{ij}}-3\bigg{]}-\frac{2}{5},$
(52)

where value $C^{(3)}_{\mathrm{b.c}}=2/5$ of border crossing term [23] is used.
In the case of $k>3$, the border crossing terms are proportional to the
negligibly small fractions $\big{(}\frac{a}{R_{\mathrm{ms}}}\big{)}^{k}$,
therefore we assume $C^{(k>3)}_{\mathrm{b.c}}\rightarrow 0$. It may be
convenient to introduce the angle dependent functions $a_{k}(\theta_{ij})$ as
follows

$\displaystyle a_{3}(\theta_{ij})$
$\displaystyle=(1/8)(-45\cos^{4}{\theta_{ij}}+36\cos^{2}{\theta_{ij}}-3),$
(53) $\displaystyle a_{6}(\theta_{ij})$
$\displaystyle=(1/4)(-36\cos^{4}{\theta_{ij}}+9\cos^{2}{\theta_{ij}}+3),$
$\displaystyle a_{8}(\theta_{ij})$
$\displaystyle=(3/4)(-30\cos^{4}{\theta_{ij}}+2\cos^{2}{\theta_{ij}}+4),$
$\displaystyle a_{9}(\theta_{ij})$
$\displaystyle=(9/8)(-33\cos^{4}{\theta_{ij}}+18\cos^{2}{\theta_{ij}}-1).$

Those functions, Eq. (53), can be used to rewrite the magnetostriction factor
as a series containing contributions from all coordination spheres $ns$ inside
the micro-sphere

$\\!\\!\\!\\!K^{\mathrm{affine}}_{\mathrm{micro}}(k)|_{\mathrm{distr.}}=\frac{4b_{k}}{3\phi}\\!\\!\sum_{r^{(0)}_{\mathrm{ns}}<R_{\mathrm{ms}}}\\!\\!\\!\\!N_{\mathrm{ns}}\Big{(}\frac{a}{r^{(0)}_{\mathrm{ns}}}\Big{)}^{k}\langle
a_{k}(\theta_{ij})\rangle_{\mathrm{ns}}^{\mathrm{distr.}},$ (54)

where the constants $b_{k}$ are equal to $1,\beta_{1},\beta_{2}$ and
$\beta^{2}_{1}$ for $k=3,6,8$ and $k=9$, respectively. The averaging in
$\langle
a_{k}(\theta_{ij})\rangle^{\mathrm{distr.}}_{\mathrm{ns}}=\frac{1}{N_{\mathrm{ns}}}\sum_{r^{(0)}_{ij}\in
r_{\mathrm{ns}}}a_{k}(\theta_{ij})$, see Eq. (54), is performed over the
$n$-th coordination sphere. Since those average values are dependent on the
distribution pattern of the particles inside an elastomer, so does the
magnetostriction factor
$K^{\mathrm{affine}}_{\mathrm{micro}}(k)|_{\mathrm{distr.}}$. In order to gain
some qualitative insight, it may be instructive to investigate the angular
dependency of the magnetostriction factor for the first coordination sphere
$r^{(0)}_{\mathrm{1s}}$. One can note that the volume fraction of the
particles $\phi\propto(\frac{a}{r^{(0)}_{\mathrm{1s}}})^{3}$, thus, the
rescaled functions, similar to those in Eq. (44), can be introduced

$m(r^{(0)}_{\mathrm{1s}};\theta)=\sum_{k=3}^{9}m_{k}(r^{(0)}_{\mathrm{1s}};\theta),$
(55)

where

$m_{k}(r^{(0)}_{\mathrm{1s}};\theta)=b_{k}\big{(}\frac{a}{r^{(0)}_{\mathrm{1s}}}\big{)}^{k-3}a_{k}(\theta).$
(56)

In the Eq. (55) we sum the first four nonzero contributions ($k=3,6,8$ and $9$
see Eq. (53)).

Figure 5: Comparison of the angular dependencies for the functions
$m(r^{(0)}_{\mathrm{1s}};\theta)$ from Eq. (55) and
$M(r^{(0)}_{\mathrm{1s}};\theta)$ from Eq. (44). The former function is
derived from the series expansion Eq. (1) truncated at $n=9$, the latter one
originates from the numerical approximation [4, 26] of the series expansion
truncated at $n=100$. Several cases for different radii of the first
coordination sphere, $r^{(0)}_{\mathrm{1s}}$, are shown.

The dependence of $m(r^{(0)}_{\mathrm{1s}};\theta)$ on the angle $\theta$ is
compared to the one considered earlier in Eq. (44). The results of the
comparison are plotted in the Fig. (5). It can be seen that, at the distance
between the nearest neighbors $r^{(0)}_{\mathrm{1s}}=5a$, both functional
forms are undistinguishable from the prediction of the dipolar potential (see
Fig. (4)). When the radius of the first coordination sphere approaches
$r^{(0)}_{\mathrm{1s}}\sim 2.2a-2.3a$, the differences between
$m(r^{(0)}_{\mathrm{1s}};\theta)$ and $M(r^{(0)}_{\mathrm{1s}};\theta)$ become
apparent. For the larger particle separations, $r^{(0)}_{\mathrm{1s}}>2.4a$,
these functions are almost identical. Interestingly, the behavior of
$M(r^{(0)}_{\mathrm{1s}};\theta)$ for $r^{(0)}_{\mathrm{1s}}=2.4a$ and
$m(r^{(0)}_{\mathrm{1s}};\theta)$ for $r^{(0)}_{\mathrm{1s}}=2.3a$ looks very
similar, especially for small values of the angle $\theta<30^{\circ}$. This
may be related to the form of the approximation formula in Eq. (2), where the
terms with $k>3$ quickly converge. Thus, the increment of the radius vector
values may effectively correspond to the approximation of the series in Eq.
(1) with smaller number of terms $n$. However, we haven’t studied this in more
detail.

Now, turning our attention to the Eq. (54), we can sum up the contribution
from the dipolar term, $k=3$, and the known higher order terms, $k=6-9$, to
obtain the complete micro-sphere magnetostriction factor

$K^{\mathrm{affine}}_{\mathrm{micro}}(k\leq
9)=K^{\mathrm{affine}}_{\mathrm{micro}}(3)+\sum_{k=6}^{9}K^{\mathrm{affine}}_{\mathrm{micro}}(k),$
(57)

where the first four non-zero terms in the series expansion Eq. (1) were used
[29]. In a such framework, the inhomogeneous magnetization effects are
accounted for in the $k>3$ part of the sum. Since those terms decay with the
increase of the distance $r^{(0)}_{\mathrm{ns}}$, we can restrict ourself to
the first few coordination spheres only. Contrary to that, the dipolar part of
the magnetostriction factor,
$K^{\mathrm{affine}}_{\mathrm{micro}}(3)|_{\mathrm{distr.}}$, exhibits complex
behavior, which requires accounting for large number of coordination spheres
$r^{(0)}_{\mathrm{ns}}$. This property of the dipolar magnetostriction factor
may originate from the interplay of the radial and angular dependent parts,
see Eq. (54). The radial part is proportional to $\propto 1/r^{3}$, which
makes dipolar interactions long-range. However, in the case of uniform
distribution, the angular dependent part turns to zero, which makes each term
in the sum from Eq. (54) equal to zero as well. The dipolar magnetostriction
factor was calculated for the simple cubic (SC), body-centered cubic (BCC) and
hexagonal close-packed (HCP) lattices inside a micro-sphere
$R_{\mathrm{ms}}\gg r^{(0)}_{\mathrm{1s}}$ in Ref. [23]. For the face-centered
cubic (FCC) lattice we performed respective calculations within the interior
of the micro-sphere, given by its radius,
$R_{\mathrm{ms}}=100r^{(0)}_{\mathrm{1s}}$. The function
$K^{\mathrm{affine}}_{\mathrm{micro}}(k\leq 9)$, given by the Eq. (57), is
plotted in the Fig. (6).

Figure 6: Dependency of the function
$K^{\mathrm{affine}}_{\mathrm{micro}}(k\leq 9)$ from Eq. (57) on the radius of
the first coordination sphere $r^{(0)}_{\mathrm{1s}}$. With the exception of
the BCC lattice, at $ns=1$, magnetostriction factor is negative. Additionally,
the cases of the uniform (UNF) particles arrangement as well as the SC, FCC
and HCP lattices are shown. With the increase of the interparticle separation,
$r^{(0)}_{\mathrm{1s}}$, the function quickly converge to the dipole-dipole
model predictions.

It can be seen that $K^{\mathrm{affine}}_{\mathrm{micro}}(k\leq 9)$ quickly
converges to the respective dipole-dipole model prediction, as the radius of
the first coordination sphere increases. The case of UNF distribution was
calculated by replacing the summation inside a micro-sphere with the
integration, under the same assumptions as in Eq. (35). It can be seen, that
for sufficiently large interparticle separation, $r^{(0)}_{\mathrm{1s}}/a>3$,
the magnetostriction factor approaches zero, as expected. For larger densities
of the particles it became negative. The BCC and FCC lattice distributions are
the most sensitive to the number of coordination spheres $ns$ accounted for in
the calculation. If only the first coordination sphere is considered, the
magnetostriction factor $K^{\mathrm{affine}}_{\mathrm{micro}}(k\leq 9)$ for
BCC lattice changes sign. It is positive for the denser systems,
$r^{(0)}_{\mathrm{1s}}/a<2.5$, and negative for less dense ones,
$r^{(0)}_{\mathrm{1s}}/a>3$. The BCC lattice arrangement leads to the negative
magnetostriction factor for all densities, if first two or three coordination
spheres are considered. Comparing this with the results for the BCC lattice in
Fig. (2), one can see that accounting for more terms ($n=100$) in the series
in Eq. (1) increases the weight of the contribution from the first
coordination sphere. Contrary to that, the magnetostriction factor for the SC
lattice is almost independent on the number of coordination spheres
considered. Some magnetostriction factors depicted in Fig. (6) are smaller
than the contribution from the initial shape discussed in Sec. 3. For example,
for the initially spherical sample one predicts
$K^{\mathrm{affine}}_{\mathrm{macro}}=2/5$. As a result, for the FCC lattice
distribution with smaller densities, $r^{(0)}_{\mathrm{1s}}/a>2.4$, the
“macro” term is dominant, predicting the net elongation of the sample, but for
higher densities, $r^{(0)}_{\mathrm{1s}}/a<2.3$, the “micro” contribution is
larger, predicting contraction of an MSE. In the case of the BCC lattice
arrangement of particles inside a spherical sample, the elastomer is expected
to elongate for all densities. Both the SC and HCP cases lead to the
contraction of the spherical sample, upon application of the external magnetic
field $\bm{H}_{0}$. The elastomers with the other initial shapes are
considered in the next section, see Fig (10).

###  4.3 Contribution from the dipolar interactions at various volume
fractions of magnetic particles

Figure 7: Ratio between the dipolar contribution and the complete micro-sphere
magnetostriction factor $R_{\mathrm{dip.}}$, given by the Eq. (58), as a
function of the first coordination sphere radius $r^{(0)}_{\mathrm{1s}}$.
Respective ratios $R_{\mathrm{dip.}}$ for various lattices are shown. For each
lattice types, two cases are shown. The first one shows $R_{\mathrm{dip.}}$
derived numerically according to Eq. (1) with $n=100$ terms used (solid lines
and symbols). The second one shows the same ratio derived from Eq. (1) where
$n=9$ terms are used explicitly.

Now, let us estimate the influence of inhomogeneous magnetization on the
magnetostriction effect. To this end, we introduce the ratio between the
dipolar contribution and the complete micro-sphere magnetostriction factor

$R_{\mathrm{dip.}}=\frac{|K^{\mathrm{affine}}_{\mathrm{micro}}(3)|}{|K^{\mathrm{affine}}_{\mathrm{micro}}(3)|+|K^{\mathrm{affine}}_{\mathrm{micro}}(k>3)|},$
(58)

that hereafter will be termed as the “dipolar ratio”. When both dipolar and
higher order terms in $K^{\mathrm{affine}}_{\mathrm{micro}}$ are of the same
sign, it can be rewritten as
$K^{\mathrm{affine}}_{\mathrm{micro}}(3)/K^{\mathrm{affine}}_{\mathrm{micro}}(k\geq
3)$. In general, dipolar ratio $R_{\mathrm{dip.}}$ is the function of the
first coordination sphere radius $r^{(0)}_{\mathrm{1s}}$, and we plot this
dependency in the Fig. (7). Two cases are considered: (i) when the dipolar
ratio is derived from the numerical approximation formula of the first $n=100$
terms in the series expansion, given by the Eq. (1); (ii) when it is derived
from the first nine ($n=9$) terms explicitly. In both cases, the first three
coordination spheres are used in the calculation of the magnetostriction
factor higher than the dipolar one,
$K^{\mathrm{affine}}_{\mathrm{micro}}(k>3)$. We should note here that the
magnetostriction factor derived from the nine ($n=9$) terms is more sensitive
to the number of coordination spheres $ns$ accounted for in calculation, which
is also evident comparing Figs. (2) and (6). As before, we use dipolar
magnetostriction factor calculated for SC, BCC, HCP lattices in Ref. [23],
while FCC lattice result we obtain by summation inside the micro-sphere with
the radius $R_{\mathrm{ms}}=100r^{(0)}_{\mathrm{1s}}$. It can be seen that, in
$n=100$ case, the inhomogeneous magnetization effects become more pronounced
for small values of $r^{(0)}_{\mathrm{1s}}$. Interestingly, the more terms
$n=100$ in the Eq. (1) lead to “more short-range” inhomogeneous magnetization
effects, than for the series truncated at $n=9$. For the case $n=100$, at
distances between the nearest neighbors $r^{(0)}_{\mathrm{1s}}>3a$, the
dipolar term comprises more than $90\%$ of magnetostriction factor. This
distances correspond to different volume fractions $\phi$ for different
lattices, as can be seen in the Table (2). The dipole ratio
$R_{\mathrm{dip.}}$ dependency on $\phi$ is shown in the Fig. (8) for the
$n=100$ case.

Figure 8: Dependency of the dipolar ratio $R_{\mathrm{dip.}}$, given by the
Eq. (58), on the volume fraction of the particles $\phi$ inside an elastomer
is shown. For each lattice, respective plots terminate where the maximal
possible packing of spherical particles is reached. In the case of the BCC
lattice, the dipolar ratio is always greater than $\approx 10\%$.

It can be seen, in Fig. (8), that for the SC, FCC, and HCP lattices, the
dipolar ratio $R_{\mathrm{dip.}}$ is negligible for dense systems. And for the
BCC lattice it is always greater than one-tenth of the total magnetostriction
factor. The maximal possible value of the volume fraction
$\phi_{\mathrm{max}}$ can be reached when the interparticle separation is
minimal $r^{(0)}_{\mathrm{1s}}=2a$. In the case of the SC and BCC lattices,
the $\phi_{\mathrm{max}}\approx 0.52$ and $\approx 0.68$, respectively, while
for both the HCP and FCC lattices $\phi_{\mathrm{max}}\approx 0.74$. In the
case of BCC lattice, the dipolar contribution comprises more than $90\%$ of
the complete micro-sphere magnetostriction factor for the volume fraction less
than $\phi<0.3$. Similar value of dipolar ratio $R_{\mathrm{dip.}}$ is reached
when $\phi<0.2$ in the case of the SC, FCC and HCP lattices. The SC and FCC
arrangements predict similar dipolar contribution to the magnetostriction
effect on the interval $\phi\in(0;0.35)$.

In the current section we showed an influence of the inhomogeneous
magnetization on the magnetostriction in the framework of PIA. The effect of
finite size of the particles become important when they are sufficiently close
to each other. For example, in Fig. (7) it can be seen that, the dipolar
contribution comprises more than $90\%$ of the complete magnetostriction
factor, when $r^{(0)}_{\mathrm{1s}}>3.2a$. This interparticle separation
corresponds roughly to the volume fraction $\phi<0.2$ for most of the lattices
considered, see Fig. (8). The exception being the BCC lattice, where the
limiting interparticle separation is $r^{(0)}_{\mathrm{1s}}>2.8a$ and the
respective volume fraction $\phi<0.3$, see Figs. (7,8). This results are
obtained in the framework involving approximation formula for $n=100$ terms of
the Laplace equation solution, given by Eq. (1). They are also compared with
the results originated from the explicit treatment of the first $n=9$ terms in
Eq. (1). It can be seen in Fig. (7), that the latter approach underestimates
dipolar contribution to magnetostriction effect for larger values of
$r^{(0)}_{\mathrm{1s}}$ (less dense distributions). It also overestimates the
importance of dipole-dipole interaction when particles become close to each
other $r^{(0)}_{\mathrm{1s}}<2.4a$. This may be related to the slow
convergence of the sum in Eq. (1), reported in Refs. [4, 26] earlier. In our
analysis we ignore the three- and higher particle interactions and non-
additive nature of inhomogeneous magnetization effects, reported for similar
electrically polarizable particles [28]. However, we speculate that the
influence on these interactions with more than two neighboring particles may
actually make the magnetization more uniform. And the range of applicability
of the dipolar model may increase to the higher volume fractions than those
roughly estimated here.

##  5 Estimation of the equilibrium strain

In this chapter we would like to estimate the equilibrium strain, given be the
Eq. (12), for the spheroidal MSE with the magnetic particles distributed in
the vertices of four different lattices. Similar studies where performed
earlier, where only dipolar interactions where accounted for [23]. Due to
inhomogeneous magnetization of the particles, the dependency of
$K^{\mathrm{affine}}_{\mathrm{micro}}$ on the volume fraction $\phi$ may
arise. This is expected to create richer magnetostriction behavior of the
MSEs. In our framework we accounted for inhomogeneous magnetization in PIA,
which may not be accurate enough in order to make confident predictions about
behavior of the MSE. The aim is, however, to show the connection between our
theoretical model and measurable macroscopic properties of the sample. And
also to demonstrate that the equilibrium strain dependency on the volume
fraction, $\varepsilon_{eq}\propto\phi^{2}$, may deviate from the quadratic
law in principle, since it can be interpreted as a manifestation of the
inhomogeneous magnetization of the particles.

Figure 9: The values of the equilibrium strain, given by the Eq. (12), for
varying strength of the magnetic field $H_{0}$ and volume fraction of
particles $\phi$. The initial aspect ratio $\gamma_{0}=1$ and the Young’s
modulus $E=2$ Mpa of the elastomer are fixed.

To this end, we use the values of the magnetostriction factors
$K^{\mathrm{affine}}_{\mathrm{micro}}$ calculated in the previous section.
Additionally, we specify some physical parameters. We would consider the
sample with fixed effective Young’s modulus $E=2$ MPa. The magnetic field
applied along the axis of the spheroidal sample is considered in the range
$H_{0}\in[80;125]$ kA/m, which is below the saturation field for carbonyl iron
particles. The volume fraction of the particles
$\phi\propto(a/r^{(0)}_{\mathrm{1s}})^{3}$ is another varying parameter. Also
we consider different aspect ratios of the initial shape of the spheroid
$\gamma_{0}\in[0.1;4]$. The quantitative measure of the magnetostriction
effect is the equilibrium strain $\varepsilon_{\mathrm{eq}}$, which the
elastomer attains when the external magnetic field is applied.

Figure 10: The values of the equilibrium strain, given by the Eq. (12), for
varying initial aspect ratio $\gamma_{0}$ and the volume fraction $\phi$. The
magnetic field strengths and Young’s modulus of the elastomer are fixed at
$H_{0}=125$ kA/m and $E=2$ MPa, respectively.

The dependencies of the equilibrium strain, $\varepsilon_{\mathrm{eq}}$, on
the volume fraction of particles $\phi$ and on the magnitude of the external
magnetic field $H_{0}$ are shown in Fig. (9). It can be seen, that the
magnitude of the effect is lower than several percents. For higher values of
the field and larger volume fractions of the particles, the equilibrium strain
for the SC and HCP lattices reaches $2\%$. The sign of the magnetostriction
effect for those two lattices is negative. For the BCC lattice, one can see
that the equilibrium strain is of smaller amplitude, reaching $0.5\%$, while
the sign is positive. In the case of the FCC lattice, the magnitude of the
equilibrium strain is $0.2\%$, and the sign of the magnetostriction effect is
different for different set of the parameters. For the intermediate volume
fractions the sign is positive, while for very densely packed particles the
effect is predicted to be negative. We want to note, that FCC lattice is
considered positioned in a way that the magnetic field is directed parallel to
the side of the cube. If we considered FCC lattice rotated to the close-packed
conformation, where direction of the field would be collinear to the diagonal
of the cube, the predicted magnetostriction effect would be more similar to
the one shown for the HCP lattice case.

The dependency of the equilibrium strain, $\varepsilon_{\mathrm{eq}}$, on the
initial aspect ratio $\gamma_{0}$ of a sample and the volume fraction $\phi$
for each lattice is shown in the Fig. (10). The effective Young’s modulus of a
sample is assumed to be fixed, and equal to $E=2$ MPa. It can be seen that the
SC and HCP lattices are similar and predict negative magnetostriction effect
with $\varepsilon_{\mathrm{eq}}$ larger than $2\%$ and $1\%$, respectively.
The equilibrium strain decreases significantly as the volume fraction of the
particles, $\phi$, becomes smaller. For the BCC lattice, our model predicts a
positive sign for the magnetostriction, and the effect is more pronounced for
the slightly oblate spheroids. In the case of the FCC lattice, the oblate
spheroids would elongate along the applied magnetic field, positive
magnetostriction. But the prolate initial shape of the elastomer would lead to
the contraction of the sample along applied field, negative magnetostriction.
The magnitude of the effect is again smaller than for the other lattices
reaching $\varepsilon_{\mathrm{eq}}\sim 0.2\%$.

##  6 Conclusions

In the current paper we start from the magnetic energy of the pair of linearly
magnetized spheres. At small interparticle separations these particles
inhomogeneously magnetize each other. In order to account for such
inhomogeneous magnetization one need to solve the Laplace equation for this
spherical particles. The approximate formula, designed to fit the solution of
the Laplace equation, given by the Eq. (1) and truncated at its hundredth
term, was reported in [26]. We use this formula in the PIA to construct the
magnetic energy density of the spheroidal elastomer. Then, we study the effect
of the magnetostriction for the MSE sample upon application of the external
magnetic field $\bm{H}_{0}$. The direction of $\bm{H}_{0}$ is chosen to be
collinear to the symmetry axis of the sample. The elastic energy density of
the elastomer is obtained assuming that the sample behaves as a spring with
Young’s modulus $E$. In our model, the magnetic energy density, as well as the
interparticle distances, depend on the strain $\varepsilon$ of the elastomer.
The equilibrium strain $\varepsilon_{\mathrm{eq}}$ is evaluated within the
linear approximation for the magnetic energy density. Then, we analyze how the
initial shape and the distribution of particles inside the MSE affect sign and
amplitude of the equilibrium strain.

The contribution to the magnetostriction effect from the initial shape of the
sample is determined by the magnetostriction factor $K_{\mathrm{macro}}$, (see
Eq. (20)). In the current study we show that only the dipolar interactions
provides a significant contribution to the $K_{\mathrm{macro}}$, as long as
MSE is much larger than the particle radius. This result is in agreement with
the continuum mechanics approach.

The magnetostriction factor inside a micro-sphere,
$K^{\mathrm{affine}}_{\mathrm{micro}}$, describes the contribution to the
magnetostriction effect from the initial distribution of the particles. We
found that the interactions of a higher order than dipolar, originated from
the effects of inhomogeneous magnetization, are short ranged. Thus, for volume
fractions of the particles, $\phi<0.2$, in the case of SC, FCC, and HCP
lattices, the dipolar term comprises more than $90\%$ of the magnetostriction
factor. The same is true for $\phi<0.3$, in the case of BCC lattice, as can be
seen Fig. (8).

At small interparticle separations, when first coordination sphere radius is
less than $r^{(0)}_{1s}\approx 3a$, the inhomogeneous magnetiation of the
particles leads to $K^{\mathrm{affine}}_{\mathrm{micro}}$ dependence on the
volume fraction of the particles, $\phi$. Since such behavior is not predicted
when only dipolar interactions are taken into account, it may be thought as
the manifestation of the effect of inhomogeneous magnetization. Consequences
of this can be observed in Figs. (9,10), where dependence of the equilibrium
strain, $\varepsilon_{\mathrm{eq}}$, on the volume fraction, $\phi$, is shown.

Since $K^{\mathrm{affine}}_{\mathrm{micro}}$ contains terms proportional to
$\cos^{4}{\theta_{ij}}$, it is dependent on the orientation of the lattice
relatively to the direction of the magnetic field $\bm{H}_{0}$. Such effect
can be illustrated by the differences between the FCC and HCP lattices, Figs.
(9,10). The orientation of FCC lattice was chosen such that the edge of the
cube is collinear with $\bm{H}_{0}$. If the field direction was collinear with
its diagonal, the differences between those lattices would have been minor.

When particles arrangement is close to a close contact, $r^{(0)}_{1s}\approx
2a$, the SC, FCC and HCP lattices predict giant negative values of
magnetostriction factor $K^{\mathrm{affine}}_{\mathrm{micro}}$, see Fig. (3).
We treat this as an artefact of our model. We assume strong particle-matrix
coupling and affine deformation of the sample. At short distances, when a
small amount of elastic component is present, this assumption is not valid any
more. Additionally, our magnetic energy density is derived under the pairwise
interaction approximation, which in the case of inhomogeneously magnetized
particles, may not be sufficiently accurate [28].

In order to make our model more realistic, it is important to consider other
laws of elastic reaction to the microscopic movement of magnetic particles.
Also, energy of multiple inhomogeneously magnetized particles should be
considered without pairwise interaction approximation. This may shed light
onto the limitations and the range of applicability of our model. Also, it
would be interesting to calculate equilibrium strain above linear
approximation of magnetic energy. This would led to additional effects, for
example, for the change of an effective Young modulus as magnetic particles
rearrange under application of the external field.

Acknowledgements
D.Y. thanks the National Academy of Sciences of Ukraine for funding this
research (the grant for research laboratories/groups of young scientists No
07/01-2022(4)).
D.Y. and J.I are grateful to the Armed Forces of Ukraine for the protection
during this research work.

## References

  * [1] J. P. Joule.  On the effects of magnetism upon the dimensions of iron and steel bars.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 30(199):76–87, 1847.
  * [2] J. M. Ginder, M. E. Nichols, L. D. Elie, and S. M. Clark.  Controllable-stiffness components based on magnetorheological elastomers.  In Norman M. Wereley, editor, Smart Structures and Materials 2000: Smart Structures and Integrated Systems, volume 3985, pages 418 – 425\. International Society for Optics and Photonics, SPIE, 2000.
  * [3] D. Ivaneyko, V. Toshchevikov, M. Saphiannikova, and G. Heinrich.  Mechanical properties of magneto-sensitive elastomers: unification of the continuum-mechanics and microscopic theoretical approaches.  Soft Matter, 10:2213–2225, 2014.
  * [4] A. Biller, O. Stolbov, and Yu. Raikher.  The forces between magnetizable particles embedded in an elastomer.  Computational Continuum Mechanics, 7(1):61–72, 2014.  URL https://doi.org/10.7242/1999-6691/2014.7.1.7.
  * [5] D. Romeis, V. Toshchevikov, and M. Saphiannikova.  Elongated micro-structures in magneto-sensitive elastomers: a dipolar mean field model.  Soft Matter, 12(46):9364–9376, 2016.
  * [6] Z. Rigbi and L. Jilkén.  The response of an elastomer filled with soft ferrite to mechanical and magnetic influences.  Journal of Magnetism and Magnetic Materials, 37(3):267–276, 1983.  ISSN 0304-8853.
  * [7] M. R. Jolly, J. D. Carlson, B. C. Muñoz, and T. A. Bullions.  The magnetoviscoelastic response of elastomer composites consisting of ferrous particles embedded in a polymer matrix.  Journal of Intelligent Material Systems and Structures, 7(6):613–622, 1996.
  * [8] J. D. Carlson and M. R. Jolly.  Mr fluid, foam and elastomer devices.  Mechatronics, 10(4):555–569, 2000.  ISSN 0957-4158.
  * [9] S. Bednarek.  The giant linear magnetostriction in elastic ferromagnetic composites within a porous matrix.  Journal of Magnetism and Magnetic Materials, 301(1):200 – 207, 2006.  ISSN 0304-8853.
  * [10] N. C. Koon, C. M. Williams, and B. N. Das.  Giant magnetostriction materials.  Journal of Magnetism and Magnetic Materials, 100(1):173–185, 1971.  ISSN 0304-8853.
  * [11] A. E. Clark and H. S. Belson.  Giant room-temperature magnetostrictions in tb${\mathrm{fe}}_{2}$ and dy${\mathrm{fe}}_{2}$.  Phys. Rev. B, 5:3642–3644, May 1972.
  * [12] N. C. Koon, C. M. Williams, and B. N. Das.  Giant magnetostriction materials.  Journal of Magnetism and Magnetic Materials, 100(1):173–185, 1991.  ISSN 0304-8853.
  * [13] R. Ahamed, S.-B. Choi, and M. M. Ferdaus.  A state of art on magneto-rheological materials and their potential applications.  Journal of Intelligent Material Systems and Structures, 29(10):2051–2095, 2018.
  * [14] Y.-C. Li, J.-C. Li, W.-H. Li, and H.-P. Du.  A state-of-the-art review on magnetorheological elastomer devices.  Smart Materials and Structures, 23(12):123001, November 2014.
  * [15] S. S. Kang, K. Choi, J.-D. Nam, and H. J. Choi.  Magnetorheological elastomers: Fabrication, characteristics, and applications.  Materials, 13(20), 2020.  ISSN 1996-1944.
  * [16] P. J. Flory.  Molecular theory of rubber elasticity.  Polymer Journal, 17(1):1–12, Jan 01, 1985.  ISSN 1349-0540.
  * [17] B. Erman and J. E. Mark.  Structures and Properties of Rubberlike Networks.  Oxford University Press on Demand, 1997.
  * [18] K. Urayama.  An experimentalist’s view of the physics of rubber elasticity.  Journal of Polymer Science Part B: Polymer Physics, 44(24):3440–3444, 2006.
  * [19] Y. Amamoto, K. Kojio, A. Takahara, Y. Masubuchi, and T. Ohnishi.  Complex network representation of the structure-mechanical property relationships in elastomers with heterogeneous connectivity.  Patterns, 1(8):100135, 2020.  ISSN 2666-3899.
  * [20] N. E. Jackson, M. A. Webb, and J. J. de Pablo.  Recent advances in machine learning towards multiscale soft materials design.  Current Opinion in Chemical Engineering, 23:106–114, 2019\.  ISSN 2211-3398.  Frontiers of Chemical Engineering: Molecular Modeling.
  * [21] D. Ivaneyko, V. Toshchevikov, M. Saphiannikova, and G. Heinrich.  Magneto-sensitive elastomers in a homogeneous magnetic field: A regular rectangular lattice model.  Macromolecular Theory and Simulations, 20(6):411–424, 2011.
  * [22] D. Ivaneyko, V. Toshchevikov, M. Saphiannikova, and G. Heinrich.  Effects of particle distribution on mechanical properties of magneto-sensitive elastomers in a homogeneous magnetic field.  Condensed Matter Physics, 15(3):33601, September 2012.
  * [23] D. Romeis, V. Toshchevikov, and M. Saphiannikova.  Effects of local rearrangement of magnetic particles on deformation in magneto-sensitive elastomers.  Soft Matter, 15(17):3552–3564, 2019.
  * [24] S. Chougale, D. Romeis, and M. Saphiannikova.  Magneto-mechanical enhancement of elastic moduli in magnetoactive elastomers with anisotropic microstructures.  Materials, 15(2), 2022.  ISSN 1996-1944.
  * [25] X. Zhang, S. Peng, W. Wen, and W. Li.  Analysis and fabrication of patterned magnetorheological elastomers.  Smart Materials and Structures, 17(4):045001, may 2008.
  * [26] A. Biller, O. Stolbov, and Yu. Raikher.  Modeling of particle interactions in magnetorheological elastomers.  Journal of Applied Physics, 116(11):114904, 2014\.  URL https://doi.org/10.1063/1.4895980.
  * [27] D. Yaremchuk, V. Toshchevikov, J. Ilnytskyi, and M. Saphiannikova.  Magnetic energy and a shape factor of magneto-sensitive elastomer beyond the point dipole approximation.  Journal of Magnetism and Magnetic Materials, 513:167069, 2020.  ISSN 0304-8853.
  * [28] H. J. H. Clercx and G. Bossis.  Many-body electrostatic interactions in electrorheological fluids.  Phys. Rev. E, 48:2721–2738, Oct 1993.
  * [29] D. J. Jeffrey.  Conduction through a random suspension of spheres.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 335(1602):355–367, November 1973.
  * [30] J. A. Osborn.  Demagnetizing factors of the general ellipsoid.  Phys. Rev., 67:351–357, Jun 1945.
  * [31] J. C. Maxwell.  A Treatise on Electricity and Magnetism, Vol. 2.  Clarendon Press, 1881.  62-64 pp.  https://archive.org/details/treatiseonelect02maxwrich.
  * [32] L. D. Landau and E. M. Lifshitz.  Electrodynamics of Continuous Media, volume 8.  Pergamon Press, 1960.
  * [33] Th. Gundermann and S. Odenbach.  Investigation of the motion of particles in magnetorheological elastomers by X-µCT.  Smart Materials and Structures, 23(10):105013, sep 2014.
  * [34] S. Zhang, Y. Zhai, R. Peng, M. Shayegannia, A. G. Flood, J. Qu, X. Liu, N.P. Kherani, and A.R. Wheeler.  Assembly of topographical micropatterns with optoelectronic tweezers.  Advanced Optical Materials, 7(20):1900669, 2019\.
  * [35] S. Goh, A. M. Menzel, R. Wittmann, and H. Löwen.  Density functional approach to elastic properties of three-dimensional dipole-spring models for magnetic gels.  The Journal of Chemical Physics, 158(5):054909, 2023.
  * [36] H. A. Lorentz.  The theory of electrons and its applications to the phenomena of light and radiant heat.  GE Stechert & Company, 1916.
  * [37] J. Bourgain, Z. Rudnick, and P. Sarnak.  Spatial statistics for lattice points on the sphere i: Individual results, 2016.
  * [38] W. Duke.  Hyperbolic distribution problems and half-integral weight maass forms.  Inventiones mathematicae, 92(1):73–90, Feb 01, 1988.  ISSN 1432-1297.

## Appendix A. Integral

We are interested in calculating the integral

$F_{\mathrm{macro}}^{(k>3)}=\beta_{1}^{p_{k}-2}\int_{0}^{1}dx\int_{R_{\mathrm{ms}}}^{r_{x}}drr^{2}\bigg{[}\frac{a^{k-3}A_{k}}{(r-aB_{k})^{k}}+\frac{a^{k-3}C_{k}}{(r-aD_{k})^{k}}x^{2}\bigg{]}.$
(A.1)

Let us consider the following integral first

$\int_{R_{\mathrm{ms}}}^{r_{x}}drr^{2}\frac{a^{k-3}A_{k}}{(r-aB_{k})^{k}}=a^{k-3}A_{k}\int_{R_{\mathrm{ms}}-aB_{k}}^{r_{x}-aB_{k}}dy\frac{(y+aB_{k})^{2}}{y^{k}}=$
$=\frac{A_{k}}{3-k}\Big{[}\frac{a^{k-3}}{(r_{x}-aB_{k})^{k-3}}-\frac{a^{k-3}}{(R_{\mathrm{ms}}-aB_{k})^{k-3}}\Big{]}+\frac{2A_{k}B_{k}}{2-k}\Big{[}\frac{a^{k-2}}{(r_{x}-aB_{k})^{k-2}}-\frac{a^{k-2}}{(R_{\mathrm{ms}}-aB_{k})^{k-2}}\Big{]}+$
$+\frac{A_{k}B^{2}_{k}}{1-k}\Big{[}\frac{a^{k-1}}{(r_{x}-aB_{k})^{k-1}}-\frac{a^{k-1}}{(R_{\mathrm{ms}}-aB_{k})^{k-1}}\Big{]}.$
(A.2)

Or

$\int_{R_{\mathrm{ms}}}^{r_{x}}drr^{2}\frac{a^{k-3}A_{k}}{(r-aB_{k})^{k}}=\sum_{l=1}^{3}2^{1-|l-2|}\frac{A_{k}B_{k}^{3-l}}{l-k}\Big{[}\frac{a^{k-l}}{(r_{x}-aB_{k})^{k-l}}-\frac{a^{k-l}}{(R_{\mathrm{ms}}-aB_{k})^{k-l}}\Big{]}.$
(A.3)

And the original integral

$F_{\mathrm{macro}}^{(k>3)}=\beta_{1}^{p_{k}-2}\sum_{l=1}^{3}\frac{2^{1-|l-2|}}{l-k}\Big{(}\int_{0}^{1}dx\Big{[}\frac{a^{k-l}A_{k}B_{k}^{3-l}}{(r_{x}-aB_{k})^{k-l}}+\frac{a^{k-l}C_{k}D_{k}^{3-l}x^{2}}{(r_{x}-aD_{k})^{k-l}}\Big{]}-\frac{a^{k-l}A_{k}B_{k}^{3-l}}{(R_{\mathrm{ms}}-aB_{k})^{k-l}}-\frac{1}{3}\frac{a^{k-l}C_{k}D_{k}^{3-l}}{(R_{\mathrm{ms}}-aD_{k})^{k-l}}\Big{)}$
(A.4)

As long as $\frac{a}{A}\ll 1$, the terms which contain $r_{x}$ can be
neglected. First, one can expand the denominator into Taylor series to get

$\int_{0}^{1}dx\frac{a^{k-l}}{(r_{x}-aB_{k})^{k-l}}\rightarrow\big{(}\frac{a}{A}\big{)}^{k-l}\int_{0}^{1}dx\big{[}(1-\gamma_{\varepsilon}^{2})x^{2}+\gamma_{\varepsilon}^{2}\big{]}^{(k-l)/2}\bigg{(}1+\frac{aB_{k}(k-l)}{A}\big{[}(1-\gamma_{\varepsilon}^{2})x^{2}+\gamma_{\varepsilon}^{2}\big{]}^{1/2}+\cdots\bigg{)}.$
(A.5)

Than, multiplying the factors $\bigg{(}\frac{a}{A}\bigg{)}^{n}$ to the unity
in the form $1=\frac{1+\gamma_{\varepsilon}^{n}}{1+\gamma_{\varepsilon}^{n}}$,
one get the sum of the terms like the following

$\bigg{[}\bigg{(}\frac{a}{A}\bigg{)}^{n}+\bigg{(}\frac{a}{B}\bigg{)}^{n}\bigg{]}\frac{1}{1+\gamma_{\varepsilon}^{n}}\int_{0}^{1}dx\big{[}(1-\gamma_{\varepsilon}^{2})x^{2}+\gamma_{\varepsilon}^{2}\big{]}^{n/2}.$
(A.6)

Similar procedure can be done for
$\int_{0}^{1}dxx^{2}\frac{a^{k-l}}{(r_{x}-aD_{k})^{k-l}}$. The corresponding
functions
$\frac{1}{1+\gamma_{\varepsilon}^{n}}\int_{0}^{1}dx\big{[}(1-\gamma_{\varepsilon}^{2})x^{2}+\gamma_{\varepsilon}^{2}\big{]}^{n/2}$
and
$\frac{1}{1+\gamma_{\varepsilon}^{n}}\int_{0}^{1}dxx^{2}\big{[}(1-\gamma_{\varepsilon}^{2})x^{2}+\gamma_{\varepsilon}^{2}\big{]}^{n/2}$,
emerging in the process, are bound for all $\gamma_{\varepsilon}$. Thus, for
macroscopic sample condition

$\bigg{[}\bigg{(}\frac{a}{A}\bigg{)}^{n}+\bigg{(}\frac{a}{B}\bigg{)}^{n}\bigg{]}\ll
1,$ (A.7)

the following integrals tend to zero

$\int_{0}^{1}dx\frac{a^{k-l}}{(r_{x}-aB_{k})^{k-l}}\rightarrow
0;\qquad\int_{0}^{1}dx\frac{a^{k-l}x^{2}}{(r_{x}-aD_{k})^{k-l}}\rightarrow 0.$
(A.8)

As can be seen from the last term in Eq. (A.4) and the assumption that radius
of the microsphere is in order of ten interparticle distances
($R_{\mathrm{ms}}\geq 10\langle r_{ij}\rangle$), one can obtain
$\frac{a}{R_{\mathrm{ms}}}<0.05$ at least, with $A_{4}=0$, $C_{4}=0.0342$ the
largest term in $F_{\mathrm{macro}}^{(k>3)}$ would be of the order $10^{-3}$
or smaller. We want to note also that the lower bound of the integral would be
the upper bound of the sum, thus
$F(\varepsilon)=F_{\mathrm{micro}}(\varepsilon)+F_{\mathrm{macro}}(\varepsilon)$
from Eq. (19) does not actually depend on the $R_{\mathrm{ms}}$ if it is large
enough to allow sum interchange with integral.

## Appendix B. Higher $K_{\mathrm{macro}}$

Higher than ordinary dipole terms in $F_{\mathrm{macro}}^{(k>3)}$ may lead to
the higher terms in the factor
$K_{\mathrm{macro}}^{(k>3)}=\big{[}\frac{\partial}{\partial\varepsilon}F_{\mathrm{macro}}^{(k>3)}\big{]}_{\varepsilon=0}$.
It turns out, that for affine deformations this higher terms are zero in the
macroscopic sample limit. To check this, one may note first, that terms
containing $R_{\mathrm{ms}}$ in the denominator in Eq. (A.4) do not depend on
the strain $\varepsilon$, and vanish after taking the derivative. And since
$\big{[}\frac{\gamma_{\varepsilon}}{\partial\varepsilon}\big{]}_{\varepsilon=0}=\frac{3}{2}\gamma_{0}$,
it is left to show that

$\gamma_{0}\frac{\partial}{\partial\gamma_{0}}\int_{0}^{1}dx\frac{a^{k-l}}{\big{(}r_{x}(\gamma_{0})-aB_{k}\big{)}^{k-l}}\rightarrow
0;\qquad\gamma_{0}\frac{\partial}{\partial\gamma_{0}}\int_{0}^{1}dx\frac{a^{k-l}x^{2}}{\big{(}r_{x}(\gamma_{0})-aD_{k}\big{)}^{k-l}}\rightarrow
0.$ (B.1)

Using $\frac{a}{A}\ll 1$ as in Eq. (A.5), one can get Taylor series to obtain

$\gamma_{0}\frac{\partial}{\partial\gamma_{0}}\int_{0}^{1}dx\frac{a^{k-l}}{(r_{x}(\gamma_{0})-aB_{k})^{k-l}}\rightarrow\sum_{n=1}C^{B_{k}}_{n}\bigg{[}\bigg{(}\frac{a}{A}\bigg{)}^{n}+\bigg{(}\frac{a}{B}\bigg{)}^{n}\bigg{]}\frac{1}{1+\gamma_{0}^{n}}\gamma_{0}\frac{\partial}{\partial\gamma_{0}}\int_{0}^{1}dx\big{[}(1-\gamma_{0}^{2})x^{2}+\gamma_{0}^{2}\big{]}^{n/2}$
(B.2)

and

$\gamma_{0}\frac{\partial}{\partial\gamma_{0}}\int_{0}^{1}dx\frac{a^{k-l}x^{2}}{(r_{x}(\gamma_{0})-aD_{k})^{k-l}}\rightarrow\sum_{n=1}C^{D_{k}}_{n}\bigg{[}\bigg{(}\frac{a}{A}\bigg{)}^{n}+\bigg{(}\frac{a}{B}\bigg{)}^{n}\bigg{]}\frac{1}{1+\gamma_{0}^{n}}\gamma_{0}\frac{\partial}{\partial\gamma_{0}}\int_{0}^{1}dxx^{2}\big{[}(1-\gamma_{0}^{2})x^{2}+\gamma_{0}^{2}\big{]}^{n/2},$
(B.3)

where $C^{B_{k}}_{n}$ and $C^{D_{k}}_{n}$ are constants and the prefacrots
$\big{(}\frac{a}{A}\big{)}^{n}$ are multiplied by
$\frac{1+\gamma_{0}^{n}}{1+\gamma_{0}^{n}}$. The last step is to show that the
functions on $\gamma_{0}$ in the right-hand side of Eqs. (B.2,B.3) are bounded
for all $\gamma_{0}\geq 0$. Let us consider the function from Eq. (B.2) as
follows

$\frac{\gamma_{0}}{1+\gamma_{0}^{n}}\frac{\partial}{\partial\gamma_{0}}\int_{0}^{1}dx\big{[}(1-\gamma_{0}^{2})x^{2}+\gamma_{0}^{2}\big{]}^{n/2}=\frac{n}{2}\frac{2\gamma_{0}^{2}}{1+\gamma_{0}^{n}}\int_{0}^{1}dx(1-x^{2})\big{[}(1-\gamma_{0}^{2})x^{2}+\gamma_{0}^{2}\big{]}^{(n-2)/2}.$
(B.4)

This function is bounded for large $\gamma_{0}$. For small $\gamma_{0}$ it is
bounded for $n\geq 2$. In the case of $n=1$ and small $\gamma_{0}$ we need to
perform integration first, to avoid the dependency on $x$ before taking
derivative and taking small $\gamma_{0}$ limit. The integral would be as
follows

$\int_{0}^{1}dx\big{[}(1-\gamma_{0}^{2})x^{2}+\gamma_{0}^{2}\big{]}^{1/2}=\frac{1}{2}+\frac{1}{2}\frac{\gamma_{0}^{2}}{\sqrt{1-\gamma_{0}^{2}}}\mathrm{Arcth}\big{(}\sqrt{\gamma_{0}^{-2}-1}\big{)}.$
(B.5)

The corresponding derivative is

$\gamma_{0}\frac{\partial}{\partial\gamma_{0}}\int_{0}^{1}dx\big{[}(1-\gamma_{0}^{2})x^{2}+\gamma_{0}^{2}\big{]}^{1/2}=\frac{\gamma_{0}^{2}}{2}\frac{2-\gamma_{0}^{2}}{(1-\gamma_{0}^{2})^{3/2}}\mathrm{Arcth}\big{(}\sqrt{\gamma_{0}^{-2}-1}\big{)}-\frac{1}{2}\frac{\gamma_{0}^{2}}{1-\gamma_{0}^{2}}.$
(B.6)

Using L’Hôpital’s rule, one get $\lim_{\gamma_{0}\rightarrow
0}\gamma^{2}_{0}\mathrm{Arcth}\big{(}\sqrt{\gamma_{0}^{-2}-1}\big{)}\rightarrow
0$, and thus

$\lim_{\gamma_{0}\rightarrow
0}\frac{1}{1+\gamma_{0}}\gamma_{0}\frac{\partial}{\partial\gamma_{0}}\int_{0}^{1}dx\big{[}(1-\gamma_{0}^{2})x^{2}+\gamma_{0}^{2}\big{]}^{1/2}\rightarrow
0.$ (B.7)

As a result, the function in Eq. (B.2) is bounded for all $\gamma_{0}\geq 0$.
Now, let us consider the function from Eq. (B.3) as follows

$\frac{\gamma_{0}}{1+\gamma_{0}^{n}}\frac{\partial}{\partial\gamma_{0}}\int_{0}^{1}dxx^{2}\big{[}(1-\gamma_{0}^{2})x^{2}+\gamma_{0}^{2}\big{]}^{n/2}=\frac{n}{2}\frac{2\gamma_{0}^{2}}{1+\gamma_{0}^{n}}\int_{0}^{1}dxx^{2}(1-x^{2})\big{[}(1-\gamma_{0}^{2})x^{2}+\gamma_{0}^{2}\big{]}^{(n-2)/2}.$
(B.8)

This function is also bounded for all $\gamma_{0}\geq 0$. Finally, for the
macroscopic sample limit
$\bigg{[}\bigg{(}\frac{a}{A}\bigg{)}^{n}+\bigg{(}\frac{a}{B}\bigg{)}^{n}\bigg{]}\ll
1$ one can write

$\bigg{[}\frac{\partial}{\partial\varepsilon}F_{\mathrm{macro}}^{(k>3)}\bigg{]}_{\varepsilon=0}\rightarrow
0.$ (B.9)

One of the consequences is that the border crossing terms defined in Ref. [23]
can be calculated from the dipole interactions while others can be ignored.
Using the same procedure as in Ref. [23], the border crossing terms explicitly
can be written as follows

$C_{\mathrm{b.c}}=\frac{3}{2}\bigg{[}\gamma_{0}\frac{\partial}{\partial\gamma_{0}}F_{\mathrm{macro}}(\gamma_{0})\bigg{]}_{\gamma_{0}=1}=-\frac{3}{2}\int_{0}^{1}dx(3x^{2}-1)(1-x^{2})=\frac{2}{5}.$
(B.10)

## Appendix C. Short note on the averaging of the cosine in the lattices

Let us consider $n$-th coordination sphere around $i$-th particle in the
lattice with $N_{\mathrm{ns}}$ particles on it. The radius-vector
$\bm{r}_{ij}$ which connects the particle in the center of the coordination
sphere and some particle laying on it would be directed to some angle
$\theta_{ij}$ with respect to the external magnetic field $\bm{H}_{0}$. Let us
chose the magnetic field to be directed along the $Ox$ axis. Then, one may
calculate the average of the $\cos^{2}{\theta_{ij}}$ on the $n$-th
coordination sphere as follows

$\langle\cos^{2}{\theta_{ij}}\rangle_{\mathrm{ns}}=\frac{1}{N_{\mathrm{ns}}}\sum_{j\in
N_{\mathrm{ns}}}\frac{(r_{ij}^{2})_{x}}{r^{2}_{ij}}.$ (C.1)

One can see that

$\frac{1}{N_{\mathrm{ns}}}\sum_{j\in
N_{\mathrm{ns}}}\frac{(r_{ij}^{2})_{x}+(r_{ij}^{2})_{y}+(r_{ij}^{2})_{z}}{r^{2}_{ij}}=1.$
(C.2)

Thus, for simple cubic (SC), body centered cubic (BCC) and face-centered cubic
(FCC) lattices, where all three directions are equivalent, the average of the
square of cosine would be as follows

$\langle\cos^{2}{\theta}\rangle_{\mathrm{ns}}=1/3.$ (C.3)

This result holds for any site $i$ and any coordination sphere $n$ of the
aforementioned lattices. This calculations are essentially the same as in the
Note $55$ in [36]. The average value of the fourth power of the cosine

$\langle\cos^{4}{\theta_{ij}}\rangle_{\mathrm{ns}}=\frac{1}{N_{\mathrm{ns}}}\sum_{j\in
N_{\mathrm{ns}}}\frac{(r_{ij}^{4})_{x}}{r^{4}_{ij}},$ (C.4)

does not have such a nice property for those lattices, and depends on the
particular coordination sphere. In the case of SC, for example, the several
first values are compared with the averaging for the homogeneous distribution
$\langle\cos^{4}{\theta_{ij}}\rangle_{\mathrm{ns}}=1/5$ in the Fig. (11).

Figure 11: Several first values of the average of $cos^{4}{\theta_{ij}}$ for
the simple cubic lattice (SC) are compared with the value of the same function
averaged for homogeneous particle distribution. The most significant
deviations are for the coordination sphere numbers of the form $n=4^{n}$,
which would happen more and more rarely as the $n$ grows. The points
corresponding to the values of $n=7$ and $n=15$ are missing, as they should,
due to Legendre’s three-squares theorem.

At first glance, the graph does not look very encouraging, but the most
diverging values correspond to $n=1,4,16$ or $n=4^{n}$ and those would happen
more and more rarely as $n$ grows. And if we understood the remarks found in
Ref. [37, 38] correctly, the points on unit sphere become uniformly
distributed for large $n\neq 0,4,7\leavevmode\nobreak\
\mathrm{mod}\leavevmode\nobreak\ 8$. Average value of the
$\cos^{4}{\theta_{ij}}$ for some lattices, up third coordination sphere can be
found in the Table 2.

Lattice type | $\langle\cos^{4}{\theta_{ij}}\rangle_{\mathrm{1s}}$ | $\langle\cos^{4}{\theta_{ij}}\rangle_{\mathrm{2s}}$ | $\langle\cos^{4}{\theta_{ij}}\rangle_{\mathrm{3s}}$ | $r_{\mathrm{2s}}$ | $r_{\mathrm{3s}}$ | $N_{\mathrm{1s}}$ | $N_{\mathrm{2s}}$ | $N_{\mathrm{3s}}$ |  $\phi$ (cubic cell)
---|---|---|---|---|---|---|---|---|---
SC | $1/3$ | $1/6$ | $1/9$ | $\sqrt{2}r_{\mathrm{1s}}$ | $\sqrt{3}r_{\mathrm{1s}}$ | 6 | 12 | 8 | $\frac{4\pi a^{3}}{3r^{3}_{\mathrm{1s}}}$
BCC | $1/9$ | $1/3$ | $1/6$ | $\sqrt{4/3}r_{\mathrm{1s}}$ | $\sqrt{8/3}r_{\mathrm{1s}}$ | 8 | 6 | 12 | $\frac{\sqrt{3}\pi a^{3}}{r^{3}_{\mathrm{1s}}}$
FCC | $1/6$ | $1/3$ | $2/9$ | $\sqrt{2}r_{\mathrm{1s}}$ | $\sqrt{3}r_{\mathrm{1s}}$ | 12 | 6 | 24 | $\frac{4\sqrt{2}\pi a^{3}}{3r^{3}_{\mathrm{1s}}}$
HCP | $2/9$ | $1/3$ | $1$ | $\sqrt{2}r_{\mathrm{1s}}$ | $\sqrt{8/3}r_{\mathrm{1s}}$ | 12 | 6 | 2 | $\frac{4\sqrt{2}\pi a^{3}}{3r^{3}_{\mathrm{1s}}}$
Table 2:

Interestingly, for the two dimensional triangular lattice the analogous result
$\langle\cos^{2}{\theta}\rangle=1/2$ can be obtained from the property of the
particle distribution on the coordination circles of this lattice. Namely, the
particles can be grouped in sets where they would be each equally spaced from
their neighbors. In other words, they would be situated in the vertices of the
regular polygon (in each group). And one can see that for the polygon with
$N_{\mathrm{ns}}\geq 3$ sides, the average would be as follows

$\langle\cos^{2}{\theta}\rangle_{\mathrm{ns}}=\frac{1}{N_{\mathrm{ns}}}\sum_{k=0}^{N_{\mathrm{ns}}-1}\cos^{2}{(\theta_{0}+2\pi
k/N_{\mathrm{ns}})},$ (C.5)

where $\theta_{0}$ is the arbitrary angle the polygon may be rotated around
its center. Using power reducing formula

$\langle\cos^{2}{\theta}\rangle_{\mathrm{ns}}=\frac{1}{2N_{\mathrm{ns}}}\sum_{k=0}^{N_{\mathrm{ns}}-1}\big{[}1+\cos{(2\theta_{0}+4\pi
k/N_{\mathrm{ns}})}\big{]},$ (C.6)

and rewriting the cosine as en exponent $\cos{\theta}=\mathrm{Re}e^{i\theta}$
one may obtain the following

$\langle\cos^{2}{\theta}\rangle_{\mathrm{ns}}=\frac{1}{2}+\frac{1}{2N_{\mathrm{ns}}}\mathrm{Re}\leavevmode\nobreak\
e^{i2\theta_{0}}\sum_{k=0}^{N_{\mathrm{ns}}-1}q^{k}=1/2.$ (C.7)

Here we use the geometric sum formula with $q=e^{4\pi/N_{\mathrm{ns}}}\neq 1$
and $q^{N_{\mathrm{ns}}}=1$. The similar procedure can be performed to
calculate the fourth power average

$\langle\cos^{4}{\theta}\rangle_{\mathrm{ns}}=\frac{3}{8}+\frac{1}{8N_{\mathrm{ns}}}\mathrm{Re}\sum_{k=0}^{N_{\mathrm{ns}}-1}\big{[}e^{i4\theta_{0}}q^{2k}+4e^{i2\theta_{0}}q^{k}\big{]},$
(C.8)

where the $q=e^{4\pi/N_{\mathrm{ns}}}$, so the last term in parenthesis gives
zero contribution for $N_{\mathrm{ns}}\geq 3$. In the case of the triangular
lattice, the first term in the parenthesis also vanishes, since $q^{2k}\neq
1$, as a result

$\langle\cos^{4}{\theta}\rangle_{\mathrm{ns}}^{\mathrm{triang.latt.}}=\frac{3}{8}.$
(C.9)

If the average in each group is the same
$\sum_{1}/N_{1}=\sum_{2}/N_{2}=\cdots=\alpha$, then it is also the average of
the all values on the circle together
$\frac{\sum_{1}+\sum_{2}+\cdots}{N_{1}+N_{2}+\cdots}=\alpha$, because
$\sum_{i}/\sum_{1}=N_{i}/N_{1}=\sigma_{i}$ and
$\frac{\sum_{1}}{N_{1}}\frac{1+\sigma_{2}+\cdots}{1+\sigma_{2}+\cdots}=\frac{\sum_{1}}{N_{1}}=\alpha$.

But more generally $q^{2}$ can be equal to one, for example in the square
lattice, when $N_{\mathrm{ns}}=4$, and calculation of the forth power average
would be coordination circle dependent, similarly as in the three dimensional
case.


The paper is temporarily withheld