Field- and concentration-dependent relaxation of magnetic nanoparticles and optimality conditions for magnetic fluid hyperthermia

The field-dependent relaxation dynamics of suspended magnetic nanoparticles continues to present a fascinating topic of basic science that at the same time is highly relevant for several technological and biomedical applications. Renewed interest in the intriguing behavior of magnetic nanoparticles in response to external fields has at least in parts be driven by rapid advances in magnetic fluid hyperthermia research. Although a wealth of experimental, theoretical, and simulation studies have been performed in this field in recent years, several contradictory findings have so far prevented the emergence of a consistent picture. Here, we present a dynamic mean-field theory together with comprehensive computer simulations of a microscopic model system to systematically discuss the influence of several key parameters on the relaxation dynamics, such as steric and dipolar interactions, the external magnetic field strength and frequency, as well as the ratio of Brownian and Néel relaxation time. We also discuss the specific and intrinsic loss power as measures of the efficiency of magnetic fluid heating and discuss optimality conditions in terms of fluid and field parameters. Our results are helpful to reconcile contradictory findings in the literature and provide an important step towards a more consistent understanding. In addition, our findings also help to select experimental conditions that optimize magnetic fluid heating applications.


Field-and concentrationdependent relaxation of magnetic nanoparticles and optimality conditions for magnetic fluid hyperthermia
Patrick Ilg 1,3* & Martin Kröger 2,3 The field-dependent relaxation dynamics of suspended magnetic nanoparticles continues to present a fascinating topic of basic science that at the same time is highly relevant for several technological and biomedical applications.Renewed interest in the intriguing behavior of magnetic nanoparticles in response to external fields has at least in parts be driven by rapid advances in magnetic fluid hyperthermia research.Although a wealth of experimental, theoretical, and simulation studies have been performed in this field in recent years, several contradictory findings have so far prevented the emergence of a consistent picture.Here, we present a dynamic mean-field theory together with comprehensive computer simulations of a microscopic model system to systematically discuss the influence of several key parameters on the relaxation dynamics, such as steric and dipolar interactions, the external magnetic field strength and frequency, as well as the ratio of Brownian and Néel relaxation time.We also discuss the specific and intrinsic loss power as measures of the efficiency of magnetic fluid heating and discuss optimality conditions in terms of fluid and field parameters.Our results are helpful to reconcile contradictory findings in the literature and provide an important step towards a more consistent understanding.In addition, our findings also help to select experimental conditions that optimize magnetic fluid heating applications.
The relaxation dynamics of Magnetic Nanoparticles (MNPs) is of great interest in solid state and colloidal science.This field of research has received even more attention in recent years when MNPs have found numerous technological, environmental, and biomedical applications 1,2 .For example, Magnetic Fluid Hyperthermia (MFH) is currently studied intensively as a promising method for cancer treatment by field-controlled local heating of tissue with the help of MNPs [3][4][5] .In this method, the dissipated heat results from magnetic losses of MNP relaxation following externally applied oscillating magnetic fields 3,6,7 .
One of the main challenges today is to increase the efficiency of this method by determining optimal parameter combinations of MNPs and applied fields 3,8 .Trial and error approaches are notoriously difficult due to the large parameter space 4 .Therefore, a better understanding of the physical mechanism underlying this method is highly desirable 6,9 .While earlier studies focused on magnetic losses of individual MNPs 5,7 , many recent works consider also the effect of interactions and MNP concentration on heating efficiency [10][11][12][13][14][15][16] .The latter is typically measured in terms of the specific loss power ( SLP ), also called specific absorption rate, which is the volumetric work done by the external field per unit cycle and magnetic mass density.
At present, there is a controversy in the literature about the concentration-dependence of the SLP .1][12] and references therein).However, Martinez-Boubeta et al. 13 instead reported non-monotonic behavior of SLP with concentra- tion, recently confirmed by Kim et al. 14 This would suggest there might be optimal conditions for MFH in terms of MNP concentration which maximize heating efficiency 15 .Some works suggested that the different findings might be due to intrinsic properties of MNPs and cluster shapes 10,11 , the role of Brownian rotational particle motion 17,18 , or dipolar interactions 19,20 .A recent simulation study 21 that found chain-formation to increase SLP seems to support some of these suggestions, whereas a corresponding theoretical work 22 arrives at opposite conclusions.To an extent, the confusing and sometimes contradictory findings might be related to the fact that these studies are difficult to compare, with some studies considering mobile MNPs, while other consider them immobile, and using different MNPs with different magnetic anisotropies or different coatings, adding to uncertainties in the effective dipolar interaction strengths.These differences can have pronounced effects on the heating efficiency 4,23 .Nevertheless, there seems to be agreement that MNP interactions play an important role on heating efficiency of MFH and that more research is needed to shed light on these unsolved issues [10][11][12][13]19 . It as been suggested 24,25 that the intrinsic loss power ( ILP ) is a preferable parameter compared to SLP when discussing the efficiency of magnetic fluid hyperthermia from ensembles of MNPs.We will follow this suggestion here and report results for ILP.
Besides MFH, relaxation dynamics of interacting MNPs in applied fields is also an interesting topic from a theoretical point of view.Starting with the famous Landau-Lifshitz-Gilbert equation of monodomain magnetics, the field-dependent relaxation of non-interacting MNPs suspended in viscous liquids has already been studied many years ago [26][27][28] and recently been revisited thanks to more efficient numerical methods [29][30][31][32] .Over the last years, considerable progress has been made to extend analytical results to moderately interacting MNPs within dynamical mean-field theory [33][34][35][36][37] .These theoretical results, however, rely on the rigid-dipole approximation, i.e. they neglect the internal, so-called Néel relaxation and consider Brownian particle rotation only (Fig. 1a).On the other hand, a number of simulation studies were performed for the opposite case of immobile MNPs where only Néel relaxation is present [38][39][40] , some studies include comparison to the corresponding experiments 41,42 .For immobile MNPs with strong magnetic anisotropy, interactions are found to lead to weaker heating efficiency, whereas the opposite is observed for MNPs with low or zero magnetic anisotropy 22,39,43 .
It is interesting to note that the case of interacting MNPs where Brownian particle rotation and internal Néel processes are both present is much less explored.It is only in the absence of external magnetic fields and particle interactions that Brownian and Néel relaxation can be considered as independent processes.In this case, the effective relaxation time τ eff is dominated by the faster process and given by 7 where τ B and τ N denote the Brownian and Néel relaxation time, respectively.While Eq. ( 1) is heavily used to analyze experimental results, it should be kept in mind that applied magnetic fields or dipolar interactions lead to coupling of Brownian and Néel processes, not only invalidating Eq. (1) but also giving raise to non-exponential relaxation 44 with corresponding non-Debye susceptibilities. 45,46In the presence of a static external field, it has been shown for non-interacting MNPs that the alternating current (AC) magnetic susceptibility can develop a bimodal shape where the maximum loss peak might not correspond to the longest relaxation time. 32Corresponding experimental studies 45 emphasized in particular the need for improved models to better described the field-dependent Néel contribution.
Here, we investigate the field-and interaction-dependent relaxation and AC susceptibility of MNPs by comprehensive computer simulation studies and dynamic mean-field theory (Fig. 1b).Thereby, we consider the case of mobile and magnetically hard MNPs where Néel processes can be considered as rare, thermally activated events occurring alongside rotational diffusion of the nanoparticles.We have verified in our earlier work 32 that the simplified diffusion-jump model provides a highly efficient and accurate description in this regime.From kinetic mean-field theory we obtain explicit expressions for the effective relaxation times, susceptibilities and SLP as well as ILP .These mean-field expressions predict computer simulation results of the underlying microscopic model (1) The model system is composed of mobile and magnetically hard spherical MNPs, each characterized by diameter σ and magnitude of the magnetic moment µ , while the direction of the moment either rotates with the particle (Brownian rotational relaxation time τ B ), or flips its direction inside the particle (Néel relaxation time τ N ).(b) We study an interacting ensemble of N such spherical particles at volume fraction φ and temperature T both theoretically and numerically in the presence of an oscillating magnetic field, H(ω) , whose strength is characterized by the time-independent Langevin parameter h and the magnitude h(ω) of the added small oscillating component.Our results will be discussed in terms of dimensionless parameters such as h, φ , dipolar interaction strength ∝ µ 2 /σ 3 k B T , the ratio q = τ B /τ N , and Langevin susceptibility χ L = 8 φ.
system rather well up to moderate dipolar interaction strengths.Our studies elucidate how dipolar interactions, the strength of a static bias field, and the ratio of Brownian and Néel relaxation time influence magnetization relaxation and the SLP and ILP .The results obtained here also shed light on the controversy about a possible non-monotonicity of SLP with concentration of MNPs and the underlying physical mechanisms.In particular, we identify different parameter regimes where SLP and ILP is either monotonically increasing or decreasing with concentration or attains a maximum.The latter allows us to identify optimal conditions of fluid and field parameters such that ILP and SLP are maximized.

