Comparative Modeling of Frequency Mixing Measurements of Magnetic Nanoparticles Using Micromagnetic Simulations and Langevin Theory

Dual frequency magnetic excitation of magnetic nanoparticles (MNP) enables enhanced biosensing applications. This was studied from an experimental and theoretical perspective: nonlinear sum-frequency components of MNP exposed to dual-frequency magnetic excitation were measured as a function of static magnetic offset field. The Langevin model in thermodynamic equilibrium was fitted to the experimental data to derive parameters of the lognormal core size distribution. These parameters were subsequently used as inputs for micromagnetic Monte-Carlo (MC)-simulations. From the hysteresis loops obtained from MC-simulations, sum-frequency components were numerically demodulated and compared with both experiment and Langevin model predictions. From the latter, we derived that approximately 90% of the frequency mixing magnetic response signal is generated by the largest 10% of MNP. We therefore suggest that small particles do not contribute to the frequency mixing signal, which is supported by MC-simulation results. Both theoretical approaches describe the experimental signal shapes well, but with notable differences between experiment and micromagnetic simulations. These deviations could result from Brownian relaxations which are, albeit experimentally inhibited, included in MC-simulation, or (yet unconsidered) cluster-effects of MNP, or inaccurately derived input for MC-simulations, because the largest particles dominate the experimental signal but concurrently do not fulfill the precondition of thermodynamic equilibrium required by Langevin theory.


Introduction
Magnetic nanoparticles (MNPs) have a plethora of applications not only in biomedical diagnostics, mainly determined by sample preparation, but also in detection [1,2]. MNPs are used as markers for biomolecules in immunoassays; in addition to the well-established techniques of AC-Susceptometry [3] and Relaxometry [4], Frequency-Mixing Magnetic Detection (FMMD) [5] has been increasingly applied during the past decade because of its high selectivity. This technique relies on simultaneously applying a low-frequency magnetic driving field, which brings the particles close to saturation, and a high-frequency excitation field, which probes the particles' susceptibility. Due to the nonlinear magnetization of the MNP, harmonics of both the individual incident frequencies and the intermodulation products of both frequencies are generated. Their signal can usually be picked up using a receive coil; however, other magnetic detectors can also be employed [6]. Due to its high sensitivity, FMMD has been successfully applied to realize magnetic immunoassays for the detection of a multitude of different analytes, for instance, viruses [7], antibiotics [8], or bacterial toxins [9]. Furthermore, FMMD enables the distinction of different types of MNP based on their different frequency mixing response spectrum [10][11][12], which opens the potential for multi-analyte detection.
In this work, we present offset-field-dependent FMMD measurements of MNPs and compare them quantitatively with two different modeling approaches: a simple Langevinfunction-based thermal equilibrium model based on a lognormal size distribution, and a micromagnetic Monte Carlo (MC)-simulation method [13]. Although particles are immobilized in the experiment and can thus only relax according to Néel relaxation, the MCsimulation also includes Brownian relaxation. With this, micromagnetic MC-simulations provide new insight into the frequency mixing excitation of MNPs by modeling their nonequilibrium dynamics. For the Langevin model, the nature of relaxation is irrelevant. The core size distribution parameters derived from the Langevin model were used as starting values for the MC-simulation. As a consistency check, the MC results were again compared against the Langevin model. Prospects and limitations of both models are discussed, and suggestions for future research are developed.

