Dynamics of interacting magnetic nanoparticles: effective behavior from competition between Brownian and Néel relaxation

Abstract

The intriguing properties of magnetic nanoparticles have sparked a growing number of theoretical studies as well as practical applications. Here, we provide the first comprehensive study of the influence of interactions on the two main relaxation mechanisms: internal (Néel) and Brownian relaxation. While non-interacting magnetic nanoparticles show Debye behavior with an effective relaxation time, many authors use this model also for the interacting case. Since Néel relaxation is typically a thermally activated process on times scales that are many orders of magnitude larger than the underlying micromagnetic times, we use extensive computer simulations employing a Brownian dynamics/Monte-Carlo algorithm to show that dipolar interactions lead to significant deviations from the Debye behavior. We find that Néel and Brownian relaxation can be considered as independent processes for short enough times until dipolar interactions lead to a coupling of these mechanisms, making the interpretation more difficult. We provide mean-field arguments that describe these short and long-time, effective relaxation times well for weak up to moderate interaction strengths. Our findings about the coupling of Brownian and Néel process and the effective relaxation time provide an important theoretical insight that will have also important consequences for the interpretation of magnetic susceptibility measurements and magnetorelaxometry analysis.

---

1 Introduction

Magnetic nanoparticles (MNPs) possess fascinating properties that can be exploited in several engineering and biomedical applications.^1–4 With typical diameters of their magnetic core on the order of 5–30 nm, MNPs are magnetic monodomain particles and therefore show “superparamagnetic” behavior.⁵ It is well-known that two basic mechanisms govern the magnetization relaxation of MNPs: (i) internal magnetization relaxation within the MNP, so-called Néel relaxation, on time scale τ_N, and (ii) so-called Brownian relaxation on time scale τ_B by rotational diffusion of the whole MNP when the particle is suspended in a viscous liquid (see Fig. 1). The relative contribution of Brownian and Néel relaxation is crucial for optimal use of MNPs in many technical as well as biomedical applications, such as hyperthermia.^3,6 Therefore, methods for determining their relative importance are currently being explored, helping to find the optimal colloids for the given application.^7,8 For ferrofluids, where MNPs are suspended in a non-magnetic carrier fluid, Rosensweig assumed that both relaxation processes occur independently so that the corresponding rates can be added up to yield an effective relaxation given by⁵Eqn (1) is fundamental for ferrofluid research and is used in numerous textbooks and research articles (see e.g.ref. 9, 10 and references therein). However, one should bear in mind that eqn (1) was originally suggested for non-interacting MNPs in the absence of external fields.

Fig. 1 Schematic of the two basic relaxation mechanisms. Brownian relaxation: entire particle including its magnetic moment μ rotates in fluid. Néel relaxation: direction of magnetic moment rotates within particle core. In the thermally activated regime considered here, deviations from the easy axis are energetically disfavored and short-lived and thus neglected in our implementation of the Néel relaxation dynamics.

Since τ_N grows almost exponentially with the magnetic volume of the nanoparticle, while τ_B grows only linearly with its hydrodynamic volume, different ratios τ_N/τ_B can be realized for core–shell particles with the same magnetic material by different sizes of the magnetic core and different thickness of the nonmagnetic shell (Section 2). The dependence of τ_B and τ_N on the strength of an external magnetic field has recently been determined experimentally by comparing the field-dependent magnetic susceptibility in the fluid and freeze-dried state for very dilute conditions.¹¹ In another set of very recent experiments, the concentration dependence of the magnetization relaxation has been measured and the Brownian and Néel contributions have been identified.⁷ While the effective Brownian relaxation time was found to increase with increasing concentration, a weaker, opposite behavior was observed for the effective Néel relaxation time. Furthermore, the corresponding dynamic magnetic susceptibility deviates strongly from the Debye law for non-interacting MNPs and the effective relaxation time (1) is not sufficient to describe the behavior. Intriguingly, some experiments suggest that Brownian effects seem to play a role under conditions where only Néel relaxation was expected.¹² In addition, hysteresis measurements pointed out the importance of dipolar interactions for magnetic losses that are relevant e.g. in hyperthermia applications.¹³

A number of computer simulation studies have investigated the magnetization dynamics when Brownian and Néel relaxation are both present. Few of these studies have thereby incorporated the stochastic Landau–Lifshitz–Gilbert (LLG) equation for a faithful representation of the internal magnetization dynamics^14–16 as suggested by Shliomis and Stepanov.^17,18 In the physically relevant regime where Néel relaxation is a rare, thermally activated process, the LLG approach is computationally very inefficient. Therefore, kinetic Monte-Carlo schemes have been used^19,20 to simulate the magnetization response of frozen multi-core magnetic particles to oscillating fields. For ferrofluids, a Monte-Carlo scheme to equilibrate the magnetic moments alongside their translational diffusion was proposed.²¹ Recently, a diffusion-jump model was proposed to efficiently model internal Néel relaxation as a thermally activated jump process alongside the Brownian rotational diffusion of the MNPs.²² For non-interacting MNPs, good agreement with the underlying model of Shliomis and Stepanov was found over a reasonable range of model parameters.

However, since the magnetic properties of interacting MNPs are much less understood,^23–25 we here address the question about the validity and form of eqn (1) for interacting MNPs. In ref. 25, concentration effects on the effective relaxation time in magnetorelaxometry were observed experimentally, but due to the lack of theoretical models were interpreted in terms of non-interacting particles. How do τ_B and τ_N depend on concentration and dipolar interaction strength? How meaningful is the distinction between Brownian and Néel relaxation in an interacting system? To answer these questions, we build on the model proposed in ref. 22 and perform extensive Brownian dynamics simulations of the translational and rotational dynamics of interacting dipolar particles coupled to thermally activated Néel relaxation processes for a broad range of concentrations and dipolar interaction strengths.

The paper is organized as follows. Section 2 provides formulation of the diffusion-jump model for interacting MNPs. Results of extensive computer simulations of this model are presented and analysed in Section 3. In particular, Sections 3.3 and 3.4 deal with the case that only Brownian and Néel relaxation is present, respectively. The corresponding models, denoted as brownian rigid dipole (RBD) model and frozen ferrofluid with Néel flip/magnetization reversal (FFMR) model, specified in Appendix B, are special cases of the full model (coupled case), for which results are presented in Section 3.4. Conclusions are offered in Section 4.

2 Model formulation

Consider a system of N interacting particles in a volume V corresponding to the number density n = N/V. Let r_i and μ_i = μu_i denote the position and magnetic moment of particle i, respectively, where μ denotes the magnitude and u_i the three-dimensional unit vector of the orientation of its magnetic moment. For simplicity we here consider monodisperse systems. To better represent experimental systems with significant polydispersity in nanoparticle sizes, generalizing the model is straightforward (e.g. along the lines of ref. 26).

The potential energy of the system in the presence of an external magnetic field H can be expressed as

where k_B and T are Boltzmann's constant and temperature, respectively. The spherically symmetric potential Φ^s models steric repulsion. Following common practice,^27,28 we employ a purely repulsive Lennard-Jones potential,with r_ij = |r_i − r_j| the distance between particles i and j, r_cut = 2^1/6σ, the interaction strength ε and the Lennard-Jones diameter σ as a measure for the spherical diameter of the particle. The dimensionless local field acting on particle i is given bywhere h = μ₀μH/k_BT denotes the dimensionless magnetic field, λ = μ₀μ²/(4πσ³k_BT) the dimensionless strength of the dipole–dipole interaction relative to thermal energy, and μ₀ the permeability of free space. The equilibrium properties of system described by eqn (2) depend strongly on λ and the volume fraction ϕ = nπσ³/6.²⁹

In principle, eqn (2) should be supplemented by the magnetic anisotropy energy Kv_m(u_i·n_i)² due to deviations of the magnetic moment direction from the particle's easy axis orientation n_i, where K is the anisotropy constant of the magnetic material and v_m the volume of the magnetic core of the nanoparticle.^10,18,30 Here, we consider MNPs that are sufficiently large with a correspondingly large magnetic anisotropy energy Kv_m that magnetic moment and easy axis can be considered to be well-aligned. For iron oxide nanoparticles, K ≈ 10⁴ J m⁻³, this means we consider magnetic core diameters larger than 12 nm, so that κ = Kv_m/k_BT ≳ 2 at room temperature, whereas for cobalt nanoparticles, K ≈ 2 × 10⁵ J m⁻³, already magnetic core diameters greater than 5 nm lead to κ ≳ 3. We note that experiments have confirmed that magnetization reversals of individual MNPs on the order of 15–30 nm are well-described by thermal activation over a single-energy barrier.³¹ We use this classical picture in the diffusion-jump model that we describe next.

Having defined the interaction potential, we also need to specify the dynamics of the system. We assume over-damped motion in a viscous carrier fluid of viscosity η_s with translational and rotational friction coefficients ξ and ξ_rot, respectively. The single particle Brownian diffusion time is given by τ_B = ξ_rot/2k_BT with ξ_rot = πη_sσ³, where η_s denotes the solvent viscosity. Let F(r̲,u̲;t) denote the probability density to find the position and magnetization orientation at time t, with the short notation r̲ = r₁,…,r_N, u̲ = u₁,…,u_N. We propose the following model for the translational and rotational dynamics