Specific and intrinsic loss power
The central quantity used to measure magnetic fluid heating efficiency is the so-called specific loss power ( SLP ), also called specific absorption rate, that measures the volumetric power dissipation P per total mass density ρ of MNPs present in the volume element in a time-varying magnetic field 5,25 , In a cyclic process with angular frequency ω , the volumetric power dissipation is given by P = ω 2π µ 0 H • dM and is therefore proportional to the area enclosed by the hysteresis curve. 7The permeability of free space is denoted by µ 0 , the external magnetic field H and the magnetization M.
We consider time-varying fields of the form H(t) = H 0 + H(t) , where a harmonically oscillating field H(t) = H cos(ωt) is applied in addition to a static field H 0 .If the amplitude H = | H| is small enough to remain within the linear response regime, the resulting magnetization becomes M(t) = M 0 + M(t) , where M is parallel to H and can be written in terms of the in-phase ( χ ′ ) and out-of-phase ( χ ′′ ) component of the susceptibility as M(t) = H[χ ′ cos(ωt) + χ ′′ sin(ωt)] , such that SLP in Eq. (2) becomes 7 where denotes the dimensionless Langevin parameter associated with the static field H 0 with µ the magnitude of the magnetic moment of a single MNP and k B T the thermal energy.Setting h = 0 in the results below corresponds to the absence of a static field, H 0 = 0 , i.e. when only a weak oscillating field H(t) is applied.In principle, the oscillating field H can be applied in different directions with respect to the static field H 0 .Here, we will focus on the case where these fields are applied parallel to each other. 8We denote the corresponding susceptibilities as χ to distinguish them from the perpendicular case χ ⊥ .In the absence of a static bias field, H 0 = 0 , the system is isotropic and χ � = χ ⊥ = χ.
Due to the appearance of the external factors ω and H in Eq. ( 3), it has been argued that instead of SLP , the intrinsic loss power ( ILP ) should be reported to quantify the efficiency of magnetic hyperthermia from ensembles of MNPs 25 .The ILP parameter is defined as 24,25 where the last equality follows from the linear response result (3).Since the dimensions of ILP are Hm 2 /kg = m 4 /(As) 2 , we introduce the dimensionless ILP by ILP * = ρILP/(πµ 0 φ) , where ρ/φ denotes the mass density of a single MNP.From the relation (5) we obtain the dimensionless i.e. the dimensionless ILP * is equal to the imaginary part of the susceptibility normalized with the volume frac- tion of magnetic material, thereby justifying the intrinsic nature of this parameter.The density of a single MNP with Fe 2 O 3 core is ρ/φ = 5175 kg/m 3 , such that we find ILP/ILP * ≈ 7.6 × 10 −10 m 4 /(As) 2 = 0.76 nHm 2 /kg.

