LiquidVapor Coexistence in a Primitive Model for a Room

9046

2009, 113, 9046–9049 Published on Web 06/17/2009

Liquid-Vapor Coexistence in a Primitive Model for a Room-Temperature Ionic Liquid Marianela Martı´n-Betancourt, Jose´ M. Romero-Enrique,* and Luis F. Rull Departamento de Fı´sica Ato´mica, Molecular y Nuclear, Area de Fı´sica Teo´rica, UniVersidad de SeVilla, Apartado de Correos 1065, 41080 SeVilla, Spain ReceiVed: April 22, 2009; ReVised Manuscript ReceiVed: June 5, 2009

We present a primitive model for a room-temperature ionic liquid, where the cation is modeled as a charged hard spherocylinder of diameter σ and length l and the anion as a charged hard sphere of diameter σ. Liquid-vapor coexistence curves and critical parameters for this model have been studied by grand-canonical Monte Carlo methods. Our results show a decrease of both the critical temperature and density as the cation length l increases. These results are in qualitative agreement with recent experimental estimates of the critical parameters. Room-temperature ionic liquids (RTILs) are a class of organic salts which present unusual low melting temperatures (usually below 100 °C). In recent years, these liquids have attracted much interest because of their peculiar properties: extremely low volatility, nonflammability, chemical and thermal stability, and interesting solvation and coordination properties.1-4 Due to these characteristics, RTILs have been suggested to be the environmentally friendly replacement of conventional media in chemical reactions. They can be tailored for specific chemical processes, due to the large number of available anions and cations.5 Despite the huge interest on these systems and the increasing amount of experimental results which can be found in the literature, the theoretical understanding of RTILs and their phase diagrams is still incomplete. From a computational point of view, most of the simulation studies of RTIL consider atomistic models,8-10 although simple models have been applied to their liquid crystalline phases11 or dynamical properties.12 This Letter will focus on the vapor-liquid phase transition of some families of RTIL, in particular the imidazolium- and pyridinium-based RTILs. From the experimental point of view, it is widely accepted that these liquids have very low volatility at room temperature. However, recent experimental studies report that some families of RTILs can be successfully distilled without decomposition at very low pressures and high temperatures.6 Previously, the critical temperatures were estimated for some RTILs by extrapolation from available experimental data.7 We will consider a model which is an extension of the simplest model for electrolytes: the so-called restricted primitive model (RPM). This model has been extensively studied by computer simulation,13-18 finding the remarkably low values in reduced units T*c ) 0.0493(3) and F*c ) 0.075;18 see eq 1. Modified versions of this model to account for charge or size asymmetry have also been studied in the literature.19-23 We will consider that the RTIL is formed by N cations, modeled as hard spherocylinders of length l and diameter σ with an embedded point charge +q placed at the center of one of the hemispherical caps, and N anions modeled as charged hard spheres of diameter σ which carry a charge -q at their center. The hard-body * To whom correspondence should be addressed. E-mail: [email protected].

10.1021/jp903709k CCC: $40.75

interactions are supposed additive. The interaction energy between two nonoverlapping ions, i and j, of charges qi and qj () (q) separated by distance rij is Uij ) qiqj/4πε0rij, where ε0 represents the vacuum dielectric constant. The primitive model of the RTIL implicitly assumes that the vapor-liquid phase transition is driven mainly by Coulombic forces, so the dispersive intermolecular interactions can be treated as a perturbation close to the critical point. Additionally, we consider the cationic alkyl chains to be rigid. We will choose σ as the length scale, and the relevant energy scale is defined by the interaction energy between two unlike charges at their closest distance q2/4πε0σ. The reduced temperature and density have the expressions13

T* ) kBT4πε0σ/q2

and

F* ) Fσ3

(1)

