FORC analysis of magnetically soft microparticles embedded in a polymeric elastic environment

First-order reversal curve (FORC) analysis allows one to investigate composite magnetic materials by decomposing the magnetic response of a whole sample into individual responses of the elementary objects comprising the sample. In this work, we apply this technique to analysing silicone elastomer composites reinforced with ferromagnetic microparticles possessing low intrinsic coercivity. Even though the material of such particles does not demonstrate significant magnetic hysteresis, the soft matrix of the elastomers allows for the translational mobility of the particles and enables their magnetomechanical hysteresis which renders into a wasp-waisted major magnetization loop of the whole sample. It is demonstrated that the FORC diagrams of the composites contain characteristic wing features arising from the collective hysteretic magnetization of the magnetically soft (MS) particles. The influence of the matrix elasticity and particle concentration on the shape of the wing feature is investigated, and an approach to interpreting experimental FORC diagrams of the MS magnetoactive elastomers is proposed. The experimental data are in qualitative agreement with the results of the simulation of the particle magnetization process obtained using a model comprised of two MS particles embedded in an elastic environment. Supplementary material for this article is available online

The unique properties of MAEs are enabled by the coupling between the interacting MS particles and mechanically soft polymeric matrix. At the microscopic scale, the coupling manifests itself in two major phenomena: (i) the soft matrix allows particles to move relative to each other under the action of the magnetostatic force, and (ii) the elastic forces produced when the moving particles deform the surrounding matrix, enable the elastomer to restore its initial shape and bring the particles to their initial relative configuration. These two phenomena additionally lead to another somewhat unexpected property of MAE samples. Even though the MS material of the particles by definition has a miniscule coercivity and, effectively, a single-valued dependence of its magnetization on the applied field, the magnetization curves of MAE composites possessing elastic modulus ∼10 kPa exhibit a distinct waspwaisted hysteresis with descending and ascending branches coinciding at the origin (see figure 1) [16][17][18]. Taking into account the immense importance of the magnetization vs. field dependence for any magnetic materials including MAEs, it is essential to investigate the phenomenon of the MAE magnetic hysteresis thoroughly and find out how it is affected by the concentration of particles, magnetic properties of the particle material, and by the elastic modulus of the polymeric matrix.
In literature, the classical loop of the magnetic hysteresis is usually characterised by the values of saturation magnetization, remanence, intrinsic coercivity, energy product, etc. In the context of the uniquely shaped magnetization loops of MAEs, some of these quantities cease to make any sense, and the hysteresis must be described using a different approach. Instead of introducing additional macroscopic parameters for the composites in which particle interaction plays a crucial role, we suggest the use of the first-order reversal curve (FORC) technique which represents a magnetic hysteresis of a sample as a result of the hysterestic magnetization of many elementary objects (e.g. separate particles, domains) dwelling in the composites under study [19][20][21]. It is required that each of such elementary objects possesses a hysteresis loop describable in terms of its half-width H c and the field-coordinate of its centre H u . The borders of a hysteresis loop are defined by the field values at which rapid changes in magnetization (jumps) take place, or, alternatively, by the field range where two different values of magnetization correspond to the same field value. If the requirements are met, the final result of the FORC analysis will be a two-dimensional contour plot commonly called FORC-diagram showing the distribution of those elementary objects over H c and H u [21][22][23].
In this work, we will investigate the hysteresis of MAEs filled with the MS particles by means of the experimental FORC technique and will compare the measured data with a simple theoretical model where the elementary system is comprised of two MS particles nested in an infinite bulk of a soft elastomer [24].

