An explicit expression for the static structure factor for a multi-Yukawa fluid: the one-component case

ABSTRACT By using the mean spherical approximation, we obtain an analytical expression for the static structure factor (SSF) for a monodisperse system of particles interacting through a potential given by a hard-sphere contribution and M Yukawa terms. This expression depends on scaling matrix Γ, which is determined by solving a set of nonlinear equations. Our theoretical results show that using three Yukawa terms in the closure relation greatly improves the accuracy when compared with hypernetted-chain closure and Monte Carlo simulation data, which display a secondary low-k peak in the SSF, due to the formation of an intermediate range order structure governed by a short-range attraction and a long-range repulsion. We discuss the appearance of such a peak in terms of the microstructure order given by the radial distribution function. Following the original proposal made by Waisman (Mol Phys. 1973; 25:45–48), we give an explicit expression that improves the structural properties of a hard-sphere system.


Introduction
The classical statistical mechanics of fluids requires an effective pair potential to adequately represent a real system and, consequently, to be able to predict the structure and thermodynamic behavior of the system. [1] A model effective potential of special interest and useful is represented by hard spheres plus Yukawa tails. This model has been applied to study a wide range of fluids as simple liquids, electrolytes, molten salts, liquid metals, dense plasma, colloidal dispersions, micelles, and microemulsions and for predicting the liquid-glass ideal transition. [2][3][4][5][6][7] Since the Ornstein-Zernike (OZ) equation can be solved for a Yukawa model in the mean spherical approximation (MSA), [8][9][10] it is possible obtain an analytic expression for the static structure factor (SSF). For multi-Yukawa fluid mixtures with factorizable coefficients, the SSF has been calculated by Ginoza et al. [11], and by Blum et al. [12] using Laplace transform. In addition, using the Baxter factorization, Cruz-Vera et al. calculated the SSF for two Yukawa terms with an arbitrary number of components, [13] while for the case of one-component two-Yukawa fluid, Liu et al. obtained the SSF applied Hoye's approach. [14,15] The SSF is an important quantity in the study of physics of atomic liquids, charged colloids, and globular proteins in solution. In particular, the theoretical structure factor of a two-Yukawa fluid, obtained by Liu et al. [15], has been used for fitting small-angle neutron-scattering (SANS) intensity distribution; it also has been applied to study lysozyme solutions. [16,17] Moreover, the analytical SSF has been used as input in mode coupling approach to determine structurally arrested states, [18] as well as in the formalism of the generalized Langevin equation. [6] The main result of this article is the determination of an analytical expression for the SSF for a one-component fluid using a model consisting of a hard-sphere contribution plus M Yukawa terms. This is done using the Baxter factorization for multi-Yukawa closure solution given by Blum et al. as a function of the Ginoza scaling parameter Γ. [10,19] In addition, we use an iterative Newton-Raphson method with finite difference to determine the physically meaningful root of the scaling diagonal matrix. The remaining of this article is organized as follows: in Section 2, the Blum solution for a one-component fluid with a multi-Yukawa closure is reviewed [10,20]; in Section 3, the explicit calculation of the SSF for this closure is shown. Some results for the analytical SSF are shown in Section 4 for a system modeled by a short-range attractive and long-range repulsive potential, represented by two or more Yukawa terms; regarding these results, we discussed the formation of intermediate range order (IRO) peak that occurs due to the presence of a short-range attractive and long-range repulsive potential. Section 5 gives our conclusions, and in the Appendices we summarize the relevant expressions that allow us to improve the properties of a hard-sphere system, and the used method to find the scaling matrix involved in the calculation of the reported results.