Previous theoretical results
The central quantity governing the loss power in the linear response regime, Eqs.(3) and ( 5), is the out-of-phase component of the AC magnetic susceptibility.For non-interacting MNPs in the absence of externally applied fields, the complex susceptibility χ * = χ ′ − iχ ′′ is given by the Debye form where χ 0 denotes the zero-frequency susceptibility and τ can be identified with the effective relaxation time τ eff defined in Eq. (1).The imaginary part χ ′′ D of the Debye susceptibility (7) entering ILP * is given by χ ′′ D = χ 0 ωτ/[1 + (ωτ ) 2 ] and achieves its maximum at ω = 1/τ where χ ′′ D ( ω) = χ 0 /2 .The corresponding (2) SLP = P ρ . ( result when inserting the Debye expression χ ′′ D into Eq.( 3) or ( 5) is routinely used to analyze magnetic heating experiments 3,7 .For non-interacting and thermally blocked MNPs, χ 0 can be identified with the Langevin sus- ceptibility χ L = nµ 0 µ 2 /(3k B T) , with n = N/V the number density of MNPs.Boltzmann's constant is denoted by k B .When internal Néel relaxation is included, χ 0 in general depends on the magnetic anisotropy and the orientation of the easy axis of the MNPs relative to the oscillating H field 47 .Some authors 45 used the empirical Havriliak-Negami model, originally proposed for dielectric spectra, to capture magnetic susceptibilities via where χ ∞ denotes the susceptibility at infinite frequency, α, β stretching exponents and χ o = χ * (0) − χ ∞ is related to the zero-frequency susceptibility.For the special case α = 0, β = 1 , χ ∞ = 0 , and χ o = χ 0 , one recov- ers the Debye model χ * D (ω) .In the general case, Eq. ( 8) describes asymmetric and broadened spectra where the maximum of the loss peak at ω is not necessarily located at 1/τ , i.e., ω depends nonlinearly on τ and the remaining Havriliak-Negami parameters.Experimental susceptibility data could be well fitted to Eq. ( 8). 45 However, due to its empirical nature, the microscopic foundation of the Havriliak-Negami model remains unclear in the present context, which makes the interpretation of the fit parameters challenging.
In the presence of a static external field, the susceptibilities and thus SLP and ILP depend on the orientation of the weak field oscillating relative to the static field.Within the classical model of non-interacting rigid dipoles, the AC susceptibilities for parallel and perpendicular orientations relative to the applied static field are obtained using the effective field approximation (EFA) 27 as with L 1 (x) = coth(x) − 1/x the Langevin function.The corresponding field-dependent relaxation times τ �,⊥ (h) are discussed in detail later.The EFA Eqs. ( 9)-( 10) hold only for non-interacting and thermally blocked MNPs.Recent works have shown that a modified mean-field approximation is able to extend these results to include dipolar interactions up to moderate strengths 35,48 .However, it remains unclear how Eqs. ( 9), ( 10) can be generalized beyond the rigid-dipole approximation, i.e. how Néel relaxation can be included alongside dipolar interactions.We found corrections to the Debye model due to Néel relaxation in terms of a diffusion-jump model 49 , but these results were restricted to zero external field.
In this context, it should be noted that the relaxation times τ are treated as fit parameters when Eqs. ( 7) and (8) are used to fit experimental data 10,45,46,50 .Below, starting from a microscopic model, we aim to derive a generalization of Eqs.(9), (10) which allows us to include the field-induced effects of Brownian as well as Néel relaxation with no adjustable parameters.

Results of dynamic mean-field theory
In the following, we study the effects of interparticle interactions and externally applied magnetic fields on the effective relaxation behavior of a system composed of N spherical particles contained in a volume V, where each particle carries a magnetic moment µu i , i ∈ {1, . . ., N} , whose orientation is encoded by the unit vector u i .Thereby, we are especially interested in the dynamics of the magnetization M which is defined as the total magnetic moment per unit volume, M = N i=1 µ�u i �/V = M sat m , where M sat = Nµ/V denotes the saturation magnetization and a mean dimensionless magnetization per particle.Here and in the following, time-dependent configurational averages are denoted by �•� .Within the many-body diffusion-jump model 49 considered here, the magnetization equation can be written as Details of the model and derivation of Eq. ( 12) are provided in the Methods section.We note that Eq. ( 12) does not represent a closed-form magnetization equation as it includes higher-order correlations due to the coupling of the local field h loc i to other particles, see Eq. (37).To make further analytical progress with Eq. ( 12), we assume the local field to be identical for all particles, h loc i = h loc .Furthermore, we assume that cross-correlations in the Néel relaxation can be neglected.Under these assumptions, Eq. ( 12) simplifies to www.nature.com/scientificreports/ We will specify to particular mean-field approximations later.For the moment, we only assume h loc = h loc ĥ with h loc = h loc (h) .Furthermore, we employ the effective field approximation (EFA) introduced by Martsenyuk et al. 51 Since EFA assumes that Brownian motion randomizes dipolar orientations sufficiently, we expect this approximation to limit the results of this section to the regime τ N ≫ τ B , and thus, q ≪ 1 , i.e. corrections to the rigid-dipole approximation.Within EFA, we find that the long-time relaxation is described by with ĥ = H 0 /|H 0 | and the field-dependent relaxation times, τ and τ ⊥ , parallel and perpendicular to the applied field, Details of the derivation can be found in the Methods section.Equations ( 15) and ( 16) show additivity of the field-dependent rates resulting from Brownian and Néel relaxation.These results involve the Langevin function L 1 defined after Eq. ( 4) and its derivative In the absence of Néel relaxation, τ N → ∞ , these expressions for the relaxation times agree with those derived earlier 35,52 .
To arrive at more explicit expressions, we need to specify h loc , i.e. choose a particular mean-field approxima- tion.In the first order modified mean-field theory (MMF1), we have 53 .Within the rigid- dipole approximation, MMF1 was found to provide accurate predictions for the relaxation dynamics in case χ L 0.5 37,54,55 .Using MMF1, explicit expressions for Eqs. ( 15) and ( 16) up to first order in χ L are obtained, see Eqs. ( 49)- (54).We verified that the Brownian contributions to these expressions agree with existing results 34 .To cover also stronger dipolar interactions, a second-order modified mean-field theory (MMF2) had been proposed with 56 Inserting this expression for h loc into Eqs.( 15) and ( 16) and expanding to second order in χ L leads to rather cumbersome expressions.For weak applied fields h, the relaxation times simplify to where τ eff = τ B /(1 + q) defined in Eq. (1) denotes the effective single-particle relaxation time in zero field.The six coefficients c �,⊥ k appearing in Eq. ( 18) are given in Eqs. ( 55)-( 57), while the h-independent coefficients are identical for parallel and perpendicular relaxation and reduce to 1 + χ L /3 + O(χ 2 L ) for q → 0 , in agreement with existing predictions. 49rom Eq. ( 18), because the c 0 's are negative, we see that initially τ decreases quadratically with increasing field strength h.We also note from Eq. ( 15) that the Néel contribution is non-monotonic in the effective field strength, i.e. the corresponding rate changes from increasing to decreasing around h loc ≈ 2 .Consequently, for fixed q and sufficiently large h loc , we find τ to be governed by the Brownian contribution and decrease corre- spondingly, while the Néel contribution freezes out.
To calculate also the dynamic magnetic susceptibility, we consider a weak, time-dependent magnetic field is externally applied in addition to the static field H 0 .Repeating the above calculations, we arrive with δm � = m � − m eq and h� = h • ĥ , at , from which we infer the two susceptibilities with amplitudes ( 14) In the present approximation, due to linearizing the relaxation dynamics ( 14) near equilibrium, the susceptibility ( 19) is of the Debye form and peaked at the relaxation times τ � , τ ⊥ , respectively.Besides the location, also the height of the peak shows a dependence on the static field strength h.For h → 0 we find isotropic behaviour with α � = α ⊥ = 1/τ eff , as expected.
The field-dependence of the Brownian contribution to the susceptibility has been calculated within EFA already some time ago 27 , at least for non-interacting MNPs.For the Néel contribution, most results in the literature are obtained for frozen particles without Brownian contribution and given (often isotropic) orientation of easy axes 27,39 .The case of mobile MNPs with Brownian and Néel contributions both present considered here has been hardly explored so far from a theoretical point of view.
With the susceptibilities at hand, we are now able to discuss the parameters SLP and ILP .As stated above, we follow Refs 24,25 .and focus first on the parameter ILP .For concreteness, we specify in Eq. ( 6) to the case where the oscillating field is parallel to the static field.Inserting the mean-field result (19) into Eq.( 6) we obtain From Eq. ( 22), we find that the frequency dependence of ILP * follows the classical Debye law (7), i.e. first increases monotonically with frequency ω of the applied field, reaching a maximum at ω = 1/τ � before decreas- ing to zero for further increasing frequencies.
For weak fields and concentrations, the amplitude (20) decreases quadratically with h and linearly with χ L , where q ≡ τ B /τ N , as before.In the same limit, the relaxation times are given by Eq. ( 18) so that we obtain from Eq. ( 22) for all q and ω and for small h and χ L , where q w serves as a useful abbreviation.From Eq. ( 24) we find that ILP * is independent of the magnetic field strength h and increases initially with the dipolar interaction strength for h ≪ 1 and χ L ≪ 1 .In the weak field regime where Eq. ( 24) holds, the dependence of ILP * on interaction strength χ L changes, i.e. develops a maxi- mum with respect to χ L at provided that q w > 0 .This condition is satisfied either (i) for any frequency if q > 3/5 ; or (ii) if q < 3/5 (i.e.τ eff /τ B > 5/8 ) and (ωτ eff ) 2 > (9 − 15q)/(23 + 15q) .Thus, we expect a maximum in ILP to occur as a function of dipolar interaction strength or concentration when the frequencies ω are large enough.Note, however, that ILP * decreases for frequencies larger than τ −1 eff , so this regime might not be very relevant for MFH applications.The other case when ILP develops a maximum with concentration occurs if τ B is as large or larger than τ N .( 22) The mean-field result (22) for the dimensionless ILP * as a function of the applied field strength h and volume fraction φ .The remaining parameters are chosen as ωτ eff = 1 , = 4 , with q = 0.01 in panel (a) and q = 1 in panel (b).The solid black line indicates φ * ≡ χ * L /8 , where χ * L is given in Eq. ( 26).(c) The meanfield prediction for the critical field strength h c versus ωτ eff and q.For h < h c , the ILP * initially increases with MNP concentration, whereas the opposite is the case for h > h c .The analytic approximation (59) reproduces the exact mean-field result for h c qualitatively, and quantitatively for ωτ eff > 1 .Contour lines equidistantly spaced at h c = 0.1 intervals up to h c = 1.7.
Figure 2a,b illustrates the mean-field result (22) for ILP * for frequency ω = 1/τ eff , dipolar interaction strength = 4 and two values for q = 0.01 and q = 1 .As expected, ILP * develops a maximum as a function of φ in cer- tain cases.We note that the predicted location of the maximum, φ * = χ * L /8 with χ * L given by Eq. ( 26), is only approximate since it was derived for weak fields and small concentration.From Fig. 2a,b we furthermore find that applying a constant bias field h leads to a decrease of ILP * .The initial quadratic decrease is worked out in the Methods section, Eq. ( 58).The dependence of ILP * on concentration is more involved, showing increasing as well as decreasing behavior.We define a critical field h c by lim χ L →0 dILP * /dχ L = 0 , i.e.where ILP * changes from increasing to decreasing trend with interaction strength χ L at small concentrations.Note that within mean- field theory, this critical field is independent of the dipolar interaction strength .Figure 2c presents a contour plot of h c , showing that it tends to decrease to zero with q and ωτ , with its maximum h c ≈ 1.7104 for q → 0 and ωτ eff → 0 .Taking into account terms up to order O(h 2 , χ L ) in ILP * , we find an explicit expression for h c given by Eq. ( 59).For the special case ωτ eff = 1 , this expression simplifies to h 2 c = 15/[2(5 + 3q)].Figure 3 shows an overview of the different qualitative behaviors of ILP * (φ) according to the mean-field prediction (22).For different values of the dipolar interaction strength and frequencies, we show the behavior of ILP * in the plane spanned by the ratio of relaxation times q and the field strength h.We distinguish between different types of increase with φ (blue colors), existence of a maximum (green), and different forms of decreas- ing with φ (yellow and red colors).We observe from Fig. 3 that the ILP * behavior shows pronounced changes with h, but varies only weakly with q.We also observe that the critical field given by Eq. ( 59) is a good predictor for the location of a local maximum in ILP * as a function of concentration, especially at high frequencies, where h c is low enough to be captured by a fourth-order approximation.For ωτ eff = 6 , for example, we find h c to be small enough such that the analytic expression (59) captures the transition from initial increase of ILP * with concentration to a decrease accurately.
We will continue the discussion of ILP * in more detail later when comparing these predictions to simulation results.But before, we first evaluate the accuracy and limits of validity of the mean-field results presented in this section via a detailed comparison to simulation results of the model system presented above.

Simulation results and comparison to theory
Dynamic properties of the system within linear response regime are encoded in the complex magnetic susceptibility χ * γ (ω) .Since dynamic properties are anisotropic in the presence of a static external field H 0 , we distin- guish parallel and perpendicular components of the susceptibility with respect to H 0 .These quantities can be computed by the formula where we introduced the auto-correlation function with the fluctuations δU γ (t) = U γ (t) − �U γ � , where U γ (t) = N −1 N j=1 u j,γ (t) denotes the instantaneous γ -component of the reduced magnetization.The dot in Eq. ( 27) denotes the time derivative.The static (or zerofrequency) susceptibility is given by χ γ (0) = 3Nχ L C γ (0).
Figure 4 shows the imaginary ( χ ′′ � ) part of the complex susceptibility χ * � (ω) = χ ′ � (ω) − iχ ′′ � (ω) parallel to the applied field.The simulation data are obtained by numerically evaluating Eq. ( 27) from the correlation functions C γ (t) for a range of frequencies ω .The correlation functions themselves are calculated by evaluating Eq. ( 28) from computer simulations in the stationary state using time averages.Up to moderate dipolar interaction strengths C γ (t) = �δU γ (t)δU γ (0)�, Figure 3. Qualitative theoretical behavior of ILP * versus Langevin parameter h and the ratio q, for a particular and for three different dimensionless frequencies ωτ eff in (a)-(c).All colors (except for circles) based on theoretical MMF2 behavior of ILP * defined in Eq. ( 6) in φ ∈ [0, 0.1] : (1) ILP * increases monotonically+convex with φ , (2) increases monotonically but neither convex nor concave, (3) increases monotonically+concave, (4) goes through a maximum, (5) decreases monotonically+concave, (6) decreases monotically but neither convex nor concave, (7) decreases monotonically convex.Region (1) is absent for the present range of parameters.The green region (4) is the one which exhibits a local maximum, for all blue regions ILP * increases with φ .Beyond a certain critical h = h c (interfacial line between green and orange), ILP * decreases with φ in any case.Black line shows our analytical approximation (59) for the critical field strength h c where the initial ILP * changes from increasing to decreasing with φ .The analytical expression works very well only for h c ≪ 1 , as it was obtained via Taylor expansion about h = 0 .Colored circles mark simulation results for q ∈ {0.01, 1}.
( ≤ 2 ), we find that the susceptibilities are well described by a Debye form.Moreover, as seen in Fig. 4a,b, the mean-field result (19) with Eqs. ( 15) and ( 20) is found to predict the simulation data quite accurately for all field strengths and over a broad range of model parameters.We emphasize that no fitting parameter is involved in the comparison between mean-field theory and simulation results.For strong dipolar interactions, the MMF2 approximation is no longer accurate (see e.g.Fig. 9 and Refs 49,54 .)and the mean-field predictions (19) are therefore less reliable.From Fig. 4c,d we indeed find not only quantitative discrepancies, but that the dynamic magnetic susceptibility can no longer be accurately described in terms of a Debye model.Instead we find that the data can be fitted to the Havriliak-Negami model, Eq. ( 8).There, the stretching exponents α and β describe the asymmetric and broadened shape of the susceptibility spectra.For = 4 and φ = 0.02 , best fits are obtained for parameter values (α, β) = (0.05, 0.77), (0.07, 0.77), (0.13, 0.78), (0.21, 0.70) for q = 0.01 and h = 0.5, 1, 2, 10 , respectively.And correspondingly for q = 1 these values are (α, β) = (0.09, 0.83), (0.14, 0.91), (0.19, 0.98), (0.26, 1.56) .Thus, the values of the stretching exponent α are found to increase with magnetic field strength h, reflecting a broad- ening of the relaxation spectrum.Interestingly, the asymmetry parameter β behaves differently for q = 0.01 and q = 1 .While β remains rather insensitive to h for q = 0.01 , a significant increase of β with h is found for q = 1 .We also note the development of a shoulder at low frequencies for strong fields that is not captured by Eq. (8).
Next, we define characteristic relaxation times τ ′′ ω = 1/ ω from the peak frequencies ω of χ ′′ � .In the parameter range investigated here, χ ′′ � shows a single peak so that we can obtain ω unambiguously from a fit to the Debye function (7) in the vicinity of the maximum.Figure 5 shows the corresponding relaxation times τ ′′ ω as a function of the dimensionless field strength h.While the field-dependent decrease of the effective relaxation times has already been studied within the rigid-dipole approximation by several authors [33][34][35][36][37] , our results extend these studies beyond the rigid-dipole approximation by including additional Néel relaxation.We observe from Fig. 5 that additional Néel contributions do not change qualitatively the decrease of the effective relaxation times with increasing field strength.In the absence of an external field ( h = 0 ), dipolar interactions increase the effective � (ω)/χ � (0) parallel to the applied field is shown as a function of reduced frequency ωτ eff .Data for several values of the constant bias field h ∈ {0.5, 1, 2, 10} are included.Other parameters are chosen as ∈ {2, 4} and q ∈ {0.01, 1} in (a)-(d).Simulation results obtained from Eq. ( 27) for selected frequencies are indicated by circles.In panels (a) and (b), solid lines show the corresponding theoretical results (19) with τ , α , and h loc calculated from Eqs. ( 15), (20), and (17).In panels (c) and (d), solid lines show fits of the simulation results to Eq. ( 8); fitting parameters are mentioned in the text part.Simulations were performed using N = 2000 ( = 2 ) and N = 10000 ( = 4 ) particles for a duration of 2 × 10 5 τ B .www.nature.com/scientificreports/relaxation times, τ ′′ ω > τ eff , as already discussed by us earlier 49 (see also Ivanov and Camp 57 for a corresponding study in the rigid-dipole limit).It is interesting to note that the field-induced decrease of τ ′′ ω /τ eff is strongest in the rigid-dipole limit q → 0 and that larger values of the ratio q = τ B /τ N lead to a more gradual decrease of the relaxation times with h.We also point out that for = 2 , the mean-field approximation (15) provides accurate predictions for τ ′′ ω except for weak fields and q = 1 where Eq. ( 15) underestimates the relaxation times.In the rigid-dipole limit, a modified Weiss mean-field approach was found to be more accurate than modified meanfield approximation for the relaxation dynamics 57 .It will be interesting for future work to explore the accuracy of such an approach for field-dependent relaxation times beyond the rigid-dipole approximation.For = 4 , the MMF2 approximation is no longer applicable and the corresponding mean-field result (15) seriously underestimates τ ′′ ω for weak fields ( h 2).The deviation of the dynamic magnetic susceptibility from the Debye form reflects the non-exponential nature of the magnetization relaxation.As we have pointed out previously 32,49 , different characteristic relaxation times can be defined in this case, where the longest relaxation time does not necessarily correspond to the maximum in χ ′′ .Therefore, Fig. 6 compares different definitions of characteristic relaxation times.In addition to τ ′′ ω defined from the peak frequency of χ ′′ � , we also show the integrated relaxation time, τ = ∞ 0 C � (t)dt , a quantity that is often discussed for non-exponential relaxation 49,58 .Furthermore, we also consider the short time relaxation time τ short obtained from the initial slope of the autocorrelation function (28), τ −1 short = −(d/dt) ln C � (t)| t=0 .Finally, we also show in Fig. 6 the long-time relaxation time τ long that governs the late stages of the relaxation.To determine τ long , we fit a double-exponential form to the autocorrelation function (28), where 0 ≤ c ≤ 1 is a weight factor.Having already determined τ short , Eq. ( 29) contains only two fit parameters due to the relation We use a Bayesian information criterion to decide between single-( c = 0 ) and double-( c = 0 ) exponential fits (see Ref. 49 for more details on this procedure).When fitting data for = 4 to the Havriliak-Negami model, there is in principle also the fit value of the effective relaxation time τ in Eq. ( 8).We find this fit value to coincide with τ ′′ ω within a factor 2-3 for these parameters.To not overburden the graph, we do not include this quantity in the comparison in Fig. 6.We observe from Fig. 6 that these different relaxation times agree with each other rather well for = 2 .Some discrepancies can be discerned for weak external fields h, where τ ≈ τ long are found to be somewhat larger than τ ′′ ω ≈ τ short .For strong fields, all four definitions become indistinguishable.For strong dipolar interactions ( = 4 ), however, significant differences are apparent.Except for very weak fields, τ is found to be the larger than the other relaxation times.In addition, for weak up to moderate field strengths, τ long and τ ′′ ω are found to be roughly similar to each other, while τ short is smaller, reflecting a faster initial relaxation.For strong fields, τ long , τ ′′ ω and τ short become identical.Thus, besides the failure of MMF2 for strong dipolar interaction strengths, linearizing the magnetization Eq. ( 14) in addition misses the important distinction between τ short and τ long .
We next turn our attention to results for the ILP relevant for hyperthermia applications.Here, we limit ourselves to the linear response regime where the intrinsic loss power is given by Eq. ( 5).Thus, in addition to the frequency dependence shown above, we now focus on the field-and concentration dependence, as well as the influence of ratio q of Brownian and Néel relaxation times.To check the accuracy and limits of validity of the mean-field prediction (22), we compare this expression to our simulation results for the susceptibilities χ ′′ � .The dimensionless ILP * , Eq. ( 6), as a function of the volume fraction φ of MNPs is shown in Fig. 7a; the dipolar interaction strength is chosen as = 2 .Results for different field strengths h and different dimensionless frequencies ωτ eff of the external magnetic field are shown.Both of these parameters are seen to have a substantial effect on SLP * .The ratio q of Brownian to Néel relaxation time, on the other hand, has only a minor impact on ILP * , as seen by a comparison of the data for q = 0.01 and q = 1 .Depending on the chosen parameters, we find ILP * increasing or decreasing with φ , in qualitative with mean-field predictions.We also find that the mean-field Stochastic simulation results at volume fraction φ = 0.02 .Comparison of different relaxation times, τ ′′ ω , τ , τ short , τ long , defined in the text normalized with τ eff defined by Eq. (1).Data are shown as a function of applied field strength h for the same parameter sets as in Fig. 4, i.e., ∈ {2, 4} and q ∈ {0.01, 1} (parameters mentioned in panels).Note the different scales on the y-axis.Simulations were performed using N = 2000 ( = 2 ) and N = 10000 ( = 4 ) particles for a duration of 2 × 10 5 τ B .
prediction (22) provides a rather accurate quantitative description of the simulation data for most parameter values.In Fig. 7b, we again show ILP * as a function of φ , but for stronger dipolar interactions, = 4 .Again, we compare stochastic simulation results to the mean-field prediction (22).We observe that the values of ILP * are significantly increased compared to the case = 2 for all concentrations, field strengths and frequencies.This increase is in agreement with the mean-field prediction (22).Although strong dipolar interactions are beyond the range of validity of the MMF2 approximation, this mean-field result still describes ILP * semi-quantitatively for small concentrations or strong enough fields as long as the frequencies are not too high, ωτ eff ≤ 1.
We want to emphasize that the mean-field predictions agree with the simulation results also on a more qualitative level, as summarized in Fig. 3, showing that ILP * increases with φ for small magnetic field strengths h, but decreases for large fields.From Fig. 7a for = 2 , ωτ eff = 0.5 and 1, we indeed find ILP * increasing with φ for h ≤ 1 and decreasing for h ≥ 2 .Similarly, from Fig. 7b for = 4 , ωτ eff = 0.5 and 1, ILP * increases with φ for h = 0.5 , goes through a maximum for h = 1 , and decreases with φ for h ≥ 2 .We note that the maxima in this regime are very shallow on the scales presented in Fig. 7.
Having presented our results on ILP , we close this section with a discussion on an alleged correlation between SLP and cluster sizes 10,11 .Based on their experimental results, a phenomenological relation between SLP and the mean chain size n of MNPs was suggested.Branquinho et al. 11 argued that strong dipolar interactions decrease SLP whereas weak interactions increase SLP , such that SLP attains a maximum for rather short chains, �n� ∼ 2 . . .6 .Note that the concentration and dipolar interaction dependence of SLP and ILP are identical since these quantities differ only by the frequency and squared amplitude of the oscillating field, see Eq. ( 5).Therefore, the proposed relation SLP( n ) should be inherited by ILP( n ).
Here, we test this hypothesis using a cluster analysis of MNP configurations obtained from our computer simulations.We adopt a commonly employed energy criterion (see e.g. 58,59) according to which two MNPs belong to the same cluster if their dipolar interaction energy satisfies � dd ij /k B T ≤ −1.5 .From analyzing many statistically independent snapshots, we calculate the average cluster size n .In agreement with earlier results, Fig. 8a finds n increasing with h, φ , and .For the chosen parameters, average cluster sizes are found to be rather small, �n� ∼ 1 . . . 2 .In Fig. 8b, we plot ILP * / parametrically versus n for different values of the magnetic field h and different dipolar interaction strengths .While both parameters, h and , have a significant influence on ILP , Fig. 8b clearly shows that their combined effect is not captured by the average cluster size n alone, since the same value of n can correspond to rather different values of ILP * .Note that the quantity ILP * / depends on n and h, but is approximately independent of in this regime.These findings demonstrate that one has to be very careful using phenomenological arguments in terms of rigid chain-like aggregates to explain the behavior of SLP or ILP in this parameter regime.22): dashed and solid lines).For both cases (a) and (b) the figure shows result at three different frequencies: ωτ eff = 0.5 (left), ωτ eff = 1 (middle), and ωτ eff = 6 (right) with τ eff from Eq. (1).

