Interfacial properties of the ionic liquid [bmim][triflate] over a wide range of temperatures

Abstract

We carried out molecular dynamics simulations of the liquid/vacuum equilibrium of the ionic liquid [bmim][triflate] in a wide range of temperatures (323.15 to 573.15 K). The results showed liquid phases with high densities even at temperatures close to the decomposition temperature of the liquid. The density and surface tension behaviors are linear across this wide range of temperatures, which is an extension of the behaviors of these systems at low temperatures, where these properties have been experimentally measured. The interfacial region shows peaks of adsorption of the ions; they are ordered, with the alkyl chains of the [bmim] cations pointing out of the liquid, and the tailing angle of the chains becomes 90° at higher temperatures. The alkyl chains are part of the outermost interfacial region, where intra- and intermolecular tangential forces are in equilibrium; thus, they do not contribute to the total surface tension. Unlike simpler organic liquids, the surface tension is composed of positive normal contributions of intermolecular interactions; these are almost in equilibrium with the negative normal contributions of intramolecular interactions, which are mainly vibrations of the distance and the angle of valence. The pressure profiles show that the molecules are in ‘crushed’ conformations internally in the bulk liquid and even more so in the normal direction at the interface. The total pressure profiles show values with physical meaning, where the tangential peaks show higher values than normal pressures and give rise to the surface tension. Short cutoff radii for the calculation of intermolecular forces (less than 16.5 Å) produce a system that is not mechanically stable in the region of the bulk liquid (confirmed by radial distribution function calculations); this produces a difference between the normal pressure and the average of the tangential pressures, which affects the calculation of the surface tension due to overestimation by up to 20% when using the global expression, which is extensively used for the calculation of surface tension. The use of a sufficiently long cutoff radius avoids these mechanical balance problems.

---

Introduction

Ionic liquids combine relatively large organic cations and organic and inorganic anions; this prevents the formation of ordered crystals over a wide range of temperatures, including room temperature. These liquids are more complex than ordinary organic liquids, and it is thought that they can be best described as fluids that self-assemble into amphiphilic nanostructures.^1–3 Simulations of molecular dynamics in the liquid/vacuum equilibrium using 368 pairs of ions with a simulation area of ∼20 nm² and 5 ns for sampling indicate the presence of amphiphilic interfaces. This is due to ordering in the interfaces because the alkyl chains of diverse ionic liquids based on imidazole leave the bulk liquid, pointing towards the vacuum; this is similar to the behavior of complexes that form water bridges with amphiphilic molecules.⁴ Other simulations, with smaller samples (216 pairs of ions), a smaller surface area (∼17 nm²), and 2 ns for sampling, produce similar results.⁵

On the other hand, the thermophysical properties of ionic liquids are unique. Their vapor pressure is extremely low, and the dependence of their surface tension, γ, on temperature is different from that of most organic liquids (γ∝T^11/9). Ionic liquids behave almost linearly with respect to temperature; therefore, the linear van der Waals–Guggenheim equation, which describes the behavior of molten metals, also accurately describes the behavior of many ionic liquids.⁶ Some exceptions to this linear behavior have been attributed to problems arising during the preparation of the liquids, such as insufficient purification.^6,7 The behavior of γ with respect to the temperature of the ionic liquid [bmim][triflate] has been studied experimentally by Součková et al. in the temperature range of 292.16 to 356.00 K,⁸ by Freire et al. in the range of 292.20 to 343.20 K,⁹ and by Cione et al. at 298.15 K.¹⁰ All the reported sets of data show linear behavior. Freire et al.⁹ used force measurements (the du Noüy ring method) to obtain their data, and the purity of the employed mixture was >98%. Meanwhile, Součková et al.⁸ also used force measurements; however, they employed two methods (the du Noüy ring and Wilhelmy plate methods). They reported that the data acquired from the du Noüy ring method showed seven times more scatter than those from the Wilhelmy plate method; the purity of their mixture was >99.5%.

The density of ionic liquids, ρ_L, unlike many organic liquids, also presents linear behavior with respect to temperature. The behavior of ρ_L with respect to the temperature of the ionic liquid [bmim][triflate] has been studied experimentally by several groups, and data have been reported for a temperature range of 288.15 to 343.15 K at ambient pressure.^11–14 The effects of pressure on ρ_L have been studied by Gardas et al.¹⁵ They found that, as with many common liquids, there is a weak dependence of ρ_L on pressure; thus, when the pressure changes from 0.10 to 10.0 MPa, the density of the ionic liquid only increases by 0.56% at 373.15 K.

Many measurements of ρ_L and γ have only been studied over a small range of temperatures and pressures; therefore, more experimental and simulation measurements are necessary in order to understand the behavior of these fluids over a wide range of temperatures and pressures. It is also necessary to obtain more information on how the arrangement of the alkyl chains of the [bmim] cations contributes to the interface, and their possible effects on the interfacial properties.

