Beyond Newton's law of cooling in evaluating magnetic hyperthermia performance: a device-independent procedure

Abstract

Accurate knowledge of the heating performance of magnetic nanoparticles (MNPs) under AC magnetic fields is critical for the development of hyperthermia-mediated applications. Usually reported in terms of the specific loss power (SLP) obtained from the temperature variation (ΔT) vs. time (t) curve, such an estimate is subjected to a huge uncertainty. Thus, very different SLP values are reported for the same particles when measured on different equipment/in different laboratories. This lack of control clearly hampers the further development of nanoparticle-mediated heat-triggered technologies. Here, we report a device-independent approach to calculate the SLP value of a suspension of magnetic nanoparticles: the SLP is obtained from the analysis of the peak at the AC magnetic field on/off switch of the ΔT(time) curve. The measurement procedure, which itself constitutes a change of paradigm within the field, is based on the heat diffusion equation, which is still valid when the assumptions of Newton's law of cooling are not applicable, as (i) it corresponds to the ideal scenario in which the temperature profiles of the system during heating and cooling are the same; and (ii) it diminishes the role of coexistence of various heat dissipation channels. Such an approach is supported by theoretical and computational calculations to increase the reliability and reproducibility of SLP determination. Furthermore, the new methodological approach is experimentally confirmed, by magnetic hyperthermia experiments performed using 3 different devices located in 3 different laboratories. Furthermore, the application of this peak analysis method (PAM) to a rapid succession of stimulus on/off switches which results in a zigzag-like ΔT(t), which we term the zigzag protocol, allows evaluation of possible variations of the SLP values with time or temperature.

---

1 Introduction

The use of nanomaterials to generate heat in a controlled manner through remote activation using external stimuli has resulted in the design and development of new therapeutic approaches for cancer treatment.¹ Increasing the temperature of tumour cells can lead to their death, thereby removing or shrinking the tumour. To achieve such an increase in temperature, magnetic hyperthermia (MH), based on the use of magnetic nanoparticles (MNPs) and alternating (AC) magnetic fields, has been largely explored.

This interesting approach, however, fails on a crucial task before its arrival in clinical practice: there is a lack of a standardized methodology to characterize the material heating properties when exposed to the AC magnetic field. One of the main reasons for this is that there have been a multitude of home-made devices generated to record the temperature variation of a suspension of nanoparticles over time when exposed to the AC magnetic field.

The heating efficiency of these materials when exposed to an AC magnetic field is quantified in terms of the specific loss power (SLP). In addition to SLP, the heating properties of the MNPs are also sometimes described by a different parameter called the intrinsic loss power (ILP). The term SAR (specific absorption rate), used very extensively in the literature, should be kept to describe power dissipation in tissues.² Nevertheless, the alternative approach presented here can be applied for calorimetric measurements of such properties, generally referred to either as SAR or SLP.

In any case, determining the SLP value of a specific material with accuracy and low uncertainty is not free of difficulties as there are many sources of error in this type of characterization (e.g. field inhomogeneity, location of the temperature probe, heat losses, failure to reach an equilibrium temperature before the beginning of the measurements, etc.).^2–4 The most extensively used approach to determine the SLP value of a magnetic nanoparticle suspension is the use of calorimetric methods using the temperature variation (ΔT) vs. time (t) curve.⁵ However, when comparing SLP values obtained by different laboratories problems start to arise, as the measurement of the heat released constitutes in itself a rather complicated task to reproduce.^6,7 In fact, a recent study done in 21 different laboratories reported large variations between laboratories in the heating capability of a single batch of particles.² Such discrepancies originate from differences in the measurement setups, the analysis techniques and the field conditions used in each laboratory.²

In the context of magnetic hyperthermia, most of the devices designed to measure the temperature variation ΔT(t) curve when the MNPs are exposed to the AC field are non-adiabatic.² Therefore, the heat losses that appear during the measurement may be significantly different depending on the design of the device. Moreover, several works have already described that different mechanisms of heat losses can coincide within a given setup⁸ and can have different timescales.⁹

As a result of all the discussed problems, researchers working in the field of magnetic hyperthermia still lack a reliable and precise method to accurately determine the SLP value of a given particle suspension. Therefore, the development of an alternative approach, less dependent on the measurement devices and able to unify how the SLP values are calculated in a precise and reproducible way, becomes critical.

The objective of the current work is to present a measurement protocol that diminishes the dependence on the specific device characteristics and environmental conditions. The work is presented in four sections. First, the theoretical framework behind the usual calorimetric methods is outlined, followed by a review of the most widely used data analysis approaches to determine SLP. This first section also includes a comparison of experimental measurements and SLP data calculated from measurements performed in three different laboratories, as an illustrative example. The second section describes the origin of the problems that affect the SLP determination using current methods, focusing on the coexistence of various heat loss mechanisms and the inhomogeneous heating of the sample. We have performed a combined experimental/theoretical effort aimed at (i) differentiating effects attributable to the particles themselves from those defined by the device thermal properties, and (ii) disentangling overlapping heat-loss effects on the determination of the heating performance, so that their roles may be understood and minimized. The third section is devoted to the description of the new protocol to determine SLP. The proposed protocol, which we refer to as the “zigzag protocol”, is based on a set of repeated short time heating–cooling cycles and the subsequent analysis of the peaks arising when the alternating magnetic field is switched off. We show that it is more beneficial to shift the SLP determination from the initial time of the heating curve, to the transition between heating and cooling. This is because the difference in losses during the heating and cooling phases is minimised close to the alternating magnetic field on/off transition, and this allows the determination and subtraction of the correct heat loss contribution. Therefore, a more precise determination of SLP values is obtained. The final section provides validation of the proposed SLP calculation methodology. The theoretical validation includes the numerical generation of a test case where SLP values vary over time, to calculate the error in the SLP value. We then validate this approach using as an example an experimental inter-laboratory comparison in which three different devices are used to characterize the same magnetic nanoparticles, showing how the differences between the SLP values estimated using the most common standard protocols are diminished when the proposed zigzag protocol is applied.

2 The basics

2.1 Newton's law of cooling

