On the statistical hydrodynamics for a binary mixture of magnetic and non-magnetic particles

Dynamic properties of a binary mixture of magnetic and non-magnetic particles are considered using the method of non-equilibrium statistical operator. The generalized hydrodynamic equations are derived and analyzed. Based on this, the hydrodynamic collective modes spectrum is calculated and expressions for sound velocity and damping coefficients of these modes are obtained. We also propose the consistent scheme for calculation of hydrodynamic time correlation functions and thus the expressions derived for them in the paramagnetic state are discussed.


Introduction
Magnetic liquids, mixtures of magnetic and non-magnetic particles in the external fields of mechanical or electromagnetic origin have already occupied their significant place in chemical, electronic and other modern technologies.The recent experimental investigations [1] have already proved the existence of ferromagnetic phase in under-cooled liquid alloy Co 80 -Pd 20 with Heisenberg exchange interaction.A similar result has been obtained previously for metallic alloy Au-Co [2][3][4].This revives investigations in the area of physics of magnetic liquids.We note that this issue remains urgent regardless of whether there exists a ferromagnetic state in magnetic liquids or not, because the study of dynamic properties in the systems with translational and spin degrees of freedom in external magnetic fields is of great importance.In particular, one faces the above mentioned problem while describing spin relax-ation in liquids [5][6][7], and while investigating the dynamics of ferrocolloid systems [8][9][10][11][12], polar liquids [13][14][15][16][17][18] as well as surface absorption phenomena [18,19].That is why the investigations of thermodynamic, structural and dynamic properties of liquid magnets are very urgent for a deeper understanding and forecasting of their behaviour [4,8,9].
The investigation of time-dependent correlation functions as well as of generalized transport coefficients of a liquid mixture of magnetic and non-magnetic particles enable us to have a deep insight into the dynamic processes in the systems with coupled degrees of freedom.One of the most interesting tasks is to investigate the behaviour of hydrodynamic collective modes, which describe the properties of heat, sound, and mass fluctuations.From the experimental point of view, the study of these phenomena as well as the study of effects of sound propagation, spin waves and their dependence on temperature, concentration and pressure is very valuable.The explicit information concerning the bulk properties of mixtures of magnetic and non-magnetic particles is also essential for understanding of their interaction with active surfaces (metals, semiconductors or dielectrics), where the processes of metallic arrangement, absorption, desorption, and chemical transformation can take place.Another important aspect of this problem is the derivation of expressions for dynamic structure factors.It is known that these functions can be extracted from neutron scattering experiments.Such a theoretical study should be based on the statistical approach, and, particularly, on the equations of generalized hydrodynamics.A similar approach was previously applied to one-component magnetic fluid [20][21][22][23][24].The hydrodynamic collective modes of a Heisenberg ferrofluid were studied in [24].
Statistical hydrodynamics for a mixture of magnetic and non-magnetic particles in an external inhomogeneous magnetic field B(r; t) was considered recently in our paper [25].Note that from the physical point of view, the models of a binary mixture are more interesting to the theoretical studies (comparing with the case of one-component magnetic fluid [20][21][22][23][24]), because of a wider range of practical applications.For example, if a ratio of masses for particles in different species is large, such models can be used for the description of a ferrocolloid suspension.In particular, in paper [25] using the non-equilibrium statistical operator method, we derived the generalized hydrodynamic equations for a mixture which are valid for a description of both strong and weak non-equilibrium states.In the latter case, magnetic and non-magnetic subsystems were characterized by individual non-local thermodynamic parameters.As a result, the generalized thermodynamic relations and the generalized hydrodynamic equations were obtained.
The goal of this paper is to study the spectrum of hydrodynamic collective modes for a binary magnetic mixture, consisting of magnetic and non-magnetic particles.We use the perturbation like theory, which permits us to obtain two sound propagating modes for the system considered and three purely diffusive hydrodynamic modes.The expression derived for sound velocity generalizes the result obtained previously for a simple ferrofluid [24].
In chapter 2 we briefly overview the main ideas of the method of non-equilibrium statistical operator [26], used for derivation of the generalized hydrodynamic equa-tions for a binary mixture of magnetic and non-magnetic particles.The description of the model as well as the definitions of the basic dynamic variables is given in the next chapter.In chapter 4 we consider static correlation functions, and some useful relations for them are derived in the limit k → 0. This enables us to calculate the frequency matrix and the matrix of memory functions in the hydrodynamic limit.The collective mode spectrum is obtained and analyzed in detail in chapter 6.Then, in chapter 7, the problem of calculation of time correlation functions is considered.For this purpose, we develop a consistent scheme which permits us to find the weight coefficients describing the partial contributions from each mode to the hydrodynamic time correlation functions.The results obtained are discussed in the case of a binary magnetic mixture in paramagnetic state.