In this paper, we report ρ_L and γ in the liquid/vacuum equilibrium region for a wide range of temperatures, from 323.15 to 573.15 K, which approach the thermal stability limit. We also explain how the molecules are organized at the interface, as well as how the components of the pressure tensor affect the surface tension.

The methods and parameters of the interactions used in this work are described in the Methodology section. In the Results section, we discuss the main results, in terms of the total density profiles, by component; the profiles of the concentration of sites in the interface and their structuring and ordering; and the profiles of pressure and surface tension. We also present an analysis of the main contributors to the surface tension. The Conclusions section summarizes our findings.

Methodology

The liquid/vacuum interface of the ionic liquid [bmim][triflate] was studied directly through the simulation of a thin layer of liquid (thickness of more than 100 Å) surrounded by a vacuum^16–18 using the molecular dynamics method at various temperatures of between 323.15 and 573.15 K. The simulation cell consisted of a parallelepiped with dimensions of 94 × 94 × 180 Å and an interfacial area of 88.36 nm² which contained 2592 pairs of ions (85 536 atoms). The initial system consisted of a cubic cell, to simulate a single liquid phase, which contained half of the ion pairs under the conditions of 298.15 K and 0.1 MPa once equilibrium was reached. This cell was repeated once in the inhomogeneous direction and was left to balance. Finally, vacuum was added in the inhomogeneous direction, and the portions of the molecules that were not entirely located at the interfaces due to the initial periodic conditions shifted. The system was initially equilibrated at 573.15 K in order to avoid the difficulty of structuring of the liquid phases at low temperatures, which requires a long simulation time to equilibrate.^19,20 Systems at lower temperatures were obtained by cooling the system at a speed of 10 K ns⁻¹. The system was simulated using the NVT assembly (number of atoms, volume, and temperature constants) using the Nosé thermostat,²¹ implemented in the open source for large-scale atomic/molecular massively parallel simulator (LAMMPS),²² with a time step of 1 fs.

The components of the ionic liquid were simulated using the models of the ions [bmim] and [triflate] developed by Cadena et al.²³ and Lopes et al.,²⁴ respectively. The molecular ion models are flexible and include intramolecular interactions due to vibrations of the bonding distance, valence, dihedrals, and improper angles. The intermolecular interactions of the original model of the [bmim] cation include truncated Lennard-Jones and shifted potential through a switching function, with an internal cutoff radius r_ci and an external cutoff radius r_ce of 10.5 and 12 Å, respectively. The original model of the [triflate] anion includes the Lennard-Jones truncated potential, with cutoff radii of 10 and 12 Å depending on the state of the system; long-range corrections were applied. The potentials from these two ions have been optimized in ionic liquids with other ions, but not for the [bmim][triflate] liquid; however, no other specific potential for this liquid has been reported in the literature, so we used these potentials as the best available approximation. Cione et al. also used this potential approximation to study how [bmim][triflate] ionic liquid wetted and deposited on silica coated with alkyl chains.¹⁰ In this study, we used the truncated and displaced potential without long-range corrections of the Lennard-Jones forces to ensure the compatibility of the two models. We studied the effects of r_ci on the calculated properties, using values for r_ci of 10.5, 12.5, …, 30.5 Å, while r_ce = r_ci + 1.5 Å. The maximum cutoff radius employed (2 × 32 Å), plus 10% for the neighbor list (64.0 + 6.4 Å), is smaller than the length of the simulation cell in the tangential direction, which allowed us to correctly sample the intermolecular interactions. The use of a large cutoff radius requires the computation of larger numbers of intermolecular interactions, an estimation which can be obtained by calculating the relationship between the volumes of spheres at different r_ce values, , which are directly related to the number of particles and, ultimately, to the number of intermolecular interactions (assuming no intramolecular interactions). For a shorter cutoff radius, the computational effort is considerably lower: . The electrostatic interactions were modeled using the 3-D particle–particle, particle–mesh method,²⁵ using the r_ce of the Lennard-Jones potential. In the inhomogeneous direction, sufficient space was added to avoid interactions between periodic liquid layers.

The pressure profiles were calculated with the tension tensor that was used to calculate the pressure profiles of Harasima,²⁶ implemented in LAMMPS,²⁷ where contributions to the profiles are distributed only in the two slabs that originate the interactions and not uniformly throughout all the slabs between those two slabs.^28,29 The previous procedure were applied to all intra- and intermolecular interactions. The normal pressure profile obtained from this definition was not uniform along the interface, which was thought to make no physical sense, even though it is not experimentally feasible to validate the measurements of the components of the pressure at the interface. Despite this, the pressure profile of Harasima²⁶ has been used in several studies of surface tension because when it is integrated to obtain the profile of the surface tension, the same result is obtained whether the interaction forces are distributed uniformly in several slabs or in only two slabs.^30–33 As contributions to the Harasima²⁶ pressure profiles are specifically located, regions that contain inhomogeneities in their density profiles are highlighted in the pressure profiles. γ was calculated through its mechanical definition:^29,34

