Stability of paramagnetic spheroid in precessing field

Stability analysis of paramagnetic prolate and oblate spheroidal particle in precessing magnetic field is studied. Bifurcation diagram is calculated analytically in dependence of magnetic field frequency and precession angle. The direction of particle in synchronous regime is calculated. Rotational dynamics and mean rotational frequency in asynchronous regime is found. Obtained theoretical model gives possibility to calculate analytically dynamics of the particle in limiting case (when motion is periodic). Theoretically obtained models were compared with experimental results of rod like particle dynamics in precessing magnetic field. Experimental results are in good agreement with proposed theory.


I. INTRODUCTION
Recently active topic of research is magnetic particle structure formation in different field configurations [1,2], which has been investigated experimentally, numerically and theoretically. Great interest in magnetic particles is due to possible applications. For example paramagnetic particles can be used to measure rheological properties of the fluid [3]. Using external magnetic field we can tune interaction between paramagnetic particles as can be seen in [4]. By tuning interaction between paramagnetic particles by changing properties of precessing magnetic field we can form chains [5] or if we use more complex particles like Janus particles we can obtain layers or tubes [6]. Magnetic Janus particles posses magnetic anisotropy and therefore can be modeled as oblate paramagnetic particles.
To fully understand these structure formation of magnetic particles, we need to start by fully understanding dynamics of individual particle. Previous studies have found experimentally the transition between synchronous regime and asynchronous regime for prolate paramagnetic particle in precessing field [7]. Numerical studies approved experimentally obtained results and expanded bifurcation diagram including oblate paramagnetic particle and wider range of field parameters [8]. Nevertheless full dynamics of the particle in asynchronous regime was not investigated.
In this work we formulate a model of the dynamics of the spheroidal paramagnetic particle which is applicable both for prolate and oblate particles, analytically obtain bifurcation diagram and calculate trajectory of the particle in stable "periodic" motion. Due to symmetry, all obtained equations can be used both for prolate and oblate particle. The only difference is the region of applicability, which is defined by bifurcation diagram. Although the equations for prolate and oblate particles are the same, due to different regions of applicability, the observable behavior differs.

II. THEORETICAL MODEL
A spheroidal paramagnetic particle has anisotropic magnetic susceptibility due to the anisotropy of its demagnetizing field factors N ,⊥ [9,10]. If we denote the symmetry axis of the particle by n, then magnetic moment m in an external magnetic field H = H h reads where V is volume of the particle and χ ,⊥ are magnetic susceptibilities, which depend on the demagnetizing field factors according to [11] χ ,⊥ = χ 0 /(1 + χ 0 N ,⊥ ), where χ 0 is the isotropic intrinsic susceptibility of the particle. The particle is in the fluid which rotates with angular velocity ω H = ω H e H . If we are interested in rotational dynamics of magnetic particle in precessing magnetic field which precesses with angular velocity ω H in stationary fluid, then this is equivalent to look at particle in stationary magnetic field and rotating fluid with angular velocity − ω H .
The balance of the viscous and magnetic torque (for low Reynolds number) reads: where ξ r is a rotational drag coefficient. ξ r for spheroidal particles can be found in [7].
If we introduce anisotropy frequency ω a = V H 2 (χ −χ ⊥ )/ξ r , then equation of the rotation of the particle can be written in the form: In general the angle between h and e H is arbitrary. For prolate particle ω a > 0, and for oblate ω a < 0.

