Computer Simulations of Dynamic Response of Ferroﬂuids on an Alternating Magnetic Field with High Amplitude

: The response of ferroﬂuids to a high-amplitude AC magnetic ﬁeld is important for several applications including magnetic hyperthermia and biodetection. In computer simulations of the dynamic susceptibility of a ferroﬂuid outside the linear response region, there are several problems associated with the fact that an increase in the frequency of the AC ﬁeld leads to the appearance of additional computational errors, which can even lead to unphysical results. In this article, we study the dependence of the computational error arising in the computer simulation of the dynamic susceptibility on the input parameters of the numerical algorithm: the length of the time step, the total number of computer simulation periods, and averaging period. Computer simulation is carried out using the Langevin dynamics method and takes Brownian rotational relaxation of magnetic particles and interparticle interactions into account. The reference theory [Yoshida T.; Enpuku K. Jap. J. Ap. Phys. 2009] is used to estimate computational error. As a result, we give practical recommendations for choosing the optimal input parameters of the numerical algorithm, which make it possible to obtain reliable results of the dynamic susceptibility of a ferroﬂuid in a high-amplitude AC ﬁeld in a wide frequency range.


Introduction
In recent years, there has been growing interest in smart materials that can be controlled by external stimuli such as light, temperature, mechanical force, electrical field, and magnetic field. These materials usually contain active fillers or microstructures in various types of soft or liquid matrices [1][2][3]. Colloidal suspensions of magnetized nanoparticles in liquid nonmagnetic carriers are known as ferrofluids. The ability of ferrofluids to respond to external magnetic fields makes them useful in many areas, ranging from engineering [4], to biomedical applications [5,6]. Such materials can be used as an active medium in devices (sealants, heat conduction media, media for hydraulic suspensions) and in materials processing (separation media, gas-fluidized beds). They can also be adapted for biomedical uses such as targeted drug delivery, diagnosis, and localized treatment by hyperthermia [7][8][9] as well as noninvasive tomographic techniques such as magnetic particle imaging (MPI) [10][11][12]. The main characteristic that describes the dynamic magnetic response of a ferrofluid to an AC field is the dynamic susceptibility. Most of the known theoretical models of the dynamic magnetic response of ferrofluids concern noninteracting samples [13][14][15][16][17][18][19][20][21]. For example, the dynamic susceptibility of a ferrofluid in a weak AC magnetic field is described by the well-known Debye formula [15,16]. Yoshida and Enpuku proposed simple analytical expressions for determining the dynamic response of a ferrofluid to an AC field with arbitrary amplitude [17]. The dynamic susceptibility of ferrofluids under magnetic fields of various configurations was considered in [18][19][20]. The aforementioned theoretical models successfully describe experimental results for diluted ferrofluids; however, they cannot describe the behavior of concentrated samples. In concentrated ferrofluids, interparticle interactions play an important role and significantly change the magnetic susceptibility of the samples [22,23]. The influence of interparticle interactions on the dynamic response of a ferrofluid was investigated in theoretical models [24][25][26][27]. In these works, two different approaches were used to take interparticle interactions into account: modified mean field theory and virial expansion. Theories [24][25][26] show a good agreement with computer simulation results for dynamic susceptibility of moderately concentrated ferrlfluids in AC field with low amplitude. In [27], a theory of the dynamic response of a ferrofluid to an AC magnetic field that includes both the effects of interparticle dipole-dipole interactions and the dependence on field amplitude is obtained.
Computer simulation of the dynamic response of a ferrofluid is usually based on linear response theory and the magnetic susceptibility is expressed in terms of the equilibrium (zero-field) magnetization autocorrelation function [28]. Using this approach, the dynamic response of a ferrofluid can predict reliably only for weak AC fields, as detailed in [29,30], for example. The response of ferrofluids to an AC magnetic field outside the linear response region is important for several applications including magnetic hyperthermia, magnetic particle imaging (MPI), and biodetection. The problem is that when we deal with the high-amplitude AC fields in computer simulations, an increase in the frequency of the AC field leads to the appearance of additional computational errors, which can even lead to physically incorrect results (for example, a negative value of the magnetic susceptibility). However, these errors must be corrected by carefully adjusting the input parameters of the numerical algorithm: decreasing the time step, increasing the total number of computational periods, complicating the calculation method, etc. In the literature, there are no studies and practical recommendations about the subject of how to choose the input parameters of the numerical algorithm to minimize computational errors and correctly determine the dynamic susceptibility in a high-frequency region. In this regard, obtaining reliable results of computer simulation of the dynamic susceptibility of a ferrofluid under a high-amplitude AC field in a wide frequency range is still a challenge.
In this article, we carry out the Langevin dynamic (LD) simulation of the dynamic susceptibility of a ferrofluid in a high-amplitude AC field using different time steps in the numerical scheme, changing the number of calculated periods. We compare the simulation results with the known Yoshida-Enpuku theory [17] in order to determine the optimal input parameters of the numerical algorithm, which make it possible to obtain reliable data of the dynamic susceptibility of a ferrofluid in a high-amplitude AC field in a wide frequency range. We also analyze the dependence of the time step inputted into the numerical algorithm on the magnitude of the computational error arising in the computer simulation of the dynamic susceptibility for various frequencies and amplitudes of the AC field and give practical recommendations for carrying out the calculations. Using these data, the dynamic susceptibility of a ferrofluid with different viscosity, concentration, and intensity of interparticle dipole-dipole interactions is calculated in a wide range of amplitudes and frequencies of the AC field.