with the Fokker–Planck operatorwith ∇_i ≡ ∂/∂r_i and the rotational operator . The first and second line describe translational and rotational diffusion subject to the interaction potential Φ and an imposed flow field with (transpose) velocity gradient κ and (one-half) vorticity Ω. Thus, the model reduces for [scr Q, script letter Q] = 0 to the “rigid dipole model” (τ_N → ∞) where Néel relaxation is fully suppressed and the magnetic moment is assumed to remain permanently fixed within the particle. The rigid dipole model is frequently used to study equilibrium and nonequilibrium dynamics³² in ferrofluids.

The term [scr Q, script letter Q] in eqn (5) describes Néel relaxation as jump processes u̲ → u̲′,

with transition rates w. Equations of the type (5) with (6) and (7) are known in the literature as “differential Chapman–Kolmogorov” equations.³³ In order to complete the model, we therefore need to specify the transition rates w in eqn (7). We here follow ref. 22 and assume sufficiently large magnetic anisotropy energies κ that the magnetic moments are well-aligned with the particle's easy axis, u_i = ±n_i, so that the transition rates w vanish unless u_i′ = −u_i. Assuming furthermore individual and statistically independent magnetization reversals, u̲⁽ⁱ⁾ ≡ u₁,…,u_i−1, −u_i, u_i+1…,u_N, we use the ansatzWe ensure that the Boltzmann equilibrium F_eq ∼ exp[−Φ/k_BT] is a stationary solution to the dynamics (5), i.e. [scr Q, script letter Q][F_eq] = 0, by the detailed balance conditionGeneralizing the model proposed in ref. 22 to include dipolar interactions we choose

Λ_i(r̲,u̲) = Λ₀ exp[−u_i·h_i,loc], (10)

which satisfies the detailed balance condition (9) identically. Through the local field h_i,loc given by eqn (4), the rates Λ_i include dipolar interactions. Choosing the constant Λ₀ = 1/(2τ_N) ensures the correct limit in the absence of external fields and dipolar interactions²² (Appendix A). Our model can also be used to study the blocking temperature in relation to the field-cooled and zero-field cooled magnetization, extending previous studies on non-interacting MNPs.³⁴ Some authors found it necessary to include distributions of magnetic anisotropy constants, representing the energy barriers of individual particles.¹² This could be incorporated in the present model by a distribution of Λ₀ values. The dependence of the single-particle Néel relaxation time τ_N on the anisotropy parameter κ was already derived by Brown and reads for κ ≳ 2, where τ₀ is typically in the order of 10⁻¹⁰–10⁻⁹ s, usually two to four orders of magnitude smaller than τ_N.¹⁰ Therefore, by eliminating the time scale τ₀, our diffusion-jump model is orders of magnitude more efficient than approaches based on the LLG equation that need to resolve the fast internal magnetization dynamics on the time scale τ₀.

To summarize, the magnetic particles are specified by the dimensionless parameters λ and τ_B/τ_N, denoting the dipolar interaction parameter and the ratio of single-particle Brown and Néel relaxation times, respectively. The many-particle system is characterized by the amount and size of such magnetic particles, described by the volume fraction ϕ, and the strength of the steric repulsion given by ε = k_BT.

3 Results and discussion

3.1 Simulation approach

To solve the coupled model described by eqn (5)–(10), we consider a time step Δt short enough so that the solution to eqn (5) can be approximated by

F(t + Δt) ≈ e^ΔtLF(t) + Δt[scr Q, script letter Q][F(t)]. (11)

For simplicity we used the short notation F(t) for the probability density F(r̲,u̲;t). The first term on the right hand side of eqn (11) represents Brownian translational and rotational motion described by the Fokker–Planck operator (6), whereas the second term describes magnetization reversals according to eqn (7)–(10). The numerical solution can therefore be broken up in two independent parts per time step: (i) Brownian dynamics simulations are employed to advance particles' positions and orientations for a time step Δt according to the translational and rotational Brownian motion (6). This step is identical to many previous simulation studies on ferrofluids that used the rigid-dipole approximation (see e.g.ref. 32 and references therein). (ii) A kinetic Monte-Carlo scheme to implement the jump processes described by [scr Q, script letter Q][F]: (a) Select a dipole i at random and calculate the effective rate Λ_i = Λ_i(r̲(t),u̲(t)). (b) Flip the orientation of the magnetic dipole, u_i → −u_i, with probability p = 1 − e^−ΔtΛ_i. (c) repeat steps (a) and (b) N times.

Parts (i) and (ii) conclude one integration step time of length Δt. More details on the algorithm are given in Appendix B. This hybrid algorithm where Brownian dynamics and kinetic Monte-Carlo steps are alternated is similar to stochastic reaction-diffusion algorithms.^33,35

The main quantity of interest in this study is the instantaneous magnetization which is defined as the total magnetic moment per unit volume,

where we have defined the saturation magnetization M_sat = Nμ/V and the unit vector parallel to the external field direction Ĥ. The reduced magnetization U is defined by U(t) = Ĥ·M(t)/M_sat.

Before discussing the validity and interpretation of eqn (1) for interacting magnetic nanoparticles, we first study the Brownian (RBD model) and Néel (FFMR model) relaxation separately in Sections 3.2 and 3.3, before addressing their combined effect (coupled model) in Section 3.4. To prevent confusion from the outset, we are going to use a unique notation (Table 1) for relaxation times occurring under different circumstances. The single particle relaxation times are the basic quantities determined by the nanoparticle and the solvent. Since dipolar interactions lead in general to non-exponential relaxation, we characterise the relaxation of the interacting system by the short-time behavior τ_short and an effective, integrated relaxation time [small tau, Greek, macron]. We also fit the relaxation to a superposition of exponential decays with τ_n the relaxation time of mode n. Experimental data can be analyzed the same way. The various relaxation times are necessary to characterize the features of the full time-dependent signal in a compact fashion.

Table 1 Different relaxation times for the Brownian RBD (subscript B), Néelian FFMR (subscript N), and coupled models

	RBD model	FFMR model	Coupled
Single particle	τ B	τ N	τ eff
Short-time	τ B	τ ^shortN	τ short
Mode n	τ B,n	τ N,n	τ n
Effective	[small tau, Greek, macron] B	[small tau, Greek, macron] N	[small tau, Greek, macron]

---

3.2 Brownian relaxation for rigid dipoles (RBD model)

Computer simulation studies of ferrofluids frequently rely on the rigid-dipole approximation (RBD model), where Néel relaxation is completely ignored (τ_N → ∞).^{29,32,36–38} In this section, we study this model of interacting magnetic nanoparticles, U_B(t), for a wide range of particle concentration and interaction strengths. Fig. 2 shows a selection of our computer simulation results for the magnetization relaxation after a strong ordering magnetic field has been switched off at time t = 0. Results for different volume fractions ϕ and dipolar interaction strengths λ are shown, corresponding to different values of the Langevin susceptibility χ_L = 8λϕ. We find from Fig. 2 that the magnetization relaxation deviates from a single-exponential behavior except for very weak interaction strengths. Therefore, we use the multi-mode ansatzto fit our simulation data, where k = 2, 4,… denotes the number of parameters, τ_B,n denotes the relaxation time of mode n and c_n the corresponding weight. Note that the relaxation spectrum {τ_B,n} can in principle be determined from the Laplace transform of the signal, where the shortest relaxation time captures the initial decay, while the longest relaxation time describes the exponential decay at long times. The short-time behavior is always given by the single-particle relaxation time, U_B(t) = 1 − t/τ_B + [scr O, script letter O](t²). Therefore, we determine c₀ and τ_B,0 from and . In addition, we use the convention τ_B,n > τ_B,n+1 for n > 0 so that τ_B,1 denotes the longest relaxation time.

Fig. 2 RBD model. Magnetization relaxation U_B(t), eqn (12) from a strongly magnetized initial state as a function of scaled time t/τ_B. The three panels correspond to Langevin susceptibilities (a) χ_L = 0.5, (b) 1 and (c) 5, respectively. In each panel, the concentration increases from top to bottom from ϕ = 0.02 to 0.4 as indicated in the legend. Measured data shown in gray, our fits using (13) with k according to Fig. 3 are colored. Fitted longest relaxation times τ_B,1 shown in Fig. 4.

Using a Bayesian information criterion (Appendix B.5 for further details), we determine the minimum number of parameters k in eqn (13) needed for a quantitative description of the magnetization relaxation data. Fig. 3 shows the result for k in the χ_L − ϕ parameter space. We observe that pronounced multi-mode relaxation is found for small concentrations and large values of χ_L, i.e. large dipolar interaction strengths λ.

Fig. 3 RBD model. The minimum number of parameters k in the magnetization relaxation eqn (13) determined by a Bayesian information criterion for different points in the χ_L − ϕ parameter plane. The shaded gray area shows the regime λ > λ_fs with λ_fs = 8, corresponding to significant finite-size effects.

We use a cluster analysis in order to investigate structure formation in this system directly. Following earlier works,^39,40 we define two particles to belong to the same cluster if their dipolar interaction energy is within 75% of the contact energy of two aligned dipoles. From the cluster analysis, we find that clusters start spanning the simulation box massively for large dipolar interaction strengths λ ≳ λ_fs with λ_fs = 8, whereas our results seem not to be significantly effected by finite-size effects for λ ≲ λ_fs.