Fabrication of MAE
The MAE samples used in this work were prepared by admixing two-component silicone rubber Elastosil RT623 (Wacker Chemie AG, Germany) with liquid silicone oil M1000 (GE Bayer Silicones, Germany) and with spherical MS particles made of carbonyl iron with the mean diameter ∼ 4 µm (carbonyl iron powder grade CC, BASF, Germany). The MS microparticle powder was pretreated with the silicon oil to provide enhanced particle compatibility with the matrix. Further process of manufacturing of composite samples consists of four main phases. Initially all prepared components of the composite are thoroughly mixed (I), then the liquid mixture is degassed (II) and poured into appropriate molds (III). The final phase is the polymerization process (IV), which was performed at room temperature. The elastic modulus of the elastomer was customized by varying the ratio between the amount of the silicon rubber and oil. Additionally, reference specimens were made with an epoxy resin matrix in which the particles are immobilised. Each sample used in the magnetic measurements was prepared in the shape of a disc with diameter ∼4.65 mm and height ∼1.1 mm, and its weight was registered with the Ohaus Explorer EX225D/AD semimicro balance. In addition to particle-filled specimens, rodlike samples consisting of a polymer matrix without magnetic particles were also fabricated. Such unfilled samples were used to evaluate the elastic modulus of the matrix itself. The geometrical parameters of the rod-like samples were registered by an electronic caliper. A typical rod-like specimens had the length l in the range of 35-45 mm and the diameter of d ∼ 4.65 mm. A set of discrete values of the particle volume concentration ϕ p was chosen from the range [5,35] vol.%. Samples with ϕ p > 40 vol.% are not stable, deteriorate and fall apart, whereas samples with ϕ p < 5 vol.% typically do not exhibit measurable magnetic hysteresis. The saturation magnetization M s of the particle material has been estimated as 1.5 MA m −1 , whereas its initial internal susceptibility χ 0 was taken equal to 10 5 [25].
The shear moduli of the elastomeric matrices without a magnetic filler were estimated utilizing a quasi-static torsion test on the MCR301 rheometer (Anton Paar GmbH, Austria). For fixation in the rheometer SCF7 clamps, the ends of the samples under study were inserted in circular plugs made of hard polyethylene. The shear modulus G of the rod under torsion is to be estimated from the torque T versus deflection angle ψ as For small strains, the curves T(ψ) are linear and the modulus G was calculated from the corresponding slopes of the curves obtained for each specimen under consideration. Due to the high resolution of the MCR301 rheometer the modulus estimation is as well very accurate (minimum torque = 0.1 µNm, torque resolution <0.1 µNm, resolution of the optical incremental encoder for measurement of the shear strain <1 µrad).
The measuring configuration as well as the advantages of the described method to evaluate shear properties of soft elastomers over plate-plate rheometry are discussed e.g. in [26]. An additional quality control was carried out on the discshaped samples intended for further magnetic measurements. In particular, using computed microtomography (TomoTU µCT, TU Dresden, Germany), the radiographs of the samples were measured in order to detect defects such as air inclusions and non-uniform distribution of the magnetic micropowder inside the bulk of the samples. Tracking and identification of individual microparticles used in this study was not possible, however, even a simple visual assessment of the grey values (brightness) distribution over the sample volume enables conclusions to be drawn about the homogeneity of the spatial powder distribution.
The parameters of the samples used in this work are presented in table 1.

Magnetometry setup
Magnetic behaviour of the MAE samples was investigated by means of the LakeShore 7407 s vibrating sample magnetometer (Lake Shore Cryotronics, Inc. USA). The MAE disc-shaped samples were oriented with the flat surface perpendicular to the external magnetic field generated by the electromagnet of the VSM as demonstrated in figure 2. To calibrate the VSM, the standard nickel sphere with a diameter of 3 mm and magnetic moment of m = 6.92 emu reached at the external field H = 400 kA m −1 (5 kOe) was utilized. All observations have been conducted at room temperature T = 23 • C. The sample averaging period used for measuring the magnetic moment was set to 1 s pt −1 . This corresponds to the noise floor of 175 nemu. The accuracy of the external field measurements of the VSM in the used field range was 8 A m −1 (0.1 Oe).
Prior to performing the FORC measurements, each specimen was trained using a cyclic magnetization process at the full range of the used magnetic field strength (between −1200 and 1200 kA m −1 ). The training was conducted in order to obtain reproducible magnetization curves regardless of the changes in the field polarity.

