The multi-density integral equation approach in the theory of complex liquids

Recent development of the multi-density integral equation approach and its application to the statistical mechanical modelling of a different type of association and clusterization in liquids are reviewed. The effects of pairing, polymerization, solvation, formation of the network bonds and selfassembling are discussed. The numerical and analytical solutions of the integral equations in the multi-density formalism are used for the description of the association phenomena in the electrolyte and polyelectrolyte solutions, water, polymers, microemulsions and other fluids. The application of the multi-density integral equation approach for the treatment of the percolation phenomena, the adsorption of fluids in porous media and the description of electronic structure of associative fluids are illustrated.


Introduction and background
For the last decades the integral equation methods have been intensively used in the modern liquid state theory.This technique is based on the analytical or numerical calculation of the pair distribution function g(r) by the solution of the Ornstein-Zernike (OZ) equation within different closures: the Percus-Yevick (PY) approximation, the hypernetted chain (HNC) one, the mean spherical approximation (MSA) and its different modifications [1].It was shown that such closures are efficient enough for the fluids having not so strong interparticle attraction but they break down with the increase of the interparticle attraction which leads to the clustering of atoms or molecules into pairs and into larger groups such as chains, networks, selfassembling aggregates etc. Due to the clustering the δ-like intraparticle distribution function appears in the distribution function which can be divided into the intra and interparticle parts g(r) = g intra (r) + g inter (r). (1) The situation will be more clear if we use the function of the corresponding running integration numbers.
which describes the average number of particles in the sphere of radius R surrounding one of them which is found in the center of this sphere, ρ is the number density of particles.Resulting from clusterization, n(R) divides into bounded (intra) and nonbounded (inter) parts n(R) = n intra (R) + n inter (R). ( Due to the saturation of bounding where the number of bonds n bond is fixed.Specifically, for pairs n bond = 1, for chains n bond = 2, for network n bond = 4 etc.The background of the traditional closure relations in the integral equations technique is connected with diagram analysis of the Mayer density expansions of the pair distribution function which are not adapted correctly for the description of the clusterization phenomena.The clusterization is caused by the attraction part of interparticle interaction and to describe their contribution it is more convenient to use the activity expansions.In particular, in order to reproduce the correct lowdensity limit for the fluids with strong clusterization, an infinite number of terms in the density expansion must be included, while only a few terms of the activity expansion are sufficient for this purpose [2,3].Consider, for example, the series in the activity Z for the pressure P and density ρ terminated at a second order terms where b 2 is the attraction part of the second virial coefficient and β = 1/(kT ).Elimination of Z from the equations (5) yields The more familiar corresponding virial equation of state is obtained from (6) by expanding in the power series of ρ For the systems with weak attraction b 2 → 0 and both expressions ( 6) and ( 7) reproduce the ideal gas equation of state βP/ρ = 1.However in the limit of strong attraction between pairs, b 2 → ∞ and correct equation of state for an ideal gas of pairs βP/ρ = 1/2 follows directly from equation (6) and in order to obtain this result from equation (7) the infinite number of terms should be included.
In a more general case where b n is the attraction part of the nth virial coefficient.Hence, the strong interparticle attraction Z 1 and in the limit After elimination of Z from equations ( 8) we obtain the equation of state which will be changed from the ideal gas equation βP/ρ = 1 to the equation of ideal gas of pairs βP/ρ = 1/2, ideal gas of trymers βP/ρ = 1/3,..., and in general the ideal gas of n-mers βP/ρ = 1/n.Moreover the summation of the infinite series in (8) leads to the possibility of the self-assembling in the system.In approximation b n ≈ b n−1 2 the expansions (8) reduce to the following form which is well known in the thermodynamical theory of micellization [4].For this reason Z = ρ 0 can be identified with the monomer density of amphiphilic molecules and the divergence point Z c = ρ c = 1/b 2 has the sense of the critical micelle concentration (CMC).Considering this, it might be expected that a theory which combines the activity and density expansions would be advantageous.The second and higher terms in the expansion (8) can be interpreted as the dimers density ρ 1 = b 2 Z 2 ,..., n-mers density ρ n−1 = b n Z n correspondly.Such interpretation suggests the possibility of the description of clusterization by introducing the multi-density formalism for this purpose.Recently a consistent integral equation theory for the description of the clusterization in liquids has been proposed [5,6].This is based on the multi-density formalism in which the description in terms of the activity and density expansions are combined.
The theory has been applied to describe clusterization of various type in the fluids modelled by hard spheres with a number of off-centre square-well bonding sites embedded into the hard-core region.The multi-density formalism was reformulated in order to treat the effects of clusterization in fluids with spherically symmetric attraction and it was applied for ionic systems [3,[7][8][9][10][11][12][13].The multi-density integral equation approach was used to describe the chain [14][15][16] and network [17] formation fluids, molecular and macromolecular liquids [18], the fluids with the multi-arm star polymerization [19] etc.The multi-density formalism was also adopted for the treatment of the percolation phenomena in the network forming fluids [20], for the adsorption of associative fluids in porous media [21,22] and for the description of electronic structure of associative fluids [23].
The recent progress in the application of the multi-density integral equation approach in the theory of complex liquids will be reviewed in this report.The general scheme of this approach in the framework of the two-density formalism including the applications to ionic systems will be considered in the second section.In the third section the possibility of the multi-density formalism for the description of polymerization and network formation are presented.Some applications and conclusions will be considered in the fourth section.

Two-density approaches: The application to ionic systems
The general idea of the multi-density formalism is connected with the separation of the potential of interparticle interaction U(r) into the bonding and non-bonding parts where U bond (r) is some short-range attractive part of interaction and includes at least the potential energy minimum of U(r).The nonbonding part U nonb (r) includes a repulsion part and a long-range tail of U(r).The diagrams appearing in the activity expansions of the one-point density are classified with respect to the number of associative bonds incident with the labelled white circle.Thus, the total number density of the system is separated into two densities, the density of nonbonding particles (monomers) ρ 0 and the density of bonding particles ρ 1 : Similarly the pair distribution function will be splitting into four terms.
where x = ρ 0 /ρ is the fraction of a nonbonding part of particles.Due to saturation of bonding the restriction only by pair formation leads to the self-consistent relation for x where The classification and topological reduction of the diagrams for pair correlation function leads to a generalized version of the OZ equation h(r 12 ) = C(r 12 ) + dr 3 C(r 13 )ρh(r 32 ), (16) where the corresponding matrices have the following form: As usual, the OZ-like equation should be supplemented by closure relations.Among them we should distinguish the associative HNC (AHNC) closure g 00 = e −βU nonb +t 00 , g 01 = g 00 t 01 , g 10 = g 00 t 10 , g 11 = g 00 [t 10 t 01 +t 11 +f as ], (18) where t αβ = h αβ − C αβ , the associative PY (APY) closure g 00 = e −βU nonb y 00 , g 01 = e −βU nonb y 01 , g 10 = e −βU nonb y 10 , g 11 = e −βU nonb (y 11 + f as y 00 ), where y αβ = g αβ − C αβ , the associative MSA (AMSA) closure where d is the diameter of particles.
The analytical solution of OZ-like equation for the dimerizing hard spheres in APY approximation was done in [24][25][26].Of special interest is the application of the two-density formalism to ionic systems.The numerical solution of AHNC for 2-2 aqueous electrolyte solution was obtained in [3].The different version of analytical solution AMSA for ionic systems is done in [8][9][10]27].It was shown that AHNC and AMSA essentially improve ordinary HNC and MSA.
In [28] in the framework of AMSA to calculate the nonbonding fraction x in the equation ( 14) an exponential approximation was introduced for where g nonb (r) is the pair distribution function of the reference system with the interparticle interaction U nonb (r), The approximation (21) for very dilute ionic solutions is consistent with the intuitive idea of Bjerum [29] to introduce the concept of ionic pair to improve the Debye-Huckel [30] treatment of ionic systems when interionic interaction increases.It was shown [28] that the electrostatic part of the excess thermodynamic properties has the same form, as in ordinary MSA, but the screening parameter Γ should be changed to a new one Γ B .The exponential approximation (21) was also used [31] for the extension of the Hoye-Stell scheme [32] of the calculation of thermodynamic properties in the framework of the two-density formalism for the fluids with arbitrary interparticle long-range interaction.Recently it was proposed in equation ( 14) for the calculation of the nonbonding fraction x together with approximation (21) for g 00 (r) for f as (r) instead of (15) to put [33] f as (r) = 1 − e βU attr (r) 1 − βU attr (r) + 1 2 where instead of the separation ( 11) we divide the total potential into the repulsion and attractive parts The choice of f as (r) in the form (23) according to the Ebeling idea [34,35] fixes the second virial coefficient of the system and gives the correct description of thermodynamic properties at least for the small concentrations.Such a scheme has been used to describe the thermodynamical properties of nonaqueous electrolyte solutions [33].
The generalization of AMSA results applied to ion-dipole systems has been considered [36,37].
The specific situation occurs in the case of the solutions of spherical polyelectrolytes, charge stabilized colloids, micelles and globular proteins which can be viewed as highly asymmetric electrolytes, with large and highly charged polyions and small ordinary counterions.Kalyuzhnyi and Vlachy [11,12] take into account the fact that due to the high asymmetry in size and charge the counterions bind to a limited number of polyions, while the polyions can bind to an arbitrary number of counterions.As a result, the density of counterions ρ c is separated into the nonbonding and bonding parts while the density of polyions ρ p is nonbonding.Similar to (13) for the pair distribution functions we have where for the fraction of nonbonding part of counterions x c we have the equations The generalized version of OZ equation in the case considered has the form similar to (16), where now only the matrices h cc , C cc and ρ c have forms similar to (17) while As for the closure relations we have forms similar to ( 18)-( 20) ), g 00 cc (r) = e −βUcc(r)+t 00 cc , g 01 cc (r) = g 00 cc (r)t 01 cp , g 11 cc (r) = g 00 cc (r)t 11 cc (r), (29) in AHNC approximation, ), g 00 cc (r) = e −βUcc(r) y 00 cc (r), g 01 cc (r) = e −βUcc(r) y 01 cc (r), g 11 cc (r) = e −βUcc(r) y 11 cc (r), (30) in the APY approximation, in the AMSA approximation.
The numerical solution of AHNC equation for highly asymmetrical ionic systems has been solved in [11,12].It was demonstrated that the AHNC yields accurate structural and thermodynamical predictions for a wide range of states, including the region where the ordinary HNC approximation does not give a convergent solution.The analytical solution of AMSA equation for highly asymmetrical ionic systems has been obtained in [13].
The formalism developed for the highly asymmetrical ionic systems can also be useful for describing the ionic solvation in ion-molecular systems [38].

