Multiparticle collision dynamics for ferrofluids

Multiparticle collision dynamics for ferrofluids Cite as: J. Chem. Phys. 156, 144905 (2022); https://doi.org/10.1063/5.0087981 Submitted: 11 February 2022 • Accepted: 28 March 2022 • Accepted Manuscript Online: 29 March 2022 • Published Online: 12 April 2022 Patrick Ilg ARTICLES YOU MAY BE INTERESTED IN Multi-particle collision dynamics with a non-ideal equation of state. II. Collective dynamics of elongated squirmer rods The Journal of Chemical Physics 155, 134904 (2021); https://doi.org/10.1063/5.0064558 Evolution of jets during drop impact on a deep liquid pool Physics of Fluids 34, 022110 (2022); https://doi.org/10.1063/5.0081064 Transient oscillation response characteristics of an electrohydrodynamic settling drop subjected to a uniform electric field Physics of Fluids 34, 043601 (2022); https://doi.org/10.1063/5.0086168 J. Chem. Phys. 156, 144905 (2022); https://doi.org/10.1063/5.0087981 156, 144905 The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp Multiparticle collision dynamics for ferrofluids Cite as: J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Submitted: 11 February 2022 • Accepted: 28 March 2022 • Published Online: 12 April 2022 Patrick Ilga) AFFILIATIONS School of Mathematical, Physical, and Computational Sciences, University of Reading, Reading RG6 6AX, United Kingdom a) Author to whom correspondence should be addressed: [email protected] ABSTRACT Detailed studies of the intriguing field-dependent dynamics and transport properties of confined flowing ferrofluids require efficient mesoscopic simulation methods that account for fluctuating ferrohydrodynamics. Here, we propose such a new mesoscopic model for the dynamics and flow of ferrofluids, where we couple the multi-particle collision dynamics method as a solver for the fluctuating hydrodynamics equations to the stochastic magnetization dynamics of suspended magnetic nanoparticles. This hybrid model is validated by reproducing the magnetoviscous effect in Poiseuille flow, obtaining the rotational viscosity in quantitative agreement with theoretical predictions. We also illustrate the new method for the benchmark problem of flow around a square cylinder. Interestingly, we observe that the length of the recirculation region is increased, whereas the drag coefficient is decreased in ferrofluids when an external magnetic field is applied compared with the field-free case at the same effective Reynolds number. The presence of thermal fluctuations and the flexibility of this particle-based mesoscopic method provide a promising tool to investigate a broad range of flow phenomena of magnetic fluids, and the method could also serve as an efficient way to simulate solvent effects when colloidal particles are immersed in ferrofluids. Published under a nonexclusive license by AIP Publishing. https://doi.org/10.1063/5.0087981 I. INTRODUCTION Colloidal suspensions of magnetic nanoparticles, also known as ferrofluids, are fascinating model systems that combine superparamagnetic and magnetoviscous effects.1,2 The classical experiments of McTague3 demonstrated that pipe flow of ferrofluids can be manipulated by external magnetic fields, showing the anisotropic response depending on the relative orientation of the magnetic field with respect to the flow direction. Martsenyuk et al.4 successfully explained these experiments with a kinetic theory of hindered particle rotations and the resulting changes in rotational viscosity. These ground-breaking works sparked numerous subsequent studies on the magnetoviscous effect, leading to more refined experiments and a better theoretical understanding as well as several novel applications of ferrofluids that rely on their field-dependent flow behavior (see, e.g., Refs. 5–9 and references therein). In addition, various computer simulation methods for ferrofluid flow have been developed in recent years. The simulations range from simple channel or pipe flows10–12 to more complicated geometries13 including free-surface flows.14 Different numerical methods were used in these studies, ranging from perturbative solutions10 to adaption of finite-element and finite-volume computational fluid dynamics codes,11,13 the lattice Boltzmann method,12 and smooth particle hydrodynamics simulations.14 Such flow simulations are helpful for various J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing applications, e.g., planning magnetic drug targeting treatments.15 In addition, micro-magnetofluidics shows promising potential for contactless mixing, separation, and trapping of particles and polymers on small scales, where flow simulations help to improve their effectiveness.16 Thermal fluctuations are important in small-scale flows which, however, are neither captured by traditional computational fluid dynamics codes nor in the lattice Boltzmann method. Instead, multi-particle collision (MPC) dynamics provides a very versatile and flexible method for simulating fluid flow including thermal fluctuations.17–19 In this off-lattice method, each particle represents a small fluid element and evolves through a sequence of streaming and simplified collision steps. By locally conserving mass, momentum, and energy, the correct hydrodynamic behavior is generated. In addition, thermal fluctuations are naturally present in this particlebased scheme. Recently, the MPC method has been extended to model anisotropic fluids such as nematic liquid crystals.20–22 Here, we use similar reasoning to extend the MPC method to model ferrofluid flow including stochastic magnetization dynamics. We note that none of the previous approaches to ferrofluid flow simulations10–14 include thermal fluctuations in the magnetization dynamics. For simplicity, we restrict ourselves to a quasitwo-dimensional system, i.e., two-dimensional translational motion with three-dimensional magnetization dynamics. We demonstrate 156, 144905-1 The Journal of Chemical Physics ARTICLE that this model is able to reproduce fluctuating ferrohydrodynamics. In particular, we verify that the field-dependent effective viscosity resulting from this hybrid method quantitatively agrees with theoretical predictions. Furthermore, we show the flexibility of this approach by simulating the flow around a square cylinder. Analyzing the flow field, we determine the length of the recirculation region and find it decreasing with increasing field strength. Moreover, we also calculate the drag coefficient and study its dependence on the magnetic field. The new hybrid method presented here is a useful tool that also allows us to study the effect of fluctuations on ferrofluid flow on a mesoscopic scale. Different models of magnetization dynamics can be incorporated in a straightforward manner. Thereby, the flexibility of the MPC method also facilitates the study of more complicated geometries. This paper is organized as follows. In Sec. II, we briefly review the ferrohydrodynamic theory. The hybrid MPC model coupled to stochastic magnetization dynamics is described in Sec. III. Section IV describes the simulation setup before results are presented and discussed for channel flows in Sec. V and for flow around square cylinder in Sec. VI. Finally, some conclusions are offered in Sec. VII. II. FERROHYDRODYNAMICS We here give a brief overview of the basic equations of ferrohydrodynamics to make the paper self-contained. Further details can be found, e.g., in Refs. 1, 23, and 24. The fluid momentum balance equation reads dv ≙ −∇p/ρ + ν∇2 v + fM , dt ∇ × H ≙ 0, ∇ ⋅ B ≙ 0, (2) (3) where B ≙ μ0 (H + M) denotes the magnetic induction and μ0 denotes the permeability of free space. Equations (1)–(3) are not closed due to the appearance of the magnetization M. Some of the previous simulation approaches have assumed quasi-equilibrium conditions and locally approximate M ≈ Meq (H).13,14 Such approaches, however, neglect relaxation phenomena, and a more accurate treatment of the magnetization dynamics is often desirable. Unfortunately, despite many efforts and long debates in the literature about the appropriate form of the magnetization equation for ferrofluids, no consensus has been reached yet (see, e.g., Refs. 8 and 23–26 and references therein). One of the advantages of the present method is that it can handle different magnetization equations rather straightforwardly. Since J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing − ζ(ωi − Ω(ri )) + μui × H(ri ) + TBi ≙ 0, (4) where μui and ωi are the magnetic moment and angular velocity of magnetic nanoparticle i, respectively. The magnitude of the magnetic moment is denoted by μ. The rotational viscous torque is proportional to the rotational friction coefficient ζ and arises when the particle rotation does not match the local vorticity, Ω ≙ (1/2)∇ × v. The magnetic torque is generated by the magnetic field H at the location of the nanoparticle. Brownian fluctuations are modeled by random torques TBi with Gaussian white noise, ⟨TBi ⟩ ≙ 0 and ⟨TBi (t)TBj (t ′ )⟩ ≙ 2kB Tζδij δ(t − t ′ )I, kB and T being Boltzmann’s constant and temperature, respectively, and I being the three-dimensional unit matrix. The fluid magnetization is calculated from local averages over the magnetic moments of the nanoparticles as M(r, t) ≙ ⟨∑i μui (t)w(r − ri (t))⟩, where w is a suitable weight function. The magnetization arising from an applied field is responsible for the backflow effect via the magnetic force density given in Eq. (2). In equilibrium, the model describes superparamagnetic behavior, Meq ∼ L(H)H/∣H∣ with L being the Langevin function.1 III. MODEL AND NUMERICAL IMPLEMENTATION where the magnetic field and the magnetization are denoted by H and M, respectively. We assume the fluid to be incompressible and non-conducting; therefore, ∇ ⋅ v ≙ 0, the scheme naturally includes thermal fluctuations in the hydrodynamic variables, we here choose to implement a mesoscopic model that also accounts for thermal fluctuations in the magnetization. In particular, we employ the classical model of ferrofluid dynamics of rigid dipoles proposed in Ref. 4 and widely studied since.8,27,28 In the free-draining limit and neglecting inter-particle interactions, the torque balance of viscous, magnetic, and Brownian torques reads (1) with v being the fluid velocity field, dv/dt denoting the material derivative, ρ and ν being the fluid density and kinematic viscosity, respectively, and p being the scalar pressure. With allowance for internal rotations, the influence of magnetization effects on fluid flow is described by the force density, 1 ρfM ≙ (M ⋅ ∇)H + ∇ × (M × H), 2 scitation.org/journal/jcp MPC is a very versatile and flexible, particle-based method to simulate fluid flow. The method has recently been extended to complex fluids19 such as polymer solutions29 and, in particular, also to anisotropic fluids such as nematic liquid crystals.20–22 Inspired by these recent studies, we here present an application to another class of anisotropic fluids, namely, ferrofluids. Within MPC, the fluid is represented by a system of N identical particles, each with mass mi ≙ m, where ri and vi denote the position and velocity of particle i, respectively, and i ≙ 1, . . . , N. Here, we follow previous studies20–22,29 and let an MPC particle represent a coarse-grained description of a fluid element, containing solvent as well as solute particles. Thus, each MPC particle also carries a magnetic moment μui , where μ denotes the magnitude and the three-dimensional unit vector ui describes the orientation of the magnetic moment of particle i. The basic idea underlying MPC is that particle dynamics can be split into a streaming and a simplified collision step; both can be computed very efficiently. For anisotropic fluids, the translational motion needs to be coupled to rotational dynamics as detailed in the following for the case of ferrofluids. In the streaming step, the positions and velocities of the particles are updated according to30 ri (t + Δt) ≙ ri (t) + Δt vi (t) + v′i (t) ≙ vi (t) + Δt Fi (t), mi Δt 2 Fi (t), 2mi (5) (6) 156, 144905-2 The Journal of Chemical Physics ARTICLE with the time step Δt. The force Fi ≙ fext + fM (ri ) acting on particle i consists of a uniform external forcing fext (e.g., due to an applied pressure gradient) and the magnetic force fM defined in Eq. (2). The collision step ensures momentum exchange between particles and can be expressed as21 vi (t + Δt) ≙ VCi (t) + βth R ⋅ ∥v′i (t) − VCi (t)∥, (7) where VCi is the center-of-mass velocity of the collision cell to which particle i belongs, VCi ≙ ∑j∈𝒞 i mj v′j /MCi . The total mass in this collision cell is given by MCi ≙ ∑j∈𝒞 i mj , and 𝒞 i labels all particles currently located in collision cell Ci . In MPC, the collision step (7) is performed simultaneously between all particles j currently residing in the same collision cell. These cells are defined by a regular square grid of length a. Collisions lead to rotation of the relative velocities, which is described by the rotation matrix R in Eq. (7). In two spatial dimensions considered here, R is completely specified by the rotational angle α. Thus, R rotates the relative velocities by an angle ±α with equal probabilities. As has been emphasized before,19 angular momentum conservation is important for the rotational dynamics of anisotropic fluids but is, in general, violated by Eq. (7). Angular momentum conservation can be restored for this scheme; however, if the rotation angle α is not fixed but chosen in collision cell Ci according to cos α ≙ A21 − A22 , A21 + A22 sin α ≙ − 2A1 A2 , A21 + A22 (8) where A1 ≙ ∑j∈𝒞 i ∥rj × ṽj ∥z and A2 ≙ ∑j∈𝒞 i rj ⋅ ṽj with relative velocities ṽi ≙ v′i − VCi .19 In flow simulations, a thermostat is generally needed in order to remove the energy input into the system. Here, we follow common practice (see, e.g., Ref. 21) and use a simple rescaling of the relative √ velocities to the bath temperature T by βth ≙ T/TCi . The instantaneous kinetic temperature in collision cell Ci is defined in two spatial dimensions by kB TCi ≙ (1/2NCi )∑j∈𝒞 i mj ṽ2j with NCi ≙ ∑j∈𝒞 i 1 being the number of particles in the collision cell. Using a fixed grid of collision cells violates Galilean invariance and might lead to spurious correlations. Therefore, following Ref. 18, in each time step, the grid is shifted by a random vector, where each component is uniformly distributed in ∥−a/2, a/2∥. The simulations discussed in the following employ two types of boundary conditions. Bulk simulations employ periodic boundary conditions in all spatial dimensions. For flow in a channel geometry, we employ periodic boundary conditions only in the flow direction, while no-slip boundary conditions on the confining walls are realized by the well-known bounce back rule, where the velocity of the particles are reversed as a result of collision with the wall.30 Since the random grid shifting can lead to underpopulated cells at the walls that would undermine the no-slip condition, additional, so-called ghost particles are added in the collision step.30 So far, we have specified the translational motion only. To simulate ferrofluid flow, we also need to specify the dynamics of the magnetic moments of the particles. From the kinematic equation for rotations, u̇i ≙ ωi × ui , and Eq. (4), a weak first-order scheme can be used, ui (t + ΔtB ) ≙ uPi /∣uPi ∣, where uPi ≙ ui (t) + Δωi × ui (t) with 1 ΔtB 1 Δωi ≙ [τB ΩCi + ui × hCi ] + √ ΔWi . 2 τB τB J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing (9) scitation.org/journal/jcp In Eq. (9), we introduced the dimensionless magnetic field, which is defined as the ratio of the local Zeeman energy over the thermal energy, hCi ≙ μH(rCi )/kB T. Furthermore, ΔWi denotes the Wiener increment over a time interval Δt B and τ B ≙ ζ/(2kB T) denotes the Brownian rotational relaxation time of a single nanoparticle. Equation (9) requires Δt B /τ B ≪ 1. For the simulations shown in the following, we use a second-order stochastic Heun algorithm. More details on the algorithm are given in Appendix A. In Eq. (9), we denote the local vorticity and local magnetic field with ΩCi and hCi , respectively, to indicate that these fields are not evaluated on the precise location of particle i but at the center of the collision cell to which particle i belongs at this time step. Similarly, we define the instantaneous magnetization of the collision cell as MCi ≙ (Msat /NCi )∑j∈𝒞 i uj with M sat ≙ nμ being the saturation magnetization and n being the number density of magnetic nanoparticles. For the mean magnetization, we recover the usual relation ⟨MCi ⟩ ≙ Msat ⟨u⟩. The magnetization dynamics (4) is coupled to the flow field via the vorticity, i.e., gradients of the velocity field, Ω ≙ (1/2)∇ × v, whereas gradients in the magnetization or magnetic field lead to backflow effects via the magnetic force density in the momentum balance, Eq. (2). In order to compute these gradients, we evaluate the velocity and magnetization field in each time step via a kernel smoothing method and determine spatial gradients from a finite-difference approximation. Further details on the numerical procedure are provided in Appendix B. IV. SIMULATION SETUP Consider a simple channel geometry of width L with infinite plates perpendicular to the y-direction, as sketched in Fig. 1. The externally applied field is denoted by H0 . We assume the exterior to be non-magnetic such that M0 ≙ 0 with the external magnetic induction B0 ≙ μ0 H0 . The magnetic field inside the channel is denoted by H. The presence of ferrofluid gives rise to a magnetization M such that the magnetic induction inside the channel is given by B ≙ μ0 (H + M). If k̂ denotes the normal vector on the interface, the continuity conditions of magnetostatics require the normal component of the magnetic induction B and the tangential component of the magnetic field H, which are continuous at the interface,1 k̂ ⋅ (B0 − B) ≙ 0, k̂ × (H0 − H) ≙ 0. (10) The internal field can be expressed as H ≙ H0 − D̄ ⋅ M where the demagnetization tensor D̄ depends on the geometry. For the simple channel geometry of Fig. 1 with infinite plates and a spatially homogeneous external field H0 , we obtain H x ≙ H 0,x and H y ≙ H 0,y − M y . The same result is obtained for an infinite needle-shaped ellipsoidal geometry.31 It is important to note that Maxwell’s equation together with the continuity conditions (10) requires for our case the magnetization M to be independent of the x-position; thus, M ≙ M(y) can depend only on the vertical position within the channel. Further details are provided in Appendix C. In terms of the dimensionless field h, the relation between the external and internal field can be written as hy ≙ h0,y − 3χ L ⟨uy ⟩ with the Langevin susceptibility χ L ≙ nμ2 /(3kB T). 156, 144905-3 The Journal of Chemical Physics ARTICLE and the time step Δt ≙ 1. We note that we limit ourselves to small forcings; otherwise, smaller values of Δt are needed. Having set the basic MPC units, Table I gives an overview of the remaining model parameters and the values or range of values used in subsequent simulations. A number of comments are in order. Temperature is 2 measured in units of Tref ≙ mvmax /kB with vmax ≙ a/Δt being the maximum propagation speed. The mean free path is defined as √ ∗ λ ≙ Δt kB T/m. We choose T ≙ T/T ref < 1 to ensure the mean free path is smaller than the grid size. The mean number of MPC particles per collision cell Q is chosen between 20 and 100. Furthermore, we choose τB∗ ≙ 100 to model slower relaxation of nanoparticles compared to the base fluid. This then allows us to set Δt B ≙ Δt, satisfying Δt B /τ B ≙ 10−2 ≪ 1. Since every MPC particle represents a small volume element of fluid, there is no contradiction of choosing n∗ independent of Q. Although particles in the original MPC scheme are point-like, introducing the rotational friction coefficient ζ in Eq. (4) allows us to associate an effective hydrodynamic diameter σ to the particles via Stokes’ formula ζ ≙ πηs σ 3 , where ηs denotes the viscosity of the solvent. Here, we refer to the pure MPC fluid as “solvent.” Therefore, it is possible to express the density ratio n∗ in terms of the volume fraction of magnetic particles ϕ as n∗ ≙ 3ϕνs∗ /τB∗ , where we used the definition of the reference pressure pref and the dimensionless kinetic viscosity νs∗ ≙ νs Δt/a2 with νs ≙ ηs /ρ. Therefore, even though every particle in the simulation carries a magnetic moment μ, the condition n∗ τB∗ /νs∗ ≪ 1 still corresponds to dilute suspensions of magnetic nanoparticles, which, strictly speaking, is needed for justifying the magnetization dynamics described by Eq. (4). FIG. 1. (a) Schematic of the channel geometry and velocity profile. (b) Channel including square cylinder in Sec. VI. V. RESULTS FOR CHANNEL FLOW In order to study Poiseuille flow, we model an applied pressure gradient as a homogeneous external forcing in the x-direction, fext ≙ f ext î. This constant force is applied to all particles; see Eq. (5). For the channel geometry of Fig. 1, the dimensionless magnetic force density f∗M ≙ (Δt 2 /a)fM from Eq. (2) becomes All simulations are started from an equilibrium initial state, where random velocities are assigned to the particles according to the Maxwell–Boltzmann distribution at temperature T. Initially, particles are placed at random positions to uniformly fill the simulation box, and the orientations ui of the particles’ magnetic moments are chosen randomly from the three-dimensional isotropic distribution. Typically, simulations are run for 105 integration steps, and averages are collected only in the second half of the run where the steady-state has been reached. For n∗ ≙ 0, backflow effects are absent and the standard translational dynamics of the MPC fluid is recovered. We tested our numerical implementation for this case against literature data, dropping angular momentum conservation and choosing a fixed rotation angle α. For equilibrium bulk simulations with periodic boundary conditions in both spatial dimensions, we recover the theoretical ⎛ 1 [(3χL ⟨uy ⟩ − h0,y ) d⟨ux ⟩ + (h0,x + 3χL ⟨ux ⟩) d⟨uy ⟩ ]⎞ dy∗ dy∗ ⎟ ⎟, ≙ −n ⎜ ⎜ ⎟ d⟨uy ⟩ 3χ ⟨u ⟩ L y ⎝ ⎠ dy∗ (11) where n∗ ≙ nkB T/pref is a density ratio and pref ≙ ρa2 /Δt 2 is a reference pressure. In Eq. (11), we used the fact that the magnetostatic fields are independent of the x-coordinate. We follow common practice in MPC simulations and choose particle masses as mi ≙ 1, the linear size of the collision cell a ≙ 1, f∗M scitation.org/journal/jcp ∗⎜ 2 TABLE I. Overview of model parameters and their numerical values used in these simulations: mean number of MPC particles ∗ per collision cell Q ≙ ⟨NCi ⟩, reduced temperature T ≙ T/T ref , dimensionless Brownian relaxation time τB∗ ≙ τB /Δt, the Langevin susceptibility χ L , the Langevin parameter h ≙ μH/k B T, the density ratio n∗ ≙ nk B T/pref , the strength of the external forcing f ∗ext , and the maximum velocity v max . ∗ Variable Values Q T τB∗ χL h n∗ f ∗ext vmax 20–100 0.05–0.4 100 0 0–5 0–0.005 0–5 × 10−5 0–0.1 J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing 156, 144905-4 The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp result for the self-diffusion coefficient, D∗ ≙ T ∗ Δt(1/2 + b/∥1 − b∥) with b ≙ 1/Q + (1 − 1/Q) cos α, where Q ≙ ρa2 /m is the average number of particles in a collision cell.18 For this test, we chose Q ≙ 20 and varied α between 70 and 120○ . As a second test, we reproduced the channel flow simulations presented in Ref. 32 for Q ≙ 35, α ≙ 90○ and T ∗ ≙ 0.4. We verified that the no-slip condition is satisfied. Moreover, from a fit of the velocity profile to the theoretical result for Poiseuille flow, vx (y) ≙ f ext y(L − y), 2ν (12) we extract the fluid viscosity ν∗ ≈ 0.089 ± 0.001 in agreement with the findings reported in Ref. 32. Having tested our numerical implementation in the uncoupled case, we now proceed to investigate the fully coupled system including angular momentum conservation and backflow effects (n∗ > 0). For the channel geometry shown in Fig. 1, we choose the channel width L ≙ 32 or L ≙ 64, length 50, and Q ≙ 100 particles per collision cell, i.e., in total 1.6 or 3.2 × 105 particles. For T ∗ ≙ 0.1 and no externally applied field, h0 ≙ 0, we find a parabolic velocity profile from which we extract the kinematic viscosity of the base fluid from Eq. (12) as νs∗ ≙ 0.114 ± 0.001. We refer to this as the solvent viscosity since the magnetic force density (2) vanishes in this case. For the same conditions but in the quiescent state, we extract a self-diffusion coefficient D ≈ 0.05 in agreement with the Q ≫ 1 limit of the theoretical expression, D∗ → T ∗ Δt/2. Therefore, the Schmidt number measuring the ratio of viscous over molecular diffusion becomes Sc ≙ νs /D ≈ 2.3. The condition Sc > 1 indicates that collisional transport dominates over kinetic transport, which is the relevant regime for fluid simulations.33 Particle-based flow simulations, in general, and MPC, in particular, do not strictly observe the incompressibility condition ∇ ⋅ v ≙ 0. For small Mach numbers, Ma ≙ vmax /cs ≪ 1, incompressibility is approximately restored. √ Here, vmax is the maximal fluid velocity and cs ≙ 5kB T/(3m) is the speed of sound. For the present conditions, T ∗ ≙ 0.1, Q ≙ 100, and f ext ≙ 5 × 10−5 , we find vmax ≈ 0.056, leading to relatively small Mach numbers, Ma ≈ 0.14. Therefore, we can consider the fluid to be approximately incompressible. The Reynolds number is defined as Re ≙ UL/νs , where U denotes a characteristic flow velocity. For U ≙ vmax , we obtain Re ≈ 15.7, well in the laminar regime for channel flow. The Weissenberg number Wi ≙ τ B U/L quantifies the ratio of elastic to viscous forces. Choosing again vmax as the characteristic flow velocity U, we find Wi ≈ 0.18, which can be considered to be still in the Newtonian limit.34 Having ensured that the simulations are performed in the proper parameter regime for laminar, viscous, near-incompressible flow, we now investigate the effect of an external magnetic field on the flow behavior. For simplicity, we choose χ L ≙ 0 in agreement with the model (4) so that demagnetization effects are absent and the internal field h coincides with the external field h0 . Other parameters are chosen as τB∗ ≙ 100 and n∗ ≙ 10−3 . Figure 2 shows the flow profile for identical external forcing f ∗ext ≙ 5 × 10−5 but different strengths of the magnetic field h applied in the gradient direction. We observe from Fig. 2 that the maximum flow velocity decreases with increasing magnetic field strength. We find that the Poiseuille profile (12) fits the velocity field very J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing FIG. 2. Velocity profiles for f ∗ext ≙ 5 × 10−5 , n∗ ≙ 10−3 , χ L ≙ 0, and various strengths of the magnetic field h. Solid lines show fits to the parabolic profile, Eq. (12). accurately. From these fits, we extract a field-dependent effective kinematic viscosity ν(h). For h ≙ 0, we obtain the solvent viscosity, ν(0) ≙ νs , as discussed above. For h > 0, we find that the effective viscosity increases with increasing magnetic field strength. This phenomenon was pioneered by McTague3 and is known in the literature as the “magnetoviscous effect.”5,8 Figure 3 shows the relative viscosity increase Δν(h)/ν(0), where Δν(h) ≙ ν(h) − ν(0) denotes the viscosity change and ν(h) FIG. 3. Dependence of the relative viscosity Δν/ν(0) on the dimensionless magnetic field strength h for n∗ ≙ 10−3 and n∗ ≙ 2 × 10−3 . Other parameters are chosen as f ∗ext ≙ 5 × 10−5 and χ L ≙ 0. Symbols show simulation results, whereas solid lines correspond to the theoretical result, Eq. (13). 156, 144905-5 The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp is obtained from fits to the velocity profile (12) as described above. We verified that we obtain identical results within numerical accuracy for weaker forcing, f ∗ext ≙ 2 × 10−5 . For weak flows, Wi ≪ 1, the model (4) can be solved analytically to give Δν(h) 3 hL2 (h) ≙ ϕ , ν(0) 2 h − L(h) (13) where L(h) ≙ coth(h) − 1/h denotes the Langevin function.4,27 As shown in Fig. 3, the numerical simulations are in very good agreement with this theoretical result. In addition, the prefactor (3/2)ϕ, giving the maximum of the relative viscosity increase, determined from fits to the numerical data (ϕ ≈ 0.27 ± 0.01 for n∗ ≙ 10−3 ) is in good agreement with the relation ϕ ≙ n∗ τB∗ /(3νs∗ ) obtained above (giving here ϕ ≈ 0.29). Before commenting on the slight deviation in the value of ϕ, we want to mention that ferrofluids are typically dilute, corresponding to smaller values of ϕ. For the purpose of demonstrating the method, however, the current choice of parameters should be sufficient to validate the correct implementation of the hybrid model combining MPC with the stochastic magnetization dynamics. Upon closer inspection, we find from the velocity profiles that the no-slip boundary condition and therefore the parabolic velocity profile is not perfectly satisfied. This effect is known in the literature, and methods have been proposed to deal with this effect.30 To investigate these deviations in more detail, we show in Fig. 4(a) the profile of the flow vorticity Ωz ≙ (1/2)∥∇ × v∥z ≙ −(1/2)dvx (y)/dy. The same conditions as used for Fig. 2 are chosen. We scale Ωz with the maximum vorticity expected for Poiseuille flow in the field-free case, Ωmax ≙ f ext L/(4νs ). While the profiles are nicely linear in the center of the channel as expected for Poiseuille flow, deviations near the channel walls are apparent. As expected, these deviations are more pronounced in the vorticity compared to the velocity profile due to spatial gradients. We performed additional simulations for wider channels (L ≙ 64) and verified that the perturbation of the profile does not grow with L but remains a boundary effect, propagating only a distance around 2a into the channel. Figure 4(b) shows the perpendicular magnetization profile M – /M sat ≙ ⟨uy ⟩ across the channel. As expected, M – vanishes in the center of the channel since the vorticity is zero there. Overall, we find a linear profile of M – that follows the linear profile of the vorticity of the Poiseuille flow shown in Fig. 4(a). However, deviations from the linear profile are apparent near the channel walls, which result from deviations in the vorticity near the walls observed in Fig. 4(a). It is worth noting that deviations from the linear magnetization profile occur over a slightly broader boundary layer compared to the deviations from the Poiseuille velocity profile. Discussing these boundary effects in more detail is beyond the scope of the present work. Here, we just want to mention that the reduced vorticity near the wall leads to a corresponding reduction in the local perpendicular magnetization component. As a consequence, the overall viscosity change Δν is slightly reduced compared to the theoretical value, as seen in the slightly reduced value of ϕ found above. J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing FIG. 4. (a) The profile of the flow vorticity Ωz /Ωmax is shown across the channel under the same conditions as in Fig. 2. The position y within the channel is scaled with the channel width L and the maximum vorticity Ωmax for h ≙ 0. Solid symbols correspond to L ≙ 32, and open symbols correspond to L ≙ 64. Solid lines show a linear fit. (b) The perpendicular magnetization component M – normalized by the saturation magnetization M sat is shown across the channel with y under the same conditions as in Fig. 2. VI. RESULTS FOR FLOW AROUND CYLINDER As a further demonstration of the flexibility of the MPC method, we here consider the two-dimensional flow of a ferrofluid around a square cylinder of diameter D inside a planar channel of width L and length Lx . No-slip boundary conditions are imposed on the channel and cylinder walls by the same bounce back algorithm as used in Sec. V. Following Ref. 32, we impose flow by prescribing the average velocity vx (y) ≙ 4vLmax 2 y(L − y) for 156, 144905-6 The Journal of Chemical Physics particles within the inlet region 0 ≤ x < 10. In addition, periodic boundary in the flow direction are imposed. The horizontal position of the cylinder center is chosen as Lx /4, which was found to sufficiently reduce the influence of inflow and outflow boundary conditions.35 Following previous studies,32,35 we choose the blockage ratio D/L ≙ 1/8. Determining the flow around a cylinder is a classical problem in fluid dynamics where different flow regimes can be distinguished according to the value of the Reynolds number Re ≙ vmax D/ν. The creeping flow regime for Re < 1 is dominated by viscous forces, and no separation is observed. For larger Reynolds numbers, the flow field separates at the downstream side of the cylinder, forming two steady, counter-rotating vortices behind the cylinder. The size of this recirculation region Lr increases with Re until the onset of the van-Kármán vortex street at a critical Reynolds number.35 We here focus on the regime 5 ≤ Re ≤ 30 where we expect to see steady recirculation regimes. In order to obtain reasonable spatial resolution, we choose L ≙ 200 and Lx ≙ 600. For the average number of particles per collision cell, we set Q ≙ 30, which means the simulations contain around 3.6 × 106 particles. The temperature is chosen as T ∗ ≙ 0.05, corresponding to a viscosity FIG. 5. Stationary velocity field for the flow around a square cylinder (gray) for h ≙ 0 (top) and h ≙ 4 (bottom). For better visibility, only a part of the velocity field is shown. The simulation parameters are chosen as T ∗ ≙ 0.05, n∗ ≙ 5 × 10−3 , χ L ≙ 0, and v max ≙ 0.1. J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing ARTICLE scitation.org/journal/jcp νs∗ ≈ 0.0913 ± 0.0001. Furthermore, to observe magnetic fieldinduced changes in the velocity field more clearly, we extend the range of concentrations to n∗ ≙ 5 × 10−3 . Choosing again χ L ≙ 0 eliminates contributions of the demagnetization field. Simulations are performed for a total of 105 integration steps and averages extracted after 5 × 104 steps. We calculate two-dimensional velocity and magnetization fields and their derivatives using kernel smoothing methods as described in Appendix B. Figure 5 shows the stationary velocity fields with and without the applied magnetic field for a particular choice of model parameters (T ∗ ≙ 0.05, n∗ ≙ 5 × 10−3 , vmax ≙ 0.1). We observe that the applied magnetic field changes the flow field and causes a smaller recirculation region. Qualitatively, such a change is expected due to the field-induced increase in viscosity. Looking at a larger portion of the velocity field, Fig. 6 shows the significant changes the external magnetic field induces not only on the wake but also on the larger-scale flow field. FIG. 6. The color-coded magnitude of the stationary velocity field is shown (h ≙ 0 in the top and h ≙ 4 in the bottom panel) under the same conditions as in Fig. 5 but for a larger part of the flow field. 156, 144905-7 The Journal of Chemical Physics To make these observations more quantitative, we extract the length of the recirculation zone, Lr , from the velocity fields. In particular, we determine Lr from the zero-crossing of the x-component of the centerline velocity at the end of the wake. We use a parabolic fit to the centerline velocity profile near the end of the wake to determine Lr and estimate error bars based on uncertainties in the fit parameters. Figure 7 shows the length of the recirculation region Lr in units of the diameter D of the cylinder as a function of the magnetic field and the Reynolds number. Increasing the magnetic field strength leads to a decrease in Lr , the effect being more pronounced for larger concentrations. This observation is consistent with the increased effective viscosity we found in Sec. V. In the absence of FIG. 7. (a) Length of the recirculation region Lr scaled with the diameter of the cylinder D as a function of the magnetic field strength h. Diamonds and circles correspond to n∗ ≙ 2 × 10−3 and 5 × 10−3 , respectively. Other parameters are chosen as T ∗ ≙ 0.05 and v max ≙ 0.1. (b) The data in panel (a) are shown as a function of the effective Reynolds number. In addition, open black squares show the result for h ≙ 0 and n∗ ≙ 5 × 10−3 but varying v max . J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing ARTICLE scitation.org/journal/jcp an external magnetic field, by varying the maximum inflow velocity vmax , we recover earlier results showing a linear increase in the length of the recirculation zone with the Reynolds number, Lr /D ≙ c0 + c1 Re, (14) with c1 ≈ 0.059 ± 0.007, indicated as the black dashed line in Fig. 7(b). Keeping instead vmax ≙ 0.1 fixed but varying the magnetic field strength h, we also find the linear relationship Eq. (14) when plotting the data in Fig. 7(a) against the effective Reynolds number Re ≙ vmax D/ν(h) with ν(h) being the field-enhanced effective kinematic viscosity discussed in Sec. V. Interestingly, however, we observe that Lr varies much less strongly with Re due to a magnetic field compared to the non-magnetic case. Indeed, we find there a slope c1 ≈ 0.021 ± 0.003, less than half the value in the field-free case when varying vmax . To further illustrate the effect of a magnetic field on ferrofluid flow, we also study the drag coefficient Cd . The drag coefficient 2 is defined in two dimensions as Cd ≙ 2F∥ /(Qvmax D), where F ∥ denotes the force in the flow direction exerted by the fluid on the cylinder.32,35 Within the MPC scheme, the force F ∥ can readily be evaluated from the parallel component of the total momentum the particles transfer to the cylinder wall when bouncing back. Figure 8 shows the drag coefficient Cd as a function of the Reynolds number. In the absence of an externally applied field (black open symbols), we observe the characteristic strong decrease in Cd with Re, approximately as an inverse power law, Cd ∼ Re−z with z ≈ 1.35 ± 0.15, in agreement with finite-volume simulations in Ref. 35. However, when a magnetic field is applied (filled colored symbols), we observe that the drag coefficient is significantly smaller compared to the non-magnetic fluid at the same Reynolds number. We observe this magnetic drag reduction in the low-Reynolds number regime 5 ≲ Re ≲ 20. FIG. 8. The dimensionless drag coefficient Cd is shown as a function of the Reynolds number under the same conditions and color coding as in Fig. 7. Stars indicate the results of Ref. 35, and the dashed line indicates a power-law fit. 156, 144905-8 The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp VII. CONCLUSIONS AUTHOR DECLARATIONS Here, we present an extension of the MPC method to simulate dynamics and flow of magnetic fluids, including fluctuation and backflow effects. In analogy to recent extensions of MPC to nematic liquid crystals,21,22 we equip each MPC particle with a magnetic moment and include stochastic magnetization dynamics via additional rotational motion of the individual particles. Fluid and magnetization dynamics are coupled to each other via velocity gradients and the magnetic force density. We successfully tested this hybrid scheme in a standard, two-dimensional channel geometry. For Poiseuille flow, we reproduce the magnetoviscous effect and recover quantitatively the theoretical result for the relative viscosity increase with magnetic field strength. Using standard methods to implement no-slip conditions in MPC simulations, we nevertheless observe small deviations of the velocity gradient from the theoretical profile very close to the wall, leading to corresponding deviations in the magnetization profile. This effect is already known in the literature, and more refined methods have been proposed to better realize the no-slip condition.30 We leave this more technical aspect for future research. We also illustrate the new MPC method for magnetic fluids for the benchmark problem of flow around a square cylinder. In the absence of an applied field, we verify known results for the length of the recirculation region behind the cylinder increasing linearly with the Reynolds number and the decrease in the drag coefficient for non-magnetic fluids. In addition, we study the dependence of these quantities on an externally magnetic field for ferrofluids. We observe that the length of the wake is increased and the drag is reduced by a magnetic field compared to the non-magnetic case at the same effective Reynolds number. Being able to reduce drag by a magnetic field in the low-Reynolds number regime might be of interest in microfluidics applications. More generally, we find that these quantities are not described by the effective Reynolds number alone even when accounting for the fieldinduced effective viscosity. This is likely due to the anisotropy of the effective viscosity induced by the magnetic field,3,36 which leads to ferrofluids showing different flow behaviors than ordinary viscous fluids. The method proposed here to simulate ferrofluid flow naturally includes thermal fluctuations that are particularly relevant at small scales. Furthermore, as a solver for fluctuating ferrohydrodynamics, the method inherits all the benefits of the MPC approach, which can easily be extended to three spatial dimension and more complicated geometries. The MPC method is also particularly well-suited as an efficient way to model solvent effects in colloidal suspensions.37 Future studies may also include the effect of the demagnetization field that we have neglected here. In addition, the method is very flexible and can straightforwardly accommodate different magnetization equations, such as, e.g., the so-called chain model38 to describe chain-forming ferrofluids or including mean–field interactions.8 Therefore, we expect this method to be a useful tool for the simulation of flow and hydrodynamic effects in magnetic fluids. Conflict of Interest ACKNOWLEDGMENTS Discussions with Anoop Varghese during a very early stage of this project are acknowledged. J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing The authors have no conflicts to disclose. DATA AVAILABILITY The data that support the findings of this study are available from the corresponding author upon reasonable request. APPENDIX A: STOCHASTIC HEUN ALGORITHM FOR ORIENTATIONS The stochastic Heun algorithm is a predictor–corrector scheme with second-order weak convergence in the Stratonovich sense.39 For rotational motion, we need to ensure that the orientations ui (t) remain unit vectors for all times. In the predictor step, new orientations ŭi are calculated from an Euler scheme as indicated in Sec. III, ŭi ≙ uPi , ∣uPi ∣ uPi ≙ ui (t) + Δωi (t) × ui (t), (A1) (A2) and Δωi (t) is given by Eq. (9) evaluated at time t. These predicted orientations ŭi are now used to computed new angular velocity changes Δω̆i , where for the latter, the right-hand side of Eq. (9) is evaluated with uPi instead of ui (t). In the corrector step, the new orientations are obtained from ui (t + Δt) ≙ uCi , ∣uCi ∣ 1 uCi ≙ ui (t) + ∥Δωi (t) × ui (t) + Δω̆i × ŭi ∥. 2 (A3) (A4) APPENDIX B: KERNEL SMOOTHING Since the evaluation of velocity and magnetostatic fields is important for the model, we give here some details on their evaluation in the simulation. We exemplify the method for the velocity component in the flow direction, vx . The other fields are evaluated in the same manner. First, we use the Nadaraya–Watson kernel regression estimator40 to find the instantaneous field as ∑i≙1 vx,i (t)K( ∣ri (t)−r∣ ) b , ∣rj (t)−r∣ N ∑j≙1 K( b ) N vx (r; t) ≙ (B1) where K(z) is known as the kernel and the parameter b is known as the bandwidth or smoothing length. While the uniform kernel is frequently used, we here employ the Epanechnikov kernel, ⎧ 3 2 ⎪ ⎪ ⎪ (1 − z ), ∣z∣ ≤ 1, K(z) ≙ ⎨ 4 ⎪ ⎪ 0, else, ⎪ ⎩ (B2) which is known to minimize the mean integrated square error. We found that choosing b equal to the linear size of the collision cells a 156, 144905-9 The Journal of Chemical Physics ARTICLE provides a good compromise between smoothing and keeping local fluctuations. Since there is a systematic attenuation bias due to smoothing, we resort to finite-difference schemes to calculate spatial gradients. Having evaluated the instantaneous velocity field vx at the centers of the collision cells from Eq. (B1), we use the central difference scheme to approximate the spatial partial derivatives as ∂vx (r; t) vx (r + ak̂; t) − vx (r − ak̂; t) ≈ ∂y 2a (B3) with the unit vector k̂ defined in Fig. 1 and correspondingly with the unit vector î in the x-direction for the partial derivative with respect to x. Only near the channel walls, we use a first-order approximation instead. We found that the central difference scheme provides good results also in comparison to using higher-order, differentiable kernel functions. APPENDIX C: MAGNETOSTATICS IN CHANNEL GEOMETRY For non-conducting fluids, Maxwell’s equations are given by Eq. (3) and a separate, medium-dependent magnetization equation. First, we consider the exterior of the system to be non-magnetic, M0 ≙ 0, B0 ≙ μ0 H0 , (C1) and we assume the fields H0 (and consequently B0 ) are spatially uniform. Thus, Maxwell’s equations (3) are identically satisfied in the exterior, ∇ × H0 ≙ 0 and ∇ ⋅ B0 ≙ 0. Let us denote the fields inside the fluid as H and B with B ≙ μ0 (H + M). (C2) Define the contribution Hff of the magnetic fluid to the internal field, H ≙ H0 + Hff . (C3) Note that Hff is often denoted in terms of a “demagnetization field,” H ≙ H0 − D ⋅ M with D being the demagnetization tensor, which depends on the shape of the sample. Note that D is symmetric, has trace one, and all diagonal elements are non-negative.41 We now specialize to the channel geometry sketched in Fig. 1. From the continuity condition on the magnetic field, Eq. (10), we find that By ≙ B0,y , Hx ≙ H0,x ⇒ Hff,x ≙ 0. (C4) (C5) Thus, we know that Hff ≙ Hff k̂ is oriented normal to the wall. Next, to satisfy Maxwell’s equation ∇ × H ≙ 0, we require ∇ × Hff ≙ 0 and therefore conclude that H ff ≙ H ff (y) must be independent of the horizontal position x along the channel. Similarly, we can conclude from ∇ ⋅ B ≙ 0 and the continuity condition that ∂Bx /∂x ≙ 0 and therefore ∂M x /∂x ≙ 0 since By is constant. Therefore, the magnetization component M x is also independent of the horizontal position along the channel. Finally, we can also write ∇ ⋅ B ≙ 0 as ∇ ⋅ (H + M) ≙ 0, which leads to ∂(H ff + M y )/∂y ≙ 0. Since we already concluded that H ff is independent of x, so must J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing scitation.org/journal/jcp M y to satisfy this condition for all points within the channel. Therefore, none of the magnetostatic fields depends on the x-position within the channel. REFERENCES 1 R. E. Rosensweig, Ferrohydrodynamics (Cambridge University Press, Cambridge, 1985). 2 Colloidal Magnetic Fluids, Lecture Notes in Physics Vol. 763, edited by S. Odenbach (Springer, Berlin, 2009). 3 J. P. McTague, J. Chem. Phys. 51, 133 (1969). 4 M. A. Martsenyuk, Y. L. Raikher, and M. I. Shliomis, Sov. Phys. J. Exp. Theor. Phys. 38, 413 (1974). 5 S. Odenbach, Magnetoviscous Effects in Ferrofluids, Lecture Notes in Physics Vol. 71 (Springer, Berlin, 2002). 6 M. Kröger, P. Ilg, and S. Hess, J. Phys.: Condens. Matter 15, S1403 (2003). 7 M. Colombo, S. Carregal-Romero, M. F. Casula, L. Gutiérrez, M. P. Morales, I. B. Böhm, J. T. Heverhagen, D. Prosperi, and W. J. Parak, Chem. Soc. Rev. 41, 4306 (2012). 8 P. Ilg and S. Odenbach, in Colloidal Magnetic Fluids: Basics, Development and Applications of Ferrofluids, Lecture Notes in Physics Vol. 763, edited by S. Odenbach (Springer, Berlin, 2008). 9 L. J. Felicia, S. Vinod, and J. Philip, J. Nanofluids 5, 1 (2016). 10 C. Rinaldi and M. Zahn, Phys. Fluids 14, 2847 (2002). 11 K. R. Schumacher, I. Sellien, G. S. Knoke, T. Cader, and B. A. Finlayson, Phys. Rev. E 67, 026308 (2003). 12 M. Hirabayashi, Y. Chen, and H. Ohashi, Phys. Rev. Lett. 87, 178301 (2001). 13 M. Sheikholeslami, M. Barzegar Gerdroodbary, S. V. Mousavi, D. D. Ganji, and R. Moradi, J. Magn. Magn. Mater. 460, 302 (2018). 14 L. Huang, T. Hädrich, and D. L. Michels, ACM Trans. Graphics 38, 93 (2019). 15 S. Kayal, D. Bandyopadhyay, T. K. Mandal, and R. V. Ramanujan, RSC Adv. 1, 238 (2011). 16 A. Munaz, M. J. A. Shiddiky, and N.-T. Nguyen, Biomicrofluidics 12, 031501 (2018). 17 A. Malevanets and R. Kapral, J. Chem. Phys. 110, 8605 (1999). 18 T. Ihle and D. M. Kroll, Phys. Rev. E 63, 020201 (2001). 19 G. Gompper, T. Ihle, D. M. Kroll, and R. G. Winkler, Advanced Computer Simulation Approaches For Soft Matter Sciences III (Springer, Berlin, Heidelberg,, 2009), Vol. 221, p. 1. 20 T. N. Shendruk and J. M. Yeomans, Soft Matter 11, 5101 (2015). 21 K.-W. Lee and M. G. Mazza, J. Chem. Phys. 142, 164110 (2015). 22 S. Mandal and M. G. Mazza, Phys. Rev. E 99, 063319 (2019). 23 R. E. Rosensweig, in Ferrofluids: Magnetically Controllable Fluids and Their Applications, Lecture Notes in Physics Vol. 594, edited by S. Odenbach (Springer, Berlin, 2002), pp. 61–84. 24 M. I. Shliomis, in Ferrofluids: Magnetically Controllable Fluids and Their Applications, Lecture Notes in Physics Vol. 594, edited by S. Odenbach (Springer, Berlin, 2002), pp. 85–111. 25 A. Leschhorn and M. Lücke, Z. Phys. Chem. 220, 219 (2006). 26 M. Liu and K. Stierstadt, in Colloidal Magnetic Fluids: Basics, Development and Applications of Ferrofluids, Lecture Notes in Physics Vol. 763, edited by S. Odenbach (Springer, Berlin, 2008). 27 P. Ilg, M. Kröger, and S. Hess, J. Chem. Phys. 116, 9078 (2002). 28 D. Soto-Aquino, D. Rosso, and C. Rinaldi, Phys. Rev. E 84, 056306 (2011). 29 B. Kowalik and R. G. Winkler, J. Chem. Phys. 138, 104903 (2013). 30 J. K. Whitmer and E. Luijten, J. Phys.: Condens. Matter 22, 104106 (2010). 31 J. A. Osborn, Phys. Rev. 67, 351 (1945). 32 A. Lamura, G. Gompper, T. Ihle, and D. M. Kroll, Europhys. Lett. 56, 319 (2001). 33 M. Ripoll, K. Mussawisade, R. Winkler, and G. Gompper, Phys. Rev. E 72, 016701 (2005). 34 P. Ilg, M. Kröger, and S. Hess, Phys. Rev. E 71, 031205 (2005). 156, 144905-10 The Journal of Chemical Physics 35 M. Breuer, J. Bernsdorf, T. Zeiser, and F. Durst, Int. J. Heat Fluid Flow 21, 186 (2000). 36 A. Sreekumari and P. Ilg, Phys. Rev. E 92, 012306 (2015). 37 A. Malevanets and R. Kapral, J. Chem. Phys. 112, 7260 (2000). 38 A. Y. Zubarev and L. Y. Iskakova, Phys. Rev. E 61, 5415 (2000). J. Chem. Phys. 156, 144905 (2022); doi: 10.1063/5.0087981 Published under a nonexclusive license by AIP Publishing ARTICLE scitation.org/journal/jcp 39 J. L. Garcia-Palacios, On the Statics and Dynamics of Magneto-Anisotropic Nanoparticles (John Wiley & Sons, 2000). 40 C. Loader, Local Regression and Likelihood, Statistics and Computing (Springer, 1999). 41 R. Moskowitz and E. Della Torre, IEEE Trans. Magn. 2, 739 (1966). 156, 144905-11