where A is the surface area, and z₁ and z₂ are the integration limits that, if they cover the entire simulation cell, require γ to be divided by two because there are two interfaces in the system. P_N and P_T are the normal and tangential pressures, respectively, and z represents the inhomogeneous direction of the simulation cell.

Results

Density profiles

Simulations of the molecular dynamics of the ionic liquid [bmim][triflate] in an inhomogeneous arrangement in order to study the liquid/vacuum equilibrium produced the density profiles shown in Fig. 1a. The profiles show the formation of a liquid layer and its interface, with thicknesses of 112 Å (373.15 K) and 129 Å (573.15 K). The profiles show liquid regions in the bulk that range between 90 (373.15 K) and 106 Å (573.15 K), with adsorption peaks in the liquid/vacuum interface that are commonly observed in binary systems.^35,36 These are highlighted in Fig. 1a as shaded areas above the adjustment to the next equation, which commonly describes simple systems with a single component or binary and ideal systems where there are no molecules in the vapor phase and where the density at the interface decays monotonically:^20,29,34where z₀ is the position of the Gibbs dividing surface and d is the thickness of the interface divided by 2.7192.¹⁷ In the obtained density profiles, none of the ions formed a vapor phase at any time during the simulation (∼60 ns). If we consider that both ions are present in the liquid/vacuum interface (which is shown below), we find that, when integrating the density profiles on the adsorption regions shown in Fig. 1a, a pair of [bmim][triflate] ions is present every 12 (373.15 K) or 34 nm² (573.15 K). These additional concentrations of ions at the interface demonstrate how different these systems are from ideal systems, where the concentration drops monotonically at the interface. An analysis of the dependence of the saturated ρ_L in bulk with respect to the r_ci used in the simulation shows that ρ_L varies slightly with r_ci (Fig. 1b), presenting an underestimation of ∼1.4% at r_ci = 10.5 Å when compared with the calculated value using r_ci = 30.5 Å. This is 7.5 times higher than the largest Lennard-Jones parameter (among all the values used) associated with the equilibrium separation between two Lennard-Jones sites. In previous studies, it has been found that this r_ci is long enough to reproduce the interactions using the total Lennard-Jones potential.³⁷

Fig. 1 (a) Total density profiles as a function of the position on the inhomogeneous axis of the ionic liquid [bmim][triflate]. Black lines represent the results of the simulation; red lines represent the results of the adjustment to eqn (2). (b) Total density of the liquid as a function of the cutoff radius. The solid line represents the best fit to an exponential decay function. (c) Total density of the liquid as a function of temperature. Red dots represent experimental values;¹⁴ crosses represent the results of this work (errors calculated as the standard deviation are smaller than the symbols). The dashed line represents values of the empirical correlation of Tokuda et al.¹⁴ The y-axes of the three graphs represent density in the same units.

A graph illustrating the results of saturated ρ_L in bulk over a temperature range of 323.15 to 573.15 K using r_ci = 30.5 Å is shown in Fig. 1c. The simulation results of this work over a wide range of temperatures show a linear dependence of ρ_L on temperature; the same trend occurred in the reduced temperature range that was studied experimentally by Tokuda et al.¹⁴ (287.15 to 313.15 K). Linear extrapolations of the experimental data of Tokuda et al.¹⁴ coincide very well with the results of this work; at 573.15 K, the simulation overestimates the correlation density data by only 0.05%. The linearity of ρ_L over the wide range of temperatures studied, which reached 40 K below the thermal stability limit (613.15 K),^12,38 suggests that these temperatures are very far from the critical temperature and that they do not influence the rapid change in saturated ρ_L with temperature.²⁹ This effect is similar to that observed for the linearity of saturated ρ_L in long linear alkanes, such as n-hexadecane.^39–41 Even so, we do not rule out additional effects, such as the nano-structuring of the system, which induces effects such as an unusual linear dependence of the viscosity on temperature in mixtures with water.⁴²

The individual density profiles of each of the ions showed ordering at the interface and, to a lesser extent, in the bulk liquid, which did not disappear when we increased r_ci (Fig. 2). This ordering has previously been observed in simulations of ionic liquids of [bmim] and other anions;⁴ it has also been experimentally reported in X-ray reflectivity measurements, from which the formation of intercalated ion layers has been inferred.⁴³ We did find when increasing r_ci that the peaks of adsorption at the interface were more pronounced with increasing r_ci, which is probably due to an increase in cohesiveness because more intermolecular forces are taken into account in the Lennard-Jones potential.

Fig. 2 Density profiles for each of the components of the ionic liquid [bmim][triflate] at a temperature of 373.15 K and at r_ci values of (a) 10.5 Å and (b) 30.5 Å. Black and red lines represent the [bmim] and [triflate] ions, respectively.

Structuration at the interface

