Three-dimensional distributions of demagnetization factors *N*_{d} within assemblies of magnetic nanoparticles have been investigated along the axes of cuboid containing vessels. From the results of a numerical polar-based model, a significant skew toward high values in the number distribution is observed and often overlooked by the assumed uniformity of the conventional analytical approach. To enable comparison with experiment, new transverse susceptibility techniques have been developed, which are also applicable to superparamagnetic assemblies that do not have the magnetization features normally required using conventional methods. Applying the two techniques to a system of ∼13 nm magnetite (Fe_{3}O_{4}) particles resulted in the difference between the in-plane and out-of-plane *N*_{d} factors of (0.21 ± 0.03) and (0.201 ± 0.009), respectively, which shows closest agreement with the simulation value for the mode of (0.19 ± 0.02). The median and mean results of the model move further away from the experimental result, yielding values of (0.17 ± 0.02) and (0.16 ± 0.02), respectively, which is consistent with the skewed distributions observed here. In all cases, the sum of the *N*_{d} factors from each orthogonal axis was equal to 1, giving further confidence in the model. The new methods allow measurements on the superparamagnetic systems often found at this scale, and the agreement with the model means that the spatial distribution of *N*_{d} factors may now be taken into account in studies on any nanoscale material that considers the whole structure as a distribution of magnetic elements.

## I. INTRODUCTION

The study of magnetic nanostructures, including discrete and fixed nano-islands and assemblies of magnetic nanoparticles (MNPs), continues to be an expanding and fundamentally important field of interest. Potential and real applications cover a wide range, which includes such examples as magnetic data storage,^{1,2} collections of skyrmions,^{3} biomedical treatments utilizing magnetic hyperthermia [heating of superparamagnetic iron oxide nanoparticles (SPIONs) under an ac field],^{4,5} or for the targeted release of drugs,^{6} and the importance of particle dipolar interactions in forming chains in such systems,^{7–9} through to room-temperature magnetic refrigeration.^{10}

^{11}The demagnetizing field

*H*_{d}is proportional to the magnetization

**that produces it with the constant of proportionality being the demagnetizing factor**

*M**N*

_{d}. The negative sign indicates that

*H*_{d}is in the opposite direction to

**,**

*M*The demagnetization factors of a sample depend on the sample’s shape. For an ellipsoid in a uniform applied field, the geometry is such that the demagnetization field is constant throughout the sample. The three orthogonal axes’ demagnetization factors are constant, and their sum is equal to 1:^{12} for a spherical sample, the demagnetization factor in each of the three axes is ⅓. For other shapes, the demagnetization fields are not constant throughout the sample and average values are typically quoted. For a homogeneous solid cuboid sample, the demagnetization factors are estimated with a dimension-dependent analytical expression.^{13}

*N*

_{z}can be estimated using the following formula:

^{14,15}

*f*defined as the volume packing fraction of the powder and

*D*

_{z}the shape demagnetization factor of the sample-containing vessel in the given direction. At full packing fraction, this equation reduces to the shape factor, and at very low packing, it tends to ⅓, and the value of the constituent spherical particles as the interactions between them tends to zero. However, in the case of interest in assemblies of

*nano*particles at intermediate packing fractions (as found in a typical sample vessel), there is little work in the literature to test the validity and/or accuracy of the estimate given in Eq. (2).

Earlier modeling work^{16} predicted that *N*_{z} would reduce as a function of decreasing thickness in thin-film assemblies of MNPs and hence the conventional approximation to an infinite sheet would no longer be valid and would need to be considered as technologies advance ever-more into the nanoscale. More recently, the effect of the packing fraction and size distribution was numerically modeled for multi-domain *macro-sized* particles^{17} and predicted close agreement with Eq. (2). A subsequent experimental study^{18} on ferrimagnetic MNPs showed just how important it is to ensure *N*_{z} is taken into account. Their study of susceptibility curves as a function of temperature showed that features described by previous researchers were, in fact, artifacts. By calculating *N*_{z} using known *f* and *D*_{z} factors in their cylindrical samples, the artifacts were removed and the curves converged for in-plane and out-of-plane measurements. For further insight into particulate systems, measurement of *N*_{z} in itself is first required for comparison with the predicted results of Eq. (2) along each axis of the containing vessel.