Conclusions
In this communication, we have developed a kinetic mean-field theory and performed extensive computer simulations of the dynamics of a microscopic model system that accounts for steric and dipolar interactions between MNPs, as well as their Brownian and Néel relaxation in the presence of external magnetic fields.Explicit expressions for the field-dependent relaxation times are obtained by a combination of modified mean-field theory and effective field approximation.The theoretical expressions are seen to provide a rather accurate prediction of the simulation results over a broad range of parameters, including weak up to moderate dipolar interaction strengths.Mean-field theory and simulations both show that the dynamic magnetic susceptibility and the SLP and ILP depend not only on the magnitude and frequency of the external magnetic field, but also on the con- centration and intrinsic properties of the MNPs and their environment.Our mean-field results allow us to map out the qualitative behavior of the SLP and ILP over a broad range of model parameters.In particular, we explore parameter regions where the ILP increases with MNP concentration, regions where ILP decreases, and where a maximum in ILP is expected.We have tested these theoretical predictions and overall obtained good, often quantitative, agreement with simulation results.Therefore, our findings help to maximize ILP by identifying optimal conditions not only for the external magnetic field strength and frequency, but also for MNP properties, concentration and environment, via dipolar interaction strength and ratio of Brownian to Néel relaxation time.Our theoretical results become inaccurate for strong dipolar interactions where modified mean-field theory is known to break down.In this regime, we also observe deviations from Debye-like susceptibilities, pointing at the importance of non-exponential relaxation behavior.Extensions of the theory to strong dipolar interactions which are able to capture these effects are left for future work.
The presented results can be used to support several experimental findings and interpretations.In the experiments reported by Piñeiro-Redondo et al. 17 , for example, SLP was reported to increase or decrease with concen- tration when MNPs with smaller or larger steric shells were used, respectively.Since the size of the steric shell affects the Brownian but not the Néel relaxation time, their observations can be rationalized with our model using different values of q.For the experiments reported by Kim et al. 14 , the ratio q is very small and frequencies higher than 1/τ eff are chosen.For these conditions, the critical field h c defined in Eq. ( 59) is small and since field strengths h > h c were applied by Kim et al. 14 , all their experiments in the dilute regime showed SLP decreasing with increasing concentration, in qualitative agreement with our mean-field result.In the opposite limit of slow Brownian relaxation, q ≫ 1 , we find the critical field to become very low for any frequency, h c ∼ O(q −1/2 ) .Although the approximation (EFA) we use is not quantitatively accurate in this regime, a number of experiments for quasi-immobilized MNPs indeed find SLP to decrease with increasing concentration 10,60 .Some authors have put forward a phenomenological model to explain the behavior of SLP in terms chain- like structures formed by MNPs 11 .We have tested this idea using a detailed cluster analysis applied to many MNP configurations obtained in our computer simulations.In agreement with earlier results, we confirm that stronger dipolar interactions and larger field strengths lead to an increase in the effective cluster size n of MNPs.However, we find that the corresponding increase in ILP or SLP can not be explained by the change in n alone, strongly suggesting that the phenomenological model is too simplistic for the parameter region investigated here.A more qualitative explanation for the non-monotonic concentration-dependence of SLP has been proposed 14 based on the relative importance of various energy contributions.For large enough concentrations, these authors also assume rigid, chain-like structures to dominate the response, leading to SLP decreasing with increasing concentration.For lower concentrations, incoherent modes are postulated to dominate since dipolar interactions are balanced by steric and magnetic field interactions.In this regime, it is argued that SLP increases with increasing concentration.Within our study, we have not detected evidence supporting this picture.Instead, our studies reveal that SLP and ILP need to be considered in the wider parameter space including the external field as well as the fluid parameters.
Steric and dipolar interaction effects are of great interest in the recent literature on relaxation and power absorption of suspended MNPs in general.In view of biomedical applications such as hyperthermia, magnetic www.nature.com/scientificreports/particle imaging or magnetic spectroscopy, this interest is in part driven by notorious agglomeration kinetics occurring in biological environments 61 .The vastly different time scales influencing MNP dynamics presents a major challenge for modelling and simulations.Although our model allows us to consider different degrees of MNP rotational mobility relative to internal Néel relaxation processes, further work is needed to better understand the interplay between MNP relaxation and agglomeration kinetics in biological environments.