A. Fixed points and their stability
Setting˙ n = 0 let us find fixed points. By simple calculations, which can be found in appendix A, cubic equation can be found, whose roots determine fixed points: where y = ( n · h) 2 , ω = ω H ω a and σ = h · e H . The last term of the cubic equation (4) is negative, that means that product of all roots is positive, therefore at least one root is positive. Since y = ( n · h) 2 , this gives that for positive y there exists fixed points. We can conclude that for every value of ω and σ there exists fixed point.
Equation (4) does not change if the sign of the ω (or ω a ) is changed. This means that fixed points for prolate (ω a > 0) and oblate (ω a < 0) particles are the same. The only difference is stability of these points, which depend on the sign of ω a .
Fixed point in situation, where magnetic field is kept constant and magnetic fluid rotates corresponds to synchronous rotation of the particle with the field frequency if fluid is stationary and magnetic field precesses (rotates).
In order to investigate stability of fixed point (synchronous regime) we introduce small perturbation ε from the synchronous state n. Since n is unit vector and its length can not change ε should be perpendicular to n and can be written in form: where ε 1 and ε 2 are small perturbations in two non-collinear directions.
After some arithmetical manipulations shown in appendix B we can find time dependence of ε: The stability of the regime is determined by eigenvalues of Jacobi matrix [12] and gives: In this situation ∆ and τ can be calculated from first and second derivative of the cubic equation (4): Therefore it is sufficient to look at the first and second derivative of cubic function f (y) at its roots f (y) = 0 to examine the stability of fixed points. The possible τ and ∆ values and consequences to the fixed point f (y) = 0 are following: • If τ < 0 (f ′′ (y) > 0) and ∆ > 0 (f ′ (y) > 0), then λ 1 ω a < 0 and λ 2 /ω a < 0 (7), therefore prolate particle (ω a > 0) is stable in this point, but oblate particle (ω a < 0) is unstable.
The coefficient at y 3 term in f (y) is positive, therefore f (y) is mainly increasing and if f (y) = 0 has only real one root then f ′ (y) > 0 at that root. The change in stability at that point is determined from f ′′ (y) which changes sign at y = 1 3 . Putting this in f (y) = 0 we get neutral stability curve: This equation has real solution only if σ < 1 √ 3 . Otherwise f (y) has one root with f ′′ (y) > 0 what means that prolate particle (ω a > 0) has stable fixed point and oblate particle (ω a < 0) has unstable.
If the equation f (y) = 0 has 3 real roots, then between these 3 real roots one must have f ′ (y) > 0 and f ′′ (y) < 0, one must have f ′ (y) > 0 and f ′′ (y) > 0 and one must have f ′ (y) < 0, which are between first two. therefore both particles prolate ω a > 0 and oblate ω a < 0 have one stable fixed point if f (y) = 0 has 3 real roots. The equation (4) has 3 real solution if: This equation has real solution if 0 < σ < 1 3 . If we put σ = 1 3 in equation (9) and (8)  The results of eq. (8) and (9) can be visualized by bifurcation diagram fig. 2, where ϑ is angle between h and e H ( h · e H = cos ϑ). In the bifurcation diagram region I II f (y) has one root, where in region I prolate particle has stable fixed point, but in region II oblate particle has stable fixed point, but in region III f (y) has three roots therefore prolate and oblate particle have stable fixed points.

B. Limit circle and rotational frequency
In the situation, where prolate particle does not have stable fixed points, there should be a limit circle. Luckily in this situation limit circle is "round" and in this situation n rotates in plane perpendicular to some stationary vector e L . As shown in appendix C, the equation for e L is the same as for n (3), but with opposite sign at ω a . This gives that limit circle for prolate particle is stable if oblate particle has stable fixed point. Corresponding fixed point of oblate particle is fixed point of the vector e L for prolate particle, which defines the plane of the rotation. The same reasoning is true for oblate particle without fixed points.
We also see that n of prolate particle is always perpendicular to e L , which is n of oblate particle. Therefore in situation where eq. (4) has three roots (region III in fig. 2) n of prolate particle is perpendicular to n of oblate particle (with equal |ω a |).
The rotation around limit circle is not uniform, but is periodic with mean angular velocity (D5):ω where y is the real root of (4).
It can be checked that mean angular velocity (10)  2, but is not zero at boundary between regions I and II fig. 2. For prolate particle if we slowly go from region III into region I than mean rotational frequency would continuously increase; same for oblate particle if it slowly goes from region III into region II. If Prolate particle slowly goes from region II into region I than there would be jump in mean rotational frequency at the boundary between regions; same for oblate particle going from region I into region II. The stable limit circle or fixed point near boundary between regions I and II is not very attractive, therefore it could take long time to reach stable state. But on the boundary of region III, the fixed point is already on the limit circle, therefore the stable state would be reached much faster.
The mean angular velocity dependence of the fluid rotation is shown in fig. 3 for prolate particle and for oblate particle it is shown in fig. 4. The place of the jump in mean angular velocity is found from (8) and the height of the jump can be found by setting y = 1 3 , which gives: If instead of rotating fluid, we will put particle in stationary fluid, but magnetic field rotates, we get that particle rotates asynchronously in this situation. And particle lags behind magnetic field with average angular velocity ω L . Therefore average rotational frequency of spheroidal magnetic particle in precessing magnetic field in asynchronous regime is