In general, to have a correct description of the temperature evolution during the hyperthermia process, it is necessary to consider the detailed temperature profile both in time and space. This requires solving the heat diffusion equation including the sample (heat source region), the container and the surrounding environment:where is the temperature, and provides the time (t) evolution of the temperature. ρ is the density, c the heat capacity, k is the thermal conductivity and S is the heat source term. The subscript “r” indicates that eqn (1) needs to be considered in each point in space.

The existing calorimetric methods used to determine the SLP value, which we will refer to as the “classical models”, are based on simplified temperature dynamics where, rather than using the fundamental heat equation (eqn (1)), the temperature evolution is defined in terms of Newton's law of cooling as:¹⁰

where a is a phenomenological heat loss coefficient of the system, S defines the heating source (≡SLP) and T_amb is the environmental temperature. Newton's law of cooling is a simplified version of the heat diffusion equation and the validity can be checked based on the Biot number.¹¹ In this case, the temperature T is considered constant inside the sample. The validity of this approach is, however, limited as will be discussed later.

Note that eqn (2) requires the temperature within the heat source region to be uniform. In other words, under this assumption, the complex temperature profile inside and around the sample (eqn (1)) is replaced by just two temperatures: one for the sample, T, and a second one for the surrounding medium, T_amb. Assuming that both the heat source (S) and the heat loss coefficient (a) are constant over time, the SLP value should be obtainable from the time evolution of the heating curve, which from eqn (2) results in:

where:

In eqn (4), ρ and c stand for the density and specific heat, respectively, of the nanoparticles (NP) and the dispersion medium (w stands for water), and c_vol corresponds to the volume fraction (%) occupied by the particles.

2.2 Usual methods to determine SLP: diversity of results

Under the standard experimental procedure, the initial sample temperature is in equilibrium with the environmental temperature. Since , and S is uncorrelated from the thermal losses of the system (namely the “a” parameter), determining SLP from the initial slope of the ΔT(t) curve (ΔT(t) = T(t) − T_amb) would allow, in principle, the intrinsic SLP value to be obtained while neglecting the role of the thermal characteristics of the measurement device. Such an approach constitutes the so-called initial slope method (ISM), where the SLP can be computed from the linear fit to the initial ΔT(t) curve. In this approach the adiabaticity assumptions imply that, at the initial time, heat losses can be neglected. Such estimates are subject to huge uncertainty due to dynamic changes in the heat loss mechanisms inside the sample and even from the measurement device/environment. In other words, both the characteristics of the experimental set up and time frame in which the slope of the ΔT(t) curve is calculated will have a strong impact on the SLP estimation, leading to unreliable values also with very high uncertainties associated.

In order to solve some of the problems associated with the ISM, other alternatives based on a more complex analysis of the initial slope have been proposed. A detailed discussion on such methods can be found in Wildeboer et al.³ Briefly, the Box–Lucas method (BLM) uses a different equation to adjust the initial slope.¹² Alternatively, the Corrected Slope Method (CSM) analyzes several time intervals of the initial slope and uses a correction for an estimated linear heat loss to extract the SLP data. To illustrate how the use of these different data analyses affects the calculation of the SLP value, we have analyzed a suspension of magnetic nanoparticles (dextran coated, 32.0 ± 6.7 nm average size, see Section 1 of the ESI† for further characterization details) using the ISM, BLM and CSM approaches; see Fig. 1. In addition, the samples were measured in three different devices operating under very similar conditions (163.3 kHz and 35 mT for device 1, 165 kHz and 35 mT for device 2, and 172.4 kHz and 35 mT for device 3, in all cases using a suspension of 1 mL with an iron concentration of 1 mg mL⁻¹). For a given device, significant differences, up to 16%, were observed in the SLP values calculated by the different methods (see specific values in Fig. 1). When comparing devices, SLP values calculated from device 2 were significantly larger than those obtained with the other two devices, the lowest obtained value (194 W g⁻¹ Fe₃O₄) being for the CSM with device 1 and the highest for the same method in device 2 (294 W g⁻¹ Fe₃O₄), which corresponds to a 34% difference in the obtained values. From this analysis, it is clear that there is a large uncertainty in the SLP value associated with both the device used for the measurements and also the method selected for the data analysis.

Fig. 1 Description of the different SLP-determination analyses, with three different measurement devices located in three different laboratories. Measurements were performed using a 1 mL suspension of the same particles (dextran coated, 32.0 ± 6.7 nm average size) at a concentration of 1 mg_Fe mL⁻¹. AC field conditions were 163.3 kHz and 35 mT for device 1, 165 kHz and 35 mT for device 2 and 172.4 kHz and 35 mT for device 3. The solid lines are the fitting curves. Details of the calculations can be found in Section 2 of the ESI.†

In addition to the data analysis methods focusing on the analysis of the initial slope of the heating curve, an alternative approach is the Decay Method (DM), that includes in the sample characterization recording both the heating curve, when the alternating magnetic field is switched ON, and the cooling curve, once the alternating magnetic field is switched OFF. Observing the different curves displayed in Fig. 1, where small changes (and even abrupt jumps) in the data may lead to significant changes in the slope, one may easily imagine that extending the fitting range would help minimising the dependence on the specific features of the curves. Furthermore, given that there are only two unknown values in eqn (1), it might seem reasonable to devise two different scenarios to fit two similar curves, so that we have two equations with two unknowns. Therefore, the decay method uses the cooling phase to obtain a characteristic time of the system cooling down and the steady state temperature. However, from an experimental point of view, this method requires longer characterization times, to ensure that the steady state temperature has been reached and also to characterize the cooling phase; therefore, this approach is much less frequently used than those described earlier. However, it is not just a matter of time: as will be explained in the next section, the fit of the cooling curve is not straightforward and, thus, the use of the decay method to obtain the SLP value also has associated systematic errors.

3 Identifying the problems – applicability of Newton's law

