Role of the Counterions in the Surface Tension of Aqueous Surfactant Solutions. A Computer Simulation Study of Alkali Dodecyl Sulfate Systems

: We have investigated the surface tension contributions of the counterions, surfactant headgroups and tails, and water molecules in aqueous alkali dodecyl sulfate (DS) solutions close to the saturated surface concentration by analyzing the lateral pressure proﬁle contribution of these components using molecular dynamics simulations. For this purpose, we have used the combination of two popular force ﬁelds, namely KBFF for the counterions and GROMOS96 for the surfactant, which are both parameterized for the SPC / E water model. Except for the system containing Na + counterions, the surface tension of the surfactant solutions has turned out to be larger rather than smaller than that of neat water, showing a severe shortcoming of the combination of the two force ﬁelds. We have traced back this failure of the potential model combination to the unphysically strong attraction of the KBFF counterions, except for Na + , to the anionic head of the surfactants. Despite this failure of the model, we have observed a clear relation between the soft / hard character (in the sense of the Hofmeister series) and the surface tension contribution of the counterions, which, given the above limitations of the model, can only be regarded as an indicative result. We emphasize that the obtained results, although in a twisted way, clearly stress the crucial role the counterions of ionic surfactants play in determining the surface tension of the aqueous surfactant solutions.


Introduction
Surface tension, the intensive counterpart of the surface area of a given phase, is a physical quantity of key importance both in pure interfacial science [1] and in its numerous industrial applications. Surface tension is related to the force that contracts the surface area of condensed phases and originates from the fact that, due to the lack or loss of interactions with the opposite phase, molecules experience an energetically less favorable environment at the vicinity of the surface than in the bulk phase. The fact that surface tension can thus be related to the properties (i.e., energy loss) of the individual molecules naturally gives rise to the question of how the individual molecules contribute to it [2].
A particularly interesting question in this respect is that of how the surface tension is distributed between the different molecules and moieties at the surface of aqueous surfactant solutions. Clearly, surfactant molecules adsorbed at an aqueous surface reduce the surface tension very efficiently; however, it is not clear how the residual surface tension is contributed by the water molecules as well as by the apolar tails and polar or ionic headgroups of the surfactant molecules and, in the case of ionic surfactants, also by the counterions. Computer simulation methods seem to be very suitable to study this problem, as in such studies an appropriately chosen model of the system of interest is seen at atomistic resolution [3]. On the other hand, the validity of the chosen model always has to be checked by comparing its relevant, measurable properties with existing experimental data. However, while numerous properties of the surface of aqueous solutions of a number of different surfactant molecules have been investigated in detail both by experimental  and computer simulation  methods in the past decades, the question of how different molecules contribute to the surface tension has scarcely been studied both for surfactant solutions [56] and for other systems [2,57].
In a recent study, we addressed this question for various surfactants at several surface densities and obtained the rather surprising result that, in the case of ionic surfactants, the counterions play a key role in this respect [56]. Even more strikingly, counterions were found to play opposing roles in the cases of two surfactants considered. More specifically, considering a saturated aqueous surface, the Na + counterion of the dodecyl sulfate (DS − ) anionic surfactant contributed 760%, while the Cl − counterion of the dodecyl trimethyl ammonium (DTA + ) surfactant cation contributed −240% to the surface tension of the corresponding system. These huge contributions were found to be partly compensated by those of the water molecules and surfactant heads, indicating that in the former and latter systems, water molecules and Cl − counterions, respectively, experience a much more favorable environment at the vicinity of the surface than in the bulk liquid phase [56].
This rather unexpected result can be explained using several different assumptions. One possible explanation is related to the asymmetry of the positively and negatively charged counterions, caused by the geometric asymmetry of the charge distribution of the water molecules. An alternative explanation assumes that the exact role of the counterions in this respect is related to their "soft/hard" character according to the Hofmeister series [58]. Thus, "hard" ions or "kosmotropes", i.e., small ions of high charge density and low polarizability, are very strongly attached to the water molecules of their first hydration shell, and hence they are effectively repelled from the aqueous surface [59][60][61][62][63][64][65]. On the other hand, larger "soft" ions of higher polarizability, also referred to as "chaotropes" are less strongly bound to their hydration shell. As a consequence, losing a part of this shell, occurring when the ion approaches the surface of the aqueous phase, corresponds to a smaller energy cost, and hence such ions are less strongly repelled from or even effectively attracted to the liquid surface [59][60][61][62][63][64][65]. In other words, hard counterions, forced to stay close to the surface by the oppositely charged layer of the surfactant heads, have to pay a large energetic penalty and hence give a huge contribution to the surface tension. On the other hand, this energetic penalty is smaller for softer counterions and might even be overcompensated by the energetic gain coming from the vicinity of the layer of oppositely charged surfactant headgroups. This explanation is certainly supported by the fact that Na + is clearly a harder ion than Cl − [59,60,62]; however, its verification requires a systematic study in which only the soft/hard character of the counterion is varied, while all other properties of the system are kept unchanged.
In this paper, we intend to investigate the validity of the second assumption, namely that the surface tension contribution of the counterions is related to their soft/hard character. For this purpose, here we report a set of computer simulations of the liquid-vapor interface of the aqueous solutions of alkali dodecyl sulfate at saturated surface density, in which the alkali counterion is systematically varied from Li + to Cs + . This way, the soft/hard character of the counterions is systematically changed from the hardest Li + to the softest Cs + , while all other characteristics of the system are left unchanged. Our above explanation predicts that upon going to larger and softer counterions, the counterion contribution to the surface tension decreases and might even turn to negative values. The paper is organized as follows. In Section 2 it is explained how the contribution of the different particles to the total surface tension of the system can be estimated from the lateral pressure profile along the surface normal. Details of the simulations performed are given in Section 3, while in Section 4 the obtained results are presented and discussed in detail. Finally, the main conclusions of this study are summarized in Section 5.

