Tuning field amplitude to minimise heat-loss variability in magnetic hyperthermia

Abstract

In this work, we theoretically investigate how shape-induced anisotropy dispersion and magnetic field amplitude jointly control both the magnitude and heterogeneity of heating in magnetite nanoparticle assemblies under AC magnetic fields. Using real-time Landau–Lifshitz–Gilbert simulations with thermal fluctuations, and a macrospin model that includes both the intrinsic cubic magnetocrystalline anisotropy and a shape-induced uniaxial contribution, we analyze shape-polydisperse systems under clinically and technologically relevant field conditions. We show that for relatively large particles, around 25 to 30 nm, the relative dispersion of local (single-particle) losses exhibits a well-defined minimum at moderate field amplitudes-between 4 to 12 mT-, hence identifying an optimal operating regime that minimizes heating heterogeneity while maintaining substantial power dissipation. The position of this optimal field depends mainly on particle size and excitation frequency, and only weakly on shape dispersion, offering practical guidelines for improving heating uniformity in realistic MFH systems.

---

1 Introduction

Magnetic fluid hyperthermia (MFH) refers to the generation of heat by magnetic nanoparticles (MNPs) subjected to an alternating (AC) magnetic field.¹ Although MFH has been primarily developed for cancer therapy,^2,3 its scope has progressively expanded to a wide range of applications, mainly biomedical,⁴ but also extending to non-biomedical fields such as catalysis⁵ and environmental remediation.⁶ For a recent review of non-oncological applications of MFH, see ref. 7.

Heat dissipation in MFH occurs at the single-particle level, predominantly through irreversible magnetization switching processes.⁸ As a consequence, the macroscopic heating response of an ensemble reflects an underlying distribution of local (single-particle) energy losses. In practical applications, large particle-to-particle variations in heat generation may lead to local over- or under-heating, and thus reduced efficiency or controllability of the process.⁹ From this perspective, MFH-based applications would benefit not only from maximizing the average heating power, but also from minimizing the dispersion of energy dissipation among individual nanoparticles. This consideration motivates the present work, which aims to explore how heating heterogeneity can be reduced in realistic MNP assemblies that are intrinsically polydisperse,¹⁰ both in size and—more critically—in magnetic anisotropy.

The influence of anisotropy dispersion on heating performance is particularly important when compared with size dispersion, because anisotropy plays a dual role in magnetic hyperthermia. Within the commonly accepted picture—where heat dissipation is mainly attributed to irreversible intraparticle magnetization processes^11–13—the particle volume V, together with the magnetic anisotropy, determines the height of the energy barrier E_B governing magnetization reversal. Size polydispersity has been shown to affect heating both locally and at the ensemble-averaged level, leading either to a reduction in energy dissipation (by about 30–50% for a standard deviation s = 0.20–0.25)¹⁴ or, in some cases, to an enhancement of heating performance.¹⁵

Magnetic anisotropy, on the other hand, also defines a characteristic field threshold for efficient heat generation. For randomly oriented assemblies of MNPs, the heating performance as a function of the field amplitude typically exhibits a sigmoidal behavior, evolving from negligible to substantial values above a threshold field proportional to the anisotropy field H_K.^5,16 As a result, the applied field amplitude emerges as a key experimental control parameter, capable of selectively activating subsets of nanoparticles depending on their anisotropy. Understanding how this field threshold is distributed in polydisperse systems is therefore central to controlling both the magnitude and the homogeneity of heat generation in MFH applications.^16,17

In this work, we investigate the interplay between particle-shape dispersity and field amplitude on local (single-particle) heating performance, focusing on magnetite nanoparticles, the most widely used material in magnetic hyperthermia. There is growing recognition that, for theoretical models to accurately reproduce experimental observations, deviations from ideal spherical particle shapes must be explicitly taken into account.^18–20 Such deviations, often described as small elongations from sphericity, give rise to a shape-induced uniaxial anisotropy that combines with the intrinsic magnetocrystalline anisotropy to determine the overall magnetic behavior.²¹ In magnetite, the intrinsic anisotropy is cubic, and neglecting this contribution may lead to an incomplete or misleading description of the heating response, particularly under low-field conditions where the competition between anisotropy terms becomes critical.²²

Here, we revisit and extend a theoretical prediction previously reported by our group, which showed that the variability of local (single-particle) heating can be reduced by appropriately tuning the field amplitude H_max as a function of the anisotropy dispersion.¹⁶ That earlier study was based on Monte Carlo simulations of systems with only uniaxial anisotropy, and did not account for real-time magnetization dynamics. Building on those results, we now investigate how shape-induced anisotropy polydispersity and field amplitude jointly control the distribution of single-particle hysteresis losses under dynamic conditions. Our objective is to identify operating regimes in which heating heterogeneity in shape-polydisperse nanoparticle assemblies is minimized, while maintaining significant overall heat generation. To this end, we employ real-time magnetization dynamics simulations and a comprehensive anisotropy model that combines uniaxial shape anisotropy with the intrinsic cubic magnetocrystalline contribution.²²

2 Physical model