Model
Ferrofluid is modeled as a system consisting of N spherical particles of mass m in a three-dimensional volume V at temperature T. We assume that the ith particle exhibits a permanent point dipole moment µ i = µμ i at its center (where µ is a length of a point dipole moment vector andμ i is a unit vector of a point dipole moment), which can freely rotate in 3D. The interaction potential between particles i and j is a sum of shot-range Weeks-Chandler-Andersen (WCA) potential u WCA ij and the long-range dipolar (d) potential u d ij . The WCA potential is where r ij is a distance between the ith and jth ferroparticles, r min = 2 1/6 σ is the minimum of the Lennard-Jones (LJ) potential, and σ are the LJ energy and range parameter, respectively. The magnetic dipole interaction is wherer ij is a unit vector of the interparticle separation vector r ij = r ijrij , µ 0 is the vacuum magnetic permeability. The interaction energy u H i between the ith particle magnetic moment and the magnetic field H = H cos(ωt)Ĥ can be written in the Zeeman form: where t denotes time,Ĥ stands for the unit vector, and ω and H are AC field angular frequency and amplitude, respectively. As usual, reduced units are employed: the reduced concentration is ρ * = ρσ 3 , where ρ = N/V; the reduced magnetic moment is µ * = µ 0 µ 2 /(4πσ 3 ); the reduced temperature is T * = k B T/ (where k B is the Boltzmann constant); the dipolar coupling constant is λ = (µ * ) 2 /T * ; and the reduced time t * = t/t 0 , which is measured in units of t 0 = √ σ 2 m/ .