Calculation of the Surface Tension Contribution of the Molecules through the Lateral Pressure
Surface tension (γ) can be calculated as an integral of the imbalance between the lateral and normal pressure components (p L and p N , respectively) along the surface normal axis, denoted here by X, i.e., for a single interface [1], It should be noted that while p L varies along the surface normal axis, the normal component of the pressure has to be constant along X due to the requirement of the mechanical stability of the system, and hence it is also equal to the (scalar) pressure of the two bulk phases. Therefore, the profile of p L is directly related to that of the surface tension, and in cases when the scalar pressure of the system is negligible compared to the values of the lateral pressure component at the vicinity of the surface, such as in the case of aqueous systems being in equilibrium with the saturated vapor at ambient temperature, it differs from the surface tension profile only in its sign.
The elements of the pressure tensor can be calculated as: where the first term accounts for the kinetics, while the second one accounts for the virial contribution to the pressure. In this equation, indices i and j run over the particles, while α and β run over the three spatial directions (or, in the case of calculating lateral pressure components, only over the two spatial directions laying in the surface plane); m, v, and f are the mass and velocity of the corresponding particle and the force acting between the corresponding particles, respectively; V is the volume of the system, δ is the Dirac distribution, C ij represents an open path connecting the two particles i and j, parametrized by the vector s, and the brackets < . . . > denote ensemble averaging. As is evident, the value resulting from Equation (2) depends on the particular choice of the integration path, C ij [66]. Fortunately, it has been shown that the use of several reasonably chosen integration contours, such as the Irving-Kirkwood [67] or the Harasima path [68], result in comparable lateral pressure profiles [69]. Among these paths, the use of the Harasima contour, connecting the interacting particles by lines parallel with the space-fixed axes, has several important advantages. Thus, this way the lateral component of the pressure can also be calculated in cases where the force is not pairwise additive [69]. The importance of this point becomes apparent considering that when the long range part of the electrostatic interaction is calculated by an Ewald summation-based method, its reciprocal space term is not pairwise additive either. The other important advantage is that this way the lateral pressure contribution coming from the interaction of a given particle pair is distributed evenly between the two particles [57]. As a consequence, lateral pressure can be treated as if it was an additive property of the individual atoms, and hence the contribution of an appropriately chosen group of atoms to the lateral pressure and, through Equation (1), the surface tension can be meaningfully discussed.