We consider assemblies of magnetite nanoparticles described within the macrospin approximation, whereby each particle is represented by a single magnetic moment µ_i = M_sVm_i with M_s the saturation magnetization, V the particle volume, and m_i a unit vector defining the magnetization direction. The total magnetic energy of each nanoparticle includes three contributions: (i) the intrinsic cubic magnetocrystalline anisotropy; (ii) a uniaxial shape-induced anisotropy, accounting for small deviations from spherical geometry; and (iii) the Zeeman interaction with the applied time-dependent magnetic field. For magnetite, the cubic anisotropy constant is taken as K_c = −1.1 × 10⁴ J m⁻³.

Shape-induced anisotropy is modeled by approximating nanoparticles as prolate ellipsoids with aspect ratio r = c/a (c the long axis, a = b the short axes).²² The corresponding uniaxial anisotropy constant is given by

where N_c and N_a are the demagnetizing factors along the long and short axes, respectively. Their analytical expressions follow standard results for prolate ellipsoids.²³

A key novelty of the present work is that we consider a distribution of shape-induced uniaxial anisotropy constants to be directly correlated with the dispersion in aspect ratios, which we assume to follow a Gaussian distribution. In the present study, we will consider, for all particle sizes, that the system always has a shape dispersity with mean 〈r〉 = 1.1 and different standard deviations σ_r, that will be set to 0.1, 0.2. An illustrative scheme of shape distributions and the direct correlation with the associated K_u values is shown in Fig. 1.

Fig. 1 (a) Examples histograms of Gaussian distributions of aspect ratio r, centered around the mean 〈r〉 = 1.1, for different standard deviations σ_r, 0.0 (monodisperse) and 0.1 and 0.2. Schematic representations of the corresponding shape-polydisperse ensembles are shown in the upper right corner for illustration purposes. (b) Dependence of the corresponding uniaxial shape anisotropy constant K_u on r. The dashed vertical line highlights the mean axial ratio 〈r〉 = 1.1 of the distributions.

The heating performance is quantified via the Specific Loss Power (SLP), obtained from the area of the hysteresis loop as

where f is the field frequency and ρ the material density. The applied magnetic field is sinusoidal, with amplitude H_max.

Interparticle interactions and Brownian rotation are not considered in the present study, so that we can isolate the effects of intrinsic anisotropy polydispersity and magnetization dynamics on local heating behavior.

3 Computational details

Simulations are performed in the same way as in ref. 5 and 22 in which the magnetic evolution of a system of magnetite nanoparticles under a time-varying magnetic field is simulated with the OOMMF micromagnetic package, which numerically integrates the Landau–Lifshitz–Gilbert (LLG) equation,²⁴ with a random field to account for thermal effects.²⁵ Each nanoparticle is modeled as a single discretized cubic cell, whose volume is taken equal to the particle volume under the macrospin approximation. All the simulations are carried out for a system of N = 1000 particles at T = 300 K, using a Gilbert damping parameter α = 0.1. To obtain the heating performance, different hysteresis loops are simulated until convergence is reached, and then the hysteresis area is evaluated to obtain the SLP through eqn (2).

The present results correspond always to size-monodisperse cases, to focus on the role of shape (anisotropy) dispersity. Yet, to assess size effects, simulations have been performed for selected particle diameters D = 15, 20, 25, and 30 nm. The aspect ratio r is varied between 1.0 and 1.29, spanning the transition from spherical to moderately elongated particles. It is important to keep in mind that in each simulation the particle volume is the same for every particles, independently of its elongation; the diameter D is reported solely as a reference, for the ideally spherical particle equivalent. The orientations of both uniaxial and cubic anisotropy axes are randomly assigned. For each set of parameters, 30 independent realizations are simulated to ensure statistical robustness. A time-dependent magnetic field H(t) = H_max sin(2πft) is applied along a fixed spatial direction, considering two frequencies, f = 100 kHz and f = 1000 kHz. The material parameters for magnetite are M_s = 4.8 × 10⁵ A m⁻¹ and ρ = 5170 kg m⁻³.

4 Results and discussion

Munoz-Menendez and colleagues¹⁶ have demonstrated that, for blocked-like particles, the variation in local (single-particle) hysteresis losses can be controlled by adjusting the maximum magnetic field intensity H_max in accordance with the anisotropy polydispersity distribution parameter, σ_r, using Monte Carlo simulations; see Fig. 2, reproduced from ref. 16.

Fig. 2 Left axis: standard deviation of the normalized local hysteresis losses σ_HL, as a function of the normalised applied field amplitude, H_max/〈H_K〉, for a distribution in anisotropy constants K_u of σ = 0.2. Right axis: percentage of this standard deviation in relation to the global normalized hysteresis losses. The H_opt value indicates the optimal magnetic field for minimising dispersion in heat losses. Adapted from ref. 16 with permission from the Royal Society of Chemistry.

The standard deviation of the local losses, σ_HL, vs. H_max (left axis in Fig. 2), might apparently suggest that the smallest fields are the more adequate to diminish the dispersion in local heating. However, since at very small fields the heating is also very small, it is more appropriate to analyse the relative importance σ_HL against the correspondent global losses, as shown by the right axis in Fig. 2. In this case, it is clear that at low fields, the distribution of local hysteresis losses is quite significant compared to the relatively low global losses. In contrast, at higher fields, the dispersion of local losses diminishes significantly, showing a minimum at ∼H_max = 0.8H_K. The key objective of the present work is to investigate the possible existence of this optimising field amplitude, under the more realistic conditions of the present model.