Model system
We here employ the same model system that we studied earlier 49 , but now include an externally applied magnetic field.In order to make the paper self-contained, we here give a brief description of the model, referring the reader to the existing publication for more details and the original Ref 62 .for a justification of the model in the dilute, non-interacting limit.

Interaction potential
Our system consists of N interacting particles (MNPs) in a volume V, corresponding to a volume fraction φ = πnσ 3 /6 with effective diameter σ .We denote with r i and µ i = µu i the position and magnetic moment of particle i, respectively, with i = 1, . . ., N .The three-dimensional unit vectors u i denote the orientation of the magnetic moment.For simplicity, we consider monodisperse systems where the magnitude µ = |µ i | of the magnetic moment is identical for all particles.Generalizations to polydisperse systems are straightforward.
Interparticle interactions as well as the interaction with the external magnetic field H are described by the potential with h = µ 0 µH/k B T the dimensionless external field, and where dd ij denotes the point-dipole interaction between particles i and j, The unit vector pointing from particle j to i is defined by rij = (r i − r j )/r ij with r ij = |r i − r j | .The dipolar inter- action parameter plays an important role since it measures the strength of dipole-dipole interactions relative to thermal energy 63 .Note, that the two dimensionless parameters h and adsorb µ 0 , µ , k B T , and H, while σ sets the length unit.
The spherically symmetric contribution s ij describes steric repulsion among the particles and is commonly modelled by a purely repulsive Lennard-Jones potential, where r c = 2 1/6 σ and �(x) = 1 if x > 0 and zero otherwise, denotes the Heaviside step function.The Lennard- Jones parameter ε controls the strength of the repulsion, whereas σ is a measure for the hydrodynamic particle diameter.In the following, we fix ǫ = k B T , such that the static properties are fully determined by the dimension- less parameters φ , , h.Static properties of dipolar systems in an applied magnetic field described by Eq. ( 30) have been studied in the literature to great extent, as also reflected by existing reviews 64,65 .
It should be noted that the potential (30) does not include contributions from the magnetic anisotropy energy that arises when the magnetisation direction u i deviates from the particle's easy axis 26,46 .Here, we assume that the anisotropy energy is sufficiently large so that magnetic moment and easy axis can be considered to be wellaligned.For cobalt and iron-oxide particles, for example, this condition is fulfilled for particles with magnetic core diameters larger than 5 and 12 nm, respectively.