In the work reported here, we develop two new experimental methods that allow the measurement of demagnetization factors in randomly oriented MNPs to be determined from AC transverse susceptibility curves as a function of an applied DC field. As these do not contain the anisotropy peaks of an oriented (textured) sample, the conventional method cannot be used. Furthermore, if the MNPs are also superparamagnetic, and hence with no hysteresis in their $MH$ magnetization curves, the alternative and less accurate “loop-closure” point method^{19} is also not applicable.

The results on randomly oriented magnetite (Fe_{3}O_{4}) nanoparticles that we pack into rectangular prisms, with the appropriate cuboid shape demagnetization factor and measured packing fraction applied, are then compared with Eq. (2) and those of a polar-based computational model we develop based on a lattice of single domain particles. From the subsequent number and spatial distributions of the demagnetization factors, we show that the experimental results are closest to those of the mode, being consistent with the significant skew observed, and so provide insight into the structure of the magnetic interactions within the assembly at this increasingly relevant scale.

## II. METHOD

### A. The model

The dipolar model is based on the work by Bissell and Cookson *et al.*, who were examining demagnetization factors in particulate recording media.^{16,20} Further details are given in the Appendix.

Our model consists of spherical particles that are set in regular rows and columns (simple cubic packing) to form the cuboid shape of the containing vessel. The size of the cuboid is defined by the number of particles in each of the orthogonal axes. For each particle in turn, the model calculates the particle’s self-demagnetization field, caused by free poles at the particle’s surface, and the combined field created by all the other constituent particles (interparticle interactions). Initial testing of the model was carried out with three axes of equal length containing 101 particles each: a perfect cube containing just over a million particles. Adjusting the size of each “length” of the cuboid allowed the aspect ratio of the cuboid to be changed as required. By treating the cuboid as the containing vessel of the powder, the container shape demagnetization factor *D*_{z} can then be found using the analytical expression derived by Aharoni.^{13}

The packing fraction of the powder *f* is set by altering the distance that the particles are apart within the sample. The same separation distance is used in all three orthogonal axes. Since the particles are in a regular aligned lattice, there is an upper limit on the packing fraction: this occurs when the particles are just touching, giving a maximum packing fraction of ∼0.52. The maximum possible packing fraction for spheres of the same dimension in any arrangement is ∼0.74,^{21} cubic close packing.

The model uses an assumption that each particle is fully magnetized in the direction of a saturating applied external magnetic field. Equation (2) states that the average demagnetization factor for a given sample consisting of packed spherical magnetic spheres is constant: it is independent of the magnetized state of the sample. By measuring or modeling the demagnetization factor at any level of magnetization, it should yield the same value. For the modeling work, the easiest state to work with is one in which the particles’ moments are fully aligned in the direction of the applied field.

The demagnetizing field that any given particle is exposed to is caused by the free poles on its own surface and on the surfaces of all the other particles within the cuboid. The model calculates this demagnetization field in the direction of the applied field at the center of each particle in two steps. First, using numerical integration, it calculates the demagnetization field that a given particle itself would generate at its own center. It then, by repeating the numerical integration, calculates the contributions that the surfaces of all the other particles will have at the center of the given particle. Further details are given in the Appendix.

### B. Transverse susceptibility curves

Transverse susceptibility, the initial susceptibility of a material measured in a transverse direction to an applied magnetic field, was originally conceived by Gans in 1909.^{22} In 1957, Aharoni *et al.*^{23} calculated the transverse susceptibility as a function of applied field for a system of identical Stonar–Wohlfarth particles^{24} with a random distribution of easy axis orientations; this function contained distinct features at the particle’s coercivity and anisotropy field values. These features were confirmed experimentally in 1987 by Pareti and Turilli^{25} but with the restriction that it was only valid for particles whose volume was below that for multiple domain formation. Samples of magnetite spherical particles, typical diameter of 13 nm, packed into a rectangular glass tube of internal dimensions: thickness (0.40 ± 0.04), width (4.0 ± 0.4), and length (10.3 ± 0.1) mm, were used in this study. The particles were synthesized using an established method. Briefly, a solution of ferric chloride and ferrous chloride (2:1 ratio) was mixed with an ammonia solution, and the reaction mixture was stirred for 1 h. The magnetite (Fe_{3}O_{4}) nanoparticles produced were washed three times with deionized water and then oven-dried for 24 h at 70° C before use. Further details of the methodology used to produce the nanoparticles can be found elsewhere.^{26} Having control of the process ensured a predominantly superparamagnetic assembly was obtained, but with the small hysteresis required to test the validity of the new transverse susceptibility method developed in Sec. II C.

