1. Introduction

Ferrofluids are colloidal suspensions of permanently magnetized particles of nanometer size dispersed in a solvent. Their rheological properties such as viscosities, storage of energy and dissipative moduli can be controlled with external magnetic fields ¹, which makes them valuable in diverse technological applications. Much of the research on these magnetic fluids relates to the class of magnetorheological suspensions which develop bulk magnetization under external magnetic field conditions ². External fields produce chain-like aggregates of the colloidal particles, an effect that can also be developed in ferrofluids in the absence of fields for sufficiently high particle’s dipolar moment or concentration. It has been observed that when the fluid is subjected to an applied increasing shear rate they change from their Non-Newtonian (nonlinear relationship of stress with the applied shear flow, which it can be produced by thermal fluctuations too) to Newtonian (the stress response of the liquid is linear through a kernel with the shear rate. With the kernel been the stress relaxation modulus) under magnetic fields. Such rheological behavior has been demonstrated to be anisotropic through the use of macroscopic rheology experiments ³, and also more recently through a no invasive technique, the passive microrheology on a probe particle diffusing in the colloidal suspension ⁴ in thermal equilibrium. In their development of this latter experimental technique by Mason et al⁵, they express “bulk relaxations have the same spectrum as the microscopic stress relaxations affecting the suspended particle motion". They assumed that the complex shear moduli can be determined from the particle’s mean squared displacement in the linear regime of low shear flow. This technique has a significant sensitivity and detects high frequencies of relaxations as compared to macroscopic rheology. Ferrofluids look dark with visible light because of its strong absorption by the ferro-colloids. Mertelj et al.⁴ were able to measure the translational diffusion of the probe particle with their use of low-intensity laser beams in dynamic light scattering experiments. In their research, they used small samples which together a careful selection of wavelengths allowed them to reduce the absorption of light by the ferrocolloids. Whereas Yendeti et al⁶ made their microrheology experiments observing a fluorescence silica colloidal particles suspended in the ferrofluid. This way they were able to use video microscopy to track the silica particle positions driven by its Brownian motion inside the ferrofluid. Thus, fine tuning of the wavelength to the order of the probe particle size allows the determination of its mean squared displacement. Therefore, the scattered light originates only from the probe silica colloidal particle since the molecules of the solvent are much smaller than the wavelength of the used light. Solvent provides the reservoir that dissipates energy from the translational and rotational movement imparted by the Brownian motion of all the colloidal particles. It was shown in Refs. ⁷ and ⁸ that dynamic light scattering and video microscopy in microrheology experiments allows the determination of the orientation variable, thus, providing the rotational diffusion of the probe. From its angular displacements, the time-dependent mean squared angular displacement was measured. Consequently, the microrheological shear modulus of the host suspension was obtained. Comprehensive theoretical knowledge of the main measured rheological properties have been reached with thermodynamic ^{9,10,11,12,13,14,15}, and kinetic theories of liquid crystal dynamics ¹⁶ that involve macroscopic constitutive magnetization’s relaxation mechanisms, and both mean field and microscopic free energies of the thermodynamic state, respectively. Thus, the experimentally observed dependence of viscosity with magnetic field was accurately confirmed with a kinetic approach from low up to moderate concentrations before the formation of chain structures by the particles ¹⁷. After the formation of chains, a so-called kinetic chain theory of rod-like aggregates ^18,19 confirms the observations on the increase of viscosity with field quantitatively. In a similar set of macroscopic experiments at low strain in a system where the finite size of the vessel containing the liquid was important ²⁰, it was determined the energy dissipated due to aggregate’s structural deformation when they are partially bound to the confining walls. It was found that the aggregates experience cohesion and movement of their inner constituting particles which generates heat through this process. Such a heating mechanism is relevant in hyperthermia applications where ferro-particles diffuse, rotate and stick to the boundaries in the multicomponent medium inside of a tissue ²¹. These findings have raised the interest to investigate the rheology of complex fluids of mixtures of added polymer solutions to the magnetorheological fluid ^20,22,23,24. Thus, the theoretical description of their observed rheology and effect of size polydispersity ^25,26,27,28, have been described with polymers motivated models ^29,30. More recently, theories based on macroscopic thermodynamics predict as in polymer theory elongational viscosity and its consequence on the ratios of the viscoelastic moduli ³⁰, and of the module to viscous dissipation time. Such predictions were adopted in the microrheology experiments of Refs. ^4,6 in order to get estimates for the values of all anisotropic properties. In the present manuscript we propose expressions for the complex shear modulus response for the rheology of an equilibrium ferrofluid derived from a statistical mechanics approach valid in the regime of linear shear flow of microrheology experiments. The viscoelastic moduli depend on the pair distribution function which captures the bulk microstructure of the ferro-particles due to their pair direct interactions. We show that for a given concentration of particles, ferrofluids at equilibrium state display rheological behavior in a similar way as it is observed to occur in non-equilibrium anisotropic magnetic fluids under external fields and shear flow. For the dispersed homogeneous phase of the colloid, we found that the dissipative modulus is dominant at low frequencies where the colloid displays a liquid character with a crossover to an elastic fluid mode at moderate frequencies where the elastic modulus becomes dominant at higher frequencies. The calculated loss modulus has good agreement with Langevin dynamic simulations up to the transition to the elastic behavior where it starts to deviate. Whereas the general shape of the elastic modulus, as obtained from simulations of the mean square translational displacement of a probe ferro-particle is qualitatively well described by the model proposed. An increase of the dipolar moment per particle or of their concentration can lead to the formation of chain like aggregates in absence of external fields. Our model theory of microrheology is useful in this case too. Since the formation of chains is taken into account properly by the pair correlation function which describes the bulk micro-structural order of the particles due to their direct interactions. We also provide the analytical expressions for these moduli when the observable is the mean square angular rotation for experiments that observe the rotational diffusion of ferro-particles. These viscoelastic properties can also be used to study in the linear regime of flow the rheology of equilibrium isotropic ferroelectric colloids ^31,32,33.

