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

Abstract

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.

1. 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 Langevin-function-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 MC-simulation also includes Brownian relaxation. With this, micromagnetic MC-simulations provide new insight into the frequency mixing excitation of MNPs by modeling their non-equilibrium 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.

2. Materials and Methods

2.1. 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\left(t\right)=H_0+H_1\sin \left(2\pi f_{1}t\right)+H_2\sin \left(2\pi f_{2}t\right)$$

(1)

where $H_0$ denotes the static magnetic offset field, $H_1=1.29$ mT/µ₀ is the magnetic field amplitude at high frequency $f_1=30,543 \mathrm{Hz}$, and $H_2=16.4$ mT/µ₀ 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₁ 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=\left(0,\dots ,24\right)$ mT/µ₀ in steps of 0.48 mT/µ₀. 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\cdot 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.

2.2. 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:

$$\left|{\mathbf{m}}_p\right|=m_p\left(d_c\right)=\frac{M_s\pi d_c^3}{6}$$

(2)

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 $\left|\mathbf{H}\right|=H$, which is governed by a Langevin function [16]:

$$ℒ \left(\xi \right)=\coth \left(\xi \right)-\frac{1}{\xi }$$

(3)

with the dimensionless magnetic field variable:

$$\xi =\frac{m_p{\mu }_{0}H}{k_{B}T}$$

(4)

and with temperature T, Boltzmann’s constant k_B = 1.38 ×⋅10^–23 J/K, and the permeability of vacuum µ₀ = 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]:

$$\mathrm{PDF}\left(d_c,d_0,\sigma \right)=\frac{1}{\sqrt{2\mathsf{\pi }}\cdot d_c\cdot \sigma }\cdot \exp \left(-\frac{{\ln }^2\left(d_c/d_0\right)}{2{\sigma }^2}\right) ,$$

(5)

with the median diameter $d_0$ and the standard deviation $\sigma$ 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:

$$m_{tot}=N_p{{\displaystyle \int }}_0^{\infty }d{d}_c·\mathrm{PDF}\left(d_c\right)·m_p\left(d_c\right)·ℒ \left(\frac{M_s\pi d_c^3}{6k_{B}T} {\mu }_{0}H\right) .$$

(6)

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:

$$\begin{array}{l}{m}_{tot}=\frac{N_p{M}_s\pi }{6}{\displaystyle \underset{0}{\overset{\infty }{\int }}}d{d}_c·\mathrm{PDF}\left(d_c\right)·d_c^3·\\ ·ℒ \left(\frac{M_s\pi d_c^3{\mu }_0}{6k_{B}T} \left[H_0+H_1\sin \left(2\pi f_{1}t\right)+H_2\sin \left(2\pi f_{2}t\right)\right]\right)\end{array}$$

(7)

As shown in [5], the nonlinearity of the magnetization curve gives rise to the emergence of intermodulation products m⋅f₁ ± n⋅f₂ of the total magnetic moment of the particle ensemble (with m and n denoting integers). In particular, the frequency mixing components f₁ + f₂, f₁ + 2⋅f₂, f₁ + 3⋅f₂ and f₁ + 4⋅f₂ appear. In the limit of small excitation amplitudes H₁ and H₂, 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₁ and H₂, 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\left(d_c\right)$ at frequency f₁ + f₂ of one particle with diameter $d_c$ is given by:

$$\begin{array}{l}{m}_1\left(d_C\right)=\frac{M_s\pi d_c^3}{6}·\frac{2}{k}{\displaystyle \sum _{i=0}^k}\cos \left[2\pi \left(f_1+f_2\right)t_i\right]· \\ ·ℒ \left(\frac{M_s\pi d_c^3{\mu }_0}{6k_{B}T} \left[H_0+H_1\sin \left(2\pi f_1{t}_i\right)+H_2\sin \left(2\pi f_2{t}_i\right)\right]\right)\end{array}$$

(8)