Fitting our numerical results for the magnetization relaxation in the RBD model, Fig. 2, to the form (13), we obtain the effective Brownian relaxation times τ_B,n for a range of concentrations and interaction strengths. In Fig. 4, we show the longest effective Brownian relaxation times τ_B,1 normalized by the single-particle τ_B as a function of χ_L. For relatively weak interactions, we find that the increase of the longest effective relaxation time τ_B,1 is well described by the so-called first order modified-mean field model,^41,42

τ^MMF_B,1 = τ_B(1 + χ_L/3) (14)

The derivation of eqn (14) within mean-field theory is provided in Appendix C. Our findings are in good agreement with those of Sindt et al.,⁴² who found that eqn (14) provides a good description up to χ_L ≲ 1 for λ ≲ 4. For stronger interactions, Fig. 4 shows that the relaxation times grow significantly stronger than predicted by the modified mean-field model. This result is in qualitative agreement with earlier findings of slow relaxation processes in strongly interacting ferrofluids using the same rigid-dipole model.³⁶ We also note that the relaxation times are no longer a function of χ_L alone, but depend on the volume fraction ϕ and λ separately. This finding is in line with predictions from cluster expansions for the magnetization⁴³ and magnetoviscosity.⁴⁴ As a word of caution, we estimate that finite-size effects for λ > λ_fs affect those data with τ_B,1 ≳ 100τ_B.

Fig. 4 RBD model. The longest effective Brownian relaxation time τ_B,1 normalized by the bare relaxation time τ_B as a function of the Langevin susceptibility χ_L (on a semilogarithmic scale). Dashed line shows the prediction of the modified mean-field model (14). Panel (a) is a zoom of panel (b) for small χ_L values.

Fig. 5 shows that the deviation of the relaxation time from the modified mean field prediction is a function of λ only and can be captured in the form

where a ≈ 40, λ_c = 4.9 ± 0.5 and z = 0.4 ± 0.1. The strong increase in the Brownian relaxation time beyond the modified mean-field prediction shown in Fig. 5 and indicated by eqn (15) suggests that major structure formation occurs in dilute systems for λ ≳ 5. Our cluster analysis confirms this conclusion. In addition, this conclusion is also in line with previous computer simulations of various dipolar model systems that have established the formation of chains and rings in the dilute, strongly interacting regime.^29,45

Fig. 5 RBD model. The longest Brownian relaxation time τ_B,1 scaled by the prediction τ^MMF_B,1, eqn (14)versus dipolar interaction strength λ. Solid black line shows the approximation (15). Same data as shown in Fig. 4. Panel (a) shows the data on a linear, panel (b) in a double-logarithmic scale.

3.3 Néel relaxation for immobile particles (FFMR model)

In order to study the Néel relaxation process in isolation, we consider in this section immobile particles by suppressing Brownian motion completely (τ_B → ∞). For a better comparison, we start with the same initial conditions as used in Section 3.2, i.e. we first perform standard BD simulations in the presence of a strong ordering external field to prepare initial configurations. Once steady state is reached, at time t = 0 we instantaneously switch off the external field and run the FFMR model. The frozen particle configuration represents a quenched disorder and we later need to average our results over independently prepared equilibrium starting configurations. Such systems and the corresponding field-cooled magnetization have attracted considerable interest in the literature on solid magnetic materials and spin glasses due to slow relaxation processes.^46,47 In view of biomedical applications, more studies on magnetic properties of immobile MNPs have appeared in recent years (see e.g. Zhao et al.¹⁶ for a simulation study on magnetization relaxation of non-interacting MNPs and Jonasson et al.²⁰ where dipolar interactions within a multi-core particle are included).

Fig. 6 shows our FFMR model result for the relaxation of the reduced magnetization U_N(t) = M(t)/M_sat averaged over independent, frozen particle configurations from a strongly magnetized initial state. For weak interactions, we observe a single-exponential decay. Upon increasing the interaction strength χ_L and depending on volume fraction ϕ, we find both, the appearance of a remanent magnetization M_r = U_∞M_sat > 0, and two slightly different relaxation times, τ^short_N < [small tau, Greek, macron]_N, where τ^short_N characterizes the initial decay of the magnetization. The appearance of a remanent magnetisation M_r in field-cooled magnetic systems at low temperatures has been observed experimentally.⁴⁸ One expects M_r > 0.5M_sat, i.e. U_∞ > 0.5, for ferromagnetic interactions and M_r < 0.5M_sat for anti-ferromagnetic interactions.⁴⁶

Fig. 6 FFMR model. Magnetization relaxation U_N(t) averaged over 10 independent, frozen particle configurations from a strongly magnetized initial state as a function of reduced time t/τ_N for different ϕ. Each panel has its own χ_L value. Measured data (black) shown together with the approximant (16). Note the different range of U_N values for the largest χ_L. A zoom into the short time window is provided in the appendix, Fig. 13.

We identify the short-time relaxation time τ^short_N (Fig. 7a) from our data (exemplarily shown in Fig. 6) by the inverse initial slope of d ln U_N/dt for times smaller than 0.05τ_N. The remaining characteristics U_∞ and [small tau, Greek, macron]_N are obtained reliably with τ^short_N at hand. Because more than two different relaxation times are not required to fit all our results with a relative deviation of less than 1% to an expression of the form (13) with k = 2 for U(t) − U_∞, we use the ansatz

U_N(t) = U_∞ + (c − U_∞)e^−t/τ_N,1 + (1 − c)e^−t/τ_N,2 (16)

where τ_N,2 = τ_N,1τ^short_N(1 − c)/[τ_N,1 + τ^short_N(U_∞ − c)], U_∞ ∈ [0,1], c ∈ [0,1] and c − U_∞ ∈ [0,1]. With this unique choice for τ_N,2, eqn (16) exhibits by construction the four required features, U_N(0) = 1, U_N(t) = 1 − t/τ^short_N + O(t²), , as well as an exponential decrease at t ≫ τ_N,1 ≥ τ^short_N that is governed by relaxation time τ_N,1 in the absence of a remanent magnetization. An effective relaxation time [small tau, Greek, macron]_N we determine uniquely, also in the presence of nonvanishing U_∞, viawith U_N(t) from eqn (16) to capture the regime of times beyond simulation time. We thus have [small tau, Greek, macron]_N = (c − U_∞)τ_N,1 + (1 − c)τ_N,2 using the parameters already obtained. Both, the resulting remanent U_∞ and effective relaxation time [small tau, Greek, macron]_N, are shown versus χ_L for various ϕ in Fig. 7b and c.

Fig. 7 FFMR model. Quantities characterizing the magnetization dynamics in the absence of Brownian relaxation. The following quantities had been evaluated viaeqn (16) from the measured data: (a) inverse short relaxation time τ^short_N, (b) remanent reduced magnetization U_∞, and (c) inverse effective Néel relaxation time [small tau, Greek, macron]_N normalized by the bare relaxation time τ_N, all as a function of the Langevin susceptibility χ_L. Solid line in (a) and dashed line in (c) highlight predictions of the mean-field estimates (18) (b_a = 1) and (19) (s = 1), respectively.

For weakly interacting systems (χ_L ≲ 0.5), we expect U_∞ = 0 when starting with sufficiently random positions,³⁰ whereas U_∞ ≈ 0.5 was found in the low-temperature regime in an earlier study where particles occupy randomly the sites of a regular lattice for volume fractions ϕ < 0.5 and λ = 1.⁴⁸ These authors found that dipolar interactions can lead to U_∞ > 0 also in the high-temperature regime. Our FFMR model results (Fig. 7b) confirm these expectations, while U_∞ quickly exceeds 0.5 with increasing χ_L. Beyond χ_L ≈ 3 all systems up to the largest volume fractions exhibit a remanent magnetization.

From Fig. 7a and c we find that the short time τ^short_N drops almost mono-exponentially with respect to χ_L, while the effective Néel relaxation time [small tau, Greek, macron]_N decreases only approximately linearly with χ_L. The relaxation times approach unique curves with increasing concentration, i.e. decreasing λ. To understand these behaviors, we here offer simple mean-field arguments for (i) the short-time relaxation τ^short_N and (ii) the effective relaxation time [small tau, Greek, macron]_N. From eqn (5)–(7), we find dU_N/dt = −τ_N⁻¹〈e^−u·h_locu_z〉_t where U_N(t) = 〈u_z〉_t and the local field h_loc defined by eqn (4). Since we are dealing with an interacting many-particle system, we do not expect to find a closed-form solution to this equation, and thus treat problem (i) and (ii) using a mean-field approach.

(i) For the short-time relaxation from a strongly ordered initial state, we can use a factorization approximation 〈e^−u·h_locu_z〉_t ≈ 〈e^−u·h_loc〉_tU_N(t). Together with a first-order cumulant expansion, we find dU_N/dt ≈ −U_N/τ^short_N with the short-time effective relaxation time τ^short_N = τ_Ne^{〈u·h_loc〉}. To make further progress, we need to approximate the local field. For weak to moderate dipolar interactions, we adopt the modified mean-field model, h_loc ≈ h + χ_L〈u〉.⁴¹ For strong interactions, we follow the chain-formation model⁴⁹ and include in the local field only the contributions from a single, perfectly straight chain. Using in addition that U_N ≈ 1 for short times, we arrive at the estimate