Computational Details
Molecular dynamics simulations of the liquid-vapor interface of five aqueous alkali dodecyl sulfate solutions, containing Li + , Na + , K + , Rb + , and Cs + counterions, were performed on the canonical (N, V, T) ensemble at the temperature T = 298 K. The X, Y, and Z edges of the basic simulation box were 265, 62.82, and 62.82 Å, respectively, with X being the surface normal axis. The basic box has contained 8000 water molecules, 192 DS − ions, and 192 alkali counterions in every case. Considering that practically all DS − ions were attached to the liquid surface during the entire course of the simulations, this composition corresponds to the surface density of 4 µmol/m 2 , representing saturated surface coverage.
The DS − ions, alkali counterions, and water molecules were described by the GROMOS96 [70,71] and KBFF [72] force fields and by the rigid, three-site SPC/E potential [73], respectively. The choice of the KBFF model of the counterions instead of GROMOS96, used in our previous publication [56], is dictated by the fact that, unlike GROMOS96, KBFF parameters are available for the full set of the alkali cations [72]. CH 2 and CH 3 groups of the DS − ions were treated as united atoms. According to these models, the total energy of the system was calculated as the sum of the interaction energy of the particle pairs and the intramolecular energy of the DS − ions, the former being the sum of the Lennard-Jones and Coulomb contributions of all atom pairs, and the latter consisting of angle bending and torsional terms. All interactions were truncated to zero beyond the group-based cut-off distance of 15.0 Å; the long-range part of the electrostatic interaction was taken into account using Ewald summation in its smooth particle mesh implementation (sPME) [74]. The geometry of the water molecules as well as the bond lengths of the DS − ions were kept unchanged by means of the SHAKE algorithm [75].
The simulations were performed using an in-house modified version [76] of the GROMACS 5.1 package [77], which allowed us to calculate the lateral pressure contribution of the individual atoms and hence also the lateral pressure profile as well as the contributions of the various molecules and moieties to it. These lateral pressure contributions, including the terms corresponding to the sPME correction of the electrostatic interaction, were calculated employing the Harasima path as described in our previous publication [78]. We recall here that using the Harasima path allows us to account for the lateral pressure contribution corresponding to the sPME correction [69], and it enables us to distribute the lateral pressure contributions among the individual atoms [57]. Equations of motion were solved using integration time steps of 2 fs. The temperature of the systems was kept constant by means of the Nosé-Hoover algorithm [79,80], using a time constant of 1 ps.
Simulations were started from the final configuration of the simulation of the Na + DS − solution reported in our previous paper [56]. The systems were further equilibrated for 10 ns. Density and lateral pressure profiles were then calculated, averaging also over the two liquid surfaces present in the basic box, using 5 × 10 4 equilibrium configurations, separated by 0.4 ps long trajectories each, which resulted from a subsequent, 20 ns long production run. Equilibrium snapshots of the surface portion of the systems consisting of Na + and Cs + counterions are shown in Figure 1a,b, respectively, as taken from our simulations.

Surface Tension and Radial Distribution Functions
The density profiles of the systems simulated along the macroscopic surface normal axis, X, as well as the contributions coming from the water molecules, sulfate headgroups and dodecyl tails of the DS − ions, and counterions to these profiles are shown in Figure 2, while the surface tension