2. Translational diffusion microrheology

Tracking microrheology experiments are useful techniques that accurately provide the linear viscoelasticity of suspensions of polystyrene particles and other complex fluids ^5,7,34,35. In these experiments it is assumed that the complex shear moduli can be determined from the particle’s mean squared displacement in the linear regime of low shear flow. The rheological properties are extracted through the equilibrium Einstein expression of the friction on a probe particle that performs a Brownian motion in the colloid. The friction due to interaction with its neighbors, in the frequency space △ζ(ω), that affects its translational movement is related to the mean square displacement <Δω(ω)2> and to the complex shear modulus by

with i=-1. a the particle radius. ζ0 the solvent friction on a particle. In what follows we show how the friction on the probe particle is obtained.

The Langevin equations of the translational and rotational Brownian motion of the probe particle that interacts with the others in the suspension are the stochastic Eqs. ³⁶

V(t) and W(t) are the translational and angular velocities of the probe. We used a space fixed frame with origin at the particle center of mass, and frame axis following the orientation of the probe main axis of symmetry. M and I are the mass and particle’s matrix of moment of inertia, respectively. We shall not consider neither hydrodynamic interactions among particles nor external magnetic fields. The first two terms in Eqs. (2) are the solvent friction force and torque. The short time free particle diagonal friction tensors ζ0, ζTR0, ζRT0, ζTR0, ζRT0 represent hydrodynamic drag forces, and torques. These friction and random forces are the only quantities that convey information on the nature of the solvent. And they ignore their molecular degrees of freedom of position and orientation coordinates and momenta. They are coupled to the thermally driven solvent random forces f⁰ and torques t⁰ by the fluctuation-dissipation theorems

