Nanoparticle Size Threshold for Magnetic Agglomeration and Associated Hyperthermia Performance

The likelihood of magnetic nanoparticles to agglomerate is usually estimated through the ratio between magnetic dipole-dipole and thermal energies, thus neglecting the fact that, depending on the magnitude of the magnetic anisotropy constant (K), the particle moment may fluctuate internally and thus undermine the agglomeration process. Based on the comparison between the involved timescales, we study in this work how the threshold size for magnetic agglomeration (daggl) varies depending on the K value. Our results suggest that small variations in K-due to, e.g., shape contribution, might shift daggl by a few nm. A comparison with the usual superparamagnetism estimation is provided, as well as with the energy competition approach. In addition, based on the key role of the anisotropy in the hyperthermia performance, we also analyse the associated heating capability, as non-agglomerated particles would be of high interest for the application.


Introduction
Based on the possibility to achieve local actuation by a harmless remote magnetic field, magnetic nanoparticles (MNPs) are very attractive candidates for novel medical applications [1,2]. Particularly iron oxides, based on their good biocompatibility [3], have been the subject of intense research in recent years, for example for magnetic hyperthermia cancer therapy [4,5] or drug release [6,7]. A key aspect defining the performance of the MNPs under external magnetic fields is their magnetic anisotropy, as it allows them to transform the absorbed electromagnetic energy into the required physical stimuli to promote specific cell behaviours, acting, in practice, as medical nanorobots [8].
Another key aspect to consider when dealing with magnetic nanoparticles for biomedical applications is the agglomeration likelihood, as it could affect not only the metabolising process but also the magnetic properties by changing the interparticle interactions [9]. Considering for example magnetic hyperthermia, it is known that the particles tend to agglomerate when internalized by the cells and such may lead to a decrease of the heating performance [10]. However, the opposite behaviour has also been reported, with an increase of the heat release if the particles form chains [11]. The problem is that while accounting for the effect of interparticle dipolar interactions is of primary importance for a successful application [12], the usual estimate of agglomeration likelihood, i.e., the ratio between the dipolar energy of parallel-aligned moments and thermal energy [13,14], in the limit case of touching particles (i.e., l cc = d), does not consider the magnetic anisotropy despite its key role in governing the magnetisation behaviour. In Equation (1), µ 0 = 1.256 × 10 −6 Tm/A is the permeability of free space, M S the saturation magnetisation, and l cc the center to center interparticle distance. In this work we suggest an approach to consider the magnetic anisotropy into the agglomeration likelihood based on the comparison of characteristic relaxation times. The complex role of the interparticle interactions often prompts researchers to the use of superparamagnetic (SPM) particles, with the idea that the rapid internal fluctuation of the particles' magnetic moments shall prevent their agglomeration. Thus, at the first approximation one could be tempted to consider that agglomeration will not occur for particles with blocking temperature (T B ) below the desired working temperature, since for T > T B the particles are in the SPM state (i.e., they behave as giant paramagnetic-like supermoments). However, it must be kept in mind that behaving SPM-like is not an absolute term, but it is defined by the experimental timescale. Thus, regarding agglomeration, a particle could be referred to as SPM if its Néel relaxation time, τ N , is smaller than the characteristic timescales that allow agglomeration, i.e., diffusion (τ di f f ) and rotation (τ B ) [15]. For the simplest case of uniaxial anisotropy, τ N can be estimated as [16] where the prefactor τ 0 usually ranges between 10 −9 and 10 −12 s, K is the uniaxial anisotropy constant, V the particle volume, and k B the Boltzmann constant. The diffusion time can be expressed as where x 2 is the mean square displacement for a translating Brownian particle [17], η the viscosity of the embedding media, and R hyd is the hydrodynamic radius, defined by the particle size plus a nonmagnetic coating of thickness t nm . For simplicity we consider spherical particles of diameter d. The rotation time τ B (also referred as Debye [18] or Brownian time [19]) is expressed as where V hyd is the hydrodynamic volume. Extensive details about the different relaxation mechanisms can be found in Coffey et al. [20]. The objective of this work is to estimate the size threshold for magnetic agglomeration, d aggl (i.e., size for which τ N > τ di f f , τ B , so that agglomeration is likely) in terms of K. Focusing on magnetite-like parameters based on its primary importance for bioapplications, we will consider different effective K values, which can be ascribed to dominance of shape anisotropy over the magnetocrystalline one [21,22]. Comparison will be made with the usual estimate of agglomeration likelihood, Equation (1). Then, the hyperthermia properties for the obtained d aggl will be studied. It must be recalled here the double role of K in the heating performance, as it determines both the maximum achievable heating [23,24] and the effectiveness in terms of field amplitude [25]; for completeness, this double role of K will also be briefly summarized. Please note that we are using "agglomeration" referring to a reversible process, distinct from the irreversible "aggregation" [26].