The factor 2/k normalizes the sum and accounts for the full-period average over sin²(..) which is ½. 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₁, i.e., Δt = t_i – t_i–1 = 0.1/f₁, yielding sufficient numerical precision. Although in our experiments, the high frequency f₁ was 485 times larger than the low frequency f₂, it was sufficient to select f₁ = 20⋅f₂ for our numerical calculations. Thus, k = 200 was used.

In a similar fashion, the response component $m_2\left(d_c\right)$ at frequency f₁ + 2⋅f₂ is obtained from:

$$\begin{array}{l}{m}_2\left(d_c\right)=\frac{M_s\pi d_c^3}{6}·\frac{2}{k}{\displaystyle \sum _{i=0}^k}\sin \left[2\pi \left(f_1+2f_2\right)t_i\right]· \\ ·ℒ \left(\frac{M_s\pi d_c^3{\mu }_0}{6k_{B}T} \left[H_0+H_1\sin \left(2\pi f_1{t}_i\right)+H_2\sin \left(2\pi ·2f_2·t_i\right)\right]\right) .\end{array}$$

(9)

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\left(d_c\right)$ and $m_4\left(d_c\right)$ are calculated similarly.

The total magnetic moment component $m_{n,tot}$ at the frequency mixing component f₁ + n⋅f₂ is then obtained by integration over the lognormally weighted particle ensemble:

$$m_{n,tot}\left(d_0,\sigma ,N_p\right)=N_p{{\displaystyle \int }}_0^{\infty }d{d}_c·\mathrm{PDF}\left(d_c,d_0,\sigma \right)·m_n\left(d_c\right).$$

(10)

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₀ and width σ, and the total number of particles $N_p$. The measured nonlinear magnetic moment components of nanoparticle samples at frequencies f₁ + n⋅f₂, n = 1, 2, 3, 4, were fitted with this model using the Levenberg–Marquardt least-squares algorithm.

2.3. 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_{\mathrm{p}}$, is governed by the Landau–Lifshitz–Gilbert equation (LLG) [18]:

$$\frac{d{\mathbf{m}}_p}{dt}=\frac{{\mu }_0\gamma }{1+{\alpha }^2}\cdot \left({\mathbf{H}}_{\mathrm{eff}}\times {\mathbf{m}}_p+\alpha {\mathbf{m}}_p\times \left({\mathbf{H}}_{\mathrm{eff}}\times {\mathbf{m}}_p\right)\right)$$

(11)

with the permeability of free space, ${\mu }_0$, the electron gyromagnetic ratio, $\gamma$, the damping parameter, $\alpha$, and the effective field $H_{\mathrm{eff}}$. The Brownian rotation of the MNP easy axis, $\mathbf{n}$, is described via the generalized torque, $\Theta$, as follows [19]:

$$\frac{d\mathbf{n}}{dt}=\frac{\mathbf{\Theta }}{6\eta V_H}\times \mathbf{n}$$

(12)

with the carrier matrix viscosity, $\eta$, and the MNP hydrodynamic volume, $V_h=\pi /6\cdot d_h^3$, in which $d_h$ is the hydrodynamic particle diameter. Néel and Brownian relaxation are coupled in the internal particle energy:

$$U=-{\mu }_0\cdot m_p\left({\mathbf{m}}_p\cdot H\right)-K\cdot V_c{\left({\mathbf{m}}_p\cdot \mathbf{n}\right)}^2+{\epsilon }_{IA}$$

(13)