III. EXPERIMENT
To test provided theory, the experiment was made. In experiment superparamagnetic rods, which were synthesized according to the method detailed in [13], was immersed in stationary fluid and subjected to precessing magnetic field. Precessing field consists of constant field with strength H ω and rotating field with strength H r = 54 Oe, which is perpendicular to stationary field. Frequency of rotating field was fixed in range 0.05 − 2 Hz and the constant field was changed in step like manner. In this experiment mixture of glycerol and water in volume fractions 3:7 were used.
Rods were tracked with image processing algorithm in Matlab. An ellipse was used as the descriptive shape of particles. The angle and length of major axis were chosen as the descriptive parameters and plotted along time axis. In cases when the ellipse was very short, the image was either considered faulty or the rod had flipped along z axis.
Experimentally obtained angle of the rod as function of time, when stationary field is changed stepwise, is shown in Fig. 5. It can be seen, that for small constant field strengths H ω rod rotates slower than field frequency, because it is in asynchronous regime, but, when the stationary field strength H ω increases above critical value, rod starts to rotate with field frequency.
If we plot rotation frequency as function of stationary field strength H ω we obtain Fig. 6.
In Fig. 6 is shown only angular velocities when particle stabilized its rotation. It is assumed that particle is stabilized its rotation after change of the field strength if fit to the straight line of ϕ(t) slope has error less then 20%. Otherwise some first points are omitted in fit. If fit has less then 10 periods and still is not stabilized, this field parameter is omitted in Fig.  6.
The value of H ω where particle switch between synchronous and asynchronous regime we denote by critical field H c . Knowing H c , we can calculate that σ 2 = H 2 r H 2 r +H 2 c , when particle changes stability. Using neutral stability curve Eq. (8) or Eq. (9), which depends on critical σ value, we can calculate anisotropy frequency ω a of the particle. Anisotropy frequency depends on magnetic field strength, fluid viscosity and form of the particle. We will use that hydrodynamic anisotropy axis coincide with magnetic anisotropy axis, therefore ξ r (for rotation about the equatorial semiaxes) is given in [7] ξ = 8πηV Γ , where η is dynamic viscosity of the particle surrounding fluid, V is volume of the particle and Γ is geometric factor which depends on spheroid semi-axis a and b: where N and N ⊥ are demagnetization factors [9]. We define formfactor F of the particle Formfactor of the particle defines magnetic and hydrodynamic properties of the particle and depends only on the particle.
Obtained H c and calculated ω a gives rotational frequency of the particle in asynchronous regime which is in good agreement to theory (11) as can bee seen in Fig. 6.

IV. CONCLUSIONS
Analytic expressions of spheroidal paramagnetic particle stable synchronous and asynchronous rotation in precessing magnetic field is found. Results are found in non-inertial coordinate system rotating with magnetic field, which synchronous rotation change to stable point but asynchronous rotation to limit circle, which can be found analytically. It is shown that numerically obtained bifurcation diagram [8]  where ε 1 and ε 2 are small perturbations in two non-collinear directions.
Giving small perturbation ε to (A1) gives equation for perturbation in the linear form: Before we can go further, some relations should be computed, which comes from (A2) and its consequences (A3), (A4) un (A5): From eqs. (B2) and (5) The stability of the fixed point is determined by eigenvalues of Jacobi matrix J, where J is seen in (B3)  The eigenvalues of Jacobi matrix are: is trace of J divided by ω a and is determinant of J divided by ω 2 a , where ( e H · n) is excluded using (A4) and (A5).

Appendix C: Limit circle
In ω and σ region where no fixed points can be found, there should be limit circle. We assume that without fixed points, n will rotate around some vector e L with time varying frequency ω L . We assume that e L ⊥ n. So we can write ω L e L = n ×˙ n. Here we introduce one more unit vector e n with propriety˙ n = ˙ n e n . This gives relation between 3 unit vectors e L = n × e n .
We need to find how e L changes: Since vectors n, e L and e n are perpendicular we can express d 2 n dt is basis of these vectors as d 2 n dt = d 2 n dt · n n + d 2 n dt · e L e L + d 2 n dt · e n e n .
This simplifies (C1) to d e L dt = − 1 ˙ n d 2 n dt 2 · e L e n .
Next we expand from (A1): Putting this in (C2) gives Multiplying (A1) by e L and requiring that d n dt is perpendicular to e L gives that: by adding and subtracting (C4) multiplied by n to eq. (C3) we get: Next we introduce unit vector k which is in the direction of the projection of the vector h onto the plane, where n rotates and angle between k and n is φ (shown in fig. 7). So we can write that n · h = n · k k · h = k · h cos φ e n · h = e n · k k · h = k · h sin φ k · h 2 = 1 − e L · h 2 which simplifies eq. (10) ω L (t) = dφ(t) dt = ω a ( k · h) 2 sin φ(t) cos φ(t) + ω H ( e H · e L ) (D2) Separation of the variables and integration over the period gives: or we can write that mean angular velocity is: Since e L fulfils the same dynamic equation as n, then stable fixed points can be found from the same cubic eq. (4). Putting solutions of (4) into mean velocity and using relations y = ( e L · h) 2 and ( e L · h)( e L · e H ) = ( e H · h) 2 = σ 2 , we get: Using eq. (4) it can be rewritten in form: where y is the real root of (4).
Solving (D2) we can also write how φ(t) changes: where t 0 is some arbitrary constant and by changing k ∈ IR we can fulfil continuity of φ(t).