Surface Tension and Radial Distribution Functions
The density profiles of the systems simulated along the macroscopic surface normal axis, X, as well as the contributions coming from the water molecules, sulfate headgroups and dodecyl tails of the DS − ions, and counterions to these profiles are shown in Figure 2, while the surface tension values obtained in the five systems simulated are collected in Table 1. For reference, the surface tension value corresponding to neat SPC/E water [81] is also included in Table 1. It should be noted that it is rather difficult to compare these surface tension values with experimental data corresponding to the same systems because, due to the small size of the systems simulated, practically no surfactant ion has penetrated into the bulk aqueous phase; instead, they are all attached to the liquid surface. As a consequence, the bulk phase concentration of the surfactant cannot be meaningfully defined in our simulations. Nevertheless, it is clear that aqueous dodecyl sulfate solutions corresponding to the more or less saturated surface concentration of the DS − ions of 4 µmol/m 2 must have considerably lower surface tension than neat water. Thus, it is rather striking that, in our simulations, even this qualitative trend is only reproduced when using Na + counterions, while in the presence of all the other counterions considered here the surface tension of the surfactant solution turns out to be considerably higher than that of neat water. It should be emphasized that the KBFF model of alkali cations was developed by testing the properties of aqueous alkali halide solutions using the SPC/E model of water, but it has not been tested in the presence of anionic surfactants [72,82]. Although the GROMOS96 force field, used here to describe the DS − ions, is supposed to be compatible also with the SPC/E water model, our results clearly point out that the KBFF model of the alkali ions is incompatible with the GROMOS96 + SPC/E force field combination. To understand the reasons of the apparent and striking failure of this model combination, we have calculated the radial distribution functions between the cations and the water oxygen atoms, sulfate S atoms, and sulfate O atoms (referred to here as g +Ow (r), g +S (r), and g +Os (r), respectively). The obtained radial distributions are plotted in Figure 3. As seen in the figure, the g +Os (r) functions of all cations, with the exception of Na + , exhibit a huge peak between about 2 and 2.5 Å, indicating the very strong tendency of these cations to form contact pairs with the sulfate group of the DS − ions. A similar picture is seen for the g +S (r) functions, with the only differences being that (i) here the first peak appears at 1-1.5 Å greater distances than in the g +Os (r) functions, reflecting the trivial fact that the O atoms rather than the central S atom of the sulfate group is directly accessible for the counterions, and that (ii) the first peak of g +S (r), with the exception of the case of Li + , is split due to the presence of two O atoms of different constitutional positions in the sulfate group. In the case of both of these radial distribution functions, the results corresponding to the Na + ion exhibit a markedly different behavior, as here the amplitude of the first peak is at least an order of magnitude smaller than for all other alkali counterions. This finding indicates the much weaker tendency of Na + than all the other alkali cations, as described by the KBFF force field, to form contact pairs with the negatively charged group of the DS − ions. Correspondingly, the first peak of the g +Ow (r) function is an order of magnitude higher in the case of Na + than for the larger alkali ions (i.e., K + , Rb + , and Cs + ). More importantly, the first peak of g +Ow (r) of all cations, with the exception of Na + , is an order of magnitude lower than that of g +Os (r), occurring at the same position, while for Na + the opposite relation is observed. This result indicates that while the Na + ion is typically fully surrounded by water molecules and prefers, as expected, to form solvent separated pairs with the DS − headgroups, all other alkali counterions described by the KBFF force field strongly prefer to form contact pairs with the DS − headgroups, or, in other words, they are simply stuck to them. This unphysical behavior of the counterions is clearly seen in the example of Cs + in the snapshot of Figure 1b. This finding is also reflected in the density profiles (see Figure 2). Indeed, the density peaks of the headgroups and counterions occur in the same position along the interface normal axis for all counterions but Na + , in the case of which these two peaks are separated by about 2.7 Å, i.e., roughly the size of a water molecule. The observed unphysical attachment of the counterions to the DS − headgroups is related to the inconsistency of the surfactant (GROMOS) and counterion (KBFF) force field parameters and leads to the unphysical increase of the surface tension of the system in the presence of the alkali dodecyl sulfate surfactant.
Colloids Interfaces 2020, 4, x 6 of 15 cannot be meaningfully defined in our simulations. Nevertheless, it is clear that aqueous dodecyl sulfate solutions corresponding to the more or less saturated surface concentration of the DS − ions of 4 μmol/m 2 must have considerably lower surface tension than neat water. Thus, it is rather striking that, in our simulations, even this qualitative trend is only reproduced when using Na + counterions, while in the presence of all the other counterions considered here the surface tension of the surfactant solution turns out to be considerably higher than that of neat water. It should be emphasized that the KBFF model of alkali cations was developed by testing the properties of aqueous alkali halide solutions using the SPC/E model of water, but it has not been tested in the presence of anionic surfactants [72,82]. Although the GROMOS96 force field, used here to describe the DS − ions, is supposed to be compatible also with the SPC/E water model, our results clearly point out that the KBFF model of the alkali ions is incompatible with the GROMOS96 + SPC/E force field combination.