The applicability of Newton's law of cooling has been investigated by Vollmer.¹¹ Vollmer states five conditions for the validity of Newton's law of cooling. Three of these relate to the nature of the cooling process: for example, in the case of radiative transfer the temperature difference must be small enough that the heat transfer can be approximated by its linearised form, also the convective heat-transfer coefficient must stay constant during the cooling process. In the context of the current work, (1) it is also assumed to be the case that the only internal energy source of the object is the stored thermal energy, which of course distinguishes the heating and cooling phases of the hyperthermia measurement, and (2) it is shown that Newton's law of cooling applies only if the object is characterized by a single temperature. It is important to note that in ref. 2 it was found that in 6 out of 17 laboratories Newton's law was demonstrably not valid. This has also been recently analysed in detail by Hanson et al.¹³ In the following we examine these factors in the context of measurements of hyperthermia. We use the general expression for heat flow, eqn (1), including a generic source term which is appropriate for both magnetic hyperthermia and photothermia. We demonstrate numerically that a measurement protocol, the Peak Analysis Method (PAM), can successfully be used to measure SLP values even when Newton's law of cooling is demonstrably not valid. Departure from Newton's law of cooling is characterised by the Biot number,¹¹ which describes the heat transfer rate inside and outside the measurement vial.

Aiming to elucidate what might be the underlying reasons for the diversity in SLP values, not only between different devices but also when applying different protocols to the same data ΔT(t), in this section we perform a detailed analysis of the underlying physical background: the assumption that Newton's law of cooling is applicable. Thus, we analyse in detail the two implicit simplifications assumed in deriving eqn (3) from eqn (2). The first is that there is only one heat-loss channel, defined by the “a” parameter. The second is that the temperature of the sample is homogeneous. The differences between the ideal situation assumed in Newton's law and the real system are schematically depicted in Fig. 2.

Fig. 2 Schematic illustration of the differences between the ideal scenario under which Newton's law of cooling would be applicable (i.e. homogeneous temperature, and one heat-exchange mechanism; left panel); and the real situation in which the sample temperature is not homogeneous, and heat can be exchanged through different channels (right panel).

3.1 Assumption 1: a single heat-loss channel

The first simplification performed to reach eqn (3) from eqn (2) assumes that there is only one heat-loss channel, defined by the “a” parameter. As explained in the previous section, the decay method uses the cooling phase (where SLP = 0, as no alternating magnetic field is applied) to obtain the phenomenological heat loss coefficient “a”, and that value is taken to obtain the SLP from the heating phase. Since Newton's law of cooling, eqn (2), corresponds to a single exponential decay, the exponential fit of the cooling part would provide the “a” value that could be used to fit the heating part. However, as discussed in the insightful work by Landi,¹⁰ it is often observed that a single exponential fit does not match the experimental cooling part. This is illustrated in Fig. 3, where a single “a” value to fit the cooling part of the ΔT(t) cannot be obtained in a measurement performed using the same particles described in the previous section. In fact, depending on how the fit of the cooling part is performed, very different “a” values are obtained (see the caption of Fig. 3 and Section 2 of the ESI† for details). The lack of a single “a” value, indicative of the presence of several heat-loss mechanisms, is more clearly emphasized when plotting vs. ΔT, which according to eqn (2) should be a straight line if a single heat loss mechanism is occurring. As the inset in Fig. 3 shows, this is clearly not the case. The heat loss mechanisms and the corresponding “a” parameter are highly influenced by the device setup and environmental conditions, which can differ significantly from one laboratory to another. The presence of several heat-loss mechanisms has also been recently reported by Iglesias et al.⁸ In summary, several heat loss mechanisms may simultaneously occur during the calorimetric methods to determine the SLP value, being a source of the variable results reported in the literature.

Fig. 3 Single heating and cooling cycle measurement of the same particles from Fig. 1. Three different approaches to fit the cooling part using a single exponential have been performed. The blue fit (dotted line) was performed using the initial time part of the cooling curve, obtaining a = 0.0020. The red fit (dashed line) was a free fit to the entire cooling part, with a = 0.0015. The green fit (continuous line) was obtained from the final part of the curve, with a = 0.0017. (Inset) Representation of the vs. ΔT of the same data. If a single heat loss mechanism was occurring, the data should be a straight line according to eqn (2). Results indicate that there is not a single heat loss mechanism. Straight lines are a guide for the eye.

3.2 Assumption 2: homogeneous temperature

The second simplification performed to reach eqn (2) from eqn (1) is that the temperature of the sample is homogeneous. Essentially this requires the temperature within the vial to be uniform, which has been observed experimentally to not be the case, as temperature gradients are clearly observed when using infrared cameras to monitor the whole volume of the sample.^14,15 This factor was investigated in ref. 11 and quantified using the ‘Biot number’ Bi, which is a dimensionless quantity, usually describing the ratio of two adjacent heat-transfer rates, in this case inside and outside the vial. In the context of the current investigation we use the definition Bi = k_wall/k_fluid.

We investigate this non-uniform heating within the sample by simulating the heating and cooling processes using a simple 1D heating model (Fig. 4). The reality of heat loss processes in practical experimental set-ups is complex, possibly involving more than one heat-loss process, and is difficult to model. Our approach was intended to illustrate the heating and cooling processes using a simple and physically transparent model, which is successful in highlighting the need for an advanced measurement protocol and evaluating its likely efficacy prior to experimental validation. We modeled the increase/decrease of temperature with the SLP driving term on/off. The conjecture was as follows: when the alternating magnetic field is ON, the SLP is localised giving rise to a rapid (uniform) heating with losses mainly through the boundary of the sample holder, whereas, when the alternating magnetic field is OFF, during the cooling phase, there will be a slow migration of heat out of the system due to small gradients within the sample holder. In 1D the time variation of the temperature (with ΔT = T − T_amb) is given by:

where α = k/(ρC) is the thermal diffusivity, and for simplicity of notation will be used. Eqn (5) is solved numerically with interface conditions between the fluid and the vial corresponding to heat loss by conduction and by convection. These conditions are given in the Methods section and in Section 3 of the ESI.†

Fig. 4 Plots of the temperature profile during heating and cooling. Vertical lines show the extent of the heated fluid. (A) High thermal conductivity of the vessel (k_wall = 0.8 W m⁻¹ K⁻¹) and (B) low thermal conductivity (k_wall = 0.01 W m⁻¹ K⁻¹) of the vial. The effect of improved insulation as expected leads to a larger temperature rise, and additionally a lower temperature gradient within the vial. The red and orange lines correspond to the heating process, with the orange line corresponding to the end of the heating process, whereas the blue colour corresponds to the cooling process. The S parameter corresponds to an SLP of 1000 W g⁻¹.