The atomic concentration profiles of key sites in the structures of the ions are shown in Fig. 3a and b at temperatures of 373.15 and 573.15 K. For the [bmim] ion, the profiles of the carbons corresponding to the methyl attached to the cycle (C₁) and the terminal carbon of the alkyl chain (C₂), are plotted, while for the [triflate] ion, the profiles of the central atoms, C and S, are plotted. The profiles show a pronounced peak of adsorption at the interface for the C₁ carbons of [bmim] at 373.15 K, surrounded by two adsorption peaks in the profile of the C₂ carbons. At 573.15 K, the magnitude of the peak of the C₁ carbons is reduced, and the inner peak of the C₂ carbons disappears. The same occurs with the [triflate] ion at 373.15 K, where a broad adsorption peak, representing the S atoms, is surrounded by adsorption peaks of the C atoms; meanwhile, at 573.15 K, the peaks of the carbons disappear, and the peak of S also almost disappears. Enlarged views of the atomic concentration profiles of the carbons of [bmim] are shown in Fig. 3c and d at the same temperature. Integration of the interfacial area shows that the concentration of C₁ sites (which is the same as that of the [bmim] ions) decreases from 4.00 to 3.33 [bmim] ions per nm² when the temperature is increased from 373.15 to 573.15 K. On the other hand, when calculating the separation between the peaks of the profiles of the C₁ and C₂ carbons, we can infer the orientations of the [bmim] ions at the interface. At 373.15 K, the [bmim] ions are oriented with their aliphatic chains forming a tailing angle of ∼45° with respect to the interface, pointing away from the bulk liquid, while at 573.15 K, the tailing angle becomes almost perpendicular to the interface. This amphiphilic behavior has been previously found in simulations of ionic liquids using cations [bmim] and other anions using r_ci = 10.5 Å.⁴ Silicon oxide coated on alkyl, alcoholic, and fluoroalkane chains show similar dependencies of their tailing angles when external forces are introduced.^44,45

Fig. 3 (a) Atomic concentration profiles of the carbons at the ends of the [bmim] ion, C₁ and C₂ (see figure), as a function of their positions on the inhomogeneous axis. Black and red lines represent the values at 373.15 and 575.15 K, respectively. Solid lines represent the C₁ profile; dotted lines represent the C₂ profile. Blue lines represent the values of the outer peaks of site C₂. (b) Atomic concentration profiles of the C and S atoms of the [triflate] ion (F₃C–SO₃)–. Solid lines represent the profile of C; dotted lines represent the profile of S. (c and d) Expanded views of the atomic profiles of the [bmim] cation at 373.15 and 575.15 K, respectively. Black and red lines represent the C₁ and C₂ carbons, respectively. Blue lines represent the position of the liquid/vacuum interface. The orientations of the [bmim] cations with respect to the interface are also depicted in the figures.

The imidazolium ring of the [bmim] cation is planar, and the angle between the imidazolium ring plane and the interfacial plane showed a distribution of angles of around 60° in a range between 0 and 150° (Fig. 4). Temperature does not have a considerable effect on the distribution of the angles; the mean value changes by only a few degrees when the temperature is increased from 373.15 to 573.15 K, which is probably due to the structuration of the liquid phase. The distribution of the angles is consistent with the conformation that the alkyl chains take in the interface; however, a direct correlation cannot be obtained because the angle ‘imidazolium ring-N–C’ at the junction between the imidazolium ring and the alkyl group can rotate at the cost of a small energetic barrier.

Fig. 4 Probability distribution of the angles between the planes of the imidazolium ring of the [bmim] cation and the interface liquid/vacuum plane. Black and red dots represent distributions at 373.15 and 573.15 K, respectively.

Cutoff radius effects on pressure profiles

The pressure profiles, calculated using eqn (1) and the original r_ci and r_ce of the potential [bmim][triflate] at 573.15 K, are shown in Fig. 5a. For these, r_ci and r_ce and P_T and P_N showed differences in the bulk liquid, which should not occur in a system in mechanical balance. The P_T values were homogeneous in the region associated with the bulk liquid (−50 to 50 Å), while P_N showed structuring. When calculating the difference between P_N and the average P_T, the structuring was maintained and, more importantly, the average density along the liquid phase in bulk was not zero (an average value of zero would be expected in a system in mechanical equilibrium); however, a positive average pressure difference of 0.92 MPa (±1.73 MPa) was observed, which implies that the system was not in mechanical equilibrium in the bulk liquid. We increased r_ci and r_ce by increments of 2 Å until we reached values of 30.5 and 32 Å, respectively; we found that from a r_ci of 16.5 Å, P_N lost its structuring, and the difference between P_N and the average P_T attained an average value of zero.

Fig. 5 (a) Tangential (black and red) and normal (green) pressure profiles, and the difference between the normal pressure profile and average tangential profiles (blue), for the [bmim][triflate] system at 573.15 K along the inhomogeneous axis with a cutoff radius of 10.5 Å. (b) Surface tension profiles (eqn (1)) with cutoff radii of 10.5 Å (black) and 16.5 Å (red).

The use of r_ci = 10.5 Å also affected the calculation of γ, not only because the Lennard-Jones interactions were truncated and their contributions to γ were not complete, but also because the artificial mechanical instability created in the bulk liquid caused the profile of γ to increase continuously in the region of the bulk liquid, which should not occur in a system in mechanical equilibrium (Fig. 5b). Starting with a r_ci of 16.5 Å, the average value of the pressure difference became zero, and the contributions to γ in the bulk liquid were zero (Fig. 5b).