Magnetic hysteresis
Several experimental and theoretical works [16][17][18] have demonstrated that the magnetic hysteresis shown in figure 1 is due to the hysteretic rearrangement of the particles inside MAE samples, which in its turn is due to the interplay between the magnetostatic and elastic forces at the microscopic scale of the MAEs under magnetization. It has been observed, that the softer the elastomeric matrix, the more pronounced the hysteresis, in other words, the ascending and descending branches of the hysteresis loops are separated with a wider gap between them and are more easily distinguished visually on the magnetization graphs of the mechanically softer samples rather than of the harder ones. Two MS particles is the simplest MS structure able to magnetize in a hysteretic way. Letters L and θ denote the length of the line segment connecting the centres of the particles and the acute angle between the segment and the field vector, respectively. Notations L 0 and θ 0 are used for the same quantities in the absence of the field (H 0 = 0). Reproduced from [24]. © IOP Publishing Ltd. All rights reserved.
The simplest structure of the MS particles able to exhibit such a hysteresis consists of only a pair of spherical particles with their centres being initially at a distance L 0 apart (see figure 3). To describe magnetization of such a system, we use the model proposed in [24] and give its summary here. According to the model, the magnetization of each particle obeys the Fröhlich-Kennelly law [27], the interaction of the particles is calculated by the dipole-dipole approximation and the elastomer is described as a linear isotropic non-compressible elastic continuum. The state of such a two-particle system (2psystem) at any particular applied field H 0 is found via the minimization of the system's potential energy W = W(H 0 , U 1 , U 2 ) comprised of the magnetostatic energy of the particles W m [24] and the elastic energy of the deformed matrix W e [28]: where µ 0 is the vacuum permeability; V stands for the particle volume; χ 0 and M s are the initial internal susceptibility and the saturation magnetization of the particle material, respectively; H is the field inside the particles; G denotes the shear modulus of the elastomer matrix; U i is the displacement of the ith particle; and a is their radius. Apart from the particle size, χ 0 , M s and G do not change during the magnetization process of a given 2p-system. By definition, the volume concentration of the particles ϕ p and the mean initial distance between themL 0 are connected asL An example of the magnetization loop of a 2psystem is shown in figure 4 (see also the GIF animation in supplementary material available online at stacks. iop.org/JPD/55/155001/mmedia). In a uniform quasi-static external field H 0 , the magnetostatic force acting on the particles is due to the quasi-static stray field generated by the polarized particles themselves whose magnetization grows with the field and turns into zero when the field is switched off. Depending on the relative location of the particles, i.e. angle θ 0 , they can either repel from or attract to each other. In general, this interaction is of a many-body nature and, therefore, is quite complicated to calculate. Nevertheless, the comprehension of the magnetomechanical hysteresis can be achieved by means of some basic assumptions. For example, it can be expected that the magnetostatic force between two particles grows faster with the change in the distance between them than the opposing it elastic stress in the deformed matrix [25,29]. If the displacement of the particles U i is small, the state of the system is defined by the equilibrium between the magnetic and elastic forces.
In the case when the magnetized particles attract to each other, after the distance between them reaches a certain critical value at a high enough field H ↑ , the magnetostatic attraction overcomes the elastic repulsion and the particles start collapsing, while their magnetization surges because it gets affected by the stray fields which are stronger in the vicinity of the particles. The collapse stops at the distance where the elastomer between the particle surfaces cannot be deformed any further or when the particles start touching each other. With the further increase in the field, the magnetization of the particles also grows until it reaches the saturation M s . Correspondingly, when the field is reduced, the magnetization of the particles and the attraction between them diminishes. At a certain field H ↓ < H ↑ , the magnetostatic force becomes lower then the elastic force striving to bring the particles to their initial positions, and the particles come apart. Since the stray fields fall down with the distance from the particles, the particle magnetization plummets at that moment. The collapse and the separation of the particles happen at different field values H ↑ ̸ = H ↓ , hence the hysteresis of their relative distance and of their combined magnetization [17,24].
Due to the fact that the particles are magnetically isotropic and do not possess any remanence, the magnetization curve is symmetrical relative to the origin, and there are four jumps of magnetization located in the first and third quadrants of figure 4. Because of that symmetry, the hysteresis of a 2psystem can be characterised only by the two positive values of the fields H ↑ and H ↓ at which magnetization jumps in the first quadrant. Alternatively, one can use the half-width of the loop H c and the x-coordinate of its centre H u which can be expressed via H ↑ and H ↓ as By introducing a distribution of particles overL 0 and θ 0 one can model magnetization of a whole MAE sample. The magnetization curve of such an ensemble of 2p-systems closely resembles that shown in figure 1. In the remainder, all the calculation is done for MAEs filled with carbonyl particles whose M s and χ 0 are given in section 2.
This simplified representation of the complex processes taking place in MAEs during magnetization serves the purpose of speeding up the modelling of the FORC diagrams being discussed in section 4.