where $m_p=\left|m_{\mathrm{p}}\right|=V_c\cdot M_S$ (cf. Equation (2)) is the magnitude of the MNP magnetic moment, and $V_c=\pi /6\cdot d_c^3$ 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, ${\epsilon }_{\mathrm{therm}}$, dominates magnetic interaction energy by two orders of magnitude in our MNP samples, so that particle–particle interaction is negligible (${\epsilon }_{IA}\ll {\epsilon }_{\mathrm{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_{\mathrm{th}}$ and ${\Theta }_{\mathrm{th}}$, which are implemented as Gaussian-distributed white noise with zero mean values ($\langle H_{th}^{\mathrm{i}}\left(t\right)\rangle =0$ and $\langle {\Theta }_{th}^{\mathrm{i}}\left(t\right)\rangle =0$) and variances $\langle H_{\mathrm{th}}^{\mathrm{i}}\left(t\right)H_{\mathrm{th}}^{\mathrm{j}}\left(t’\right)\rangle =\left(2k_B T\cdot \left(1+{\alpha }^2\right)\right)/\left(\gamma m_p\alpha \right)\cdot {\delta }_{ij}\delta \left(t-t^{\prime }\right)$ and $\langle {\Theta }_{\mathrm{th}}^{\mathrm{i}}\left(t\right){\Theta }_{\mathrm{th}}^{\mathrm{j}}\left(t’\right)\rangle =12k_{B}T\eta V_H\cdot {\delta }_{ij}\delta \left(t-t^{\prime }\right)$, respectively. Here $T$ is the global temperature in the system. With this, the effective field and generalized torque read:

$$H_{\mathrm{eff}}=-\frac{1}{m_p\cdot {\mu }_0}\cdot \frac{\partial U}{\partial m}+H_{th}=H+\frac{2K\cdot V_c}{m_p\cdot {\mu }_0}\cdot \left(m_{\mathrm{p}}\cdot n\right)n+H_{th}$$

(14)

$$\Theta =\frac{\partial U}{\partial n}\times n+{\Theta }_{th}=-2K\cdot V_c\left(m_{\mathrm{p}}\cdot n\right)\left(m_{\mathrm{p}}\times n\right)+{\Theta }_{th}$$

(15)

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. Simulation parameters applied in the micromagnetic simulation model.

Effective Anisotropy Constant	Saturation Magnetization 1 MS	Mass Density of Magnetite	Viscosity of Surrounding (Water)	Temperature
11 kJ/m3	476 kA/m	5.2 g/cm3	$8.9\times {10}^{-4}$Pa$\cdot$s	300 K

¹ The literature value for bulk magnetite from [24] was used.

with the MNP properties matching the experimentally determined values from

Table 2. Material properties of MNPs from fitting the experimental data with the Langevin model.

Core Diameter ${\mathbf{d}}_{\mathbf{c}}$	Log-Normal Distribution Width $\mathbf{\sigma }$	Polydispersity Index (PDI)	Hydrodynamic Diameter 2 ${{\mathbf{d}}_{\mathbf{h}}}^{ }$	Concentration 1
7.81 nm	0.346	0.127	20 nm	2.4 mg(Fe)/mL

¹ 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. ² The hydrodynamic diameter d_h is taken from the datasheet of the manufacturer.

. Furthermore, the damping parameter $\alpha$ 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/µ₀, $f_1=40,000$ Hz, $H_2=16$ mT/µ₀, $f_2=2000$ Hz, and the static field was varied from $H_0=\left(0,\dots ,24\right)$ mT/µ₀ with a step size of $1$ mT/µ₀.

3. Results

3.1. 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 least-square 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\cdot 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₁ + f₂ at high offset field did the simulation predict slightly higher values than measurement, whereas for f₁ + 2⋅f₂ and f₁ + 3⋅f₂, the simulation slightly underestimated measurements.

Fitting yielded the MNP material properties of median core diameter, $d_c$, and its log-normal distribution width, $\sigma$, 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.

3.2. 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.

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>\left(H_1+H_2\right)$, e.g., for $H_0=20$ mT/µ₀ 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\cdot 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₁ + f₂ that MC-simulation consistently underestimates the measurement results. However, the peak at $H_0=15$ mT/µ₀ coincides. For the mixing term f₁ + 2⋅f₂, MC-simulation underestimates the measurement results for $H_0<18$ mT/µ₀ and $H_0>20$ mT/µ₀ and shows a mismatch of over 50% for the peak value at $H_0=15$ mT/µ₀. MC-simulation and measurement of component f₁ + 3⋅f₂ coincide for values up to $H_0=10$ mT/µ₀ and again for $H_0\ge 20$ mT/µ₀. However, around the peak at $H_0=16$ mT/µ₀, MC-simulations overestimate experimental data by approx. 70%. For mixing term f₁ + 2⋅f₂, MC-simulations match experimental data, except around the peaks at $H_0=$11 mT/µ₀ and $H_0=16$ mT/µ₀.

3.3. 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/µ₀, $f_1=40,000$ Hz, $H_2=16$ mT/µ₀, $f_2=2000$ Hz, and fixed the mean core diameter with $d_c=7.81$ nm and variable distribution parameter $\sigma$ 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 $\sigma =1.466$. This is significantly different than the input parameters to MC-simulations ($\sigma =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.

4. 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\approx 7.8$ nm and $\sigma \approx 0.35$; cf. Table 2) are furthermore in good agreement with literature values reporting a core diameter between $d_c\approx 7$ nm [28] and $d_c\approx 11$ nm [29] for nanomag^®-D SPIO. The broad size distribution (with width $\sigma >0.3$) reflects a heterogeneous particle size with an effectively range of $d_c=\left(3-25\right)$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 ~ {\tau }^{-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\cdot V_C$.

This can be assessed by assuming an AC-field-dependent Brownian relaxation time of [30]:

$${\tau }_{\mathrm{B}}\left(H_{\mathrm{AC}}\right)=\frac{{\tau }_{\mathrm{B}}}{\sqrt{1+0.21\cdot \xi {\left(H_{AC}\right)}^2}}$$

(16)

and an AC-field dependent Néel relaxation time of [16]:

$${\tau }_{\mathrm{N}}\left(H_{AC}\right)={\tau }_0\cdot \exp \left(\frac{K\cdot V_{\mathrm{C}}}{k_{\mathrm{B}} T}\cdot {\left(1-\frac{H_{AC}}{H_K}\right)}^2\right)$$

(17)

where $H_{AC}$ denotes the AMF amplitude, ${\tau }_B=3\eta V_{\mathrm{H}}/\left(k_{\mathrm{B}}T\right)$ the field-independent Brownian relaxation time, $\xi \left(H_{AC}\right)=m{H}_{AC}/\left(k_{\mathrm{B}}T\right)$ the reduced field parameter (cf. Equation (4)), ${\tau }_0={10}^{-9}$ s the time constant, and $H_K=2K/{\mu }_0{M}_S$ the (uniaxial) anisotropy field strength. From Equations (16) and (17), the effective relaxation time follows with:

$$\tau =\frac{1}{{\tau }_B^{-1}+{\tau }_N^{-1}}.$$

(18)

Using the (mean) values from Table 1 and Table 2, and the experimental setup parameters with $H_{AC}=H_1+H_2 ~ 17$mT/µ₀, one calculates $\tau ~ {\tau }_0 ~ {10}^{-9}$ s so that for $f ~ \left(30-40\right)$ kHz the condition for the onset of hysteresis, $f ~ {\tau }^{-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 \ge 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³ [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 ~ {\tau }^{-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, Figure 2 and Figure 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 \ge 12.1$ nm. Similarly, 99% of the signal stems from $d_c \ge 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, \sigma , 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 \ge 20$nm at frequencies $f ~ \left(30-40\right)$ 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 $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 Figure 4 and Figure 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.

5. 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:

• Both the Langevin model and MC-simulations reproduce the shape of the experimental signal satisfactorily (see Figure 1 and Figure 3). However, punctual deviations between experimental data and MC-simulations are observed (see Figure 3).
• 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.
• 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).
• 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 particle–particle 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 dual-frequency 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 multi-parameter 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.