Simulation Details
The translational and rotational movement of the ith particle is determined by the Langevin equations of motions [31]: where I is the particle moment of inertia, and v i and Ω i are its translational and angular velocities, correspondingly. The resulting forces acting on the ith particle include LJ force F LJ ij , the force of dipole-dipole interaction F d ij , and particle-field force µ 0 µ i × H. The symbols Γ T and Γ R stand for the translational and rotational friction constants, respectively. Gaussian distributed random forces ξ T i and torques ξ R i have the properties that the first moments vanish while the second moments satisfy where φ ∈ {x, y, z} and ψ ∈ {x, y, z}, {x, y, z} are components in the Cartesian coordinates, t and s denote time moments, δ ij and δ ϕψ are the Kronecker delta, and δ(t − s) is the Dirac delta function.
In actual LD simulations, we use the reduced variables where v * i is the reduced translational velocity of the ith particle, F * i is the reduced force of the ith particle, Γ * T is the reduced translational friction constant, ξ T * i is the reduced Gaussian distributed random force for the ith particle, I * is the reduced particle moment of inertia, Ω * i is the reduced angular velocity of the ith particle, τ * i is the reduced torque of the ith particle, Γ * R is the reduced rotational friction constant, and ξ R * i is the reduced Gaussian distributed random torque for the ith particle. Therefore, Equations (4) and (5) can be rewritten in reduced form For the computer simulation of ferrofluids, Equations (6) and (7) are augmented by an equation of the motion of magnetic moments It should be noted that the model (6)-(8) corresponds to a monodisperse ferrofluid. The direct application of this model to experimental results is complicated by the polydispersity of real ferrofluids. Nevertheless, the extension of this model to the polydisperse case is possible by using the discrete ferroparticle diameter distribution. In general, the set of Equations (6)-(8) can be solved using various numerical methods. For translational equations of motion (6), the velocity Verlet algorithm is widely used, which is an example of a second-order integrator [32]. The problems with solving Equations (7) and (8) are that the direct application of integration methods involves using Euler angles to represent rotational degrees of freedom, which leads to singularities when the azimuthal angle of the particle is 0 or π [31,33]. Another problem is that the value of the magnetic moment must be kept constant during numerical integration. To overcome the difficulties mentioned above, the velocity-Verlet-like algorithm is used for the quaternions [34]. For a system of noninteracting magnetic particles, the above-mentioned algorithms for solving Equations (6)- (8) are implemented in ESPResSo [35], which we used for LD simulations. We included the interparticle dipole-dipole interactions in the numerical scheme of ESPResSo with help of the dipolar-P 3 M algorithm [36].
The rotational motion of the particles can be described by two time scales. They are the inertial time τ I = I/ Γ R and the Brownian time τ B = Γ R / 2k B T. First, τ I describes how long a particle resists any change in its angular velocity. Second, τ B is connected with thermal fluctuations. For ferrofluids, the ratio τ I τ B must be fulfilled. To satisfy this condition, we set I * = 1, T * = 1 and Γ * R from the range 20 to 80. A cubic simulation box with periodic boundary conditions was used for simulations. The volume fraction ϕ was varied from 0.01 to 0.125. The strength of the dipole-dipole interaction was defined by λ = 1.
The reduced normalized magnetization is defined as M * (t) was calculated until the moment of time that corresponds to 22 periods of the external magnetic field. The results were averaged over 18 periods, while the first four periods were omitted. The average value of the magnetization M * H (t) in the direction of the field was used to calculate the dynamic susceptibility. The real χ (ω) and imaginary χ (ω) parts of the dynamic susceptibility are defined as the first term in the Fourier series of M * H (t): where χ L = µ 0 µ 2 ρ/3k B T is the Langevin susceptibility and α = mH/k B T is the Langevin parameter characterizing the AC field amplitude. χ L is a complex characteristic of ferrofluid density and the intensity of interparticle dipole-dipole interactions. Integrals (9) are calculated numerically using the trapezoidal rule.

Theory
In order to validate the LD simulation methodology, we first consider a system of noninteracting particles. For this system, computer simulation results are compared with the Yoshida-Enpuku theory [17], which is a modification of the Debye theory for the high field amplitudes where χ id (0) is the susceptibility in the quasistatic case (when the frequency approaches zero), τ id is the effective relaxation time, which can be obtained as τ id = 1/ω p (ω p is the frequency at which the imaginary part of the dynamic susceptibility χ id has a peak), and k id denotes the factor. Subscript id denotes the ideal system, i.e., the system without dipolar interactions between magnetic particles. Expressions (10) were obtained by approximating the numerical solutions of the Fokker-Planck equation [37] for a single dipole particle in an AC magnetic field. Real and imaginary parts of the susceptibility (10) reliably describe the susceptibility of a noninteracting particles, at least for α ≤ 10 [17]. Figure 1 shows dynamic magnetization curves for the systems of noninteracting particles in the AC field with amplitude α = 5 and angular frequencies ωτ B = 0.1 (Figure 1a) and ωτ B = 20.0 (Figure 1b). Each magnetization curve contains 22 periods of the AC field. For the low frequency ωτ B = 0.1, magnetization immediately goes into a steady state and the values of the magnetization is similar in each period, whereas for the high frequency ωτ B = 20.0, the magnetization reaches the steady state slowly. Therefore, in the high-frequency region, it is necessary to omit several initial periods to obtain good statistical data. It is worth to note that the number of periods that should be omitted for each system is different and depends on the time when the system achieves the steady state. In this work, to calculate the dynamic susceptibility using Equation (9), the magnetization was averaged over 18 periods (the first four periods were omitted). This is the minimum number of calculation periods that made it possible to obtain the best agreement between the computer simulation results and the reference Yoshida-Enpuku theory (10) in the whole range of considered parameters ωτ B ≤ 10 2 , 5 ≤ α ≤ 10. The deviation of the computer simulation results from the theoretical data (10) for the real and imaginary parts of the susceptibility depending on the number of computational periods is presented in Figure 2. The vertical axis shows the number of periods of the AC field over which the averaging of the magnetization was carried out for calculating the susceptibility; the horizontal axis is the number of the first periods of the calculation, which are omitted during averaging. Colors show the absolute difference between the results of the computer simulation and the Yoshida-Enpuku theory for the system with ωτ B = 20.0, α = 5, Γ * R = 40. When the total number of calculation periods is less than 20, the averaging result is sensitive to the number of omitted periods; there is an intermittent color change in this region, which is more pronounced for the imaginary part. An increase in calculation periods over 20 leads to a smoother change in colors on the diagram, and the dependence on omitted periods becomes weak, that is, the results of computer simulations are stabilized. Nevertheless, in order to exclude the influence of the unstable behavior of magnetization at small t * (Figure 1b), we omitted the first four periods when calculating the susceptibilities.