Theory
In the theory of simple liquids, the direct correlation function, c(r), is defined in terms of the total correlation function, h(|r|), through the OZ equation [1]: where h r ð Þ ¼ gðrÞ À 1, g(r) is the radial distribution function (RDF), which is proportional to the probability density of finding a particle at a distance r from a particle situated at the origin of coordinates, and ρ is the number density of the hard-core particles. The solution of the OZ equation for a closure of M Yukawa terms with factorizable coefficients, together with the exact excluded-volume condition, h r ð Þ ¼ À1 for r < σ, was first found by Blum et al. in terms of the Ginoza scaling parameter, [10,19] using the Weiner-Hopf-Baxter factorization method. [21] In Equation (2), K (n) and z n are the closure parameters calculated for the system under consideration by means of the MSA in such a way that c r ð Þ ¼ Àβu r ð Þ for r > σ, where u(r) is the pair potential, σ is the effective diameter of the hard-core particles, β À1 ¼ k B T, k B is Boltzmann's constant, and T is the absolute temperature. Now, in the Generalization of the MSA (GMSA), K (n) and z n are determined by imposing a thermodynamic consistency condition through the calculation of some known state property. [8,22] The factors δ (n) are quantities arising in the solution for the case of mixtures [10] and, in this work, they will be taken only as auxiliary mathematical parameters and, in general, it will be assumed that δ n ð Þ ¼ 1 for all elements n = 1, 2, . . ., M. Equation (1), written in Fourier space and using Baxter's factorization, [10,21] is: where S(k) is the SSF, andĥ k ð Þ andĉ k ð Þ are the Fourier transforms of h(r) and c(r), respectively. The Baxter factorQ k ð Þ is given bỹ andQ Àk ð Þ is the complex conjugate ofQ k ð Þ. The solution forQ k ð Þ or Q(r) has been discussed in [10,20] and has been found as where Q 0 r ð Þ ¼ 0 for r > σ, and C n ð Þ e Àz n r À e Àz n σ ð Þ ; for r < σ; also D n ð Þ ¼ Àδ n ð Þ a n ð Þ e z n σ : The coefficients appearing in the above equations may be written as: ! a n ð Þ e z n σ ; The quantities labeled with the superscript 'HS' are the contributions of the hard-sphere repulsion and have the form with where η ¼ πρσ 3 =6 is the reduced density or volume fraction.
and X n ð Þ ¼ δ n ð Þ þB n ð Þ z n 1 À e Àz n σ ð ÞþσΔ n ð Þ : The elements of the matrix M appearing in the Equation (11) are given by Now, in the case of the closure of M Yukawa terms with factorizable coefficients, the symmetry conditions for the direct correlation function imply that Q(r) must be symmetric in the origin and, therefore, that Å n ð Þ must be proportional to X (n) . The most convenient scaling relation which may be used to obtain this proportionality comes from Equation (11) [10] where Γ mn are the elements of the scaling matrix Γ, which has dimensions M Â M. Hence, from Equations (15) and (18), X (n) can now be written as [10,20] where N À1 nm are elements of the inverse matrix N , and this matrix has elements with and Using the expression for Q(r) in the Yukawa closure (see Equation (2)), it may be seen that the closure condition is equivalent to the following equation [20]: The expression for the scaling matrix Γ may be obtained from Equation (26) when it complemented by the set of M(M − 1) equations obtained from the symmetry conditions in the origin for the direct correlation function and also from the symmetry conditions for the contact pair correlation function. These symmetry conditions, in matrix notation, are given by [10] where I is the identity matrix and z is a diagonal matrix whose diagonal elements are given by z n . Within the Yukawa closure, the termB n ð Þ denotes the excess internal energy [10] and, consequently, the contact pair correlation function may be written as [9,10] and the inverse of the compressibility, via the fluctuation theorem, as [9] β ρ κ À1 ¼ @βp @ρ where p is the pressure and g HS ð Þ σ þ ð Þ is the contact pair correlation function in the PY approximation, given by [23] 3. Static structure factor in the Yukawa closure: general case for a one-component fluid In this section, we obtain an analytical expression for the SSF of a one-component fluid in the general case of a Yukawa closure of M terms. For convenience, from here onwards, we will take dimensionless quantity for wave number k;kσ and for the parameter z n ;z n σ, and we will write the number density as ρ ¼ 6η=ðπσ 3 Þ. The explicit expression for the SSF is obtained by solving the integrals defined in Equations (4) and (5). These results are substituted in Equation (3), and after some mathematical manipulation, we obtain that the SSF may be written as If we neglect all the quantities in the right-hand side of Equation (32) except those having the subindex 'HS', which are given bỹ we will be left with the expression for the SSF in the Percus-Yevick (PY) approximation for a hard-sphere system. [23] This result may be obtained in the limit K n ð Þ ! 0 for all n. In addition, χ ¼ k=2.
On the other hand, the two remaining terms in Equation (32), expressed by the summation symbols, contain information corresponding to the Yukawa closure of n = 1, 2, . . ., M terms. There, The coefficient Q n ð Þ χ ð Þ appearing in the above equations is defined as Also, where j 0 (x) and j 1 (x) are the spherical Bessel functions of order zero and one, respectively: Our explicit result for the SSF is more general than that obtained by Liu et al., [15] because they considered only two Yukawas in the closure relation. Furthermore, our result depends on the solution for the scaling matrix Γ given in Equation (26) and the symmetry relations (Equations (27) and (28)). In Appendix B, we will provide a brief explanation of the algorithm used to obtain this scaling parameter. It is worth mentioning that, at least in principle, the general expression for the SSF of mixtures with two Yukawa-type terms given by Cruz-Vera et al. [13] should be reduced to our result if one takes a simple liquid with M = 2.