FORC analysis
Information about the magnetic hysteresis of separate elementary objects constituting a system under investigation can be obtained by means of the FORC analysis. This technique was first introduced by Mayergoyz as a supplementary experimental tool enabling one to extract distribution of hysterons for describing the hysteresis of a sample in accordance with the Preisach model [19]. A hysteron is a model of a ferromagnetic particle which possesses a rectangular hysteresis loop with the saturation magnetization M s (see figure 5). The rectangular loop can be characterized by the fields H ↑ and H ↓ or H c , H u , which are analogous to those introduced in expression (5).
The process of acquiring data for the FORC analysis involves measuring a set of the FORCs (hence the name FORC) as schematically illustrated in figure 6. The measurement of each such curve is preceded by ramping the applied field H 0 to some value H max 0 high enough to saturate the sample. After that, the field is decreased to H r which is a defining parameter of a specific FORC because it is the starting value of the field at which the measurement of that curve begins. The data points of a single curve are the values of magnetization M(H r , H f ) measured at different H f taking on the values in range [H r , H max 0 ] with a certain step ∆H 0 . The values of H r of any two consecutively measured curves are also separated by ∆H 0 and lie in (−H max 0 , H max 0 ). After measuring (or simulating) all the curves, the so-called FORC distribution is calculated:  Process of acquiring data for the FORC analysis. The blue lines denote the curves that have been already measured (simulated); the red curve is being currently measured. All the curves for which Hr ≤ 0 coincide with the bold blue branch of hysteresis at H f ≥ 0. The dashed curve denotes the major hysteresis loop, exaggerated for clarity. Reproduced from [24]. © IOP Publishing Ltd. All rights reserved.
There exist several interpretations of this quantity [21][22][23]. The geometrical interpretation would be that ρ(H r , H f ) shows how the slope of two neighbouring curves at the same field H f changes while going from the curve defined by H r to the curve defined by H r + ∆H 0 . If the slopes of the curves are equal everywhere except for the field values where the magnetization jumps, ρ(H r , H f ) will be non-zero only at H r and H f corresponding to those jumps. In the case of a rectangular hysteron, ρ(H r , H f ) is non-zero only at H ↓ and H ↑ where the slope of the curves changes from the horizontal orientation to the vertical one.
Mayergoyz demonstrated that the FORC distribution built in coordinates actually represents the distribution of hysterons constituting a sample under investigation over fields H c and H u [19]. Even though expressions (5) and (7) look similar, the meaning of the variables on their right-hand sides is different. The values of H ↑ and H ↓ reflect the shape of individual model hysterons, whereas H f and H r are the field values directly measured during the FORC experiment. From section 3, it can be concluded that there is a strong similarity between the hysteresis loops of a hysteron and a 2p-system. This fact encourages one to use the FORC technique for describing hysteretic magnetomechanics of a MAE sample by representing it as the result of the magnetization of all the 2p-systems in that sample. However, the slopes of the two branches of the 2p-system loop are different, and this is the major source of discrepancy between the FORC distribution of a 2p-system and that of a hysteron. This was explored in detail in [24].
The standard way to study the FORC distributions of different samples is by building and comparing the contour plots of ρ(H c , H u ) or ρ(H r , H f ). In the case of the graph of ρ(H c , H u ), the horizontal and vertical coordinates of positive zones gives an idea of the values of half-width and location of the centres of the elementary system loops, respectively. The negative regions of the diagrams bear additional information about relation between neighbouring FORCs and of the sample hysteresis as a whole [30][31][32]. For instance, figure 7(a) shows the FORC diagram of the 2p-system whose hysteresis is illustrated in figure 4.
Since the landscape of the FORC-diagrams in either (H f , H r ) or (H c , H u ) planes is by and large determined by the field values where the jumps of magnetization take place, it is therefore convenient to sketch the shape of a modelled diagram not by calculating and drawing the full FORCdistribution (see expression (6)) but by plotting markers in the same coordinate system at the points with coordinates corresponding to the switching fields of the loops of the modelled 2p-systems. For instance, the peaks in figure 7(a) and markers in figure 7(b) are located at H c = 98 kA m −1 , H u = ±220 kA m −1 , that is, at the coordinates calculated using expression (5) and the values of H ↑ = 318 kA m −1 and H ↓ = 122 kA m −1 obtained from the loop in figure 4.
Previously it was shown that MAEs can exhibit the socalled training effect manifesting itself in the discrepancy between the shapes of the first several consecutively measured magnetization loops [33]. Even though this phenomenon can be investigated by the FORC analysis, the FORC diagrams of the MAEs studied in this work were measured after the samples were trained beforehand by magnetizing them in ten field cycles between H max 0 and −H max 0 . It ensured that, while the FORCs were measured, the particles kept returning to the same spacial locations every time the field was equal to zero, and that no irreversible change in the structure of the elastomer took place.