Dynamics
We follow common practice in colloidal science and model the dynamics of the nanoparticles as overdamped Brownian motion in a viscous carrier with translational and rotational friction coefficients ξ and ξ rot , respectively 63,66 .With the N-particle, time-dependent probability density F N (r, u; t) , where r = {r 1 , . . ., r N } , u = {u 1 , . . ., u N } , the diffusion-jump model proposed by us can be written as 49 The Fokker-Planck-Smoluchowski operator describes translational and rotational Brownian motion subject to the interaction potential as 66 , www.nature.com/scientificreports/Here, we introduced ∇ i ≡ ∂/∂r i and the rotational operator L i ≡ u i × ∂/∂u i .The model corresponding to QF N = 0 , known as "rigid-dipole approximation", neglects internal Néel relaxation processes and has been studied extensively in the literature. 35,48,54,63We note that the model (35) corresponds to the "free-draining" approximation where hydrodynamic interactions between particles are neglected.While this approximation is routinely adopted, notable exceptions including these effects have been reported 21,67 .
As indicated above, we here go beyond the rigid-dipole approximation and include Néel relaxation for the case of large magnetic anisotropy energies, where magnetisation reversals u i → −u i within individual nanoparticles can be considered as thermally activated events.Assuming statistically independent reversals and employing the detailed balance condition to ensure the Boltzmann equilibrium F eq ∼ exp [−�/k B T] is a stationary solution to the diffusion-jump dynamics (34), we find 49 with u (i) = {u 1 , . . ., u i−1 , −u i , u i+1 , . . ., u N } denoting the orientation state with the magnetic moment of particle i reversed.The hereby introduced dimensionless local field h loc i acting on particle i, combining dipolar interactions as well as the Zeeman energy, is given by Note that other (e.g.Glauber-type) rates are in principle also admissible 62 but that we found the Arrhenius rates used in Eq. ( 36) to be more appropriate to recover field-dependent Néel relaxation in limiting cases 32 .
We treat the bare, single-particle Néel relaxation time τ N as input parameter of the model.Together with the single-particle Brownian rotational diffusion time, τ B = ξ rot /(2k B T) , they form the basic time scales of our model.The microscopic time scale τ D related to the attempt frequency does not appear in the model since Néel processes are modelled in Eq. ( 36) as thermally activated events.While τ D ∼ 10 −10 s, τ B and τ N are typically many orders of magnitude larger.Therefore, the present scheme is much more efficient than directly coupling Brownian particle motion to the Landau-Lifshitz-Gilbert equation for modeling internal relaxation. 68rom Eq. ( 34), we can write the magnetization dynamics as Inserting the explicit form of the operators L and Q from Eqs. (35) and (36) into Eq.( 38) and performing partial integration, we arrive at Eq. (12).