4.1 Global (entire system) heating: role of aspect ratio

We first consider shape-monodisperse assemblies (σ_r = 0.0; part 4.1.1), in which all nanoparticles have the same axial ratio r, to establish a direct relationship between shape anisotropy and heating performance. We then extend the analysis to shape-polydisperse systems (σ_r > 1; part 4.1.2), where a distribution of aspect ratios is introduced to account for realistic nanoparticle ensembles. In this second stage we shall always restrict to the case 〈r〉 = 1.1 as a reference sample of what is usually considered (quasi)homogeneous experimentally, as ideal perfectly spherical (or cubic) particles cannot be synthesised yet.²⁶

4.1.1 Shape-monodisperse system. The influence of the aspect ratio on the magnetic response is first assessed through the hysteresis loops obtained under alternating magnetic fields. Representative loops for monodisperse systems with different aspect ratios r are shown in Fig. 3 for a 30 mT sinusoidal field at f = 100 and 1000 kHz.

Fig. 3 Hysteresis loops for D = 25 nm particles with different aspect ratios r. Solid lines correspond to f = 100 kHz and dashed lines to f = 1000 kHz, in both cases for µ₀H_max = 30 mT.

Two main features are observed in Fig. 3. First, increasing r changes the general shape of the curves, from very narrow and quasi-saturated at r = 1.0, widening and deviating from saturation with growing r. Second, increasing f does not significantly change the shape of each r case, only slightly opening the loop. These changes in the magnetic hysteresis directly translate into variations of the heating efficiency. The dependence of the SLP on the applied field amplitude H_max for different aspect ratios and frequencies is presented in Fig. 4. For clarity, the results are reported as SLP/f, which allows a direct comparison between the f = 100 and 1000 kHz cases, as estimates of the heating performance per cycle.

Fig. 4 SLP/f vs. µ₀H_max for D = 25 nm particles, for systems with different aspect ratios r; for both f = 100 kHz (solid lines) and f = 1000 kHz (dashed lines).

It is can be seen that variations in r lead to markedly different heating performances, with clearly distinct behaviour in the low and high field regimes. At high fields, a general trend of increasing SLP/f values with increasing r is observed. In contrast, the low-field behaviouris more intricate: significant heating occurs only for r = 1.1, while it remains negligible for both r = 1.0 and r = 1.2. Overall, increasing the frequency results in non-negligible losses for the ideal case r = 1.0 and accentuates the sigmoidal character of the curves for r > 1. A detailed explanation of these features can be found in ref. 5 and 22 in terms of the relative contributions of the intrinsic cubic magnetocrystalline anisotropy and the shape-uniaxial one to the energy barriers.

Briefly, at r = 1.1 the uniaxial shape-anisotropy contribution becomes comparable to the cubic one, increasing the effective barrier and hindering the ability of the magnetization to follow the alternating field. Consequently, hysteresis losses increase and measurable SLP/f values appear already at low field amplitudes (about 4 mT at 100 kHz). The crossover to uniaxial shape-anisotropy dominance occurs at r ≃ 1.22;²² therefore, at r = 1.22 the effective barrier is substantially larger (12 times higher than the cubic barrier), and appreciable losses are observed only above higher amplitudes (around 10 mT at 100 kHz). Beyond ∼15 mT, SLP/f tends to increase more moderately. Increasing the frequency further enhances dynamic losses, yielding the largest SLP/f for r = 1.2 at 1000 kHz.

4.1.2 Shape-polydisperse system. We now consider polydispersity in the uniaxial anisotropy values to account for shape-induced variability, as presented in Fig. 1. Since we assume non-interacting particles, to construct polydisperse ensembles we can simply weight the relative contribution of each fraction of particles with a given r according to the desired σ_r distribution. Thus, in the simulations we followed the same procedure described in Sec. 4.1.1, for different values of r. Specifically, we carried out 30 independent simulations for each σ_r value, varying r from 1.0 to 1.29 in steps of 0.01. For each r, hysteresis loops and the corresponding SLP values were computed by sweeping µ₀H_max from 1 to 10 mT in 1 mT increments, and from 10 to 40 mT in 2 mT increments. This procedure was repeated for the different nanoparticle diameters D. It is important to recall that for all cases 〈r〉 = 1.1, as stated earlier.

For a given degree of shape polydispersity, the effective SLP of the ensemble was obtained by weighting the individual SLP contributions corresponding to each axial ratio according to the Gaussian distribution defined by the chosen value of σ_r. In this way, the resulting SLP captures the combined response of nanoparticles with different shapes within the ensemble. To quantify the dispersion, we define the standard deviation of the local SLP values as

where P denotes the number of particle categories, N is the total number of nanoparticles in the ensemble, and SLP_j corresponds to the SLP associated with the j-th category. As an illustrative example, Fig. 5 shows the SLP obtained for the case σ_r = 0.2 as a function of µ₀H_max. The shaded region represents the standard deviation of the SLP values.