Size Threshold for Magnetic Agglomeration, d Aggl
To estimate d aggl we followed the same approach as we did in Ref. [15]: to compare the characteristic Néel, diffusion, and rotation times, to obtain d aggl as the size for which τ N > τ di f f , τ B . In Equations (3) and (4) we have at first set t nm = 0, and used η = 0.00235 kg/m·s, as in Ref. [19], which is comparable to that of HeLa cells for nm-scale dimensions [27]. We considered three cases for Equation (2): K = 8, 11, and 15 kJ/m 3 , i.e., values of the order found in the literature for magnetite particles [15,28,29]. The diffusion distance in Equation (3) is set as the interparticle distance at which the magnetostatic energy dominates over the thermal one, i.e., Γ > 1 [15], so that: Note that while we have chosen Γ = 1 to have a well defined criterion, agglomeration usually requires higher Γ values [30]. That is to say, we are searching for the lower d aggl boundary. With the same spirit, in Equation (1) we used M S = 4.8 × 10 5 A/m, i.e., the upper value for magnetite so that the interaction is, most likely, overestimated. The relaxation times as a function of the particle size are shown in Figure 1. In Figure 1 it is clearly observed how increasing K leads to more stable moments, thus favouring agglomeration at smaller sizes (from d aggl ∼ 25 nm for K = 8 kJ/m 3 , to d aggl ∼ 20 nm for K = 15 kJ/m 3 ). The inset shows the size dependence of Γ, which i) does not distinguish among particle characteristics (in terms of K, as previously mentioned), and; ii) predicts the dominance of the dipolar energy for much smaller particle sizes, with d aggl ∼ 7 nm. It is worth noting that the threshold value obtained for the K = 11 kJ/m 3 case, d aggl ≈ 22, is slightly bigger than the one previously reported in Ref. [15], for which d aggl ≈ 21 nm. This is due to the larger M S value used here, which enhances the diffusion time (through the diffusion distance, Equation (5)). Nevertheless, the great similarity despite the different M S values emphasizes the key role of the anisotropy in the agglomeration likelihood. The fact that so far we are not considering a nonmagnetic coating has a minor effect, as discussed next.
While we considered t nm = 0 in order to determine the boundary where clustering might appear, biomedical applications will always require a biocompatible nonmagnetic coating and therefore it is important to consider its role. That being said, the analysis shows that including a non-magnetic coating does not significantly modify the obtained threshold values: if considering t nm = 5 nm, d aggl increases just by ∼0.2 nm; and by ∼0.5 nm if t nm = 20 nm. This is illustrated in Figure 2A. A slightly larger influence is that of the viscosity of the embedding media, as illustrated in Figure 2B. Considering for example that of water, η = 0.001 kg/m·s, it is observed a 0.6 nm decrease from the average size. This value of viscosity is very significant because of being very similar to that of the cell's cytoplasm, although it must be kept in mind that large variations can be observed within the same cell type and among different types of cells [31]. A much higher viscosity would have a more significant effect, as illustrated for example with the macroscopic value of HeLa cells, η = 0.044 kg/m·s; nevertheless this values would be unrealistically high for the current particles, as such large η would correspond to much bigger sizes (over ∼86 nm for HeLa cells) because of the size-dependent viscosity at the microscale [27].
It is important to note that for the anisotropy values considered here, in all cases the size threshold d aggl is always defined by the competition between diffusion and Néel times, as τ B < τ di f f for all cases shown in Figure 2.
Next we will compare the predictions from the relaxation times with those obtained from zero field cooling/field cooling (ZFC/FC) measurements, the common way to estimate SPM behaviour (and thus likely non-agglomeration). Thus, if associating the onset of SPM behaviour to the blocking temperature, estimated as T B = KV/25k B [32], the corresponding threshold size, d T B , is readily obtained. The comparison between the agglomeration thresholds predicted by both approaches at room temperature (i.e., setting T B = 300 K) is summarized in Table 1.  Table 1 shows that, on average, the ZFC/FC approach predicts agglomeration to occur for sizes ∼4.2 nm bigger than the ones predicted by the relaxation times approach. In fact, the obtained d T B values correspond to a lower boundary, as they were estimated considering the limit case of no applied field, which is not possible in real ZFC/FC experiments. In general, applying the field during the measurements will result in lower T B [33][34][35], which would correspond to larger d T B (at least for the monodisperse case considered here; polydispersity might result in more complex scenarios [36][37][38]).

