Surface tension at the liquid-vapor interface of screened ionic mixtures

M.González-Melchor 1 , C.Tapia-Medina 2 , L.Mier-y-Terán 1 , J.Alejandre 3 ∗ 1 Depto. de Fı́sica, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186, Col. Vicentina, 09340 México D.F., México 2 Depto. de Energı́a, Universidad Autónoma Metropolitana-Azcapotzalco, Av. San Pablo 180, Col. Reynosa, 02200 México D.F., México 3 Depto. de Quı́mica. Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186, Col. Vicentina, 09340 México D.F., México


Introduction
Phase behavior and interfacial properties are relevant in many technological applications.Simulation results concerning interfacial properties of complex and simple fluids have been performed in the last decade, particularly, at the liquid-liquid and liquid-vapor interfaces.These studies include the surface tension of simple fluids, molten salts, hydrocarbons, water, bilayers and monolayers of surfactant molecules among other systems [1][2][3][4][5][6].
The coexistence and interfacial properties can be obtained from the direct simulation of interfaces using molecular dynamics (MD) and Monte Carlo (MC) techniques.These tools give us insight about how molecular interactions affect interfacial properties.
In this work we have undertaken the study of interfacial properties of systems with screened electrostatic interactions.This physical situation is often found in complex systems, where different length and time scales are involved and charged particles are present.Ionic fluids are systems where the properties are strongly determined by electrostatic interactions.Examples of these systems are colloidal suspensions of charged particles, electrolyte solutions and molten salts.The Restricted Primitive Model (RPM) has been the starting point to describe properties of charged particles using theory and simulation methods [7,8].The interaction between particles in RPM is purely electrostatic and its range cannot be modified.However, there are systems like colloidal suspensions of charged particles, where the interaction range can be changed by the counterions in solution and by addition of salt [9,10].In this case the relevant forces between particles can be well described as screened electrostatic interactions and a Yukawa-type potential model is suitable.This model contains repulsive and attractive interactions and has been used to study simple fluids and charged particles in solution.Figure 1 (a) shows an all components ionic solution while a simplified model of an ionic mixture is shown in figure 1 (b).
Colloidal suspensions of charged particles have been extensively studied and a "liquid-vapor" phase transition has been observed experimentally [11].The phase behavior exhibited by mixtures of oppositely charged colloidal particles in suspension can be described, at a basic level, by a model containing screened attractive/repulsive interactions [9,10].The simplest model consists of charged hard sphere particles interacting with the screened Coulombic potential and it is called the Yukawa restricted primitive model (YRPM).
The Yukawa model with a hard-core repulsive interaction has been very useful in obtaining thermodynamic properties of one component and mixtures of charged particles using integral equations theory [7,8,12,13].The liquid-vapor equilibrium of YRPM has been studied by the mean spherical approximation (MSA) equation [8,13].One of the main conclusions obtained in those works is that the critical temperature (T c ) increases and has a maximum as the potential decreases its range of interaction.This behavior is opposite to that found in simple fluids for the attractive hard-core Yukawa model (AHCYM) where T c decreases as the potential is short ranged [14,15].As the range of interaction decreases (large λ) the binary ionic fluid should behave as a binary mixture of hard sphere particles without stable liquid-vapor equilibrium.We have recently confirmed in binary ionic fluids by using simulation techniques that T c increases by decreasing the range of interaction [16].The surface tension of the attractive HCY has been obtained previously by using molecular dynamics (MD) and Monte Carlo techniques [15].To the best of our knowledge, there are not calculations of interfacial properties for the YRPM.In this work we used two types of MD simulations to obtain the surface tension of binary ionic fluids.First we applied a hybrid molecular dynamics (HMD), where the hard sphere and continuous forces are evaluated during simulations.In the second case, the hard-sphere interaction is replaced by a soft repulsive potential and standard continuous MD is used.
We have shown that when the hard sphere interaction in the RPM is replaced by this soft repulsive potential, the coexistence densities and interfacial properties of this soft primitive model are essentially the same as those predicted by the true RPM [17,18].On the other hand, MD simulations of continuous potentials are faster than those performed with HMD.This is particularly important in binary ionic fluids where fast frequency collisions between hard spheres is found due to molecular association.In addition, to check the equivalence between the soft and hard core models we also studied one component systems using the AHCYM.The interfacial properties for the AHCYM have been reported before [15] and we used them to validate the results from this work.
The rest of this work is organized as follows: in section 2 the potential models employed are presented.Simulation details and definition of properties are given in section 3. Section 4 contains the interfacial properties of one component and binary ionic fluids interacting with screened Coulomb potential.Finally, conclusions are presented.