where F ) 2N/V is the total ionic number density. As in the case of the primitive models for electrolytes, the low temperatures and densities of interest here (2F*c < 0.2) imply that the clustering is controlled by the mimimum collision distance between unlike charges σ, so eq 1 provides the natural definition for the reduced density. Finally, the reduced simulation box length and spherocylinder length are defined similarly via L* ) L/σ and l* ) l/σ. We have performed neutral grand-canonical Monte Carlo (MC) simulations (characterized by a temperature, T, and a chemical potential for a pair of unlike ions, µ) on cubic boxes of length L, under periodic boundary conditions. The basic MC steps are biased attempts to insert/delete a pair of unlike ions which are accepted/rejected using a modified Metropolis criterion, as described in ref 13. To handle the long-ranged character of the electrostatic interactions, we used the Ewald summation technique with the conducting boundary condition at infinity, 518 Fourier-space wave vectors, and a real-space damping parameter of κ ) 5. With this choice of Ewald parameters, the relative error due to the infinite sums truncation in the electrostatic energy was shown to be less than 10-5 for charged systems at random configurations in small systems.24 Two different sets of simulations have been considered. In  2009 American Chemical Society

Letters

Figure 1. Schematic representation of two unlike ions of the RTIL primitive model (in two dimensions), as considered by the finediscretization approximation. The shaded areas correspond to their hard cores. The different length scales of the system are highlighted.

addition to the usual MC simulations in the continuum (following ref 13), we considered a fine-discretization approach, which is computationally more efficient than the continuum counterparts25 (see Figure 1). In this approach, the positions available to each point charge are the sites of a simple cubic lattice of spacing a. In our simulations, we choose the lattice refinement parameter to be ζ ) σ/a ) 10. In this way, we can precompute the electrostatic interaction between two ions placed at any pair of sites, reducing the computational cost of the simulation with respect to the MC simulations in the continuum. To place the spherocylinder, we allow the center of the noncharged hemispherical cap to be located in any lattice site which is at a distance between l and l + δl from the site where the +q charge is. The value of δl is chosen large enough to have a significant sample of orientations distributed uniformly. Typically, we choose δl/σ ) 0.02 for l* g 0.5, so the number of possible orientations is 134, 144, and 318 for l* ) 0.5, 0.75, and 1.0, respectively. The case l* ) 0.25 presented the problem of nonexistence of available spherocylinder orientations for δl ) 0.02σ. Instead, we considered in the fine-discretization approach the case l* ) 0.2 with δl/σ ) 0.1, so the average spherocylinder length is 0.25, with 96 available orientations. Finally, in order to make the MC sampling more efficient, we choose the bias function for the pair creation to be proportional to the Boltzmann factor associated with the total pair potential between the two inserted unlike ions (including the hard-core contribution).26 Histogram reweighting techniques were used to obtain the vapor-liquid envelopes up to T < 0.98Tc.27,28 Effective critical points for given L* and their extrapolations in the thermodynamic limit were estimated using mixed-field finite-size scaling methods,29 assuming Ising-type criticality, which has been shown to be the universality class for the RPM.17 This approach should be satisfactory to discern a systematic dependence on l*, although the inclusion of the pressure in the field-mixing31 may have an effect on our results, especially in the predicted values of the critical density.32 The phase diagrams from the continuum and fine-discretized simulation schemes were calculated using L* ) 15 data and are shown in Figure 2. We note that both procedures yield quite similar curves, so configurational space discretization does not have a dramatic effect on the results. The critical point parameters obtained from the continuum simulations are listed in Table 1. These results show that there is a marked downward shift in T*c as l* increases. The critical density F*c also decreases with l* but at a lesser rate. As the MC acceptance ratios decrease with larger values

J. Phys. Chem. B, Vol. 113, No. 27, 2009 9047

Figure 2. Coexistence curves of the RTIL primitive model for l* ) 0.0 (circles), l* ) 0.25 (squares), l* ) 0.5 (diamonds), l* ) 0.75 (downside triangles), and l* ) 1.0 (upside triangles). The filled symbols correspond to the results from the MC simulations in the continuum, and the open symbols, from the MC simulations in the fine-discretization approximation.

TABLE 1: Continuum MC Critical Parameters l* 0

0.25

0.50

0.75