To save computational effort, γ is commonly calculated through the global expression

where L_z represents the dimension of the simulation cell along the inhomogeneous axis z, 〈·〉 represents temporary averages over periods on the ns scale, and the number 2 is due to the existence of 2 liquid/vacuum interfaces. In equilibrium systems, the value obtained from this global expression should be the same as that obtained by integrating the mechanical expression along the entire inhomogeneous axis (eqn (1)); however, in the global expression, it cannot be deduced whether the system is in mechanical equilibrium in the bulk liquid. Values of the global expression (eqn (3)) for the calculation of γ appear in Fig. 6, as well as values of γ using its mechanical definition (eqn (1)), integrating from −∞ to −50 Å (internal limit of the left interface) using a r_ci < 16.5 Å. Starting at r_ci = 16.5 Å, the global expression and the mechanical definition produce the same values, and only one dataset is shown. For r_ci = 10.5 Å, we observed differences in the calculation of γ of up to 20% between the global expression and the integration only on the left interface, without including the bulk liquid. The values obtained from integrating only on the left interface (r_ci < 16.5 Å) show a continuous profile with the values obtained from integrating on both interfaces (r_ci ≥ 16.5 Å), showing a maximum for r_ci ∼ 24.5 Å; later, the values show the beginning of what can be considered to be an oscillation zone, but with small differences between the peaks and valleys of ∼1.6 mN m⁻¹. The oscillating behavior of the surface tension of simpler models of ionic liquids has been previously studied and has been associated with the size of the simulation cell, which can be indirectly associated with the cutoff radius.^46,47

Fig. 6 Surface tension as a function of the cutoff radius. Open circles represent values calculated using their mechanical definition, integrating the pressure difference profile only in the left interfacial region (eqn (1)), while the closed circles represent integration over the entire system. The crosses represent values calculated using the global expression (eqn (3)).

Structuration of the bulk liquid

We calculated the radial distribution functions (RDF) of two key sites of the [bmim] ions in the bulk liquid at 573.15 K, created in the space between −50 and 50 Å in the inhomogeneous direction (Fig. 7). The RDF of one of the nitrogens in the imidazolium ring, which is the most hindered atom of the [bmim] ion, and to which the alkyl chain is attached, showed two prominent peaks at ∼8.1 and ∼15.2 Å. The second peak is larger than the cutoff radius employed in the original optimization of the potential parameters of the [bmim] and [triflate] ions; this indicates that this specific combination of ions requires a longer cutoff radius, as shown in the calculation of the pressure tensor. The RDF of the C₂ carbon located at the end of the alkyl chain, which is the most accessible atom (excluding the hydrogens), showed one prominent peak at ∼4.2 Å and a small peak at ∼10.5 Å, which are both within the cutoff radius of the original parameterization. To verify if the larger cutoff radius employed in this work affected the results of the homogeneous liquid phase, the RDF values of the two key components were calculated using cutoff radii of r_ci = 10.5 Å and r_ci = 30.5 Å. We found that the positions of the peaks did not change, and the peaks had almost the same height; this indicates that the use of a larger cutoff radius did not affect the structural properties of the bulk liquid phase.

Fig. 7 Radial distribution functions of two key sites in the [bmim] ion in the [bmim][triflate] ionic liquid at 573.15 K, which correspond to (a) one of the nitrogens of the imidazolium ring where the alkyl chain is attached, and (b) the terminal carbon of the alkyl chain, C₂.

Cutoff radius effects on surface tension

Fig. 8 shows the results of the calculation of γ at temperatures of 373.15, 473.15, and 573.15 K and r_ci = 30.5 Å, which agree well with linear extrapolations of the experimental results of Součková et al. at lower temperatures.⁸ Linear extrapolations of other experimental data have produced considerable deviations from our simulation results, probably due to the methods employed and the purity of the mixtures.⁹ Linear extrapolations were used in the comparisons because experimental measurements of γ with respect to temperature in other ionic liquids have shown good fit to the linear van der Waals–Guggenheim equation,⁴⁸ which also describes the behavior of molten metals. At all studied temperatures, the simulations showed a maximum overestimation of ∼0.56 mN m⁻¹ (573.15 K), which is equivalent to a positive deviation of ∼2.26%. This indicates that the saturated ρ_L obtained in the simulations represents an overestimation because γ is a strong function of the saturated ρ_L. The linearity of γ with respect to temperature could also be assumed to be a consequence of the linearity of the saturated ρ_L that has been reported experimentally and in the simulation data acquired here. This is ultimately a consequence of the structuring of the system not only in the interface, but in the bulk liquid phase,^42,43 although a long separation from the critical temperature may influence the linearity of this property with temperature. We also studied the properties of the system at a temperature of 298.15 K. Even though the density profiles showed a liquid phase in equilibrium, we failed to produce symmetric pressure tensors using simulation periods of up to 60 ns. At this low temperature, the liquid tends to be more structured, and we probably need more time to equilibrate the system or design a strategy to cross the local equilibrium shown by one of the interfaces.

