Compressibility of ferrofluids: Towards a better understanding of structural properties

This paper addresses a computational method aimed at obtaining the isothermal compressibility of ferrofluids by means of molecular dynamics (MD) simulations. We model ferrofluids as a system of dipolar soft spheres and carry out MD simulations in the NPT ensemble. The obtained isothermal compressibility computed via volume fluctuations provides us with a strong evidence that dipolar interactions lead to a higher compressibility of dipolar soft sphere systems: the stronger the dipolar interactions, the bigger is the deviation of the compressibility from the one of a system with no dipoles. Furthermore, we use the isothermal compressibility to calculate the structure factor of ferrofluids at low values of wave vectors, i.e. in the range where it is difficult to predict its behaviour because of a problem with accounting for long-range particle correlations that give the main contribution to the structure factor in this range. Our approach based on the interpolation of the structure factor and the computed isothermal compressibility allows us to obtain the smooth structure factor in the range of low wave vectors and the reliable fractal dimension of the clusters formed in the system.


Introduction
Nowadays many applications of soft matter heavily rely on the structural behaviour of their compounds.In turn, the structural behaviour of these materials is usually very sensitive to both external stimuli and intrinsic characteristics of their compounds.One of the classical examples of dipolar soft matter materials with a large variety of phase and structural transitions is a ferrofluid, also addressed sometimes as a magnetic fluid [1].Ferrofluids represent a class of fluids consisting of single-domain magnetic particles suspended in a nonmagnetic liquid carrier.Magnetic and steric interaction between particles in ferrofluids might result in their self-assembly even in the absence of an external magnetic field.
The two past decade can be fairly called as a boom in studying the microstructure of ferrofluids and the major contributions aimed at understanding microstructural properties were published .The experimental technique usually used for this purpose is the small angle neutron scattering (SANS) [1,2,5,10,13,14].The scattering pattern obtained in these experiments transforms into a Contribution to the Topical Issue "Advances in Computational Methods for Soft Matter Systems" edited by Lorenzo Rovigatti, Flavio Romano, John Russo.
a e-mail: elena.minina@univie.ac.at structure factor which carries information about the sample microstructure.In order to be consistent with experimental measurements, many following computational and theoretical works put an effort to obtain structure factors of ferrofluids as this property can be compared to real experiments [15][16][17][18][19][20]22].The structure factor S(q) with the wave vector q is usually analytically calculated by evaluating a Fourier transform of the radial distribution function g(r).Nowadays, there are two known theoretical approaches to calculate g(r) for ferrofluids.The first one is based on the diagram expansion [19], which gives accurate predictions for the pair distribution function of ferrofluids at low density and dipolar interactions not strong enough to lead to self-assembly.The second approach takes into account the self-assembly occurring at moderate dipolar interactions by considering the system as an ideal gas of chain clusters of different lengths [17].Despite both theoretical approaches are proven to be reliable for predicting the structure factor first peak, they fail at low wave vectors [17][18][19][20]22].Here, the inaccurate behaviour of the structure factor is related to the fact that the method of diagram expansion neglects high-order correlations which are important at long distances.In turn, the chain model lacks interchain interactions which are important for concentrated systems.In computer simulations due to the finite size of the system, the structure factor obtained for low values of q is also not reliable [15-18, 20, 22].
The knowledge about the microstructure of ferrofluids is not only valuable by itself, but also is directly linked to thermodynamic properties such as pressure, compressibility factors and heat capacity which alter together with the microstructure [20,[24][25][26][27][28][29][30].To our knowledge, nobody has studied isothermal compressibility of ferrofluids by means of computer simulations despite it has a strong implication for thermodynamics and dynamics of ferrofluids [31].The aim of this paper is to fill this gap and to explain how dipolar interactions influence the isothermal compressibility.Using molecular dynamics computer simulations, we compute the isothermal compressibility χ T of dipolar hard sphere systems in an isothermal-isobaric (NP T ) ensemble as a measure of volume fluctuations.This method captures well the qualitative behaviour of χ T which clearly displays the influence of dipolar correlations.Moreover, isothermal compressibility has a direct relevance to the structural properties of ferrofluids, since it measures the structure factor at the zero limit of the wave vector length as follows [32]: Here ρ is the number density of dipolar particles, k B and T denote the Boltzmann constant and temperature of the system, respectively.The notation q signifies the magnitude of the wave vector q (q = |q|).The range q → 0 is exactly the range of wave vectors where conventional methods experience difficulties in predicting the structure factor.Note that the structure factor in this range holds an important information regarding long-range particle correlations and its scaling contains information about the fractal dimension of a sample that characterizes its microstructures [15].In order to obtain the structure factor in this important range of wave vectors, we combine two simulation approaches.We compute the structure factor at large wave vectors in a conventional way by means of molecular dynamics simulations of dipolar hard sphere systems in the canonical NV T ensemble.At low q values, we have only one point of the structure factor obtained from NP T simulations according to eq. ( 1), where the pressure P is the equilibrium one of the NV T ensemble.
In the next step, we join the results of two simulations techniques to interpolate the structure factor in the range of small values of wave vectors.This way, we estimate the structure factor of ferrofluids at low wave vectors and verify the range of validity of the aforementioned analytical approaches.
The paper is organized as follows.Section 2 describes the details of modelling the systems of dipolar spheres and explains both methods of computer simulations used in this paper.The simulation results and their analysis are given in sect.3 and summarized in sect.4.

Simulation method
We consider a system of N spherical particles of diameter σ.Each particle possesses a central magnetic moment μ whose rotation is related to the particle rotation.Particles interact with each other via two potentials: the Weeks-Chandler-Andersen (WCA) potential [33] to account for the particle steric repulsion (2) and the potential of dipole-dipole interactions Here r ij is the distance between particles i and j and r ij is the displacement vector connecting the centers of these particles.Vectors μ i and μ j denote their dipolar moments of magnitude μ.The constant μ 0 stands for the vacuum permeability and the parameter is the strength of the WCA potential.The strength of the dipolar interactions can be described by the parameter λ = βμ 0 μ 2 /4πσ 3 .The inversed product of temperature T and the Boltzmann constant k B is denoted as All simulations mentioned in this paper were carried out in ESPResSo [34,35] at the following system parameters: number density ρ is in the range between 0.01 and 0.1 with step 0.01, and the dipolar interaction strength varied from λ = 1.0 up to λ = 4.0 with step 1.0.
The first molecular dynamics simulations of dipolar soft spheres were done in the canonical ensemble (NV T simulations).We placed 512 dipolar particles into a cubic simulation box of length l = 3 N/ρ.The initial positions of the particles and orientations of their magnetic moments were randomly chosen.We use the periodic boundary conditions to model the system in a bulk [36].The dimensionless temperature T = 1 of the system was kept constant by means of the Langevin thermostat [37] at the dimensionless friction coefficient γ = 1.We also employed the dipolar-P3M algorithm to compute the longrange dipolar interactions between particles in metallic periodic boundary conditions [38].The system was evolved using the Velocity-Verlet algorithm with the time step Δt = 0.01.The equilibrium of the system was reached during 2 × 10 4 time steps and the statistics was accumulated during subsequent 6 × 10 6 time steps.The target observables in these simulations were the radial distribution function (RDF) g(r) and the pressure P .The RDF is related to the structure factor of the system as follows [32]: Note, the isothermal compressibility χ T can be also calculated via RDF as expressed below: The pressure P obtained in NV T simulations was used to set up the next molecular dynamics simulations in the isothermal-isobaric ensemble (NP T simulations).This time, we kept the pressure constant and used the volume of the simulation box for this purpose.We changed the volume of the simulation box according to the acceptance rule for making a volume move [36], which has the following form: Here, U (V o ) and U (V n ) denote the energies corresponding to the systems with volume V o before the volume move and V n after the volume move was made.The attempt to change the volume from U (V o ) to U (V n ) was made every 10 4 time steps by making a trial volume move The volume was changed equally along all three directions.The remaining routine of the simulation protocol stayed the same as for NV T simulations, including the Langevin thermostat to preserve the temperature, the dipolar-P3M method to compute dipolar interactions and the Velocity-Verlet scheme to propagate the system.We increased the size of the system to N = 2048 particles and the time spent to collect statistics to 10 7 .We are interested in volume fluctuations which are related to the isothermal compressibility as expressed below: Note that we also considered smaller systems consisting of N = 512 particles.The comparison of data obtained by means of NP T simulations of smaller systems to those of bigger systems did not show any dependence of the isothermal compressibility of the system size.
In the next section, we compare the isothermal compressibility obtained using eqs.( 5) and (7).

Results and discussion
Figure 1 shows how the isothermal compressibility χ T of the dipolar soft sphere systems alters with number density ρ.The compressibilities obtained by two different methods (using eq. ( 5) and eq.( 7)) demonstrate the same qualitative behaviour: the compressibility decays as the system becomes denser and therefore less compressible.Moreover, the compressibilities calculated via two different methods even coincide at λ = 1 and 2. For higher λ, we observe the discrepancy between the methods which grows with increasing particle concentration: the compressibility calculated via RDF is consistently lower than its counterpart at a given density because of the finite size effect of the system in NV T simulations which has been enhanced by integration in eq. ( 5).In the following, we will only talk about the compressibility calculated in NP T ensemble.
To estimate the impact of dipolar correlations on the compressibility, we compare it to the compressibility of pure hard sphere systems.The latter can be analytically calculated by differentiating volume V with respect to pressure P at constant temperature T [32]: Using the Carnahan-Starling approximation for the compressibility factor Z of hard spheres [28] where ϕ = ρv 0 is the particle volume fraction and v 0 = πσ 3 /6 is the particle volume, the isothermal compressibility for hard spheres writes as follows: The isothermal compressibility obtained using eq.( 10) is depicted as a black dashed curve in fig. 1. Comparing data for dipolar and nondipolar systems, we notice that the compressibility of the system of dipolar spheres resembles that of a nondipolar one.Moreover, they almost coincide at λ = 1 (see fig. 1  is due to the effective interparticle attraction caused by dipolar interactions and self-assembly of dipolar particles.Both these phenomena lead to a higher compression without the need of exerting any additional pressure on the system.It turns out that the compressibility can be well predicted the power law f (ρ) = Aρ B .The result of the fitting of the compressibility by the power law is depicted in fig. 1 as a solid red curve and the obtained coefficients A and B are given in table 1.
Since the isothermal compressibility carries the information about the structure factor at the zero limit of the wave vectors q → 0, the obtained above results can help us to estimate the structure factor at low wave vectors q and therefore add to understanding the structural properties of ferrofluids.
As has been mentioned in the introduction, the structure factor obtained by means of computer simulations is very accurate for values of q corresponding to short interparticle distances in the real space.Therefore, we plot the structure factor obtained via NV T simulations using eq.( 4) only for qσ > 1.5 as solid curves (see fig. 2).Different colours correspond to different number densities ρ.In the range qσ > 3, the structure factor qualitatively coincides with the one obtained in previous works [15][16][17][18][19]:  the increase of the density leads to a growth of the first peak as the number of particle pairs in the sample becomes greater (compare the curves within each subfigure); the position of the first peak slightly shifts towards larger wave vectors with the enhancement of dipolar interactions as it results in the effective attraction diminishing the interparticle distance (track the shift by looking at the subplots from (a) to (d) consequently).Circles at q = 0 in fig.2(a)-(d) represent the SF obtained via isothermal compressibility using eq.( 7).In order to get the missing part of the structure factor at 0 < qσ < 1.5, we interpolate the points of the structure factor at q = 0 and 1.5 ≤ qσ ≤ 2.5 by a quadratic polynomial The interpolation coefficients can be found in table 2. The part of structure factor obtained by the interpolation is depicted as dashed curves of the colour corresponding to the figure colour scheme.
It is worth noticing that the structure factor at zero q decays with increasing particle concentration at λ = 1.At λ = 2, these points almost converge to one value, whereas for high λ = 3 and 4, the structure factor at zero wave vectors increases with the density of the system.This crossover is happening because dipolar particles tend to form chains at λ > 2 and these chains are captured by the structure factor also at distances L σ, i.e. when q → 0. The scaling of the structure factor in this range by the function S(q) ∼ (qσ) −D allows us to characterize the formed structures because the meaning of the coefficient D is the fractal dimension.To find the scaling, we use the obtained by interpolation values of the structure factor at 0.8 ≤ qσ ≤ 2.0.This is the range where the structure factor has a rapid linear decay in log-log scale.The results of the scaling are provided in fig. 3 as black lines and the obtained values for the coefficient D are given in table 3.At the highest λ = 4 considered here, the coefficient D ∼ 1.This confirms the chain formation and additionally justifies the accuracy of compressibility computed in this work.These results are consistent with scalings obtained in previous works on structure factor calculations [15][16][17].However, this coefficient has not reached unity for systems at λ = 3 at any particle density.We attribute this to the fact that although the chains can already be found in these systems, they are mostly short, i.e. dimers-like, with few longer chains.Therefore, they give low impact in long-range chainlike spatial correlations that results in such a small value for D. Now having the reliable structure factor even at low values of q, we have the opportunity to compare it to the two known theoretical predictions through the complete range of the wave vectors.Figure 3 shows structure factors for systems with different λ and ρ obtained in several ways.Blue curves in fig. 3 correspond to the analytical method, employing the diagram expansion for the calcu-lation of g(r) [18,19].This method is advised to be used for very dilute systems with weak dipolar interactions.In case of higher values of dipolar interactions, self-assembly starts playing an important role.To this respect, another method to calculate g(r) was proposed which considers the system as an ideal gas of chains formed by magnetic particles [17].The equilibrium distribution of chains by their lengths is obtained at given values of λ and ρ by minimizing the free energy density functional [6].However, this method neglects the interchain interactions that make the model not applicable for concentrated systems.Black curves in fig. 4 represent the structure factor obtained via this model.Figure 4 compares the black and blue curves to the structure factors obtained using our approach based on simulation data and their interpolation (red circles).In the following, we call the latter as simulation data for convenience.At low lambda (λ ≤ 2), all three methods are in good agreement in the range of qσ > 3 (see fig. 4(a) and (b)).However, the simulation data are closer to the blue curve at ρ = 0.1 and slightly deviate from the prediction made by the chain model.However, the chain model gives better predictions of the structure factor from λ = 3 and higher, since this approach accurately accounts for the distances between particles inside chain structures that is reflected in the structure factor first peak.This is an expected result obtained in the former works [17][18][19].However, we would like to draw your attention to the range of low values of q.One can note that the structure factor of dilute systems (ρ < 0.05) at λ ≤ 2 obtained by the diagrammatic method coincides with the simulation data for low values of q as there are almost no long-distance correlations in these systems.The chain model becomes less accurate for systems with higher ρ (see fig. 4 at ρ = 0.1 at small values of q), albeit it seems to work for low concentrations and λ = 3.Finally, at λ > 3, we cannot rely on any of the two models when analysing the structure factor for low values of q.

Conclusion
This paper is aimed at computing the isothermal compressibility of dipolar soft sphere systems using computer simulations.We employed molecular dynamic simulations which we carried out for two ensembles.Firstly, simulations in the NV T ensemble, which have provided us with the pressure of the system to start the next simulations in the NP T ensemble.We computed the isothermal compressibility using data obtained in both simulations: via integration of the RDF obtained in NV T simulations and via volume fluctuations of the simulation box in NP T simulations.The isothermal compressibility obtained by these methods show the same qualitative behaviour, though the first method is less reliable because of the finite size effects.The isothermal compressibility decreases with an increasing particle concentration.However, dipolar interactions slow down the observed decay.The difference between the obtained isothermal compressibilities of dipolar soft sphere systems and those of nondipolar ones increases with We increase the number density ρ in each column going from left to right so that ρ has the following values: 0.01, 0.04 and 0.1 accordingly.Circles represent data obtained using two simulation approaches and interpolation.Black and blue curves correspond to theoretical predictions of the structure factors using the chain model and the model based on the virial expansion, respectively.
the strength of dipolar interactions.Apparently, the enhanced interparticle attraction due to dipolar interactions, which also results in chain formation at high enough λ, facilitates the compressibility of such systems.The obtained result estimating an impact of dipolar correlations to compressibility is also relevant to other kinds of dipolar systems, e.g., magnetic gels and elastomers, for which compression is an essential basis for any possible application.
The high importance of the isothermal compressibility is also determined by its relevance to the structural property.Namely, the compressibility is the structure factor at the zero limit of the wave vector q.This is exactly the range where both theories as well as computer simulations usually experience difficulties in predicting the structure factor.In this paper, we offer the method to predict the structure factor at low wave vectors which is based on the interpolation of the obtained compressibility at q = 0 and the structure factor in the reliable range of q to the range of small wave vectors.The advantage of this method is that the structure factor obtained by interpolation has a smooth behaviour in the range of small q and therefore allows us to find an accurate scaling of the fractal dimension characterizing structures formed in the sample which is in agreement with previous results.

Fig. 1 .
Fig. 1.Isothermal compressibility χT as a function of the number density ρ obtained in NP T simulations using eq.(7) (depicted as red circles) and NV T simulations using eq.(5) (depicted as blue diamonds).The black dashed curve is the isothermal compressibility for the pure hard sphere system.Red curves are the result of the fitting of NP T simulation data by the power law Aρ B , where the coefficients A and B are given intable1.Each subfigure corresponds to different λ -strength of dipolar interactions: (a) λ = 1.0;(b) λ = 2.0; (c) λ = 3.0;(d) λ = 4.0.

Table 2 .
Interpolation coefficients a, b and c obtained by interpolating the structure factor of a dipolar hard sphere system by a polynomial of the second order F (q) = ax + bx + cx 2 at different strengths of the dipolar interaction λ.

Fig. 4 .
Fig. 4. Structure factor S(q) obtained for dipolar hard sphere systems at different densities ρ depicted in different colours.Each subfigure is plotted for a different strength of dipolar interactions: (a) λ = 1.0;(b) λ = 2.0; (c) λ = 3.0;(d) λ = 4.0.We increase the number density ρ in each column going from left to right so that ρ has the following values: 0.01, 0.04 and 0.1 accordingly.Circles represent data obtained using two simulation approaches and interpolation.Black and blue curves correspond to theoretical predictions of the structure factors using the chain model and the model based on the virial expansion, respectively.

Table 1 .
Coefficients A and B obtained by fitting the isothermal compressibility by the function f (ρ) = Aρ B .
(a)) as dipolar interactions are very weak and give almost no contribution to the isothermal compressibility for such dilute systems.With increasing values of λ (look at fig.1(a) to (d) consequently), the curves start deviating towards higher compressibility.This