Measured and modelled FORC diagrams
The experimental data presented in this section were measured in accordance with the process discussed in section 4 using the field step ∆H 0 ≈ 3.8 kA m −1 . Such a step is small enough to unveil the important features of the diagrams, and is large enough to keep the time of the FORC measurement decently short. The FORC distribution was calculated using the locally weighted regression algorithm (LOESS) described in [34,35]. The fraction of the neighbouring points used for the smoothing of every value of M was found using the smoothing factor optimization approach as suggested in [36] and was equal to 0.32% that corresponds to the effective smoothing factor of the truncated algorithm [20] SF = 10 (441 points).
The FORC diagram of the powder, which has been used in this work, immobilised in epoxy is shown in figure 8(a). The particles are made of carbonyl iron that is a rather MS but nonetheless ferromagnetic material yielding a high positive ridge near H c = 0 axis. This is the only major feature of the diagram of the immovable particles. In contrast, the diagram of a MAE sample shown in figure 8(b) contains additional positivevalued areas extending from the origin towards the first and fourth quadrants and mostly located inside the area limited by the vertical H u axis and two diagonals H c = H u and H c = −H u . The visual similarity of that feature with the wings of an insect allows us to refer to it as 'wing feature' in the remainder.
Since neither epoxy nor polymeric elastomer introduces any significant signal to the magnetic measurements, one can conclude by comparing the two diagrams, that the wing feature must be caused by the magnetomechanical hysteresis of the MS particles, and its size can serve as an estimation of the particle mobility in any given matrix.
Taking into account that the same sort of the magnetic powder was used in all the samples presented in this work, there are only two major factors able to significantly affect the shape of the wing feature: (i) the elasticity of the elastomer matrix and (ii) the concentration of the particles inside a sample. In terms of the model, those two factors are represented by the shear modulus G and the average initial interparticle distanceL 0 , respectively. Figure 9 illustrates the differences between the FORC diagrams of MAEs possessing different matrix shear moduli. Frames 9(a)-(c) contain experimental diagrams, whereas frame 9(d) demonstrates the expected location of the positive features of the diagrams calculated using 2p-systems at three different G and whose orientation angles θ 0 are evenly spaced between 0 • and 90 • ; the initial distance between the pairs of the particlesL 0 was calculated using expression (4). All the points of the same colour represent the locus of only the central ridge of the corresponding wing feature. The unique wing shape of the experimental features is due to the fact that the magnetic properties of the particles as well as the mechanical properties of the matrix in the vicinity of those particles can significantly vary in real samples.
The experimental diagrams differ from each other in three major ways. First, the overall vividness and the magnitude, in particular, of the positive features in the diagrams decrease with growing G. This phenomenon indicates that the particles get less mobile as the matrices get harder, and the hysteresis becomes less pronounced than that shown in figure 1. Second and third, the wing features get more elongated and move towards axis H u with the increase in G, implying that the hysteresis of 2p-systems becomes narrower (lower H c ) and it takes higher fields to force the particles to collapse in the first instance (higher H u ). The presence of the wing features in all the experimental FORC diagrams evidently shows that even in our least compliant MAE sample s1, some of the iron particles experience magnetomechanical hysteresis. In contrast, it would be difficult to make such a conclusion by only examining the traditional magnetization curve of sample s1, inasmuch as it would be impossible to tell the difference between the visually coinciding ascending and descending branches of the hysteresis.
In a qualitative way, the last two noted differences between the measured FORC diagrams are supported by the scatter plot of figure 9(d), where each graph was built by simulating 19 2p-systems, with θ 0 evenly spaced in range [0 • , 90 • ]. One can see that the angle between the ridges and axis H u decreases, whereas the average H u coordinate of the points representing samples with a certain G increases with the growth in the matrix shear modulus. According to the model, in matrices with G > 10 kPa and particle concentration ϕ p ≈ 5 vol. %, the magnetomechanical hysteresis of 2p-systems is impossible: the particles are located too far from each other and the matrix is not compliant enough for them to collapse. That is why the graphs in figure 9(d) were built only for G ≤ 10 kPa, and the number of markers in each of them is different.
Consequently, the vanishing but still remaining parts of the wing features in figures 9(b) and (c) indicate that some particles in the corresponding samples were located closer to each other than atL 0 estimated by expression (4), or the local matrix elasticity in the vicinity of particular particles differed from that measured on a macroscopic sample of a pure matrix. This result is not surprising: it is expected that the concentration of the particle and the mechanical properties of the matrix are not absolutely uniform in the bulk of real MAE samples. The essential message here is that in contrast to the conventional magnetization curves, the FORC diagrams are indeed able to unveil the details of magnetization of separate elementary systems comprising a sample.
The coordinates of the wing-features in the FORC diagrams and of the points in the scatter plot gives an idea how the magnetization curves of separate 2p-systems look like. In softer matrices, the 2p-systems exhibit hysteresis of different width with the loops close to the origin (high range of H c , low H u ), whereas in harder matrices their hysteresis tends to be narrow and shifted away from the origin (narrow range of H c , high H u ). Figure 10 illustrates how the concentration of particles ϕ p affects the FORC diagrams of MAEs possessing similar matrices. As in figure 9, frames 10(a)-(c) contain experimental diagrams and frame 10(d) illustrates where the wing feature would appear on modelled diagrams. The simulation was run using three groups of 19 evenly oriented 2p-systems each, at the same G but different values ofL 0 (i.e. ϕ p ).
From comparing the experimental diagrams, it is seen that with growing ϕ p the wing features become more elongated in the direction away from the origin and move closer to the H u axis as if the shear modulus of the matrices was also growing, similar to what was happening in figures 9(a)-(c). This can be explained by the dependence of the sample elastic modulus on the particles concentration [37][38][39][40], as well as by the so called magnetorheological effect [41] when the effective elastic modulus of an MAE sample grows with the magnetic field: the higher the particle concentration, the more pronounced the phenomenon. Because of the increase in the particle concentration, the magnetization of the samples also increases-this, in its turn, affects the height of the features in the FORC diagrams, which is in full agreement with expression (6). The increasing intensity of the wing features leads to the increase of their area if the diagrams are built using the same colourbar limits.
In the case of the growing concentration, the model can provide a general idea of how the wing-features could be oriented in FORC diagrams. Since the 2p-systems do not allow for the magnetorheological effect, the angle between the ridges of the scatter plots in figure 10(d) and axis H u does not vary much. The number of the markers comprising the scatter plots grows with ϕ p reflecting the fact that when the particles are located closer to each other, there are more orientations of the particles (represented by the angle θ in the model) at which they can magnetize in a hysteric way. Similarly to what we discussed regarding figure 9(d), the difference in number of points between the graphs is caused by the fact that not all of the 19 simulated 2p-systems magnetize hysteretically.
The most obvious difference between the model and the experimental data in figures 9 and 10 is the discrepancy between the exact coordinates of the scatter points and the coordinates of the wing features. This can certainly be explained by the general simplicity of the 2p-system, contrasting with the rigorous approach of the models with many particles and treating the matrix deformation in a more thorough way. However, the calculation of the MAE magnetization model based on the 2p-systems is not at all timeconsuming and, therefore, it allows one to quickly obtain the approximate shape and location of the wing features and to acquire an idea of how the FORC diagrams of MAE samples should behave under a change in sample parameters. Figure 11 demonstrates both qualitative and to a certain extent quantitative agreement between the model and experiment for particles concentration values ϕ p = 20 and 35 vol. % (the case of ϕ p = 5 vol % has already been discussed in the context of figure 9). In the figure, the locations of the central ridges of the wing features calculated using the 2p-system are juxtaposed with the experimental diagrams of the samples with different G. The scatter plots were obtained by simulating the magnetization of 100 2p-systems randomly oriented in space, that is, whose distribution of orientation angles θ 0 corresponds to the continuous uniform distribution of the solid angle.