1.00

L*

T*c × 102

-µ*c

F*c × 102

-s

12 15 18 ∞ 12 15 18 ∞ 12 15 18 ∞ 12 15 18 18 ∞ 12 15 18 ∞

4.94(2) 4.92(2) 4.93(1) 4.91(2) 4.58(1) 4.54(1) 4.55(1) 4.53(2) 4.21(2) 4.21(3) 4.17(1) 4.17(3) 3.93(3) 3.89(4) 3.87(1) 3.87(1) 3.83(4) 3.57(5) 3.58(4) 3.57(3) 3.57(5)

1.3383(3) 1.3378(1) 1.3379(1)

7.7(1) 7.6(2) 7.7(1) 7.6(1) 7.0(1) 6.9(2) 6.8(2) 6.5(2) 6.5(4) 6.3(2) 6.1(2) 5.6(4) 5.8(3) 5.6(2) 5.7(1) 5.7(1) 5.6(3) 5.7(5) 5.1(1) 5.0(2) 4.0(5)

0.744(4) 0.750(6) 0.744(4)

1.2957(1) 1.2950(6) 1.2952(5) 1.2768(9) 1.2768(9) 1.2764(3) 1.2673(9) 1.2668(9) 1.2666(3) 1.2666(3) 1.2604(9) 1.2605(9) 1.2603(7)

0.745(1) 0.770(8) 0.763(9) 0.759(4) 0.757(5) 0.773(8) 0.754(7) 0.766(3) 0.770(6) 0.770(6) 0.745(9) 0.764(9) 0.764(5)

of l*, we could not obtain reliable values of the critical parameters for l* > 1 due to the increasing sampling difficulties. We argue that the shift toward lower temperatures of the vapor-liquid coexistence has an entropic origin, as the presence of the noncharged tails reduces the number of available energetically favorable ionic configurations. We check this hypothesis by noting that the average cohesive energy at a given density in the liquid phase is roughly independent of the spherocylinder length. The reason for this is twofold: the main contribution to the cohesive energy comes from the interaction between the cation and its associated anion, and on the other hand, the low densities involved in the liquid phase imply that there is still free space to accommodate the ions in allowed energetically favorable configurations. In fact, the reduced cohesive energy per ion pair in the liquid phase near coexistence lies in a range between -1.2 and -1.3 for all of the considered cases. Thus, the reduction of the critical temperature must come from the weakening of the effective interactions between the relevant ionic clusters (mostly dimers) due to the reduction of available energetically favorable configurations. We have performed a cluster analysis of the vapor phase at a corresponding state of T ) 0.95Tc and F ) Fc/3. The fraction of (i,j) clusters composed by i cations and j anions is shown in

9048

J. Phys. Chem. B, Vol. 113, No. 27, 2009

Figure 3. Fraction f(i,j) of (i,j) clusters composed by i cations and i ions in the vapor phase at the corresponding state T ) 0.95Tc and F ) Fc/3 for l* ) 0.0 (circles), l* ) 0.5 (diamonds), and l* ) 1.0 (triangles). The values (i,j) are highlighted for the neutral most populated clusters.

Figure 4. Snapshot of an ionic configuration of the RTIL primitive model for l* ) 1.0 at T* ) 0.0351 and instantaneous density F* ) 0.00233. The blue particles are cations (the noncharged tails are in a lighter color than the charged cap), and the anions are in red. The reduced box size is L* ) 30.

Figure 3 for different values of l*. We have used Gillan’s definition33 based on the position of the charges, with a cutoff distance of Rc ) 1.5σ. First, we note that the ions cluster mainly into neutral aggregates as in the RPM.34,35 For l* e 0.5, the dependence of the cluster fractions on the number of ions which compose them is quite similar, with the dimers being the dominant aggregates of the vapor phase. This observation is in agreement with a recent experimental analysis of the RTIL vapor phase.36 For l* > 0.5, larger clusters are more often observed in the vapor phase, and even for l* ) 1.0, there are more tetramers than dimers. Visual inspection of typical snapshots of the vapor phase for l* ) 1.0 (see Figure 4) shows that both ringlike and more compact clusters are observed, where the noncharged cationic tails orient outward from the cluster. We argue that the cationic tails are excluded from the cluster core to disturb to a lesser extent their energetically favorable conformations. From an experimental point of view, it is not possible to directly get the critical parameters as the RTIL molecules decompose. However, they were estimated from empirical