and k_B Boltzmann constant, T = 300K the absolute temperature. The total force F_TOT and torque T_TOT on the probe by the other colloidal particles which are at the concentration n(r,Ω,t)=∑i=1Nδ(r-ri(t))δ(Ω-Ωi(t)) are given by

with Ω=(θ,φ) and θ,φ being the polar angles, ψ(r,Ω,Ω') is the pair potential, and ∇Ω=n^×∂∂n^ the angular gradient operator. It should be noted that F_TOT (t) and T_TOT (t) are calculated using the probe’s dipole located along the Z axis direction of the its local frame, that is Ω'=(0,0). The unitary cartesian vector n^ of the orientation of any other particle is in the direction of that particle’s axis of symmetry.

Both Langevin equations can be written compactly as

where ∇↔=(∇,r×∇+u×d/du) is the general gradient operator. r denotes the vector joining the centers of the probe and another particle. u is a unitary vector with orientation θ,ϕ. The generalized velocity V⃡=V,W, M⃡ij=Mδij (i,j=1,2,3),M⃡ij=δijIi-3 (i, j = 4,5,6), with I ₁, I ₂, I ₃ being the principal moments of inertia of the tracer. The friction tensor ζ0⃡ is a diagonal matrix which have the nonzero components with elements ζ110=ζ220=ζ⊥0, ζ330=ζ∥0, ζ440=ζ550=ζR0, and ζ660=0 which are external inputs to this theory and are provided by experiment or an external theory. For spherical particles ζ⊥0=ζ∥0=ζ0. Such coefficients describe the particle hydrodynamic friction with the solvent when its diffusion occurs along its dipole orientation (∥), perpendicular (⊥), or performs rotational motion R. Equation (5) can be written to first order in concentration fluctuations δn(r,Ω;t)=n(r,Ω;t)-neq(r,Ω), with the profile distribution of host particles in the probe’s field given as neq(r,Ω)=⟨n(r,Ω,t)⟩ been an equilibrium ensemble average. Thus, the Langevin equation takes on the final form

where neq(r,Ω) does not contribute to the total force and torque. Similarly, a stochastic evolution equation for δn(t) is derived with help of linear irreversible theory of fluctuations ³⁷

with ⟨∇⃡⋅j(r,Ω;t)∇⃡⋅j†(r',Ω';0)⟩≡L(r,r',Ω,Ω';t), and ∇⃡⋅j(t) been a random diffusive flux. This is the more general form of the diffusion equation of the particles. Its solution is

The diffusion relaxation of the host particles around the probe is provided by χ(t) that fullfils

which has initial condition χ(r,r',Ω,Ω';t=0)=δ(r-r')δ(Ω-Ω'). Equation (9) can be converted into a dynamical equation for the Van Hove function of fluctuations in the concentration of particles C(r,r',Ω,Ω';t) with respect to its equilibrium value

which determines the relaxation modes of the cage of particles surrounding the probe. Its initial condition σ=Ct=0=⟨δn(0)δn(0)⟩ is the inhomogeneous static correlation function with inverse given by

Substituting the above solution for δn(t) into the Langevin equation (6) leads to

where

is a fluctuating generalized force arising from the spontaneous departures from zero of the net direct forces exerted by the other particles on the tracer. It groups a random force and torque on the tracer with zero mean value, and time dependent correlation function given by ⟨F⃡(t)F⃡†(0)⟩=kBTΔζ⃡(t), and the time-dependent friction function on the probe (tracer) particle is

Using the Wertheim-Lovett’s relation ³⁸

we derive other useful forms of the friction function Δζ⃡(t) as

with β=1/kBT and C=χ∘σ , ∘≡∫drdΩ, and † is transpose. The angular average in (15) (division by Ω2) is necessary since this is an experimental observable. In Eq. (15) neither hydrodynamic interactions nor external magnetic fields were included. It can also be extended to describe mixtures of species of anisotropic particles with axial symmetry.

The time-dependent memory function Δζ⃡(t) contains the dissipative friction effects derived from the direct interactions of the probe particle with the particles around it. This memory function defines the relaxation time τI≫τB(τB) which is the relaxation time of the momenta of the particles) for the particles to diffuse a mean distance among them. Thus, in the diffusive regime t≫τB, long time overdamped regime means t≫τI. At this time scale the momenta of the colloidal particles have already relaxed and the only relevant dynamic variables are their positions and orientations (r(t),Ω(t)).