Results from this analysis can be found in Fig. 4, which shows plots of the temperature profile during heating and cooling. Vertical lines show the extent of the heated fluid. Data are shown for high (k_wall = 0.8 W m⁻¹ K⁻¹) and low (k_wall = 0.01 W m⁻¹ K⁻¹) thermal conductivity of the vial. As expected, it can be observed that the temperature gradients occurring within the liquid sample differ significantly between the heating and cooling phases as a consequence of the heat flow and heat loss processes. Moreover, the impact of the thermal conductivity values of the vial has also been tested (see Fig. 4). In the case of high thermal conductivity of the vial the temperature gradient varies strongly across the heated fluid and also varies considerably with time. As expected the increased insulation leads to a higher temperature rise. It also leads to profiles with rather less curvature. These results are consistent with the analysis of ref. 11. Specifically, in panel A of Fig. 4 the corresponding Biot number is Bi = 1, which is consistent with the observed strong temperature variation, whereas for panel B the Biot number is Bi = 0.0125, which is in the region where weak temperature variations are expected.

These results illustrate the fact that heat losses and temperature gradients within the sample can be strongly dependent on the details of the experiment and in particular the thermal properties of the container. This clearly indicates that the validity of Newton's law of cooling cannot be generally assumed to be true. However, in the following we show, using the general heat flow equation, that it is possible to define a peak analysis method, under which protocol it is possible to estimate heat losses and precisely compensate for the calculated SLP even when Newton's law of cooling does not apply. The protocol is validated by experiments in 3 different measurements in 3 laboratories.

4 A new approach to calculate the SLP

In order to solve the problems associated with the presence of temperature gradients and the co-existence of simultaneous heat loss mechanisms described in the previous section, we have generated an alternative protocol to perform the calorimetric experiments and the data analysis using which it is possible to significantly diminish the role of the heat-loss mechanisms associated with the different devices. This protocol is based on the following principles. First, there is a shift in the data that is going to be analyzed, from the initial slope measured via the classical methods, to the peak generated when the field is switched off. Second, instead of doing a single measurement, a series of on/off switches of the alternating magnetic field are performed resulting in a zigzag shape ΔT(t) curve. The reasons behind these two modifications are explained in detail below.

4.1 Change of paradigm: the peak analysis

The main reason to shift the data analysis from the initial slope to the peak generated by the alternating magnetic field on/off switch is that if we were able to measure the cooling curve immediately after heating, the effect of losses at the transition between heating and cooling phases and also the coexistence of longer relaxation time loss mechanisms on SLP determination would be minimised. Further, as we show later, by measuring the decay curve it becomes possible to correct accurately for the heat losses thereby removing a significant systematic error in SLP measurements. This is illustrated in Fig. 5, where the 1D model described in the previous section has been used to plot the gradient along the sample at some specific time points during the heating and the cooling parts during the ΔT(t) measurement. We show that, when we are looking at data far away from the peak, the temperature profiles are very different at a given temperature during the heating and cooling phases. In contrast, when we are closer to the peak, these gradients are similar and the form of the temperature profile is essentially unchanged around the peak temperature.

Fig. 5 Plots of the temperature profile during a single heating and cooling experiment. The vertical lines show the extent of the heated fluid. Near the peak, the gradient in the heating curve and the cooling curve is similar. Further away from the peak, when the temperature of the suspension is 9.5 °C, larger differences in the temperature gradients inside the vessel are observed.

In fact, we can write the heat diffusion equation for the two processes: (1) heating when the alternating magnetic field is ON:

and (2) cooling, when the alternating magnetic field is OFF:

A reasonable assumption, verified by the simulations (see the discussion around Fig. 5 and 6), is that at the transition between heating and cooling the spatial derivatives will be very similar, and thus subtracting eqn (6) from eqn (7) leads to as precise as possible a determination of S and implicitly SLP:

Fig. 6 Comparison between single heating/cooling cycle (A) and repeated heating–cooling cycles (B). The plots show the curvature of the temperature profile inside the heated fluid. The negative sign indicates the concave profile of temperature and the value indicates the degree of homogeneity of the temperature profile. The inset in (a) is the curvature as a function of temperature for alternating magnetic field ON (red: heating) and alternating magnetic field OFF (blue: cooling). This shows a complex temperature evolution inside the fluid both during the heating and cooling process, which is not symmetric. The results correspond to case A from Fig. 4.

Therefore, shifting the data analysis from the initial slope to the peak seems a feasible way of avoiding the problems associated with the temperature gradients that originated within the sample during the heating and cooling processes. Moreover, this methodology will compensate for uncertainties associated with the temperature probe position.

4.2 Repeated alternating magnetic field ON/OFF switches

The peak analysis proposed in the previous section could be done for one heating–cooling cycle or multiple heating–cooling cycles. Other groups have already proposed the use of a “stepped heating procedure”¹⁶ to acquire heating and cooling data at different temperatures or to verify the absence of parasitic signals from the thermocouple.¹⁷ In addition, the idea of applying the alternating magnetic field intermittently has also been recently described as a way to control the global temperature reached in tissue phantoms.^18,19

Our proposal is that, given that the data needed for the peak analysis described in the previous section corresponds only to the values closest to the peak, the temperature profile remains essentially constant and the SLP value can be determined using eqn (8). It follows that there is no need to wait for the temperature to reach an equilibrium (neither in the heating, nor the cooling part), and a sequence of fast cycles of ON/OFF alternating magnetic field can be performed. This approach will provide two main advantages. First, the calculation of SLP values in several “peaks” will allow calculation of the error associated with the SLP determination in repeated measurements faster than repeating “classical methods” several times. Furthermore, it will allow tracking of any possible changes in SLP values due to differences in the global temperature, also providing a tool to calculate the SLP at body temperature. To verify the use of this pulsed alternating magnetic field approach, we have quantified the degree of temperature homogeneity inside the sample (as illustrated in Fig. 4). The curvature of the temperature profiles as a function of temperature is presented in Fig. 6. We calculate a radius of curvature as given in the Methods section. The characteristic radius of curvature is calculated at the center of the heated vial.