Theoretical framework: Non-equilibrium statistical operator method
Let us start with the Liouville equation: where classical part of i L is determined as a Poisson brackets of the function ρ with a classical part of the Hamiltonian of the system and the quantum one is determined as a commutator with its quantum part.Statistical operator ρ(x N ) is a function of phase variables x N = {r, p, ŝ} N , where N is a total number of particles.Following Zubarev's method of non-equilibrium statistical operator [26], we can rewrite equation (1) in the form: where ε → 0 after the thermodynamic limit, ρ q is the so-called quasi-equilibrium statistical operator.Nonzero right-hand side of the equation ( 2) imposes the boundary conditions which destroy the time reversal symmetry of the Liouville equation and permits one to choose the retarded solutions.
In hydrodynamic region we can restrict our consideration to the set of most slow dynamical quantities { Pα } which are thought to determine a weak non-equilibrium state.According to the ideas of abbreviated description of the system, these variables are naturally connected with the set of conserved quantities.The Gibbs-like form for statistical operator ρ q can be found from the extremum principle for information entropy with fixed parameters of abbreviated description Pα t , preserving the normalizing condition Sp ρ q (t) = 1.This gives where the parameters {F α (t)} are defined from the self-consistency conditions: or Sp Pα ρ(x N , t) = Sp Pα ρ q (x N , t) .
The index α = {i, k} in (3) and ( 4) denotes a combination of discrete index i, labelling the variables, and wave vector k, so that this summation means: Taking into account the properties of projection operators, the formal solution of equation ( 2) can be written [26] as follows where Îα (t) = (1 − P(t)) Ṗ α are the generalized fluxes, Ṗ α ≡ i L Pα , and is the operator of time evolution with the Kawasaki-Gunton projection operator P q (t).Operator P q (t) is simply connected with the generalized Mori projection operator P(t) We note that P(t) acts on the dynamical variables as follows and has the following properties: Using the solution (5) for the non-equilibrium statistical operator ρ(t) one can derive the generalized transport equations in the form [26] where are the generalized memory functions, or generalized transport kernels.The equations (4), (8), and (9) make up the basic set of nonlinear equations which describes both strong and weak non-equilibrium states.
Let us consider now in more detail a weak non-equilibrium case.Then, the deviations δ Pα (t) = Pα t − Pα 0 (10) of the averages Pα t = Sp Pα ρ(t) from the equilibrium values Pα 0 = Sp Pα ρ 0 (x N ), where is the equilibrium statistical operator, as well as the deviations of the intensive quantities δF n (t) = F n (t) − F 0 n can be considered as small quantities.The relations between them follow from the self-consistency conditions (4).This gives where with ∆ Pi ≡ Pi − Pi and ( Â, B) denotes a correlation function defined as follows ( Â, B) = Note that ( Â, B) ≡ Â B 0 in the case of classical treatment.For small deviations from an equilibrium we can rewrite (8) in the form where the projection operator is given by Using ( 13) and the equality we get the generalized transport equations in the form where iΩ = ( ˙P , ∆ P + )(∆ P , are the frequency matrix and the matrix of memory functions, and δ P (ω) denotes the Fourier transform of δ P (t).The matrix equation for the Laplace transforms (∆ P , ∆ P + ) z of time correlation functions (∆ P , ∆ P + ) t , (∆ P , ∆ P + ) z = ∆ P , 1 has the structure similar to (15), namely, We note that the retarded correlation Green functions Hence, the poles of the Laplace transforms of G

Model Hamiltonian and hydrodynamic variables
Let us consider the system consisting of N 1 non-magnetic and N 2 magnetic particles posed in an external magnetic field.Taking into account the interaction between subsystems, the Hamiltonian of such a system can be written as follows [29,30] Here and further the subscripts 1 , 2 or superscripts in parentheses (1) , (2) indicate the non-magnetic and magnetic subsystems, respectively.Thus, H 1 and Ĥ2 denote the Hamiltonians of isolated non-magnetic and magnetic subsystems, H int describes interaction between them, and Ĥex is the energy of the spin interaction with an external magnetic field.The Hamiltonian H 1 of non-magnetic subsystem can be taken in the classical form where p (1) j is the momentum of jth non-magnetic particle; V (11) (r jl ) denotes the potential of interaction between two non-magnetic particles j and l; and m 1 is a mass of nonmagnetic particle.The term Ĥ2 in (20), consists of classical "liquid" part H 2L , having the same structure as the term H 1 , and quantum part, which describes the spin subsystem and can be taken in Heisenberglike form.Other terms in ( 20) could be written as follows where V (12) (r ij ) is the potential of interaction between i-th non-magnetic and j-th magnetic particles, and B(r i , t) describes an external magnetic field.The Liouville operator, which corresponds to Hamiltonian (20) can be written in the form where V (11) , and i Ls is a purely quantum part of Liouville operator, determined by the commutator It is well-known that in the hydrodynamic region the dynamics near an equilibrium is mainly determined by the slowest processes, which are directly associated with the densities of the conserved quantities.Hence, for the model considered five parameters of abbreviated description { Pi }, i = 1 . . . 5 can be chosen.They are: the partial densities of particles' numbers n1 (k) and n2 (k), the density of total momentum p(k), the spin (or magnetization) density m(k), and the density of total energy ε(k), where the index α indicates now the spatial α-component of the corresponding vector, and ε(1 (1) V (11) (r ij ) + 1 2 (2) For the set of dynamical variables Pi , the microscopic equations of motion could be written in the following form (see Appendix) where Pi (k) and P α l (k) denote the scalar [n 1 (k), n2 (k), ε(k)] and vector [ p(k), m(k)] variables, respectively.The terms R i (k) and R α l (k) appear in (30) mainly due to an effect of inhomogeneous external magnetic field.If the applied field B(r, t) does not depend on r [i.e., B(r, t) = B(t)], these terms disappear and the corresponding dynamic variables become conservative.

Static correlation functions
Let us assume that the external magnetic field is static, homogeneous, and directed along the '0z' axis, so that B(r, t) = {0, 0, b}.In this case only the z-th component of spin density mz (k) ≡ m(k) is conservative (see (A.10)).This means that the basic set of dynamic variables includes four scalar dynamic variables and one vector-like dynamic variable p(k).Moreover, because the transverse and longitudinal dynamics can be considered separately, for the study of the longitudinal dynamics we deal in fact with five scalar variables P (k) = {n 1 (k), n2 (k), p(k), m(k), ε(k)}, where p(k) denotes the longitudinal component of p(k).We also note that in this sense the transverse dynamics of the "liquid" subsystem (for the model ( 20)) is equivalent to the transverse dynamics of a simple liquid.
In order to calculate the frequency matrix ( 16) and the matrix of memory functions (17) we have to analyze in more detail the properties of static correlation functions constructed on the dynamic variables P (k).
Let us define now a static correlation function (â, b) as an average, constructed on the deviations â and b from its equilibrium values: In order to have the relation with thermodynamics we consider the equilibrium statistical operator ρ 0 (x N ) as a Gibbs distribution for the grand canonical (µ 1 , µ 2 , V, T, b)-ensemble, so that where Ω(µ, b, T ) is the thermodynamic potential; β is the inverse temperature; n1 , . . ., ε are the variables (31), taken at and b denotes an external magnetic field.
For an arbitrary operator â and parameter γ (â does not depend on γ) it is easy to prove the equality where the averaging is performed with the operator (33).If we put, for example, â = n1 and γ = µ 1 in (35), one gets Here and further the quantity, written by a capital letter, denotes the average value of a corresponding operator, written in a small letter, for instance, N 1 = n1 0 .In the same manner using (35) we obtain Considering (32) as the definition of scalar product, it is seen in ( 36)-(38), that the set of Pi (k = 0) is not orthogonal in the sense that the non-diagonal elements of matrix ( P , P + ) are nonzero.For practical needs it is often more convenient to use orthogonalized dynamic variables.Such a new set can be found using Gram-Schmidt orthogonalization procedure.Let us introduce new dynamic variables being mutually orthogonal with one exception: the first two variables in (32), namely, n1 and n2 , form the non-orthogonal (2×2) sub-block of matrix ( P , P + ) with the non-diagonal elements (n i , nj ), i = j.In order to proceed further, we define the projection operator where (n, n) −1 ij denotes the {ij}th element of the inverse matrix (n, n+ ) −1 , and consider the new variable defined by It is obvious that ŝ is simply connected with the spin density m and orthogonal to n1 and n2 in the sense that (ŝ, n1 ) = (ŝ, n2 ) = 0.
Note that the density of momentum p(k) is orthogonal to all the variables by the definition, and where mN = N 1 m 1 + N 2 m 2 is a total mass of all particles in the mixture considered.
It is important to stress that the projection using the operators (39) and (41) permits us to perform the transition from one ensemble to another.For example, the projection, given by (39), describes, in fact, the transition from (µ, V, T, b) to (N, V, T, b) ensemble.Thus, the magnetic susceptibility in (N, V, T )-ensemble can be defined on the 'projected' variable ŝ, namely, Using (35) we can prove another equality for the entropy in grand canonical ensemble so that an expression for the specific heat in (µ, V, T, b)-ensemble can be written as follows C µ,b = T (∂S/∂T ) = β 2 (ω, ω), and performing the transition to (N, M, V, T )-ensemble, as it was done in (44), we obtain C N,M = β 2 ((1 − P n − P ŝ) ω , ω) = β 2 ( ĥ, ĥ).
Taking into account the relations (40) and (42), it is seen that the set of variables {n 1 , n2 , p, ŝ, ĥ} is orthogonal in the sense mentioned above.Generalizing the results obtained, one can introduce the new dynamic k-dependent variables which are mutually orthogonal.One exception is for the variables n1 (k) and n2 (k).
For the dynamic variables ŝ(k) and ĥ(k) one has with the projection operators defined as follows The correlation functions, constructed on the variables (46), can be considered as the generalization of the thermodynamic derivatives, derived above (see ( 36)-( 38), (44), and (45)), for nonzero values of k.Hence, ( ĥ(k), ĥ(−k)) = 1 where S ij (k) are the so-called partial structure factors, and χ T,N (k) and C N,M (k) denote the generalized magnetic susceptibility and specific heat for the mixture of magnetic and non-magnetic particles, respectively.

Frequency matrix and matrix of memory functions
Taking into account the symmetry of the correlation functions under time inversion and spatial symmetry operations, or performing the direct calculations, one can prove that Ẏ i (k), Ŷj (−k) = 0, if ( Ŷi = p and Ŷj = p) or ( Ŷi = p and Ŷj = p). (54) Moreover, because of the equality , it can be shown that to calculate the frequency matrix we have to find only the correlation functions which involve the operator of momentum density: In the limit k → 0 one can prove that (see [24]) where U(r N ) is the total potential energy of the system Let us remind now the expression for pressure P , which follows from the equilibrium treatment (see [24]) Comparing ( 55) and ( 56), one obtains From ( 57) one can conclude that in homogeneous external field ∂ B(r, t)/∂r = 0 the pressure of a system with the 'isotropic' interactions, which depend only on a distance between particles, is expressed via the average of the stress tensor J αβ p (k).In particular, using (57) we get It is seen in ( 58) that only the longitudinal component p(k) of the current density p(k) is coupled with the density and thermal fluctuations.Taking into account this result and using (57), the following relations for static correlation functions can be derived where k = |k|.For the orthogonalized variables we get Using Gibbs-Duhem equation the right-hand sides of (59) can be rewritten in another form where are the coefficients of magnetostriction, isothermal compressibility and isobaric thermal expansion for the mixture of magnetic and non-magnetic particles, respectively.The derivative (∂V /∂P ) in the second equation of ( 60) is defined at constant magnetization M.After some algebra using Gibbs-Duhem and Maxwell relations, one has , Using the results obtained above and performing the generalization for k = 0 we can write some elements of the frequency matrix iΩ as follows Let us consider now the elements iΩ p,n i of the frequency matrix iΩ, which involve the densities of particles' numbers nν The elements of matrix (n, n+ ) in the limit k→0 are the derivatives β −1 ∂N i /∂µ j (see ( 36)), therefore it can be proved that in the canonical (N 1 , N 2 , V, T, b)-ensemble one has The right-hand side of this equation can be rewritten by performing the transition to the ensemble with constant pressure where v ν = (∂V /∂N ν ) T P b is the partial molar volume per particle in the νth species.Taking into account that and using Gibbs-Duhem equation, we find for isobaric processes This allows us to express the matrix elements (63) via the partial molar volume v ν .
In general case k = 0 one has [27] iΩ p,n where the generalized k-dependent functions v ν (k) are introduced.The matrix of time-dependent longitudinal memory functions for the set of dynamic variables (46) is defined as follows where the k-dependence of the hydrodynamic projection operator P(k), ŝ + . . ., ĥ ĥ, ĥ −1 ĥ, is omitted for the sake of simplicity.For the generalized fluxes in the limit k → 0 one has where f i are regular longitudinal parts of fluxes.Using (68) and (70), it is seen that all the matrix elements of φ(k, t) are proportional at least to k 2 .We also note that in the hydrodynamic limit the Markovian approximation for the memory functions is in fact explicit and, therefore, can be used for the subsequent calculations.Summarizing the results obtained above one gets the frequency matrix in the form It can be proved that due to the symmetry properties, the matrix of memory functions in the hydrodynamic limit k → 0, ω → 0 has an opposite structure to iΩ, namely, All the other elements in (72) could be neglected because they are proportional to the higher powers of k.Thus, in the hydrodynamic limit, the Laplace transforms of memory functions can be written as follows where L ij are the kinetic coefficients, defined by Green-Kubo-like formulas

Hydrodynamic collective modes
The spectrum of collective modes can be found from the equation ( 19) and gives in fact the set of eigenvalues for the matrix Considering δ ≡ ik as a small parameter and taking into account that iΩ ∼ δ and φ ∼ δ 2 , we can develop the perturbation-like theory which permits us to solve the equation (75 being the sound velocity.Three other eigenvalues z 3,4,5 describe the diffusive modes and vanish (z 3,4,5 = 0) in this approximation.Using the results, obtained in the previous section, the expression (78) can be rewritten as follows where the mass density ρ = mN/V is introduced.To find the next approximation let us consider z in (75) as a series over δ, so that Putting (80) into (75), the new set of algebraic equations for the coefficients z 0 and D can be derived.In such a way we find for the propagating sound modes while for purely diffusive hydrodynamic modes one has where D * and D i are the corresponding damping coefficients.For sound modes this results in the expression and damping coefficients D i for the purely diffusive modes can be found as the roots of an algebraic equation of third order with coefficients p i , defined via the elements of matrix κ ≡ ω + φ with an additional constraint κ pp = 0. Namely, it can be shown that

Hydrodynamic time correlation functions
The matrix equation (18) for the Laplace transforms of time correlation functions can be rewritten as follows where Ã(k, z) = z• 1 − T (k), the matrix T is defined by (76), and F (k, z) = ( Ŷ , Ŷ + ) z and F0 (k) are the matrices of Laplace transforms of time correlation functions and static correlation functions, respectively.The formal solution of (85) is given by with Ã−1 (k, z) = (z 1 − T ) −1 being the matrix inverse to Ã(k, z).The ijth element ( Ã−1 ) ij of this matrix can be written as an algebraic adjunct (Ad [A ij ]) of element A ij , divided by the determinant ∆( Ã), so that one has For Ad [A ij ] we have where the dependence on z is omitted for the sake of simplicity.The determinant of matrix T can be presented as a product of its eigenvalues {z γ }, ∆( T ) ≡ 5 γ=1 z γ .
Therefore, one gets The eigenvalues {z γ } were found in the previous section up to the second order with respect to the small parameter δ (see ( 81) and ( 82)).Note also that, because the matrices ω and φ have the opposite structure (see ( 71) and ( 72)), we can write Hence, using (88), one gets where are the so-called 'weight coefficients', describing the contribution from corresponding collective modes.Because of {z k } are found as series in δ (see (81) and (82)), the expression for G γ ij can be written as follows It should be stressed that the expression (94) is valid only in the hydrodynamic limit.A more general result can be found from (86) using simple matrix relations.Let us introduce the notation X for the matrix of eigenvectors of T [see (76)], so that T • X = X• Z, with Z = ||z γ δ γγ ′ || being the diagonal matrix of eigenvalues {z γ }.It is easy to show that the matrix Ã(z) = z• 1 − T has the same set of eigenvectors X with the eigenvalues {z − z γ }, therefore we get Taking into account that the matrix (z 1 − Z) is diagonal, the last expression can be rewritten in the following form Hence, comparing (92) and (95), we conclude that This expression can be used to calculate the weight coefficient for an arbitrary value of wavenumber k.In particular, the result (94) can be easily reproduced from (96), when k is small.Note also that for practical applications, the expressions (94) can be more convenient if analytical results are needed.

Discussion and concluding remarks
In this paper the spectrum of collective hydrodynamic modes for a binary ferromagnetic mixture is studied based on the rigorous statistical treatment.We have found two complex-conjugated collective modes, which are responsible for propagating the sound excitations.All the other modes are purely diffusive and describe the processes connected with relaxations of temperature-, mass-and spin-fluctuations.The results, obtained in the hydrodynamic region, are valid for an arbitrary binary mixture under the same external conditions, namely, when the external magnetic field is homogeneous and the interactions are isotropic.In particular, it was shown that the sound velocity is inversely proportional to the adiabatic compressibility in the ensemble with constant magnetization.This conclusion coincides with the result found previously for a simple magnetic liquid [24].We also derived the expression for the damping coefficient of sound excitations [see (83)] and the equation ( 84), which permits one to calculate the damping coefficients of purely diffusive collective modes.
For zero magnetic field (b = 0), if we consider paramagnetic phase, the spin and liquid subsystems are decoupled and the damping coefficient of spin mode D m can be easily found from (84).This results in D m = ϕ ss .The expressions for two other damping coefficients, describing the hydrodynamic thermal and concentration excitations, coincide formally with the results known for a binary mixture of simple fluids.The same can be said about the expression for hydrodynamic time correlation functions.For instance, the Laplace transform of "spin density-spin density" time correlation function can be written in the form used widely in the theory of solid magnets.In a more general case, when b = 0, an additional coupling between spin and translational degrees of freedom appears and the expressions for time correlation In a similar way, the corresponding expressions for the spin-density current can be derived.We find (1) [V (11) ] ′ (r ij ) (r ij p J α ε,2 has the same as J α ε,1 with the substitution '1'→'2' and implying difference between ε(1) i and ε(2) i ; and l( =i =j) J(r ij )J(r il )(ŝ i , [ŝ j , ŝl ])r α ij P k (r ij ) e ikr i , (A.14) i e ikr i (2) j [V (12) ] ′ (r ij ) ij (t), which give the spectrum of collective excitations, are determined by the matrix equation det |z − iΩ + φ ε (z)| = 0. (19)