Associated Heating Performance
Similar to its importance on the agglomeration likelihood, the anisotropy plays a principal role in defining the hyperthermia performance. On the one hand, it defines the maximum energy that can be dissipated [4,39]: it is easy to see that for aligned easy axes the maximum hysteresis losses per loop are 8K [40] (2K for the random easy axes distribution [24]). On the other hand, it settles the response to the applied field (of amplitude H max ) through the anisotropy field, defined as H K = 2K/µ 0 M S [25,39]. This double key-role is illustrated in Figure 3, where the heating performance is reported in terms of the usual Specific Absorption Rate parameter, SAR, as SAR = A * f , where A stands for the area of the loop (hysteresis losses), and f is the frequency of the AC field. The simulations were performed in the same way as in Ref. [15]: we considered a random dispersion of monodisperse non-interacting nanoparticles (with the easy axes directions also randomly distributed), and simulated their response under a time varying magnetic field by using the standard Landau-Lifshitz-Gilbert equation of motion within the OOMMF software package [41]; for the random thermal noise (to account for finite temperature) we used the extension module thetaevolve [42].
The results displayed in Figure 3 show how, similar to how the apparently different hysteresis loops (A panel) are scaled by the anisotropy field (B panel), the apparently different SAR vs. H max trends scale if plotting SAR/(2K * f ) vs. H max /H K (the 2 f factor is just for normalisation). Note, however, that those results correspond to the Stoner-Wohlfarthlike case at T = 0 K [43]. In real systems with finite temperature, K also defines-as previously discussed-the stability of the magnetization within the particle. Thus, the ideal T = 0 K situation may vary significantly due to the effect of thermal fluctuations, as shown by the open symbols in Figure 3D, which correspond to the T = 300 K case for the two particle types considered. It is clearly observed how the strict H max ∼ 0.5H K threshold does not hold, and that the SAR is much smaller than the maximum possible.
The results shown in Figure 3 illustrate well the double role of the anisotropy on the heating performance. What is more, it must be kept in mind that the magnetic anisotropy is the only reason why small particles, such as the ones considered here of typical hyperthermia experiments (well described by the macrospin approximation) release heat under the AC field: if no anisotropy were to exist, there would be no heating (at least not for the frequencies and fields considered). This applies both to Néel and Brown heating, as with no anisotropy the magnetization would not transfer torque to the particle for its physical reorientation. Of course, larger sizes could display different heating mechanisms (due to non-coherent magnetization behaviour [44] or even eddy currents [45]), but that is not the present case.
We will analyse now the hyperthermia properties of the obtained threshold sizes for the different K values. Since the roles of surface coating and media viscosity are not very significant in relation to d aggl , we have focused, for simplicity, on the K − d aggl pairs summarized on Table 1, which would set an ideal limit. Thus, we simulated the dynamic hysteresis loops for the three cases considered, to then evaluate the heating capability. Some representative hysteresis loops are shown in Figure 4, for a small (205 kHz) and a large (765 kHz) frequency, as reported in experimental works [11,39,46]. Note that for simplicity we have considered a random easy-axes distribution, but depending on the specific experimental conditions (particle shape and properties of the embedding media, mainly), it might occur as an easy-axes reorientation leading to different hyperthermia performance [46][47][48][49].   (Table 1).
The results displayed in Figure 4 show large differences depending on the value of H max , illustrative of the minor-major loops competition [24,25]. This is further emphasized by the fact that a higher frequency results in narrower loops for the small fields, but wider for the larger ones. The differences between the different K cases are due to the different H max /H K ratios, as discussed in Figure 4. This is systematically analysed through the associated SAR values, shown in Figure 5. The results plotted in Figure 5 nicely fit within the general scenario discussed previously discussed (Figure 3): larger K allows higher SAR, provided enough field amplitude is reached (see corresponding 0.5H K values-vertical dashed lines-for reference); for small H max values, however, it may occur that smaller-K particles result in higher SAR due to the minor/major loops conditions, as discussed elsewhere [25]. This is an important aspect to consider regarding the variation in local heating due to size and/or anisotropy polydispersity [25,50]), as the difference between blocked and SPM particles would be the highest and thus also the locally released heat [25,51]. The results are also clearly divergent from the linear response theory model [19], for which SAR ∝ H max 2 ; this is not surprising as we are far from its applicability conditions (see e.g., Refs. [52,53] for a detailed discussion).
The predicted SAR values are quite large, implying that those particles would make efficient heat mediators. However, it is important to recall here that, so far, we made no considerations to the role of sample concentration. While this may appear reasonable as an initial approach, the fact is that the sample concentration is a key parameter to determine: first, because it defines the amount of deliverable heat; and second, because interparticle interactions (even without agglomeration) may significantly change the heating performance [11,24,29,39,54]. To provide some hint on how the sample concentration, c (% volume fraction), relates to the assumptions made, we can consider it through the nearest-neighbors interparticle distance, l NN . Following Tewari and Gokhale [55], for a random distribution of monodisperse particles we can approximate l NN as [24] l NN = (d + 2 · t nm ) · 0.4465 c 1/3 1 + 1.02625 c 0.64 2 3 .
Thus, by equating l NN to the diffusion distance x (Equation (5)) of the different d aggl values, we can obtain the related sample concentration threshold, c aggl . This is shown in Figure 6. The results shown in Figure 6 indicate that for bare particles (t nm = 0 nm) the applicability of the discussed arguments would be limited to very small concentrations, with c aggl ∼ 0.2% for the K = 1.1 × 10 4 J/m 3 case. However, the presence of a nonmagnetic coating significantly enlarges c aggl , as illustrated in the main panel for the cases of t nm = 5 and 10 nm. This trend is systematically summarized within the inset for the different values of K. It is observed that a coating of a few nanometers allows extending the applicability of our arguments within the 1-10% range. It is interesting to notice how with higher K this trend occurs with thinner t nm , as expected due to the smaller d aggl sizes. At this point it is worth noting that for iron oxides it has been reported the existence of an essentially non-interacting regime at low concentrations [56,57], characteristics that are very attractive for the application viewpoint as it would allow discarding the complex role of interparticle interactions. Nevertheless, it would be clearly interesting to consider the role of interparticle interactions in a more accurate way (e.g., by considering their role on the Néel relaxation time [58]), but unfortunately such an approach is difficult to carry out for a randomly distributed system.