Results and Discussion
The magnetization curves averaged over 18 periods are used to calculate the dynamic susceptibility (9). Figure 3 shows the comparison of the Yoshida-Enpuku theory (10) with LD simulation results for dynamic susceptibility of noninteracting particles in the AC field with amplitudes α = 5 ( Figure 3a) and α = 10 ( Figure 3b). Lines are from the Yoshida-Enpuku theory (10); symbols show the computer simulation data. Filled symbols correspond to the real part of a dynamic susceptibility; unfilled ones are the imaginary part. Computer simulations were performed with different time steps h t , so squares are chosen for h t = 0.01, circles are used for h t = 0.0025, rhombi correspond to h t = 0.001, and triangles stand for h t = 0.00025.
The maximum deviation of the LD simulation results from the theoretical data is observed in the frequency range ωτ B ∈ (2, 100) at h t = 0.01, where the real part of the susceptibility in computer simulation demonstrates nonphysical negative values. Note that at the same time step, h t = 0.01 for the frequencies ωτ B ∈ (10 −2 , 2), the results of the computer simulations are in satisfactory agreement with the theory. Decreasing the time step up to h t = 0.00025 allows us to achieve the quantitative agreement between theoretical curves (solid and dashed lines) and computer simulation results (triangles) in the wide frequency range ωτ B ≤ 10 2 .    In Figure 4, the difference between theoretical and computer simulation results is shown for different time steps. For all time steps h t presented in Figure 4, the computer simulation results χ CS and χ CS have the same deviation from theoretical data χ T and χ T for frequencies ωτ B 1. The error begins to behave differently with increasing frequency, which is most noticeable in the range of ωτ B > 10. As it is seen, the maximum difference is observed for the biggest time step h t = 0.01. For the high frequency ωτ B , the smallest difference is reached for the smallest time step h t = 0.00025. In other words, decreasing the time step in the high-frequency region allows us to reduce statistical errors.  The step size h t is not the only parameter that affects the accuracy of the obtained computer simulation data of the dynamic susceptibility. Another characteristic is related to the calculation of the integral Equation (9). This integral was calculated using the trapezoid rule [38] where {t k :t 0 = 0,t k =t k−1 + h int ∀k ∈ {1, .., K − 1},t K = 2π/ω}. Discrete time points {t k } from the numerical solution of Equations (6) and (7) and the points {t k } for calculating Equation (9) do not have to coincide surely. If {t k }={t k }, then the saving of the magnetization value should be at each step of solving Equations (6) and (7). This approach is more accurate, but it is slower, and it leads to increased memory consumption for storing data. We performed calculations for h int = Sh t , where S is the aspect ratio between these sets of points. It is worth to note that in Figures 2 and 3, h int = 5h t . The influence h int on the accuracy of the computer simulation results is clearly pronounced in Figure 5 for the time step h t = 0.001. The results are presented for different magnetic field amplitudes α = 5 in Figure 5a and α = 10 in Figure 5b. Blue squares are chosen for h int = 0.1; green circles are for h int = 0.005. It is clearly seen that the difference between the computer simulation data for different h int can reach the value 0.1, which corresponds to the order of h int .
The absolute deviation is shown by horizontal dashed lines. Susceptibility values correspond to the system from Figure 3b, Thus, Figures 4 and 6 show a detailed analysis of the computational error versus a field frequency and time steps. Since LD simulations are time-consuming, based on Figures 4 and 6, the following practical recommendations can be given for calculating the dynamic susceptibility of a ferrofluid in a high-amplitude AC field.

1.
The Langevin dynamics equations should be solved until the time moment that corresponds to at least 22 periods of the AC magnetic field. Obtained dynamic magnetization should be averaged over all calculated periods, excluding at least the first four ones.