The multi-density formalism
For the particles having more than one bonding state, the formation of chains, rings, networks and more complex aggregates is possible.Such aggregates can be considered as a collection of monomers (segment) bonded at asymmetric attraction sites.The multi-density formalism is needed for the description of such fluids.In general for the particles with M bonding sites the density is separated into 2 M densities of different bonding states.The diagram analysis leads to the generalized version of the OZ equation which has the form similar to (16) where in the general case h, C, ρ are the matrices 2 M × 2 M .In general 2 M − 1 self-consistent relations are needed instead of the relation (14) for the pairing case since the total number density of the system is separated into 2 M densities of bonding states.Some simplification can be connected with the approximation that the bond creation between two molecules is independent of the existence of other bonds.As a result, the fraction of the molecules which have n < M bonded neighbours can be given by the binomial distribution where p = 1 − α, α is the fraction of particles nonbonded by one fixed site.
For the molecules with two attractive sites A and B (one is donor, the other is acceptor) the assumption (33) leads to the ideal chain approximation (ICA) [14] and α = 1/m, where m can be considered as the mean chain length.For this case the matrices h, C and ρ have the following forms where the lower indices α denote the unbonded (α = 0), singly bonded (α = G) and doubly bonded (α = Γ) states of the corresponding particles.Similar to (13) the pair distribution function where in ICA approximation g 0Γ (r) = g AΓ (r) = g BΓ (r) = g ΓΓ (r) = 0.The analytical solution of the OZ-like equation in polymer PY (PPY) approximation for the chain forming fluids in ICA approximation was obtained and discussed in [14,16].In the ICA approximation the formation of the ring polymers is neglected.This approximation can be used to describe a system of chain polymers, polydisperse in length, that is characterized by a prescribed mean chain length.Some other approach in the framework of PPY approximation is used in [15,39].These studies are focused on the investigation of the polymer fluids represented by the system of monomers that have completely associated into polymers of a fixed size, usually into linear freely-joined chain polymers of a fixed length.Such an approach makes it possible to consider the system at all the degrees of polymerization, including the limit of completely dissociated and completely associated systems.The latter limit enables us to describe molecular and macromolecular liquids in the multidensity approach [18].In [40,41] the multi-density formalism was applied for the polyelectrolytes, with the flexible linear polyion chains.It was shown again that the electrostatic part of the excess thermodynamic properties has the same form, as the ordinary MSA, but screening parameter Γ should be changed to the parameter Γ B for the chain polyelectrolytes [42].The combination of polymer multi-density approach with the approach that was used in the theory of spherical polyelectrolytes (see equations (25)(26)(27)(28)) enables us to describe the fluids of multi-arm polymer star molecules [19].
For the molecules with four attractive sites A, B, C and D (the two are donors and the two are acceptors) the assumption (33) leads to the ideal network approximation (INA) [17].The matrices h, C and ρ can be written in the form similar to (34) with the difference that every matrix element is a submatrice.For example The analytical solution of the OZ-like equation for the network forming fluids in network PY (NPY) and INA approximation was done in [17,43].It was shown [16,17] that the calculated structure factors for chain and network systems exhibits a peculiarity (a so-called pre-peak) at small wave numbers connected with the forming of relatively large molecular aggregates.It was shown that such a peculiarity plays an important role in describing the structure factor of liquid sulfur [44].The possibility of the existence of a prepeak in the structure factor of water is discussed based on the combination of computer simulation and integral equation theory [45].Lately, such a peculiarity has been detected in many polymer and glassing substances and in ionic liquids [46], where, due to the correlations between relatively large clusters, the ordering in mesoscopic scale can appear.
The general analytical solution of the OZ-like equation in the multi-density formalism in PY approximation for the mixture of associating hard sphere was done in [HP].This result [47] was generalized for the ionic fluids [48].
The number of the bonding states of molecules can principally change the thermodynamic properties of fluids.For illustration we present here the generalization of Van der Waals equation of state for the fluids of molecules with M bonding sites.The expression for pressure P for nonionic fluids can be represented in the following form where βP HS /ρ is the hard sphere contribution, βP MAL /ρ is the contribution connected with clusterization (the mass action law), βP atr /ρ is the contribution connected with Van de Waals attraction.To describe βP MAL /ρ contribution we make use of the expression obtained by thermodynamic perturbation theory [6], where y HS (1) is the contact value of the cavity distribution function y(r) = e βU (r) g(r) for the hard sphere model.We put the term βP atr /ρ in Van der Waals form where η = 1 6 πρd 3 is the hard-sphere packing fraction, W is the constant characterized Van der Wals attraction.
After the density expansion If we suppose that this series creates the geometrical progression, the equation of state can be represented in the form For the pairing case M = 1 and in the limit α → 0 For the chain formation case M = 2, α = m −1 and with increasing chain length T c increases and η c decreases.Such a behaviour is qualitatively in agreement with n-alkanes and n-perfluoroalkanes [49].In the limit α → 0 The specific situation appears for the network formation fluids (M=4).We can see that in the case M > 2 the limit α → 0 in (42) is not possible.It means that some new phenomena appear for M > 2 which is connected with the possibility of forming the infinite cluster in case M > 2. Hence, a new mechanism of the criticality appears in which the "ideal" term plays a principal role and has sense only for α 1 − 2/M.The criticality of the system is connected with the network formation and in this regime the Van der Waals term (39) has no influence.More sophisticated treatment [50] does not change the described behaviour qualitatively.The careful analysis shows that the critical value of the compressibility factor Z c = β c P c /ρ c changes from Z c ≈ 1/3 for the Van der Waals fluid to Z c ≈ 1/5 for the network formation fluid which is very close to the value Z c ≈ 0.229 for the real water [51].
The multi-density integral equation theory was reformulated for studying connectedness properties [52] in order to understand such a peculiarity of the network forming fluids.The division of the potential of interparticle interaction into the blocked and connectedness parts leads to the similar separation also for the pair and direct correlation functions The connectedness pair and direct correlation functions satisfy the OZ equation similar to (16).The mean cluster size is given by As the percolation transition is approached S increases and becomes infinite at the percolation threshold.The connectedness version of the OZ equation supplemented by the NPY-like closure and INA approximation was solved analytically [52].In these approximations It is seen that S → ∞ when α → 1 3 .From the analysis of the spinodal and percolation curves it was shown that liquid phase including the critical point is inside the percolation region.Such a conclusion is in accord with the recent computer simulation for supercritical water [53].