τ_N/τ^short_N ≈ e^−b_aχ_L (18)

where b_a = 1 for weak and moderate dipolar interactions and b_a = ζ/(4ϕ) in the chain-forming regime, where ζ defined in eqn (44) depends weakly on the mean chain length (more details on the derivation can be found in the Appendix D). We observe that eqn (18) describes the FFMR data in Fig. 7a accurately for weak interactions and also explains qualitatively the behavior for stronger interactions.

(ii) Estimating the effective relaxation time [small tau, Greek, macron]_N is more challenging. For weak dipolar interactions, we again use the modified mean-field approximation for h_loc. We solve the differential equation for U_N to first order in χ_L to arrive at the effective relaxation time

τ_N/[small tau, Greek, macron]_N ≈ 1 − sχ_L (19)

with s = 5/9. The corresponding calculation in the chain-formation regime leads to a more complicated expression, but might be approximated by (19) with s = 5ζ/(36ϕ). Details of the calculation are given in Appendix D. We emphasize that these theoretical estimates assume U_∞ = 0 and therefore apply only for moderate to strong but not too strong interactions. Comparing the prediction (19) to Fig. 7c, we find that the mean-field predictions are qualitatively correct. However, they are quantitatively less accurate than for the short-time regime.

Finally a word of caution on interpreting [small tau, Greek, macron]_N, eqn (17), seems in order. While [small tau, Greek, macron]_N is well-defined also when U_∞ > 0, its interpretation as an effective relaxation time becomes more and more questionable as U_∞ increases. Therefore, we have omitted in Fig. 7c those data for which U_∞ > 0.

3.4 Effective relaxation in coupled system

Having studied the effect of inter-particle interactions on Brownian and Néel relaxation via the RBD and FFMR models separately in Sections 3.2 and 3.3, we now investigate the fully coupled system where Brownian translational and rotational motion occurs simultaneously with Néel processes. As above, we equilibrate a number of statistically independent configurations in the presence of a strong external magnetic field. At t = 0, the field is turned off instantaneously and we follow the relaxation dynamics as described by the differential Chapman–Kolmogorov eqn (5). Snapshots of some particle configuration during relaxation are shown in Fig. 8.

Fig. 8 Coupled model. Snapshots of particle configurations at different times during relaxation for N = 4000 particles, χ_L = 0.9, ϕ = 0.1, τ_N/τ_B = 0.1. The time t and corresponding reduced magnetization U(t) are indicated in each panel. The colour scheme encodes the projection of the magnetic moments on the direction of the orienting field (now switched off).

Fig. 9 Coupled model. The minimum number of parameters k in the magnetization relaxation eqn (21) determined by a Bayesian information criterion. Different points in the χ_L–ϕ–(τ_B/τ_N) parameter space are shown, where τ_B/τ_N – four values for each (χ_L,ϕ)-pair – increases from front to back.

The resulting relaxation behavior of the reduced magnetization U(t) as a function of time for selected parameter values is shown in Fig. 10. It shows some similarities with the corresponding Fig. 7 for the FFMR model. We now carefully analyse the data and perform a thorough comparison. From the inverse initial slope of ln U(t) versus time t, we determine the effective short-time relaxation time τ_short. The values obtained for τ_short are shown in Fig. 11 as τ_eff/τ_short versus χ_L for different volume fractions ϕ and different ratios τ_N/τ_B. For vanishing interactions, χ_L → 0, τ_short becomes identical to the single-particle effective relaxation time τ_eff defined in eqn (1). Fig. 11 shows that dipolar interactions lead to a significant increase of τ_short compared to τ_eff. We find that the mean-field prediction

describes this increase quite accurately for ϕ ≳ 0.1. Details on the derivation of eqn (20) are given in Appendix E. For more dilute systems, the mean-field prediction is accurate only for weak enough dipolar interactions but underpredicts τ_short for larger values of λ. It is interesting to note that the Néel contribution to 1/τ_short vanishes for large χ_L so that 1/τ_short → 1/τ_B. These observations are very similar to the FFMR case when only Néel relaxation is present. Since we know from the RBD model that dipolar interactions do not affect the short-time relaxation when Néel processes are absent, we conclude that Brownian and Néel relaxation can be considered as independent processes at short times. Therefore, the effective short-time relaxation rate 1/τ_short is to a very good approximation the sum of the Brownian and Néel short-time relaxation rates.

Fig. 10 Coupled model. The magnetization relaxation U(t) as a function of reduced time t/τ_eff for τ_B/τ_N = 1 (each panel has a different χ_L value). Measured data (black) together with our fits to eqn (21) using a Bayesian information criterion, see Appendix (B.5) and Fig. 9. Different panels correspond to different values of χ_L, with different volume fractions ϕ color-coded as indicated in the legend.

Fig. 11 Coupled model. The scaled, (a) inverse short-time relaxation time τ_eff/τ_short and (b) inverse effective relaxation time τ_eff/[small tau, Greek, macron] are shown versus the Langevin susceptibility χ_L for various volume fractions ϕ. In each panel τ_N increases from bottom to top: yellow, τ_N/τ_B = 0.1; green, τ_N/τ_B = 0.25; blue, τ_N/τ_B = 1.0; purple, τ_N/τ_B = 2.5. Dashed lines indicate the mean-field prediction (20) for τ_short and (23) for [small tau, Greek, macron].

Analysing the relaxation curves in Fig. 10 in more detail, we find that they can be described very accurately by a superposition of decaying exponentials,

in the same form as in the RBD case, eqn (13). Indeed, the data (in black) are indistinguishable on the scale of the plot from the fit (in color) in Fig. 10. Note that for weak interactions, χ_L = 0.1, results for different volume fractions ϕ superimpose on the same curve. For larger χ_L, however, this is no longer the case and smaller ϕ corresponding to larger λ lead to slower relaxation. We determine the minimum number of parameters k needed in the fits shown in Fig. 10 from a Bayesian information criterion (Appendix B.5). In Fig. 9, we show the results for the minimum number of parameters k needed in the fits for the coupled system. For sufficiently weak interactions, we find k = 0, i.e. a single-exponential relaxation. In this case, the relaxation time has already been identified above from the short-time behavior. We find that additional modes of relaxation are needed for stronger interactions in a very similar manner as for the RBD model (Fig. 3). Interestingly, for most points in the χ_L–ϕ plane, different values of τ_B/τ_N do not change the number of modes needed for an accurate fit.

Since the interpretation of the longest relaxation time in a multi-mode fit can be unclear when the corresponding weight is small, we define an effective relaxation time for the coupled system analogue to eqn (17) as

Note that Brownian relaxation enforces U_∞ = 0. Definition (22) of [small tau, Greek, macron] we find to be very useful and robust. Fig. 11b shows the inverse of [small tau, Greek, macron] scaled with the single-particle effective relaxation time τ_eff, eqn (1). Since our fits to U(t) are very accurate, we evaluate [small tau, Greek, macron] from eqn (22) as using the fit values of τ_n and c_n. We observe from Fig. 11 a that dipolar interactions lead to a strong increase in [small tau, Greek, macron] that, contrary to the short-time relaxation, does not level off at large χ_L. We find from Fig. 11b that the mean-field prediction for [small tau, Greek, macron],is not as accurate as the one for τ_short, but provides a good description for weak and moderate interaction strengths. Eqn (23) holds for χ_L < 3 and (χ_L − 1)τ_B < τ_N and we have defined τ^MF = τ_eff/(1 − χ_L/3) and b₁ = χ_L(τ_N − 2τ_B)/[(3 − χ_L)(τ_N + τ_B)]. Details on the derivation of eqn (23) are given in Appendix E. For stronger interactions, mean-field theory predicts a vanishing rate 1/[small tau, Greek, macron] at some critical value of χ_L, whereas the simulations show a gradual decrease, deviating from the mean-field predictions earlier the larger λ.

The effective relaxation time [small tau, Greek, macron] increases with χ_L as shown in Fig. 12a on a linear scale. Interestingly, a stronger increase of [small tau, Greek, macron] relative to the non-interacting τ_eff is seen when compared to the relative increase of τ_B,1 relative to τ_B obtained using the RBD model (Fig. 4). For small enough χ_L, this increase of [small tau, Greek, macron] is well described by the mean-field result (23). For smaller ϕ corresponding to larger dipolar interaction strength λ, the data start to diverge earlier from the mean-field result. In order to better analyse deviations from the mean-field result, Fig. 12b shows the effective relaxation time [small tau, Greek, macron], eqn (22), scaled with the mean-field prediction (61)versus the dipolar interaction parameter λ. As mentioned above, the mean-field result (23) predicts a divergence of [small tau, Greek, macron] at sufficiently strong interactions, but (61) provides a good approximation for χ_L < 1. Therefore, we consider in Fig. 12b the mean-field result (61) which does not diverge but is accurate up to second order in χ_L only. From Fig. 12b, we observe an S-shaped curve, where the mean-field result holds quite well for small λ, but seriously underestimates the simulation results for larger λ. We also note that data collapse is encouraging, even though not as good as in the case where only Brownian relaxation is present, Fig. 5. The onset of the strong increase of [small tau, Greek, macron]/[small tau, Greek, macron]_MF occurs around λ ≈ 5, very similar to our RBD model observations. Therefore, an improved theory for Brownian relaxation in strongly interacting systems is very likely to lead to correspondingly improved predictions also in the coupled case.