Experimental Setup for Frequency Mixing Magnetic Detection
A custom-built measurement setup was employed for simultaneous sample excitation and demodulation signal detection. This comprises two excitation and two pick-up coils in a measurement head unit, allowing digital frequency demodulation directly from the detected signal. Details on the setup can be found in earlier works [5,14,15]. The applied alternating magnetic excitation field (AMF) was: H(t) = H 0 + H 1 sin(2π f 1 t) + H 2 sin(2π f 2 t) (1) where H 0 denotes the static magnetic offset field, H 1 = 1.29 mT/µ 0 is the magnetic field amplitude at high frequency f 1 = 30, 543 Hz, and H 2 = 16.4 mT/µ 0 is the amplitude at low frequency f 2 = 62.95 Hz. With this low-frequency amplitude, the particles are driven well into the nonlinear regime of the magnetization curve. The almost 400-fold higher f 1 yields well-detectable voltages at the mixing frequencies. ABICAP ® columns from Senova GmbH (Weimar, Germany) containing polyethylene filters with pore sizes of approximately 50 µm, a height of 5 mm, and a diameter of 5 mm were primed with ethanol in a desiccator to remove air bubbles from the filter. After washing the column twice with 500 µL of distilled water, 450 µL of nanomag ® -D SPIO (Prod.#: 79-00-201; Micromod Partikeltechnologie GmbH, Rostock, Germany) was flushed through the column, followed by another washing step to ensure homogeneous MNP distribution and to remove unbound MNPs. The sample was dried at ambient conditions. Then the columns with the immobilized MNPs were measured at varying static field strength H 0 = (0, . . . , 24) mT/µ 0 in steps of 0.48 mT/µ 0 . Due to coil-current resistive heating, the temperature T in the measurement head was approximately 318 K. The first four nonlinear magnetic moment demodulation components from frequency mixing (both real and imaginary part) were stored: f 1 + n · f 2 with n = 1, 2, 3, 4. Background subtraction was performed using data from reference measurements without a sample, and phase correction for frequency-dependent phase shift inside the induction coil and amplification chain was performed.

Thermodynamic Langevin Model of a Magnetic Nanoparticle Ensemble
In the classical thermodynamic model description, the MNP sample can be described as an ensemble of noninteracting particles. (The validity of this assumption for our sample is assessed in Appendix A by an estimation of the dipole-dipole energy.) Neglecting surface effects, the saturation magnetic moment of a spherical particle with a core diameter d c is given by: with M s denoting the saturation magnetization of the MNP. For simplicity, all particles are assumed to be spheres. The total magnetic moment of the particle ensemble in thermodynamic equilibrium is calculated by averaging over a Boltzmann distribution of the orientations of the individual particle moments, yielding a dependence on the amplitude of the applied magnetic field |H| = H, which is governed by a Langevin function [16]: with the dimensionless magnetic field variable: and with temperature T, Boltzmann's constant k B = 1.38 ×·10 -23 J/K, and the permeability of vacuum µ 0 = 4π ×·10 -7 Vs/Am [17].
In the average over the particle ensemble, it has to be considered that the saturation magnetic moment m p of each particle depends on its diameter d c . Usually, particle ensembles exhibit a lognormal size distribution with a probability density function PDF(d c ) given by [17]: with the median diameter d 0 and the standard deviation σ of the diameters' natural logarithm.
The total magnetic moment of the ensemble of N p particles is then calculated by integrating over the lognormal distribution: In our FMMD scheme [5], the particle ensemble is exposed to a two-frequency excitation with static offset magnetic field (see Equation (1)). Inserting Equations (1) and (2) into (6) yields: As shown in [5], the nonlinearity of the magnetization curve gives rise to the emergence of intermodulation products m·f 1 ± n·f 2 of the total magnetic moment of the particle ensemble (with m and n denoting integers). In particular, the frequency mixing components f 1 + f 2 , f 1 + 2·f 2 , f 1 + 3·f 2 and f 1 + 4·f 2 appear. In the limit of small excitation amplitudes H 1 and H 2 , these frequency mixing responses can be calculated with a Taylor expansion of Equation (3), yielding offset field dependencies of the mixing components proportional to the higher order derivatives of the Langevin function [5]. In the case of larger excitation amplitudes H 1 and H 2 , the Taylor approximation is no longer valid and has to be replaced by the respective Fourier components of Equation (7). For instance, the average nonlinear moment response m 1 (d c ) at frequency f 1 + f 2 of one particle with diameter d c is given by: The factor 2/k normalizes the sum and accounts for the full-period average over sin 2 (..) which is 1 2 . The sampling time steps t i should be chosen such that a sufficient number of samples is taken in one high frequency period. In our calculations, we took 10 steps in a period 1/f 1 , i.e., ∆t = t i -t i-1 = 0.1/f 1 , yielding sufficient numerical precision. Although in our experiments, the high frequency f 1 was 485 times larger than the low frequency f 2 , it was sufficient to select f 1 = 20·f 2 for our numerical calculations. Thus, k = 200 was used.
In a similar fashion, the response component m 2 (d c ) at frequency f 1 + 2·f 2 is obtained from: Note the cos[..]/sin[..] alternation in the reference frequency term behind the sum symbol in Equations (8) and (9), which is due to the fact that with increasing order, the frequency mixing responses are alternately uneven (point-symmetric) and even (axisymmetric) functions. Components m 3 (d c ) and m 4 (d c ) are calculated similarly.
The total magnetic moment component m n,tot at the frequency mixing component f 1 + n·f 2 is then obtained by integration over the lognormally weighted particle ensemble: Equation (10) constitutes our forward model for calculating the frequency mixing signals. The model contains just three fitting parameters, the lognormal distribution characteristics median diameter d 0 and width σ, and the total number of particles N p . The measured nonlinear magnetic moment components of nanoparticle samples at frequencies f 1 + n·f 2 , n = 1, 2, 3, 4, were fitted with this model using the Levenberg-Marquardt least-squares algorithm.

Micromagnetic Monte Carlo (MC-)Simulation
The nonlinear particle relaxation dynamics in nonequilibrium conditions under the influence of an applied AMF, H, can be described by combined Néel-Brownian relaxation [13]. The Néel relaxation of the direction of the magnetic moment of a single MNP, m p , is governed by the Landau-Lifshitz-Gilbert equation (LLG) [18]: with the permeability of free space, µ 0 , the electron gyromagnetic ratio, γ, the damping parameter, α, and the effective field H eff . The Brownian rotation of the MNP easy axis, n, is described via the generalized torque, Θ, as follows [19]: with the carrier matrix viscosity, η, and the MNP hydrodynamic volume, V h = π/6 · d 3 h , in which d h is the hydrodynamic particle diameter. Néel and Brownian relaxation are coupled in the internal particle energy: where m p = m p = V c · M S (cf. Equation (2)) is the magnitude of the MNP magnetic moment, and V c = π/6 · d 3 c is the MNP core volume. The first term describes the Zeeman energy including the applied field H. The second term incorporates the anisotropy energy via the effective anisotropy constant, K, and under the assumption of uniaxial anisotropy and spherically shaped particles. The third term in Equation (13) includes magnetic dipole-dipole interaction. However, thermal energy, ε therm , dominates magnetic interaction energy by two orders of magnitude in our MNP samples, so that particleparticle interaction is negligible (ε I A ε therm ). Please refer to Appendix A for a detailed estimation of the effect of magnetic dipole-dipole interaction energy that corroborates our assumption. Thermal fluctuations are included by adding the terms H th and Θ th , which are implemented as Gaussian-distributed white noise with zero mean values ( H i th (t) = 0 and Θ i th (t) = 0) and variances H i Here T is the global temperature in the system. With this, the effective field and generalized torque read: We apply the Stratonovic-Heun scheme to solve the system of coupled stochastic differential Equations (11) through (15) and implement a Monte Carlo method routine as described in our previous works [20][21][22]. The full source code is available as listed in the Data Availability Statement and its results are denoted as MC-simulations henceforth. Simulation parameters were chosen as listed in Table 1 with the MNP properties matching the experimentally determined values from Table 2. Furthermore, the damping parameter α was set to unity [23]. One thousand particles were simulated simultaneously and initialized with randomized directions of magnetization and easy axes for each MNP. The MNPs were then thermalized for one-fifth of the total number of time steps, N, before the AMF was applied. We used N = 50,000 and averaged the magnetization over five independent simulation runs to achieve a good compromise between accuracy and acceptable computation time. The time step sizes were then 10 ns. The simulations were performed with the open-access Python code referenced in the Data Availability Statement. Calculations were carried out on a PC cluster consisting of 2 × CPU Intel Xeon E5-2687W, 3.1/3.8 GHz, with 8 clusters each and RAM 64 GB each. The typical calculation time for one offset field value was approx. 53 h.
To approximate experimental data for comparison, the excitation field parameters were chosen as H 1 = 1 mT/µ 0 , f 1 = 40,000 Hz, H 2 = 16 mT/µ 0 , f 2 = 2000 Hz, and the static field was varied from H 0 = (0, . . . , 24) mT/µ 0 with a step size of 1 mT/µ 0 .   1 The concentration c is taken from the datasheet of the manufacturer. The concentration in the filter might be smaller due to unbound particles being washed out undetected. 2 The hydrodynamic diameter d h is taken from the datasheet of the manufacturer.

Experimental Results and Thermodynamic Langevin Model Fitting
The Langevin model (Equation (10), with inputs (8) and (9), see Section 2.2) was fitted to the measured real parts of the experimental data using the Levenberg-Marquardt leastsquare algorithm routine, whose results are plotted in Figure 1. Although the immobilized MNPs in the experimental setup are blocked in Brownian rotation, this step is justified because Langevin theory approximates the magnetization of the entire ensemble of MNP independently of the underlying mechanism of relaxation. The imaginary parts of the response signal were found to be two orders of magnitude weaker and were therefore disregarded. For all four demodulation components f 1 + n · f 2 with n = 1, 2, 3, 4, a very good agreement between experimental and simulated results was observed, confirmed by a coefficient of determination of R 2 > 0.98. Only for component f 1 + f 2 at high offset field did the simulation predict slightly higher values than measurement, whereas for f 1 + 2·f 2 and f 1 + 3·f 2 , the simulation slightly underestimated measurements.  1 The concentration c is taken from the datasheet of the manufacturer. The concentration in the filter might be smaller due to unbound particles being washed out undetected. 2 The hydrodynamic diameter dh is taken from the datasheet of the manufacturer.

Experimental Results and Thermodynamic Langevin Model Fitting
The Langevin model (Equation (10), with inputs (8) and (9), see Section 2.2) was fitted to the measured real parts of the experimental data using the Levenberg-Marquardt leastsquare algorithm routine, whose results are plotted in Figure 1. Although the immobilized MNPs in the experimental setup are blocked in Brownian rotation, this step is justified because Langevin theory approximates the magnetization of the entire ensemble of MNP independently of the underlying mechanism of relaxation. The imaginary parts of the response signal were found to be two orders of magnitude weaker and were therefore disregarded. For all four demodulation components + ⋅ with = 1,2,3,4, a very good agreement between experimental and simulated results was observed, confirmed by a coefficient of determination of > 0.98. Only for component f1 + f2 at high offset field did the simulation predict slightly higher values than measurement, whereas for f1 + 2⋅f2 and f1 + 3⋅f2, the simulation slightly underestimated measurements. Fitting yielded the MNP material properties of median core diameter, d c , and its log-normal distribution width, σ, which are listed in Table 2. The hydrodynamic diameter, d h , and the concentration of the MNP were taken from the datasheet of the manufacturer.

Micromagnetic MC-Simulation Results
We used the MNP material properties derived from fitting the Langevin model to the experimental data as described in the previous Section 3.1 (see Table 2) as input parameters for the micromagnetic MC-simulations. These simulations yield magnetization curves (M(H)-loops), which are shown for exemplary static offset fields in Figure 2.
Fitting yielded the MNP material properties of median core diameter, , and its log normal distribution width, , which are listed in Table 2. The hydrodynamic diamete , and the concentration of the MNP were taken from the datasheet of the manufacture

Micromagnetic MC-Simulation Results
We used the MNP material properties derived from fitting the Langevin model to th experimental data as described in the previous section (3.1, see Table 2) as input param ters for the micromagnetic MC-simulations. These simulations yield magnetization curve (M(H)-loops), which are shown for exemplary static offset fields in Figure 2. The magnetization curves are (almost) identically overlapping, but shift towards th magnetization saturation of MNP as the static offset field, , increases. The M(H)-loop are very slightly opened, revealing minor hysteresis area (see Figure 2, inset). For value of > ( + ), e.g., for = 20 mT/µ0 in Figure 2, the applied field is constant keeping the MNPs in (almost) saturation and no hysteresis area is present.
From the M(H)-loops, the first four demodulation components + ⋅ were ca culated following Equations (8) and (9). For comparison, the experimental measuremen data and the MC-simulated data were normalized to their respective highest value, and plotted alongside each other in Figure 3. The agreement between experimental an simulated data is acceptable, as confirmed by a coefficient of determination of > 0.7 for all four demodulation components. We observe for component f1 + f2 that MC-simul tion consistently underestimates the measurement results. However, the peak at = 1 mT/µ0 coincides. For the mixing term f1 + 2⋅f2, MC-simulation underestimates the mea urement results for  The magnetization curves are (almost) identically overlapping, but shift towards the magnetization saturation of MNP as the static offset field, H 0 , increases. The M(H)-loops are very slightly opened, revealing minor hysteresis area (see Figure 2, inset). For values of H 0 > (H 1 + H 2 ), e.g., for H 0 = 20 mT/µ 0 in Figure 2, the applied field is constantly keeping the MNPs in (almost) saturation and no hysteresis area is present.
From the M(H)-loops, the first four demodulation components f 1 + n · f 2 were calculated following Equations (8) and (9). For comparison, the experimental measurement data and the MC-simulated data were normalized to their respective highest value, M max , and plotted alongside each other in Figure 3. The agreement between experimental and simulated data is acceptable, as confirmed by a coefficient of determination of R 2 > 0.76 for all four demodulation components. We observe for component f 1 + f 2 that MC-simulation consistently underestimates the measurement results. However, the peak at H 0 = 15 mT/µ 0 coincides. For the mixing term f 1 + 2·f 2 , MC-simulation underestimates the measurement results for H 0 < 18 mT/µ 0 and H 0 > 20 mT/µ 0 and shows a mismatch of over 50% for the peak value at H 0 = 15 mT/µ 0 . MC-simulation and measurement of component f 1 + 3·f 2 coincide for values up to H 0 = 10 mT/µ 0 and again for H 0 ≥ 20 mT/µ 0 . However, around the peak at H 0 = 16 mT/µ 0 , MC-simulations overestimate experimental data by approx. 70%. For mixing term f 1 + 2·f 2 , MC-simulations match experimental data, except around the peaks at H 0 =11 mT/µ 0 and H 0 = 16 mT/µ 0 .

Comparing Micromagnetic MC-Simulation Results and Thermodynamic Langevin Model Fitting
To test whether predictions from the Langevin model and MC-simulations show reproducible results, we fitted the Langevin model directly to the results from MC-simulation. We used the same input parameters for MC-simulations and Langevin model fitting of = 1 mT/µ0, = 40,000 Hz, = 16 mT/µ0, = 2,000 Hz, and fixed the mean core diameter with = 7.81 nm and variable distribution parameter for the fitting. The results are plotted in Figure 4, confirming overall good agreement with a coefficient of determination of > 0.989 for all four demodulation components. From the qualitative comparison, one sees that Langevin model fitting and MC-simulations coincide, except for the secondary peaks for = 3 and = 4 (cf. Figure 4). The fitting yields a distribution width of = 1.466. This is significantly different than the input parameters to MC-simulations ( = 0.346), questioning our hypothesis that assumes identical outputs from the Langevin model and MC-simulation. As we discuss in detail in the next section, we suspect the reasons for this lie with MNP properties that are not (yet) accurately represented in the modeling.

Comparing Micromagnetic MC-Simulation Results and Thermodynamic Langevin Model Fitting
To test whether predictions from the Langevin model and MC-simulations show reproducible results, we fitted the Langevin model directly to the results from MC-simulation. We used the same input parameters for MC-simulations and Langevin model fitting of H 1 = 1 mT/µ 0 , f 1 = 40, 000 Hz, H 2 = 16 mT/µ 0 , f 2 = 2000 Hz, and fixed the mean core diameter with d c = 7.81 nm and variable distribution parameter σ for the fitting. The results are plotted in Figure 4, confirming overall good agreement with a coefficient of determination of R 2 > 0.989 for all four demodulation components. From the qualitative comparison, one sees that Langevin model fitting and MC-simulations coincide, except for the secondary peaks for n = 3 and n = 4 (cf. Figure 4). The fitting yields a distribution width of σ = 1.466. This is significantly different than the input parameters to MC-simulations (σ = 0.346), questioning our hypothesis that assumes identical outputs from the Langevin model and MC-simulation. As we discuss in detail in the next section, we suspect the reasons for this lie with MNP properties that are not (yet) accurately represented in the modeling.

Discussion
The Langevin model has been extensively applied to describe MNP magnetization for various applications such as diagnostic biosensing [25], Magnetic Particle Imaging (MPI) [26], and therapeutic Magnetic Fluid Hyperthermia (MFH) [27]. Its application to frequency mixing excitation seems therefore naturally reasonable. This assumption is fos tered by our results reporting very good agreement between experimental data and Langevin model fitting ( > 0.98; cf. Figure 1). The resulting fitting parameters ( ≈ 7.8 nm and ≈ 0.35; cf. Table 2) are furthermore in good agreement with literature values reporting a core diameter between ≈ 7 nm [28] and ≈ 11 nm [29] for nanomag ® -D SPIO. The broad size distribution (with width > 0.3) reflects a heterogeneous particle size with an effectively range of = (3 − 25) nm, cf. Figure 5.
Nevertheless, the Langevin model is only valid for noninteracting, single particle en sembles with uniaxial anisotropy in thermodynamic equilibrium and (quasi)static mag netic fields [16]. Therein lies its major limitation, because the Langevin model canno model an opening of M(H)-loops (hysteresis), which are, however, expected for AMF ex citations at applied frequencies that approximate the inverse MNP relaxation time  Figure 2), which is inaccessible via equilibrium Langevin theory. We attribute this minor hysteresis to the small portion of large particles, whose magnetic relaxation is (thermally) blocked by the volume (and therefore size cubed) dependent anisotropy bar rier: ⋅ . This can be assessed by assuming an AC-field-dependent Brownian relaxation time of [30]:

