# Interfacial tension behavior of binary and ternary mixtures of partially miscible Lennard-Jones fluids: A molecular dynamics simulation.

## I Introduction

The investigation of interfacial tension behavior in complex fluids is of relevance from the theoretical as well as from the practical point of view . In particular, the nature of the liquid-vapor interface has been extensively studied both analytically and numerically, using Monte Carlo and Molecular Dynamics  algorithms. On the other hand, the properties of the liquid-liquid interface, which are of importance in different biological systems as well as in various technological applications , have received significantly less attention. One reason may be the complexity of the topology of the phase diagram for liquid mixtures. It is known that the intermolecular properties of a single component fluid have no important effect on the topology of its thermodynamic phase diagram leading to the concept of corresponding states. By contrast, for liquid mixtures the interplay between the molecular properties of the components as well as the type of the intermolecular potential between different species leads to a highly non-trivial topology of the thermodynamic phase diagram. Recently, the properties of liquid-liquid planar interfaces have begun to be studied by means of molecular dynamics simulations using Lennard-Jones (LJ) interactions, by Monte Carlo simulations of benzene-water interface and hexanol-water  as well as by density functional theory. In addition, there have also been molecular dynamics simulations for the interface of water-alcohols using more complicated intermolecular potentials that include the relevant molecular degrees of freedom. . The density and pressure profiles of a symmetrical binary mixture of LJ fluids at the reduced temperature $`T^{}=0.827`$ was studied in a region of the phase diagram where the system is separated into two phases. It was found that as the system becomes more miscible there is more diffusion of particles from one phase to the other and, as a consequence, the interfacial tension decreases. Some dynamic properties of this system were investigated at two reduced densities, $`\rho ^{}=`$0.8 and 1.37 for several temperatures from which a crude sketch of the phase diagram was suggested. A third type of particle representing a simple model of amphiphilic molecule was introduced finding that for low concentrations (of the order of 6%) there is no important effect on the binary fluid interface. In recent papers the structural properties of the above mentioned binary mixture was investigated at low temperatures and high pressures finding a stable oscillatory behavior in the density profiles close to the interface. These oscillations decrease as the temperature increases. It is also suggested that as pressure grows the interfacial tension increases. A recent density functional investigation of a binary fluid mixture shows that the interfacial tension as a function of temperature can have a maximum.

With the purpose of getting a more complete understanding of the thermodynamic and structural properties of this liquid-liquid interface we have performed extensive equilibrium molecular dynamics simulations. The main results found in this work are: (i) the non-monotonic behavior of the interfacial tension as a function of temperature and (ii) the monotonic decay of this quantity as the concentration of a third type of particle, amphiphilic-like, increases. The layout of this paper is as follows: In section II we introduce the model and details of the simulations. In section III we present the results for the thermodynamic as well as structural properties of the binary and ternary systems. Finally, in section IV the conclusions of this investigation are presented.

## II Model and simulations.

We consider a symmetrical binary mixture of partially miscible Lennard-Jones fluids. The intermolecular interactions $`F_{_{AA}},F_{_{BB}}`$ and $`F_{_{AB}}`$ between particles A-A, B-B and A-B are described by modified LJ potentials that yield the forces,