Fig. 12 Coupled model. (a) Reduced effective relaxation time [small tau, Greek, macron](22)versus χ_L for various ϕ and τ_N/τ_B = 1/4, up to χ_L = 1. Results are quite insensitive to the precise value of τ_N. Both the mean field result (23) (thick solid black) and its expansion (61) (dashed black) are shown for comparison. The mean field expression is not able to capture data for χ_L > 1 well. (b) The effective relaxation time [small tau, Greek, macron] scaled with the second-order mean-field prediction (61)versus the dimensionless dipolar interaction strength λ. Data for all available χ_L, ϕ and τ_N/τ_B values are shown. Black dashed line indicates the mean-field prediction (61).

3.5 Magnetic susceptibility for coupled system

Measurements of the dynamic magnetic susceptibility [small chi, Greek, tilde](ω) are more common than magnetization relaxation experiments U(t). From a theoretical point of view, these quantities are related by [small chi, Greek, tilde](ω) = χ₀[1 − iωC̃(iω)], where χ₀ denotes the static susceptibility and the Laplace-transform of the magnetization fluctuation C(t) = 〈u(t)·u(0)〉. Thanks to the fluctuation–dissipation theorem, we can identify C(t) with the relaxation function U(t). We note that by definition C̃(0) = [small tau, Greek, macron] and therefore the low-frequency behavior is given by [small chi, Greek, tilde](ω) = χ₀[1 − iω[small tau, Greek, macron] + [scr O, script letter O](ω²)], i.e. governed by the effective relaxation time defined in eqn (22).

For the coupled system, we have seen in Section 3.4 that the relaxation is very well described by a weighted sum of exponentials, eqn (21). Using this form, we evaluate the real and imaginary part of the dynamic susceptibility, [small chi, Greek, tilde] = χ′ − iχ′′ as

From Fig. 9, we find only a small parameter region with k = 0, corresponding to the classical Debye model. In the rest of the parameter space, however, dipolar interactions lead to a multi-mode relaxation and corresponding deviations from the simple Debye behavior. The multi-mode nature of the magnetic susceptibility prohibits the unique definition of a single characteristic relaxation time from χ′′. Similar difficulties are encountered in interpreting the non-exponential magnetization relaxation. While defining an effective relaxation time from relaxation measurements viaeqn (22) is common, many authors identify relaxation times with peaks in χ′′. While this approach has been used to identify Brownian and Néel relaxation times when those processes happen on well separated time scales,^7,8 our results show that such a separation runs into problems when dipolar interactions are present and time scales are not well separated, since the individual τ_n already contain Brownian and Néel contributions.

4 Conclusions and outlook

With the current work, we provide a comprehensive study of the effective magnetization dynamics of interacting magnetic nanoparticles in solution when Brownian and Néel processes are both present. We perform extensive computer simulation studies using Brownian Dynamics coupled with Monte-Carlo methods of magnetization reversals to simulate the magnetization relaxation after a strong ordering field is switched off.

First, we study Brownian (RBD model, Section 3.2) and Néel (FFMR model, Section 3.3) contributions separately, before considering the coupled case when both processes are active (Section 3.4). In all three cases, we find that dipolar interactions lead to deviations from a single-exponential decay. Therefore, a careful analysis of the data is required. For the RBD model we find that the short-time relaxation is unaffected by dipolar interactions, whereas in the other cases, the effective short-time relaxation time increases. We propose a mean-field argument that captures this increase rather well. The long-time or integrated effective relaxation time is found to increase due to dipolar interactions in all our simulations. Using the RBD model we confirm earlier findings that the increase in the long-time relaxation time is well described by the modified mean-field theory for weak to moderate interaction strengths.⁴² We also show that corrections to this result for stronger interactions should be sought in terms of the dipolar interaction strength λ, rather than χ_L. For Néel relaxation in the absence of Brownian motion (FFMR model), we propose a mean-field argument to explain the increase of the effective relaxation time. We note that an earlier study found that dipolar interactions reduce rather than increase the Néel relaxation time.⁵⁰ However, in the latter work, the authors consider CoPt core–shell nanoparticles that are compressed into dense conglomerates. Therefore, it would be interesting to extend our studies to higher volume fractions to find out conditions under which the effective Néel relaxation would be enhanced due to dipolar interactions.

When Brownian and Néel processes are both present (coupled model), we find that the corresponding rates from the RBD and FFMR model are in general not additive for interacting systems. Accordingly, equations of the form (1) hold only for very weak dipolar interactions. A noteworthy exception is the short-time relaxation, ln U(t) = −t/τ_short + [scr O, script letter O](t²) for t → 0, where the coupled rate 1/τ_short is well described as the sum of the short time RBD and FFMR rates. Therefore, in interacting systems, Brownian and Néel relaxation can be considered as independent processes for very short times only. Since our findings apply also to magnetic susceptibility spectra, their interpretation is challenging when the single particle Brownian and Néel relaxation times are comparable.

We note that our study explored the regime 0.1 ≤ τ_N/τ_B ≤ 10. For larger values of τ_N/τ_B, the rigid-dipole approximation becomes more and more accurate. In the opposite limit of dominant Néel relaxation, τ_N ≪ τ_B, pronounced two-step relaxation was observed experimentally in ref. 7. This regime is left for future research.

Finally, we want to emphasize that the coupled model studied here of combining Brownian translational and rotational dynamics with Néel-type magnetization reversals can be simulated very efficiently. The efficiency relies on modeling the Néel process as thermally activated magnetization reversals via Monte-Carlo methods. This physically appealing picture underlies also Brown's classical treatment (see e.g.ref. 10 and references therein). It breaks down, however, for particles with sufficiently small magnetic core, where the energy barriers for deviations of the particle's magnetic moment from their easy axis (anisotropy energy) become equal to or even smaller than the thermal energy. In the latter case, the internal magnetization dynamics should rather be simulated with the stochastic LLG equation. While the stochastic approach is considered to provide a more faithful representation of the internal magnetization dynamics, this model becomes very inefficient for ferrofluids containing magnetic nanoparticles with large anisotropy energies. Therefore, we propose our model to susbstitute the LLG approach in the physically relevant regime of large anisotropy barriers where magnetization reversals can be modelled as thermally activated events. It is well-suited for studying further, e.g. field-dependent properties of interacting ferrofluids.

Conflicts of interest

There are no conflicts to declare.

A Néel contribution to magnetization relaxation

In terms of the probability density F(r̲,u̲;t), the time-dependent magnetization is defined by eqn (12). Using eqn (5), the magnetization dynamics can be expressed aswhere and the first and second term describe the contribution due to Brownian dynamics and Néel relaxation, respectively. For the present modeling of Néel relaxation by eqn (7) and (8), we find for the contribution of dipole iwhere we made use of a change of integration variables from u̲ to u̲^(j) and u̲⁽ⁱ⁾ for j ≠ i and j = i, respectively. In the absence of external fields and dipolar interactions (of particular relevance for the FFMR model), we find from eqn (10) that Λ_i reduces to Λ₀ = 1/(2τ_N) and eqn (27) leads to a contribution −M/τ_N to the magnetization relaxation.

B Implementation of RBD, FFMR, and coupled model

All three parts of our algorithm to be described below (orientation-BD, translation-BD, flip-MC) can be independently switched on or off, which allows us to prepare the conditions required for the various ‘experiments’, and to define three models. The frozen ferrofluid with Néel flip/magnetization reversal (FFMR model) uses the flip-MC part only. The Brownian rigid dipole (RBD model) uses orientation and translation BD only. The full BD + MC algorithm including all parts is denoted as coupled model.

For a sphere of diameter σ, the Brownian translational and rotational friction coefficients are respectively given by ξ = 3πη_sσ and ξ_rot = πη_sσ³ with solvent viscosity η_s. We use the Brownian rotational relaxation time τ_B = ξ_rot/2k_BT to introduce reduced units. For a system with given ϕ and r_RF and periodic box size ≥2r_RF = 16σ for our choice of r_RF, the number of particles must exceed 48ϕr_RF³/πσ³ ≈ 7823ϕ. In this section we omit all units σ, ε and τ_B, as they do not appear in the algorithm.

B.1 Orientation-BD

The orientational Brownian dynamics corresponding to u̇_i = (ω_i + Ẇ) × u_i with and Wiener process W reads, using the abbreviation a_i = h^loc_i,y cos ϕ_i − h^loc_i,x sin ϕ_i,with a time step Δt, and reflecting boundary conditions for z_i ∈ [−1,1]. The ΔW's are independent random numbers with 〈ΔW²〉 = Δt. The constraints |u_i| = 1 are taken care of automatically, as we switched to spherical coordinates to represent u_i by an angle ϕ_i and its z_i-component. We use equally distributed random numbers ζ ∈ [0,1] and to generate them. Calculations are done using the reaction-field approximation with metallic boundary conditions, which amounts to replacing h^loc_i by in eqn (28). Several test of the above algorithm were performed, checking e.g. the correct stationary magnetization and rotational diffusion in the non-interacting case.

B.2 Translation-BD

