Methods of investigating the demagnetization factors within assemblies of superparamagnetic nanoparticles

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,³ 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,⁶ and the importance of particle dipolar interactions in forming chains in such systems,^7–9 through to room-temperature magnetic refrigeration.¹⁰ 

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:¹² 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.¹³ 

Earlier modeling work¹⁶ 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¹⁷ and predicted close agreement with Eq. (2). A subsequent experimental study¹⁸ 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¹⁹ is also not applicable.

The results on randomly oriented magnetite (Fe₃O₄) 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.

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.¹³ 

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,²¹ 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.

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.²² In 1957, Aharoni et al.²³ calculated the transverse susceptibility as a function of applied field for a system of identical Stonar–Wohlfarth particles²⁴ 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²⁵ 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₃O₄) 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.²⁶ 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.

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.²⁸ This means that another feature needs to be identified that can then be used as a “marker” between the two measurements.

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.

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.²⁹ 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.,¹⁹ 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.

The initial test of the model was to simulate a cubic sample (shape factor D_z value of ⅓¹³), 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.

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.¹³ 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.

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−f/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.¹¹ 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.

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.

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.²⁵ 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⁻¹ 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_zM, where N_zM 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.

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³³ 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⁻¹. From the measured mass of the sample and assuming all of it was magnetite with a density of 5.18 gcm⁻³ (5180 kgm⁻³), 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⁻¹ 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¹⁹ 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.

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¹³ 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.

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.

The authors have no conflicts to disclose.

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

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

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 2ma in the x and z directions and by a distance 2mb 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³ 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³ (2n − 1)³, 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³ matrix, which is equal in size to the sample (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.