Pressure and compressibility factor of bidisperse magnetic fluids

In this work, we investigate the pressure and compressibility factors of bidisperse magnetic fluids with relatively weak dipolar interactions and different granulometric compositions. In order to study these properties, we employ the method of diagram expansion, taking into account two possible scenarios: (1) dipolar particles repel each other as hard spheres; (2) the polymer shell on the surface of the particles is modelled through a soft-sphere approximation. The theoretical predictions of the pressure and compressibility factors of bidisperse ferrofluids at different granulometric compositions are supported by data obtained by means of molecular dynamics computer simulations, which we also carried out for these systems. Both theory and simulations reveal that the pressure and compressibility factors decrease with growing dipolar correlations in the system, namely with an increasing fraction of large particles. We also demonstrate that even if dipolar interactions are too weak for any self-assembly to take place, the interparticle correlations lead to a qualitative change in the behaviour of the compressibility factors when compared to that of non-dipolar spheres, making the dependence monotonic.


Introduction
Magnetic fluids or ferrofluids are suspensions of nano-sized magnetic dipolar particles in a non-magnetic carrier. These fluids show a strong response to an external magnetic field, however, remain in a liquid state. The remarkable combination of these properties distinguishes magnetic fluids from other magnetic materials and is successfully used in a variety of applications ranging from mechanical engineering (e.g. magnetorheological fluid dampers) to medical applications (e.g. cancer treatment by hypothermia) [1][2][3][4][5][6][7][8][9].
Conventional ferrofluids used in applications are usually polydisperse with a wide distribution of particle size. Polydispersity inevitably influences properties of ferrofluids. The first attempts to demonstrate this effect were made for bidisperse ferrofluids-the simplest approximation of polydisperse ferrofluids-by splitting the particles into two fractions of approximated characteristic sizes. The early works on bidisperse dipolar systems reveal the striking change in microscopic properties due to bidispersity in the absence of an external magnetic field. For example, bidispersity causes the poisoning effect that does not allow particles to form long chains [10][11][12][13][14] in contrast to monodisperse ferrofluids where long chain formation was confirmed theoretically as well as by computer simulations [15][16][17]. Later, Novak and co-authors managed to compute the partial radial distribution functions and structure factors of bidisperse dipolar spheres for different granulometric compositions, both theoretically and in comp uter simulations [18].
Along with the interest in the understanding of microstructure, the thermodynamic macroscopic properties of ferrofluids and the influence of polydispersity on them have also been studied partially [19][20][21][22][23]. In particular, one of the main issues crucial for applications of ferrofluids is to understand how In this work, we investigate the pressure and compressibility factors of bidisperse magnetic fluids with relatively weak dipolar interactions and different granulometric compositions. In order to study these properties, we employ the method of diagram expansion, taking into account two possible scenarios: (1) dipolar particles repel each other as hard spheres; (2) the polymer shell on the surface of the particles is modelled through a soft-sphere approximation. The theoretical predictions of the pressure and compressibility factors of bidisperse ferrofluids at different granulometric compositions are supported by data obtained by means of molecular dynamics computer simulations, which we also carried out for these systems. Both theory and simulations reveal that the pressure and compressibility factors decrease with growing dipolar correlations in the system, namely with an increasing fraction of large particles. We also demonstrate that even if dipolar interactions are too weak for any self-assembly to take place, the interparticle correlations lead to a qualitative change in the behaviour of the compressibility factors when compared to that of non-dipolar spheres, making the dependence monotonic.
Keywords: magnetic fluid, pressure, compressibility factor, computer simulations (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. much magnetic fluids are compressible without the presence of any external factors changing their properties. To answer this question, monodisperse dipolar systems were examined. Experimentally, it has been shown that dipolar interactions lower the pressure induced by the particles even at relatively weak dipolar strength [24]. Later, this was confirmed theoretically as well as in computer simulations [25,26]. However, little is known about the influence of polydispersity on such properties as pressure and the compressibility factor. In one of the studies we are aware of, these properties were computed for systems of bidisperse dipolar hard spheres using Monte Carlo simulations [27]. The systems considered in this work are nearly monodisperse with weak dipolar interactions and contain only a small portion of particles of larger size and with stronger dipolar interactions (up to 20% of the total volume fraction of particles). Therefore, the full spectrum of granulometric compositions has not been tested yet, and the behaviour of the pressure and compressibility factors has not been entirely understood. We would like to underline the importance of this lacking information because the compressibility factors, as it turned out, reveal an interesting dependence on the granulometric compositions.
The present work is aimed at filling this gap by clarifying how the granulometric composition alters pressure and compressibility factors of bidisperse dipolar systems with weak dipolar interactions. We compute the pressure and the compressibility factors of such systems using molecular dynamics simulations. We vary the fraction of small and large particles so that we consider various granulometric compositions including monodisperse systems. While changing the granulometry, we keep the total volume fraction of particles constant to preserve of the total volume of magnetic material for all considered granulometric compositions. Therefore, despite different granulometry, these systems have the same satur ation magnetization. In addition, we calculate pressure and compressibility factors theoretically extending the previously developed theory for calculating the radial distribution functions [18] which is based on the method of diagram expansions. Considering two type of systems (systems of dipolar hard spheres and systems of dipolar soft spheres), we determine the range where the theoretical predictions can be made based on the conventional model of dipolar hard spheres, where the softness of the particles plays a crucial role and where this theoretical approach should not be trusted at all. This paper is organized as follows: the theoretical approach for calculating pressure and the compressibility factor of systems of bidisperse dipolar spheres is presented in section 2. The details of molecular dynamics simulations for modelling such systems and computing the desired properties are given in section 3. The results are presented and discussed in details in section 4 and summarized in section 5.

Theoretical approach
We consider a bidisperse magnetic fluid comprised of N dipolar particles which are split into N s small particles of diameter d s and volume v s = πd 3 s /6 and N l large particles of diameter d l and volume v l = πd 3 l /6 (d s < d l ). Each particle possesses a central magnetic moment represented as a vector (see figure 1) of a corresponding magnitude: µ s for small particles and µ l for large particles. We consider a 3D system, in which the rotation of dipoles is associated with particle rotation and translational motion with the displacement of the particle's center of mass. The positions of the particles and directions of their magnetic moments are determined by the steric repulsions and dipole-dipole interactions. The latter can be written as: where β = 1/k B T is the inverted product of the Boltzmann constant k B and temperature T, µ 0 is the vacuum permeability and r ij is the vector connecting the centers of the particles i and j. The strength of this potential is characterized by the parameter which is the ratio of the dipolar energy of two particles at side-by-side close contact with parallel dipole moments, to the thermal energy. Here d ij = (d i + d j )/2 is the distance between the particle centers at close contact. In our systems, three parameters of dipolar interaction have to be distinguished: the interaction strength between two small particles, that of two large particles, and the one between small and large particles. Another important detail, concerning modelling a magnetic fluid, is the particle steric interactions. Traditionally, there are two type of steric potentials that are used to describe steric repulsions between particles of ferrofluids: the hard sphere potential [28] βU and the soft sphere or Weeks-Chandler-Andersen (WCA) potential [29] βU 2d ij is the cut-off radius and ε is the interaction strength. As it has previously been demonstrated, the choice of the steric potential plays a big role in theoretical modelling. On one hand, the hard sphere potential is much more convenient to use for analytical calculations and one can get a direct result without involving numerical methods [30]. More crucial is that some properties like the structure factor basically do not depend on the steric potential for at least moderate concentrations [18]. On the other hand, the softness of the steric potential can influence the precision of theoretical calculations leading to inaccurate results, because, at some conditions, the real steric interactions are underestimated by the hard sphere repulsion. In these cases, the WCA potential has to be used to model the soft steric interactions correctly [18,25,26]. In the present work, we considered both potentials for calculating the pressure P and the compressibility factor Z of bidisperse ferrofluids in order to understand what potential gives an accurate theoretical prediction of these quantities.
The pressure of bidisperse systems can be written as a combination of partial pressures: where n = N s /(N s + N l ) and (1 − n) = N l /(N s + N l ) are the fractions of small and large dipolar particles comprising the sample and P ss , P ll are the partial pressures triggered by small particles only and large particles only, and P sl is the cross-term due to the presence of both types of particles. According to the equation of state [28], the pressure P and the compressibility factor Z are linked by the following relation: where ρ is the number density of particles. Therefore, in order to calculate both the total pressure and the compressibility factor, we have to find the partial pressures. The latter can be obtained through the corresponding partial radial distribution functions g αγ (i, j) (where α, γ ∈ {s, l} throughout the paper) as where the product of the partial RDF and derivatives of the potentials are averaged over all orientations of magnetic moments Ω i and Ω j . Without loss of generality, this averaging for any index i can be written as Here the angles ω i and ξ i denote polar and azimuthal angles of the magnetic moment µ i , correspondingly. The partial RDFs for systems of bidisperse dipolar hard-and soft-sphere particles have recently been calculated for different granulometric compositions using the method of diagram expansion [18]. According to this method, a RDF can be presented as a virial expansion over the particle volume fraction ϕ up to the 2nd order, i.e. truncated with B 3 (r ij )-the third virial-like coefficient: The coefficients B 2 (r ij ) and B 3 (r ij ) are the contributions of two-particle and three-particle correlations and have the following form: with the Mayer function which has been expanded over λ ij up to the 4th order. More detailed information concerning the RDF calculation including the diagrams taken into account can be found in [18]. Then using properties of the diagrammatic expansion, we obtain the final expressions for both pressure and compressibility factors: where Z * denotes the total compressibility factor of binary hard sphere mixtures. The next terms are the partial pressures but without steric contribution β −1 Z * αγ ρ αγ : . This steric contribution has already been taken into account in Z * . For more details of how the expression (11) is derived, see appendix A. Dividing the total pressure (equation (11)) by the particle number density ρ gives us the total compressibility factor of the system.
Note that Z * has previously been obtained by Mansoori and co-authors [31], who solved the Percus-Yevick integral equation to find the RDF of such system, and has the following form: Here ϕ = ϕ s + ϕ l is the total volume fraction of the system (ϕ s = ρ s v s ; ϕ l = ρ l v l ). Equation (12) allows us to calculate steric contributions more accurately by taking into account not only two-and three-particle correlations as it is assumed in the method of diagram expansion, but also the correlations of the higher order. It is important to realize that dipolar particles of ferrofluids, in reality, are soft due to the surface stabilization. This softness is usually modelled by the WCA potential. For our purpose that particular potential has one clear disadvantage, i.e. it does not allow us to perform analytical integrations for obtaining direct analytical results. The common way to avoid this problem is to use hard spheres with effective diameters d * . This effective diameter mimics the softness of particles in order to switch from hard sphere particles to the soft ones and is approximated as [28]: In section 4, we discuss whether the simple model of hard spheres with effective diameter gives us an accurate approximation for the pressure and compressibility factors or the softness of the particles has to be explicitly taken into account.

Computer simulations
To study the desired macroscopic properties of magnetic fluids computationally, we performed molecular dynamics simulations in the simulation package ESPResSo [32,33]. To model our system, we randomly throw N = N s + N l particle into a cubic simulation box, which size was scaled to the required volume fraction of particles ϕ. The diameter of the small particles is set to unity in reduced units σ s = 1. Large particles of diameter σ l are bigger than small particles by the factor of σ l /σ s = d l /d s . Each particle possesses a magn etic moment of the following magnitudes in dimensionless units: µ s = λ ss σ 3 s for small particles and µ l = λ ll σ 3 l for large particles. The set of parameters considered in this work is chosen to be consistent with real systems and is shown in table 1. The simulation box with the particles, whose number varied from 514 to 5791 for different total volume fractions ϕ, was replicated using metallic periodic boundary conditions [34].
The steric repulsion between the particles was modeled by the WCA potential with the cut-off radius r c = 2 −5/6 (σ i + σ j ), where σ i and σ j are the diameters of the ith and jth particles. Dipolar interactions of the particles were computed using the dipolar-P3M method [35]. The system has been propagated using a velocity Verlet integrator with the fixed time step ∆t = 0.01. The temperature of the system was kept constant at T = 1 by means of a Langevin thermostat with dimensionless friction 1. The system was first equilibrated for 2.5 × 10 5 time steps. After equilibration, we calculated the pressure using the virial equation [34] collecting the statistics in the next 7.5 × 10 6 time steps.
In order to understand the influence of granulometric composition on the pressure, we varied the volume fraction ϕ s and ϕ l of small and large particles respectively, keeping the total volume fraction ϕ = ϕ s + ϕ l constant. All the granulometric compositions of small and large particles considered in this work are presented in table 2. Note that among them there are compositions with zero volume fraction of one of the particle types, i.e. the monodisperse cases, where all the particles are of one size (See GC1 and GC11 in table 2). Figure 2 shows the pressure of bidisperse dipolar spheres for systems of type S1 obtained in molecular dynamics simulations. Three different types of symbols correspond to the simulation data obtained for three different total volume fractions of particles. Each point corresponds to one granulometric composition. For better understanding, it is convenient to look at figure 2(a) where the pressure P is plotted as a function of the volume fraction of small particles ϕ s . Going from left to right, we start from a purely monodisperse system consisting of large particles only (CG1: ϕ s = 0, ϕ l = ϕ). Then we remove some portion of large particles and add a portion of small particles according to table 2, keeping the total volume fraction of particles constant, so that at the end, we are left with a pure monodisperse system comprising of small particles only (CG11: ϕ s = ϕ, ϕ l = 0). This order, where the leftmost and rightmost points correspond to monodisperse  cases with respective granulometric compositions GC1 and GC11 and with growing volume fraction of small particles ϕ s between these points, is valid for all figures in this section. Figure 2(b) shows the same data for pressure as figure 2(a), but as a function of the total number density of the particles ρ = 6/π ϕ s /σ 3 s + (ϕ − ϕ s ) /σ 3 l . Both figures 2(a) and (b) demonstrate that an increase of the volume fraction of particles leads to a pressure growth. However, the pressure is also sensitive to the granulometric composition: the higher the volume fraction of small particles at ϕ = const , the larger the total number density and, therefore, the pressure is higher as it has previously been confirmed for monodisperse systems [25]. The pressure of the systems of the type S2 displays the same trend, though the data are not shown here. Figure 3 compares the simulation data (symbols) with our theoretical predictions for the pressure P of dipolar soft spheres represented by green solid curves and dipolar hard spheres with effective diameters depicted as red curves, both given by equation (11) in section 2. Plots in the upper row correspond to the systems of type S1, whereas plots in the  Pressure of bidisperse dipolar spheres (plots in the upper row are for systems of the type S1 with λ ll = 1.81, plots in the lower row are for systems of the type S2 with λ ll = 2.34) as a function of the total number density ρ. Each column corresponds to different total volume fractions of particles: ϕ = 0.05, 0.09, 0.13. The curves depict theoretical predictions: the black dashed curves are reference curves for pure hard sphere system; the red curves are for systems of dipolar hard spheres with effective diameters; the green curves are for dipolar soft sphere systems. Symbols represent simulation data for all granulometric compositions from table 2 starting from GC1 (leftmost circles) to CG11 (rightmost circles). lower row are obtained for systems of type S2. In order to examine the significance of the contribution of dipolar interactions, we also plot the pressure of non-dipolar hard spheres obtained from the product of the number density and the compressibility factor given by equation (12) without effective diameters. All theoretical predictions show a linear growth in the pressure with the particle number density ρ, however, the dipolar interactions significantly lower the pressure. Whereas the rightmost points on the curves for dipolar and non-dipolar systems, which correspond to the systems comprising mainly of small particles, almost coincide, the leftmost points already differ noticeably, since in this case, the main dipolar contribution comes from large particles, whose dipolar interactions are stronger than those of their smaller counterparts. They enhance the dipolar correlations that lead to an effective attraction of particles and, therefore, lower the pressure.

Results and discussions
The pressures of systems of the type S1 with weaker dipolar interactions converge to the simulation data noticeably better than the ones of the systems S2 especially for the systems with low number density. This tendency is observed because the smaller the number density is, the more large particles the system contains. Such particles provide a higher impact of the dipolar interactions. Furthermore, it is already known that the method of diagram expansion becomes less accurate when λ > 2 as it has recently been reported in the work on pressure calculation of monodisperse dipolar spheres [25,26]. The discrepancy between the theory and simulation data increases towards the smaller number density especially for the systems of dipolar hard spheres, where not only the dipolar contributions to the pressure are underestimated, but in addition, steric repulsion is not modelled accurately. Apparently, it is not sufficient to model the softness of the steric repulsion via a hard sphere potential and by assigning each particle an effective diameter, especially not for computing pressures of dense systems. However, our theory based on the method of diagram expansion remains reliable in its prediction of the pressure of bidisperse dipolar systems at all possible granulometric compositions under the condition that the dipolar interactions of both particle fractions are fairly weak λ ss(ll,sl) < 2. In this case, the theory is relatively accurate even for systems with a rather high particle volume fraction (see figure 3(c) for S1 at ϕ = 0.13).
Despite the fact that dipolar correlations do not alter the qualitative behaviour of the pressure for bidisperse fluids, compressibility factors demonstrate the unexpectedly different behaviour. Figure 4 depicts the compressibility factors Z of bidisperse dipolar systems as a function of the total number density ρ. The structure of figure 4 is the same as that of figure 3 with red curves representing theoretical predictions for bidisperse hard sphere systems with effective diameters and green curves depicting theoretical predictions for bidisperse soft sphere systems. The symbols are the compressibility factors obtained in simulations. We also plot the compressibility factors of pure hard sphere systems (black dashed lines). There are two main facts to be mentioned about the compressibility factor. Firstly, the method of diagram expansions gives less appropriate theoretical predictions of the compressibility factors. Multiplying pressure by the total number The compressibility factor of bidisperse dipolar spheres (plots in the upper row are for systems of the type S1 with λ ll = 1.81, plots in the lower row are for systems of the type S2 with λ ll = 2.34) as a function of the total number density ρ. Each column corresponds to different total volume fractions of particles: ϕ = 0.05, 0.09, 0.13. Symbols represent simulation data and the curves depict theoretical predictions: the black dashed curves are reference curves for pure hard sphere systems; the red curves are for systems of dipolar hard spheres with effective diameters; the green curves are for dipolar soft sphere systems. density, which is of the order of 10 −2 − 10 −1 , inevitably magnifies the discrepancy between the theory and simulation data. The discrepancy diminishes for less concentrated systems, i.e. consisting of large particles mostly with the higher impact of dipolar correlations. Plotting Z also uncover the deviations of the soft-sphere model from the simulation data, which are basically invisible in figure 3. The more important fact, however, is the difference in qualitative behaviour of the compressibility factors of dipolar and non-dipolar systems. As one can see in figure 4, the compressibility factor of non-dipolar hard spheres does not monotonically depend on the number density of particles and has a minimum. However, this is no longer valid for dipolar systems. As it is clearly seen in figure 4, the compressibility factor of bidisperse dipolar spheres grows monotonously with the total number density of particles. Apparently, dipolar correlations play a significant role leading to the deviation from the compressibility factor of a pure hard sphere system. Comparing the plots in each column in figure 4, one can observe that the difference between the compressibility factors of dipolar and non-dipolar systems grows increasing large particle fraction due to the enhancement of dipolar interactions in these systems.

Conclusion
In this work, we studied the behavior of the pressure and compressibility factors of bidisperse dipolar systems in the absence of an external magnetic field using both molecular dynamics simulations and a theoretical approach. In the latter case, we employed the method of diagram expansion to calculate these properties for bidisperse hard spheres with an effective diameter to model their softness and bidisperse soft spheres. We also considered different granulometric compositions in order to observe how it affects the pressure and the compressibility factor at fixed total volume fraction of particles. We found that the pressure decays gradually at a subsequent replacement of small particles by large ones at fixed total volume fractions. By adding large particles the impact of dipolar interactions grows prompting an effective attraction between particles and as a result a decrease in the pressure. The theoretical predictions made for both dipolar hard and soft spheres reproduce this tendency. They do, however, deviate from the simulation data when large particles with the strength of the dipolar interactions greater than 2 are dominating in the systems.
Despite such a good convergence of the theoretical predictions for pressure to its simulation data, the reliable agreement for compressibility has been found only for dilute systems ϕ = 0.05. One of the most important observations is that the dipolar correlations change the qualitative behaviour of the compressibility factors. Namely, dipolar correlations not only lower the compressibility factor but lead to a monotonous behaviour when adding large particles with stronger λ, whereas for non-dipolar systems the compressibility factor depends on the granulometric composition in a non-monotonous fashion. .
(A.1) Taking the derivative of the magnetic dipole-dipole interaction potential U dd (i,j) with respect to the distance r ij between the particles i and j, we find that this derivative can be expressed through the potential of dipolar interactions itself as follows: In order to proceed further with the calculations, we have to make a brief revision of how the diagrams are constructed according to the method of diagram expansion. In particular, we draw attention to the dashed lines in figure A1 or figure 1 in [18]. Here m dashed lines represent the expression [−βU dd (i, j)] m /m!. Since the RDF g(i, j) is a sum of the diagrams, the multiplication of the RDF by βU dd (i, j) leads to an appearance of one additional dashed line between the particles i and j on each diagram with the correction factor −m!/(m − 1)! = −m. This way, we renew the set of diagrams. However, not every renewed diagram gives a non-zero contrib ution to the RDF. As it was pointed out in the earlier work on the RDF calculations for monodisperse magnetic fluids [30], if a diagram has at least one vertex with an odd number of dashed line entries, it gives a zero contribution to the RDF. Taking this fact into account, we can consider only the diagrams depicted in figure A1 for calculating the pressure. Note, the multipliers in front of the diagrams are the aforementioned correction factors. Calculating contributions of the diagrams to RDFs requires averaging over the orientations of magnetic moments, therefore, equation (A.1) can be rewritten as: with the following integrand of the second integral obtained by substituting the partial RDFs g αγ (r) (see equations (18) and (19) in [18]) and without the Jacobian r 2 : • for small or large particles only (α = γ ) (A.5) • for small-large particles (α = s and γ = l) βU dd (i, j)g αγ (r) = −2I 2b (r, d α , d γ , λ αγ ) − 4I 2c (r, d α , d γ , λ αγ ) (A.6) The exact formulas for computing contributions from each diagram can be found in appendix B. Now we split equation (A.4) into two terms. The first term contains the steric contribution only which is expressed though the partial compressibility factors Z * αγ . The second term includes dipolar correlations. Therefore, g(r) αγ in the first integral of equation (A.4) is the partial RDF, excluding those diagrams containing only steric contributions, i.e. g αγ (r) = g αγ (r) − I 2a − I 3a (see equations (7) and (10) of [18] for two-and three-particle diagrams I 2a and I 3a ). This leads to the following expression for the total pressure: where Z * denotes the total compressibility factor of binary hard sphere mixtures, and the next two terms in equation (11) are the partial pressures P * αγ = P αγ − β −1 Z * αγ ρ αγ (α, γ = {s, l}) excluding the steric contribution β −1 Z * αγ ρ αγ which has already been taken into account in Z * .