Longest relaxation time versus maximum loss peak in the field-dependent longitudinal dynamics of suspended magnetic nanoparticles versus maximum loss peak in the field-dependent longitudinal dynamics of

Magnetic nanoparticles in suspensions provide fascinating model systems to study ﬁeld-induced effects. Their response to external ﬁelds also opens up promising new applications, e.g., in hyperthermia. Despite signiﬁcant research efforts, several basic questions regarding the inﬂuence of external ﬁelds on the magnetization dynamics are still open. Here we revisit the classical model of a suspended magnetic nanoparticle with combined internal and Brownian dynamics in the presence of an external ﬁeld and discuss the ﬁeld-dependent longitudinal relaxation. While internal and Brownian dynamics are independent in the ﬁeld-free case, the coupling of both processes when an external ﬁeld is present leads to richer and more complicated behavior. Using a highly efﬁcient and accurate solver to the underlying Fokker-Planck equation allows us to study a broad parameter range. We identify different dynamical regimes and study their respective properties. In particular, we discuss corrections to the popular rigid-dipole approximation which are captured in terms of a simpliﬁed diffusion-jump model in the Brownian-dominated regime with rare Néel relaxation events. In addition, we discover a regime with surprising mode-coupling effects for magnetically soft nanoparticles. We explain our ﬁndings with the help of a perturbation theory, showing that in this regime the magnetization relaxation at late times is slaved by the slow Brownian motion of the nanoparticle. We discuss consequences of these ﬁndings such as the discrepancy of the longest relaxation time and the inverse frequency of the loss peak of the magnetic susceptibility.