$$F_{_{XY}}(r)=\{\begin{array}{cc}\frac{24}{r}ϵ\left[2\left(\frac{\sigma }{r}\right)^{12}\alpha _{_{XY}}\left(\frac{\sigma }{r}\right)^6\right]\hfill & \text{ if }rR_c\hfill \\ 0\hfill & \text{ if }r>R_c\hfill \end{array}$$
(1)   
where $`ϵ`$ and $`\sigma `$ are the same for all interactions and the parameter $`\alpha _{_{XY}}`$ controls the miscibility of the two fluids. For a partially miscible binary fluid we choose $`\alpha _{_{AA}}=\alpha _{_{BB}}=1`$ and $`0\alpha _{_{AB}}<1`$. The interaction between particles of the same type is energetically more favourable, stronger binding, than the interaction A-B between different particles. For this reason one would expect a separation of the species at lower temperture when the entropy looses versus potential energy. With the aim of understanding the process by which a “surfactant” weakens the interfacial tension and leads to the mixed state, a third species C is introduced to try to emulate a surfactant-like particle. In this case we consider the same parameters in Eq. 1 and $`\alpha _{_{CC}}=\alpha _{_{AC}}=\alpha _{_{BC}}=1`$.

We have carried out extensive equilibrium molecular dynamics (MD) simulations using the (N,V,T) ensemble for the particular miscibility parameter value $`\alpha _{_{AB}}=0.5`$. The range of the inter-molecular potential was set equal to three times the particle diameter $`\sigma `$ unless otherwise stated. In most of the simulations a total of N=$`1728`$ particles were used. Nonetheless, to check for finite size effects we also carried out simulations with N=$`2592`$ particles. These particles are placed in a paralelepiped of volume $`L_x\times L_y\times L_z`$, with $`L_x=L_y`$ and $`L_z=2L_x`$, applying periodic boundary conditions in the $`x,y`$ and $`z`$ directions. The particles were initially placed on the sites of an FCC lattice forming a perfect planar interface, that is, all particles of type A are on the left side of the box while those of type B are on the opposite side. In this way one can obtain a minimum of two interfaces due to periodic boundary conditions. If one starts with a statistical mixture usually more than two demixed regions develop giving rise to more than two interfaces.

The three component system initial configuration was chosen from an equilibrated separated binary system picking at random $`N_c/2`$ particles from type A and $`N_c/2`$ particles from type B and replacing them by particles of type C. This way of putting the third species in the system makes the total density to remain constant. It is customary to carry out the simulations using the following reduced units for the distance $`r^{}=\frac{r}{\sigma }`$, particle linear momentum $`p^{}=\frac{p}{\sqrt{mϵ}}`$ and time $`t^{}=\frac{t}{\sigma }\sqrt{ϵ/m}`$. In these definitions $`m`$ is the mass of each particle, which is taken to be the same for all particles, $`\sigma `$ is the particle diameter and $`ϵ`$ is the depth of the LJ potential. Similarly, one can define reduced thermodynamic quantities as follows: $`T^{}=\frac{k__BT}{ϵ}`$ that represents the reduced temperature with $`k__B`$=Boltzmann’s constant and $`\rho ^{}=\rho \sigma ^3`$ for the reduced density, with $`\rho =N/V`$. The equations of motion were integrated using a fourth order predictor-corrector algorithm with an integration step-size of $`\mathrm{\Delta }t^{}0.005`$, which in standard units is of the order of $`10^5`$ nanoseconds in the scale of argon. The particles initial velocities were assigned from a Boltzmann distribution. The equilibration times for most of the simulations were of the order of $`10^4`$ time steps. Thermodynamic quantities were measured every 50 time-step iterations up to a total of 5$`\times 10^5`$ to $`10^6`$ measurements from which averages were evaluated. This amounts to a simulation time between 5-10 nanoseconds in the scale of argon. At the start the reduced homogeneous density in the simulation cell was set equal to $`\rho ^{}=0.844`$, which is close to density of the triple point of argon. Due to the inhomogeneity which develops in the system around the interface, the densities of the bulk phases later are slightly higher than this starting density. The bulk densities are evaluated when the system has been equilibrated and they depend on the conditions of the system and on $`T^{}`$.

## III Result and discussion

### A Binary mixture

Since our interest is to study the structural and thermodynamic properties of the interface, we have carried out a set of simulations for a sequence of temperatures below the critical demixing temperature for two system sizes, namely, $`N=1728`$ and $`N=2592`$ particles. We have considered these values of $`N`$ to find out about possible finite size effects. All the quantities studied show qualitatively the same tendency for the two system sizes. The interfacial behavior is investigated by calculating the density profiles, the pressure, and the interfacial tension at different temperatures. The relevant parameters of the investigated systems are summarized in Tables I and II where the bulk-values of the density and the pressure are specified. To emphasize the separated nature of the system, the values of the reduced total density in the bulk A-rich phase $`\rho ^{\mathrm{bulk}\mathrm{A}}`$, the reduced density of particles A, $`\rho __A^{\mathrm{bulk}\mathrm{A}}`$, and particles B, $`\rho __B^{\mathrm{bulk}\mathrm{A}}`$, are given. Due to the symmetry of the interactions, the B-rich phase is symmetric to the A-rich phase.

At sufficiently low temperatures the reduced density of B particles in the phase of A is very small indicating that the system is fully separated. As temperature increases the total density of the bulk region decreases slightly because the diffusion of particles through the interface increases, making the inhomogeneity smaller and driving the system towards a mixed state. This can be seen in Figures 1 and 2 where the density profiles along the $`z`$ direction (longer side of the box) are shown for the $`N=1728`$ and 2592 systems, respectively. The continuous line corresponds to a temperature of $`T^{}=0.827`$, which can be considered low, (systems 2 in Tables I and II) while the dashed line corresponds to a higher temperature $`T^{}=1.4`$ (systems 4 in Tables I and II).

When looking at the density profiles for the two system sizes studied, Figures 1 and 2, we observe that at $`T^{}=0.827`$ these quantities exhibit an oscillatory structure at the interfaces. This reminds us about the behavior of the density profile in front of a hard wall, as discussed in reference . This kind of structure is also found in the liquid-vapor system, but the oscillations are significantly stronger in the liquid-liquid interface. As expected, one also sees from these figures that the oscillations in the bulk region are less pronounced for the larger system.

The values of the reduced normal pressure $`P_n^{}`$ shown in column sixth of Tables I and II follow an almost linear behavior as a function of temperature in the range studied. The normal and transversal pressures profiles $`P_n(z)`$, $`P_t(z)`$ were calculated using the definition of the Irving-Kirkwood pressure tensor, which for a planar interface is given by 

$`P_n(z)`$ $`=`$ $`\rho (z)k__BT`$ (2)   
$``$ $`{\displaystyle \frac{1}{2A}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{z_{_{ij}}^2u^{}(r_{_{ij}})}{r_{_{ij}}|z_{_{ij}}|}}\theta \left({\displaystyle \frac{zz__i}{z_{_{ij}}}}\right)\theta \left({\displaystyle \frac{z__jz}{z_{_{ij}}}}\right),`$ (3)   
$`P_t(z)`$ $`=`$ $`\rho (z)k__BT`$ (4)   
$``$ $`{\displaystyle \frac{1}{4A}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{[x_{_{ij}}^2+y_{_{ij}}^2]u^{}(r_{_{ij}})}{r_{_{ij}}|z_{_{ij}}|}}\theta \left({\displaystyle \frac{zz__i}{z_{_{ij}}}}\right)\theta \left({\displaystyle \frac{z__jz}{z_{_{ij}}}}\right).`$ (5)   
In these expressions the first term corresponds to the ideal gas contribution while the second comes from the intermolecular forces. For a dense system and long ranged intermolecular potentials the latter term yields the larger contribution. On the left side of Figure 3 we show the normal and tangential components of the pressure tensor as functions of $`z`$ for $`T^{}=0.827`$ and $`T^{}=1.4`$. The first feature to notice is that $`P_n^{}(z)`$ remains constant for both temperatures as one would expect for a system in thermodynamic equilibrium. The pressures $`P_n^{}`$ range from 0.25 kbar to 5.18 kbar in the scale of argon for the systems of Tables I and II. On the right side of the same figure the ideal as well as the configurational parts of $`P_n^{}(z)`$ and $`P_t^{}(z)`$ are plotted. At the lower temperature $`T^{}=0.827`$ the pressure profiles clearly show the oscillations related to the oscillatory behavior of the density profiles. One can learn that the configurational contribution is much larger than the ideal gas part since we are dealing with a dense system. These profiles are needed to evaluate the interfacial tension by means of the mechanical definition.

$$\gamma =_{\mathrm{bulk}_\mathrm{A}}^{\mathrm{bulk}_\mathrm{B}}[P_n(z)P_t(z)]𝑑z.$$
(6)   
Also, a straightforward evaluation of $`\gamma `$ can be done using the Kirkwood-Buff formula

$$\gamma =\frac{1}{4A}\underset{i<j}{}\left(1\frac{3z_{_{ij}}^2}{r_{_{ij}}^2}\right)r_{_{ij}}u^{}(r_{_{ij}}).$$
(7)   
In this latter equation there appears an additional factor of $`\frac{1}{2}`$ which comes from having two interfaces in the system. From the calculational point of view it is better to use the Kirkwood-Buff formula rather than the mechanical definition since fluctuations in $`P_n(z)`$ and $`P_t(z)`$ may introduce important inaccuracies in the evaluation of $`\gamma `$ according to Eq. 6. As a matter of fact, we calculated the interfacial tension using both expressions and found consistency in the values within the statistics of the simulations. The range of temperatures, in reduced units, in which $`\gamma `$ has been studied is $`0.6T^{}3.0`$, that in the scale of argon corresponds to $`72KT359.5K`$.

Unlike for the liquid-vapor interface we find a non-monotonic behavior of $`\gamma (T)`$ for the liquid-liquid interface. In Figure 4 we show the behavior of $`\gamma ^{}`$ as a function of $`T^{}`$ for systems with $`N=1728`$ and $`N=2592`$ particles. The main feature of this graph is the maximum of the interfacial tension at about $`T^{}1.1`$. Such a maximum has been reported recently in the context of density functional theory. It can be understood physically as follows: A stronger mixing near the interface introduces weaker A-B bonds and raises the potential energy. This leads to an increase of $`\gamma `$, the free energy of the interface, until at high temperatures the entropy contribution drags it down. We have calculated $`\gamma (T)`$ for these two systems to chek for possible finite size effects. In doing so we have to assure that the bulk thermodynamic quantities are approximately the same for both systems. In particular, we monitored very closely the bulk density since small variations in this quantity lead to significant changes in the pressure of the system at low temperatures. So, the important key is to have approximately the same pressure for both systems. This can be achieved by adjusting the side $`L_z`$ of the simulational box for the larger system, mantaining its cross section constant. In this way we eliminate possible variations of $`\gamma (T)`$ due to variations in the cross section area. We found that the optimal average value of $`L_z^{}`$ was 30 for the temperature range studied. In Figure 4 and tables I and II one also sees that at low temperatures, the values of $`\gamma ^{}(T^{})`$ (between the two system sizes) are at the most 4% different, while at higher temperatures ($`T^{}>1.1`$) the values are much closer to each other. One also observes that small variations in the bulk density lead to important changes in the pressure at low temperatures. However, at higher temperatures the variations in pressure are smaller. This is the expected behavior for a single LJ fluid. Therefore, differences in $`\gamma (T)`$ are due to the variations in pressure, although the former are much smaller than the latter. This happens because the interfacial tension is obtained from the difference between the normal and tangential components of the pressure tensor. This behavior explains why the differences and consistencies in the values of $`\gamma ^{}(T^{})`$.

### B Ternary mixture

In this section we investigate the role that a third species plays in the interface properties of a demixed binary system. A third species — C-particles — is introduced in the system and the behavior of the interfacial tension is studied as a function of concentration. The interactions of the C-particles $`F_{_{CC}}=F_{_{CA}}=F_{_{CB}}`$ are all assumed equal to the strong interactions $`F_{_{AA}}=F_{_{BB}}=F_{_{CC}}`$. When placed between the A and B, C particles avoid the weak A-B bonds and lower the potential energy. It is found that as $`N__C`$ increases $`\gamma `$ decays monotonically. This is reminiscent of surfactant-like behavior in some ternary systems. In Figure 5 we plot the density profiles of the ternary system as well as the reduced total density for different concentrations of C particles at $`T^{}=1.4`$. These profiles yield a clear evidence that C-particles like to be at the interface position rather than in the bulk phases. In this way the C-particles diminish (screen) the energetically unfavourable interactions between particles A-B. In Table III a summary of the reduced bulk densities of particles A, B and C in the A-rich phase for $`T^{}=0.827`$ and $`T^{}=1.4`$ is given. As $`N__C`$ increases the density of B in the A-rich phase also increases. This means that the C particles help the B particles to diffuse into the A-phase because the B particles can be solvated by C and therefore avoid the weak A-B interactions. As $`N__C`$ becomes sufficiently large this mechanism drives the system to a mixed state. Therefore, the C particles act like a surfactant or emulgator. This mechanism may have potential technological implications when trying to design compatibilizers. In Figure 6 the reduced interfacial tension is plotted as a function of $`N__C`$ for two different temperatures. With increasing bulk concentrations of C-particles we see a monotonic decay. This low concentration behavior is consistent with that found in reference where amphiphile like-particles C are introduced. In such a model $`\gamma ^{}`$ as a function of $`N__C`$ shows a linear decay in the small amphiphile concentration range.

In this same figure we also observe that the curve $`\gamma ^{}(N__C)`$ at $`T^{}=1.4`$ is slightly above the curve of $`\gamma ^{}(N__C)`$ at $`T^{}=0.827`$. This is an unexpected result since usually $`\gamma ^{}`$ decreases when $`T^{}`$ increases. Nonetheless, we should recall that these two points lie in the region where $`\gamma ^{}`$ has its maximum as a function of $`T^{}`$ for the binary system. This fact might also be responsible for the shoulder in the curve for $`\gamma ^{}`$ at $`T^{}=1.4`$ in the region $`200<N__C<300`$. We would like to emphasize that this structure is outside of the statistics of the simulations since the error bars are smaller than the difference between the points. We have also evaluated the excess $`\mathrm{\Gamma }__C`$ of C-particles given by Equation 8 and plot this quantity as a function of $`N__C`$ in the inset of Figure 6.

$$\mathrm{\Gamma }__C=_{\mathrm{bulk}_\mathrm{A}}^{\mathrm{bulk}_\mathrm{B}}\left(\rho __C(z)\rho __C^^B\right)𝑑z$$
(8)   
We find that for all values of $`N__C`$, $`\mathrm{\Gamma }__C(T^{}=0.827)>\mathrm{\Gamma }__C(T^{}=1.4)`$ which means that there are more particles of type C located at the interface at the lower temperature and therefore they screen more effectively the energetically unfavourable interactions between A and B particles, diminishing even more the interfacial tension. In the range $`200N__C400`$ the excess of C particles at the interface at $`T^{}=1.4`$ increases slower than at $`T^{}=0.827`$ giving rise to a shoulder in $`\gamma ^{}`$ as shown in figure 6.

Since in the above results not the pressure but the overall density $`\rho ^{}=0.844`$ was kept constant, in Figure 7 we plot the behavior of $`P_n^{}`$ as a function of $`N__C`$ for the two temperatures indicated above. One notes that the hydrostatic pressure decreases slowly as $`N__C`$ increases. This is due to the fact that $`N__C/2`$ particles from fluid A and the same amount from fluid B were switched to particles C, replacing the weak A-B attraction by the stronger C-A and C-B attraction. This increases the cohesion and reduces the pressure. One could be tempted to obtain the upper curve from the lower by simply adding the ideal gas term $`\rho ^{}\mathrm{\Delta }T^{}`$, however, this is not the case because, as noted above with regard to Figure 3, the contribution of the configurational part to the pressure is the more relevant one.

## IV Conclusions

Molecular dynamics simulations for a mixture of two kinds of molecules A and B lead to a phase separation with an A-rich and B-rich phase and an interface between the phases. This happens if the attractive potential between A-B is weaker than that between A-A and B-B and if density and temperature are in the two phase region of demixing. The interface between the demixed phases has also been investigated. We have calculated density profiles at different temperatures, the components of the pressure tensor across the interface and have evaluated the interfacial tension $`\gamma `$ for a series of temperatures. The important result of the first part of this investigation is the maximum in $`\gamma (T)`$ (Figure 4). We explain this result as a consequence of a strong mixing near the interface that yield weaker A-B bonds thus raising the potential energy, that in turn leads to an increase in $`\gamma (T)`$, until at higher temperatures the interfacial free energy decreases due to an increase of the entropy term.

Then we have added a third kind of particle, “surfactant-like”, and have studied the interfacial tension as a function of its concentration. We found that these particles reduce the interfacial tension $`\gamma `$ as its concentration increases. This happens because these particles screen the unfavorable A-B interactions and lead to a mixing state of the system.

## Acknowledgments

We thank S. Iatsevitch and M. Swiderek for helpful discussions. F.F. would like to thank the Foundation Sandoval Vallarta for financial support and the colleagues at the Physics Department of UAM-I for the grand hospitality during his visit. The calculations were carried out on the Power-Challenge computer at the UAM-I and at DGSCA-UNAM. This work is supported by CONACYT research grants Nos. L0080-E and 25298-E.