The translational Brownian dynamics corresponding to is governed bywithand r_ik = r_i − r_k, r_ik = |r_ik|, r̂ = r/r, where the components of ΔW are all independent random numbers with 〈ΔW²〉 = Δt, and generated as in Section B.1. The sum in (29) does not run over all k ≠ i. For the first term λG_ik it runs over all k with r_ik ≤ r_RF. The second LJ part involves only terms with r_ik ≤ r_cut = 2^1/6. We verified the correct translational diffusion in the non-interacting case.

B.3 Flip-MC

As has been noted in ref. 22, the ansatz (10) can be generalized to

Λ_i(r̲,u̲) = Λ₀r(|x_i|)exp[−x_i] (31)

with x_i = u_i·h_i,loc which satisfies the detailed balance condition (9) identically for all choices of functions r with r > 0 and r(0) = 1. Note that the choice r = 1 lead to Arrhenius-like expression (10), whereas r(x) = sech(x) leads to Glauber-like rates of the form Λ_i = (2τ_N)⁻¹[1 + tanh(−x)]. Also the phenomenological expression r(x) = cosh(ax) has been considered²² with parameter 0 ≤ a ≤ 1.

We use a Monte-Carlo scheme to implement this jump process. At each time step Δt, every magnetic moment is flipped with its individual probability p_i = 1 − exp(−ΔtΛ_i) with Λ_i = Λ₀r(|x_i|)exp(−x_i) and x_i = u_i·h^loc_i, where h^loc_i is the same local field appearing in eqn (28), and thus already available at no additional computation cost. If not otherwise stated, r(x) = 1 is used. For the non-interacting case, we verified the equilibrium fluctuations satisfy 〈u₁(t)·u₁(0)〉 = e^−t/τ_N.

B.4 Orientation relaxation ‘experiments’

To measure the effective orientational relaxation time, starting from a state of fully saturated magnetization M = μe_z that is realized in the limit h → ∞, we prepared equilibrated systems with initial ∀_iz_i = 1 in the presence of translation-BD only. At time t = 0 we set h = 0 and either continue running using the FFMR, the RBD, or the coupled model. In each case we measure the decorrelation function 〈〈u_i·e_z〉〉(t) = M(t)·e_z/M_sat ≡ U(t), where the average is taken not only over all particles i, but also an ensemble of equilibrated initial configurations. The saturation magnetization is defined by M_sat = nμ with n = N/V the number density.

The limiting case of noninteracting particles, or χ_L → 0, we study using many particles upon switching off the translation-BD and setting λ = 0. Under these conditions each of the three cases (i)–(iii) is characterized by a monoexponential decay (U_∞ = 0), with the following relaxation times (i) [small tau, Greek, macron]_N(Λ₀) = τ_N, (ii) [small tau, Greek, macron]_B = τ_B, and (iii) we confirm that U(t) = exp(−t/τ_eff) with τ_eff defined in eqn (1). The result can be considered trivial as the rotation-BD and flip-MC are completely independent in the absence of interactions, and because there is no reason for a trapped orientation.

For the remaining results to presented we investigated 300 different systems with all combinations of χ_L ∈ {0.1, 0.2,…, 0.9, 1, 1.5, 2, 3, 4, 5}, ϕ ∈ {0.02, 0.05, 0.1, 0.2, 0.4}, and Λ₀ ∈ {0.2, 1, 2, 5} or τ_N ∈ {2.5, 1, 0.25, 0.1}, implying λ's to reside within the range λ ∈ [0.03, 31]. For each set of parameters we run all three models, and evaluated U_N(t), U_B(t), U(t), as well as the corresponding characteristic times [small tau, Greek, macron]_N, τ_B,1, τ_B,2, τ₁, τ₂, and coefficients c_B, c, and U_∞.

B.5 Multi-exponential fitting of the relaxation function

In this section, we deal with the problem of how to characterize the region in the three-dimensional (ϕ − χ_L − τ_N/τ_B) parameter space where relaxation is monoexponential to a good approximation. For each system we measured N data points (t_i,U_i) with uncertainties σ_i and we have competing fits U_fit(t) with several exponentials. For the FFMR model, cf.eqn (16), containing k = 4 (or less) parameters: τ_N,1, τ_N,2, U_∞, c. For the RBD and coupled model, cf.eqn (13) and (21) with parameters {τ_B,n,c_n} and {τ_n,c_n}, respectively. For each of those fit functions, we define the corresponding quantitieswhere U_∞ ≠ 0 occurs only in case of the FFMR model. Since uncertainties σ_i in our data do not vary significantly, they do not appear in eqn (32).

We generally expect fits with larger number of parameters k to be more accurate than with lower k. But to decide whether the double- or multi-exponential fit is preferable or rather overfitting the data, we employ the Bayesian information criterion (BIC) defined by

BIC_k = Nχ_k² + k ln N (33)

We conclude that the k-parametric fit is appropriate if BIC_k is the smallest of all BIC's.

We find that for the RBD model the short-time behavior is always given by the single particle τ_B. Therefore, dipolar interactions lead to a multi-exponential decay with two (k = 2), three (k = 4) or four (k = 6) modes, see Fig. 3. For the FFMR model, the decay is monoexponential for k = 1, monoexponential with offset for k = 2, double-exponential without offset for k = 3, and k = 4 is found to be irrelevant by its BIC value. In the coupled case, we find a single-exponential regime (k = 0) as well as superposition of two (k = 2) and three (k = 4) exponentials, see Fig. 9.

C Mean-field arguments for RBD relaxation

C.1 Short-time RBD dynamics

Consider the magnetization dynamics (26) in the absence of Néel relaxation,Using the modified mean-field approximation,⁴¹ we approximate the local field by h_loc ≈ h + χ_L〈u〉. Let n denote the orientation of the ordering field. We assume that the magnetization remains parallel to n during the relaxation where the field has been switched off, 〈u〉 = U_Bn. Assuming furthermore uniaxial symmetry, 〈uu〉 = S₂nn + (1 − S₂)I/3, where I denotes the three-dimensional unit matrix, we can rewrite eqn (34) aswhere S₂ = 〈P₂(u·n)〉 and P₂(x) = (3x² − 1)/2 denotes the second Legendre polynomial. Starting from a strongly oriented initial state, U_B ≈ S₂ ≈ 1, we find from eqn (35) that the short-time relaxation is given by the bare single-particle Brownian relaxation time, τ_B, i.e. remains unaffected by interparticle interactions.

C.2 Effective RBD relaxation time

For longer times, eqn (35) predicts a non-exponential decay due to the time-dependence of S₂. Using a simple closure approximation S₂ ≈ U_B², we arrive at the ordinary differential equation τ_bU̇_B + U_B = −bU_B³, whereand b = χ_L/[3 − χ_L]. The solution of this differential equation to the initial condition U_B(0) = 1 readsFrom eqn (37), we recover the short-time behavior U_B = 1 − t/τ_B + [scr O, script letter O](t²) for t ≪ τ_b and an exponential decay for long times governed by τ_b. Note that eqn (37) agrees to first order in χ_L with the mean-field result for the RBD model, τ^MMF_B,1 = τ_B(1 + χ_L/3), eqn (14). Our mean-field arguments therefore confirm the result (14) derived earlier. It should be noted that no expansion in χ_L was used in the derivation of eqn (37) and therefore the result is expected to hold in the whole regime where the modified mean-field approximation and the closure relation S₂ ≈ U_B² apply.

Defining an effective relaxation time , we find from (37) for small χ_L

Since the linear term is quite close to the mean-field result for the long-time relaxation, τ^MMF_B,1, the longest and the effective relaxation times agree quite well for small χ_L.

D Mean-field arguments for FFMR relaxation

Let n denote the direction of the external field with which the system was prepared and U_N(t) = 〈u·n〉_t the time-dependent scaled magnetization. Point of departure is the magnetization dynamics (26) and (27), in the absence of Brownian motion,

D.1 Short-time FFMR dynamics

Fig. 13 shows the initial magnetization relaxation U_N(t) when Brownian relaxation is absent. For the short-time dynamics, the system remains in the vicinity of the highly oriented initial state. Therefore, we use a factorization approximationInserting eqn (40) into (39) and using a cumulant expansion, we arrive atwith τ_N/τ^short_N = 〈e^−u·h_loc〉_t = e^{−κ₁+κ₂/2+…}, where the first cumulants are defined by κ₁ = 〈u·h_loc〉_t, κ₂ = 〈(u·h_loc)²〉_t − 〈u·h_loc〉_t².

Fig. 13 FFMR model. Initial magnetization relaxation U_N(t) averaged over 10 independent, frozen particle configurations from a strongly magnetized initial state as a function of reduced time t/τ_N for different χ_L and ϕ. Measured data (black) shown together with the approximant (16). Note the different range of U_N values for the largest χ_L.

In order to evaluate the cumulants, we need to specify the local field. In the regime of weak up to moderate dipolar interactions, we use the modified mean-field model⁴¹ where h_loc ≈ h + χ_L〈u〉. After the field is switched off, this expression becomes h_loc = χ_LU_Nn. With this, the cumulants can be evaluated, κ₁ = χ_LU_N² and κ₂ = χ_L²U_N²[2S₂ − U_N + 1]/3, where S₂ = 〈P₂(u·n)〉 and P₂ denotes the second Legendre polynomial. For short times where the magnetic moments are still strongly aligned, U_N ≈ S₂ ≈ 1, we find κ₁ ≈ χ_L and κ₂ ≈ 0 and therefore