The established method for determination of a sample’s demagnetization factors requires two transverse susceptibility measurements across the width and thickness, respectively, as a function of an applied DC field.^{19,27} In both cases, a small AC field is applied transverse to the DC field direction along the sample length as shown in Fig. 1. This perturbation of the sample’s magnetization by the AC field, the linear transverse susceptibility, *χ*_{t}, is then measured in the AC field direction using pickup coils.

*H*

_{k}. The sample is immersed in an

*effective*field due to the opposing demagnetization field of Eq. (1) in any given direction, and so the peaks will be shifted along the DC field axis accordingly. Hence, using Eq. (1), the difference in the position of the anisotropy peaks between the two measurements is used as a measure of the difference in the demagnetization fields between the two orientations by

*z*and

*x*representing the directions across the thickness and width, respectively. The common magnetization point,

*M*, is found by cross-referencing the field strengths with standard magnetization curves in the two orientations, and then, the difference in the demagnetization factors may be determined. In summary,

*H*

_{z}and

*H*

_{x}are the anisotropy fields (applied field by the susceptometer) for the same sample measured twice at two different orientations on the susceptometer. The magnetization,

*M*, is measured using a vibrating sample magnetometer, duplicating one of the sample orientations and respective applied field settings from the susceptometer.

### C. Development of the transverse susceptibility method

The presence of an anisotropy peak in the susceptometer output is dependent on the orientation of the sample’s particles, with randomly packed powders commonly suppressing this peak.^{28} This means that another feature needs to be identified that can then be used as a “marker” between the two measurements.

*et al.*

^{23}derived an equation that describes how the transverse susceptibility

*χ*

_{t}of a Stoner–Wohlfarth particle is dependent on the reduced field

*h*, defined as the ratio of the applied magnetic field and the particle’s anisotropy field,

^{24}

*M*

_{s}defined as the particle’s saturation magnetization and

*K*being the anisotropy constant; taking the

*z*axis as the direction of the applied DC field,

*ϕ*

_{K}and

*θ*

_{K}are the spherical polar angles between the applied DC field and the easy axis, illustrated in Fig. 2, and

*θ*

_{M}is the angle between the applied field and the magnetization.

*x*–

*y*plane, the transverse susceptibility is found by integrating Eq. (4) over ϕ

_{k}(0 to 2π) with θ

_{k}equal to π/2 to give

*θ*

_{M}equals zero) and h ≫ 1, Eq. (5) simplifies to

It should be noted that this methodology is still based on the shift in *H* of a single common feature as used in the established textured measurement. With no anisotropy features expected in these randomly packed powders, the only difference anticipated between the curves in the two orientations will be a widening of the overall peak centered on zero applied field when the demagnetization factor is higher: in this case across the thickness. This provides a difference of *ΔH* between the two curves at any given *χ*_{t} that is also investigated here as an alternative method of determining the demagnetization factors.

### D. Comparison with loop closure points

Samples of iron oxide MNPs sometimes show a small hysteresis (very narrow but distinct near-closed loops) despite direct size measurements indicating they are expected to be within the superparamagnetic regime. This is often attributed to magnetic interaction effects, due to agglomeration, in the dry powdered state.^{29} Selection of such a sample was done in this study as a means of comparing the new *χ*_{t} method with that of another accepted technique, as described by Bissell *et al.*,^{19} that is independent of particle orientation effects. Here, the hysteresis curve of the sample is measured using a standard Vibrating Sample Magnetometer (VSM), once with the applied magnetic field through the thickness and once with it through the width of the sample. The difference in the applied magnetic field needed to bring the hysteresis loop to closure between the two orientations is equal to the difference in the demagnetization fields. Hence, taking the common magnetization level that the closure point occurs at, the difference in the demagnetization factors can be calculated using Eq. (3). In this manner, a good comparison with the new *χ*_{t} method gives confidence in its applicability, in general, to that of randomly oriented nanoparticles, including extension to closed loops with no measurable hysteresis.

## III. RESULTS AND DISCUSSION

### A. Testing the model