This is a general expression for the friction contribution on the tracer due to direct interactions with the particles about it. It depends on the microstructural inhomogeneous total correlation function h(r,Ω)=neq(r,Ω)/ρ-1 of the host suspension of N particles at concentration ρ around the field of the probe, and of the free friction constants ζ0⃡ through χ(t). We now introduce the homogeneity approximation which amounts to ignore the tracer’s field on the properties σ, χ(t) or equivalently on C(t), which then can be determined in the bulk solution. Thus, σ(r,r',Ω,Ω')≈σ(r=∣r-r'∣,Ω,Ω',Ω∣r-r'∣). It is also adopted the Fick's diffusion approximation for χ(t)=exp(-tL∘σ-1), with L(r-r',Ω,Ω',Ωr-r')=ρ[D*0∇2+DR*0∇Ω2]δ(r-r')δ(Ω-Ω'). In this approximation D*0=D0+Dother0, DR*0=DR0+DR,other0, (The case Dγ,other0≪Dγ0, γ=‖, ⊥, R was made in ³⁶). The short-time diffusion coefficients Dother0=D0,DR,other0=DR0 of the other particles are approximated by those of the tracer, and D0=kBTζ0,DR0=kBT/ζR0. For spherical solid particles, ζ0=3πηsold,ηsol is the viscosity of the pure solvent, d is the diameter of the Brownian particle, and ζR0=πηsold3.

Using the above approximation for χ(t) and C=χ∘σ, Eq. (9) can be written as

This last equation governs the diffusive relaxation of C(t), as described from the tracer’s reference frame. In this manner, we have obtained a closed approximate expression for Δζ↔(t) in terms only of the static properties ψ, σ and of the phenomenological quantities D⁰ and DR0.

From (15) the 6x6 diagonal friction matrix, where it is ignored translational and rotational coupling is

with Ω=4π. In this theory the orientation of the probe particle’s main axis of symmetry was defined with respect to its body fixed frame. It has its dipole oriented along the main shaft u = (0, 0, 1). Thus, another particle’s orientation u=u(Ω) is denoted in this frame by the polar angles of its orientation Ωj=(θj,ϕj) for j = 1,…, N.

Thus, the friction matrix (17) has only five non zero diagonal elements; △ζ11(t)=△ζ22(t)=△ζ⊥(t) is the friction on the tracer for perpendicular motion to the dipole or of its anisotropy vector of orientation. △ζ33(t)=△ζ∥(t) describes the parallel motion along the main axis of symmetry. The contribution to the rotational friction is △ζ44(t)=△ζ55(t)=△ζR(t). In Eq. (17) β=1/kBT. The general expressions for the friction function in (17) can be used also to study tracer diffusion on rod-shaped particles suspensions ³⁹ such as fd viruses on which the experimental techniques of birefringence ⁴⁰, and forced Rayleigh scattering measure the diffusion properties ⁴¹. On the other hand, the propagator in (17) C(t):=⟨δn(t)δn(0)⟩ governs the collective relaxation of the particle’s configuration variables (r(t),Ω(t)) due to thermal fluctuations in their local concentration δn(r,Ω,t)=n(r,Ω,t)-neq(r,Ω) about the equilibrium neq=⟨n(r,Ω,t)⟩ of its local instantaneous value n(t)=∑i=1Nδ(r-ri(t))δ(u-ui(t)) ³⁸. It has initial condition C(t=0)=ρS(r,Ω;r',Ω'):=⟨δn(r,Ω;0)δn(r',Ω';0)⟩ where S is the structure factor of particles with numerical density ρ=N/V. The temporal evolution of C(t) is given by Fick’s diffusion law which amounts to an exponential relaxation time and thus corresponds to a Maxwell viscoelastic model.

The colloidal ferrofluid is contained in a volume V. It consist of a carrier fluid which molecular nature is not taken into account, plus the monodisperse system of N spherical particles where each one has a permanent dipolar moment μ.ψij is the direct pairwise interaction potential energy between particles i, j. Because ferro-particles can not overlap their interaction is modeled by a Lennard-Jones (LJ) short range (sr) repulsive, and long-range dipolar (d) potentials ψ12=usr+ϕd

with ϵ the strength, and the dipolar part

μ0 is the magnetic permeability of vacuum. D(Ω1,Ω2,Ωr):=3(r^12⋅u1)(r^12⋅u2)-(u1⋅u2), where r^=r/r a unitary vector with orientation Ωr. Due to the spatial symmetry of the dipolar potential the pair correlation function g(r,Ω1,Ω2) has the fixed space rotational invariant expansion of Wertheim in the form ³⁸

with △:=u1⋅u2.D(Ω1,Ω2,Ωr):=3(r^12⋅u1)(r^12⋅u2)-(u1⋅u2), where r^=r/r a unitary vector with orientation Ωr. u is a unitary vector with orientation u=u(θ,ϕ). Another equivalent laboratory rotational invariant expansion for g(r,Ω1,Ω2) was given by Blum ³⁸

with the identification g000:=gr, h110:=-h△(r)/(3), and h112:=hD(r)(10/3) dictated by the symmetry of the dipolar fluid.

By taking the Fourier-Bessel

⁴² with j_l being the spherical Bessel function, and Laplace transform

of (16) (See Ref. ³⁶) we find

where i=(-1) and

⁴². Using the above values of mnl = 000, 110, 112, then for dipolar liquids α=0,±1, and using the approximation D*0=2D0, DR*0=2DR0, ³⁶ we get for (22)

And the inverse relation holds ⁴²

Thus

The structure factor S=σ/ρ and total correlation function h are related by

Using the invariant spherical harmonic expansion in Eq. (17) we obtain the friction function on the probe (tracer) particle

with the single time dependent relaxation function Cα11(x,t*)=ρSα11(x)e-t*(x2+3/2)/Sα11(x), α=0,1 ³⁶. The first term on the right hand side of (28) was first derived by Nägele et al⁴³ for spherically symmetric interacting particles and the second term in the present manuscript.

3. Microrheology under linear shear flow and no external magnetic field

We assume now that the viscoelastic response of the ferrofluid to variations of strain rate γ˙(t) is given by a linear constitutive equation for the shear stress ⁵

with shear relaxation modulus G(t). Because the shear stress divided by the shear rate has dimensions of viscosity, it is defined the complex viscosity η(t)=G(t)/γ˙(t). In frequency space, ω, σ(ω)=G(ω)iωγ(ω). Mason and Weitz ⁵ proposed that the bulk stress temporal relaxation scale of G(t) is the same time scale response as the microscopic stress relaxation of the complex viscosity η(t) that affects the particle motion. Thus for spherical particles undergoing translational diffusion they found G(ω)=iωη(ω)=ζ(ω)/(3πd) ⁵, whereas for rotational Brownian movement G(ω)=ζR(ω)/(πd3) ³⁴. This approximations has resulted in a very useful experimental technique called microrheology. Thus, here we make the same approximations of Mason and Weitz in order to obtain, in the case of translational movement, the linear relationships for the rheology of the ferrofluid and described by the complex frequency ω dependent shear modulus

where ı=(-1), ∆ζω : =Re∆ζω-ıIm∆ ζω, with Re, Im the real and imaginary parts of the friction function given by the average △ζ(t)=(1/3)(2△ζ⊥(t)+△ζ∥(t)). Using the Laplace transform

we find

Because for dipolar particles △ζ∥(t)=(4/3)△ζ⊥(t) ³⁶, then from Eq. (28) we determined △ζ(t) from these contributions as the real and imaginary parts of the complex friction function given by (31)

The complex viscosity of a viscoelastic fluid is defined as η=η'-iη″, with components η'=G″/ω, η″=G'/ω. Using (30) and (31) we find that

The normalized frequency ω*=ωt0 where the Brownian time t0=d2/D0 for a particle to diffuse its diameter. The short time particle diffusion coefficient D0=kBT/ζ0 and the friction ζ0=3πηsold with ηsol the solvent (sol) viscosity. These equations are valid for a concentrated ferrofluid under a stationary shear flow. They constitute the main results of this manuscript. Finally the storage modulus is

and the dissipative modulus

with h000=g000-1. For a hard ferrofluid suspension, the obtained expression of the complex modulus can be rewritten as

where we introduced the auxiliary definitions

The imaginary part is

Given the known material values of d,μ,ρ,ηsol,T for a ferrofluid, Eqs. (35-37) do not depend on adjustable parameters. The dimensionless quantities τs,α*=τs,α/t0=S,α11(x)/[8πx2+6] are short (s) relaxation times associated to the ferrofluid’s fluctuation in concentration in transversal (α=1) and longitudinal (α=0) modes to the wave vector. And τs,iso*=τs,iso/t0=1/[x21+1S,000x] arising from the isotropic (iso) interaction of particles. In Eqs. (35) and (37) the first and second terms arise respectively from the radially symmetric part of the potential (either this been LJ, Yukawa or hardcore type) ^36,43. The second term comes from the dipolar interaction derived in this manuscript. The above equations show that there are three temporal relaxation times of fluctuations in concentration. j₂(x) is the spherical Bessel function of order 2. μ*2=μ0βμ2/(4πd3), x = kd, k is the wave number. The reduced density ρ*=ρd3 and the projections of isotropic, longitudinal and transversal structure factor to the wave vector, are respectively

The Fourier-Bessel transform of the projections of the total correlation function h000:=g(r)-1, h△(r), and hD(r) were calculated according to the method of Ref. ⁴². These structural information results from the equilibrium position and orientation of the particles and they are determined from Brownian dynamic simulations performed with the Lammps package ⁴⁴, which allows the calculation of the averages

In Lammps simulations we used dimensionless units for, the particle energy at room temperature, length and mass which were obtained with respect to those of the ferrofluid Fe₂O₃, given by the values ²; ϵ0=5.453110-21J, d = 10^-8 m, and m ₀ = 2.70710^-21 Kg, respectively. They are related by the time scale t'=dm0/ϵ0 that measures the step size in the simulations.

Figure 1 is a log-log plot that depicts the behavior of τs,α,iso* versus reduced wavenumber x for a system with dipole moment μ*=1, and density ρ*=0.3. This figure shows that τs,iso*>τs,1*≥τs,0*, that is, isotropic fluctuations of concentration decay more slowly than longitudinal ones. Thus, transversal modes (τs,1-1) decay more slowly than longitudinal ones. This phenomenon has also been observed by other authors in ferrofluids and dipolar colloids ^46,47. In the long wave length limit kd≪1, it results τs,α=t0S,α11(x)/48π, which yields a relationship between the Brownian time t₀ and the microscopic time τs,α. A similar relationship between the Brownian and microscopic times, through colloid microstructure, was found by Bagchi et al⁴⁷. Figure 1 also implies that the short time τ,000 governing the collective isotropic density fluctuations has a larger contribution in both moduli than the longitudinal and transversal modes.