Fig. 8 Surface tension of the [bmim][triflate] system as a function of temperature. Blue and red circles represent the experimental results of Součková et al.⁸ and Freire et al.,⁹ respectively, while the crosses represent the results of this work. Dashed lines represent linear extrapolations of the experimental data.

Contributions to surface tension

Fig. 9a shows the total pressure profiles (including all intra- and intermolecular, and kinetic contributions) of P_N and the two P_T values as well as the difference between P_N and the average P_T using r_ci = 30.5 Å at 573.15 K. Using this long r_ci, no structuration appeared in the bulk liquid phase, and the profile of the pressure differences in the bulk liquid showed an average value of zero (mechanical stability). Fig. 9d shows the integration profiles using the mechanical definition of γ (eqn (1)) for the total contributions as well as for the contributions of the inter- and intramolecular interactions, and the contribution due to the kinetic part of the pressure profile. The intermolecular contribution of γ shows large positive values (net contribution of ∼170 mN m⁻¹). This would commonly be expected from this type of interaction, which is very dense, because the profile of differences between normal and average tangential pressures results in a profile with peaks with positive values at the interface and in the bulk liquid; however, it is interesting that the intermolecular profile, P_N, at the interface showed peaks with values greater than the values of P_T. The opposite would be expected (tangential intermolecular forces greater than normal forces, resulting in tension forces at the interface); thus, in the normal direction, the molecules in the interface experience stronger forces than those that occur in the bulk liquid, while the tangential forces are smaller. The fact that the pressures in all directions are positive in the regions of the bulk liquid and at the interface (Fig. 9b) indicates that, on average, the intermolecular separations are shorter than the separation corresponding to the minimum energy in the intermolecular potential (Lennard-Jones + coulombic). Thus, the intermolecular pressures are positive (repulsive), unlike what occurs in a simpler fluid, such as ethane, where the intermolecular pressures are negative (attractive) in all directions.²⁹

Fig. 9 Tangential (black and red) and normal (green) pressure profiles, and the difference between the normal pressure profile and the average tangential profile (blue), for the [bmim][triflate] system at 573.15 K over the inhomogeneous axis with a cutoff radius of 30.5 Å and the contributions to the (a) total, (b) intermolecular, and (c) intramolecular pressure profiles. (d) Surface tension profiles with a cutoff radius of 30.5 Å for the total (black), intermolecular (red), intramolecular (green), and kinetic (blue) contributions. The vertical orange lines represent the positions where the different profiles begin to grow. (e) Surface tension profiles with a cut-off radius of 30.5 Å for the distance vibration (black), valence angle (red), and dihedral (blue) contributions.

For a simple fluid such as ethane, which is modeled as two flexibly-bonded Lennard-Jones sites, it can be concluded that the intermolecular separations are, on average, longer (attractive) than the separation corresponding to the minimum energy of the intermolecular potential (Lennard-Jones alone); however, the ions of the ionic liquid [bmim][triflate] are more complex because they also include coulombic interactions, and each interaction site has different parameters. Therefore, it is complicated to infer an average trend of how the intermolecular separations of the different sites in the interface produce this unexpected behavior.

A similar analysis of the intramolecular contributions (Fig. 9d) shows that these resulted in a negative contribution to γ (net contribution of about −148 mN m⁻¹). The pressure profiles are negative (Fig. 9c), and they also show an inverse behavior (P_N > P_T in the interface); thus, we can infer that the molecules in the interface are in more extended states with respect to their values in equilibrium in the normal direction, which may be a consequence of what happens with the intermolecular pressures in that direction. Fig. 9e shows the γ profiles of the intramolecular interactions; it was observed that the main contributors to the negative value obtained are the vibration interactions of the bond distance and the valence angle. Together with the pressure profiles, we noted that the conformations of the ions in the interface have extended bonds and more open bond angles than the equilibrium values in states in which it would seem that the molecules are ‘crushed’ in the bulk liquid, as well as at the interface. This may be a consequence of the high density of the ionic liquid.

The kinetic contribution is minimal, contributing ∼0.22 mN m⁻¹ to the total γ. The sum of all the contributions produces profiles of P_N and P_T that have physical meaning, where the tangential forces are greater than the normal forces in the interface, giving rise to the surface tension of the ionic liquid. The integration profiles of γ show that the inter- and intramolecular interactions begin to increase around −71 Å, while the total profile starts to increase around −62 Å; this indicates that the inter- and intramolecular contributions are almost completely canceled out in this outer interfacial region between −71 and −62 Å, which extends for ∼9 Å. As we saw in the atomic concentration profiles, this region is composed mainly of the alkyl chains of [bmim] cations and free [triflate] anions.

Conclusions