Discussion
The Langevin model has been extensively applied to describe MNP magnetization for various applications such as diagnostic biosensing [25], Magnetic Particle Imaging (MPI) [26], and therapeutic Magnetic Fluid Hyperthermia (MFH) [27]. Its application to frequency mixing excitation seems therefore naturally reasonable. This assumption is fostered by our results reporting very good agreement between experimental data and Langevin model fitting (R 2 > 0.98; cf. Figure 1). The resulting fitting parameters (d c ≈ 7.8 nm and σ ≈ 0.35; cf. Table 2) are furthermore in good agreement with literature values reporting a core diameter between d c ≈ 7 nm [28] and d c ≈ 11 nm [29] for nanomag ® -D SPIO. The broad size distribution (with width σ > 0.3) reflects a heterogeneous particle size with an effectively range of d c = (3 − 25) nm, cf. Figure 5.
Nevertheless, the Langevin model is only valid for noninteracting, single particle ensembles with uniaxial anisotropy in thermodynamic equilibrium and (quasi)static magnetic fields [16]. Therein lies its major limitation, because the Langevin model cannot model an opening of M(H)-loops (hysteresis), which are, however, expected for AMF excitations at applied frequencies that approximate the inverse MNP relaxation time, f ∼ τ −1 , as the magnetic moment of each MNP begins to lag behind the driving AMF [2]. In contrast, micromagnetic MC-simulations provide new insight in the frequency mixing excitation of MNP by modeling the non-equilibrium dynamics of MNP by including thermal fluctuations: The MC-simulations reveal a minor hysteresis in the M(H)-loop of the MNP (cf. Figure 2), which is inaccessible via equilibrium Langevin theory. We attribute this minor hysteresis to the small portion of large particles, whose magnetic relaxation is (thermally) blocked by the volume (and therefore size cubed) dependent anisotropy barrier: K · V C . This can be assessed by assuming an AC-field-dependent Brownian relaxation time of [30]: (16) and an AC-field dependent Néel relaxation time of [16]: where H AC denotes the AMF amplitude, τ B = 3ηV H /(k B T) the field-independent Brownian relaxation time, ξ(H AC ) = mH AC /(k B T) the reduced field parameter (cf. Equation (4)), τ 0 = 10 −9 s the time constant, and H K = 2K/µ 0 M S the (uniaxial) anisotropy field strength. From Equations (16) and (17), the effective relaxation time follows with: Using the (mean) values from Tables 1 and 2, and the experimental setup parameters with H AC = H 1 + H 2 ∼ 17 mT/µ 0 , one calculates τ ∼ τ 0 ∼ 10 −9 s so that for f ∼ (30 − 40) kHz the condition for the onset of hysteresis, f ∼ τ −1 , is not fulfilled. Because the measured imaginary parts of the mixing components are hundredfold weaker than the real parts, we can conclude that dissipation is indeed negligible. However, for the small portion of large particles with d C ≥ 20 nm and naturally increasing the effective hydrodynamic diameter d h > 20 nm, and presumably decreasing effective anisotropy constant for these larger core particles, K ∼ 5 kJ/m 3 [31,32], the effective relaxation time increases by several orders of magnitude. The relaxation time for these large particles is dominated by the Brownian relaxation mechanism, fulfilling the condition f ∼ τ −1 . Consequently, the minor opening of the M(H)-loop observed in Figure 2 from MC-simulations is a direct contribution from these large size particles relaxing with Brownian rotation. This assumption could, however, not be verified experimentally in our current study, because the MNPs were immobilized due to sample preparation, blocking Brownian rotation. Nevertheless, Figures 2 and 4 demonstrate the potential insight to be gained on the micromagnetic level in the underlying mechanisms in dual frequency excitation responses of MNP (see discussion regarding future experiments below).
Interestingly, these large size particles apparently also contribute most dominantly to the Langevin model fitting to MC-simulation data (cf. Figure 4): When we deliberately restrict the particle size range for Langevin model calculation to the largest portion of core sizes (i.e., limiting the minimum core size d c ), we find that 90% of the calculated signal is contributed from the particles with d c ≥ 12.1 nm. Similarly, 99% of the signal stems from d c ≥ 9.2 nm. With the reverse of the cumulative distribution function (CDF) of the lognormal distribution (see Equation (A3) in Appendix B), the corresponding size fractions (quantiles) of the distribution are calculated, beginning from the large-sized tail. This finding is visualized in Figure 5. The largest 10.3% of the particles contribute 90% of the FMMD signal. Furthermore, 99% of the signal is produced by the largest 31.8% of particles. In other words, this means that almost all of the FMMD signal originates from the large-sized tail of the size distribution.
This could also explain why the Langevin model can be fitted with very good agreement to MC-simulation data (cf. Figure 4; R > 0.989), even though it is physically unable to model nonequilibrium dynamic relaxation processes: The mathematical fitting routine can force the fitting parameters (d 0 , σ, N p ) to values outside the model's range of validity (for further details see Appendix B). As the opening M(H)-loop and our relaxation time approximations (see above) suggest, thermodynamic equilibrium is not valid anymore for large particles d C ≥ 20 nm at frequencies f ∼ (30 − 40) kHz. Therefore, fitting the model to the experimental data could lead to unreliable distribution parameters. This could also explain why the Langevin model can be fitted with very good agreement to MC-simulation data (cf. Figure 4; > 0.989), even though it is physically unable to model nonequilibrium dynamic relaxation processes: The mathematical fitting routine can force the fitting parameters ( , , ) to values outside the model's range of validity (for further details see Appendix B). As the opening M(H)-loop and our relaxation time approximations (see above) suggest, thermodynamic equilibrium is not valid anymore for large particles ≥ 20 nm at frequencies ~ (30 − 40) kHz. Therefore, fitting the model to the experimental data could lead to unreliable distribution parameters.
For other frequency-dependent biomedical applications of iron-oxide MNPs, such as MPI [33,34] or MFH [32,35], an optimal core size of ~ 25 nm has been suggested from experiment and MC-simulation, stressing the importance of open hysteresis loops for signal generation. Our results from Figures 4 and 5 could therefore indicate that an opening of the M(H)-loop is also favorable for generating strong signals in mixed frequency measurements. As mentioned above, our current sample preparation does not allow experimental verification of this assumption, because Brownian rotation is blocked for the immobilized MNP. Even more so, this assumption remains to be tested in future experiments using MNP with larger core sizes, ideally suspended and freely rotatable in solution. In future experiments, the question of whether the Brownian mechanism dominates the MNP relaxation mechanism could be further probed by suspending MNP in gelated matrices (e.g., agarose or poly-acrylamid gels), in which Brownian relaxation has been shown to be controllably blocked in MFH measurements with particles of the same size as used in the present study [36].
MC-simulations deviated quantitatively in intensity at the peaks by up to 70% from the experimental values (cf. Figure 3) when directly compared. The first reason for this could be that Brownian relaxation mechanism is blocked for immobilized MNPs in the experimental sample, while it adds to the signal in MC-simulations (see paragraph above). A second potential reason for deviations is that the parameters used for MC-simulation input slightly differ from the experimental parameters. Another reason could be that MNPs form clusters of a few particles in the experimental setting, which Dennis et al. For other frequency-dependent biomedical applications of iron-oxide MNPs, such as MPI [33,34] or MFH [32,35], an optimal core size of d c ∼ 25 nm has been suggested from experiment and MC-simulation, stressing the importance of open hysteresis loops for signal generation. Our results from Figures 4 and 5 could therefore indicate that an opening of the M(H)-loop is also favorable for generating strong signals in mixed frequency measurements. As mentioned above, our current sample preparation does not allow experimental verification of this assumption, because Brownian rotation is blocked for the immobilized MNP. Even more so, this assumption remains to be tested in future experiments using MNP with larger core sizes, ideally suspended and freely rotatable in solution. In future experiments, the question of whether the Brownian mechanism dominates the MNP relaxation mechanism could be further probed by suspending MNP in gelated matrices (e.g., agarose or poly-acrylamid gels), in which Brownian relaxation has been shown to be controllably blocked in MFH measurements with particles of the same size as used in the present study [36].
MC-simulations deviated quantitatively in intensity at the peaks by up to 70% from the experimental values (cf. Figure 3) when directly compared. The first reason for this could be that Brownian relaxation mechanism is blocked for immobilized MNPs in the experimental sample, while it adds to the signal in MC-simulations (see paragraph above). A second potential reason for deviations is that the parameters used for MC-simulation input slightly differ from the experimental parameters. Another reason could be that MNPs form clusters of a few particles in the experimental setting, which Dennis et al. consistently found for nanomag-D SPIO from small-angle neutron scattering (SANS) and transmission electron microscopy [37]. Clustering of MNPs has been suggested to lead to increased particle-particle interaction, changed effective anisotropy, and restriction of MNP rotatability [38][39][40]. Although the interplay of these clustering effects on MNP relaxation behavior is the subject of ongoing discussions [39][40][41][42], it is overall assumed to diminish the experimental signal intensity [43]. Both the rise of particle-particle interaction-although purposefully excluded from this study-and the influence of varying K could be investigated in future MC-simulations studies (e.g., by diffusion-limited colloidal aggregation [42]) to advance the understanding of its influence on mixed frequency excitation signal generation.
Finally, we are aware that our current study is limited by the extraction of input parameters for the MC-simulations from Langevin model fitting to experimental data (cf. Table 2). Future studies will complement these assumptions by experimental methods, e.g., measuring MNP core sizes via transmission electron microscopy (TEM) or the MNP magnetic core size distribution from magnetic measurement [44], in addition to determining the hydrodynamic diameters from dynamic light scattering (DLS) experiments [20], which will serve as input parameters to further validate our fitting and modelling.