The initial test of the model was to simulate a cubic sample (shape factor *D*_{z} value of ⅓^{13}), with the volume packing fraction held constant at 0.2, a typical packing fraction of the physical samples used in the experiments. The model was executed several times, each time using different side lengths so that the impact of the total number of particles involved in the simulation could be assessed. Using Eq. (2), this cubic system should give a sample demagnetization factor *N*_{z} of ⅓. The results for the mean average from the complete assemblies with cubic sides of 41, 51, 81, 101, and 151 particles, respectively, all produced the expected value as did the median and mode. Figure 3 shows slices taken through the model’s output for a cuboid of side length 101 particles, illustrating the spatial distribution of the constituent particles’ demagnetization factors. With the applied field in the *z* direction as shown, the largest *N*_{z} values are on the top and bottom surfaces, while the lowest values are observed on the side surfaces as expected. For clarity, only two “slices” through the 3D structure are shown and illustrate the tendency to converge on *N*_{z} values of ⅓ to the center particle. This is consistent with the dominance of the cuboid shape at the surfaces (and their associated free poles) reducing to that of a single isolated particle.

As the model was executed on a standard desktop computer, it was important to maximize its efficiency to obtain an acceptable run time. The details of this are given in the Appendix.

#### 1. Impact of packing fraction

A further test of the model’s consistency with the average value of the analytical model of Eq. (2) can be obtained by holding the shape factor constant while varying the packing fraction. This was executed for packing fractions between 0.1 and 0.5 in steps of 0.1 with a typical response for the average values shown in Fig. 4: in this case, a cuboid of length 201 particles along both the *x* and *y* axes and of 41 particles in the z axis was used, giving a shape demagnetization factor, *D*_{z}, of 0.6942.^{13} The mean of the particles’ demagnetization factors in the applied direction compares well with the solid line of Eq. (2) with the close fit and high degree of linearity shown. From a linear regression on the model data points, the extrapolation to a packing fraction of zero resulted in an *N*_{z} value of (0.3334 ± 0.0001) compared to the theoretical value of ⅓. These relationships were observed when other shape demagnetization factors were tested, including a perfect cube of factor ⅓.

Deviations from Eq. (2) for the median and modal points both show an increase above the mean that remains linear as illustrated by the regression curve fits and with the modal points showing the largest increase. This indicates a skew toward high *N*_{z} values in the number distribution consistent with greater interactions between the particles’ free surface poles as a function of the particles being packed closer together.

#### 2. Impact of shape demagnetization factor

By fixing the packing fraction while changing the aspect ratio of the model’s cuboid containing vessel, the impact of its shape demagnetization factor, *D*_{z}, on the average particle’s demagnetization factor may be investigated and compared with the output of Eq. (2). Figure 5 shows this relationship for a packing fraction of 0.2. Again, the mean from the computational model agrees with Eq. (2), being linear and with a regression fit yielding a gradient of (0.1992 ± 0.0005) compared to the expected value of 0.20 from the packing fraction *f*, and crosses the ordinate at (0.2667 ± 0.0003) compared to the value of 0.267 expected from $1\u2212f/3$. Similar relationships occurred for other tested packing fractions. Again, this gives confidence in the validity of the computational model, especially in terms of the numerical and spatial distribution of each *N*_{ijk} throughout a 3D nanostructure representing that of a real system. Deviations from Eq. (2) for the median and model outputs are non-linear as a function of vessel shape factor *D*_{z}. The values converge on a *D*_{z} of ⅓ as is expected by definition and consistency with Eq. (2), with values above the mean at *D*_{z} > ⅓ and below the mean at *D*_{z} < ⅓.

Insight into the reason for the median and modal curves being above or below the mean around the pivot point of *D*_{z} at ⅓ is given by consideration of the number distribution within the assembly. While the mean average is always in close agreement with Eq. (2), as illustrated in both Figs. 4 and 5, examination of the distribution of *N* also raises a question about if this is the best single measure? This is depicted by the example of Fig. 6 for a sample with packing fraction 0.2, cuboid lengths of 41 particles in the thickness and 401 particles in the width and length. The two distributions show the particles’ demagnetization factors when the saturating applied field is through the thickness and then through the length/width. In addition, the mean, median, and mode averages for the demagnetizing factors through the thickness are 0.4273, 0.4413, and 0.4480, respectively. Conversely, for the length/width, they are 0.2863, 0.2794, and 0.2762. Adding up the values in the three orthogonal axes gives a sum of ∼1 as expected.^{11} Looking at the skewness of the distributions qualitatively, there are a greater number of particles with higher values in the distribution above the pivot point of ⅓ in Fig. 6 than there is below for the out-of-plane and in-plane results: this is consistent with the corresponding *N*_{z} values above and below the mean of Fig. 5. This was investigated quantitatively by comparison with the results from the three experimental methodologies and included testing the validity of the proposed new transverse susceptibility techniques.

