Density anomaly of charged hard spheres of different diameters in a mixture with core-softened model solvent. Monte Carlo simulation results

Very recently the effect of equisized charged hard sphere solutes in a mixture with core-softened fluid model on the structural and thermodynamic anomalies of the system has been explored in detail by using Monte Carlo simulations and integral equations theory [J. Chem. Phys., 2012, 137, 244502]. Our objective of the present short work is to complement this study by considering univalent ions of unequal diameters in a mixture with the same soft-core fluid model. Specifically, we are interested in the analysis of changes of the temperature of maximum density (TMD) lines with ion concentration for three model salt solutes, namely sodium chloride, potassium chloride and rubidium chloride models. We resort to Monte Carlo simulations for this purpose. Our discussion also involves the dependences of the pair contribution to excess entropy and of constant volume heat capacity on the temperature of maximum density line. Some examples of the microscopic structure of mixtures in question in terms of pair distribution functions are given in addition.


Introduction
It is our honor and pleasure to make contribution to this special issue dedicated to brilliant scientist Prof. Myroslav Holovko on the occasion of his 70th birthday. During last years we were influenced by his developments and profound results along different lines of research in the theory of liquids and mixtures, interfacial phenomena and phase transitions. We appreciate his mentorship in some projects, many useful pieces of advice and friendly discussions concerning topics of our common interest. During many years Myroslav Holovko explored electrolyte solutions in very detail, see e.g., the classical monograph [1] written by him together with his teacher Prof. I. Yukhnovskyi that summarizes their effort prior to 1980. Later results can be found in many articles and reviews.
The principal objective of the present work is to complete one extension of our very recent study of solutions comprising solute ions and a soft-core model fluid as a solvent [2]. During last decade many studies in the theory of fluids have been focused on single-component, isotropic core-softened models in which the repulsive part of the interparticle interaction potential is peculiar ("softer") at "intermediate" interparticle separations while a hard-core repulsion acts at short distances. This region of peculiar dependence of the pair interparticle interaction is often described by linear or nonlinear ramp [3][4][5][6][7][8], by a shoulder [9,10] or Gaussian term [11][12][13], or another more sophisticated functional dependence, see e.g., recent review [14] and references therein. A more complete list of references has been given in [2] as well.
Unfortunately, the complex shape of the interaction in this type of models requires application of solely numerical procedures for their solution. However, peculiar structural, thermodynamic and dynamic properties relevant to certain real liquids can be reached. Models of single-component fluids with isotropic core-softened potential have been studied by computer simulations in numerous works, see the list of references in [2]. For purposes of the present work it is necessary to mention that the fluids with spherically symmetric core-softened potentials exhibit anomalies in thermodynamic and dynamic properties in close similarity to those observed in liquids with directional inter-particle interactions, such as water, silica, phosphorous, and some liquid metals. One of the examples of such anomalous behavior is the presence of the fluid density maximum as a function of temperature. Moreover, some of the softcore fluid models exhibit puzzling liquid-liquid coexistence envelope terminating at the critical point, in addition to common liquid-vapor phase transition. Thermodynamic and other anomalies can be related to this additional phase transition, see e.g., references [15][16][17]. The soft-core models by design miss the orientational interparticle correlations. Thus, their correspondence to, for example, water or aqueous solutions is restricted to thermodynamic properties. Nevertheless, the soft-core models are simple and permit to capture the peculiarities of thermodynamic properties quite adequately [14].
Theoretical studies of solutions in which a model solvent exhibits thermodynamic properties with anomalies have been attempted recently by different authors, see e.g., [18][19][20][21][22][23][24]. Several important particular issues were explored, namely the solubility of apolar solutes in such solvents, the effect of solute species on the temperature of maximum density line, on heat capacity extrema and others. On the other hand, following our previous theoretical and simulation developments [25,26], we applied similar modelling, but involved ionic solutes rather than nonpolar species [2]. We examined the effect of ionic solutes on the microscopic structure and thermodynamics of solutions with a soft-core solvent that possesses anomalies in one-phase region. Experiments [27][28][29][30] and computer simulations [31][32][33] have shown that certain anomalies of thermodynamic properties are preserved in ionic aqueous solutions for low and up to moderate ion concentrations. They become weaker and finally disappear, if the concentration of ions increases but remains below the solubility limit.
In contrast to the initial stage of this project [2], we consider a mixture of singly charged cations and anions modelled as charged hard spheres and solvent particles described by a soft-core model, all of them having nonequal diameters. The system is studied by the canonical Monte Carlo simulations. The following discussion of the model and methods is brief in order to avoid unnecessary repetitions of [2], only the essential elements are given to benefit the reader.