Model parameters
With the thermal energy k B T setting the energy scale, the particle diameter σ the unit length, and the Brownian rotation diffusion τ B setting the time scale, we here summarize the dimensionless parameters that fully charac- terize the model under investigation.First, we measure the dimensionless MNP concentration in terms of the volume fraction φ = πnσ 3 /6 , with n the particle number density.The dipolar interaction strength is quantified by the dimensionless parameter defined by Eq. ( 32), while the dimensionless external magnetic field strength is given by the Langevin parameter (4).Since the model contains two basic time scales, τ B and τ N , their dimensionless ratio, is another important parameter of our model, so that τ B /τ eff = 1 + q .Below, we consider time-dependent magnetic fields oscillating with frequency ω .Therefore, the final dimensionless parameter in this study is the dimensionless frequency ωτ eff .To summarize, our model is defined by five dimensionless quantities: φ , , h, q, and ωτ eff .
Within mean-field theory, the parameters φ and do not appear separately.Instead, the role of interactions is approximately captured by their combined effect via the Langevin susceptibility χ L which can be written as

Derivation of effective relaxation times within dynamic mean-field theory
Starting point is the magnetization Eq. (13).The yet unspecified local field h loc is assumed to be of the form h loc = h loc ĥ , with h loc = y(h) and where the function y(h) depends on the particular form of the mean-field approximation and thus, in particular also on χ L .( 35) To arrive at closed-form expressions, we use the effective field approximation (EFA) that has been found to provide quite accurate approximations for dynamic properties under various circumstances 27,52,63 .The basic assumption of EFA is that nonequilibrium moments can be calculated with the equilibrium probability density, where the external field h is replaced by an effective field ζ e n, where n is a unit vector that is not necessarily parallel to h and h loc .However, in equilibrium, ζ e = h , n = ĥ and Eq. ( 41) reduces to the mean-field equilibrium distribution with y(h) = h loc the local field.
Evaluating the right hand side of Eq. ( 13) with the help of Eq. ( 41) we find where S j = �P j (u • n)� = L j (y(ζ e )) are the nonequilibrium orientational order parameters, P j (x) the jth order Legendre polynomial, and L 1 (x) = coth(x) − 1/x the Langevin function.In addition, we have defined ν = y(ζ e )n − h loc .Equation ( 42) provides a closed but nonlinear magnetization equation that can be rewritten as a time evolution equation for the effective field ζ e n.
Here, instead, we are interested in the long-time relaxation and therefore linearize Eq. ( 42) around the stationary state for constant h.To do this, we introduce the dimensionless deviation of the effective field from the applied field, = ζ e n − h , as well as its parallel, � = • ĥ , and perpendicu- lar, ⊥ = − � ĥ , component.To first order in we find ζ e = h + � � and n = ĥ + h −1 ⊥ .Thus, S j = L j (y(h ) , w h e r e S eq j = L j (h loc ) .Fu r t h e r m o r e , ν = (dh loc /dh)� � ĥ + (h loc /h)� ⊥ + O(� 2 ).
Inserting these expressions into Eq.( 42) and linearizing in , we find that the relaxation equations for the magnetization component parallel ( m � = m • ĥ ) and perpendicular ( m ⊥ = m − m � ĥ ) to the applied field separate, Finally, we express , and thus the effective field ζ e n = + h , in terms of the magnetization components, � � = (m � − S eq 1 )/[L ′ 1 (h loc )dh loc /dh] , ⊥ = (h/S eq 1 )m ⊥ , such that we arrive at the linearized magnetization relaxation Eq. ( 14) with m eq = S eq 1 = L 1 (h loc ) and effective relaxation times Using the relation S eq 2 = L 2 (h loc ) with L 2 (x) = 1 − 3L 1 (x)/x , these expressions agree with Eqs. ( 15) and (16) given above.
Note that the form of h loc has not been specified yet and therefore these expressions hold for arbitrary h loc in terms of h and χ L , in particular.

First-order modified mean-field approximation
Specifying to the first-order modified mean-field approximation (MMF1), h loc = h + χ L L 1 (h) , Eqs. ( 15) and ( 16) read to first order in χ L where is the classical EFA result for the parallel relaxation time in the non-interacting limit, using the rigid-dipole approximation 51 .The coefficient  www.nature.com/scientificreports/where we used the short-hand notation F N (t) = F N (r, u; t) .From Eq. ( 60) we find that the dynamics can be integrated forward over a short time step t by a succession of a diffusion step e t L and a jump process t Q .Since diffusion is described by the Fokker-Planck operator (35), we use Brownian Dynamics simulations to propagate the configurations {r(t), u(t)} by a time step t , which requires �t ≪ τ B .
To simulate the jump process t Q for a time interval t , we consider the probability p i that the moment u i has not reversed its orientation between t = 0 and t = t .According to Eq. ( 36), this probability is given by p i = e −r i t with rate r i = e −u i •h loc i /(2τ N ) .Thus, the probability that the reversal of moment u i takes place in the time interval t is 1 − e −r i t , a well-known property of Poisson processes.Therefore, calculating the probabilities p i (t) from the current configuration and resulting h loc i (t) , we reverse the orientations of the magnetic moments u i (t) → −u i (t) in the time step t with probability 1 − p i (t).
We note that the time step t should be small enough so that (i) r i �t < 1 , i.e. at most a single reversal per magnetic moment within t .In addition, we also require (ii) that the local field h loc i can be considered constant over this time interval.This implies that the fraction x flip of magnetic moments that reverse their orientation during t must be relatively small (we ensured that x flip does not exceed 0.05).For τ N ≫ τ B (i.e.considering small corrections to rigid-dipole approximation), the time step from stochastic simulations should in general be small enough to ensure conditions (i) and (ii) above.For τ N ≈ τ B or even τ N < τ B , condition (ii) might not be satisfied.In this case, one might re-calculate the probabilities p i whenever x flip is exceeded or perform k successive jump steps, each using a reduced time step �t/k to integrate the Q part of the dynamics.In the extreme case τ N ≪ τ B (corresponding to basically immobile particles, not considered here), such schemes become inefficient and one resorts to kinetic Monte-Carlo simulations. 39,40e consider systems of N particles in a cubic simulation box of volume V.The corresponding volume fraction is defined as φ = Nπσ 3 /(6V ) , where the hydrodynamic diameter σ was introduced in Eq. (33).To simulate bulk systems we use periodic boundary conditions in all three spatial dimensions.To deal with the long-range nature of dipolar interactions, we use the so-called reaction field method 69 .In this method, a cut-off radius r RF is introduced and all interactions for particles closer than r RF are calculated explicitly according to Eq. ( 31).The far-field contribution is approximated as a magnetic continuum, leading to an additional reaction-field for every particle of the form h RF i = (σ/r RF ) 3 j∈R i u j , where the sum extends over all particles within a sphere of radius r RF centered at particle i (and includes particle i).We have verified (see e.g.Fig. 9 below) that cut-off values r RF = 5σ and 8σ are sufficient to achieve accurate results for the simulation parameters used in our study for volume fractions φ ≥ 0.05 and φ < 0.05 , respectively.The number N of MNPs ranges between N = 1000 for ≤ 2 and N = 10000 for ≥ 4 , and is generally mentioned in the figure captions.A rectangular simulation box is used which is uniaxially extended in the direction of the applied magnetic field.A time step of �t = 0.001τ B is chosen.We ensure proper equilibration before data are collected.Within the stationary regime, data are typically collected over a time interval of 10 4 τ eff .
We have performed several tests to validate our implementation of the model system.For example, we verified that results are independent of the initial configurations and there are no finite-size effects, i.e. results do not change when using larger N than reported here.We also reproduced results on cluster sizes available from the literature and several theoretical results in the limit φ → 0.
We also verified that the equilibrium magnetization M eq = M sat �u • ĥ� eq does not depend on q.This result is expected -since equilibrium properties should not depend on relaxation times -and provides a further check for proper equilibration of our samples.Moreover, our results for M eq agree quantitatively with earlier simulation www.nature.com/scientificreports/results 59 for φ = 0.157 with = 1 and = 2 (not shown).Figure 9 shows M eq normalized with the saturation magnetization M sat as a function of the dimensionless field strength h for different concentrations and dipolar interaction strengths.As reported previously 59 , dipolar interactions increase M eq compared to the non-interacting values.For = 1 and 2, this increase is very accurately described by MMF2, where the local field h loc is given by Eq. ( 17) and L 1 (x) denotes the Langevin function.Also in agreement with previous observations 48,59 , MMF2 is very accurate up to moderate interaction strengths, but is found to underpredict the magnetization for strong dipolar interactions (here = 4 ).For strong dipolar interactions, the flexible chain model could offer a promising alternative approach. 70 https://doi.org/10.1038/s41598-023-43140-8