We simulated the liquid/vacuum equilibrium of the ionic liquid [bmim][triflate] over a wide range of temperatures, from near its freezing point to 40 K below its point of thermal instability. The liquid presented linear dependencies of the density of the saturated liquid and its surface tension with respect to the temperature of the system.

The alkyl chains of the [bmim] cations were structured in, and pointing out of, the interface, with a tailing angle that became straighter as the temperature increased. Alkyl chains were highly compressed at the interface, causing an additional accumulation at the interface of a pair of [bmim][triflate] ions every 34 nm² at 573.15 K. The outer region of the interface, composed of the structured alkyl chains and free [triflate] anions, did not contribute to the total surface tension; however, this region is likely to contribute significantly when pollutants or moisture accumulate within it, interacting with the alkyl chains and creating additional tension forces. The systems did not show any vaporized molecules, which may be due to this singular outer region of the interface, which is characterized as a region in balance between its intra- and intermolecular forces. Probably, only the intervention of other molecules (surfactants) or the presence of external forces would facilitate the vaporization (mobilization) of a pair of ions, creating an imbalance between the canceled forces present in that region.⁴⁹

The use of a short internal cutoff radius (<16.5 Å) in the [bmim][triflate] system caused structuring of the bulk liquid phase in the normal direction of the interface at all simulated temperatures; additionally, it created mechanical instability due to the existence of a difference between the tangential and normal pressures in the bulk liquid, resulting in inadequate calculation of the surface tension when the global expression was used. The use of a sufficiently long cutoff radius (≥16.5 Å) for the [bmim][triflate] system avoided these problems. Any calculation of the surface tension should involve an analysis of the pressure profiles to avoid artefacts caused by a limited cutoff radius. The radial distribution functions of key sites in the [bmim] structure also show peaks at positions larger than the cutoff radius employed in the original parameterizations; this necessitates a larger cutoff radius to properly account for these characteristic regions, where the radial distribution is larger than the value distribution.

The positive contributions of the intermolecular interactions to the surface tension reached values of up to 6.8 times the value of the total surface tension. Because the liquid is a very dense system, the intermolecular pressures showed positive values (repulsive) in the bulk liquid and the interface; this presented unexpected behavior, where the normal pressures were greater than the tangential pressures. The intermolecular contributions to surface tension were compensated by negative contributions, mainly due to the intramolecular interactions of the vibrations of the bond and the angle of valence. This showed that at the interface, the molecules were more extended (crushed) than in the bulk liquid, which is probably a consequence of repulsive intermolecular separations. This may explain the almost null value of the vapor pressure, because the molecules under these normal repulsive pressures were strong and only the intermediation of the intramolecular forces kept them attached to the interface. The total pressure profiles show values that have physical meaning, where the tangential pressures were stronger than the normal pressure and were the origin of the surface tension.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

JLR and PNS gratefully thank CONACYT (México) for infrastructure fellowships. JLR and MRE gratefully thank UMSNH for research funds.

References