Model and method
The system in question consists of ionic solute (salt) and solvent. The ionic subsystem is a mixture of positively and negatively charged hard spheres with arbitrary diameters, σ + and σ − , the charge is located at the particle center. Ions and solvent particles are immersed in a dielectric background characterized by dielectric constant ε. The interaction potential between ions i and j is, where e is the elementary charge, z i = |z − | = 1 are the valencies of ions, ε 0 is the permittivity of free space, and r is the inter-particle distance. As in [2], we use the core-softened model by Barros de Oliveira et al. [11][12][13] that possesses thermodynamic, structural, and dynamic anomalies. The potential is a sum of the Lennard-Jones and a Gaussian terms, with ε ss -the solvent-solvent interaction energy and σ s -the diameter of a solvent particle. Parameters a, c, and r 0 are taken as in the original work [11] (a = 5, c = 1, and r 0 /σ s = 0.7). With this choice of parameters, the inter-particle attraction is strongly depressed and thus there is no liquid-gas transition in the solvent subsystem [34].

43607-2
The pressure-temperature curves calculated for the pure solvent at certain densities have a minimum (i.e., (∂P /∂T ) ρ = 0) and consequently (∂ρ/∂T ) P vanishes for some thermodynamic states. The region of density anomaly corresponds to the state points for which (∂ρ/∂T ) P > 0 and is bounded by the locus of points for which the thermal expansion coefficient equals zero.
There is some room for the choice of ion-solvent interaction. However, to be consistent with our previous work [2], the interaction between an ion (i ), and a solvent (s) particle is taken simply in the form of Lennard-Jones (12-6) potential, where ε is is the ion-solvent interaction energy and σ is = (σ i + σ s )/2. The reduced temperature for the solvent subsystem is defined as T * = k B T /ε s , where k B is the Boltzmann constant and T is the absolute temperature.
In order to describe the energetic aspects of interactions between all species, it is necessary to define the parameter α = T * /T * el as the ratio between the reduced temperature for the solvent subsystem (T * ), and the reduced temperature of the ionic subsystem, T * is the Bjerrum length, λ B = 7.14 Å. Similar to the previous study [2], we choose α = 0.25. Thus, if the solvent is even at a rather low temperature, the ionic system remains far from the possible vaporliquid type criticality. The mixing rule for ion-solvent potential is applied such that: We performed Monte Carlo computer simulations in the canonical ensemble and in the cubic simulation box according to common algorithms [35,36]. The simulation was started from a configuration of randomly distributed particles. The number of ions in the box, N ion (N ion = N + + N − ) was chosen to provide their reasonable statistics and along with the desired ion concentration, where N s is the number of solvent particles resulting from N ion and x ion . In addition, we picked up the desired total dimensionless density of the system, ρ * , ρ * = (N + σ 3 where L is the length of the box edge resulting from previously given parameters. Coulomb interactions were taken into account by using the Ewald summation [35,36] with 337 wave vectors and the inverse length equal to 5/L).
In this work, we explore three model salts in the soft-core solvent, namely the NaCl (σ + = 2.21 Å, σ − = 4.42 Å), KCl (σ + = 4.74 Å, σ − = 4.42 Å) and RbCl (σ + = 5.27 Å, σ − = 4.42 Å) see e.g., [37,38]. In all cases, we assume σ s = 3.166 Å in close similarity to the SPCE model of water. Each of the systems was first equilibrated over 10 8 attempted configurations followed by a series of production runs of at least 5 × 10 7 attempted particle displacements to obtain the relevant average properties of the system under investigation (excess internal energy, compressibility factor, heat capacity at constant volume, pair excess entropy, and radial distribution functions). The number of ions present in the simulation box was in the range between 150 and 400 (dependent on the desired density). The excess internal energy of the system was obtained directly as an ensemble average. The equation of state was calculated by [39] PV where U tot is the total excess internal energy of the system and g i j (σ i j ) is the contact value of the appropriate ion-ion radial distribution function. Pair excess entropy was determined from the expression [40] S 2,ex where x i is the number fraction of species i , and g i j (r ) are the pair distribution functions of species i and j . The heat capacity at constant volume was determined from the energy fluctuations as 43607-3