#### 3. Mean, median, and mode output

*p*,

*q*,

*r*,

*s*,

*t*, and

*u*being coefficients unique to either the median or mode. A non-linear regression fit

^{30–32}on the model data, over a range of packing fractions and shape demagnetization factors, i.e., fitted over the 3D surface, gave the coefficients and respective uncertainties: these are listed in Table I. The mean, median, and mode demagnetization values can then be calculated using Eqs. (2) and (7) for any powdered sample in a cuboid containing vessel and compared to experiment. The model data were extrapolated to extend beyond a packing faction of 0.52. A packing fraction of 0.52 was the uppermost value for this model, since at this point, the particles would be “physically” touching. This would be a physical limit, not a mathematical one.

Coeff. . | Median . | Mode . |
---|---|---|

p | 0.333 51 ± 0.000 05 | 0.333 26 ± 0.000 03 |

q | 0.664 7 ± 0.000 8 | 2.436 ± 0.004 |

r | −2.215 ± 0.015 | −6.576 ± 0.008 |

s | 2.001 ± 0.010 | 5.246 ± 0.005 |

t | 0.549 4 ± 0.002 3 | −0.108 3 ± 0.001 2 |

u | −0.333 59 ± 0.000 17 | −0.333 23 ± 0.000 09 |

Coeff. . | Median . | Mode . |
---|---|---|

p | 0.333 51 ± 0.000 05 | 0.333 26 ± 0.000 03 |

q | 0.664 7 ± 0.000 8 | 2.436 ± 0.004 |

r | −2.215 ± 0.015 | −6.576 ± 0.008 |

s | 2.001 ± 0.010 | 5.246 ± 0.005 |

t | 0.549 4 ± 0.002 3 | −0.108 3 ± 0.001 2 |

u | −0.333 59 ± 0.000 17 | −0.333 23 ± 0.000 09 |

It should be noted in Table I that the values for *p* and *u* agree within error and tend to their expected values of +⅓ and −⅓, respectively, thereby providing a consistency check with the calculations of the other coefficients over the whole range of packing fraction and shape demagnetization factors.

As expected for skewed distributions, the largest difference, Δ*N*_{z}, was between the mode and mean values (as compared to the median and mean), and this is shown in the 3D surface plot of Fig. 7. The pivot axis about Δ*N*_{z} = 0 along the packing fraction axis at a container demagnetization factor of ⅓ is clearly visible in the surface projection onto the 2D contour map. Significant differences typically occur at high packing fractions with a container shape demagnetization factor of about 0.70 for the positive peak and 0.15 for the negative trough. The good fit and associated small errors given in Table I mean that *N*_{d} can now be interpolated at any point in the structure without the need to run the model at a higher resolution and so becomes a powerful tool for further investigations.

### B. Magnetic characterization

#### 1. Transverse susceptibility measurements

Figure 8 shows the signal voltage output from the susceptometer when measuring the magnetite sample with its applied DC field across both the width and then the thickness of the sample. The signal voltage is proportional to the susceptibility of the sample.^{25} The feature of interest in these measurements is the point at which the curve descends with increasing magnetic field and becomes near linear. From Eq. (6), it was expected that at high applied fields, as the sample approaches saturation, the output should show an inverse relationship and would be expected to tend to linear behavior as it approaches the horizontal asymptote. The onset of linearity is shown in the inset of Fig. 8 to give field strengths of (492 ± 5) and (552 ± 5) kAm^{−1} across the width and thickness, respectively. Since this is occurring near to saturation magnetization, *N*_{z} − *N*_{x} has a value of (0.21 ± 0.03), determined using Eq. (3).

