A Cascading Mean-Field Approach to the Calculation of Magnetization Fields in Magnetoactive Elastomers

We consider magnetoactive elastomer samples based on the elastic matrix and magnetizable particle inclusions. The application of an external magnetic field to such composite samples causes the magnetization of particles, which start to interact with each other. This interaction is determined by the magnetization field, generated not only by the external magnetic field but also by the magnetic fields arising in the surroundings of interacting particles. Due to the scale invariance of magnetic interactions (O(r−3) in d=3 dimensions), a comprehensive description of the local as well as of the global effects requires a knowledge about the magnetization fields within individual particles and in mesoscopic portions of the composite material. Accordingly, any precise calculation becomes technically infeasible for a specimen comprising billions of particles arranged within macroscopic sample boundaries. Here, we show a way out of this problem by presenting a greatly simplified, but accurate approximation approach for the computation of magnetization fields in the composite samples. Based on the dipole model to magnetic interactions, we introduce the cascading mean-field description of the magnetization field by separating it into three contributions on the micro-, meso-, and macroscale. It is revealed that the contributions are nested into each other, as in the Matryoshka’s toy. Such a description accompanied by an appropriate linearization scheme allows for an efficient and transparent analysis of magnetoactive elastomers under rather general conditions.


Introduction
Magnetoactive elastomers (MAEs), also often denoted as magneto-rheological or magneto-sensitive elastomers, represent a very promising class of field-controllable functional polymer materials. Their magneto-mechanical properties can undergo plenty of different changes and effects like magnetodeformation and anisotropic changes in dynamic and static mechanical moduli depending on strength, orientation, or modulation of the external magnetic field [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. The specific constitution of an MAE sample itself can also be quite versatile, ranging from the usage of magnetically hard or magnetically soft filler particles to mixtures of both types [17][18][19]. In the following we consider the filler particles as magnetizable, i.e., they display no magnetization in the absence of a magnetic field, being produced from an ideal magnetically soft material. Additionally, different compositions of the polymeric material forming the elastic matrix can have substantial influence on the composite behavior under the applied field. Finally, it is well known that not only the bare amount of magnetic filler content, but also its arrangement into isotropic or anisotropic structures as well as the macroscopic form of the sample itself is of great importance for the control of the effective material behavior. Quite a variety of technical implementations and applications have been proposed and implemented so far. For example, MAEs can be used as a working part in actuators and sensors [20][21][22], energy harvesting devices [23][24][25], micro-robots [26] and -pumps [27], in prosthetic and orthotic devices [28], as well as in ophthalmologic magnetic fixators [29][30][31].
where µ 0 is the vacuum permeability and V S is the volume of the sample. The magnetization field M(r) within a particle inclusion located at position r is assumed to be induced by the local magnetic field H(r), whereas the externally applied magnetic field is described by the constant vector H 0 . To be able to calculate the magnetic energy, two basic inputs are necessary: (1) What kind of the particle microstructure is embedded inside the sample volume, and (2) which magnetization function describes the behavior of inclusions. As we mentioned in the introduction, an ideal magnetically soft behavior will be considered in the present work, following the most common experimental situation [30,39,40]. Having this basic information, the task of computing the magnetic energy stays highly nontrivial due to the presence of two unknown fields M and H in Equation (1), which should be calculated self-consistently. In particular, non-linear magnetization behavior, quite relevant in most experimental situations, makes the computation of magnetic energy and thus a prediction of magneto-mechanical properties technically infeasible. This will become clear from the considerations presented in the following session.

The Dipole Model
In the following let N be the total number of magnetizable particles distributed somehow within the macroscopic MAE sample. We assume a constant magnetization field within each of the inclusions, such that magnetic interactions among all the particles are described in terms of dipole fields. This assumption represents the so-called dipole model [33,[41][42][43][44][45][46][47][48][49][50][51]. It is best suited for inclusions with the shape close to a spherical one, because a homogeneously magnetized sphere generates a dipole field in its exterior exactly. Neighboring inclusion particles should also remain at some distance from each other to assure a homogeneous field over the extent of each particle. This minimum center-tocenter distance was found to be around 3r p , with r p being the particle radius [38]. This corresponds to a maximum of about 20-30% of the volume fraction of well-dispersed particles in an MAE sample [33]. However, for any arbitrarily shaped inclusions or for particles coming closer to each other the dipole approach can likewise be motivated as a first approximation. In the most general form the calculation of the local magnetization field within any inclusion a ∈ N of the sample follows from: Here, L(·) denotes the material magnetization function of the inclusions and will be specified in the following. For instance in case of isotropic linear magnetization, it simply reads M a = χH a , with magnetic susceptibility χ. In the framework of the dipole model, the local magnetic field H a is obtained via: The externally applied magnetic field H 0 is assumed to be homogeneous over the extent of the MAE sample. For an arbitrarily-shaped magnetizable particle the tensor parameter ν a represents the self-demagnetization factor of inclusion a. In Equation (3) we defined the tensorial dipole operator between the particles a and b as: 3r ab r ab − r 2 ab I r 5 ab . (4) The vector r ab = r b − r a describes the distance vector a ↔ b, I denotes the second order unit tensor, and v b accounts for the volume of inclusion b.
The fact that the sum in Equation (3) runs over all inclusions in the entire sample is due to the long-range nature of the magnetic interactions, here approximated in terms of the dipole field, Equation (4). Accordingly, even within the (simplifying) dipole approximation the calculation of the magnetic energy in an MAE sample requires a self-consistent solution of a macroscopic number (here 3N in 3 dimensions) of coupled equations, Equations (2) and (3) ∀a ∈ N . We note, that already in the elementary case of linear magnetization (M = χH), and a respective set of linear equations, a straightforward solution for any macroscopic sample (N ∼ O(10 9...12 )) is technically not feasible as it requires the calculation, storage, and inversion of a (3N ×3N )-Matrix.

Introducing the Cascading Mean Field Approach ('Matryoshka' Scheme)
For specially selected and simplified situations the problem of the equivalence of short-and long-range magnetic interactions has been solved in various previous publications [32][33][34]36,42,52]. The main strategy is based on the understanding that beyond a certain distance from any inclusion the precise local particle arrangement is no more resolvable. It appears as a homogeneously smeared continuous distribution at larger distances. Usually this 'critical' resolution range is around 10 times the average nearest neighbor distance [32]. This behavior simply resembles the fact that from a sufficiently far distance any microstructure appears homogeneous [32][33][34][35]. The idea of decoupling the micro-and macro-effect has been outlined previously [32,33] and it is sketched on the left side in Figure 1. We emphasize that in the macroscopic limit such a decomposition is inherently exact. (Right) Formal discretization of sample volume V S into mesoscopic portions V α , α ∈ [1, N]. On such scales any particle microstructure appears a homogeneous continuous distribution.
In a first step we divide the macroscopic MAE sample into N mesoscopic partial volumes V α , α = 1, N. The length scale of such mesoscopic portions is defined as very small compared to the macroscopic sample size, but very large compared to typical distances among neighboring particles, i.e., the local microstructure. The macroscopic limit assures, so to say, that on a mesoscopic scale any sample volume V S may be perfectly resolved upon increasing N to an arbitrary large number, whereas any particle microstructure appears completely diffuse. On the right in Figure 1 a formal portioning of an MAE sample is sketched. We emphasize that each mesoscopic volume V α appears from any other volume V β as a homogeneous material point containing the locally averaged volume fraction φ α of magnetizable material.
In the following we denote the limiting range around any reference particle where the microstructure starts to appear homogeneous as the mesoscopic sphere. One prominent consequence of decoupling the micro-and macro-effect is the capability to provide even analytically representable expressions for the magnetic energy of an entire macroscopic MAE sample. It can be shown how the macroscopic effect, following from the given form of the MAE sample, and the microscopic effect, following from the local particle structure, are intertwined and how this interplay affects the magneto-mechanical properties of MAEs [32,33]. Nevertheless, although the virtual implementation of the mesoscopic sphere is generally valid, in previous works we introduced additional simplifying assumptions in order to obtain a (partially) analytical tractable form. In particular, we want to mention that we were able to obtain self-consistently the locally varying magnetization field, arising due to the particle microstructure, and this was done independent of the actual macroscopic aspect ratio of the sample. Such a decomposed representation of the problem can only be achieved for the special case of a linear magnetization scheme and assuming an ellipsoidal sample form. Whenever one considers a more general, i.e., non-linear, magnetization function in Equation (2) a decomposition is no longer achievable.
In the present work we want to further generalize our existing approach towards less restrictive assumptions and thereby also towards more realistic situations as likely relevant for practical purposes. Adapting the idea of a mesoscopic sphere to the general case, we may rewrite the formal relation in Equation (3) as: The notation introduced in Equation (5) strictly separates micro-structural components, here with lower indices (i, j), from mesoscopic form components, here upper indices (α, β). Accordingly, parameter N limiting the first sum in Equation (5) corresponds to the number of mesoscopic volume elements (see Figure 1) to which the macroscopic sample has been formally discretized. Since ∑ N α=1 V α = V S the total, or macroscopically averaged, volume fraction of magnetizable inclusions in the sample, φ, amounts to: Furthermore, the field quantities follow a hierarchical order, where (·) α i denotes the corresponding field in the particle located mesoscopically at position α within the macroscopic sample and which is microscopically identified as the i − th particle within the local microstructure, see Figure 1. Regularly, in the macroscopic (or thermodynamic) limit the index α is a continuous variable scanning the macroscopic form of the sample. Thus, the sum ∑ N β=1 represents effectively an integral over the sample volume V S . For convenience we keep here its discrete form and note that an explicit numerical evaluation of arbitrarily shaped samples requires a corresponding volume discretization anyways. The tensors g in Equation (5) are adapted accordingly: In Equation (5) we present the relation in a very general form. So we allow each inclusion to be of a different size, v α i , and to have an individual self-demagnetization factor, ν α i . The average local particle density may vary, φ α , and therefore the number of particles representing an actual local microstructure, n α , as well as the micro-structural dipole operator, g α ij , may be different. In principle, one could here even allow for a locally varying external magnetic field, H 0 → H α 0 .
The present notation greatly generalizes our previous studies. Nevertheless, we introduce several appropriate simplifications since we do not account for all possible situations here. In the following the external magnetic field is strictly homogeneous. The particles immersed in the MAE are identical spheres, they are monodisperse in size with v α i = v and ν α i = ν = 1 3 I ; ∀i, α. Furthermore, we consider a homogeneous particle density on a mesoscopic scale, φ α = φ and n α = n ; ∀α. Finally, for simplicity and to keep the derivation possibly short, we constrain to g α ij = g ij ; ∀α. Thus, we consider only one single type of microstructure present in the entire volume of an actual sample. With the aforementioned simplifications we finally come to the following relation: Actually, the cases of varying mesoscopic densities φ α and/or of varying microstructures g α ij in different sections α of the sample appear practically relevant and also quite interesting for an even more comprehensive theoretical description of MAE in external magnetic fields. Perhaps this will be for future consideration. Nevertheless, for a freshly synthesized MAE sample it is reasonable to assume that the particles are similarly distributed all over its volume and thus φ α = φ and g α ij ≈ g ij ; ∀α. Upon applying an external magnetic field, the local field deforming the specimen may vary on macroscopic scales resulting in deflection, buckling, or any kind of asymmetric macroscopic distortion. Accordingly, the microstructure may rearrange differently in the different parts of the sample. In the present work we neglect such effects and assume uniform rearrangements on a macroscopic scale. We also want to note that in the course of applying an external magnetic field the variation of a particle microstructure within an individual sample has been often disregarded in literature [40,53,54].
Still, Equation (9) for a local magnetic field should be calculated self-consistently when inserted in Equation (2). We make use of what we call a Cascading Mean-Field description and implement the following notation: where represents the macroscopic magnetization averaged over all particles in the sample and denotes the mesoscopic magnetization averaged over the particles found in the mesoscopic volume element V α . Note, that the deviations ∆M α and ∆M α i are introduced such that ∑ N α ∆M α = 0 and ∑ n i ∆M α i = 0. The cascading scheme is sketched in Figure 2. It visualizes the three contributions on the micro-, meso-, and macroscale and how the contributions to the magnetization field are nested into each other, similarly to how it is performed in the Matryoshka's toy. Additionally, we introduce new tensorial operators: and We define the average microscopic structure factor G micro and the average macroscopic sample form factor G macro following our notation via: and That way we may formalize the Gs analogous to Equation (10) as: and In addition, from Equation (18), resp. Equation (17), follows that the corresponding sum over α of ∆ G α , resp. over i of ∆ G i , is equal to zero.
Using the notation convention, given by Equations (10)-(17), we rewrite Equation (9): The mean value of Equation (19) over all particles gives the average magnetic field H. Noting that the tensorial dipole operators are symmetric with respect to their indices, i.e., g αβ = g βα and g ij = g ji , the sum 1 Nn ∑ N α ∑ n i of Equation (19) reads: Here, ∆M i S denotes the mean deviation of the individual particle magnetization, ∆M α i , from the mesoscopic magnetization, M α , as averaged over the sample volume: At this point we shortly want to discuss the role of the magnetization function L(·) in Equation (2). Presuming an isotropic linear magnetization behavior would allow us to keep following a mathematically exact derivation. The reason is that then the average magnetization M and the average magnetic field H are related linearly allowing a direct insertion of Equation (20) in the magnetization function, M = χH. Subsequently, it would be possible to decouple to a great extent the (still exact) relations on the different length scales, i.e., achieving a partition into two almost independent sets of self-consistent equations: One for the macroscopic sample form and one for the microscopic particle structure. The case of isotropic linear magnetization will be discussed in detail in Section 2.4.

General Magnetization and Linearization
To proceed from Equation (20), the central difficulty is due to the inequality M = L(H) for a general, i.e., non-linear, magnetization function. In order to profit from the implications of a linear relation we deploy a classical Taylor expansion, which we restrict to first order to yield a subsequent linearization of the magnetization function. This expansion shall be developed around an approximate average magnetic field A, which we define as follows: The approximate A only accounts for the first two terms of the exact H in Equation (20). We motivate the neglect of the last two terms in Equation (20) by the following reason. One immediately notes that these terms represent the microscopic, resp. mesoscopic, corrections to the macroscopic sample averages and involve the solution/knowledge of a huge system of equations, i.e., since they are in fact coupled: O((3N × 3N)(3n × 3n)). In contrast, the approximate solution for M, resp. A, in Equation (22) only involves 3 equations, one for each spatial direction of the global M. From the computational point of view this represents an enormous benefit. To neglect any correction due to alterations ∆ G and/or ∆M represents the leading order approximation in our cascading mean-field approach.
Beyond the leading order we apply a Taylor expansion around A and approximate the magnetization M α i = L(H α i ) in the particle i of mesoscopic domain α as: Taking on both sides in Equation (23) the mesoscopic average over the microstructure, 1 n ∑ n i , we get: Finally, the mean magnetization field in the entire sample is obtained upon taking on both sides in Equation (24) the macroscopic average over the sample form, 1 N ∑ N α . This yields: Here, we make use of the definition ∆M j S introduced in Equation (21). The equations above represent the linearized approximation for a general magnetization function L(·). Note, that in case of an isotropic linear magnetization, M = χH, as well as for anisotropic linear behavior, M = χ · H, the relations are exact.

The Local Magnetization within an Individual Particle
Subtracting both sides of Equation (24) from the corresponding sides of Equation (23), we end up with the following relation for the ∆M α i : The form of Equation (26) immediately implies that within the present linearized approach ∆M α i is proportional to M α for all i. This fact greatly reduces the computational complexity of the problem. Note that the general solution for ∆M α i can be obtained fully independent of the actual mesoscopic position α within a given sample form. We adopt a more compact notation by defining: The tensors κ ij are exclusively determined by the micro-structural arrangement of the particles, i.e., by g ij . Furthermore, following the common notation for linear magnetization scheme, we may construct the prefactor in the r.h.s. of Equation (26) as a generalized susceptibility tensor: and introduce thereby also an effective generalized susceptibility via: Accordingly, we rewrite Equation (26): Since each ∆M α i consists of 3 components, Equation (30) represents a system of 3n linear equations. Traditionally, the solution of such a system may by obtained via a direct procedure, e.g., inverting the problem via Gaussian elimination. Let us introduce some generalized n × n matrices via: Note, that the entries of the generalized matrices are second-order tensors and correspondingly they effectively represent 3n × 3n matrices. Then, the solution may be expressed formally as: Alternatively, the form of Equation (30) also allows for an iterative solution procedure starting with ∆M α i ≈ χ A eff · ∆ G i · M α as the first estimate. In terms of the formalization in Equation (31) this yields: Note, the series in Equation (33) may be obtained upon expanding Equation (31c) in a Taylor series around I = (I) −1 . The iteration scheme is only meaningful if the increasing powers (c) k rapidly diminish and, correspondingly, the series in Equation (33) converges. Indeed, for any distribution of particles this is true for two reasons. Firstly, since particles do not overlap we have r ij ≥ d p , where d p denotes the particle diameter with v = πd 3 p /6, in Equation (8) and thus the entries of g ij , with i = j, as well as those of ∆ G j /n, in Equation (27) are correspondingly small. Secondly, it is important to emphasize that it is the self-demagnetization term ∝ − ν which represents by far the major contribution in the right-hand side of Equation (26). Hence, upon introducing the generalized effective susceptibility χ A eff in Equation (29) and passing over to Equation (30) this self-demagnetization is implicitly, and exactly, taken into account. Furthermore, by definition, see Equation (29), the components of χ A eff are bound to moderate values even though the components of the generalized susceptibility itself, , may become very large. In contrast to Equation (26), it is now χ A eff which is forming the prefactor in Equation (30), resp. in Equation (31a).
It is important to note that the computation of the generalized matrix s, either via some direct method or iteratively, using Equation (33), is completely independent of the mesoscopic position α within the sample. Nevertheless, we also note a fundamental drawback for magnetization functions L(·) beyond the 'true' linear magnetization scheme. Only in this case where M = χH, or general M = χ · H, the here defined χ A , Equation (28), and accordingly also the effective χ A eff , Equation (29), are pure material parameters. They simply become χ A = χ, resp. χ A eff = χ eff . Beyond linear magnetization functions these parameters depend on the approximate average magnetic field A, given by Equation (22). Accordingly, the matrix c, and consequently s in Equation (31) depend on A. Thus, whenever for example the particle content φ is altered, or the average macroscopic sample form factor G macro changes (likewise if G micro changes), or, most prominently, whenever just a different external magnetic field H 0 is applied, the computation, for example in Equation (32), must be repeated. This is not surprising, since within such a highly coupled system as an MAE sample in the applied field, the externally applied field becomes effectively a 'material' or internal parameter itself. To cut the matter short, each H 0 provokes a different material behavior and beyond the linear regime its role can not be decoupled, or resolved, from the other system parameters.
As a reminder, the leading order approximation is defined upon reducing Equations (20)- (22). In a first step beyond this leading order we will neglect the series in the brackets in Equation (33) and assume s ≈ I. With this approximation we may offer up some accuracy but at the same time we yet gain great generality of the approach, where most parameters are effectively decoupled allowing an efficient and transparent analysis of a very wide range of conditions and situations.

The Mesoscopic Magnetization
Let us proceed with an equivalent derivation of the mesoscopic ∆M α . Subtracting both sides of Equation (25) from the corresponding sides of Equation (24) we find: From Equation (32), we know that ∆M α j ∝ M α = M + ∆M α and accordingly ∆M j S ∝ M, with identical prefactors. Hence, the last line in Equation (34) transforms to: where the tensor ∆ G micro is defined as: Analogue to the derivations for ∆M α i we turn to a more compact notation upon defining a tensor: and via introducing here another, quasi mesoscopic, effective susceptibility: We will further see that χ B eff reflects somehow an effective bulk susceptibility of the composite material in the sample. Thus, we may rewrite Equation (34) in similar form as for ∆M α i in Equation (30): Again, we end up with a system of linear equations, where apparently In the current representation the solution may be again obtained via a direct or iterative procedure. Adopting the notation of generalized matrices the solution is formally computed as: Here, the corresponding generalized matrices are defined as: Analogue to Equation (33), the S is expanded in a series: As discussed below, Equation (33), applying our cascading mean-field approach in a first step beyond the leading order we will approximate S ≈ I. Note that the series in Equation (42) includes the factor φ. The volume fraction of magnetic/magnetizable inclusions is usually φ ≤ 0.3 for practically relevant MAEs. Thus, we expect the corrections from the series in Equation (42) to be considerably smaller compared to I.

The Average Magnetization
Substituting the above derived expressions for ∆M α and ∆M α i , resp. ∆M i S , into Equation (25) the average magnetization among all inclusions in the sample, M, becomes: Here, analogue to Equation (36) we defined a tensor ∆ G macro via: Although the relation in Equation (43) effectively represents 'just' a system of 3 nonlinear equations for the vector M, it is not feasible to solve it directly. Note, especially, that as well as s and S not only depend on A, and thus on M, but are also non-trivially convoluted with each other. Before applying the present approach with some general saturating magnetization behavior we first want to turn to the special case of linear magnetization behavior where most of the principal quantities are decoupled.

The Case of Isotropic Linear Magnetization
In case that we can assume a perfect linear magnetization behavior the present derivation benefits from several simplifications. Most prominently, to yield the magnetic energy in the sample it is fully sufficient to compute only the average magnetization, as we will see below. In contrast, for the general non-linear case the local deviations from the sample mean, or to say the actual magnetization in each point / inclusion, need to be known in order to calculate the magnetic energy of the sample according to Equation (1).

General Relations
In the following, we consider the frequently studied case of an isotropic linear magnetization for each inclusion a: Then, the magnetic energy in Equation (1) reduces to: We remind that the total amount/volume of magnetizable material in the sample equals V S φ. The average magnetization among all inclusions, M, needs to be calculated self-consistently. Noting, that now χ A = χ I and since ν = ν I, with ν = 1 3 for spherical inclusions, we have: Accordingly, we may obtain ∆ G micro in Equation (36) completely independent of any macroscopic or external parameters like G macro or H 0 . The tensor ∆ G micro is here, just as it is also the average G micro in a general case, an entire and exclusive parameter of the local microstructure. Thus it can be immediately obtained exactly, or to any accurate precision desired.
From Equation (43) with L(A) = χA or, equivalently, from Equation (20) with M = χH we get the exact relation: Then, we explicitly find: and the magnetic energy in the sample is given upon insertion into Equation (46). This result requires some remarks. Most importantly, all terms in Equation (48), resp. Equation (49), are completely independent of each other, except for ∆ G macro . Whenever an individual parameter of a given sample changes we may adjust that parameter correspondingly without the need to recalculate the other terms. The only coupling in case of a linear magnetization behavior is entering via ∆ G macro . The calculation of this tensor is explicitly dependent on the microstructure in the form of the prefactor χ B eff in Equation (44). For linear isotropic magnetization this effective susceptibility reads: If required, for example in microscopic simulation models, the calculation of the local magnetization field M α i in each inclusion can also be formulated yet exactly in the case of a linear magnetization scheme. Substituting Equations (32) and (40) into Equation (10) we find: Although it is remarkable that the results, Equations (48)-(51), can be expressed in such a compact and yet exact form for linear magnetization behavior we should be aware of the non-trivial coupling entering via Equation (50). Clearly, a different local arrangement of the particles changes the micro-structural contributions G micro and ∆ G micro . Accordingly, χ B eff also needs to be updated, which requires a recalculation of the macroscopical ∆ G macro . Similar implications arise in the exact form for the individual M α i in Equation (51). In addition, in the case of linear magnetization it is therefore worth considering an approximation which allows for a computationally more efficient calculation. Implementing the leading order approximation as introduced around Equation (22) we would simply neglect ∆ G macro and ∆ G micro in Equation (49) yielding an approximate average magnetization field which we denote M (1) . As suggested above, beyond the leading order we apply s ≈ I and S ≈ I in the following. Then, ∆ G micro in Equation (36) is estimated as: Accordingly, we obtain an approximation for χ B eff in Equation (50) and the local magnetization fields M α i in Equation (51) are approximated via: In addition, we can estimate ∆ G macro because with S ≈ I Equation (44) reduces to: Furthermore, Equations (52) and (54) can be inserted in Equation (49) yielding an improved average magnetization. Whether, or how, to approximate χ B eff , resp. ∆ G micro and ∆ G macro , depends on the characteristics of the actual MAE sample, the desired degree of precision, and the available computation resources.

Alternative Formulation of the Role of the Effective Mesoscopic Susceptibility
An alternative way for consideration of the microstructure contribution will become very enlightening here. We consider an isotropic linear magnetization scheme in which Equations (24) and (32) are exact with L(A) = χA and ∇ A L(A) = χ I. With Equation (22) and the definition of ∆ G micro in Equation (36) we can rewrite Equation (24): Using here the definition in Equation (14) as well as that of χ B eff for isotropic linear magnetization in Equation (50), we finally arrive at: This is a quite compact and remarkable result. Although it is exclusively determined for the case of an (isotropic) linear magnetization behavior, it clearly shines a bright light on the interpretation of effective susceptibility χ B eff . Apparently, Equation (56) exactly represents the dipole calculation scheme of the magnetization distribution within a body containing a volume fraction φ of magnetizable material, which obeys itself a bulk magnetization behavior of the composite as reflected by the susceptibility χ B eff . Generally, and depending on the actual microscopic particle distribution, this new effective bulk susceptibility χ B eff is anisotropic, although the original bulk susceptibility of the inclusion material is isotropic χ = χ I. Accordingly, once the χ B eff for a given microstructure, via calculating G micro and ∆ G micro , is known, it is also possible to apply alternative solution approaches like well-established Finite-Element Methods for an anisotropic magnetization behavior in a given sample body. The form of the microstructure determines the effective magnetization behavior on meso-, resp. macroscale, just like it does χ eff = χ/(1 + χν) on a particle scale. For an isotropic [33], as well as for an 'isotropic-like' [32,34], particle distribution the micro-contributions vanish, G micro → 0 and ∆ G micro → 0, and thus χ B eff = χ eff I. An anisotropic microstructure likewise gives rise to an anisotropic χ B eff .
Beside a better understanding of how the different length scales in an MAE may be separated and at the same time also influence each other, the important message here is how alternative, and computationally advantageous, approaches can enter the formalism to obtain a comprehensive, detailed, and accurate, but at the same time effort-or complexityreduced description of the magnetization fields in an MAE sample. The above-derived relations and interpretations concern the case of a linear magnetization behavior. In the following we want to consider the corresponding implications for the general case of a more realistic magnetization scheme.

Cascading Approximation Scheme for Saturating Magnetization
To be able to provide an accurate approximate description of the magnetization fields in the case of non-linear magnetization behavior, it is necessary to adopt a cascading scheme and start with the leading order in the form: For some given composition, i.e., microstructure, sample shape or volume fraction of magnetizable particles, all terms in Equation (57) are computed independently. In the simplest, or leading order, we may assign each particle in the sample this average magnetization M (1) as obtained upon solving just a single 3-dimensional equation, i.e., Equation (57). Beyond the simple leading order (M α i ≈ M (1) , ∀i, α) we consider the following procedure.
From the solution of Equation (57) we not only find a first approximation for the average magnetization among all particles, but also an approximation for the average magnetic field A. Using Equations (28) and (29) we obtain the tensors χ A and χ A eff . Furthermore, one can then proceed to calculate χ B eff via Equations (36) and (38). Analogously to the case of linear magnetization in Section 2.3 we use s ≈ I and S ≈ I. The local magnetization M α i of an individual particle is then evaluated in the next approximation step as: This result is similar to the approximation in the case of a linear magnetization behavior in Equation (53). Note however a small difference, the tensorial representation of χ A eff , since it may become anisotropic in the case of a saturating magnetization behavior.
We want to highlight that the present approximation approach is completely determined if the tensors G α and G i are evaluated for a given sample form and microstructure.
From these, one directly derives the averages G macro and G micro , as well as the deviations ∆ G α and ∆ G i . With Equation (57) the calculation of magnetization fields requires only one single three dimensional self-consistent equation. The rest of the computation are straightforward addition and multiplication functions. In the following we will work out the here developed approximation scheme for some examples to demonstrate its applicability in practice.

Application to Selected Examples
In this section we want to apply the present approach to some example calculations. In the first situation we consider the distribution of a finite number n of particles obeying some specific saturating magnetization behavior under the application of an external field H 0 . We compare the approximation results to the exact computation of the 3n × 3n system of self-consistent non-linear equations. In a second example we derive tensors G macro and ∆ G α for the meso-and macroscopic effects for a sample of cylindrical form. To prove the accuracy of our approach against an exact calculation for an entire macroscopic sample with, at the same time, a fully resolved particle structure on a local scale in 3 dimensions would require a tremendous computational effort for the 'true' or precise description. For our approximation scheme such calculation is straightly feasible upon combining the subsequent example schemes for micro-and macroscale in Equations (57) and (58). We may set up a 2-dimensional model in the future where an exact calculation of ∼ O(10 6−7 ) particles defining a 2D quasi-macroscopic sample is feasible.

Microstructure Calculation for Finite Number of Inclusions
To compare exact results from solving the full non-linear set of equations self-consistently and the approximation from the here developed approach for an actual particle configuration we introduce the following properties. We characterize the magnetization behavior of the particles via an isotropic saturating magnetization of the form: This function represents a Langevin-type behavior, and using the parameters M sat = 868 kA/m and ζ = 0.0218 m/kA it was shown to adequately reproduce the magnetization behavior of micron-sized carbonyl iron particles [39,40]. The function in Equation (59) is plotted in Figure 3 together with its linear approximation, i.e., M = χ H in the limit H → 0. Apparently, the linear approximation is only adequate as far as H < 100 kA/m. Hence, it is crucial to account for saturating effects and consider a more realistic form of the magnetization behavior, i.e., beyond the linear approximation, to reasonably model authentic experimental situations. Nevertheless, because of its great computation benefits and to obtain a first approximate description, the linear magnetization assumption represents a substantial modeling strategy in theory [16,32,33,44,48,49,55].  Upon solving a system of non-linear equations of the form as presented in Equations (2)-(4) the exact magnetization field within the dipole model is found. In the following we consider some selected particle arrangements of a small number n of identical spherical inclusions, with the same particle volume v and magnetization property according to Equation (59). Under the application of an external field H 0 , the induced magnetization field M i in particle i located at r i is given through: This system of 3n × 3n equations we solve via the Newton-Raphson technique [56]. The average of these individual magnetization fields M i over all n particles defines the exact average magnetization: The application of the approximation scheme to a finite number of particles is straightforward. First, we obtain the G i ∀ i ∈ n according to Equation (13) and calculate its average G micro via Equation (15). The leading order approximation Equation (57) with pure microstructure effect is then calculated via the single 3-dimensional self-consistent equation: With ∆ G i = G i − G macro , the magnetization field M i in each particle i beyond the leading order is adopted from Equation (58) and approximated as: The effective susceptibility tensor χ A eff results from Equations (28) and (29) with (1) . In the present case of an isotropic saturating magnetization, the susceptibility tensor χ A is obtained via: We compare the solution of Equation (60) for all i ∈ n to our approximation approach Equation (63). In order to quantify the precision of our scheme, we use the relative absolute error δM i in the following form: and find both, the maximum, max{δM i , i ∈ n}, and the average error among all particles: For example, in the case of a simulation or any other explicit treatment of an MAE sample it is necessary to describe the magnetization fields in the corresponding manyparticle system. Here, we want to consider a system of 100 particles arranged randomly at a volume fraction φ = 0.3 inside some fictitious boundaries assuring that the particles do not overlap each other. In the first case the sample boundaries are specified by a cube, and, thus, it should more or less represent a part of an isotropic microstructure. In the second case we set the boundaries to an elongated cylinder imitating a part of a rather chain-like particle structure. The main axis of the cylinder is aligned with the x-axis. The randomly generated particle configurations in both cases are visualized in Figure 4. In the following the diameter of one spherical particle defines the unit length in the system.
In the first case, see top of Figure 4, we set each side of the boundaries an equal length of l x = l y = l z ≈ 5.59. Accordingly, the volume fraction of magnetizable particles in the cell becomes φ ≈ 0.3. For the arrangement shown in the top of Figure 4 we find the following average microstructure tensor expressed in Cartesian coordinates: Using this G micro in Equation (62), we find an approximation for the average magnetization among all particles, M (1) . The exact magnetization fields are calculated according to Equation (60). In Table 1 we compare precise and approximate calculations. Exemplary, we apply two different external fields H 0 . In the first line in Table 1 we summarize the results for a moderately low H 0 aligned along the x-axis. In the second line, a moderately large field inclined with respect to the cubic cell is applied.
The comparison between exact M and approximate M (1) shows a quite reasonable agreement for both H 0 . In addition to it, both the average magnetization is obtained in good accordance and the individual magnetization field within each particle is found to be rather accurate. On average the deviations of the approximate M (1) i from the precise M i are less than 1%. The maximal error here is found to be just 2.2%. In the Supplemen-tary Materials we provide a full list of all M i and M (1) i (i ∈ [1, 100]) to demonstrate the remarkable accuracy of our approximation scheme. We highlight that this high degree of precision is achieved already in the quite simple form provided via Equations (62) and (63). The computation benefit compared to the full solution in Equation (60)    The comparison between exact M and approximate M (1) shows a quite reasonable agreement for both H 0 . In addition to it, not only the average magnetization is obtained in good accordance, also the individual magnetization field within each particle is found to be rather accurate. On average the deviations of the approximate M (1) i from the precise M i are less than 1%. The maximal error here is found to be just 2.2%. In the supplementary material we provide a full list of all M i and M (1) i (i ∈ [1, 100]) to demonstrate the remarkable accuracy of our approximation scheme. We highlight that this high degree of precision is achieved already in the quite simple form provided via Equation (62) and Equation (63). The computation benefit compared to the full solution in Equation (60) on the other hand is striking as the approximation only requires the self-consistent solution of a single 3-dimensional equation.
Let us consider now the case where the particles are arranged inside a small but elongated cylinder with its main axis aligned along the x-direction. The aspect ratio Γ of length to diameter is chosen to be Γ = 6 and the volume is adjusted such that the n = 100 spherical particles occupy again a volume fraction of φ ≈ 0.3. The particle structure is visualized in the bottom of figure 4. The tensor G micro we found for this case reads: (68)  Table 1. Comparison between exact calculation and approximate results for the 100-particle configuration inside cubic boundaries, see top image in Figure 4, upon applying two different external fields H 0 . The corresponding field components (x, y, z) are grouped in vectorial notation. The field quantities H 0 , M and M (1) are denoted in A/m. Let us consider now the case where the particles are arranged inside a small but elongated cylinder with its main axis aligned along the x-direction. The aspect ratio Γ of length to diameter is chosen to be Γ = 6 and the volume is adjusted such that the n = 100 spherical particles occupy again a volume fraction of φ ≈ 0.3. The particle structure is visualized in the bottom of Figure 4. The tensor G micro we found for this case reads: Apparently, and in contrast to the cubic cell in Equation (67), the tensor in Equation (68) exhibits much larger absolute values in the diagonal elements. This property reveals the rather anisotropic structure compared to the cubic boundaries. Clearly, G micro xx adopts the largest value here because the elongated structure is aligned along x.
Please note the following: Ideally, for an isotropic distribution of particles inside cubic boundaries G micro vanishes. When averaging over several randomly generated particle distributions inside a cubic cell, all entries of G micro in Equation (67) equally diminish. In contrast, averaging over several randomly generated distributions inside elongated cylindrical boundaries result in only the non-diagonal elements in Equation (68) to vanish, whereas the diagonal elements converge towards a finite value of O(10 −2 ). Since G micro is traceless and due to cylindrical symmetry then G micro yy = G micro zz = − 1 2 G micro xx on average. In addition, with increasing particle number n the absolute values of all elements in Equation (67), resp. the non-diagonal elements in Equation (68), diminish. For a single, explicit structure of n = 100 particles we usually find that such 'negligible' elements are of order O(10 −3 ) in maximum absolute value.
Again, upon using Equation (68) in Equation (62) we find an approximation for the average magnetization among all particles, M (1) and the individual M (1) i are found via Equation (63). The precise solution is calculated according to Equation (60). We consider the same two H 0 as in the case of cubic boundaries. In Table 2 we compare the corresponding results. Table 2. Comparison between exact calculation and approximate results for the 100-particle configuration inside elongated cylindrical boundaries, see bottom image in Figure 4, upon applying the same two external fields H 0 as in the cubic case. The corresponding field components (x, y, z) are grouped in vectorial notation. Magnetic and magnetization fields are denoted in A/m. Again, we find very reasonable agreement between the full and the approximate calculation. On average, the absolute deviations δM i are clearly below 1%. The maximum error found for an individual particle is just ∼1.44%. The full list of all M i and M (1) i is provided in the Supplementary Materials. Clearly, upon comparing the results in Table 1 and 2, we note that upon applying identical H 0 the particles arranged in an elongated structure along the x-axis feature considerably larger components of the magnetization field in x-directions as compared to particles arranged inside a cubic cell.

Macroscopic Shape Effect for a Cylindrical Sample
According to our approach the relevant quantity describing the effect on the magnetization field due to the macroscopic shape of a sample is given by the tensorial operator defined in Equation (14). There, for the sake of consistency with the microscopic operator in Equation (13), it is expressed in a discrete form. Formally, Equation (14) represents the discretization of the volume integration over the entire sample body. With Equation (7) in Equation (14) and the number of volume elements, N → ∞ the mesoscopic volume elements become infinitesimal, V β → dV. Thus, turning to a continuous representation α → r we have to calculate the following integral equation: Function Θ(r) denotes the Heavyside function, which is a defined unity for r > 0 and zero otherwise, in order to omit the pole at r = r , i.e., α = β in Equation (14). The integration in Equation (69) is taken over the volume of the cylindrical sample V C and shall be performed for each component of G. In the following we denote the symmetry, or main, axis of the cylinder as the x-direction and the perpendicular circular plane shall be generated by y and z, see Figure 5. The origin of the coordinate system is defined as the center of the cylinder's mass.
Due to cylindrical symmetry of the sample and the form of the integrand in Equation (69), e.g., G is symmetric with tr( G) = 0, one can conclude that G yy = G zz = − 1 2 G xx and G yz = G zy = 0. Furthermore, the two remaining cross-terms G xy = G yx and G xz = G zx must be equivalently centrosymmetric with respect to the cylinder main axis. Thus, expressed in Cartesian coordinates we find: The cross-term we denote here G xρ as it describes the interrelation between the symmetry axis along e x and the lateral direction along e ρ , with ρ = y 2 + z 2 . Finally, due to symmetry considerations G xx and G xρ can only depend on the 'height' position x and the lateral coordinate ρ. The functions G xx (x, ρ) and G xρ (x, ρ) are given in Appendix A.
The tensor G macro is found upon averaging G over the sample volume: Apparently, the cross-terms immediately vanish due to the centrosymmetric form and we find: Here, f macro describes the volume average of G xx according to Equation (71). We denote the height, or length, of the cylindrical sample as L and its cross section diameter as D. Then, for being a scale-free quantity, f macro can only be a single parametric function of the aspect ratio Γ = L D . In Ref. [36] we plotted the dependency f macro (Γ) for a cylinder in comparison to the corresponding results for a spheroid and a rectangular rod. There, we also show that for φ → 1, i.e., hypothetically pure magnetizable material, we obtain the relation 1/3 − f macro = N , with N the demagnetization factor of a bulk magnetic cylinder homogeneously magnetized along its main axis. The function N (Γ) = 1/3 − f macro (Γ) we derive by the present method exactly resembles previously reported findings [58,59]. In addition to the compliance with the demagnetizing factor, our here derived formalism features a great flexibility in predicting the magnetization field inside a composite sample in an approximation description.
Let us shortly consider the approximation scheme introduced in Section 2.5. With some prescribed microstructure, i.e., G micro in the compound MAE and given G macro , e.g., Equation (72) we immediately may predict the average magnetization among all magnetizable inclusions in approximate, but general form via Equation (57). In the most elementary case of a homogeneous-like, or random isotropic, particle structure we have simply G micro = 0 [33,34,52]. Neglecting the microscopic contribution, Equation (57) characterizes primarily the sample shape effect: Just via the solution of this 3-dimensional equation we find a first approximation for the average magnetization field M ≈ M (1) without the necessity to reconsider the full sample if the magnetization property L(H), the applied field H 0 , or the amount φ of inclusion particles is modified. Furthermore, in the case of a homogeneous-like, or random isotropic, particle structure Equation (58) is adopted to read: (1) . (74) Thus, we immediately note that ∆ G(r) = G(r) − G macro , describes the mesoscopic deviations in the magnetization field due to the macroscopic shape effects. In Figures 5 and 6 we display our results for the diagonal component ∆G xx = G xx − f macro for different aspect ratios of the cylinder (Γ = 1.0, Γ = 0.1, and Γ = 5.0). Due to symmetry considerations it is obvious that ∆G xx = ∆G xx (x, ρ) can only depend on the height position x and the radial position ρ. We introduce dimensionless units X = x L and P = ρ D , to represent ∆G xx (X, P) with X ∈ [− 1 2 , 1 2 ] and P ∈ [0, 1 2 ]. The angular coordinate ϕ = arcsin( y ρ ) with ϕ ∈ [0, 2π) is not relevant. Furthermore, since from Equation (69) we note the symmetry x → −x for G xx we also find ∆G xx (−X, P) = ∆G xx (X, P) and display our results for the cut X ∈ [0, 1 2 ] and P ∈ [0, 1 2 ], see left sketch in Figure 5. On the right in Figure 5, ∆G xx (X, P) is plotted via a 2-dimensional color map for a symmetric cylinder, i.e., Γ = 1. It clearly indicates, that by applying an external field H 0 along the x-direction, the magnetization field into the same direction x (accordingly G xx ) is reduced at the top, and correspondingly also at the bottom, of the cylinder compared to the average magnetization field. In contrast, at the lateral, or outer, fringe the magnetization field is enhanced compared to the average magnetization. These variations are due to the specific cylindrical form of the sample.
In Figure 6 the corresponding ∆G xx (X, P) for an oblate (Γ = 0.1) and a prolate (Γ = 5) cylinder are shown. The results are somehow similar to the case Γ = 1, but also a systematic difference is visible. Clearly, in the case of an oblate cylinder the region with reduced magnetization (compared to the average M) is largely grown. For Γ = 0.1 it spans an entire part inside the sample (for all positions with p < 0.39) from the top to the bottom. The reduction effect in that region is relatively small (G xx ≈ −0.1). In contrast, the part in the lateral fringe, where an enhancement effect occurs is reduced in size but the enhancement effect itself is increased instead (G xx ≈ 0.35). In the case of a prolate cylinder, with Γ = 5, the situation is exactly reversed. Here, the enhancement region is clearly widened (spanning through the entire cylinder where X ∈ [−0.34, 0.34]), but quantitatively weakened. And the reduction effect is now of a large absolute magnitude, but concentrated to a small slice on the top (and the bottom) of the cylinder. To illustrate the oblate, resp. prolate, form of the cylinders we stretched the axes.
Regarding the cross-term ∆G xρ (X, P) = G xρ (X, P), we plot the corresponding 2dimensional color map for a symmetric cylinder Γ = 1 on the left in Figure 7. Since this function is antisymmetric with respect to X → −X the full length of the cylinder is displayed here. Obviously, in the largest part of the cylindrical body the cross-term nearly vanishes. Only right at the edges does it adopt quite large absolute values, positive on the top (x ∼ L 2 ) and negative on the bottom (x ∼ − L 2 ). In fact, due to the perfectly shaped edge in the geometrical representation of a cylinder we encounter a discontinuity exactly at ρ = D 2 and x = ± L 2 . Any real body can not obey such mathematically shaped edges and thus we consider such discontinuity as a modeling artifact here. Nevertheless, the discontinuity is fortunately only of a logarithmic order and we expect our function G xρ to be yet reasonable if ρ < D 2 or |x| < L 2 . In order to give an idea about the 'full' centrosymmetric property of the cross-terms in Equation (70) we sketch it in vectorial form G xρ ye y +ze z y 2 +z 2 on the right in Figure 7. The application of an external field H 0 along the positive x-direction imposes, additionally to the magnetization field along x, as well as components in the lateral direction. At the edge on the bottom of the cylinder this lateral magnetization field is oriented 'inwardly' to the main axis. In contrast, on the top of the cylinder it is oriented 'outwardly' to the edge.
For completeness, we plot in Figure 8 the analogue results for ∆G xρ (X, P) in the case of an oblate cylinder, Γ = 0.1 on the left, and a prolate cylinder, Γ = 5 on the right. Apparently, the cross-term is almost identical for variations of the cylinder aspect ratio. Yet one systematic difference may be identified. In an oblate cylinder the major contribution to ∆G xρ (X, P) is found at the edges along the x-axis. In prolate cylinders, the principal part exists at the edges along the radial axis. To illustrate the oblate, resp. prolate, form of the cylinders we stretched the axes.

Conclusions
To understand the behavior of an MAE sample it is vital to capture the full extent of the magnetic interactions, which drive the changes of material properties under the application of an external magnetic field. From fundamental magneto-statics one immediately notes that these interactions are energetically equivalent on different length scales. Consequently, any comprehensive description of MAE samples is sophisticated as several orders of magnitude in length scale must be bridged. In the present work we developed an approximation approach which allows to greatly decouple the short-and long-range effects to the magnetization field in MAE composites. This is achieved upon introducing a Cascading Mean-Field scheme where we split the magnetization of an individual particle into three contributions: The macroscopic average magnetization field among all particles in the sample, the deviations due to the mesoscopic position within the sample volume, and the deviations resulting from the microscopic location with respect to the surrounding particle structure. In order to cover practically relevant situations, the derivation is performed in a possibly general manner. The intrinsic magnetization behavior of the inclusion material, the shape of the macroscopic sample, the actual form of the local microstructure, the overall amount of the inclusion particles, as well as the externally applied field can be adjusted according to the actual conditions without the need to restart the full computation all over again. This flexibility exemplifies the potential of our Cascading Mean-Field approach.
As a demonstration example for the mesoscopically varying magnetization fields due to a non-ellipsoidal shape of a macroscopic sample we considered a cylindrical form. Our results reveal that deviations from the macroscopic average are found most notably near the sample boundaries and strongly depend on the aspect ratio of the cylinder, see Figure 6. In particular, the cross-term ∆G xρ , coupling lateral to parallel components, displays quite large values, up to order unity near the sample edges, see Figure 7. Thus, for a sample with volume fraction φ of magnetizable particles, the relative mesoscopic deviations from the macroscopic magnetization average can reach O(φ). In order to illustrate the application of our approach for a resolved particle microstructure we considered an ensemble of 100 inclusions and compared the approximation result to the full self-consistent solution of the problem. These comparisons, as summarized in Tables 1 and 2, clearly demonstrate the high accuracy of our approach and the relative errors usually range in the order of just ∼1%.
The explicit particle arrangements are usually considered in computer simulations [37,40,[45][46][47][48][49]53,[60][61][62], where it is essential to calculate the magnetization fields in individual particles to account for the corresponding magnetic interactions. The present approach allows for an approximate but very compact and yet accurate estimate of the local fields in individual particles situated in an arbitrary mesoscopic portion of the sample. In this way, we believe that our results represent a substantial progress and will help to develop more refined hybrid models to bridge the gap among the several length scales necessary to characterize MAE samples comprehensively. Currently, our approach is based on the dipole model to account for magnetic interactions. This adaptation represents a commonly applied modeling strategy towards MAEs. The here presented scheme, where different effects on different length scales are considered as nested into each other like in the famous Matryoshka toy, can be extendable in future works to account for a more adequate description of magnetic interactions among closely-located particles where the dipole approach reaches its limits [63]. As a prospective extension of the Matryoshka toy, upon introducing a fourth element which accounts for dipole corrections, it should be additionally helpful that such modifications are bound to very short ranges.
Beyond the computation benefit for application in simulations and/or hybrid models we want to emphasize that the cascading mean-field approach also provides a systematic framework to understand and to study the interplay between different effects directly. This is possible as our scheme provides a compact analytic, or at least semi-analytic, notation describing the leading effects on different length scales. Note, such analysis is achieved through the understanding that modifications in magnetization fields due to particle (re)arrangements, sample deformations, or altered external fields induce changes in the magnetic energy. Correspondingly, investigating the energetic minimum provides systematic insights to the behavior of an MAE. In a previous work [36] we used the leading order approximation of our approach for an efficient combination of theory and experiment to clarify some field-induced effects in MAEs. Furthermore, our scheme allows to study in detail the role of potential particle rearrangements and macroscopic deformations, as well as to consider subsequently the influence of different mechanical couplings and analyze the consequences. Accordingly, the current Cascading Mean-Field approach represents a consequent continuation of our previous work [32][33][34]52] towards a fundamental generalization of the idea to establish a unified theoretical description of the behavior of MAEs under an applied magnetic field. Data Availability Statement: On inquiry, the data presented in this study is available from the authors.

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

Appendix A
The calculation of the 3-dimensional integral in Equation (69) is not straightforward, especially due the necessity of cutting off the pole at r = r . In the following we derive partially analytic closed forms, i.e., we may leave some parts of Equation (69) as numerically derived integrations. There may be yet several more mathematically involved solution methods allowing for representations with more or full analytically closed terms. However for the present example we consider partially numerical forms as sufficient.
In Figure A1 we sketch our decomposition scheme of the integral Equation (69) into three parts. The first part refers to the innermost sub-cylinder, red colored in Figure A1, and accounts for the contribution right above, resp. below, the position r. The second part, blue colored in Figure A1, captures the encapsulating cylinder around the red one where yet a full angular integration, i.e., ϕ ∈ [0, 2π] is viable. Finally, the third part describes the non-shaded rest in Figure A1. With such decomposition we find for G xx : Here, we use dimensionless cylindrical coordinates X = x L and P = ρ D with L the length, or height, of the cylinder and D its diameter. The aspect ratio of the cylinder is then Γ = L D . In the first line in Equation (A1) the analytically closed form for the first and the second part of the decomposition (red-and blue-colored sub-cylinder in Figure A1) are given. In the second line the remaining part is condensed in function F (W, P) which reads here: This 1-dimensional integral can be solved numerically.
For the cross-term G xρ , we find the following form: The function A(V) is defined as: where 3 F 2 denotes the generalized hypergeometric function, E represents the complete elliptic integral of the second kind, and K represents the complete elliptic integral of the first kind.