Fig. 5 SLP vs. µ₀H_max for shape-polydisperse MNPs with σ_r = 0.2, f = 1000 kHz, and D = 25 nm. The symbols represent the ensemble-averaged SLP, while the shaded region stands for the standard deviation of SLP, σ_SLP, arising from the distribution of particle shapes. The inset shows the SLP value as a function of µ₀H_max for other σ_r cases.

Two main features are observed in Fig. 5. First, the overall shape of the SLP vs. µ₀H_max curve closely resembles that of the shape-monodisperse systems with r > 1.0. Second, the width of the σ_SLP value increases significantly as the SLP tends to saturation. Both features originate from the progressive contribution of particles with larger shape anisotropy as H_max increases.

The first feature reflects the dominant role of the more elongated particles in determining the ensemble response. Accordingly, the polydisperse curves closely follow the monodisperse r = 1.1 and r = 1.2 cases, and the different σ_r curves nearly overlap (inset of Fig. 5); note that the σ_r = 0 curve in Fig. 5 corresponds to the monodisperse r = 1.1 case shown in Fig. 4.

The pronounced increase in σ_SLP near the saturation regime arises from the progressive activation of particles with higher uniaxial shape anisotropy as H_max increases. As discussed in ref. 16, particles with larger K_u dissipate more energy but require larger fields to reach major-loop conditions. Increasing H_max therefore enables these particles to contribute, simultaneously enhancing the global SLP and broadening the local losses.

A simplified, more intuitive explanation of this threshold at the saturation part of the curve follows from the crossover between the r = 1.1 and r = 1.2 curves in Fig. 4, which occurs at µ₀H_max ∼14 mT: before this field, the highest contribution to the SLP is from the r = 1.1 case, and the deviation in losses is more moderate (both more and less elongated particles give very small SLP), whereas above it the r = 1.2 contributes more strongly, raising the average SLP while simultaneously broadening the distribution in local SLP values.

The results displayed in Fig. 5 (similar trends are obtained for other σ_r and D, not shown) indicate that global SLP and its dispersion are not directly correlated. While the average SLP increases monotonically with H_max, σ_SLP also grows, so that the latter alone does not provide a reliable measure of the relevance of local heating effects. Identifying a more appropriate descriptor is the subject of the following subsection.

4.2 Dispersion in local (single-particle) heat production

To analyse local heating effects in a physically meaningful way, we focus on the dispersion normalized by the ensemble-averaged SLP, σ_SLP/〈SLP〉. In contrast to the σ_SLP value alone, this dimensionless quantity directly quantifies the relative heterogeneity of single-particle heat production within the ensemble, and therefore provides an appropriate descriptor for assessing the relevance of local losses. This approach follows the same rationale adopted by Muñoz-Menéndez et al.,¹⁶ who analysed the dispersion of local hysteresis losses relative to the corresponding global losses. The difference is that in our present work we can consider real-time dynamics through the LLG equation, which were not accessible through the Metropolis Monte Carlo method of ref. 16.

4.2.1 Comparison with previous work. We first compare the results of our dynamic model against the main result reported in Fig. 2. Fig. 6 shows the analogue of the original plot for the same dispersion parameter, σ_r = 0.2, now formulated in terms of the SLP instead of the hysteresis losses HL. Note that in ref. 16 no particle size was specified, implicitly assuming that the model applies to blocked-like particles, as reported elsewhere for the same technique.²⁷ Blocked-like behaviour can be safely expected for magnetite particles with D = 25 nm at f = 1000 kHz, as the case shown in Fig. 5. To enable a more direct comparison with Fig. 2, H_max has also been normalised by the average anisotropy field of the particles, 〈H_K〉; for σ_r = 0.2 and 〈r〉 = 1.1, µ₀〈H_K〉 = 22.38 mT.

Fig. 6 Absolute values of heat-losses dispersion, σ_SLP (black circles, left axis) and normalised correspondent, σ_SLP/〈SLP〉 (red squares, right axis), as a function of H_max/〈H_K〉, for D = 25 nm, f = 1000 kHz, and σ_r = 0.2. The H_opt value indicates the optimal field for minimising SLP dispersion.

Fig. 6 shows an overall increase of σ_SLP vs. H_max, as expected (note that this data corresponds to the shaded area in Fig. 5). Remarkably, a secondary minimum is observed at H_max ∼0.4〈H_K〉, which was not evident from Fig. 5. Importantly, such feature is also present in the σ_HL data of Fig. 2, which also follows an overall increase with the presence of a secondary minimum at intermediate fields. Notable differences are that in the original case an initial constant (nearly-zero) range was followed by an absolute maximum, whereas in the present work no constant range is observed, and the maximum is secondary. In addition, the (normalised) field values at which the minimum appears are significantly smaller for the σ_SLP case.

In contrast, the normalized quantity σ_SLP/〈SLP〉 exhibits a markedly different behaviour. It is largest at low fields (indicating strong heterogeneity in local heating) and decreases rapidly until reaching a minimum at a optimal field H_opt ∼0.4〈H_K〉, where heat production becomes more homogeneous. A subsequent weak increase is observed for higher H_max, followed by a plateau for H_max ≳ 0.6〈H_K〉. This confirms that, while σ_SLP alone is not a suitable descriptor (as anticipated from Fig. 5), the normalized quantity σ_SLP/〈SLP〉 captures the existence of an optimal field regime, in direct analogy with ref. 16. A key difference, however, is that in the original case a similar trend emerged only after an initial constant regime, which is entirely absent in the present data.