The other feature of interest in the curves is the clear widening of the peak width when the *H*_{DC} orientation is switched to be across the thickness instead of the width. This is expected due to *N*_{z} *> N*_{x} requiring a greater *H*_{DC} to be applied for the sample to be in the same *effective* field *H*_{eff} = *H*_{DC} − *N*_{z}*M*, where *N*_{z}*M* is the demagnetizing field of Eq. (1) at the common magnetization point *M*. As these two field points are at the same *H*_{eff}, the *χ*_{t} values are also common. At first sight, the variation of *ΔH* between the two curves as a function of *χ*_{t} does not seem to offer any further information. However, when compared to Eq. (3), it is apparent that a simple linear relationship with a gradient of *N*_{z} *− N*_{x} should be expected. This is shown in Fig. 9, where Δ*H* points from the *χ*_{t} curve are plotted as a function of magnetization values. In this case, normalized magnetization has been used and so *N*_{z} *− N*_{x} is found by multiplying the gradient by the saturation value. Figure 9 does show a linear relationship, but it is starting to break down at high magnetization levels. This is consistent with the difficultly in measuring the difference between the applied fields for the two orientations as the two curves start to converge at high *H*_{DC}, with only small changes in the *χ*_{t} signal as each curve approaches the horizontal. From the regression fit over the linear region, *N*_{z} *− N*_{x} was found to be (0.201 ± 0.009) with a significantly lower experimental error of 4% compared to the error of 14% found from the linear onset point.

#### 2. Magnetization curves

For comparison of the new transverse susceptibility results with those of an accepted method, magnetization curves measured in the same two orientations as the susceptometer are shown in Fig. 10. From TEM images (not shown), the particle size distribution was predominantly at 13 nm in a range (10 > d > 16) nm from a sample set of 100 measurements. At this small size, the behavior of these particles tends to superparamagnetic^{33} with near-closed loops. From the internal dimensions of the sample holder and the measured moment, the saturation magnetization of the magnetite powder was determined to be (286 ± 8) kAm^{−1}. From the measured mass of the sample and assuming all of it was magnetite with a density of 5.18 gcm^{−3} (5180 kgm^{−3}), the volume packing fraction was determined to be (0.21 ± 0.03).

The closure point occurs for the two orientations at (110 ± 10) and (160 ± 10) kAm^{−1} at an approximate magnetization of (0.80 ± 0.05) of the saturation magnetization, resulting in an *N*_{z} *− N*_{x} value of (0.22 ± 0.06). As can be seen in the inset of Fig. 10, it is difficult assessing the closure point: an issue raised in the original investigation^{19} and made more difficult here with the near-closed loops and hence the large error. However, comparison of this accepted method and its result with the value from our new and significantly lower error susceptibility technique is reasonable and, thus, gives confidence in using transverse susceptibility to characterize demagnetization factors in powders of randomly oriented nanoparticles. The inset around the coercivity points of Fig. 10 shows they are independent of the sample orientation, crossing the abscissa at the same values, as is expected due to *N*_{d} being zero at this juncture.

As fully closed loops do not allow the standard VSM method to be used, these results provide a method of extending characterization to the fully closed loops often found in measurements on superparamagnetic nanoparticles.

#### 3. Comparison of the model with experiment

VSM and susceptometer experimental techniques used to measure demagnetization factors examine the difference in the demagnetizing fields between two orientations; in our example, the sample was tested across its thickness and width. The difference in the applied field needed to drive the sample to a “feature” point in its characterization would be equal to the difference in the demagnetization fields between the two orientations. A direct comparison can then be drawn between the empirically determined value of N_{z} − N_{x} (the difference between the demagnetization factors for the two tested orientations) and that calculated using the model.

Comparisons of the model with experimental results are given in Table II. From the previously determined packing faction, *f*, of (0.21 ± 0.03) and sample vessel dimensions, it has a containing shape demagnetization factor^{13} along its thickness, *D*_{z}, of (0.86 ± 0.02) and along its width, *D*_{x}, of (0.11 ± 0.02). Equation (3) from the analytical model gives a value of N_{z} − N_{x} (0.16 ± 0.02). Evaluating Eq. (7) from the computational model for both median and mode gives values for N_{z} − N_{x} of (0.17 ± 0.02) and (0.19 ± 0.02), respectively. All the experimental results are larger than the theoretical values, the closest being the mode value, which was expected to be the largest theoretical value due to the skewed distribution evident in Fig. 6. Consequently, this value agrees with all three experimental values within error, although this must be viewed with caution when comparing to the VSM result with its large uncertainty. The two transverse susceptibility results are closer to the theoretical values, have smaller errors, and agree within their uncertainty.

Origin . | N_{z} − N_{x} value
. |
---|---|