Letters

Figure 5. Comparison of the MC critical temperatures (filled triangles) with their estimates from ref 7 for the [BF4] (circles), [PF6] (squares), and [Ntf2] (diamonds) families of 1-alkyl-3-methyl-imidazolium-based RTILs. Open symbols are obtained using the Guggenheim equation, and filled symbols, using the Eo¨tvos equation.

correlations with surface tension and density data.7 The reduction of the critical temperature as the aliphatic tail increases in a series of 1-alkyl-3-methyl-imidazolium-based salts (with different anions) is an indication that the dispersive forces are not the sole driving force of the vapor-liquid phase transition. If this were the case, we should expect an enhancement of dispersive forces for longer cations, inducing higher critical temperatures. The critical temperature drop might suggest that as the alkyl chain length is increased the entropic contribution to the cohesive (free) energy diminishes faster than the dispersive contribution to the cohesive energy is enhanced. In order to compare our simulations with these estimations, we used q ) e, the absolute value of the electron charge. The value of σ is obtained for 1-alkyl-3-methyl-imidazolium-based RTILs by correlating the positions of the maxima and minima of the radial distribution functions obtained from our simulations and those corresponding to atomistic models,8,37,38 as we noted that these maxima and minima locations are quite insensitive to the density. From this comparison, we estimate σ ≈ 4 Å. Finally, the dependence on the number of carbons n of the spherocylinder l is estimated from available force fields in the literature, as l ≈ 1.3n Å. In the latter, we assumed that the aliphatic tail of the cations is rigid and completely stretched. Figure 5 shows the comparison between the estimations of the critical temperature from ref 7 (using the Guggenheim and Eo¨tvos empirical equations) and our MC results. Although our MC values are above, they extrapolate for larger chains to the estimations from the experimental data. Note that as the number of carbons increases we anticipate a shallower decay in the critical temperature, as flexibility will play a role in the determination of the effective cation length. In summary, we have performed a MC simulation study of a primitive model of RTILs. We observe a clear decrease of the critical temperature with the length of the cation, in agreement with some experimental estimations. We also observe a higher clustering for large l. Our results confirm that the vapor-liquid transition is mainly driven by the Coulombic forces, and that the decrease of the critical temperature can be understood from the interplay between electrostatics and hard-core interactions. However, asymmetry between the cation and anion radii, dispersive interactions, and the flexibility of the cationic tail may have an influence. Further work is needed to elucidate this point. Acknowledgment. We acknowledge financial support from Ministerio de Ciencia e Innovacio´n (Spain) through grant