We have compared the temperature profile when performing a single heating and cooling cycle and the zigzag protocol using the simple 1D model outlined earlier. Because the radius of curvature is infinite at time t = 0 we choose to characterise the curvature as the inverse of the radius of curvature calculated as given in the Materials and methods section (eqn (14)). In the single heating/cooling cycle, it can be observed that there is a significant variation of the curvature during the heating process (red line in Fig. 6). Additionally, a more complex variation of the curvature with time occurs when the field is turned off (the cooling phase). We can observe that the behaviour is not symmetric for heating and cooling, supporting our hypothesis that the heat-loss mechanisms are complex and cannot be assumed to be the same during the two phases. As expected the curvature converges at the transition between heating and cooling, thereby supporting the use of the peak analysis. Indeed, a more linear behaviour of the curvature is observed during the simulation of the zigzag protocol, where the convergence between the heating and cooling parts of the curve at the peak is clear. These results support the idea of using the peak analysis when applying a pulsed alternating magnetic field to the sample.

4.3 The zigzag protocol

We have termed our proposal for the new approach to determine the SLP value the zigzag protocol (Fig. 7). This proposal comprises two complementary ideas. The first is shifting the data analysis from the initial slope to the peak. As described above, the peak analysis will allow the ideal limit to be studied where the heat loss mechanisms in the heating and cooling phases of the ΔT(t) curve approach a single value, providing a scenario that allows the determination of an SLP value independent of the device. Moreover, performing a zigzag type of measurement, where the alternating magnetic field is intermittently switched ON and OFF, it is possible to determine the error associated with the SLP value using different fast repetitions of the peak analysis, assuming a constant value of the SLP. A Standard Operating Procedure to perform this analysis has been developed and is provided in Section 4 of the ESI.†

Fig. 7 Scheme depicting the advantages of the use of the peak analysis method and the zigzag protocol. Shifting the data analysis to the peak resulting from the ON/OFF alternating magnetic field switch allows assuming similar heat loss mechanisms during the heating marked in red in the top panel and the cooling phases marked in blue in the top panel. Repeating the on/off switches several times in a zigzag way allows replicates to be quickly obtained and the evolution of the heating capacity of the particles to be tracked over time and temperature.

4.4 Time dependent SLP

We have shown theoretically that the PAM method provides determination of the SLP with reduced systematic errors introduced by thermal transport processes in comparison with the existing methods. In addition to reducing the errors in the SLP calculation, the use of the zigzag protocol has other interesting advantages. In the first instance, under the assumption of a constant (time invariant) SLP, the extra data provided by the zigzag method allows further averaging and a reduction in the statistical error. Further, it could also be used to accurately determine a time varying SLP. This might arise from temperature variations in the intrinsic nanoparticle properties; for example, in the case of magnetic hyperthermia, the use of low Curie temperature materials to control their heat output would give rise to a time (temperature) varying SLP.^20,21 A further example is the possibility of time dependent chaining as observed by Mille et al.²² when using magnetic nanoparticles. To illustrate this aspect, we simulated a case when SLP increases exponentially in the form SLP = SLP₀[1 − exp(−t/t₀)]. The results for SLP₀ = 1000 W g⁻¹ and t₀ = 100 s are shown in Fig. 8. As the transition between heating and cooling happens over a short period of time, we can apply the PAM method for the zig-zag protocol to obtain the SLP at each time corresponding to the transition (as shown at the bottom of Fig. 7) meaning that we can obtain the SLP at various time points as illustrated in Fig. 8. As the temperature evolution was simulated, the exact SLP time dependence is known, and this allowed us to calculate the absolute error of the PAM method at different time points. These results prove the capability of the zigzag method to accurately determine SLP, not just in the case of constant SLP but also for the case where SLP varies during the MH measurement. We believe that this feature represents a significant advantage of this protocol. In this way, any changes in the SLP due to particle clustering, or temperature dependence of magnetic properties could be tracked over time. However, it should be noted that this will require high accuracy (low noise) measurements which should be possible by using high resolution temperature probes and improved averaging.

Fig. 8 Example of time dependent SLP (black line, input SLP) and the values extracted from the zig-zag method (gold dots). The results based on the error given in the inset indicate a very accurate determination of the input SLP.

In summary, the zigzag protocol not only provides a better estimation of SLP, but can also (1) provide the variation of SLP during the heating protocol and/or (2) the SLP value at the desired operation temperature, which is generally different from the ambient temperature at which the “classical methods” based on the initial slope are generally applied.

4.5 Experimental validation of the “zigzag protocol”

While the theoretical calculations provide insight into the applicability of the proposed methodology, the final step is to proceed to the experimental validation of the protocol. In particular, we use the same particles described in the first section (see also Section 1 of the ESI†). Suspensions of the dextran coated particles are characterized using three different commercial devices for the SLP analysis and a representative measurement is depicted in Fig. 9a and b. In all devices, a suspension of 1 mL with an iron concentration of 1 mg mL⁻¹ prepared from the same batch is used. In all devices, very similar AC field conditions are used (163.3 kHz and 35 mT for device 1 (D5 Series from nB nanoScale Biomagnetics), 165 kHz and 35 mT for device 2 (Fives Celes, MP 6 kW) and 172.4 kHz and 35 mT for device 3 (Nanotherics)). The sample was measured in each device using two experimental protocols: a single heating/cooling cycle and the zigzag protocol, where several cycles of faster heating/cooling cycles are sequentially performed following the Standard Operating Procedure provided in Section 4 of the ESI.† Measurements were repeated twice for the single heating/cooling measurements. Repetitions of the zig/zag measurements were also performed to assess the reproducibility of the data obtained in selected cases (see Section 5 of the ESI†). Additionally, water samples were also tested to evaluate the possible heat contributions from the coils (see Section 6 of the ESI†).