Lateral Pressure Profiles and Surface Tension Contributions
Despite the apparent failure of the combination of the KBFF and GROMOS96 force fields in properly describing the behavior of the alkali DSsolutions, we have calculated the lateral pressure profiles and surface tension contributions of the different molecules and moieties in these systems. The main reason for this is that, regardless of the unphysical attachment of the oppositely charged ions, the qualitative trend of the surface tension contributions of the counterions of different "hardness" might still be reproduced. In this case, the present results can still be indicative (although by no means decisive) of whether the soft/hard character of the counterions might possibly be related to their surprisingly large surface tension contribution. It should be emphasized that in our previous publication, only the Na + counterion was studied in the systems containing the DS − surfactant ion; furthermore, besides KBFF, the GROMOS96 model of Na + was also considered, and both models led to a proper decrease of the surface tension [56]. Thus, the failure of the KBFF cation model in this respect, reported here, does not invalidate our previous findings at all.
The obtained lateral pressure profile and surface tension contributions of the various molecules and moieties are shown in Figures 4,5, respectively; the surface tension contribution values are also collected in Table 2. The marked oscillation of the headgroup and counterion lateral pressure profiles in the cases of K + , Rb + , and Cs + is probably a consequence of the ordering effect of the observed unphysically strong ion pairing involving these large counterions. Therefore, here we refrain from interpreting this feature.