Conclusions
With the current work, we present the first insights into dual-frequency magnetic excitation of MNP from interpreting experimental results in terms of both the thermodynamic Langevin model and nonequilibrium dynamic MC-simulations. In summary, we draw the following conclusions from our comparative study:

1.
Both the Langevin model and MC-simulations reproduce the shape of the experimental signal satisfactorily (see Figures 1 and 3). However, punctual deviations between experimental data and MC-simulations are observed (see Figure 3).

2.
MC-simulations allow the investigation of the dynamic hysteresis (M(H))-loop during AMF excitation, revealing a minor opening (cf. Figure 2). This opening is attributed to the small portion of large, thermally blocked particles. 3.
Langevin model fitting suggests that 90% of the experimentally detected FMMD signal intensity is generated by the largest 10% of the particles (cf. Figure 5).

4.
For the large particles (d c > 20 nm) which dominate the FMMD signal, relaxation cannot be neglected. However, this effect is not included in our Langevin model. We suspect this is the reason for observed deviations between the Langevin model and MC-simulations.
We are aware that these findings raise more questions and could serve as the starting point for further investigations. Based on these conclusions we propose the following respective future research directions:

•
Complementary experimental methods should be used to derive MNP properties for more accurate input in the MC-simulations (to address and remedy point 1, above). From this, we also plan to further increase our predictive accuracy in the future by coordinating simulation and experimental parameters to be identical (both materials properties and field parameters). • Experimental sample preparation should be expanded to allow-ideally gradually controllable-Brownian rotation of MNPs (e.g., in poly-acrylamide gels) in order to precisely analyze the role of the Brownian relaxation mechanism for signal generation (see also next point below). • Furthermore, MC-simulations could be expanded by including magnetic particleparticle interaction effects and the inclusion of clustering-effects of MNPs. This should be investigated with special regard to point 3 (above), because magnetic interaction scales with increasing particle core size. Additionally, the influence of effective anisotropy could be studied systematically with MC-simulations. Finally, the dominant relaxation mechanism in dual frequency excitation could be further investigated by weighting Néel and Brownian relaxation mechanisms systematically in MC-simulations, and comparing the results to experimental findings. • Point 4 (above) suggests the necessity for multi-theory approaches to interpret dualfrequency MNP excitation responses. Therefore, systematic parameter studies, e.g., varying MNP properties such as effective anisotropy, median core sizes and shell sizes, and AMF parameters (H 1 H 2 , f 1 , f 2 ), should be performed to generate a multiparameter repository from MC-simulations. This could serve as a basis for a unified model to explain FMMD signal generation in the future. However, such a study must be well designed to ensure acceptable computational effort.  From the intersection of and the corridor, the range of interparticle distances in which magnetic dipole-dipole interaction energy is of the same order of magnitude as the thermal excitation energy can be estimated with = (7.1 − 10.2) nm (see mark in Figure A1). In other words, for average interparticle distances > 10.2 nm, thermal ex- Figure A1. Estimated corridor for magnetic dipole-dipole interaction energy between neighboring MNPs (shaded) versus the average interparticle distance, estimated from Equation (A2). Thermal excitation energy is shown for reference.
From the intersection of ε therm and the corridor, the range of interparticle distances in which magnetic dipole-dipole interaction energy is of the same order of magnitude as the thermal excitation energy can be estimated with d = (7.1 − 10.2) nm (see mark in Figure A1). In other words, for average interparticle distances d > 10.2 nm, thermal excitation can be assumed to dominate magnetic interaction.
The average interparticle distance d avg can be estimated from the number of particles in solution, z, by taking the inverse third root: d avg = z −1/3 . With z = 4.3 · 10 15 1 cm 3 (taken from nanomag-D SPIO manufacturer's datasheet) we estimate an average interparticle distance of d avg ≈ 61.5 nm. From Figure A1 one can see that for such an average interparticle distance the thermal excitation energy is two orders of magnitude larger than magnetic dipole-dipole interaction energy, ε therm ε I A , corroborating our assumption to neglect magnetic interaction between MNP for our calculations.

Appendix B. Log-Normal Distribution Probability Density Function (PDF)
The log-normal distribution probability density function (PDF) for a general particle size d is defined as shown in Equation (5) with the median diameter d 0 and the standard deviation σ of the diameters' natural logarithm. Integration over PDF from 0 to d c yields the cumulative distribution function (CDF) [17]: Standard CDF cumulates particles starting from the small-sized tail of the distribution. For FMMD, large particles contribute most of the signal. To quantify their contribution, it is advantageous to consider the reverse cumulative distribution starting from the large size tail, i.e., 1-CDF. The exemplary log-normal distribution PDF(d c ,d 0 ,σ) and its reverse CDF using the parameters from Table 2 are plotted in Figure 5.
By fitting the Langevin model with fixed d c = 7.81 nm to the MC-simulation results, a significantly larger width parameter σ = 1.466 was obtained (see Section 3.3). The corresponding lognormal distribution and its reverse CDF are plotted in Figure A2. The FMMD signal calculated with this distribution is plotted in Figure 4. Analogously to the quantile analysis described in Section 4, the quantiles yielding 90%, 99%, and 99.9% of the FMMD signal were also determined for this wider distribution, as shown in Figure A2. analysis described in Section 4, the quantiles yielding 90%, 99%, and 99.9% of the FMMD signal were also determined for this wider distribution, as shown in Figure A2.  For this wider distribution, it was found that the largest 3.0% of the particles (with d c > 18.8 nm) contribute 90% of the FMMD signal. The largest 11.4% of particles (with d c > 13.7 nm) contribute 99% of signal. Almost all of the signal (99.9%) is generated by the largest 24.3% of the particles above 10.8 nm core size.