τ_N/τ^short_N ≈ e^−χ_L, λ ≲ 1 (42)

For stronger dipolar interactions, the modified mean-field model fails due to significant chain formation. In this regime, we follow the chain-formation model⁴⁹ and assume that the local field can be approximated by the contributions from neighboring particles within a single, perfectly straight chain,

Assuming furthermore that all particles in the chain are still oriented in direction of the ordering field (which is now switched off), u_j = n, we find that we can approximate the local field contribution as u·h_loc ≈ 2λζ(u·n)², whereand the average is performed over all particles i in the chain. We note that ζ is monotonously increasing with increasing chain length and reaches its asymptotic value 2ζ(3) ≈ 2.404 for infinite chain lengths. With this approximation for the chain-formation regime, the cumulants can be evaluated to give κ₁ = 2λζ[1 + 2S₂]/3 andFor short times where the magnetic moments are still strongly aligned, S₄ ≈ S₂ ≈ 1, we find κ₁ ≈ 2λζ and κ₂ ≈ 0 so that

τ_N/τ^short_N ≈ e^−2λζ, 1 ≲ λ ≲ 5 (46)

With the definition of the Langevin susceptibility χ_L = 8λϕ, we can rewrite this expression in the form (18) with b_a = ζ/(4ϕ).

D.2 Effective FFMR relaxation time

D.2.1 Moderate dipolar interactions. In order to estimate the effective relaxation time (17) of the FFMR model, we want to first derive its magnetization relaxation U_N(t) as predicted with the help of the modified mean-field approximation. Our starting point is the magnetization relaxation eqn (39) for weak dipolar interactions,Truncating at first order in χ_L and assuming the simple closure approximation S₂(t) ≈ U_N²(t), we arrive at the nonlinear differential equation τ_aU̇_N + U_N = aU_N³, with the long-time relaxation time τ_a = τ_N/(1 − χ_L/3) and a = 2χ_L/(3 − χ_L). Obviously, the range of validity of the approximation (47) is restricted to χ_L < 3. The solution to this differential equation for the initial condition U_N(0) = 1 readsFor long times, t ≫ (τ_a/2), we indeed find an exponential decay governed by τ_a, U_N ∼ e^−t/τ_a. For short times, t ≪ τ_a, we find from eqn (48)
with the short-time relaxation τ^short_N = τ_N/(1 − χ_L). From eqn (48), we find that the magnetization relaxation deviates stronger from a single-exponential decay the larger χ_L becomes.

With U_N(t) from eqn (48) at hand, we calculate the effective relaxation time from eqn (17), , to find

withand where the last equality holds when a ∈ [0,1]. Since φ(a) = 1 + a/3 + a²/5 + [scr O, script letter O](a³) for a → 0, we find that for weak dipolar interactions

D.2.2 Chain-formation regime. In order to estimate the effective relaxation time in the chain-formation regime of the FFMR model, we assume that the factorization approximation (40) remains valid during the whole relaxation process. Therefore, we considerwhere we truncated the cumulant expansion after the first term withInserting κ₁ into (52) and using the same closure approximation S₂ ≈ U_N² as we did for moderate dipolar interactions, we arrive at a closed ordinary differential equation for U_N, τ_α(dU_N/dt) = −exp(−αU_N²)U_N, where τ_α = τ_N exp[(2/3)λζ] and α = 4λζ/3. The solution for U_N(t) can only be given in implicit form aswhere Ei(x) is the exponential integral. From eqn (54), we find the effective relaxation time by a change of integration variables aswhere erfi(x) denotes the imaginary error function. For very strong dipolar interactions, λ → ∞, this expression becomes [small tau, Greek, macron]_N ≈ τ_N3e^2λζ/(8λζ), whereas in the opposite limit . Thus, we might approximate the corresponding relaxation rate in the form of eqn (19) with (10/9)ζλ = sχ_L or s = 5ζ/(36ϕ). Note, however, that the chain-formation model does not apply for too small values of λ and therefore we expect s = 5ζ/(36ϕ) to overpredict the observed relaxation rate.

E Mean-field arguments for relaxation in the coupled system

E.1 Short-time dynamics

Consider now the magnetization dynamics (26) when both, Brownian and Néel relaxation is present. Combining the corresponding expressions for the RBD and FFMR models (Appendices C and D) for the regime of weak to moderate interactions, where the modified mean-field approximation can be applied, we arrive atFrom eqn (56), we find that the short-time relaxation (where S₂ ≈ 1) is governed by the weighted sum of the RBD and FFMR individual particle rates,Eqn (57) was derived for the case of weak to moderate interactions. We have seen above that eqn (42) provides a better description for the short-time FFMR contribution for larger values of χ_L. Assuming independent Brownian and Néel contributions to the short-time relaxation, we propose to generalise eqn (57) to

E.2 Effective relaxation time

In order to describe the full time evolution, we once again employ the simple closure approximation S₂ ≈ U². With this approximation, eqn (56) can be brought into the form τ^MFU̇ + U = −b₁U³, whereand τ_eff denotes the effective relaxation time for isolated particles, 1/τ_eff = 1/τ_B + 1/τ_N. The coefficient b₁ is given by b₁ = χ_Lτ^MF(1/τ_B − 2/τ_N)/3. Since the time evolution equation is formally identical to the one encountered for the RBD model, we find the same functional form for the magnetisation relaxation,This U(t) exhibits an exponential relaxation with τ^MF at long times (t ≫ τ^MF) and a short-time regime with U(t) = 1 − t/τ^MF_short + [scr O, script letter O](t²). Eqn (60) implies, for χ_L < 3 and (χ_L − 1)τ_B < τ_N, that the effective relaxation time defined in eqn (22) is given by eqn (23), where for negative b₁ the expression can be rewritten with the identity tan⁻¹(ix) = i tanh⁻¹(x). For weak interactions we find from eqn (23)with u₁ = [(2/9)τ_N + (5/9)τ_B]/(τ_N + τ_B) and u₂ = (8τ_N² + 28τ_Nτ_B + 47τ_B²)/[135(τ_N + τ_B)²]. Therefore, to second order in χ_L we findNote that equation reduces for τ_B → ∞ to the FFMR model (61), and for τ_N → ∞ to the result (38) we obtained for the RBD model.

F Dynamic magnetic susceptibility

In order to calculate the dynamic magnetic susceptibility, we first need to find the Laplace-transform of the relaxation function. Using the mean-field result for the coupled system (60), or similarly (48) for the FFMR model, we findwith the hypergeometric function ₂F₁ in the 2nd line. Inserting this expression into [small chi, Greek, tilde](ω) = χ₀[1 − iωC̃(iω)], we find [small chi, Greek, tilde] = χ′ − iχ′′ withup to second order in b₁. For ease of notation, we defined the reduced frequency y = ωτ^MF. It is interesting to note that the mean-field model predicts deviations from the Debye model not in terms of additional Debye modes, but in the above form, with τ^MF, eqn (59) as basic time scale. The whole [small chi, Greek, tilde](ω) behavior is available from the second line of eqn (63), including its low and high frequency characteristics, but the coefficients of an expansion in terms of ω cannot be obtained in closed analytic form; they are available numerically. From eqn (61), we verify that the low-frequency behavior χ′′(ω) = χ₀[small tau, Greek, macron]^MFω + [scr O, script letter O](ω³) is governed by the effective relaxation time [small tau, Greek, macron]^MF, consistent with its definition.

Acknowledgements

M. K. acknowledges support by the Swiss National Science Foundation through grant 200021L_185052 and we thank the Swiss National Supercomputing Centre for providing computational resources.

Notes and references