Potential model
The model studied in this work is a binary mixture of charged particles with equal size σ.Half of them are negatively charged and the other half are positively charged.The interaction between two particles is given by the soft Yukawa model (SYM), where r is the distance between particles, λ is the inverse of the reduced screening length, q i = ±e is the charge of particle i and e is the electron charge.The variables 0 and r are the vacuum permittivity and the dielectric constant of the medium, respectively.The potential u ij (r) is zero for distances greater than the cutoff distance, R c .The parameters A, n and σ s define the short-range repulsive interaction.
A second model was also used to calculate the surface tension, the YRPM, in this case the first term in equation ( 1) was replaced by a hard sphere potential.The advantage of using the SYM is that it is a continuous function and standard MD can be applied while in the YRPM model a HMD has to be used, where the hard sphere and continuous forces are calculated during the dynamics.To emulate the hard sphere interaction we choose A/ε = 20000, σ s = 0.95σ and n = 225.Figure 2 shows the potential of the SYM for ionic fluids with λ = 6.The one component system has only the attractive part.Hereafter we use AHCYM and attractive soft Yukawa model (ASYM) to denote the one component systems with the hard sphere and soft interactions, respectively.The calculations were performed in reduced units.The distance is reduced with σ and the energy with ε = e 2 /(4π 0 r σ).So, the reduced temperature is where k B is the Boltzmann's constant and T is the absolute temperature.

Simulation details and definition of properties
The simulations were performed in a rectangular box elongated along the z axis having dimensions L x = L y , L z > L x and volume V = L x L y L z .The box length in z direction was at least 5 times longer than in x and y directions.A dense liquid slab of N particles surrounded by vacuum was placed at the center of the simulation cell.Full periodic boundary conditions were used in all directions.In order to minimize surface area effects on the surface tension introduced by periodic boundary conditions we considered N = 4116 particles and L x = 14σ.These values are large enough to avoid finite size effects in the calculation of surface tension.The cutoff distance was set to R c = 3.5 σ.There was used a multiple time step integration scheme [19] to improve the efficiency of simulation programs.Three different reduced time steps, t * = t( /mσ 2 ) 1/2 , were used in HMD and in continuous MD.The smaller one was t * = 0.001 and the largest was t * = 0.012.The temperature was controlled using a Nosé-Hoover thermostat.At least 50000 large time steps were performed to equilibrate the system.The average properties were evaluated for additional 150000 large time steps, i.e, 1.5 × 10 6 small time steps.
The surface tension of one planar interface was computed through the pressure tensor route, where P αβ (α, β = x, y, z) are the pressure tensor components.The component P xx for systems that contain the hard sphere interaction, AHCYM and YRPM potentials, is defined as, where all particles have the same mass m, v x i is the velocity of particle i in the x-direction, x ij = x i − x j is the relative position in the x-direction of atoms i and j and f x ij is the x-force component between particles.In equation 3 the first term is the kinetic contribution, the second is the continuous contribution of intermolecular interactions, the third takes into account the hard sphere interaction and the fourth evaluates the pressure at the cutoff distance.In the third term t sim is the total simulation time, ∆v x i is the change of velocity that a particle has got before and after a collision and N c is the number of collisions between hard-spheres.In the fourth term, δ(s) is a δ-function and was calculated as in [15].To simulate the ASYM or SYM the third term in equation 3 is changed by the corresponding continuous part of the soft potential.The reduced surface tension is given by γ * = γσ 2 / .
The density profile, ρ(z) was calculated as where ∆N is the average number of particles with coordinates between z and z + ∆z, and ∆V = L x L y ∆z is the volume of the corresponding slab.A value of ∆z = 0.02σ was used for all simulations.
At the end of simulations, the average density profiles were fitted to a hyperbolic tangent function, where ρ l and ρ v are the liquid and vapor densities, respectively, z 0 is the Gibbs dividing surface and δ is a measure of the width of the interface.
Figure 3 shows a typical density profile where the so-called 10 − 90 width of the interface δ is shown.The liquid and vapor densities and width of the interface were obtained using the fitted density profiles.The coexistence curves were build up as a function of temperature and range of interaction.The critical parameters were calculated using a rectilinear diameter law with a critical exponent β = 0.32.

Results
Although the coexistence densities and surface tension of the AHCYM have previously been reported [15] for one component systems, we re-evaluate them in the first subsection of this work but using the ASYM to see how well the results from this potential are compared with those of the AHCYM.The analysis of coexistence densities and critical properties of screened binary ionic mixtures defined by the YRPM is the subject of a future work [16].At the end of this section we report, for the first time, the surface tension of screened binary ionic mixtures using the SYM and YRPM potentials.