I. INTRODUCTION
Magnetic nanoparticles (MNPs) are not only of great theoretical interest as a physical model system, they have also found numerous practical applications in different areas, ranging from engineering to biomedicine [1][2][3][4]. A particularly interesting property of MNPs that is also relevant in many of these applications is their magnetization dynamics which can readily be manipulated by external fields. For MNPs that are suspended in a viscous liquid, two main magnetization relaxation mechanisms are known: first, rotational Brownian motion of the particle. While this mechanism occurs for any colloidal particle suspended in a viscous liquid, for MNPs there is in addition also the internal magnetization dynamics within the nanoparticle [5,6]. The two corresponding characteristic relaxation times are (i) the Brownian rotational diffusion time τ B and (ii) the Néel relaxation time for internal magnetization reversals τ N . Considering both processes as independent in the absence of an external field (h = 0), Rosensweig defined the effective relaxation time τ eff from the * p.ilg@reading.ac.uk † mk@mat.ethz.ch sum of the corresponding rates [6] 1 Equation (1) is of central importance, both theoretically and for applications, as it allows to easily distinguish different regimes where either Brownian or Néel processes dominate, depending on which relaxation is faster. In relation to hyperthermia applications, for example, Eq. (1) was employed to quantify Brownian and Néel contributions [7] and to interpret the strong size dependence of the specific power absorption of MNPs [8]. Several engineering applications rely on transmitting torques from an external magnetic field to the liquid.
Consequently, for such applications, MNP sizes and magnetic material are chosen such that Brownian processes are dominant [4]. In other applications such as hyperthermia, MNPs are carefully designed to maximize magnetic losses [9]. Some recent experiments that have been performed with the aim to distinguish Brownian and Néel contributions to the magnetization relaxation are therefore of great importance. Clear differences in magnetization hysteresis and AC susceptibility have been observed when measured in fluid and freeze-dried states [10,11], which allows the authors to infer the importance of Brownian relaxation in these cases. Similarly, Brownian and Néel contributions were separated via different frequencies of the oscillating field [12]. It is interesting to note that in magnetorelaxometry, the different magnetization relaxation mechanisms are used to detect binding kinetics due to which Brownian rotation of MNPs becomes suppressed and the magnetization signal changes correspondingly [13]. However, it has been pointed out in the literature that the simple picture underlying Eq. (1) might break down in nonequilibrium situations [14]. Of particular relevance for this study is the finding that the effective relaxation times are strongly dependent on an external magnetic field [15,16]. In particular, the Néel relaxation time can decrease dramatically and therefore can become comparable to Brownian relaxation even when Néel processes would be negligible in the field-free case [17]. Therefore, it is clearly desirable to keep both Brown and internal relaxation processes in the theoretical modeling when considering field-dependent phenomena. Indeed, the corresponding model ("egg" model) is well established [14,18,19] and describes the coupling of both processes, contrary to recent claims [20,21]. Despite having been proposed around 30 years ago, the egg model is not widely used in applied sciences these days [14]. Some noteworthy exceptions in the recent literature include the simulation of hysteresis curves [22], nonlinear response to oscillating magnetic fields [23], frequency mixing [24], and magnetization relaxation [25,26]. We note in particular that simulations of the egg model show pronounced two-step relaxation in multicore magnetic nanoparticles [25] that are confirmed in experiments [26].
Since the egg model describes the internal magnetization dynamics in terms of the stochastic Landau-Lifshitz-Gilbert equation, the basic timescale of this model is often very short compared to typical Brownian relaxation times. This difference in timescales makes the egg model impractical to be used directly for suspended MNPs in most situations [27]. Further comprehensive theoretical work on combined relaxation, coupling Brownian and Néel processes, is still needed [14]. Very recently, more efficient solutions to the egg model based on the corresponding Fokker-Planck (FP) equation have been proposed and explored [28][29][30]. In addition, some heuristic attempts have been made over the last years to eliminate the fast internal dynamics on the microscopic timescale τ D and to model Néel relaxation as thermally activated magnetization reversals, that can be implemented numerically very efficiently via Monte Carlo methods [27,[31][32][33][34]. Very recently, we studied the combined magnetization relaxation for interacting system in the absence of external magnetic fields [35]. There, we found that Brownian and Néel relaxation can be considered as independent processes only for short enough times, whereas coupling of these processes leads to long-time relaxation that can not be described in the form (1), at least for sufficiently strong dipolar interactions.
Here, we study the field-dependent longitudinal magnetization dynamics resulting from combined Brownian and internal relaxation processes for noninteracting MNPs. We use a highly efficient and accurate solver for the FP equation of the egg model to explore a wide range of parameter values. From analyzing the spectrum of eigenvalues, which correspond to the inverse characteristic relaxation times, and magnetic susceptibilities, we identify different dynamical regimes with different characteristic behavior. For magnetically hard nanoparticles with predominant Brownian relaxation and rare, thermally activated internal (Néel) processes, we analyze corrections to the rigid-dipole approximation. Since the egg model in this experimentally relevant regime with τ B , τ N τ D is very inefficient, we compare our results to a simplified diffusion-jump model, recently proposed by one of us [34]. For magnetically soft nanoparticles with long Brownian relaxation times and not too strong fields, instead, we find mode-coupling effects leading to multistep relaxation behavior with slow long-time relaxation. We develop an approximate theory to explain our findings, which gives rather accurate results and demonstrates the interplay between magnetization and particle orientation.
The paper is organized as follows. The formulation of the main model is given in Sec. II. In particular, the egg model of Shliomis and Stepanov is briefly reviewed in Sec. II A. Expressions for the resulting Néel relaxation time are given in Sec. II B. Efficient solution of the FP equation corresponding to the egg model is presented in Sec. III. Next, the approximate diffusion-jump model is briefly reviewed in Sec. IV. In Sec. V, we present a reduced moment approximation that we later use to interpret some of our findings. Solutions to the egg model are presented in Sec. VI together with thorough discussions and analysis in terms of simplified theories. Finally, conclusions are offered in Sec. VII.

A. Microscopic egg model
In order to describe the coupled dynamics of internal magnetization relaxation and particle rotation, Shliomis and Stepanov proposed the so-called "egg" model [19] that we briefly describe here. Let the three-dimensional unit vectors e and n denote the orientation of the particle's magnetic moment and its easy axis, respectively. In the presence of an external field H, the Zeeman energy contributes −me · H to the potential energy, where m denotes the magnitude of the magnetic moment of the nanoparticle. Misalignment of the magnetization direction e from the particle's easy axis n increases the potential energy by the anisotropy contribution −Kv m (e · n) 2 , where K denotes the anisotropy constant of the magnetic material and v m the volume of the magnetic core of the particle.
In the following, it will be convenient to work with dimensionless quantities such as the dimensionless field h = mH/k B T , with h = |h| the Langevin parameter and the dimensionless anisotropy constant κ = Kv m /k B T . Thus, the potential energy of the magnetic nanoparticle relative to the thermal energy k B T can be expressed as where e = e ·ĥ is the component of the normalized magnetic moment parallel to the direction of the applied field h = h/h. Equilibrium properties of magnetic nanoparticles can be inferred from the Boltzmann distribution F eq (e, n) ∼ exp [he + κ (e · n) 2 ] [5]. For dynamical properties we need to consider the time-dependent joint probability density F (e, n; t ). Combining the Landau-Lifshitz-Gilbert model for the internal magnetization dynamics with Brownian rotation of the particle, Shliomis The rotational operators appearing in Eq. (3) are defined as L e = e × ∂/∂e, L n = n × ∂/∂n, and L = L e + L n . For immobile MNPs, τ B → ∞, Eq. (3) reduces to the Landau-Lifshitz-Gilbert magnetization equation, whereas the rigid-dipole approximation of Brownian rotational dynamics is recovered by identifying e and n, i.e., by dropping the internal relaxation and correspondingly neglect the magnetic anisotropy contribution to the interaction. The two basic timescales appearing in Eq. (3) are (i) the Brownian relaxation time τ B = 3ηv h /k B T of rotational diffusion of a spherical particle with v h the hydrodynamic volume in a liquid with viscosity η; and (ii) the internal timescale τ D = m/(2αγ k B T ) of magnetization dynamics with the dimensionless damping parameter α and the gyromagnetic ratio γ . The Larmor frequency is defined by ω L = −(γ /m)∂U/∂e and is given here by ω L = γ H + 2γ KM −1 s (e · n)n, where M s denotes the spontaneous magnetization of the magnetic material M s = m/v m . Note that for the frequently used core-shell particles [4], the presence of a nonmagnetic stabilizing shell implies v m < v h .
The anisotropy constant K for the frequently used magnetite nanoparticles is around K ≈ 10 4 J/m 3 [36], leading to κ ≈ 0.9 and 10.2 for magnetic core diameters of 9 and 20 nm, respectively. The internal timescale τ D , sometimes interpreted as inverse internal attempt frequency, is typically around τ D ∼ 10 −9 . . . 10 −10 s [4,5]. Brownian relaxation times, on the other hand, in typical solvents are several orders of magnitude larger, τ B ∼ 10 −6 . . . 10 −3 s [4,11]. While the egg model provides a physically sound basis to describe coupled Brownian and Néel relaxation [14], in the absence of closed-form analytical solutions, the huge separation in timescales τ D τ B makes the model computationally very demanding to use in many applications. Therefore, the efficient numerical simulation we present in Sec. III as well as the analytical approximations in Secs. IV and V can be very helpful.

B. Néel relaxation time
For sufficiently large anisotropy energies κ 1, the internal magnetization dynamics consists of fast vibrations on the timescale τ D around the particle's easy axis n and slow magnetization reversals once the anisotropy energy barrier is overcome [37]. The latter is governed by the Néel relaxation time τ N , which was calculated approximately by Brown in the absence of a magnetic field and in the limit κ 1 as [18] increasing steeply with κ. The Arrhenius-type formula (4) has been used and tested experimentally [37,38], including detailed investigations on the prefactor (attempt frequency) [39]. For weak and moderate anisotropies κ 2, the expression was derived in Refs. [18,40] with S 2 = P 2 (e · n) eq , P 2 (x) the second Legendre polynomial. For the Boltzmann equilibrium, S 2 can be evaluated explicitly in terms of κ [cf. Eq. (A7) in Appendix A 2]. In the rigid-dipole limit, lim κ→∞ S 2 = 1 and thus e n. For very small κ, the first terms of the expansion of Eq. (5) read as τ N /τ D ≈ 1 + (2/5)κ + (16/175)κ 2 + O(κ 3 ), according to Eq. (A10). A phenomenological approximation that extends Eq. (4) to medium and lower anisotropies is [41] which reduces to Eq. (4) for κ 1 and also provides a good approximation to Eq. (5) for small κ. Strictly speaking, Eq. (6) was derived for longitudinal relaxation when a weak external field is switched off. Figure 1 shows the Néel relaxation time versus the anisotropy parameter κ on a semilogarithmic scale. As mentioned in Refs. [41], Brown's expression (4) breaks down for κ < 2, but Eq. (6) provides a good description over the whole range of anisotropies. For anisotropies in the range κ ∼ 10, we find from Fig. 1 that the difference between Eqs. (4) and (6) is still noticeable, with Eq. (4) systematically underpredicting τ N . Therefore, we here use Eq. (6) rather than the much more popular Eq. (4) to calculate the reference effective relaxation time from Eq. (1) via τ N for given κ.

III. SOLUTION OF FOKKER-PLANCK EQUATION
The time-dependent probability density F (e, n; t ) can be expanded into the complete orthonormal basis set of bipolar harmonics, 134433-3 where K = LM l 1 l 2 is a multi-index and b K the expansion coefficients [28,30]. Admissible multi-indices are all those that satisfy the four inequalities l 1 , l 2 0, |l 1 − l 2 | L l 1 + l 2 , and |M| L simultaneously. The bipolar harmonics are defined in terms of spherical harmonics Y lm and Clebsch-Gordan coefficients C LM l 1 m 1 l 2 m 2 as With the help of the expansion (7), the partial differential equation (3) can be transformed into an infinite system of coupled ordinary differential equations for the coefficients b K , which can be written in vector form as with a positive-definite square matrix A. The additional term d arises due to the conservation of the normalization of F , which corresponds to b 00 00 = 1. For details on the derivation of Eq. (9), see Refs. [28,30]. In the latter reference, we calculate A and d and we also provide the code for implementing Eq. (9). Following Ref. [28], the implementation neglects the Larmor precession term, which should be a good approximation except at very high frequencies.
With the eigenvalues μ of the matrix A, we can write the solution to the linear system (9) as where we specified to isotropic initial condition b K (0) = 0. The weights w Kμ of the eigenvalue μ contributing to b K in this case are w Kμ = ν (S −1 ) Kμ −1 μ S μν d ν , where S −1 contains the eigenvectors of A in its columns. We are going to reserve 1 for the smallest eigenvalue, and use sorted eigenvalues j j+1 when presenting their values. None of the results to be presented below require the explicit calculation of S −1 , as this would be extremely time consuming. Instead, Eq. (9) can be integrated directly using a small time step τ D , the susceptibility requires solving a system of linear equations (22), smallest eigenvalues are obtained using variants of the power method, and stationary solutions of Eq. (9) can be calculated analytically.
Note that the representation (7) can be interpreted as an expansion around the isotropic state. Therefore, in order to accurately capture strongly oriented states which occur for large anisotropies and/or strong magnetic fields, a truncation of the sum at low orders O is inappropriate. The order O of the expansion is defined as the maximum of l 1 + l 2 retained in the sum l 1 + l 2 O. However, the number of terms I (O) in the expansion grows very strongly [like ( 3 14 )O 4 ] with increasing order O, limiting the approach in practice to O ≈ 10 . . . 20, corresponding to 10 4 . . . 10 5 terms.
Here, to study also more strongly ordered states and to be able to cover a wide range of parameters, we use symmetry arguments to reduce the number of expansion coefficients in Eq. (7) considerably. First, we note that the easy axis is described by the director n, which means that n and −n are equivalent. Using the property of the spherical harmonics Y l 2 m 2 (−n) = (−1) l 2 Y l 2 m 2 (n), we conclude that the physically admissible distribution functions F correspond to representations with even l 2 . Although this symmetry reduces the number of expansion coefficients by a factor 2, we need a more substantial reduction to cover a broader parameter range. For initial states that depend only on the difference in the azimuthal angles ϕ e − ϕ n , which corresponds to M = 0 in the bipolar series, we show in Appendix B that the property M = 0 is preserved by the dynamics of the egg model. Given l 2 even and M = 0, we can furthermore conclude that also L − l 1 must be even for the probability density F to be a real function. A detailed discussion of these arguments is provided in Appendix B.
To summarize, the multi-indices in the bipolar expansion (7) to be used here to treat the general case of nonzero external field, and κ = 0, starting from an isotropic or stationary situation, is As a result, if we keep 10 4 or 10 5 terms in the expansion (7), we can now reach, using (B1), O ≈ 42 and 92, respectively, compared to O ≈ 10 without reduction.

IV. DIFFUSION-JUMP MODEL: CORRECTIONS TO THE RIGID-DIPOLE APPROXIMATION
We here focus on the regime of slow Brownian relaxation τ B τ D and of large magnetic anisotropies κ for which τ N τ D . To address this huge timescale gap between the internal timescale τ D and Brownian and Néel relaxation, the egg model presented in Sec. II A needs to be suitably coarse grained.
As mentioned above, for sufficiently large anisotropy energies κ, Néel relaxation can be described as rare, thermally activated magnetization reversals with relaxation time τ N approximately given by Eq. (4). Assuming that deviations of the magnetization direction from the particle's easy axis are negligible on timescales of τ B and τ N , the anisotropy contribution to the internal energy is no longer relevant. Then, the state of individual MNPs can be described by the probability density of the magnetic moment alone, f (e; t ). Keeping Brownian rotational diffusion as in the egg model (3) and modeling Néel relaxation as thermally activated magnetization reversals e → −e on the timescale τ N , one of us proposed to model the corresponding magnetization dynamics in terms of the following diffusion-jump model [34]: involving a yet unspecified rate function r = r(|he |). We emphasize that Eq. (12) was not derived systematically from the egg model, but only plausibility arguments are offered. In the limit τ N → ∞, Néel relaxation can be ignored and Eq. (12) reduces to the well-known rigid-dipole model introduced by Martsenyuk and coworkers [42]. For large but finite τ N , Eq. (12) extends the rigid-dipole model by including Néel processes in an approximate manner. Assuming that Néel relaxation can be modeled as a Poisson process and imposing Boltzmann equilibrium f eq (e) ∼ exp(he ) as stationary state 134433-4 specifies Eq. (12) only up to the unknown rate function r(x), with r(0) = 1 and r 0. The condition r(0) = 1 ensures Eq. (12) leads to the correct relaxation time, Eq. (1), in the field-free case. In Ref. [34], some choices for this rate function were discussed, e.g., r(x) = 1 corresponding to Arrhenius activation and r(x) = 1 + tanh(x) leading to Glauber dynamics. We can use the diffusion-jump model to derive approximate equations for the magnetization dynamics. First, define the time-dependent normalized magnetization m(t ) as m = e f (e; t )de = e , where the integration is performed over the three-dimensional unit sphere. Multiplying the kinetic equation (12) by e and subsequent integration, we arrive at [34] Equation (13) still contains averages over the unknown probability density f and therefore does not provide a closed-form magnetization equation. However, we can solve Eq. (13) approximately for the late stages of relaxation, leading to the relaxation times parallel and perpendicular to the applied field 1 where L 1 (h) = coth(h) − 1/h denotes the Langevin function and L 1 its first derivative. Details of the derivation are provided in Appendix C, where also expressions for I 2 , R 2 are given.
For τ N → ∞, we recover from Eq. (14) the known results for the rigid-dipole approximation of the relaxation times within the effective field approximation [42]. For large but finite τ N , the terms proportional to τ −1 N in Eq. (14) describe corrections to the rigid-dipole approximation. As discussed in Appendix C, these contributions depend on the particular form of rate function r appearing in Eq. (12).

V. REDUCED MOMENT APPROXIMATION
We here propose and analyze an approximate theory to capture and explain unexpected mode-coupling effects in the relaxation phenomena of magnetically soft nanoparticles. For the relaxation processes of interest, we can limit ourselves to small deviations from equilibrium, F = F eq (1 + ) with | | 1. Crucially, we assume that these small deviations can be expressed as a linear combination of N given functions X j = X j (e, n). Introducing the short notation δ X j = X j − X j eq , we can derive closed moment equations where the matrices B = XX eq − X eq X eq and M = (2τ B ) −1 (LX) T · LX eq + (2τ D ) −1 (L e X) T · L e X eq are derived using component notation in Appendix E, Eqs. (E3) and (E5), respectively. Both matrices can be evaluated with the help of the known F eq mentioned after Eq. (2). The crucial step is to choose a suitable set of functions X 1 , X 2 , . . . , X N . In the simplest approach, we include only the magnetization component in the field direction, i.e., N = 1 and X 1 = e . In this case, Eq. (16) describes a single exponential relaxation with time constant τ 10 For h = 0, τ 10 reduces to τ 0 , whereas for negligible internal relaxation we recover the field-dependent longitudinal relaxation in the rigid-dipole approximation [42]. It is interesting to note that in this model, considering uncoupled magnetization relaxation implies ignoring anisotropy effects.
To obtain a better approximation and motivated by earlier works on the rotational viscosity by Shliomis and Stepanov [43], we suggest to consider N = 3 with the following three functions to describe the coupling of the magnetization to the easy axis: Note that due to the symmetry n → −n, the easy-axis orientation enters only quadratically. Equation (16) represents a closed system of N-coupled time evolution equations, where the closure is obtained by representing as a linear combination of X j [see Eq. (E1)]. Compared to the infinite system of Eq. (9) for the expansion coefficients b K of bipolar harmonics, Eq. (16) represents a very drastic simplification. The accuracy and limits of validity of this drastically reduced set of moments will be discussed below. Although Eq. (16) is formally identical to Eq. (9) where M · B −1 plays the role of the matrix A and M · B −1 · X eq takes the role of d, it is worth noting that the reduced moment approximation (RMA) can be seen as a low-order perturbation around the Boltzmann equilibrium state F eq , whereas the series of bipolar harmonics (7) is a systematic expansion around the fully isotropic state. Note that X j can alternatively be expressed in terms of bipolar harmonics Y K , but for convenience we will use the representation given in Eq. (17). Upon diagonalizing the matrix M · B −1 , we can determine the corresponding three eigenvalues. Explicit analytic expression for these eigenvalues in the limit of weak fields h and small anisotropies κ are obtained in Appendix E using secondorder perturbation theory. In the absence of a magnetic field, these eigenvalues reduce to and Note that for immobile particles, τ B → ∞, the eigenvalue λ 2 for h = 0 describes the internal relaxation time [Eq. While we are using the convention that the eigenvalues j obtained via the FP approach are sorted, j j+1 , the numbering of the eigenvalues λ j for the RMA corresponds to the choice of moments (17) to allow for a more direct interpretation. For τ B > 2τ D , λ 3 is the smallest of the three 134433-5 RMA eigenvalues at zero field, and the RMA eigenvalues may cross upon variation of h and κ, as we will find later.
For weak magnetic fields, we find that first-order contributions vanish so that the eigenvalues vary quadratically with h. For λ 1 and λ 2 , we find a quadratic increase, corresponding to the familiar decrease of relaxation times with increasing field strength. This, however, is not true for λ 3 , which instead decreases [Eq. (F43)], implying an increasing relaxation time with increasing h. We discuss this surprising prediction and its consequences in Sec. VI C.

A. Field-free case
In the absence of an external magnetic field, internal and Brownian relaxation decouple and the faster of the two dominates the effective relaxation time τ eff defined by the well-established relation (1), as proposed in Refs. [6,19]. Figure 2(a) shows τ eff /τ D versus the anisotropy parameter κ, where we define τ eff = −1 1 (h = 0) with 1 (h = 0) the smallest eigenvalue of the matrix A, evaluated for zero field that has nonzero weight. The results shown are obtained for order O = 40 of the bipolar expansion presented in Eq. (7), but are indistinguishable when larger O are used. From Fig. 2 we find that Eq. (1) provides an excellent description, confirming the decoupling of the relaxation processes at zero field. The same conclusions were reached in Ref. [44] with the help of Brownian dynamics simulations of the egg model and extracting τ eff by fitting an exponential decay for the longtime relaxation. It is reassuring to validate that the smallest eigenvalue corresponds to the long-time relaxation observed in these simulations.
As a side comment, we remark that the agreement of the lowest eigenvalue with Eq. (1) is nearly perfect when the interpolation formula (6) is used for the Néel relaxation time, whereas some mild discrepancies at intermediate values of κ can be seen in Fig. 2(a) for Brown's formula (4).
In Fig. 2(b) we scale the data differently and show τ eff /τ B versus the ratio of the relaxation times τ N /τ B , where τ N is calculated from Eq. (6). We observe an almost perfect data collapse on a master curve, again in excellent agreement with Eq. (1). Therefore, relaxation in the field-free case can be considered to be well understood. As we will see in the following, the situation is much more complicated when an external field is present.

B. Dynamical phase diagram for field-dependent relaxation
In the presence of an external magnetic field, the dynamics of the egg model presented in Sec. II A is determined by the three dimensionless parameters κ, τ B /τ D , and h. This parameter space is too large for a brute-force exploration. Therefore, we start our analysis by considering some representative choices of parameter values. Figure 3 shows the magnetization relaxation parallel to an applied field (i.e., longitudinal) for different values of κ, τ B /τ D , and h. High-order solutions (O = 100) to Eq. (9) are compared to results from Brownian dynamics simulations of the egg model. The latter are described in detail, e.g., in (1) and τ N from Eq. (6). Dashed line is same as solid but using Eq. (4) instead of (6) for τ N . Note: eigenvalues corresponding to eigenvectors with vanishing weight w Kμ in Eq. (10) are not taken into account. This is equivalent with considering the subset |l 1 − l 2 | = L and L ∈ {0, 1} of (11), according to [30].
Refs. [30,34]. We solve Eq. (9) for isotropic initial condition and plot the transient magnetization e (t ) in Fig. 3(a) and 1 − e (t )/m eq in 3(b), where m eq = L 1 (h) denotes the equilibrium magnetization. Brownian dynamics simulations are performed for the same conditions, with ensemble averages performed over 10 6 independent realizations. The overall good agreement between both methods serves as a further test for the correct implementation of the numerical schemes. Upon closer inspection, some deviations between solutions of the FP equation and Brownian dynamics simulations are seen at short times. These deviations are stronger for large values of κ and h and are likely due to the Larmor precession term being neglected in the FP solution. Due to noise inherent in Brownian dynamics simulations, we find that the high-order implementation of the bipolar harmonics expansion allows us to extend the observable range and explore much longer timescales where the magnetization is merely a fraction 10 −4 away from the stationary state. In all cases, we find from Fig. 3 the approach to the equilibrium value proceeding in multiple stages. The first stage is always a very rapid decay on timescales shorter than τ D . For small values of the anisotropy constant κ, the initial decay gives way to the terminal relaxation for times greater than τ D . This regime is therefore characterized by a two-step relaxation. For larger values of κ, the situation is more complicated, due to the occurrence of an intermediate relaxation regime around τ D , before the late-stage relaxation set in for times greater than 10 . . . 100τ D . Therefore, strictly speaking, this regime shows a three-step relaxation. For practical purposes, however, the short-time relaxation is often negligible and/or difficult to observe. Similar comments apply to the late-stage relaxation in cases where the magnetization has relaxed almost fully, before entering the terminal regime. Therefore, we next consider the susceptibility spectra to find an alternative and possibly more practical distinction between different dynamical regimes.
Magnetization relaxation can conveniently be studied experimentally in terms of the AC magnetic susceptibility spectra, spanning a broad range of frequencies [11,15,45,46]. The complex magnetic susceptibility χ * (ω) is defined as the response function to an oscillating external field. More precisely, we here consider the normalized longitudinal susceptibility χ * defined via e (ω) = (3χ L ) −1 χ * (ω)h 1 (ω), with e (ω) the Fourier transform of the normalized magnetization and χ L the Langevin susceptibility. The amplitude of the oscillating field is weak, |h 1 | 1, to ensure the system is in the linear response regime. Note that χ * quantifies the response to a weak oscillating field with frequency ω in addition to a collinear static field with magnitude h. Therefore, χ * depends not only on the frequency ω and the model parameters κ, τ B /τ D , but also on h.
In terms of the bipolar harmonics expansion presented in Sec. III, we find where δb K = b K − b eq K denotes the deviation from the equilibrium value. Inserting h(t ) = h + h 1 e iωt into Eq. (9) and performing a first-order expansion in terms of h 1 allows us to calculate the component δb 10 10 from with the subscript zero denoting the quantities for h 1 = 0. Further details on the calculation of the complex susceptibility are given in Ref. [30]. Figure 4 shows the real and imaginary parts of the complex susceptibility, χ and χ , respectively, as a function of frequency ω of the applied field. To better compare the susceptibilities at different field strengths, we normalize χ and χ with the static (zero-frequency) susceptibility χ (0) = 3χ L L 1 (h). For small κ, we find a Debye-type behavior with a pronounced loss peak for frequencies around τ D or larger, in agreement with earlier works [28]. Applying a magnetic field in the range h 5 does not change the behavior qualitatively. Figure 3 shows that for these parameters, the late stage relaxation sets in only after the magnetization has almost reached its equilibrium value. Thus, the weight of the longest relaxation time is small and so is its contribution to the susceptibility, where a shoulder at lower frequencies is barely visible. Therefore, we observe a near-Debye behavior in this regime, despite the two-step relaxation of the magnetization. For larger κ, however, χ shows a pronounced double-peak structure, with the high-frequency peak indicating the shorttime relaxation on timescales τ D , whereas the low-frequency peak mirrors the long-time relaxation. In this case, the doublepeak structure mirrors the multistep relaxation behavior seen in Fig. 3, where the long-time relaxation carries a significant weight and sets in when the magnetization is still noticeably different from its stationary value. As the field strength h increases, the location of the low-frequency peak moves to higher frequency, indicating a speeding up of the relaxation.
Before analyzing the relaxation further in more detail in Secs. VI C and VI D, we first classify the different relaxation behavior based on the forms of the susceptibility spectra. More specifically, we study the corresponding Cole-Cole plots, where χ (ω) is shown parametrically versus χ (ω) [15]. Depending on the values of the model parameters, we observe three different characteristic shapes that we denote by o, Oo, and oO. An overview of these shapes depending on the value of the anisotropy constant κ and field strength h is shown in Fig. 5. The plot was created for τ B /τ D = 100, but looks qualitatively identical for smaller τ B as long as τ B /τ D 3. Note that the colored curves in Fig. 5 correspond to the same colored curves in Fig. 4. In the regime of weak magnetic anisotropy, i.e., for κ < κ * (denoted as "o"), the Cole-Cole plot resembles a semicircle, which is characteristic of Debye-type relaxation. Note that around κ * occurs the transition between a more diffusive relaxation [Eq. (5)] and thermally activated Néel-type relaxation [Eq. (4)]. In the regime of Néel relaxation, i.e., for κ > κ * , we find significant deviations from the semicircle shape in the Cole-Cole plots showing two maxima. For small enough fields, h κ − κ * , we find that the dominant loss peak occurs at low frequencies.
In this regime (denoted as "oO"), we expect the dominant loss peak to coincide with the smallest eigenvalue and thus the long-time relaxation time. Finally, for strong anisotropies (κ > κ * ) and strong fields (h > κ − κ * ), we observe that the shape of the Cole-Cole plots change in that the dominant loss peak moves to high frequencies. In this regime (denoted by "Oo"), the long-time relaxation corresponds to the more or less pronounced secondary peak, whereas the main loss peak is determined by the short-time relaxation. This surprising finding is important for correctly interpreting magnetic susceptibility spectra and associating loss peaks with relaxation times.
Figures 6(a) and 6(b) show contour plots of the smallest eigenvalue 1 and the ratio ω * / 1 , respectively, where ω * denotes the frequency of the largest loss peak. The eigenvalue 1 is scaled with τ B /3, the inverse of the zero-field eigenvalue (20). Comparing Fig. 6(a) to Fig. 5, we find the small-κ regime ("o") to roughly coincide with the parameter region in which the smallest eigenvalue 1 ≈ 3/τ B . For larger values of κ, the regions Oo and oO identified in Fig. 5 correspond roughly to regimes where 1 > 3/τ B and 1 < 3/τ B , respectively. In the regime of large κ and small enough fields (oO region), we find from Fig. 6(b) that the peak frequency ω * is close to the the smallest eigenvalue 1 . When comparing Fig. 6(b) to Fig. 5, we also note that the transition from the Oo to oO region is mirrored by a large change in the characteristic frequency ω * , whereas 1 changes only mildly. While ω * remains near 1 in the oO region, ω * jumps to much higher frequencies upon entering the Oo region. In the regime of small κ, as discussed above, the smallest eigenvalue carries only a small weight, so that the Debye-type susceptibility is peaked at much higher frequencies. This behavior is evident in Fig. 6(b) where we also notice a strong increase in ω * with increasing field strength h for κ < κ * , while 1 stays approximately constant.

C. Mode-coupling effects
In this section, we aim to better understand the regime of magnetically weak nanoparticles with moderate values of the anisotropy constant κ. For such nanoparticles, there is no major energy barrier for deviations of the magnetization direction e from the magnetocrystalline easy axis n. Therefore, internal relaxation is mainly diffusive and the corresponding relaxation times in zero field [Eq. (5)] increase only mildly with κ. Figure 3 shows exemplary relaxation curves for this regime with κ = 1 < κ * . As discussed in Sec. VI B, the relaxation 134433-8 is characterized by a fast initial relaxation, that brings the magnetization very close to its equilibrium value, before the slow, long-time relaxation sets in. The inherent noise makes this long-time relaxation difficult to study with Brownian dynamics simulations. However, with our high-order solution of the FP equation, the long-time relaxation can be investigated in great detail.
Although the ultimate relaxation is exponential and allows to define the long-time relaxation time, it is obvious from Fig. 3 that the fast initial relaxation is very significant in this regime. Therefore, defining a single effective relaxation time is of limited use here. Instead, we aim to better understand the reason for the two-step relaxation and to predict the corresponding two characteristic relaxation times using the reduced moment approximation presented in Sec. V. Figure 7 shows the smallest eigenvalues j obtained from high-order solutions to the FP equation. More precisely, the eigenvalues are obtained from diagonalization of the matrix A in Eq. (9) for order O = 50 and then sorted according to their values. We observe that the second and third smallest eigenvalues increase with increasing field strength h, corresponding to the familiar decrease of relaxation times due to applied magnetic field. The lowest eigenvalue, however, remains almost unaffected by the field and tends to slightly decrease as h increases. Furthermore, we observe that the lowest eigenvalue is approximately independent of κ, whereas the second and third smallest eigenvalues increase and decrease with increasing κ, respectively. For τ B = 5τ D and h = 1, we observe crossings of the smallest eigenvalues around κ = 1.5 and 2.5. Finally, for fixed κ and h, the smallest eigenvalues decrease with increasing τ B /τ D indicating the slowing down of the relaxation.
For the parameter regime investigated here, the RMA despite its simplicity provides surprisingly accurate predictions of the three lowest eigenvalues over a considerable range of parameters. Due to the degeneracy of the smallest eigenvalues, the analytical expressions for the RMA for weak fields and anisotropies are not defined for the special case τ B = 2τ D (see Appendix F), while the full solution of the RMA is not effected [cf. Fig. 7(c)]. The second lowest eigenvalue corresponding to the magnetization and the third lowest eigenvalue associated with X 2 both increase with increasing h. The lowest eigenvalue corresponds to the easy-axis orientation n 2 . Contrary to the other two, this eigenvalue slightly decreases with increasing h and is approximately independent of κ. See also Eq. (F43) where we worked out the leading quadratic order. Therefore, we can conclude that the two-step relaxation found in Fig. 3 arises due to the coupling of the magnetization dynamics to the relaxation of the easy-axis orientation. In the regime of weak magnetic anisotropies and large enough Brownian relaxation times, the easy-axis relaxation is the limiting process occurring at long times. Depending on the parameter values and sensitivity of observation, the longtime behavior might be difficult to detect in practice in the magnetization relaxation or magnetic susceptibility spectra. Since the slow process is linked to the dynamics of the easy axis, measurements detecting particle rotations, e.g., due to attached molecules, should be able to directly observe the slow relaxation.
In Fig. 8 we plot the smallest eigenvalues versus the ratio τ N /τ B and scale the eigenvalues with the effective relaxation time τ eff in zero field [Eq. (1)]. The remaining parameters are chosen as τ B = 100τ D and h = 0.1 (a) and h = 1 (b), respectively. Note that τ N /τ B τ D /τ B since τ N τ D . The gray shaded area indicates the regime λτ eff < 1 or also τ eff < 1. For τ N /τ B < 1 2 , we indeed find the smallest eigenvalue 3/τ B [Eq. (20)] to be smaller than 1/τ eff , meaning the longest relaxation time in the presence of a weak applied field exceeds τ eff . We note that this behavior is contrary to the field-induced speeding up of magnetization relaxation typically observed for magnetically hard nanoparticles [42]. Although the parameters chosen in Fig. 8 are outside the range of validity of the RMA, the model provides still a good description of the smallest eigenvalue up to τ N /τ B ≈ 1 2 . For τ N /τ B > 1 2 , the lowest eigenvalue obtained from the FP approach is given by field-induced corrections to 1/τ eff , leading to the familiar scenario of field-assisted accelerated magnetization relaxation. It 134433-9 is interesting to note that the crossing of the lowest eigenvalues for h = 0.1 is avoided for h = 1.

D. Corrections to the rigid-dipole approximation
Magnetically hard nanoparticles are characterized by strong magnetic anisotropies κ 1, so that the magnetic moment can be considered to be thermally blocked near the easy-axis direction. For not too strong fields h κ, we encounter the classical Néel scenario of internal relaxations as rare, thermally activated magnetization reversals [18] where relaxation times increase very strongly with κ, approximately given by Eq. (4).
Since Brownian motion of the nanoparticles is often slow compared to the internal relaxation, this regime is characterized by a strong timescale separation, τ N , τ B τ D . Consequently, we expect a two-step relaxation, where a very The gray shaded area indicates the regime τ eff < 1 or also λτ eff < 1. Within the shown regime, the RMA eigenvalues, derived in Appendix F, are actually sorted as λ 1 λ 2 λ 3 .
short and fast initial process is followed by a pronounced long-time relaxation with corresponding relaxation time τ . Such a behavior is indeed observed in the regime oO of Fig. 5 and can exemplarily be seen in Fig. 3 for τ B = 100τ D , κ = 5, h = 1. In this regime, the inverse of the smallest eigenvalue −1 1 agrees well with τ extracted from an exponential fit to the long-time relaxation. Figure 9(a) shows the long-time relaxation time −1 1 as a function of the applied field h for τ B = 100τ D and various values of κ. As expected, we find that τ increases with κ due to the increasing Néel relaxation time. It is interesting to note that deviations from the rigid-dipole approximations are clearly visible even for κ = 10 where τ N is significantly larger than τ B . We observe that these deviations are well described by the diffusion-jump model of Sec. IV [Eq. (14)], down to κ 8 where τ N ≈ τ B . The good comparison seen in Fig. 9(a) is obtained for Arrhenius rates r = 1, where the relaxation time is given by In Appendix D we argue that r = 1 is indeed the appropriate choice for large magnetic anisotropies that recovers Brown's result for the Néel relaxation of immobile particles. In the parameter range investigated here, we find the familiar decrease of τ with increasing field strength h. In particular, Eq. (23) predicts Brownian processes to dominate with a decrease τ ≈ τ B /h for large fields, in agreement with arguments put forward in Ref. [47]. We want to emphasize that this finding results from the coupling of internal and Brownian relaxation since otherwise the Néel relaxation time, Eq. (D6) or (D5), would dominate as long as h κ. While the agreement between lowest eigenvalues and Eq. (23) is excellent in this regime up to moderate field strengths, noticeable differences occur for stronger fields, where the inverse of the smallest eigenvalue even exceeds the corresponding rigid-dipole result. The transition to this strong-field regime occurs at higher and higher field strengths as κ increases (see Fig. 5). For smaller values of κ, a qualitatively different behavior occurs, where the relaxation times first increase with h before decreasing at 134433-10 stronger field strengths, as visible in Fig. 9 for κ = 6, 7. This behavior is far from the rigid-dipole limit and is not captured by Eq. (23) since Néel processes dominate Brownian rotation in this regime. We did not explore whether Eq. (14) with a different choice of the rate function r might be able to describe this effect.
In Fig. 9(b), we again show τ /τ D for τ B = 100τ D but versus the ratio τ N /τ B . For fixed value of the field strength h, we observe an increase of τ with τ N /τ B , slowly approaching the rigid-dipole limit. Again, the prediction (23) of the diffusion-jump model provides a very good description for sufficiently large τ N and up to moderate field strengths. For stronger fields, Eq. (23) somewhat underestimates the relaxation time. An important result is the very slow approach of τ to the rigid-dipole limit, where strong deviations at low field strengths persist even for large τ N . It is interesting to note that in the regime τ N < τ B , for which the diffusion-jump model was not designed originally, Eq. (23) still captures the second lowest eigenvalue rather well.

VII. DISCUSSION
We present an efficient implementation of the classical egg model [19], describing the coupled magnetization and particle rotational dynamics of suspended magnetic nanoparticles. Exploiting several symmetries, we are able to retain a sufficiently large number of terms in the bipolar harmonics expansion (7) to study the longitudinal dynamics for a broad range of parameters, including strongly oriented states arising for large magnetic anisotropies κ and strong external fields h.
In the field-free case, we verify that Eq. (1) provides an excellent description of the lowest eigenvalue of the Fokker-Planck equation, confirming earlier results [44] for the longtime relaxation time. From the Fokker-Planck equation (3), one finds that internal dynamics and particle rotation decouple for h = 0, justifying their independent contribution to Eq. (1).
Internal and Brownian relaxation are no longer independent if an applied magnetic field is present, leading to richer but also more complicated behavior. Studying the egg model over a broad range of parameters, we identify three different dynamical regimes, which we denote as o, oO, and Oo according to their respective Cole-Cole plots. For large enough Brownian relaxation times, τ B 3τ D , we find that these regimes are mostly determined by the anisotropy parameter κ and the dimensionless field strength h, and do not change significantly upon increasing τ B .
The regime "o" comprises magnetically soft nanoparticles with κ < κ * and κ * ≈ 2. Soft magnetic nanoparticles are very promising due to their efficient energy dissipation [48,49]. In this regime, the AC magnetic susceptibility is near Debye, showing the characteristic semicircle in the Cole-Cole plot. In this regime, we find surprising, field-induced mode-coupling effects with different dependencies of the eigenvalues on the model parameters, leading to crossings and avoided crossings of eigenvalues as, e.g., the field strength h or the anisotropy parameter κ is varied. We develop a simplified theory (RMA), that allows us to understand these effects quantitatively. In addition, RMA enables us to conclude that in this regime, the easy-axis orientation corresponds to the smallest eigenvalue and couples to the magnetization relaxation, leading to slow, long-time relaxation. In practice, the slowest mode might be difficult to detect in the magnetization dynamics due to its small weight. Measurements involving particle rotation, however, e.g., via markers attached to the nanoparticle surface, rotational friction, or orientation effects of nonspherical nanoparticles, should reveal the slow relaxation mode. The effect should already be visible at low field strengths. For interpreting AC magnetic spectra, it is important to note that the frequency of the loss peak ω * is much higher than the smallest eigenvalue 1 .
For magnetically hard particles κ > κ * and not too strong fields h < κ − κ * , we enter the regime "oO" that is characterized by predominant Brownian particle rotation with rare additional, thermally activated Néel relaxation events. In this regime, the smallest eigenvalue closely corresponds to  the dominant loss peak. A recently proposed diffusion-jump model [34] predicts the corresponding longest relaxation time rather well for κ 8 [Eq. (23)], capturing the corrections to the rigid-dipole approximation. It is important to stress that corrections to the rigid-dipole limit are clearly discernible in weak fields even for large anisotropies. These encouraging results together with the high efficiency of the diffusion-jump model open new possibilities for many-body Brownian dynamics simulations of ferrofluids, e.g., to study their structure formation in oscillating fields, which is found to depend on the Brownian as well as the Néel relaxation time [50]. For intermediate values of the anisotropy constant κ * κ 7, the diffusion-jump model breaks down and we discover an unexpected increase of the longest relaxation time with increasing field strength, reaching a maximum near h ≈ 2, before decreasing for stronger fields.
Finally, for magnetically hard nanoparticles in strong fields, κ > κ * and h > κ − κ * , which we term the "Oo" regime, the smallest eigenvalue no longer corresponds to the dominant loss peak, which is now located at much higher frequencies. While a discrepancy between the smallest eigenvalue 1 and the peak frequency ω * was already found for the special case of immobile MNPs [51,52], we here quantify these discrepancies and provide the dynamical phase diagram when Brownian relaxation is also present. Knowledge of these discrepancies is also very useful for the correct interpretation of AC spectra and designing suitable magnetic nanoparticles for hyperthermia applications.
We emphasize that the significant reduction in the terms retained in the bipolar harmonics expansion (7) used here (see Appendix B) is limited to studying longitudinal relaxation. Perpendicular magnetization components require the use of a larger set of bipolar harmonics. For this more general case, one can resort to our previous work [30].
We want to conclude with another potential application of this work to the promising new technique of magnetic particle imaging [53]. For improving the spatial resolution, the dynamical properties of magnetic nanoparticles and Néel relaxation at high frequencies are crucial [54]. For ultrahigh frequencies corresponding to very short times, however, the egg model needs to be suitably extended. One step in this direction could be including inertia effects [26], as has already been proposed in the original paper [19]. Note that for such short times, details of the phenomenological Gilbert damping become important [55], highlighting the general problem of extending coarse-grained models to shorter times, where additional processes (such as deviations from homogeneous internal magnetization) might come into play. This interesting field of research is left as a future challenge.

APPENDIX A: ORDER PARAMETERS
For frequent use, we here collect a number of properties of orientational order parameters.

Higher-order Langevin function
The order parameters L j are defined by the equilibrium averages with the jth-order Legendre polynomials P j (z). They quantify the orientational ordering of the magnetization parallel to the direction of the applied field e = e ·ĥ. Note that the functions L j depend only on h and are independent of κ. From the recursion relation of the Legendre polynomials, the jth-order Langevin function L j is defined recursively as [56] with L 0 = 1 and L 1 (h) = coth(h) − h −1 being the classical Langevin function. The jth-order Langevin function exhibits the following asymptotic properties at small and large (semipositive) h, respectively, which can be derived with [57] at hand. All L j (h) monotonically increase with increasing h for h 0. Starting from L j (0) = 0, they all approach unity in the h → ∞ limit. Expressions for, and recursive relationships between, the derivatives of the higher-order Langevin functions and their behavior can be derived based on Eqs. (A2)-(A4). In particular, we use the relation and L j (0) = δ j,0 /3, according to Eq. (A3).

Alignment order parameters
Order parameters S are related to the th Legendre polynomials as equilibrium averages of the magnetization component parallel to n: Because n is a director, S is zero for odd . For the Boltzmann equilibrium, S 2 can be evaluated explicitly to give while S 0 = 1 and where Da is the Dawson integral Da(x) = e −x 2 x 0 e −y 2 dy = −(i/2) √ π e −x 2 erf(ix) [30,34]. Higher-order parameters can be recursively obtained from S 0 and S 2 via [30] Hence, for example, with S 2 from (A7). Note that S depend only on κ and are independent of h. For even the asymptotic behaviors are given, with the help of lim x→∞ 4x Da(x) = 2, by For 4 the expansion (A10) simplifies to S 2 = 2κ/15 + O(κ 2 ) and

APPENDIX B: BIPOLAR HARMONICS EXPANSION
Here, we give more details on the expansion of the joint probability density into bipolar harmonics. In Sec. III we have already argued that we can restrict the expansion (7) to terms with even l 2 due to the n → −n symmetry.
To further reduce the number of expansion coefficients, we make additional, more restrictive assumptions. First, we consider the assumption proposed by Titov et al. [29] that the joint probability density F (e, n; t ) does not depend on the azimuthal angles ϕ e and ϕ n separately, but only on their difference ϕ e − ϕ n . In terms of the bipolar harmonic expansion (7), since Y l 1 m 1 (e) ∼ e im 1 ϕ e , Y l 2 m 2 (n) ∼ e im 2 ϕ n , this implies a restriction to terms with m 2 = −m 1 , i.e., M = 0 only. In Ref. [29], it is argued that this property of the interaction potential should somehow be inherited by F , but no justification for the claim was given. In Ref. [30], we have already shown that terms b LM l 1 l 2 corresponding to different values of M do not couple due to the dynamics of the egg model. Therefore, if initial conditions are such that moments with M = 0 are not present at time zero, they will remain zero also for later times (see Sec. 4.2 in [30]). This demonstrates that the assumption M = 0 is indeed justified if the initial state also corresponds to M = 0. Such cases include isotropic initial conditions or stationary states (as they can be produced from isotropic initial conditions). The number of moments with M = 0 and even l 2 at given order O is already given in Eq. (E.1) of [30].
Upon closer inspection, we find that a further reduction in the number of relevant terms is possible by the condition that L − l 1 is even. This condition can be inferred from the symmetry properties of the Clebsch-Gordan coefficients C LM l 1 m 1 l 2 m 2 = (−1) l 1 +l 2 −L C L,−M l 1 ,−m 1 ,l 2 ,−m 2 , or more specifically for our case with M = 0 this becomes C L0 l 1 ml 2 ,−m = (−1) l 1 +l 2 −L C L,0 l 1 ,−m,l 2 ,m . These symmetry properties are conveniently derived by converting the Clebsch-Gordan coefficients to Wigner 3 j symbols, and using basic symmetry properties of the 3 j symbols [58]. Since we already concluded that l 2 is even, the condition L − l 1 even ensures via Eq. (8) that F is a real function where the dependence on the azimuthal angles is given by cos[m(ϕ e − ϕ n )]. With these reductions, only a limited set of base functions needs to be considered, significantly reducing the number of components I for a given order O.
Since moments with odd L − l 1 are not excited as well, we conclude that all stationary moments with odd L − l 1 , nonzero M, and odd l 2 vanish, giving rise to the set (11). The remaining number of moments (including b 00 00 ) is Compared to the corresponding expression in Ref. [30],

APPENDIX C: RELAXATION TIMES OF DIFFUSION-JUMP MODEL
The terminal relaxation times we are interested in govern the late-stage approach to equilibrium. Therefore, to calculate effective long-time relaxation times, we consider small deviations from equilibrium and make the ansatz where a is unknown and f eq = Z −1 0 exp(he ) denotes the equilibrium density with Z 0 = 4π sinh(h)/h. We introduce the short notation A eq = A f eq de to indicate averages over the equilibrium density. The ansatz (C1) satisfies f (e; t )de = 1, but is valid only for weak deviations from equilibrium, i.e., |a| 1. We note that replacing Eq. (C1) by the effective field approximation (EFA) allows to derive closed-form magnetization equations also far from equilibrium [42]. We have used EFA to derive the magnetization equation for this model in Ref. [34] with an emphasis on the perpendicular relaxation. However, since EFA involves stronger assumptions and reduces to Eq. (C1) near equilibrium, we here proceed without using EFA.
With the ansatz (C1), all moments can be expressed in terms of the unknown vector a, A = A eq + ( Ae eq − A eq e eq ) · a. (C2) In particular, the nonequilibrium magnetization can be expressed as where the functions L j = L j (h) are defined in Appendix A 1.
We have also introduced a = a ·ĥ as the component of a parallel to the applied field. To find closed-form expression for the magnetization dynamics (13), we also need to express the second moment in terms of a, ee = [L 2 + (L 3 − L 1 L 2 )a ]ĥĥ as well as, for r = r(|he |), r e e −he = r ee e −he eq · a = R 1 a ĥ + R 2 a, where R 1 = (3I 2 − I 0 )/2, R 2 = (I 0 − I 2 )/2 with I n = r e n e −he eq = 1 2h n sinh(h) h −h y n r(y)dy. (C6) Plugging these expressions into Eq. (13) we arrive at with with relaxation times given by Eq. (14); they depend on the rate function r. Using Eq. (C6), we find for Arrhenius rates (r = 1) I 2 = R 2 = h/ [3 sinh(h)]. More generally, we use the ansatz r(x) = cosh(αx) with parameter 0 α 1. For α = 0 this reduces to Arrhenius rate r(x) = 1. For general α we find

APPENDIX D: NÉEL RELAXATION FROM DIFFUSION-JUMP MODEL
We here focus on Néel relaxation for large magnetic anisotropies and consider immobile MNPs, where Brownian rotations of the particles are suppressed. Within the diffusionjump model (12), only magnetization reversals e ↔ −e are possible in this case. Therefore, we use the ansatz for the probability density f (e; t ) = p + (t )δ(e − n) + p − (t )δ(e + n), where n denotes direction of the particle's easy axis. Normalization of the probability density requires p + (t ) + p − (t ) = 1. Note that we here assume all MNPs are frozen with identical easy-axis direction n. We denote with hn = h · n the projection of the dimensionless magnetic field onto the easy-axis direction. With the probability density (D1), the time-dependent magnetization can be expressed as and therefore d dt e = 2ṗ + n. On the other hand, the Néel contribution to the magnetization relaxation is obtained from Eq. (12) in the general form Evaluating the right-hand side of (D3) with (D1), we arrive at describing a single-exponential relaxation of the magnetization with the effective relaxation time strongly decreasing with increasing magnetic field strength.
In order to compare Eq. (D5) to Brown's classical result, we consider all MNPs easy axes to be perfectly aligned with the external field, n = 1, so that the magnetization component parallel to the field m = e = e · n obeys the same single-exponential relaxation, with τ eff N given by Eq. (D5) with n = 1.
Let us contrast Eq. (D5) with Brown's result, obtained for large anisotropy barriers [18] For sufficiently large anisotropies κ 1 and not too large fields h/κ 1, Brown's result (D6) simplifies to τ N (h) ≈ τ N / cosh(h), identical to Eq. (D5) for r = 1. Therefore, we conclude that the diffusion-jump model recovers the field dependence of the Néel relaxation in this limit, assuming Arrhenius activation law, which corresponds to r = 1.

APPENDIX E: REDUCED MOMENT APPROXIMATION (RMA)
The FP equation (3) for the egg model can be written in the generic form ∂F/∂t =MδS/δF with the Boltzmann-Gibbs-Shannon entropy S[F ] = − F ln(F/F eq )de dn and F eq the Boltzmann equilibrium distribution [59]. The operatorM is identified from the first bracket in the second line of Eq (3).
For the relaxation processes of interest here, we consider small deviations from equilibrium, F = F eq (1 + ), with | | 1. From the normalization of F we require eq = 0, where we have introduced the short notation for equilibrium averages • eq = • F eq de dn. With this ansatz, the entropy becomes to lowest order S ≈ −(1/2) 2 eq + · · · . Next, we make the assumption that can be represented as a linear combination = n j=1 a j (X j − X j eq ) ( E 1 ) with given functions X j = X j (e, n) and unknown coefficients a j (t ). By construction eq = 0 as required. With the ansatz 134433-14 (E1), the entropy can be expressed as a quadratic form with δ X j = X j − X j eq and the matrix B i j defined as The quadratic entropy (E2) corresponds to Einstein's fluctuation theory. With Eq. (E2) we can write the time evolution equations for the moments X i in vector notation as where the elements M i j = X iM X j de dn of the matrix M can be expressed as Note that Eq. (E5) is in a form familiar from irreversible thermodynamics with the friction matrix M and the entropy gradient B −1 · δ X [59]. Since we consider only small deviations from equilibrium, we are in the regime of linear irreversible thermodynamics where the matrix C = M · B −1 is independent of the variables X . In the following, our aim is to find the matrix C and its eigenvalues for the special choice (17) of the variables X.

Matrix coefficients
For N = 3 and the chosen set X 1 = e , X 2 = n (e · n), X 3 = n 2 , the coefficients of the matrix B in Eq. (E3) can be evaluated as where L j denotes the jth-order Langevin function (Appendix A 1), and S the th alignment order parameter (Appendix A 2). Evaluating the coefficients of the matrix M from Eq. (E5) for the current choice of X j leads to Note that the matrix B as well as M are manifestly symmetric, while C = M · B −1 is not symmetric. With B and M at hand, C and its eigenvalues can be calculated numerically for arbitrary values of the model parameters τ B , τ D , κ, and h.

APPENDIX F: RMA PERTUBATION THEORY
We want to find analytical expressions for the eigenvalues in the case of weak fields h 1. Therefore, we expand the matrix C = M · B −1 around h = 0. Perturbation theory aims at solving the corresponding eigenvalue problem considering an expansion of C, λ k , and v k in terms of h. Because we will find that first-order corrections to the eigenvalues are zero, we start with an expansion up to O(h 3 ): Since eigenvectors are determined only up to an overall constant, we choose the common normalization |v k | = 1. With the ansatz (F2), it follows that |v (0) k | = 1 and v (0) k · v (1) k = 0, i.e., the first-order correction should be perpendicular to the zeroth-order eigenvector. Using this ansatz in the eigenvalue equation (F1) and insisting that each order o ∈ {0, 1, 2} must be satisfied separately, the recursive equations determining the contribution λ (o) k to the eigenvalues λ k , as well as the corresponding vectors v (o) k read as The perturbation approach in its simplest form is only applicable for nondegenerate C (0) , as long as the eigenvalues of C (0) do not cross.

Zeroth order
In zeroth order o = 0, Eq. (F3) reduces to the eigenvalue problem in the absence of a magnetic field (C (0) − λ (0) k I) · v (0) k = 0. To this end we need to determine the matrix C h = 0 (see Appendix A 1), we obtain with The eigensystem of C (0) is where we have mentioned the left eigenvectors u (0) k of C (0) for later use [60]. The decoupling of the magnetization and easyaxis orientation in the absence of a field is evident from the orthogonality of the eigenvectors v (0) 1 and v (0) 3 . For the present case, Eqs. (F9)-(F11) imply Note that the h-independent eigenvalues of C (0) are degenerated under certain conditions In the limit of small κ 1 (cf. Appendix A 2), the zerothorder eigenvalues behave as while λ (0) 3 does not depend on κ. All three λ (0) k are thus identical for τ B = 2τ D and κ = 0 because S 2 = 0 at κ = 0. Depending on the values for τ B and τ D , there is in any case a κ value for which two of the eigenvalues coincide since (τ B − 2τ D )/(τ B + 4τ D ) ∈ [− 1 2 , 1] and S 2 ∈ [0, 1]. For example, if τ B /τ D = 3, the eigenvalues λ (0) 2 and λ (0) 3 are identical at κ ≈ 1, and if τ B /τ D = 1 2 , the eigenvalues λ (0) 1 and λ (0) 3 are identical at κ ≈ 5.

First order
To find the first-order effect of a magnetic field on the eigenvalues λ (1) k , we multiply Eq. (F3) for o = 1 by the left eigenvectors u (0) k of C (0) . This yields u (0) k · (C (1) − λ (1) k ) · v (0) k = 0 or, equivalently, Note that Eq. (F18) is known as first-order perturbation theory for the eigenvalues of nonsymmetric matrices [60]. Next, we use the ansatz B −1 =B 0 + hB (1) for the left inverse of the matrix B = B (0) + hB (1) Inserting these expressions into (F18), using the eigensystem of C (0) provided by Eqs. (F9)-(F11), gives Therefore, we find that there are no first-order contributions of the magnetic field to the eigenvalues. Using this latter result, the first-order version of Eq. (F3) simplifies to This equation will be used to calculate v (1) k required for the second-order corrections to λ k .

Second order
For the second-order contributions to the eigenvalues, we scalar multiply the second order o = 2 version of (F3) with u (0) k , use u (0) k · C (0) = λ (0) k C (0) as before, to find where we made use of λ (1) k = 0 according to Eq. (F28). Extending Eq. (F18), we can interpret Eq. (F30) as second-order perturbation theory for the eigenvalues of nonsymmetric matrices for the special case of vanishing first-order corrections. To evaluate this expression, we need the second-order matrix C (2) . Adapting the strategy we used to obtain C (1) , or by Taylor expanding the full expression for C up to h 2 , it is seen to be of the form The coefficients of C (2) where we recall that S = S (κ ) is the th alignment order parameter (Appendix A 2). The additional contribution in Eq. (F30) involves the firstorder correction to the eigenvectors v (1) k . They can be obtained from Eq. (F29), where the right-hand side reads as Requiring in addition that v (1) k · v (0) k = 0, there are as many unknowns as independent equations upon inserting Eq. (F36) into Eq. (F29). The solution reads as v (1) 1 = 0, 0, − with v (1) 3,x given by Note that all components of v (1) k required in the expression (F32)-(F34) for λ (2) k turn out to be nonzero, while all others vanish. Note also that the second-order expansion results are only applicable as long as the zeroth-order eigenvalues do not cross. The crossing conditions were stated in Eqs. (F13)-(F15).
If we now insert the v (1) k from Eqs. (F37)-(F39) into Eqs. (F32)-(F34), making use of the coefficients of the matrices C (1) and C (2) , Eqs. (F24)-(F27) and (F35), we end up with relatively long but explicit expressions for λ (2) k . The λ (2) k does not apply for τ B = 2τ D and κ = 0, and thus cannot be Taylor expanded about κ = 0, as discussed in Appendix F 1, except if we stay away from this region. Assuming τ B 2τ D , 134433-17 up to terms O(κ 2 ), again with the help of Appendix A 2. These second-order corrections to λ k = λ (0) k + h 2 λ (2) k serve to explain the qualitative behavior of the eigenvalues with respect to all parameters at small h and κ. All but λ (2) 3 at τ B τ D tend to be positive, in agreement with both the numerical RMA and the solution of the FP equation.