1. I. Torres-Diaz and C. Rinaldi, Recent progress in ferrofluids research: novel applications of magnetically controllable and tunable fluids, Soft Matter, 2014, 10, 8584–8602.
2. L. Mohammed, H. G. Gomaa, D. Ragab and J. Zhu, Magnetic nanoparticles for environmental and biomedical applications: A review, Particuology, 2017, 30, 1–14.
3. Q. A. Pankhurst, N. T. K. Thanh, S. K. Jones and J. Dobson, Progress in applications of magnetic particles in biomedicine, J. Phys. D: Appl. Phys., 2009, 42, 224001.
4. V. Uskoković, S. Tang and V. M. Wu, Targeted magnetic separation of biomolecules and cells using earthicle-based ferrofluids, Nanoscale, 2019, 11, 11236–11253.
5. R. E. Rosensweig, Ferrohydrodynamics, Cambridge University Press, Cambridge, 1985.
6. D. Shi, M. E. Sadat, A. W. Dunn and D. B. Mast, Photo-fluorescent and magnetic properties of iron oxide nanoparticles for biomedical applications, Nanoscale, 2015, 7, 8209–8232.
7. S. Ota and Y. Takemura, Characterization of Néel and Brownian Relaxations Isolated from Complex Dynamics Influenced by Dipole Interactions in Magnetic Nanoparticles, J. Phys. Chem. C, 2019, 123, 28859–28866.
8. L. Maldonado-Camargo, I. Torres-Diaz, A. Chiu-Lam, M. Hernández and C. Rinaldi, Estimating the contribution of Brownian and Néel relaxation in a magnetic fluid through dynamic magnetic susceptibility measurements, J. Magn. Magn. Mater., 2016, 412, 223–233.
9. Colloidal Magnetic Fluids, ed. S. Odenbach, Springer, Berlin, 2009, vol. 763.
10. E. Blums, A. Cebers and M. M. Maiorov, Magnetic Fluids, de Gruyter, Berlin, 1997.
11. J. Dieckhoff, D. Eberbeck, M. Schilling and F. Ludwig, Magnetic-field dependence of Brownian and Néel relaxation times, J. Appl. Phys., 2016, 119, 043903.
12. S. A. Shah, D. B. Reeves, R. M. Ferguson, J. B. Weaver and K. M. Krishnan, Mixed Brownian alignment and Néel rotations in superparamagnetic iron oxide nanoparticle suspensions driven by an ac field, Phys. Rev. B, 2015, 92, 094438.
13. S. Ota, T. Yamada and Y. Takemura, Dipole-dipole interaction and its concentration dependence of magnetic fluid evaluated by alternating current hysteresis measurement, J. Appl. Phys., 2015, 117, 17D713.
14. H. Rogge, M. Erbe, T. M. Buzug and K. Lüdtke-Buzug, Simulation of the magnetization dynamics of diluted ferrofluids in medical applications, Biomed. Eng., 2013, 58, 1–9.
15. P. Ilg, Equilibrium magnetization and magnetization relaxation of multicore magnetic nanoparticles, Phys. Rev. B, 2017, 95, 214427.
16. Z. Zhao, N. Garraud, D. P. Arnold and C. Rinaldi, Effects of particle diameter and magnetocrystalline anisotropy on magnetic relaxation and magnetic particle imaging performance of magnetic nanoparticles, Phys. Med. Biol., 2020, 65, 025014.
17. M. Shliomis and V. Stepanov, Rotational viscosity of magnetic fluids: contribution of the Brownian and Néel relaxational processes, J. Magn. Magn. Mater., 1993, 122, 196–199.
18. M. Shliomis and V. Stepanov, in Relaxation Phenomena in Condensed Matter, ed. W. Coffey, John Wiley & Sons, 1994, vol. LXXXVII of Advances in Chemical Physics, pp. 1–30.
19. S. Ruta, R. Chantrell and O. Hovorka, Unified model of hyperthermia via hysteresis heating in systems of interacting magnetic nanoparticles, Sci. Rep., 2015, 5, 9090–9097.
20. C. Jonasson, V. Schaller, L. Zeng, E. Olsson, C. Frandsen, A. Castro, L. Nilsson, L. K. Bogart, P. Southern, Q. A. Pankhurst, M. Puerto Morales and C. Johansson, Modelling the effect of different core sizes and magnetic interactions inside magnetic nanoparticles on hyperthermia performance, J. Magn. Magn. Mater., 2019, 477, 198–202.
21. D. V. Berkov, N. L. Gorn and D. Stock, Combined Langevin dynamics/Monte-Carlo simulations of the non-equilibrium ferrofluid remagnetization, J. Magn. Magn. Mater., 2004, 272, E1281–E1283.
22. P. Ilg, Diffusion-jump model for the combined Brownian and Néel relaxation dynamics of ferrofluids in the presence of external fields and flow, Phys. Rev. E, 2019, 100, 022608.
23. C. Iacovita, J. Hurst, G. Manfredi, P. A. Hervieux, B. Donnio, J. L. Gallani and M. V. Rastei, Magnetic force fields of isolated small nanoparticle clusters, Nanoscale, 2020, 12, 1842–1851.
24. K. Wu, D. Su, R. Saha, J. Liu and J.-P. Wang, Investigating the effect of magnetic dipole–dipole interaction on magnetic particle spectroscopy: implications for magnetic nanoparticle-based bioassays and magnetic particle imaging, J. Phys. D: Appl. Phys., 2019, 52, 335002.
25. D. Eberbeck and L. Trahms, Experimental investigation of dipolar interaction in suspensions of magnetic nanoparticles, J. Magn. Magn. Mater., 2011, 323, 1228–1232.
26. P. J. Camp, E. A. Elfimova and A. O. Ivanov, The effects of polydispersity on the initial susceptibilities of ferrofluids, J. Phys.: Condens. Matter, 2014, 26, 456002.
27. J. D. Weeks, D. Chandler and H. Andersen, Role of repulsive forces in determining the equilibrium structure of simple liquids, J. Chem. Phys., 1971, 54, 5237–5247.
28. S. Hess, M. Kröger and H. Voigt, Thermo-mechanical properties of the WCA-Lennard-Jones model system in its fluid and solid states, Physica A, 1998, 250, 58–82.
29. C. Holm and J. J. Weis, The structure of ferrofluids: A status report, Curr. Opin. Colloid Interface Sci., 2005, 10, 133–140.
30. E. A. Elfimova, A. O. Ivanov and P. J. Camp, Static magnetization of immobilized, weakly interacting, superparamagnetic nanoparticles, Nanoscale, 2019, 11, 21834–21846.
31. W. Wernsdorfer, E. B. Orozcho, K. Hasselbach, A. Benoit, B. Barbara, N. Demoncy, A. Loiseau, H. Pascard and D. Mailly, Experimental Evidence of the Néel-Brown Model of Magnetization Reversal, Phys. Rev. Lett., 1997, 78, 1791–1794.
32. P. Ilg and S. Odenbach, in Colloidal Magnetic Fluids: Basics, Development and Applications of Ferrofluids, ed. S. Odenbach, Springer, Berlin, 2008, vol. 763 of Lecture Notes in Physics, pp. 245–315.
33. C. Gardiner, Stochastic Methods, Springer-Verlag, Berlin, Heidelberg, 2009.
34. K. L. Livesey, S. Ruta, N. R. Anderson, D. Baldomir, R. W. Chantrell and D. Serantes, Beyond the blocking model to fit nanoparticle ZFC/FC magnetisation curves, Sci. Rep., 2018, 8, 1–9.
35. R. Erban, J. Chapman and P. Maini, A practical guide to stochastic simulations of reaction-diffusion processes, 2007, arxiv.org, arXiv:0704.1908v2.
36. A. Sreekumari and P. Ilg, Slow relaxation in structure-forming ferrofluids, Phys. Rev. E, 2013, 88, 042315.
37. G. Mériguet, M. Jardat and P. Turq, Brownian dynamics investigation of magnetization and birefringence relaxations in ferrofluids, J. Chem. Phys., 2005, 123, 144915.
38. J. Jordanovic, S. Jaeger and S. H. L. Klapp, Crossover from Normal to Anomalous Diffusion in Systems of Field-Aligned Dipolar Particles, Phys. Rev. Lett., 2011, 106, 038301.
39. Z. Wang, C. Holm and H. W. Müller, Molecular dynamics study on the equilibrium magnetization properties and structure of ferrofluids, Phys. Rev. E, 2002, 66, 021405.
40. P. Ilg and M. Kröger, Anisotropic Self-Diffusion in Ferrofluids studied via Brownian Dynamics Simulations, Phys. Rev. E, 2005, 72, 031504.
41. A. F. Pshenichnikov, V. V. Mekhonoshin and A. Lebedev, Magneto-granulometric analysis of concentrated ferrocolloids, J. Magn. Magn. Mater., 1996, 161, 94–102.
42. J. O. Sindt, P. J. Camp, S. S. Kantorovich, E. A. Elfimova and A. O. Ivanov, Influence of dipolar interactions on the magnetic susceptibility spectra of ferrofluids, Phys. Rev. E, 2016, 93, 063117.
43. B. Huke and M. Lücke, Magnetization of ferrofluids with dipolar interactions: A Born-Mayer expansion, Phys. Rev. E, 2000, 62, 6875–6890.
44. P. Ilg and S. Hess, Nonequilibrium Dynamics and Magnetoviscosity of Moderately Concentrated Magnetic Liquids: A dynamic Mean-field Study, Z. Naturforsch., 2003, 58a, 589–600.
45. L. Rovigatti, J. Russo and F. Sciortino, Structural properties of the dipolar hard-sphere fluid at low temperatures and densities, Soft Matter, 2012, 8, 6310–6319.
46. W. Luo, S. R. Nagel, T. F. Rosenbaum and R. E. Rosensweig, Dipole interactions with random anisotropy in a frozen ferrofluid, Phys. Rev. Lett., 1991, 67, 2721–2724.
47. S. Bedanta, S. Sahoo, X. Chen, W. Kleemann, D. Sudfeld, K. Wojczykowski and A. Hütten, Intra- and interparticle interaction in a dense frozen ferrofluid, Phase Transform., 2005, 78, 121–129.
48. D. Kechrakos and K. N. Trohidou, Magnetic properties of dipolar interacting single-domain particles, Phys. Rev. B, 1998, 58, 12169–12177.
49. A. Y. Zubarev and L. Y. Iskakova, Effect of chainlike aggregates on dynamical properties of magnetic liquids, Phys. Rev. E, 2000, 61, 5415–5421.
50. H. Kesserwan, G. Manfredi, J. Y. Bigot and P. A. Hervieux, Magnetization reversal in isolated and interacting single-domain nanoparticles, Phys. Rev. B, 2011, 84, 172407.

---