Conclusions
We have presented an estimation of the threshold sizes for magnetic agglomeration of magnetite-like nanoparticles, depending on their magnetic anisotropy. Our approach was based on the consideration that K determines the stability of the particle magnetization and thus the likelihood of magnetic agglomeration, which involves physical translation and rotation of the particles themselves. By comparing the associated timescales, we have obtained that magnetite particles with usual anisotropy values should be relatively stable against agglomeration up to sizes in the range ∼20-25 nm in diameter. Then, we evaluated the associated hyperthermia performance, and found it to be relatively large (hundreds to thousands of W/g) for usual field/frequency conditions. The role of the nonmagnetic surface coating and that of the media viscosity appears secondary in determining the threshold sizes for agglomeration.
The initial estimates were made with no considerations about sample concentration, despite being a critical parameter for the application. In this regard, simple estimates indicate that the assumptions would be strictly valid only for very diluted conditions. However, the presence of a nonmagnetic coating might significantly extend the validity of the approximations to higher concentrations (up to about 10% volume fraction), showing that in this sense the nonmagnetic coating would play a key role.
It is important to recall that we have focused here on purely magnetic agglomeration, i.e., an ideal assumption which does not consider the complex situation often found experimentally, where other forces-of electrostatic nature-often play a central role in the agglomeration [59][60][61] and lead to agglomeration at smaller sizes [12]. Including those falls, however, is out of the scope of the present work, as it would result in a scenario that is too complicated. We have neither considered other important system characteristics such as polydispersity in size (both regarding aggregation [15] and heating [50]), and in anisotropy. The latter is expected to play a key role based on its primary importance both for agglomeration and heating, as discussed here. However, to the best of our knowledge its role has only been investigated regarding heating performance [25], but not regarding agglomeration likelihood. Considering the combined influence of those parameters clearly constitutes a challenging task for future works.
Finally, it is necessary to recall the conceptual character of the present work: while we have considered magnetite-like values for K and M S as a representative example, for simplicity those were taken as independent of size and temperature. However, it is well known that those may vary significantly within the size range of interest [62], and therefore the accurate determination of the agglomeration likelihood and hyperthermia performance would require including also those dependencies, together with the role of the nonmagnetic coating [63].

Data Availability Statement:
The reported data is available upon reasonable request to the authors.