The main differences between our present results and those reported in ref. 16 can be attributed to both the different physical model and the computational approach. In particular, in the low-field regime the key difference arises from the inclusion of cubic anisotropy in the present work, which enables dissipation at much smaller fields.²² Thus, with increasing H_max, there is a progressive contribution of the total losses from particles with higher uniaxial-shape anisotropy, which can undergo magnetisation reversal via the effective reduction of the uniaxial energy barrier by the cubic contribution. In contrast, the uniaxial-only anisotropy particles considered in ref. 16 required substantially larger fields to dissipate energy, resulting in negligible and nearly constant losses at small H_max values. In addition, part of the difference arises from the computational methodology itself. While ref. 16 employed a Metropolis Monte Carlo scheme based on energy minimization and thus probing quasi-static configurations, our approach relies on the numerical integration of the Landau–Lifshitz–Gilbert equation and explicitly accounts for real-time magnetization dynamics. This allows us to capture frequency-dependent and transient effects in the losses, which become particularly relevant in the low-field regime and further contribute to the absence of an initial constant-loss region in the present results.

4.2.2 Determination of optimum field conditions. A key outcome of the previous analysis is that the more comprehensive model employed in the present work reproduces the occurrence of the optimal field H_opt reported in ref. 16. This provides an independent validation of those findings using a more robust computational procedure and a more sophisticated physical description. In the following, we extend that analysis to consider the role of particle size (D), frequency (f), and shape polydispersity (σ_r).

We begin with a representative case, D = 25 nm and f = 1000 kHz, and analyse the behaviour of σ_SLP/〈SLP〉 as a function of the field amplitude µ₀H_max for several values of σ_r (Fig. 7).

Fig. 7 σ_SLP/〈SLP〉 vs. µ₀H_max for f = 1000 kHz and D = 25 nm. The inset highlights the corresponding optimal fields H_opt.

Fig. 7 shows a qualitatively similar behaviour for all shape-polydisperse samples. In all cases, σ_SLP/〈SLP〉 decreases from a pronounced maximum at low field amplitudes to an absolute minimum around 7–9 mT, followed by a slight increase and a broad plateau at higher fields. The main effect of increasing σ_r is a progressive upward shift of the curves, i.e., larger values of σ_SLP/〈SLP〉 across the full field range, together with a displacement of the minimum towards larger µ₀H_max. In particular, this implies a shift of H_opt towards larger values as shape dispersion increases.

The value of H_opt can be intuitively interpreted as the field threshold beyond which adding additional particle categories characterized by their respective K_u values no longer results in a meaningful increase of the average SLP. Although including more categories initially enhances the average SLP by incorporating particles that undergo irreversible magnetization processes, it simultaneously increases the dispersion of local (single-particle) heating. When the contribution from particles with larger K_u becomes dominant, the average SLP enters a saturation regime. In this situation, further increasing H_max leads only to a modest rise in SLP, despite the persistently large σ_SLP. It is important to stress that this optimum field should not be interpreted as the threshold separating minor from major hysteresis loops, nor as the crossover from linear-response to nonlinear behavior. Such an interpretation is incompatible with the fact that no H_opt is observed in the monodisperse system, even though the minor-to-major loop crossover still occurs in that case.

Although not shown here, a similar trend in the σ_SLP/〈SLP〉 vs. σ_r is observed for the same particle size at f = 100 kHz, and the overall behaviour remains comparable for D = 30 nm at both frequencies. However, for smaller particles (D = 15 and 20 nm), the σ_SLP/〈SLP〉 vs. H_max curves change significantly.

This is illustrated in Fig. 8 for D = 20 nm, where an apparently very different behaviour with respect to Fig. 7 is observed, as no minimum can be identified. For f = 1000 kHz, the curves display a well-defined maximum within the same field range where the D = 25 nm case exhibited a minimum. For f = 100 kHz, the dispersion shows only a weak decrease with increasing H_max; it is worth noting that such behaviour does not change if extending the field range up to 40 mT (not shown). In both cases, σ_SLP/〈SLP〉 remains substantially larger (by roughly a factor of 2–3) than the minimum values obtained for the larger nanoparticles. The results for the D = 15 nm case (not shown) are very similar to the f = 100 case of Fig. 8 at both frequencies, exhibiting a smooth decrease with increasing H_max for both frequencies. Likewise, for ideal monodisperse (σ_r = 0) and even nearly monodisperse systems (σ_r = 0.01), the dispersion does not exhibit a clear minimum, preventing the definition of a optimal magnetic field.

Fig. 8 Variation of σ_SLP/〈SLP〉 as a function of σ_r, for D = 20 nm, respectively. Empty symbols (dashed lines) correspond to f = 1000 kHz (a), and filled symbols (solid lines) to f = 100 kHz (b).