Reduction to the case of one Yukawa
In the case of one Yukawa, we take M = 1 as the upper limit of the summations in Equation (32). It may be noticed that some of the coefficients will be appreciably simplified. For example, the scaling relation, Equation (20), reduces to where we have used that Γ 11 ¼ γ=σ and X ¼ X ð1Þ . The equation defining the coefficient a (1) now reads as Additionally, the Yukawa closure, given by Equation (26), reduces to which is equal to the relation given by Ginoza for the one-Yukawa closure, [24] with K ¼ K ð1Þ =σ, z ¼ z 1 σ, and where Other quantities appearing in the expression for the SSF reduce to the following: As it may be noticed, the SSF for one Yukawa has the same algebraic structure as in the general case, but some of the quantities appearing in its expression are simplified considerably and now depend on the solution for the scaling parameter γ satisfying Equation (42). The physical solution for this parameter is quickly obtained by the fixed-point iterative method assuming γ ¼ 0 (if zÞ0) as initial value for the iteration. The result obtained for the SSF by analytical reduction to the case of one Yukawa is equivalent to the expressions reported by Herrera et al. [25] and by Arai. [26] 4. Results and discussion for the static structure factor Let us start by mentioning that the MSA is asymptotically exact at high densities; however, it does not satisfy the asymptotic limit of low densities. This deficiency can be corrected by Yukawa closure through a thermodynamic consistency criterion, such correction is called the GMSA. As a first attempt, a method to improve the HS model within the PY approximation [23] is due to Waisman parametrization. [8] Therefore, one of the initial motivations of this work has been to find the parameters that correct such a deficiency. Waisman scheme is based on the fact that the Carnahan-Starling equation of state describes appropriately the data simulation for the hardsphere model. Following Waisman, we will illustrate the correction to the system of hard spheres by using the equations for one Yukawa (see Appendix A for details). The corresponding explicit parameters K η ð Þ and z η ð Þ of the single Yukawa are calculated according to Equations (A5) and (A6), respectively. Figure 1 shows our results for the RDF compared to Monte Carlo (MC) simulation data [27,28] and with integral equation theories. To obtain the RDF we use Equation (A1) given in Ref. [29]. Figure 1, we see that there is good agreement with MC simulation and the Rogers-Young (RY) closure. [30] In general, we observe that the behavior of the RDF obtained through of the parameters that correct the deficiencies of the PY approximation, are in good agreement with Barker and Henderson's simulation data [31] and the RY closure. It should be noted that for the hard-sphere model, the GMSA is more accurate than the results of the RY when both are compared with simulation results. This is because in the first, the equation of state and the contact value of the RDF of the MSA match the corresponding equations of Carnahan-Starling (CS), while the second is thermodynamically self-consistent, but excluding the CS equation.
The idea to obtain a reliable reference system for the hard-sphere model is to improve the prediction at high densities considering some perturbative theory, or also to calculate the bridge function for the HS system in the liquid state integral equation theory. Alternatively, it may be used in the liquid-glass transition of the self-consistent generalized Langevin equation theory of colloidal dynamics. [32] In recent years, there has been a considerable number of investigations of colloidal dispersion having both a short-range attraction and long-range screened Coulomb repulsion. [33][34][35][36] Such type of interaction can be found in many charged colloidal system, such as proteins, globular proteins, and weakly charged colloids in a low dielectric constant solvent with added nonadsorbing polymers to induce an attraction depletion mechanism. Recently, SANS and small angle X-ray scattering experiments reveal new peaks in the structure of colloidal and lysozyme systems. [15,17,[35][36][37][38] For these systems, has been shown the existence of equilibrium cluster, where the low-k peak was attributed to a cluster-cluster interaction. [37,38] However, the presence of a scattering peak at a low-k value does not directly indicate the presence of a cluster state as commonly defined. Recently, this low-k peak has been related with the formation of IRO between scattering centers in a sample, of which, a structured fluid consisting of clusters of particles is only one specific scenario, [35] and therefore, the IRO peak not merely has the feature to indicate the correlation between clusters particles, rather, contains information on the correlation between monomers and the particles in cluster as well as their cross-correlation, [33] and has no direct relation with the formation of equilibrium clusters in a solution. [33,35,36] In the integral equation theory of liquid state, the SSF in the MSA with two Yukawa terms slightly underestimate the height of the monomers peak and the IRO peak [34]; therefore, it has been increasingly evident that double Yukawa MSA closure cannot account for these observations, [34] and accordingly three or more terms are needed. Broccio et al. [34] used the two-Yukawa fluid to model the effective interaction among lyophobic colloid and globular proteins in solution. The two-Yukawa fluid is composed of particles interacting with the potential where K (1) and K (2) are the strengths of attraction and repulsion, respectively. The two-Yukawa model can be used to simulate the classic Derjaguin-Landau-Verwey-Overbeek (DLVO) potential, where the short-range attraction usually represents an approximation of the van der Waals force or entropic forces, whereas the long-range repulsion, to infinite dilution, represents the repulsive part the DLVO potential. We use a third Yukawa-type term in Equation (2) to be accurate with the hypernetted-chain (HNC) closure [34,39] or MC simulation data [34]; this third term does not represent another  Figure 1. RDF for the hard-sphere liquid for a reduced density η ¼ 0:484. The solid black curve represents the results improving the PY solution, and stars stand for RY closure data. The solid circles stand for MC simulation data. [27] In the inset, we plot the RDF for a reduce density η ¼ 0:4, here the solid circle are the MC simulation data from Ref. [28].
term in the pair potential given by Equation (47). Therefore, we will take M = 3 in Equation (2), and the potential parameters values given by Broccio et al. at reduced density η ¼ 0:15, where the third term in the Yukawa closure is used to modulate the height of the IRO peak, and slightly the height of the monomers peak. Figure 2 shows our results for the SSF compared to HNC closure data. This figure shows only the case in which the double Yukawa fails to reproduce quantitatively the low-k peak (see Figures 3(c), 3(d) and 5(a) of Ref. [34]). The used parameters from bottom-up are z 1 ¼ 8, K 1 ð Þ =σ ¼ 2:5, z 2 ¼ 0:5, K ð2Þ =σ ¼ 0:625, and the third term is z 3 ¼ 4:5, K 3 ð Þ =σ ¼ 1:25 refers to attractive strength. The same parameters for the inverse range of attraction and repulsion are used for the curve shown in the middle of Figure 2, except that the strength are K 1 ð Þ =σ ¼ 3, K ð2Þ =σ ¼ 0:75, and K ð3Þ =σ ¼ 1:5. For the last curve at the top, the parameters are z 1 ¼ 2, K ð1Þ =σ ¼ 1, z 2 ¼ 0:5, K ð2Þ =σ ¼ 0:1, and the third term with parameters z 3 ¼ 4:5, K ð3Þ =σ ¼ 0:4, refers to attractive strength. We have noticed that the inclusion of the third term does not change the position of the IRO peak for different values of the volume fraction. The values of the third Yukawa term in the MSA solution have been fixed to achieve a good agreement with HNC solution; therefore, K (3) and z 3 are treated as free parameters and can be adjusted to predict the structure of more accurate closure, e.g. reversed the hybridized mean spherical approximation, Zerah-Hansen, and Bemont-Bretonnet. [40][41][42] In the optimal case, the parameters of the third Yukawa term, in Equation (2), can be obtained numerically by imposing the thermodynamic selfconsistency of the theory.
We used a third Yukawa term to improve the deficiency that is displayed by double Yukawa MSA closure to predict the IRO peak. The values of third Yukawa term have been set so as to achieve an optimal fit with computer simulation data given by Godfrin et al. [33] Our SSF compared with these data is shown in Figure 3. We considered, in particular, the potential parameters given by these authors: z 1 ¼ 10, K ð1Þ =σ ¼ 4:444, z 2 ¼ 0:5, K ð2Þ =σ ¼ 0:4444 for the left panel and z 1 ¼ 10, K 1 ð Þ =σ ¼ 2:415, z 2 ¼ 0:5, K ð2Þ =σ ¼ 0:2415 for the right panel. In all the cases shown in this figure, the inverse range of the third term in the Yukawa closure (Equation (2)) is z 3 ¼ 15, while the values of the reduced density and strengths are expressed in this figure. It may be seen from the figure that the three-Yukawa fluid model performs well.
Godfrin et al. have concluded that the location of the IRO peak is, in general, not an accurate representation of inter-cluster spacing and thus from the preferred cluster size. [33] In particular, the Monte Carlo simulation of dispersions identifies four clearly different microstructures that exhibit an IRO peak: monomer dominated states, clustered, clustered percolated, and random percolated system. Thus, when the concentration is increased and the system is at a reduced temperature of T Ã ¼ 0:25, the magnitude of the IRO peak is an indication of the transition from a molecular fluid to a cluster fluid and then to a cluster percolated state (shown in the left panel of Figure 3 from bottom to top, respectively). However, when the system is at a T Ã ¼ 0:46 and the concentration increases, the fluid goes from monomer-dominated molecular fluid to random percolated states (shown in the right panel of Figure 3 from bottom to top, respectively). To understand better a system with IRO, it is very useful to calculate the potential of mean force βw r ð Þ ¼ À ln g r ð Þ ½ : (48) Figure 4 shows our results for the potential of mean force and the RDF. As can be seen from the figure, the mean force between the particles is strongly attractive at short separations (first local maximum from the left), which becomes less important with increasing the temperature (states B, F, and E). On the other hand, the potential of mean force decreases in magnitude and becomes repulsive at larger separations, as a consequence of the low value of S k ! 0 ð Þshown in Figure 3. The results of MC simulation given by Godfrin et al. have pointed that the IRO of states A and B arises from the localization of monomers; this conclusion, to some extent, can be supported by the behavior of the potential of mean force or the RDF. We note that the structural order of state B is lacking at r=σ>2, in contrast, the structural order of state A extends slightly beyond these distances. At short distances, we also note that both the magnitude and range of depletion attractive of state B are smaller than those of state A. These differences help us to distinguish the magnitude of the IRO peak between both states, while repulsion provides a stabilizing mechanism at larger separations; the particles of state A that can overcome the barriers at r=σ % 2:2 and at r=σ % 1:4 and fall into an attractive potential well are more bound than those of state B, because the short-range attractions between monomers particles become more dominant at reduced temperature of T Ã ¼ 0:25. Therefore, this mechanism explains why the maximum of the depletion region moves to higher distances, and also explains the enhancement of values of the RDF at contact, for state A.
In the case of states E and F, Godfrin et al. have concluded that the IRO peak is due to random percolate structure from the aggregation of monomers and small polydisperse clusters that eventually span the system. The RDFs of these states are almost similar to RDF of state B, except by broad peak at r=σ % 2 that indicates that there is a slight depletion zone as shown in the potential of mean force, probably due to the increase of the volume fraction, forcing all particles staying close to each other or to formation of a cluster. The magnitude of the IRO peak is smaller than the other equilibrium states due to high temperature and increasing volume fraction, causing a greater structural disorder on the IRO length scale.
The RDF of states C and D shows a structural order quantitatively different from other states. This structural order are shown by a pronounced minimum near of r=σ ¼ 1:35 and r=σ ¼ 2:25, increased considerably the repulsive barrier in the potential of mean force that resulting in depletion zones, due that the system has a low temperature and an increasing the volume fraction. Thus, the magnitude of IRO peak is increased, denoting that the depletion zone at short separation and at r=σ % 2 increases at high volume fraction. In Ref. [33], the authors have mentioned that the states C and D are equilibrium cluster; the clusters are formed at low temperature and become more abundant with increasing volume fraction. They hypothesize that the existence of clustered states can be characterized by an IRO peak of magnitude higher than three. Also, they concluded that state C is a cluster state, while state D is a cluster percolate state. [33] As it may be seen, the potential of mean force (or the RDF) for both states shows differences, particularly, at short separations the state C shows a double depletion attraction region, extending to a distance beyond the diameter of the reference particle. Therefore, cluster results from the delicate balance of a short-range attraction and a long-range repulsion. The attraction among particles favors the formation of cluster, while bulk aggregation is prevented by the long-range repulsion which increases with cluster size due to the increasing charge of the cluster (repulsive part) that provides a stabilization mechanism. [36] However, the increasing concentration causes transition of the cluster state (C) to percolated state (D). [33] The potential of mean force of state D (cluster percolated state) has the feature of being the most repulsive at short and middle distances. Therefore,  Figure 4. Potential of the mean force and the RDF obtained from MSA analytical expression using three Yukawa terms. The parameters and the symbols are same as in Figure 3.
particles or the cluster particles of state D can overcome the repulsive barriers and fall into an attractive potential well to stay together. Cluster particles within this well can only hardly go out, and presumably, at this lower temperature, the long-range repulsion leads to a preferred length scale that sets the microstructure of the percolate state in a system with IRO. Therefore, the structure with large oscillation as indicated by the first maximum of the RDF is responsible for the formation the IRO peak in the states C and D. As a result, cluster-cluster correlations become more and more important and eventually they induce an arrest transition in materials like gels. [33,36]

Conclusions
In this work, we have obtained an analytical expression for the SSF, Equation (32), of a onecomponent fluid using a closure of M Yukawa terms within the MSA. This expression is a generalization of the SSF obtained by Liu et al. [15], but unlike these authors, our expression is written in terms of a scaling parameter matrix. The numerical values of this matrix are determined through an efficient iterative method for solving the set of nonlinear equations given by the closure condition, Equation (26); see Appendix B.
From our expression for the SSF, the solution of the two-Yukawa fluid progressively fails to reproduce the low-k portion when increasing the attraction strength. [34] By using three Yukawa tails, we have improved these deficiencies and we have found good agreement with the HNC closure and Monte Carlo simulation data. Therefore, another important objective of this article was to present a formulation explicit and simple, such that the Yukawa parameters of the closure relation allow us to go beyond the MSA approximation, and thus represent some kind of system, such as lyophobic colloids or globular proteins in solution.
In summary, we believe that the combination of the MSA with appropriate parameters is able to be compared with another closure or computer simulation, or (in some cases) with experiments in certain regions of the thermodynamic state space. Therefore, the Yukawa closure may be taken as a convenient reference model in perturbative studies in the theory of liquids, as well as in the study of liquid metals or in the ideal glass transition, which we will discuss in a future work.

Disclosure statement
No potential conflict of interest was reported by the authors.