A. One component attractive soft Yukawa fluid
The coexistence curves of the one component ASYM for different ranges of interaction are shown in figure 4 and given in table 1.The main effect of reducing the range of interaction is to decrease the critical temperature and to increase the difference between the liquid and vapor densities.The ASYM results are compared with those from the AHCYM for λ = 1.8, 3 and 4 taken from [15] and shown with stars.The agreement between both models is reasonably good.We performed a systematic study of interfacial properties in the region where the AHCYM has a stable equilibrium.It has been found for the AHCYM that the liquid-vapor equilibrium is not stable for λ > 6 [14].After this value the system has a fluid-solid equilibrium.
It was not our interest to search the stability of the ASYM but to analyze how well the interfacial properties of the AHCYM are reproduced.The critical temperature and density for the ASYM are given in table 2 and compared with those of AHCYM obtained by HMD simulation of interfaces.The agreement in critical density between both models is excellent and the critical temperature of the ASYM is 3% higher than that of the AHCYM with λ = 3.The critical temperature decreases as the potential is short-ranged because the interaction between molecules is less attractive.
The surface tension of one component fluids interacting with the ASYM are shown in figure 5 for different ranges of interaction parameters and also shown in table 1.The surface tension decreases with temperature and with λ.The agreement between the ASYM and the AHCYM is excellent for λ = 3 and 4. The surface tension seems to be less sensitive to the potential details than the coexistence densities, as can be seen for results using λ = 3.The surface tension results for the ASYM were fitted to γ = γ 0 (1 − T /T c ) µ (where γ 0 and µ are constants) according to the mean field theory [20] with µ = 1.26 as an alternative to obtain the critical temperature.The critical temperature results obtained by this fitting procedure are also given in table 2 and show a good agreement with those obtained directly from the rectilinear diameter law.0.306 AHCYM [15]  2.0 1.028 1.030 0.329 ASYM 2.5 0.868 0.881 0.340 ASYM 3.0 0.725 0.351 AHCYM [15]  3.0 0.748 0.741 0.347 ASYM 4.0 0.593 0.361 AHCYM [15]  4.0 0.599 0.599 0.363 ASYM 5.0 0.526 0.525 0.378 ASYM

B. Binary screened ionic mixtures
Simulation results show that the critical temperature increases when the interaction between molecules is short-ranged [16], in agreement with predictions of the mean spherical approximation [8].In this work, the liquid-vapor interface of a binary ionic mixture is analyzed in terms of the surface tension, shown in figure 6.As expected the surface tension decreases with temperature.An excellent agreement is observed for results between SYM and YRPM potentials using λ = 2, i.e, the SYM can be safely used to study interfacial properties of screened ionic fluids.As the interaction between particles becomes short-ranged the surface tension increases, in contrast with the AHCYM of simple fluids where the surface tension increases with the range of interaction.It is interesting that, although the RPM (λ = 0) is a long ranged potential, it has a lower surface tension than those of YRPM with screening parameters λ > 0. This is because ionic fluids interacting with a screened poten-tial form clusters of different size due to molecular association [16].The region of temperature for a given interaction range, where the liquid-vapor interface is stable, decreases with λ.This suggests that at some λ the liquid-vapor equilibrium will not be stable and the system would shift into a fluid-solid transition.Table 3 resumes surface tension of the screened binary ionic systems using the SYM as a function of λ.It is not easy to compare the width of the interface for different ranges of interaction because the region of temperature where the interface is stable is not the same.However we can observe that for λ = 6 and 10 at T * = 0.16 (see table 3) the width of the interface is smaller for λ = 10 because these systems have larger surface tensions.

Figure 1 .
Figure 1.Molecular models shown schematically.(a) All the components are considered: solvent, colloidal particles and ions.(b) Simplified ionic model with screened electrostatic interactions in a continuous medium with dielectric constant r .

Figure 2 .
Figure 2. The SYM potential for a binary ionic fluid with screening parameter λ = 6.0.

Figure 3 .
Figure 3. Hyperbolic tangent function fitted to a density profile obtained from MD simulations.The width of the interface δ, is shown schematically.

Figure 4 .Figure 5 .
Figure 4. Liquid-vapor coexistence curve for the ASYM as a function of λ.Stars are AHCYM results from [15] for λ = 1.8, 3 and 4, from top to bottom.The critical properties are also shown.

Figure 6 .
Figure 6.Reduced surface tension of binary mixtures using the SYM as function of the reduced temperature T * and different ranges of interaction.The results from this work for the YRPM are shown with stars for λ = 2.The continuous line without symbols are HMD results for the RPM taken from [17].

Table 1 .
Coexistence densities and surface tension of the ASYM as a function of the screening parameter.The number 0.65 5 means 0.65 ± 0.05.

Table 2 .
Critical parameters for the ASYM for different ranges of interaction.The critical temperature obtained by fitting the surface tension from MD simulations to a function predicted by mean field theory, T * MF

Table 3 .
Surface tension as function of temperature for binary ionic mixtures using the SYM for different ranges of interaction.The results for the YRPM with λ = 2 are also given.