Some applications and conclusions
The characteristic features of numerous complex liquids are connected with associating the molecules into different clusters caused by the strong interparticle attraction.The starting point of the theory of such liquids is the combined cluster expansions for pair correlation function in which the activity expansions are used to describe the contribution of the bonding part of the interparticle interactions while the usual density expansion is used to describe a nonbonding part of interactions.The diagram analysis of these cluster expansions leads to the multi-density integral equation approach which is flexible enough to treat different associative features of liquids such as dimerization, polymerization, network formation, solvation, self-assembling etc.
The possibilities of the theory developed are tested by comparing with computer simulations.It is shown that the multi-density approach essentially improves the integral equation theory for ionic systems.The theory yields very good agreement with computer simulations including the region of the parameters where usual approximations like HNC haven't got a convergent solution.In the limited small ionic concentration, the result of AMSA is consistent with the original Bjerrum theory for a weakly associative case and takes into account the electrostatic contributions from ionic pairs which are usually neglected in the Bjerrum theory and are especially important for a strongly associative case.The multi-density integral equations are also solved for polymerization and network formation cases and are reformulated for studying the connectedness properties of network forming fluids.The gas-liquid critical point is predicted to exist for network case, and a region of liquid state is inside the percolation region.
The developed multi-density formalism was applied for describing complex liquids at the solid-associative fluid interfaces.To this end, the associative analogue of the Henderson-Abraham-Barker approach [54] was formulated [55] and it was applied to describe dymerizing fluids [55][56][57], polymers [58,59], the network forming fluids [60] close to a hard and adhesive surfaces and ionic fluids close to a charged wall [61].The developed approach was also applied to describe the adsorption of polymers [58,62,63] and network forming fluids [64] onto crystalline surfaces.It was shown that the intramolecular correlations in polymer chains lead to the cooperative increase in the surface coverage which takes place under the conditions of a strongly diluted case.This effect was recently observed experimentally using contact-angle measurements of polymer adsorption [63].For the network forming fluids a new type of the cooperative adsorption was discovered which is related to a bridging due to tree-like clusters adsorption.To describe the adsorption of associating fluids in porous media the associative analogues of the replica OZ equations [65] were formulated [21,22].Since in the presence of porous media the value of liquid-gas critical temperature increases and a coexistence curve is narrower compared to the pure fluids [66] the role of association effects in the presence of porous media increases and it can be very important for different fluids, including simple liquids.
The multi-density formalism was also adopted to describe electronic structure of associative fluids [23].The approach developed is closely connected with the single superchain-effective medium approximation [67], in which the effects of association are explicitly built-in.The theory has been tested in a simple-minded model of alkali metal, namely a hard sphere fluid with a one-level tight-binding Hamiltonian with transfer matrix elements modeled by Yukawa terms.It was shown that this treatment correctly accounts for the extreme modifications that the association induces on the band structure.If one assigns one electron per atom basis function, like one would have in alkali metals, the dimerization process induces a metal-insulator transition that is correctly described by the theory.
In [68,69] a pairing sticky hard sphere model was introduced to describe the influence of proteins on the physical properties of reverse micelles.In this model it was assumed that the empty micelles are described by the sticky hard sphere model and the presence of proteins is represented by an association parameter that characterized the pair formation micelles.The association parameter was determined by fitting the experimental structure factors and from this parameter the decrease in the percolation threshold and the change in the critical parameters induced by proteins have been estimated.Thus, it has been possible to reproduce the main effects induced by the solubilization of proteins in reverse micelles.
The multi-density integral equation approach is also useful for modelling the effects of cation hydrolysis and polynuclear ion formation in the theory of aqueous electrolyte solutions [70].The cation hydrolysis is caused by the increase of bonding between cation and oxygen of water molecules with the decrease of ionic radius and/or the increase of their charges.At the same time the repulsion between cation and hydrogen of water molecules increases and some part of water molecules in ionic hydration shell can dissociate.It means that instead of hydrated ions [M(H 2 O) n ] Z+ , there appear the hydrated complexes [MO n H 2n−h ] (Z−h)+ , where h is a molar ration of the hydrolysis.The effect of cation hydrolysis was studied recently [71] in the framework of the reference HNC approach and using the central force model for water.The ionic hydration number for oxygen and hydrogen was calculated ( The molar ration of the hydrolysis was defined as For rigid molecules h = 0.It was shown that with the increase of ionic charge the constant h increases.The possibility of computer simulation of hydrolysis effects is discussed in [72].The effects of polynuclear ion formation in aqueous solutions of metal salts in the framework of ionic approach was mimiced by the intercationic associative interaction [70].The theory predicts a significant influence of the polynuclear ion formation on the thermodynamics and structure properties of electrolyte solutions in a wide range of concentrations including a very dilute ecologically important concentration region.

PACS: 05.20.-y
generalized form of the Van der Vaals equation.The simple estimation shows that in result clusterization the critical temperature T c increases and the critical density η c decreases compared with the corresponding values T V dW c and η V dW c at the lack of clusterization as