One might hypothesize that what is observed for the D = 20 nm case in Fig. 8 is not fundamentally different from that of the 25 nm case of Fig. 7, and that the curves would eventually show the same minima, just at much larger fields. However, the comparison between the D = 25 and D = 30 nm cases does not support that hypothesis, as will be shown later in Fig. 9. Furthermore, even if this were the case, it is clear that size polydispersity plays a key role in determining the dispersion of local heating, as discussed elsewhere.²⁸

Fig. 9 Optimal magnetic field H_opt as a function of σ_r, for D = 30 nm (blue squares) and D = 25 nm (green circles), respectively. Empty symbols (dashed lines) correspond to f = 1000 kHz, and filled symbols (solid lines) to f = 100 kHz.

Nevertheless, since studying the role of polydispersity lies outside the scope of the present work, we focus in the following on the cases where H_opt is well defined, namely D = 25 and 30 nm for sufficiently broad shape dispersions. The corresponding values of H_opt as a function of σ_r for f = 100 and 1000 kHz are summarized in Fig. 9. At f = 1000 kHz, H_opt lies in the range 11–13 mT for D = 30 nm, whereas for D = 25 nm it is shifted to lower values, around 8–9 mT, indicating a clear dependence on particle volume. Reducing the particle diameter progressively shifts the minimum towards lower field amplitudes and ultimately suppresses its appearance. Moreover, decreasing the excitation frequency further lowers the optimal field: at f = 100 kHz, H_opt is found around 4–6 mT for D = 25 nm and 8–10 mT for D = 30 nm. Overall, these results show that the emergence of H_opt, associated with minimal hysteresis-loss dispersion, is primarily governed by nanoparticle diameter and excitation frequency, while the influence of σ_r remains secondary.

The results displayed in Fig. 9 define an optimization scenario for MFH, which can be practically accessed by tuning the field amplitude. Remarkably, the predicted optimal fields fall within the biologically acceptable window for in vivo biomedical applications, usually expressed in terms of the H_max × f constraint. Taking the reference value reported by the Atkinson and collaborators,²⁹ f × H_max ≤ 4.8 × 10⁸ A (m⁻¹ s⁻¹) (often referred to as Brezovich limit³⁰), for f = 100 kHz the corresponding maximum field amplitude is µ₀H_max ∼6.1 mT, i.e., very close to the H_opt values obtained for D = 30 and 25 nm particles. It is worth noting that Pankhurst and collaborators³¹ suggested that reducing the duty cycle could allow the safe threshold to increase to approximately twice the Brezovich limit. In that case, again for f = 100 kHz, the upper field limit would be µ₀H_max ∼12.5 mT, which also lies within the range of optimal fields shown in Fig. 9. Beyond biomedical applications, the results shown in Fig. 9 may also have significant implications for other non-biomedical fields. For example, they are directly relevant to sequential catalytic reactions requiring well-separated triggering temperatures.⁵ In such cases, the absence of physiological constraints on f × H_max would, in principle, allow independent tuning of the frequency to control the target SLP, while still operating within a regime of minimized local-heating dispersion.

It is important to emphasize that, although Fig. 9 identifies an optimal field amplitude for given D and f conditions that appears to be only weakly dependent on σ_r, a finite dispersion in heat production persists and must be taken into account in any practical application. To quantify its magnitude, Fig. 10 shows the corresponding values of σ_SLP/〈SLP〉 associated with the data presented in Fig. 9.In Fig. 10 it is clearly observed that, while the position of H_opt is only marginally affected by shape polydispersity, the relative dispersion of the SLP increases systematically with increasing σ_r for the different particle sizes and excitation frequencies. This behaviour demonstrates that shape-induced polydispersity substantially enhances local heating fluctuations, even at the field at which hysteresis losses are minimized. Consequently, more polydisperse systems exhibit more pronounced spatial variations in heat dissipation, whereas highly monodisperse ensembles show markedly reduced local losses at H_opt. These results underscore the importance of achieving high monodispersity to ensure spatially homogeneous and efficient heating in MFH applications. These trends further indicate that, while field tuning provides an effective means to reduce local heating variability, it cannot fully compensate for strong shape disorder. Therefore, achieving highly homogeneous heating in MFH requires a combined strategy involving both optimized excitation conditions and careful control of nanoparticle shape distributions during synthesis.

Fig. 10 Variation of σ_SLP/〈SLP〉 as a function of σ_r, for D = 30 nm (blue squares) and D = 25 nm (green circles), respectively. Empty symbols (dashed lines) correspond to f = 1000 kHz, and filled symbols (solid lines) to f = 100 kHz. The field amplitudes considered here correspond to the optimal field values, H_opt, shown in Fig. 9.

5 Conclusions

We have presented a detailed numerical study of the interplay between shape-induced anisotropy polydispersity and applied magnetic field amplitude on both global and local heating performance in magnetic fluid hyperthermia. By employing real-time magnetization dynamics simulations and explicitly accounting for the coexistence of intrinsic cubic magnetocrystalline anisotropy and shape-induced uniaxial anisotropy in magnetite nanoparticles, we extend previous theoretical descriptions¹⁶ toward a more realistic representation of experimental systems.