Mean: Model [Eq. (2)] | 0.16 ± 0.02 |

Median: Model [Eq. (7)] | 0.17 ± 0.02 |

Mode: Model [Eq. (7)] | 0.19 ± 0.02 |

Hysteresis closure (VSM) | 0.22 ± 0.06 |

Linear onset (susceptometer) | 0.21 ± 0.03 |

Field difference (susceptometer) | 0.201 ± 0.009 |

*ζ*,

^{34,35}give an indication on whether two results

*x*

_{1}and

*x*

_{2}with respective uncertainties

*u*

_{1}and

*u*

_{2}are in agreement and are determined by

*ζ*≤ 0.55 being significantly lower than those of their median and mean counterparts in the range (0.79 ≥

*ζ*≥ 1.87). This is consistent with the single value of the experimental results being closer to the mode value of the model. Furthermore, it is consistent with the skewed number distributions, such as that of Fig. 6, as the

*N*

_{z}

*− N*

_{x}values for the mode will be greater than those of the median and mean, bringing them closer to the empirical measurement as would be expected for a majority of particles in the assembly being grouped around this value.

Experimental technique . | Zeta scores . | ||
---|---|---|---|

Mean . | Median . | Mode . | |

Hysteresis closure (VSM) | 0.95 | 0.79 | 0.47 |

Linear onset (suscept.) | 1.39 | 1.11 | 0.55 |

Field difference (suscept.) | 1.87 | 1.41 | 0.50 |

Experimental technique . | Zeta scores . | ||
---|---|---|---|

Mean . | Median . | Mode . | |

Hysteresis closure (VSM) | 0.95 | 0.79 | 0.47 |

Linear onset (suscept.) | 1.39 | 1.11 | 0.55 |

Field difference (suscept.) | 1.87 | 1.41 | 0.50 |

## IV. CONCLUSIONS

A computational demagnetization factor model has been developed based on particle surface poles in an assembly of nanoparticles and their interactions. By reducing the calculations to that of geometry-only, this model may be run on a standard personal computer (PC) and thus provide an effective and efficient technique to aid in the interpretation of experimental data.

Comparison with experiment using two new transverse susceptibility techniques showed the best agreement was with the Mode value of the model, resulting in zeta scores of ≤0.55. This is consistent with the skewed numerical distribution of factors in the model weighting a significant proportion of the experimental factors around the Mode and a better understanding of the single value measured.

With interest in superparamagnetic nanoparticles of continuing and growing interest, the susceptibility technique had to be extended to include such measurements. Particles with near-closed magnetization curves were used to test the validity of the two new susceptometer methods against an accepted measurement that utilizes loop closure points. All results agreed within error, giving confidence in the two new techniques. Of these, the multiple feature result was closest to the Mode value of the model and with the smallest experimental uncertainties.

Mapping the model demagnetization factors as a function of their position within the array allows the spatial distribution to be determined and provides insight into how the particles interact in either a global or local field. As this is essentially the same as considering magnetic islands or elements in other nanostructures, this will be increasingly important to understand in terms of the potential applications of 3D nanomagnetic materials that are currently being studied and developed. Preliminary modeling on changing the shape of particles within a given sample-containing vessel has illustrated the limitations of the spherical assumption inherent in Eq. (2) and along with the effects of agglomeration is the subject of future work.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: FURTHER MODEL DETAIL AND IMPLEMENTATION

The model used in this report is based on the work by Cookson^{20} in which he was studying the average demagnetization factor for tape recording media consisting of prolate spheroid particles, as shown in Fig. 11.

*o*,

*p*,

*q*). The particle has a semi-minor axis length of

*a*in both the

*x*(

*o*) and

*z*(

*q*) directions and a semi-major axis length of

*b*in the

*y*(

*p*) direction. The particle is fully magnetized in the

*+q*(

*z*) direction. By integrating the magnetic field contributions generated at the particle’s center by isolated poles on infinitesimally small sections of the surface, and dividing through by the particle’s magnetization, the demagnetization factor,

*N*

_{d}, of the particle can then be calculated as follows:

*q*given by

*p*by