Fig. 9 SLP calculation from experimental calorimetry measurements using three different devices. (A) Representative measurement of the classical single heating/cooling cycle. (B) Representative measurement of the repeated heating–cooling cycles (zigzag measurements). (C) The ISM, CSM, and BLM approaches have been used to calculate the SLP values from single heating/cooling measurements, as the ones shown in panel (A). SLP analysis has been performed in the same time range shown in Fig. 1. Values correspond to the mean and standard deviation calculated from two replicas. The Peak Analysis Method (PAM) has been applied to each of the peaks from the zigzag measurements shown in panel (B). (D) Comparison between the extracted SLP values for the three devices when using different methods for the data analysis. Note that SLP values for the PAM method correspond to the average value for all the peaks analyzed independently and shown in panel (C). Also note that the sample size is different for the ISM, CSM, and BLM approaches (2 replicas) than for the zigzag measurements (10 peaks).

When considering the single heating–cooling cycle, measurements performed with the three devices revealed a fairly similar initial slope. It is important to note that although the initial slope may look similar, when the initial time was analyzed, significant differences appear, as shown in Fig. 1. However, differences in the cooling phases associated with different degrees of insulation and therefore different thermal losses were evident (Fig. 9a). Devices 2 and 3 have a more uniform cooling process, whereas device 1 has clearly at least two timescales associated with the cooling. This different degree of insulation was also observed when measuring the water samples (see Section 6 of the ESI†) as external heating arising from the non-adiabaticity of the devices was observed.

SLP values were calculated using the classical methods applied to the first 30–40 s range, as shown in Fig. 1. As described in the first section, striking differences in the SLP values, up to 35%, are observed in the calculations performed using the different classical methods applied to the initial slope analysis (see Fig. 1). A repetition of the measurements was performed in order to test the variations associated with the SLP value calculation depicted in Fig. 9c. Two measurements were performed in each device for the classical single heating/cooling measurements and each measurement was analyzed using different “standard models”. These results clearly show the large uncertainty in the SLP value determination of a material that arises from both the device used for the measurements and also the method selected for the data analysis. Section 2 of the ESI† shows the individual SLP values calculated for each measurement and type of data analysis.

Then, the zigzag protocol is applied to measurements with the three different devices (Fig. 9b). At least ten peaks are analyzed for each device (see individual SLP data obtained for each peak in Section 2 of the ESI†). Average values, standard deviations and relative standard deviations (RSD) are also provided in that section of the ESI.†

Fig. 9c shows the average data of the SLP values obtained from several measurements of single heating/cooling cycles and calculated using the classical data analysis approaches (ISM, BLM and CSM) applied to the initial slope of the results. When comparing relative standard deviations obtained for the different calculation methodologies, interesting differences are observed. In the case of the single heating/cooling measurements, relative standard deviations are much higher for the SLP values calculated for device 1 (RSD between 15 and 26%) than for the other devices (RSD between 2 and 11% for devices 2 and 3) independently of the calculation method used. However, given the large variance of some of the data sets, the statistical analysis, comparing if any device was providing a significantly different value than the others, rendered negative results.

SLP values obtained from the single heating/cooling measurements are compared to the individual SLP values obtained using the peak analysis method for each of the sequential peaks generated in the zigzag protocol over time (see Sections 2 and 7 of the ESI† for the individual values). Although not used in this particular case, it can be noted that, as shown earlier, this approach allows the possibility of tracking SLP variations over time. It also allows the characterization of the SLP value of a material occurring at a given temperature. In general, results from the peak analysis show a much smaller variation of the SLP values obtained in all the peaks for device 2 when compared to the other two devices. We think the reason behind this observation is the better temperature resolution (0.01 K) of the thermal probe of device 2, as can be observed also from the temperature data in Fig. 9a. Nevertheless, the zigzag protocol can detect this aspect, which can be taken into account when comparing results between different devices and laboratories.

The average SLP value calculated from all the analyzed peaks in the zigzag protocol is shown in Fig. 9d and compared with the values obtained from the classical approach that focuses on the initial slope. A good consistency among the obtained results from the three different devices is obtained for the zigzag protocol, especially compared with the results obtained from the single heating/cooling approaches. Here it is important to note that the sample size for the classical methods is different to that of the zigzag approach, as the ten peaks were averaged.

One of the advantages of the zigzag method is that it allows the relatively fast measurement of a great number of SLP values, helping to reduce uncertainties. Indeed, when the average SLP was calculated from the analysis of at least 10 peaks using the zigzag protocol, uncertainties were highly reduced for device 1, the device showing the larger standard deviations in the classical methods. It should be noted that to perform a complete single heating and cooling measurement and to allow the sample to cool down completely, in this case, around 50 min are needed. In contrast, using the zigzag method, between 10 and 14 peaks can be analyzed in the same time frame. This means that the time needed to achieve a high enough number of replicates to reduce uncertainties can be reduced dramatically (10 fold).

Although the absolute error on the average SLP value cannot be determined given the absence of reference materials or verified methods, the similarity of the results obtained when measuring the SLP in three very different devices is a promising result. This good consistency of the results is especially relevant provided the disparity of results obtained for the single heating/cooling measurements. It is especially interesting to consider the case of device 3. This device presented the lowest relative standard deviation in the SLP values in the data analysis using the classical methods (RSD < 5%); however, the average SLP values obtained from such methods are lower (230–240 W g⁻¹ with 5% RSD) than those obtained from the zigzag protocol (270 W g⁻¹ with 8% RSD). Nevertheless, although these results are promising, future work should be performed to further validate this approach by research groups with adiabatic devices.²³

Overall, these experimental results support the idea of using the zigzag protocol as an interesting alternative for the SLP analysis independently of the device being used.

5 Conclusion

We have analysed the contribution to systematic deviations between measurement systems arising from the different thermal properties of the systems and sample holders. A simple 1D conduction model was used to illustrate the central factor, specifically the difference between the temperature profile during heating and cooling. This essentially arises from the fact that the heating is rather uniform, resulting in a relatively flat temperature profile in the vial during heating which becomes more strongly curved during cooling. Consequently, Newton's law of cooling is not generally obeyed. However, around the peak temperature, when the alternating field is switched off, the temperature profile remains essentially constant allowing for exact calculation of SLP by compensating for the heat losses obtained from the cooling curve. Therefore, an alternative way of calculating the SLP value was proposed, called the Peak Analysis Method (PAM). This method was further extended into a ‘zigzag protocol’ by carrying out repeated heating/cooling cycles. The method was verified using the simple 1D temperature diffusion model and then experimentally validated using 3 different measurement systems in different laboratories. It was shown experimentally to give significant improvement to achieve good reproducibility within a given time frame over the ‘classical’ methods of SLP determination. Finally, we note that the zigzag method in principle allows the determination of a time dependent SLP as can arise, for example from the effects of chaining.²² However, this would require more accurate measurements of heating curves. Nonetheless, the PAM technique and the zigzag protocol have shown a remarkable agreement between different measurement systems and laboratories and is recommended as a significant improvement over current measurement techniques.