2.
It is necessary to use different time steps in different frequency ranges : h t = 0.0025 if The step size h int has to be the same order as h t -for example, h int = 5h t .
These recommendations with a time step decrement starting from a certain frequency were used to determine the dynamic susceptibility of an ensemble of noninteracting magnetic particles distributed in a medium with different viscosities of a carrier liquid. The viscosity is modeled by friction constant Γ * R in Equations (6)- (8). The results for different friction constants are shown in Figure 7. The results for α = 5 can be found in Figure 7a; Figure 7b shows the data for α = 10. Lines correspond to the theoretical results obtained using Equation (10); symbols are chosen for the computer simulation data. Blue squares and lines correspond to the real and imaginary susceptibility for the friction constant Γ * R = 20. For Γ * R = 40, we used green circles and lines, and red rhombi and lines are chosen for the real and imaginary susceptibility for Γ * R = 80. Increasing viscosity of the medium leads to a shift in the maximum of the imaginary part of the susceptibility to the low-frequency region; the value of the maximum does not depend on the medium viscosity. The plateau in the region of low frequencies of the real part of the susceptibility is longer for systems with lower viscosity. An increase of the amplitude of the AC field leads to a decrease in the values of the real and imaginary parts of the susceptibility, and hence to a decrease in the dynamic magnetic response of the system. The latter is confirmed by experimental data [39]. In Figure 8, the dynamic susceptibility of an ensemble of interacting magnetic particles is shown. Computer simulation results plotted by symbols are shown in comparison with the reference Yoshida-Enpuku theory (10) for noninteracting particles plotted by lines. Blue circles correspond to the real part of the susceptibility; orange squares are chosen for the imaginary part. Solid and dashed black lines are the real and imaginary parts of the susceptibility obtained in the reference theory correspondingly. Figure 8a contains the data for the Langevin parameter α = 5; Figure 8b is for α = 10. In Figure 8b, it is seen that the dipole-dipole interaction is suppressed by the dipole-field energy, and the susceptibilities of the interacting and noninteracting systems almost coincide. A similar behavior of the magnetic susceptibility at large amplitudes was described in the work for a system of immobilized interacting superparamagnetic particles [40]. With an increase of the amplitude of the AC field, dipole-dipole interactions begin to appear; that leads to an increase of the susceptibility, which is demonstrated in Figure 8a.

Conclusions
In this work, we have shown that the arising computational errors in the LD computer simulation of the dynamic susceptibility of a ferrofluid in a high-amplitude alternating magnetic field at high frequencies can be minimized by carefully adjusting the input parameters of the numerical algorithm, namely the time step h t for the solution stochastic Equations (6) and (7), the time step h int for calculation of integrals (9), and the number of calculation periods. Despite the fact that the exact values of input parameters depend on the system under consideration, we have given practical recommendations that make it possible to obtain reliable results of the dynamic susceptibility of a ferrofluid in an AC field of high amplitude 5 ≤ α ≤ 10 in a wide frequency range, at least ωτ B ≤ 10 −2 . As a criterion for assessing the reliability of the obtained computer simulation results, we chose a good agreement between the computer simulation data of dynamic susceptibility and the Yoshida-Enpuku theory (10). Recommendations are described in detail in the Results and Discussion sections and include (1) using different time steps in different frequency ranges for calculating dynamic magnetization, (2) solving the Langevin dynamics equations until the moment of time that corresponds to at least 22 periods of the external magnetic field, (3) using the averaged magnetization over all calculated periods, excluding the first few periods, and (4) increasing the integration step h int in comparison with the step h t . All of these recommendations allow us to accumulate reliable statistics. Using these recommendations, we simulate the dynamic susceptibility of a ferrofluid with different viscosity, concentration, and intensity of interparticle dipole-dipole interactions in a wide range of amplitudes and frequencies of the AC field. Acknowledgments: Computer simulations were performed using the Ural Federal University Scientific Cluster.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: The reduced temperature k B The Boltzmann constant λ The dipolar coupling constant t * The reduced time t 0 The units for time I The particle moment of inertia v i The particle translational velocity Ω i The particle angular velocity F

LJ ij
The Lennard-Jones force The Gaussian distributed random force for the ith particle The Gaussian distributed random torque for the ith particle The Cartesian coordinate φ ∈ {x, y, z} of the Gaussian distributed random force for the ith particle The Cartesian coordinate ψ ∈ {x, y, z} of the Gaussian distributed random torque for the ith particle s The time moment v * i The reduced translational velocity of the ith particle F * i The reduced force of the ith particle Γ * T The reduced translational friction constant The reduced Gaussian distributed random force for the ith particle I * The reduced particle moment of inertia Ω * i The reduced angular velocity of the ith particle The reduced torque of the ith particle Γ * The time step for calculating integral (9) {t k } The time points for calculating integral (9) K The number of segments for calculating integral (9) S The aspect ratio between sets of points {t k } and {t k }