A central result of this study is the demonstration that ensemble-averaged heating efficiency, 〈SLP〉, and local heating homogeneity, quantified in terms of standard deviation of SLP, σ_SLP are not trivially correlated. While the global specific loss power increases monotonically with increasing field amplitude, the dispersion of single-particle heat losses shows a non-monotonic behavior when evaluated relative to the ensemble average. In particular, the analysis of the normalized standard deviation, σ_SLP/〈SLP〉, reveals a local optimal field, H_opt, at which local hysteresis losses are minimized relative to global losses. This confirms, using a more comprehensive physical model and real-time dynamics, the existence of an optimal operating field previously predicted.¹⁶

The range of H_opt is governed primarily by particle diameter and field frequency, while the influence of polydispersity (σ_r) remains minor. A well-defined H_opt is observed for larger particles (D = 25–30 nm), whereas smaller (D = 15–20 nm) or monodisperse particles do not exhibit such behavior. These results underscore the key roles of particle size, shape-anisotropy distribution, and field frequency in achieving uniform and efficient heating in magnetic hyperthermia.

Regarding the direct translation of these theoretical results to practical implementations, it is important to keep in mind some key factors. First, the results were specifically obtained for magnetite nanoparticles with an average aspect ratio dispersion of 〈r〉 = 1.1. Other compositions and/or elongations would give very different uniaxial contributions.³² Second, we considered a size-monodisperse system, while it is known that size polydispersity is intrinsic to a real system, and has a strong influence on local (single-particle) heating.⁹ Third, the simulations were run for non-interacting conditions, which would need to be taken into account if aggregates are formed.^21,33

Another aspect that deserves special attention is that, while the easy axes were assumed to be randomly oriented, the average heating performance was computed for a given r category without considering the actual particle elongation direction relative to the applied field. However, the orientation of the anisotropy axes with respect to the field can have a significant impact on heating performance.³⁴ Moreover, the present treatment implicitly corresponds to a textitfrozen ferrofluid scenario. In contrast, in a viscous environment where particles are free to rotate, a dynamical reorientation of the anisotropy axes may occur,^35,36 leading to additional modifications of the heating performance.^37,38

Taken together, these considerations indicate that identifying optimal H_opt conditions for a given nanoparticle sample may require a dedicated analysis tailored to the specific physical conditions. While the present results provide general guidelines for minimizing heating heterogeneity, extending them to more complex and realistic scenarios involving orientational effects and particle mobility will require further investigation.

Author contributions

Conceptualization: D. S. Data curation: N. Ç., I. L. V. Formal analysis: N. Ç., I. L. V., D. S., Ò. I. Funding acquisition: D. S. Investigation: N. Ç., I. L. V., D. S., Ò. I. Methodology: N. Ç., I. L. V., D. S., Ò. I. Project administration: D. S., Ò. I. Software: N. Ç., I. L. V. Supervision: D. S., Ò. I. Validation: N. Ç., I. L. V., D. S., Ò. I. Visualization: N. Ç., I. L. V. Writing – original draft: N. Ç. Writing – review & editing: N. Ç., I. L. V., D. S., Ò. I.

Conflicts of interest

There are no conflicts to declare.

Data availability

All the simulation results presented in this article have been obtained using the micromagnetic simulator OOMMF, available from the NIST website, https://math.nist.gov/oommf/software.html.

Acknowledgements

We acknowledge financial support by Spanish Ministerio de Ciencia, Innovación y Universidades through projects PID2024-157172NB-I00 and CNS2024-154574, “ERDF A way of making Europe”, by the “European Union”, and Catalan DURSI (2021SGR0032). Xunta de Galicia is acknowledged for projects ED431F 2022/005 and ED431B 2023/055. AEI is also acknowledged for the Ramón y Cajal grant RYC2020-029822-I that supports the work of D.S. We acknowledge the Centro de Supercomputacion de Galicia (CESGA) for computational resources.

References