6 Materials and methods section

6.1 Magnetic nanoparticle synthesis and characterization

MNPs were synthesized using an oxidative precipitation aqueous route previously described²⁴ with slight modifications. Briefly, a 1 M solution of FeSO₄, prepared in 50 mL of 0.01 M H₂SO₄, was quickly added to a basic solution prepared with 4.25 g of NaNO₃ and 4.22 g of NaOH in a mixture of 137 mL of water and 63 mL of 96% vol ethanol. The green rust suspension obtained was stirred for 15 min and poured in a jacketed flask previously thermalized to 90 °C with a thermostatic bath. The MNPs were left to grow at that temperature for 6 hours. MNPs prepared by this route were subjected to an acid treatment and then coated with dextran.²⁵

Size and morphology were studied by transmission electron microscopy (TEM). A drop of the diluted sample was deposited on a carbon coated grid, allowing it to dry at room temperature. Micrographs were acquired in a Tecnai G2 TEM (FEI) operated at 200 kV. The particle size was defined considering the largest internal dimension of the nanoparticles. A total of 175 nanoparticles were manually measured and the histogram obtained was fitted with a probability density function.

Magnetic characterization of the sample was carried out in a Quantum Design MPMS-XL SQUID magnetometer. The liquid suspension was placed in a cotton piece allowing it to dry. Then, this piece of cotton was placed in a gelatin capsule for the magnetic measurements. Field dependent magnetization of the sample was recorded at 300 K at a maximum field of 1600 kA m⁻¹.

The magnetic hyperthermia measurements were performed using three different devices. Device 1 is a commercial equipment (D5 Series from nB nanoScale Biomagnetics) with a G model closed coil and a fiber optic sensor. This device uses, as sample holders, 2 mL disposable glass vials. Device 2 is an AMF produced by a Fives Celes 12, 118 M01 generator. This device is composed of a combination of a CELES MP 6 kW generator capable of generating resonant frequencies in the range 100–400 kHz (tunable with an ALU CU type capacitor box) and a 71 mm i.d. DT25901A chilled coil. Temperature was measured with an OSENSA fiber optic probe model PRB-G40-02 M-STM-MRI. This device uses, as sample holders, plastic microtubes. Device 3 is a MagneTherm from Nanotherics that uses 1 mL Nunc® CryoTubes® as sample holders. Two types of measurements (classical single heating–cooling cycle and a zigzag measurement) were performed in the three devices using a suspension of 1 mL with an iron concentration of 1 mg mL⁻¹ and similar AC field conditions (163.3 kHz and 35 mT for device 1, 165 kHz and 35 mT for device 2 and 172.4 kHz and 35 mT for device 3). For the single heating/cooling measurements, 1 mL of the particle suspension was placed into a specific container and located at the center of the magnetic induction coil inside an isolating holder. When the sample temperature was stable, the AC magnetic field was applied for 110–120 s. The sample temperature during the heating and cooling time was measured in each device at least 2 different times. For the zigzag heating–cooling measurements, the suspension was placed in the center of the magnetic induction coil inside an insulating support. The sample temperature was stabilized before starting the measurement. The AC magnetic field was applied until the sample temperature increased 2–3 °C and then the AC field was switched off. The magnetic field was turned on again when the sample temperature decreased by 2 °C. This process was repeated several times. Even though the measurements were performed by different operators, data analysis was performed by the same person. Furthermore, several quality control tests were performed. Reproducibility was assessed by repeated measurements of the same sample in each device (see Section 5 of the ESI†). The negligible contribution from the coils to the temperature of the sample was verified by measurements of water samples of the same volume in each device (see Section 6 of the ESI†).

6.2 Simple (1D) model of heating

The aim is to model the increase/decrease of temperature with the SLP driving term on/off. The conjecture is as follows: the SLP is localised giving rise to a rapid (uniform) heating with losses mainly through the boundary of the sample holder, whereas during the cooling phase there will be a slow migration of heat out of the system due to small gradients within the sample holder. In 1D the time variation of the temperature (with T = T − T_ambient) is given bywhere the SLP S is taken as a constant, ρ is the density, c the heat capacity and k is the thermal conductivity. α = k/(ρC) is the diffusivity. The vial is uniformly heated with SLP S for a time t_heat after which s is set to zero for the cooling phase. At the edge of the simulation region is assumed that the temperature equals the environment.

6.2.1 Boundary conditions. We now look at the experimental case of heating in a vial. At the boundary we have to impose continuity of the temperature and the heat flux, with the interface condition

In the discrete approximation, continuity of T and eqn (10) lead to

T_i = rT_i−1 + T_i+1/(1 + r), (11)

where T_i is the interface temperature, T_i−1 and T_i+1 are temperatures immediately inside the fluid and boundary respectively, and the ratio r = k_wall/k_fluid.

Next we consider the heat transfer to the surroundings via conduction through the vial and by convection from the upper surface of the fluid.

The convection BC is

with h being a constant. This leads to interface temperature at the surface of the fluidwhere Δx is the numerical spatial discretisation.

6.3 The radius of curvature

We characterise the curvature based on a radius of curvature as follows:

A “−” sign indicates the concave temperature profile and the absolute value quantifies how large the non-uniformity is, with zero meaning constant temperature. Note that because the initial radius of curvature at time t = 0 is infinite we chose to use r⁻¹ as a measure of curvature.

Author contributions

Sergiu Ruta: conceptualization; formal analysis; investigation; methodology; software; validation; writing – original draft. Yilian Fernández-Afonso: conceptualization; formal analysis; investigation; validation; visualization; writing – original draft. Samuel E. Rannala: methodology; software. M. Puerto Morales: funding acquisition; investigation; writing – review and editing. Sabino Veintemillas-Verdaguer: investigation. Carlton Jones: investigation. Lucía Gutiérrez: conceptualization; funding acquisition; project administration; supervision; visualization; writing – original draft. Roy W Chantrell: conceptualization; investigation; methodology; project administration; software; supervision; writing – original draft. David Serantes: conceptualization; funding acquisition; methodology; project administration; software; supervision; visualization; writing – original draft.