Conclusion
In this work, for the first time, MAE samples filled with the MS particles made of carbonyl iron have been studied by means of the FORC analysis. It has been demonstrated that the FORC diagrams of mechanically soft MAEs contain a distinct positive-valued feature visually resembling wings of an insect. The presence of the wing feature in the FORC diagrams indicates that the magnetic particles experience magnetomechanical hysteresis during magnetization. This feature is not clearly observable in the diagrams of the samples with G > 100 kPa, and that makes the FORC analysis especially useful in testing the local elasticity of the matrix: if the mechanical properties of the matrix are non-uniform, the particles will move more easily and contribute to the strength of the wing feature in the corresponding FORC diagrams. In contrast, it might be difficult to observe any hysteresis at all on the major magnetization curves of such samples because the gap between the ascending and descending branches can be extremely difficult to notice. Higher concentration of the particles in samples produce the FORC distribution (expression (6)) of higher magnitude and affect the shape of the wing features by indirectly increasing the elastic modulus of MAE samples.
The experimental diagrams have been compared with the results of calculations based on the 2p-system model comprised of two MS particles separated by an isotropic, linearly elastic and non-compressible media. Qualitatively, the model is in agreement with the experiment and gives clear understanding of how the wing features of samples possessing similar amount of magnetic powder but various matrix properties must differ from each other. Even taking into account the simplicity of the 2p-system representation, any significant discrepancy between the modelled and experimental diagrams can serve as an indication of the local elasticity of the matrix in the vicinity of the particles being different from the macroscopic one, and of the non-uniformity of the particle concentration in the samples under study.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.