1. E. A. Perigo, G. Hemery, O. Sandre, D. Ortega, E. Garaio, F. Plazaola and F. J. Teran, Appl. Phys. Rev., 2015, 2, 041302.
2. R. Ivkov, Int. J. Hyperthermia, 2013, 29, 703–705.
3. M. Latorre and C. Rinaldi, Puert. Rico Health Sci. J., 2009, 28, 123–129.
4. A. Rivera-Rodriguez and C. M. Rinaldi-Ramos, Annu. Rev. Chem. Biomol. Eng., 2021, 12, 163–185.
5. J. G. Ovejero, I. Armenia, D. Serantes, S. Veintemillas-Verdaguer, N. Zeballos, F. López-Gallego, C. Gruttner, J. M. de la Fuente, M. d. P. Morales and V. Grazu, Nano Lett., 2021, 21, 7213–7220.
6. B. Baral, A. Altaee, K. Simeonidis and A. K. Samal, Front. Chem., 2024, 12, 1362033.
7. H. Gavilán, A. Gallo-Córdova, M. Rábade-Chediak, A. Páez-Rodríguez, M. d. P. Morales and L. Gutiérrez, Nanoscale, 2025, 17, 27734–27761.
8. F. Soetaert, P. Korangath, D. Serantes, S. Fiering and R. Ivkov, Adv. Drug Deliv. Rev., 2020, 163, 65–83.
9. C. Munoz-Menendez, I. Conde-Leboran, D. Serantes, R. Chantrell, O. Chubykalo-Fesenko and D. Baldomir, Soft Matter, 2016, 12, 8815–8818.
10. R. E. Rosensweig, J. Magn. Magn. Mater., 2002, 252, 370–374.
11. D. Serantes, R. Chantrell, H. Gavilán, M. d. P. Morales, O. Chubykalo-Fesenko, D. Baldomir and A. Satoh, Phys. Chem. Chem. Phys., 2018, 20, 30445–30454.
12. R. Di Corato, A. Espinosa, L. Lartigue, M. Tharaud, S. Chat, T. Pellegrino, C. Ménager, F. Gazeau and C. Wilhelm, Biomaterials, 2014, 35, 6400–6411.
13. D. Soukup, S. Moise, E. Céspedes, J. Dobson and N. D. Telling, ACS Nano, 2015, 9, 231–240.
14. A. P. Khandhar, R. M. Ferguson and K. M. Krishnan, J. Appl. Phys., 2011, 109, 07B310.
15. M. Jeun, S. Bae, A. Tomitaka, Y. Takemura, K. H. Park, S. H. Paek and K. W. Chung, Appl. Phys. Lett., 2009, 95, 083108.
16. C. Munoz-Menendez, D. Serantes, J. M. Ruso and D. Baldomir, Phys. Chem. Chem. Phys., 2017, 19, 14527–14532.
17. G. Vallejo-Fernandez and K. O'Grady, Appl. Phys. Lett., 2013, 103, 142417.
18. A. G. Roca, L. Gutiérrez, H. Gavilán, M. E. F. Brollo, S. Veintemillas-Verdaguer and M. del Puerto Morales, Adv. Drug Deliv. Rev., 2019, 138, 68–104.
19. G. Salazar-Alvarez, J. Qin, V. Sepelak, I. Bergmann, M. Vasilakaki, K. Trohidou, J. Ardisson, W. Macedo, M. Mikhaylova, M. Muhammed, M. D. Baró and J. Nogués, J. Am. Chem. Soc., 2008, 130, 13234–13239.
20. J. M. Paez-Muñoz, F. Gámez, Y. Fernández-Afonso, R. Gallardo, M. P. Leal, L. Gutiérrez, J. M. de la Fuente, C. Caro and M. L. García-Martín, J. Mater. Chem. B, 2023, 11, 11110–11120.
21. N. A. Usov, M. S. Nesmeyanov, E. M. Gubanova and N. B. Epshtein, Beilstein J. Nanotechnol., 2019, 10, 305–314.
22. D. Faílde, V. Ocampo-Zalvide, D. Serantes and Ò. Iglesias, Nanoscale, 2024, 16, 14319–14329.
23. G. Bertotti, Hysteresis in Magnetism: for Physicists, Materials Scientists, and Engineers, Academic Press, 1998.
24. M. J. Donahue and D. G. Porter, OOMMF User's Guide: Version 1.0, 1999.
25. O. Lemcke, Thetaevolve Module for OOMMF, 2004.
26. I. López-Vázquez, D. Serantes and Ò. Iglesias, J. Magn. Magn. Mater., 2025, 173751.
27. D. Serantes, D. Baldomir, C. Martinez-Boubeta, K. Simeonidis, M. Angelakeris, E. Natividad, M. Castro, A. Mediano, D.-X. Chen and A. Sanchez, et al., J. Appl. Phys., 2010, 108, 073918.
28. C. Munoz-Menendez, I. Conde-Leboran, D. Baldomir, O. Chubykalo-Fesenko and D. Serantes, Phys. Chem. Chem. Phys., 2015, 17, 27812–27820.
29. W. J. Atkinson, I. A. Brezovich and D. P. Chakraborty, IEEE Trans. Biomed. Eng., 2007, 70–75.
30. I. A. Brezovich, Psychol. Med. Monogr. Suppl., 1988, 16, 82–111.
31. M. K. Kwok, C. C. Maley, A. Dworkin, S. Hattersley, P. Southern and Q. A. Pankhurst, Appl. Phys. Lett., 2023, 122, 240502.
32. N. A. Usov and J. M. Barandiarán, J. Appl. Phys., 2012, 112, 054901.
33. V. Russier, C. De-Montferrand, Y. Lalatonne and L. Motte, J. Appl. Phys., 2013, 114, 143902.
34. I. Conde-Leborán, D. Serantes and D. Baldomir, J. Magn. Magn. Mater., 2015, 380, 321–324.
35. N. A. Usov, J. Appl. Phys., 2010, 107, 123916.
36. S. Ota and Y. Takemura, Appl. Phys. Express, 2017, 10, 085001.
37. H. Mamiya and B. Jeyadevan, Sci. Rep., 2011, 1, 157.
38. K. Simeonidis, M. P. Morales, M. Marciello, M. Angelakeris, P. de La Presa, A. Lazaro-Carrillo, A. Tabero, A. Villanueva, O. Chubykalo-Fesenko and D. Serantes, Sci. Rep., 2016, 6, 38382.

---