Letters ENE2007-68040-C03-02 and Junta de Andalucı´a (Spain) through grant P06-FQM-01869. References and Notes (1) Welton, T. Chem. ReV. 1999, 99, 2071. (2) Seddon, K. R. Chem. Eng. 2002, 730, 33. (3) Rogers, R. D.; Seddon, K. R. Science 2003, 302, 792. (4) Ionic Liquids in Synthesis, 2nd ed.; Wassercheid, P., Welton, T., Eds.; Wiley-VCH: Weinheim, Germany, 2008. (5) Holbrey, J. D.; Seddon, K. R. Clean Prod. Proc. 1999, 1, 223. (6) Earle, M. J.; Esperanca, J. M. S. S.; Gilea, M. A.; Lopes, J. N. C.; Rebelo, L. P. N.; Magee, J. W.; Seddon, K. R.; Widegren, J. A. Nature 2006, 439, 831. (7) Rebelo, L. P. N.; Lopes, J. N. C.; Esperanca, J. M. S. S.; Filipe, E. J. Phys. Chem. B 2005, 109, 6040. (8) Hanke, C. G.; Price, S. L.; Lynden-Bell, R. M. Mol. Phys. 2001, 99, 801. (9) Margulis, C. J.; Stern, H. A.; Berne, B. J. J. Phys. Chem. B 2002, 106, 12017. (10) Morrow, T. I.; Maginn, E. J. J. Phys. Chem. B 2002, 106, 12807. (11) Avendan˜o, C.; Gil-Villegas, A. Mol. Phys. 2006, 104, 1475. (12) Malvaldi, M.; Chiappe, C. J. Phys.: Condens. Matter 2008, 20, 035108. (13) Orkoulas, G.; Panagiotopoulos, A. Z. J. Chem. Phys. 1994, 101, 1452. (14) Caillol, J. M.; Levesque, D.; Weis, J. J. J. Chem. Phys. 1997, 107, 1565. (15) Orkoulas, G.; Panagiotopoulos, A. Z. J. Chem. Phys. 1999, 110, 1581. (16) Yan, Q. L.; de Pablo, J. J. J. Chem. Phys. 1999, 111, 9509. (17) Luijten, E.; Fisher, M. E.; Panagiotopoulos, A. Z. Phys. ReV. Lett. 2002, 88, 185701. (18) Kim, Y. C.; Fisher, M. E. Phys. ReV. Lett. 2004, 92, 185703. (19) Romero-Enrique, J. M.; Orkoulas, G.; Panagiotopoulos, A. Z.; Fisher, M. E. Phys. ReV. Lett. 2000, 85, 4558.

J. Phys. Chem. B, Vol. 113, No. 27, 2009 9049 (20) Romero-Enrique, J. M.; Rull, L. F.; Panagiotopoulos, A. Z. Phys. ReV. E 2002, 66, 041204. (21) Panagiotopoulos, A. Z.; Fisher, M. E. Phys. ReV. Lett. 2002, 88, 045701. (22) Hynninen, A. P.; Dijkstra, M.; Panagiotopoulos, A. Z. J. Chem. Phys. 2005, 123, 084903. (23) Kim, Y. C.; Fisher, M. E.; Panagiotopoulos, A. Z. Phys. ReV. Lett. 2005, 95, 195703. (24) Hummer, G. Chem. Phys. Lett. 1995, 235, 297. (25) Panagiotopoulos, A. Z.; Kumar, S. K. Phys. ReV. Lett. 1999, 83, 2981. (26) Martı´n-Betancourt, M. M.Sc. Thesis, Universidad de Sevilla, 2008. Martı´n-Betancourt, M.; Romero-Enrique, J. M.; Rull, L. F. Work in preparation. (27) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1988, 61, 2635. (28) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1989, 63, 1195. (29) Wilding, N. B.; Bruce, A. D. J. Phys.: Condens. Matter 1992, 4, 3087. (30) Fisher, M. E.; Orkoulas, G. Phys. ReV. Lett. 2000, 85, 696. (31) Kim, Y. C.; Fisher, M. E.; Orkoulas, G. Phys. ReV. E 2003, 67, 061506. (32) Kim, Y. C.; Fisher, M. E. J. Phys. Chem. B 2004, 108, 6750. (33) Gillan, M. J. Mol. Phys. 1983, 49, 421. (34) Bresme, F.; Lomba, E.; Weis, J. J.; Abascal, J. L. F. Phys. ReV. E 1995, 51, 289. (35) Caillol, J. M.; Weis, J. J. J. Chem. Phys. 1995, 102, 7610. (36) Leal, J. P.; Esperanca, J. M. S. S.; da Piedade, M. E. M.; Lopes, J. N. C.; Rebelo, L. P. N.; Seddon, K. R. J. Phys. Chem. A 2007, 111, 6176. (37) Shah, J. K.; Brennecke, J. F.; Maginn, E. J. Green Chem. 2002, 4, 112. (38) del Po´polo, M. G.; Voth, G. A. J. Phys. Chem. B 2004, 108, 1744.

JP903709K