Conflicts of interest

C. Jones works for Nanotherics, a company that produces devices to characterize the heating performance of magnetic nanoparticles.

Acknowledgements

Authors would like to acknowledge financial support from the following projects: Project PID2021-122508NB-I00 funded by MICIU/AEI/10.13039/501100011033 and FEDER, UE, Projects PID2019-109514RJ-100 and PID2020-13480RB-I00 funded by MICIU/AEI/10.13039/501100011033 and project CNS2023-144321 funded by MICIU/AEI/10.13039/501100011033 and NextGenerationEU/PRTR. Xunta de Galicia is acknowledged for project ED431F 2022/005 (to D.S.). AEI is also acknowledged for the Ramón y Cajal grant RYC2020-029822-I to D.S. The authors would like to acknowledge Fondo Social del Gobierno de Aragón (grupo DGA E15-23R) and the use of the Advanced Microscopy Laboratory (Universidad de Zaragoza), for access to their instrumentation and expertise. We acknowledge the Centro de Supercomputacion de Galicia (CESGA) for computational resources. SR would like to acknowledge the Sheffield Hallam University Beowulf Cluster for computational resources. Authors would also like to thank Lise Grüner Hanson and Cathrine Frandsen for fruitful discussions about this work.

Notes and references

1. R. M. Fratila and J. M. De La Fuente, Nanomaterials for Magnetic and Optical Hyperthermia Applications, Elsevier, 2018.
2. J. Wells, D. Ortega, U. Steinhoff, S. Dutz, E. Garaio, O. Sandre, E. Natividad, M. M. Cruz, F. Brero and P. Southern, et al. , Int. J. Hyperthermia, 2021, 38, 447–460.
3. R. R. Wildeboer, P. Southern and Q. A. Pankhurst, J. Phys. D: Appl. Phys., 2014, 47, 495003.
4. S.-Y. Wang, S. Huang and D.-A. Borca-Tasciuc, IEEE Trans. Magn., 2012, 49, 255–262.
5. E. Garaio, J. M. Collantes, J. A. Garcia, F. Plazaola, S. Mornet, F. Couillaud and O. Sandre, J. Magn. Magn. Mater., 2014, 368, 432–437.
6. P. Lemal, C. Geers, B. Rothen-Rutishauser, M. Lattuada and A. Petri-Fink, Mater. Today: Proc., 2017, 4, S107–S117.
7. F. Soetaert, S. K. Kandala, A. Bakuzis and R. Ivkov, Sci. Rep., 2017, 7, 1–15.
8. C. Iglesias, J. de Araújo, J. Xavier, R. Anders, J. de Araújo, R. da Silva, J. Soares, E. Brito, L. Streck and J. Fonseca, et al. , Sci. Rep., 2021, 11, 1–13.
9. M. Coïsson, G. Barrera, F. Celegato, L. Martino, F. Vinai, P. Martino, G. Ferraro and P. Tiberto, J. Magn. Magn. Mater., 2016, 415, 2–7.
10. G. Landi, J. Magn. Magn. Mater., 2013, 326, 14–21.
11. M. Vollmer, Eur. J. Phys., 2009, 30, 1063.
12. O. L. Lanier, O. I. Korotych, A. G. Monsalve, D. Wable, S. Savliwala, N. W. F. Grooms, C. Nacea, O. R. Tuitt and J. Dobson, Int. J. Hyperthermia, 2019, 36, 686–700.
13. L. G. Hanson, B. L. Hansen, T. Veile, M. Zambach, N. B. Christensen and C. Frandsen, IEEE Magn. Lett., 2023, 14, 6100505.
14. D. Valdés, T. Torres, A. M. Maldonado, G. Urretavizcaya, M. Nadal, M. V. Mansilla, R. Zysler, G. Goya, E. De Biasi and E. Lima Jr, Phys. Rev. Appl., 2023, 19, 014042.
15. A. Espinosa, J. Kolosnjaj-Tabi, A. Abou-Hassan, A. Plan Sangnier, A. Curcio, A. K. Silva, R. Di Corato, S. Neveu, T. Pellegrino and L. M. Liz-Marzán, et al. , Adv. Funct. Mater., 2018, 28, 1803660.
16. N. Iacob, G. Schinteie, P. Palade, C. Ticos and V. Kuncser, Eur. Phys. J. E: Soft Matter Biol. Phys., 2015, 38, 1–6.
17. N. Liu, A. Pyatakov, A. Saletsky, M. Zharkov, N. Pyataev, G. Sukhorukov, Y. Gun’ko and A. Tishin, J. Magn. Magn. Mater., 2022, 555, 169379.
18. A.-R. Tsiapla, A.-A. Kalimeri, N. Maniotis, E. Myrovali, T. Samaras, M. Angelakeris and O. Kalogirou, Int. J. Hyperthermia, 2021, 38, 511–522.
19. H. Carlton and R. Ivkov, J. Appl. Phys., 2023, 133, 044302.
20. I. A. Brezovich, W. J. Atkinson and D. P. Chakraborty, Med. Phys., 1984, 11, 145–152.
21. E. Natividad, M. Castro and A. Mediano, Appl. Phys. Lett., 2011, 98, 243119.
22. N. Mille, D. De Masi, S. Faure, J. M. Asensio, B. Chaudret and J. Carrey, Appl. Phys. Lett., 2021, 119, 022407.
23. E. Natividad, M. Castro and A. Mediano, Appl. Phys. Lett., 2008, 92, 093116.
24. M. A. Vergés, R. Costo, A. Roca, J. Marco, G. Goya, C. Serna and M. d. P. Morales, J. Phys. D: Appl. Phys., 2008, 41, 134003.
25. J. Teller, F. Westphal and C. Gruettner, Magnetic nanoparticles having improved magnetic properties, US Pat., 7,691,285, 2010.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4na00383g

‡ These authors contributed equally.

---