Results and discussion
First series of our calculations concern the NaCl salt model (σ + = 2.21 Å, σ − = 4.42 Å) in the softcore solvent. The cation is rather small and its diameter is smaller compared to the anion and to the diameter of solvent particles. The results for P * (P * = P σ 3 s /ε s ) as a function of T * for a set of densities of the mixture at fixed ion concentration x ion = 0.1 are given by short-dashed lines in figure 1 (a). Each line has a minimum (it is difficult to see it on the scale of such combined figure) marked by a circle. By joining the circles, one can construct the temperature of maximum density line at this particular ion concentration, x ion = 0.1. The P * (T * ) lines at a given ρ * are built by joining points obtained as an average from several runs. Actually, it is quite a tedious procedure. We perform a series of calculations to construct the TMD lines for different x ion , up to x ion = 0.3. Each of the lines has an extremal point that shifts to a higher pressure and lower temperature with an increasing x ion . Moreover, the TMD lines shrink with an increasing ion concentration. These trends were already discussed by us for the model in which all the species involved are of equal diameter [2]. Slightly different view on the changes of the TMD lines is given in figure 1 (b). The shift of the TMD lines to higher ρ * with a growing x ion can be explained by augmenting the excluded volume effects in the system, if ions are added to the solvent. On the other hand, the extremum point of the TMD along the temperature axis goes down due to the overall growth of attraction, principally between oppositely charged ions. In the case of a lower density and lower ion concentration (figure 2), we can attribute the shape of the ion-ion distribution functions to the formation of Na-S-Cl, Cl-S-Cl and Na-S-Na structures or complexes, i.e., of ions separated by a single solvent molecule, corresponding to the local maxima at (σ Na + σ Cl )/2 + w 1 , σ Cl + w 1 and σ Na /2 + w 1 distances, respectively. Here, w 1 is the distance corresponding to the first maximum of the distribution of solvent species, g w w (r ) at 3.63 Å, approximately. Similarly, one can observe trends of formation of clusters at a distance (σ Na +σ Cl )/2+ w 2 , where w 2 is the distance corresponding to the position of the principal maximum of the distribution of solvent species at 7 Å, approximately. Due to a stronger repulsion between small cations, the formation of Na-S-Na structure is less favorable compared to Cl-S-Cl structure at these conditions.
On the other hand, if one considers the model at x ion = 0.3 at a point of the TMD line with slightly higher total density but substantially lower temperature (figure 3) compared to the case described in figure 2, the inter-ion electrostatic correlations are expected to be stronger. Nevertheless, screening of interactions between ions is expected to be more pronounced. Under the conditions given in figure 3, one can see that the contact value of g NaCl (r ) is practically of the same order as in figure 2. However, the height of the shoulder of this function at r ≈ 6.9 Å describing the correlation between oppositely charged ions separated by a single solvent particle is lower than in figure 2. On the other hand, the height of the first maximum of g ClCl is higher as an evidence of weaker repulsion between equally charged anions.
Possible formation of Na-w 1 -Na is also more pronounced, cf. figure 2. The cation-solvent correlations are enhanced in the present case (at this low temperature) compared to anion-solvent correlation and compared to the case discussed in figure 2. We now again return to thermodynamic manifestations of the existence of the TMD lines. A quantitative cumulative measure of the distribution of particles of different species is commonly given in terms of the pair contribution to the excess internal energy per particle as a function of the total density for different compositions of the mixture and at different reduced temperatures. In previous studies of the single-component soft-core model fluids it was suggested that the relation between the density and diffusion anomalies and the microscopic structure of the system can be explored in terms of the excessentropy-based framework, see e.g., [41]. It has been shown that the presence of the density and the dif- fusion anomalies is related to the effect of excess entropy on the fluid density. Our discussion concerns S 2,ex defined for a mixture by equation (5). However, in contrast to previous works, we have calculated the dependence of S 2,ex as a function of the temperature of maximum density rather than the total density of the model. In this respect, our results given in figure 4 (a) correspond to the TMD lines presented in figure 1. For example, the left hand part of the TMD line for x ion = 0.1 in figure 1 (b) that is characterized by dT * /dρ * > 0 corresponds to the upper part (lower absolute values) of the corresponding line for S 2,ex on T * tmd . What can we learn from these functions? First observation is that the S 2,ex (T * tmd ) curves shift to lower values of S 2,ex and lower values for T * tmd with an increasing concentration of ionic solutes. Moreover, these curves shrink with increasing ion concentration. At low values of x ion (and at vanishing x ion as well) the curve S 2,ex (T * tmd ) has a maximum such that dS 2,ex /d(T * tmd ) < 0 in the small region of T * tmd after the point dS 2,ex /d(T * tmd ) = 0. This point does not coincide with the maximum temperature of the TMD line. Discontinuity of this derivative is expected, however, around the point at which T * tmd reaches its maximum. This behavior is much less pronounced at high values for ion concentration in a mixture. The existence of such type of peculiar behavior seems to be worth to investigate in more detail for other soft-core solvent models and for say nonpolar solutes in various solvents for the sake of comparison with the case of ionic solutes.
A rather similar discussion can be presented while considering the dependence of the reduced constant volume heat capacity, C * V on T * tmd , figure 4 The heat capacity was "measured" in simulations through the fluctuations of internal energy according to equation (6). It is quite difficult to evaluate the heat capacity precisely due to fluctuations, especially at low temperatures.
One way to suppress the fluctuations of the C * V values is to consider bigger boxes with a bigger number of particles, but then the simulations become quite lengthy in time. Shift and shrinking of the C * V (T * tmd ) curves with an increasing ion concentration are well pronounced. However, the behavior of C * V (T * tmd ) around the extremum point of the TMD line must be studied in more detail. After rather comprehensively examining the properties of the TMD lines of the model NaCl salt in the soft-core solvent in question at different ion concentrations, we would like to proceed to the results obtained for other two salts with bigger cations, namely the KCl and RbCl. The assumed diameters of cations were mentioned in the previous section, but we repeat the numbers for the sake of convenience of the reader KCl (σ + = 4.74 Å) and RbCl (σ + = 5.27 Å). In the former case, the diameter of cation does not differ much from the diameter of the anion σ Cl = 4.42 Å. From the analysis of the T * −ρ * projections of the TMD lines given in figure 5 (a), we would like to conclude that in the present formulation of the interactions in the model, the principal factor determining the location of the TMD lines seems to be the excluded volume effects. Namely, at x ion = 0.1, the difference between TMD lines is reasonably big, if one changes from NaCl salt to KCl salt. The difference between KCl and RbCl cases is small as expected, because the cation diameter does not change much. It seems that the shift of the TMD line to higher densities is determined by the amount of excluded volume added to the system. On the other hand, the enhanced excluded volume effects also lead to a lower extremal point along the temperature axis of the TMD lines. However, this shift of TMD lines on temperature axis is to a great extent effected by the electrostatic interactions between ions. If the ion concentration increases for each salt model, the lines substantially shift to lower temperatures. Unfortunately, the solvent in question (due to a set of parameters assumed) does not exhibit vaporliquid criticality. Therefore, we do not have a natural lower boundary for the TMD lines in, for example, the T * − ρ * projection. Our preliminary calculations have shown that quite short TMD lines still exist for x ion = 0.31 and for 0.32 for NaCl model. The shifts of the TMD lines along the density and temperature axes given in figure 5 (a) are also reflected by the behavior of excess pair entropy and of the reduced heat capacity, see figures 5 (b) and 5 (c), respectively. Shrinking of these curves with an increasing ion concentration is quite evident. Trends of the behavior of S * 2,exc and of C * V are similar to those discussed for the model NaCl salt in the soft-core solvent just above.
To conclude, in the present short communication we considered one extension of our very recent work [2] concerning a simple model mixture of ions and soft-core solvent fluid by using Monte Carlo computer simulations. Our extension consists in considering ions of non-equal diameters and different from the diameter of solvent particles. We were interested to study how the presence of ion solutes affects thermodynamic properties of the system compared to the single-component fluid, specifically the temperature of maximum density line. For this purpose, we studied three model solute salts, namely the model for NaCl, KCl and RbCl. Our results confirm that the mixtures based on a soft-core solvent in question preserve this anomalous behavior at x ion = 0.1 0.2, 0.25, and 0.3, if the cation is small, e.g., for NaCl model. However, for the model with a big cation, e.g., RbCl the density anomaly disappears at x ion around 0.27.
The TMD lines were constructed from the simulation data for a set of isochores. They shift to lower temperatures and higher densities with an increasing ion concentration in the T * − ρ * plane. The TMD line shrinks and moves to lower values of the reduced temperature with an increasing ion concentration for each model salt. These observations agree with the results of molecular dynamics simulations [32] for much more elaborate model for aqueous solutions. We have seemingly made interesting observations about the dependence of the excess pair entropy and constant volume heat capacity along the TMD line at different ion concentrations. Still it is necessary to complement our study by simulations of uncharged solutes with different characteristics mixed with applied or another soft-core model for the solvent subsystem. It is known that in the presence of uncharged, non-polar solutes, the curves describing the anomalies shift as well, see e.g. [18]. However, for example, the TMD line changes essentially depend on the solute concentration and strength of solute-solvent interactions. In order to compare our present observations with the models of uncharged solutes in soft-core solvents, additional simulations are necessary. This project is in progress in our laboratory at present. Intuitively, one should expect the most pronounced differences between the effects of ions and non-polar solutes if the ion concentration is low or, in other words, when the screening of ion-ion interaction is not well pronounced. There remains much room for extensions of the present work. In order to validate the results concerning the trends for anomalies with increasing ion concentrations, it is worth to extend the calculations and establish the solubility limits. It is necessary to explore how the modification of the cross solute-solvent interaction potentials would affect the TMD lines. Another important issue is to investigate various properties of ionic solutions and their peculiarities by using a solvent model characterized by the presence of liquidvapor phase transition and possibly liquid-liquid separation, e.g., [42]. Theoretical developments are in fact desirable in many of these aspects.