For the contributions of the surface free poles of a neighboring particle, Eq. (A1) can be modified by simply adding the coordinates (*x*, *y*, *z*) of the center of the neighboring particle to (*o*, *p*, *q*). To simplify the coordinate system, Cookson used a scaling factor, *m* (with *m* = 1 representing touching particles). Neighboring particles in the simple lattice that he used were separated by a distance 2*ma* in the *x* and *z* directions and by a distance 2*mb* in the *y* direction. Therefore, he was able to calculate the contribution to the demagnetization field at a particle centered at the origin generated by any particle, centered at (*x*, *y*, *z*) for any given packing fraction.

For a sample consisting of a 3D lattice of evenly spread particles in a simple cubic packed arrangement, the model calculates the self-demagnetization field and the contributions of all the other particles for a particle positioned at the vertex (0, 0, 0). To illustrate this, Fig. 12 shows an example for a cubic sample consisting of spherical particles in a 3^{3} assembly. Spheres have been chosen to match that of the experimental study reported here and are simply implemented in Cookson’s model by setting *a* and *b* of equal length.

The demagnetizing factor of the particle at the (0, 0, 0) vertex can be found by summing all the contributions from the 27 particles (itself included). This surface integration process would then need to be repeated for all 27 particles in the system to be able to calculate the average demagnetization factor of the sample. As this does not take advantage of symmetries within the system, it was found to be computationally inefficient and was subject to further development in this study.

To reduce the number of times the surface integration must be performed, a new matrix is generated, created by making reflections of the matrix in Fig. 12 about the *xy*, *xz*, and *yz* planes: this new matrix has size 5^{3} (2n − 1)^{3}, where *n* is the length of the side of the cube. A representation is shown in Fig. 13. For simplicity, it only shows the *xz* plane at *y* equals 2. The particle at the origin of the original matrix, Fig. 12, is now at the center of the new, in our example, at position (2, 2, 2).

This new matrix is a demagnetization contribution map for our sample with the particle that we want to find its demagnetization factor at the center. If we think about stepping through each particle in turn in our sample, starting with the particle at (0, 0, 0), the demagnetization of each particle can be found by stepping a sub-matrix through the 5^{3} matrix, which is equal in size to the sample (3^{3} in this example), starting with highlighting elements (2, 2, 2) to (4, 4, 4). By summing the contents of this sub-matrix, this will yield the particle’s demagnetization field. So, for the particle at position (1, 0, 0) in the sample, the sub-matrix will be shifted to highlight elements (1, 2, 2) to (3, 4, 4), again summing up its contents to yield the demagnetization factor. This stepping process is repeated until all the particles’ demagnetization factors have been found.

In both stages, applying the surface integration and in the summation of demagnetization contributions, the actual numbers of integrations and summations are reduced due to several symmetries in the sample, further reducing the computational time.

Even with the above technique, there is still a limitation due to the computational speed: the model being executed on a standard desktop computer. The initial setup of the model required one numerical integration to be carried out per particle and additions equal to the square of the total number of particles. Ideally, since demagnetizing fields are dominated by the sample shape, the requirement for the model is that it contains the minimum number of particles without degradation of its output. This was tested by examining a cubic sample repeatedly with different edge lengths per repeat: Fig. 3 shows an example of such a setup. The mean, median, and mode average for all the edge lengths was ⅓. Figure 14 shows a comparison of the distribution of the constituent particle’s demagnetization factors for the largest (151 particle edge length) and the smallest (41 particle edge length) of the cubic simulations. The comparison shows that the distributions are similar in shape. The effect of decreasing the edge length in the simulation seems to be reduction in the “smoothness” of the distribution and with a noticeable shift of high valued demagnetization factor particles to clump and form a secondary high-end peak. A two-sample Kolmogorov–Smirnov test was performed with null hypothesis that the two datasets are drawn from the same distribution.^{36–38} The test indicated that the null hypothesis was not rejected at the 5% significance level, with a *p*-value of 0.8938.

To produce high *D*_{z} values requires a high aspect ratio between the lengths of the cuboid edges. The smallest length used in the model was 41 particles to keep distortion of the distributions of particle demagnetization factors to a minimum, but still allowing a model execution time that was acceptable.

## REFERENCES

_{(1−x)}La

_{x}Gd

_{y}Sm

_{z}Fe

_{(12−(z+y))}O

_{19}for recording media applications: Experimental and theoretical investigations,”

*In vitro*cellular study,”

_{5}Si

_{4},”

*Introduction to Magnetic Materials*

*Nonlinear Regression*

*Computing and Graphics in Statistics*

_{3}O

_{4}nanoparticles,”