Lateral Pressure Profiles and Surface Tension Contributions
Despite the apparent failure of the combination of the KBFF and GROMOS96 force fields in properly describing the behavior of the alkali DS − solutions, we have calculated the lateral pressure profiles and surface tension contributions of the different molecules and moieties in these systems. The main reason for this is that, regardless of the unphysical attachment of the oppositely charged ions, the qualitative trend of the surface tension contributions of the counterions of different "hardness" might still be reproduced. In this case, the present results can still be indicative (although by no means decisive) of whether the soft/hard character of the counterions might possibly be related to their surprisingly large surface tension contribution. It should be emphasized that in our previous publication, only the Na + counterion was studied in the systems containing the DS − surfactant ion; furthermore, besides KBFF, the GROMOS96 model of Na + was also considered, and both models led to a proper decrease of the surface tension [56]. Thus, the failure of the KBFF cation model in this respect, reported here, does not invalidate our previous findings at all.
The obtained lateral pressure profile and surface tension contributions of the various molecules and moieties are shown in Figures 4 and 5, respectively; the surface tension contribution values are also collected in Table 2. The marked oscillation of the headgroup and counterion lateral pressure profiles in the cases of K + , Rb + , and Cs + is probably a consequence of the ordering effect of the observed unphysically strong ion pairing involving these large counterions. Therefore, here we refrain from interpreting this feature.     Nevertheless, it is clear that the amplitude of the lateral pressure profile of the counterions decreases gradually and even changes sign by increasing size (and, hence, softness) of the counterion. This trend is accompanied by an opposite change of the lateral pressure profiles of the headgroups. Figure 5 and Table 2 show that the surface tension contribution of the alkali counterions gradually decreases from about 200 mN/m (in the case of the hardest Li + and Na + ions) to almost zero, even indicating a negative contribution in the case of the softest counterion (i.e., Cs + ). At the same time, the opposite trend is seen for the surface tension contributions of the headgroups and, to a smaller extent, also for those of the water molecules. It also turns out that in the presence of the softest counterion, i.e., Cs + , the system behaves qualitatively in the same way (i.e., positive surface tension contributions of the headgroups and water molecules, negative contribution of the counterions) as was previously observed for DTA + Cl − , i.e., a positively charged surfactant and a Nevertheless, it is clear that the amplitude of the lateral pressure profile of the counterions decreases gradually and even changes sign by increasing size (and, hence, softness) of the counterion. This trend is accompanied by an opposite change of the lateral pressure profiles of the headgroups. Figure 5 and Table 2 show that the surface tension contribution of the alkali counterions gradually decreases from about 200 mN/m (in the case of the hardest Li + and Na + ions) to almost zero, even indicating a negative contribution in the case of the softest counterion (i.e., Cs + ). At the same time, the opposite trend is seen for the surface tension contributions of the headgroups and, to a smaller extent, also for those of the water molecules. It also turns out that in the presence of the softest counterion, i.e., Cs + , the system behaves qualitatively in the same way (i.e., positive surface tension contributions of the headgroups and water molecules, negative contribution of the counterions) as was previously observed for DTA + Cl − , i.e., a positively charged surfactant and a negatively charged counterion (see Figure 9 of [56]). It should be emphasized that the system consisting of Na + counterions, which behaves in a qualitatively different (i.e., correct) way than the others, does not deviate substantially from these trends, suggesting that the failure of the model in reproducing the surface tension of the system might indeed not substantially affect the relation of the soft/hard character of the counterions and the surface tension contribution of the different molecules and moieties. The only exception in this respect is the relatively large positive contribution of the dodecyl tails to the surface tension in the presence of Na + . Nevertheless, due to the failure of the model discussed above, the present findings cannot be regarded more than being only qualitative results and need to be confirmed later by using a more reliable potential model that can, at least qualitatively, reproduce the decrease of the surface tension in the presence of surfactants.

Conclusions
In this paper, the relation of the surface tension contribution of the various molecules and moieties with the soft/hard character of the counterions has been intended to be investigated by computer simulations of the liquid-vapor interface of various aqueous alkali dodecyl sulfate solutions. The main and unexpected result of this study is, however, the striking failure of the used model combination of GROMOS96 surfactant ions, KBFF counterions, and SPC/E water in reproducing, at least qualitatively, the decrease of the surface tension in the presence of the ionic surfactant. In this respect, only the system consisting of Na + counterions has been found to behave in a physically reasonable way. It should be emphasized that, ironically, exactly the observed failure of the model combination used for all alkali counterions but Na + clearly proves the extremely important role of the counterions in determining the surface tension of the aqueous solutions of ionic surfactants. The reason for the observed failure has been traced back to the unphysically strong contact pair formation (i.e., attachment) of the counterions and the negatively charged sulfate groups of the surfactants.
The evaluation of the surface tension contribution of the different molecules and moieties has shown a clear trend between these contributions and the soft/hard character of the counterions. More specifically, with increasing size and, hence, softness of the alkali counterions, the surface tension contributions of the counterions gradually decreases, while those of the surfactant headgroups and water molecules increase and even change sign. Thus, in respect to the surface tension contributions, the system consisting of Cs + DS − has been found to behave qualitatively in the same way as was found earlier for the cationic surfactant DTA + Cl − [56]. The fact that these trends are more or less followed in the presence of the Na + counterion, i.e., the only one that does not behave in an unphysical way, suggests that the observed trend might not be strongly affected by the aforementioned failure of the model. Nevertheless, the present findings need to be confirmed by additional simulations, performed with potential models that are able to describe the qualitative behavior of the surface tension of the system upon addition of surfactants. This work might well involve various all-atom [83] and united atom [84] force fields of the DS − ion as well as of the alkali counterions [85,86]. Work in this direction is currently in progress.