1. R. Hayes, G. G. Warr and R. Atkin, Chem. Rev., 2015, 115, 6357–6426.
2. C. R. López-Barrón, D. Li, L. DeRita, M. G. Basavaraj and N. J. Wagner, J. Am. Chem. Soc., 2012, 134, 20728–20732.
3. J. Bowers, M. C. Vergara-Gutierrez and J. R. P. Webster, Langmuir, 2004, 20, 309–312.
4. M. E. Perez-Blanco and E. J. Maginn, J. Phys. Chem. B, 2010, 114, 11827–11837.
5. R. M. Lynden-Bell and M. Del Pópolo, Phys. Chem. Chem. Phys., 2006, 8, 949–954.
6. M. Tariq, M. G. Freire, B. Saramago, J. A. P. Coutinho, J. N. C. Lopes and L. P. N. Rebelo, Chem. Soc. Rev., 2012, 41, 829–868.
7. I. Bandrés, B. Giner, H. Artigas, C. Lafuente and F. M. Royo, J. Chem. Eng. Data, 2009, 54, 236–240.
8. M. Součková, J. Klomfar and J. Pátek, Fluid Phase Equilib., 2011, 303, 184–190.
9. M. G. Freire, P. J. Carvalho, A. M. Fernandes, I. M. Marrucho, A. J. Queimada and J. A. P. Coutinho, J. Colloid Interface Sci., 2007, 314, 621–630.
10. A. M. Cione, O. A. Mazyar, B. D. Booth, C. McCabe and G. K. Jennings, J. Phys. Chem. C, 2009, 113, 2384–2392.
11. G. García-Miaja, J. Troncoso and L. Romaní, J. Chem. Thermodyn., 2009, 41, 161–166.
12. C. P. Fredlake, J. M. Crosthwaite, D. G. Hert, S. N. V. K. Aki and J. F. Brennecke, J. Chem. Eng. Data, 2004, 49, 954–964.
13. M. L. Ge, R. S. Zhao, Y. F. Yi, Q. Zhang and L. S. Wang, J. Chem. Eng. Data, 2008, 53, 2408–2411.
14. H. Tokuda, K. Hayamizu, K. Ishii, M. A. B. H. Susan and M. Watanabe, J. Phys. Chem. B, 2004, 108, 16593–16600.
15. R. L. Gardas, M. G. Freire, P. J. Caryalho, I. M. Marrucho, I. M. A. Fonseca, A. G. M. Ferreira and J. A. P. Coutinho, J. Chem. Eng. Data, 2007, 52, 80–88.
16. J. L. Rivera, C. McCabe and P. T. Cummings, Phys. Rev. E, 2003, 67, 11603.
17. J. L. Rivera, M. Predota, A. A. Chialvo and P. T. Cummings, Chem. Phys. Lett., 2002, 357, 189–194.
18. J. L. Rivera, C. McCabe and P. T. Cummings, Nano Lett., 2002, 2, 1427–1431.
19. J. Alejandre, J. L. Rivera, M. A. Mora and V. de la Garza, J. Phys. Chem. B, 2000, 104, 1332–1337.
20. F. M. Juárez-Guerra, J. L. Rivera, A. Zúñiga-Moreno, L. A. Galicia-Luna, J. L. Rico and J. Lara, Sep. Sci. Technol., 2006, 41, 261–281.
21. S. Nosé, J. Chem. Phys., 1984, 81, 511.
22. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
23. C. Cadena and E. J. Maginn, J. Phys. Chem. B, 2006, 110, 18026–18039.
24. J. N. C. Lopes and A. A. H. Pádua, J. Phys. Chem. B, 2004, 108, 16893–16898.
25. R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles, Taylor & Francis Group, New York, 1988.
26. A. Harasima, Adv. Chem. Phys., 1958, 1, 203–237.
27. T. W. Sirk, S. Moore and E. F. Brown, J. Chem. Phys., 2013, 138, 64505.
28. J. L. Rivera, F. W. Starr, P. Paricaud and P. T. Cummings, J. Chem. Phys., 2006, 125, 94712.
29. J. L. Rivera, H. Nicanor-Guzman and R. Guerra-Gonzalez, Adv. Condens. Matter Phys., 2015, 2015, 258601.
30. T. Dreher, C. Lemarchand, L. Soulard, E. Bourasseau, P. Malfreyt and N. Pineau, J. Chem. Phys., 2018, 148, 34702.
31. M. Sega, B. Fábián, A. R. Imre and P. Jedlovszky, J. Phys. Chem. C, 2017, 121, 12214–12219.
32. A. P. Sgouros, G. G. Vogiatzis, G. Kritikos, A. Boziki, A. Nikolakopoulou, D. Liveris and D. N. Theodorou, Macromolecules, 2017, 50, 8827–8844.
33. M. Orsi and J. W. Essex, Faraday Discuss., 2013, 161, 249–272.
34. J. L. Rivera, F. W. Starr, P. Paricaud and P. T. Cummings, J. Chem. Phys., 2006, 125, 94712.
35. A. V. Lukyanov and A. E. Likhtman, Langmuir, 2013, 29, 13996–14000.
36. J. L. Rivera, J. Alejandre, S. K. Nath and J. J. De Pablo, Mol. Phys., 2000, 98, 43–55.
37. A. Trokhymchuk and J. Alejandre, J. Chem. Phys., 1999, 111, 8510–8523.
38. C. Maton, N. De Vos and C. V Stevens, Chem. Soc. Rev., 2013, 42, 5963–5977.
39. S. Outcalt, A. Laesecke and T. J. Fortin, J. Chem. Thermodyn., 2010, 42, 700–706.
40. T. J. Ford, Ind. Eng. Chem. Fundam., 1986, 25, 240–243.
41. D. J. Rosenthal and A. S. Teja, AIChE J., 1989, 35, 1829–1834.
42. R. Nanda, Phys. Chem. Chem. Phys., 2016, 18, 25801–25805.
43. N. Nishi, Y. Yasui, T. Uruga, H. Tanida, T. Yamada, S. I. Nakayama, H. Matsuoka and T. Kakiuchi, J. Chem. Phys., 2010, 132, 164705.
44. J. Ben Lewis, S. G. Vilt, J. L. Rivera, G. K. Jennings and C. McCabe, Langmuir, 2012, 28, 14218–14226.
45. J. L. Rivera, G. K. Jennings and C. McCabe, J. Chem. Phys., 2012, 136, 244701.
46. M. González-Melchor, F. Bresme and J. Alejandre, J. Chem. Phys., 2005, 122, 104710.
47. P. Orea, J. Lopez-Lemus and J. Alejandre, J. Chem. Phys., 2005, 123, 114702.
48. A. W. Adamson and A. P. Gast, Physical Chemistry of Surfaces, Wiley-Interscience, New York, 6th edn, 1997.
49. N. C. Osti, B. Dyatkin, M. W. Thompson, F. Tiet, P. Zhang, S. Dai, M. Tyagi, P. T. Cummings, Y. Gogotsi, D. J. Wesolowski and E. Mamontov, Phys. Rev. Mater., 2017, 1, 35402.

---