Surface of Half-Neutralized Diamine Triflate Ionic Liquids. A Molecular Dynamics Study of Structure, Thermodynamics, and Kinetics of Water Absorption and Evaporation

Surface properties of room temperature ionic liquids (RTILs) consisting of half neutralized diamine cations (H2N−(CH2)n−NH3, n = 2, 4) and triflate anions have been investigated by molecular dynamics simulations, based on an empirical atomistic force field. Planar slabs periodically repeated in 2D have been considered, and the temperature range 260 ≤ T ≤ 360 K has been covered, extending from below the melting and glass point to the equilibrium liquid range of the diamine compounds under investigation. Addition of water at 1% weight concentration allowed us to investigate the kinetics of water absorption through the RTIL surface, and to characterize the structural and dynamical properties of subsurface water. Animations of the simulation trajectory highlight the quick absorption of water molecules, progressing downhill in free energy and taking place without apparent intermediate kinetic stages. To verify and quantify these observations, a variant of the umbrella sampling algorithm has been applied to compute the variation of excess free energy upon displacing a water molecule along the normal to the surface, from the center of the slab to the vapor phase. The results provide a comprehensive picture of the thermodynamic properties underlying the kinetics of water absorption and evaporation through the surface, and they also provide the ratio of the equilibrium density of water in the vapor and liquid phase at the average concentration considered by simulations. A variety of properties such as the surface energy, the 90−10% width of the profile, the layering of different species at the interface, and the electrostatic double layer at the surface are computed and discussed, focusing on the effect of water contamination on all of them.


I. INTRODUCTION
Organic ionic compounds of low melting temperature (T m ≤ 373 K) known as room-temperature ionic liquids 1,2 (RTILs) summarize in themselves properties such as the low volatility and high thermal and chemical stability of ionic systems, together with the wide range of sizes, modular structure, and low melting point of organic molecules. 3,4 At least in their most popular manifestations, moreover, RTILs present low toxicity and moderate environmental impact, 5 and their selective affinity for biomolecules 6 opens the way to future applications in pharmacology 7 and biomedicine. 8 Last but certainly not least, their ionic character makes RTILs able to conduct an ionic current, which is an additional property that could foster applications. In most cases, however, the high viscosity and correspondingly low diffusion constants of RTILs, 9 together with a fairly high degree of ion association, 10 limit their electric conductivity to relatively low values.
Protic ionic liquids 11−13 (PILs) represent a special class of RTILs that can be synthesized by neutralizing a Brønsted acid with a Lewis base. PILs, in general, share the positive properties as well as the limitations of other classes of ionic liquids, but the acidic proton that gives them their name sometimes displays enhanced mobility, due to proton-specific mechanisms such as Grotthuss conductivity. 14 The improved ion conductivity, in turn, makes selected PILs suitable for applications in electrochemistry, 15 favored also by the wide electrochemical window of these compounds.
These considerations together with promising experimental results 16 motivate an intensive effort devoted to the investigation of protic ionic liquids, especially in view of applications in power generation, energy conversion, and storage devices such as batteries, fuel cell, accumulators, and supercapacitors. 17 Performance and efficiency considerations often impose relatively high operational temperatures for these devices, usually exceeding room temperature and reaching up to 400 K.
Ionic liquid tend to be hygroscopic, 18 and this is particularly true of protic ionic liquids, since water molecules find ample opportunities to bind and diffuse in their hydrogen bond network. Contamination by water enhances the ions and proton mobility and thus electric conductivity, but water turns into vapor at high temperature, giving origin to high pressures that might endanger the device integrity.
To investigate all these issues, a paradigmatic role could be attributed to a family of PIL compounds made of aliphatic diamine cations and the trifluoromethanesulfonic (triflate) anion, whose chemical physics properties, and proton conductivity in particular, have been investigated by experiments 16 and by simulation. 19 More precisely, cations in these compounds consist of a conjugated pair of amino (−NH 3 + )/ amine (−NH 2 ) terminations, joined by an aliphatic segment [−(CH 2 ) n −], and their charge is compensated by an equal number of triflate (CF 3 SO 3 − , [Tf] − ) anions. In our study we consider the case n = 2, ethylenediamine triflate, [DAEt] [Tf], and, to a lesser extent, the case n = 4, butyldiamine triflate, [DABu] [Tf]. The experimental measurements of ref 16 show that [DAEt][Tf] melts at T m = 351 K, and it is chemically stable up to T = 500 K. Corresponding data for [DABu] [Tf] are not reported in the same experimental paper, but they are expected to be similar. Molecular dynamics simulations 19 show a glass transition on cooling at T g = 310 K in [DAEt][Tf] and document the formation and role of mesophases, due to the interplay of Coulomb forces, packing effects, and hydrogen bonding. NMR measurements of the diffusion coefficients in [DAEt] [Tf] reveal that protons diffuse faster than cations, implying that vehicular proton transport is supplemented by proton jumps among cations, corresponding to the well-known Grotthuss proton transport in water. Both experiments and simulations show that, even at low concentration, water enhances the mobility of all species, i.e., anions, whole cations, and single protons. As already mentioned, however, water contamination gives raise to other, potentially more serious problems.
In the present study, by molecular dynamics simulations we investigate water absorption and desorption through the free planar surface of liquid [DAEt][Tf] and [DABu] [Tf] samples, aiming at quantifying the thermodynamic and kinetic parameters underlying these properties. Simulations have been carried out for temperatures 280 K ≤ T ≤ 360 K, closely corresponding to the temperature range of liquid water and covering both the experimental melting temperature of [DAEt][Tf], at T m = 351 K, and its glass transition that, according to our model, occurs at T g = 310 K.
To provide a term of comparison, we report also results for the surface properties of dry samples. In this case, simulations predict a surface energy of about 110 mN/m over the entire 280 K ≤ T ≤ 360 K range, revealing a relatively sharp surface, whose 90−10% width (see the definition below) is comparable with or lower than the ion sizes. Since the system temperature is relatively low, close or below its melting and glass temperatures, the density profile oscillates in proximity of the surface, highlighting a fair degree of ordering at the surface. The amplitude of the density oscillations, however, is modest and decays quickly moving inside the slab at all temperatures T ≥ 300 K. In contrast, the surface layering of the constitutive chemical groups (−CF 3 , −SO 3 , −NH 3 , −NH 2 , −CH 2 −) is much more marked, driving a separation of charge that causes the electrification of the surface, with a drop in the electrostatic potential from the liquid to the vacuum of ∼0.25 V.
Water addition at 1% weight concentration changes only quantitatively these properties, leaving unchanged the qualitative picture. Analysis of simulation trajectories shows that the incorporation of water into the slab is a fast process, in which the crossing of the first ionic layers takes place within a few ns time scale. The estimate of the water density in the vapor phase and well within the slab shows that water molecules inside enjoy a gain in the excess free energy of about 25 kJ/mol with respect to those far from the slab. The computation of the free energy profile experienced by a water molecule across the liquid−vapor interface confirms that no intermediate state can be identified along the absorption− desorption process.
Similar computations for the slightly larger diamine species [DABu] + (H 2 N−(CH 2 ) 4 −NH 3 + ) provides the same qualitative picture, with only quantitative differences in most properties.
This and similar studies could provide the insight required to optimize protic liquids with respect to water sorption, perhaps modifying their surface layer by a third species that could represent a barrier to water absorption and desorption into the RTIL.
More in general, the present study contributes to the accumulated knowledge of the surface of dry and water contaminated molecular ionic liquids. Over the years, these properties have been investigated many times by experiments and by simulations. Experiments on nominally dry samples, in particular, determined the surface tension, surface energy, and surface entropy as a function of temperature for broad RTIL families, in most cases involving compounds based on the alkyl-imidazolium cation. 20−23 The coverage of structural properties is less exhaustive, partly because of the wide variety of RTIL compounds. The effect of water contamination on the RTIL surface has not been studied in much detail, but a relevant exception is represented by the study of ref 24, investigating the effect of water on hydrophilic and hydrophobic RTILs and concluding that the effect is larger for hydrophobic salts.
Molecular dynamics has been used, again in the case of dry RTILs, to investigate primarily structural and dynamical properties at surfaces (see, for instance, refs 25−28 and references therein), while thermodynamic properties beyond the surface energy have been studied less extensively. Needless to say, the distinction between experimental and computational studies is only a matter of (incomplete) classification, and many studies of RTIL surfaces involve both experimental and simulation approaches. 29 Once again, the surface of wet RTILs has been sparingly investigated. The simulation study closest to the present one might be represented by ref 30, which, however, considered primarily systems of high water concentration.
Given the vast number of RTIL compounds, the studies reported until now cover a tiny fraction of the whole RTIL chemical space, even before considering the additional concentration variable of water contaminated systems. For this reason, the present investigation of compounds chemically and structurally different from the popular alkyl imidazolium salts might contribute to a better view of the variety of properties and trends at the surface of dry and wet RTIL samples. The Journal of Physical Chemistry B Article using an empirical, atomistic force field model described in detail in ref 19. The functional form corresponds to the popular OPLS 31 /Amber 32 model, representing molecules in terms of particles and covalent bonds, expressing the potential energy as a function of atomic coordinates by the sum of bonded U b and nonbonded U nb contributions:

II. METHODS
where N is the number of particles (atoms, in our case). U b accounts for stretching, bending, and torsion energies: (2) with: (5) where i, j k, and l are atoms joined by 1, 2, and 3 consecutive bonds, K s , K b , and K t are force constants and r ̅ , θ̅ , and ϕ̅ are reference values for bond stretching, bending, and torsion, respectively. The geometry of covalent bonds in the cations and anions considered in our simulations is schematically shown in Figure 1.
Nonbonded terms account for dispersion interactions, described by Lennard-Jones pair potentials, and electrostatic contributions. They are given by with ϵ ij and σ ij satisfying the Berthelot's rule: ii jj (7) ϵ = ϵ ϵ ij ii jj (8) The charges {q i } and the {σ ii }, {ϵ ii }, therefore, are the only independent parameters representing the nonbonded potential. The prime on the sum indicates that nonbonded terms between atoms separated by up to two consecutive bonds vanish, while those between atoms separated by three consecutive bonds are divided by two.
The many parameters defining both bonded and nonbonded interactions are listed in the Supporting Information section of ref 19. Water molecules, in particular, are modeled through the TIP4P potential. 33 To represent slabs of finite thickness, the sample is enclosed in an orthorhombic simulation box of fixed shape and volume, with periodicity L z = 100 Å along the direction normal to the slab surface approximately three times the slab thickness.
The short-range part of Coulomb and dispersion energies are computed within a cutoff of 10 Å. Long range Coulomb interactions are accounted for using the smooth particle mesh Ewald (SPME) algorithm. 34 Long range corrections to dispersion are applied to the computed energy, according to the Gromacs implementation. Since the implemented approach assumes the system to be homogeneous, the correction is likely to be somewhat overestimated at the free surface and underestimated in the condensed phase. The value of the long-range dispersion correction computed by Gromacs, however, is on the order of the percent of the cohesive energy, and the slight distortion is unimportant.
Constant temperature conditions have been enforced using the Nose−Hoover thermostat, with a relaxation time constant of 2 ps. The free surface of the slab and the zero or negligible density in the vapor phase imply that the pressure on the simulated sample virtually vanishes. Therefore, simulations have been carried out in the NVT ensemble, using the Gromacs package 35 with an integration time step of 1 fs.
A crucial part of our study is the determination of the free energy profile for the motion of water molecules along the direction z perpendicular to the free surface. To this aim, the water molecule closest to the center of mass (COM) of the slab has been selected and pulled along the z direction by steered MD. Intermediate positions separated by 0.1 nm from each other provide 31 initial configurations in which the COM position of the selected water molecule has been restrained using a harmonic force constant of 1000 kJ/mol/nm 2 . After equilibration, restrained simulations are used to collect umbrella sampling statistics. The free energy profile has been obtained by using the WHAM method as implemented in Gromacs 5.1.2. 36 The method is known to have drawbacks and limitations. 37 However, the consistency of the free energy profiles among themselves and with other data obtained from the simulations supports the validity and reliability of the method in the present case.
II.A. Sample Preparation and Simulation Setup. Slabs consisting of 216 [DAEt][Tf] molecules have been prepared starting from a homogeneous liquid sample equilibrated at T = 440 K (see ref 19), increasing the periodicity L z along z in discontinuous steps, while keeping constant the periodicity along xy. At first the change of periodicity has to be slow, and each increase in L z (by ∼1 Å) is followed by minimization of the energy at constant volume, and by brief (50 ps) The Journal of Physical Chemistry B Article equilibration at NVT conditions back to 440 K. Soon a density gap forms, separating two adjacent surfaces, and it becomes possible to extend L z by large (∼10 Å) steps without perturbing much the system. The final result is a liquid slab of thickness ∼36 Å, exposing two (on average planar) surfaces of 37.81 × 37.81 Å 2 cross section, in an orthorhombic cell periodically repeated in 3D, of size 37.81 × 37.81 × 100 Å 3 . During all simulations, the periodicity along xy is fixed at the size of the original homogeneous and cubic sample at T = 440 K; hence, the width of the slab shrinks slightly with decreasing T to account for the thermal expansion of the sample.
Although we did not estimate this quantity, we assume that the interaction of globally neutral surfaces separated by more than 60 Å of empty space is negligible. The leading correction to energy and forces is determined by the instantaneous dipole moment of the whole simulation cell. 38 Since, on average, the slab is symmetric, this correction term depends on the fluctuations of the dipole, whose average value vanishes. These same fluctuations are investigated below to extract the frequency dependent dielectric function of the slab, and this is found to be small. A more relevant problem is the residual interaction of the two surfaces that occurs through the condensed phase, due to the limited thickness of the slab. The presence of this spurious effect is discussed in the Results.
As stated in the Introduction, our investigation focuses on the 280 K ≤ T ≤ 360 K temperature range. Therefore, the first [DAEt][Tf] sample at T = 440 K has been annealed down to T = 360 K during 100 ns. The temperature interval of interest has been covered by decreasing the temperature in discontinuous steps of 10 K, each time equilibrating the sample during 60 ns. Statistics have been accumulated every 20 K, i.e., at T = 360, 340, 320, 300, and 280 K, each time running at NVT conditions during 200 ns.
A simulation snapshot of the [DAEt][Tf] slab at T = 300 is shown in Figure 2. sample at the same temperature and vanishing pressure allows us to compute the surface energy u(T), reported in Table 1 for selected temperatures. Somewhat unexpectedly, the surface energy decreases, although slowly, with increasing T. The analytical interpolation (not explicitly carried out) of u(T) by a low order polynomial or Padéfunction would allow the estimation of the temperature dependence of the surface entropy s, according to where s is expressed in units of k B and all quantities refer to the same P = 0 condition. The surface energy decreases slightly with increasing T and du(T)/dT < 0; hence, the surface entropy decreases with increasing T. Since the surface entropy is the difference of entropy between the sample with and without a surface (or interface), there is no strict constraint on the sign of ds(T)/dT or even of s(T) itself. Moreover, the temperature dependence of u is apparently moderate; hence, the role of entropy is modest, and in this case the surface energy could be considered as a good approximation to the surface tension γ. A semiempirical approach to compute the surface tension and entropy from u(T) is discussed in the Supporting Information. The low value of the surface entropy agrees with the findings of previous studies on different families of ionic liquids. 23 III.A.2. Structural Properties. The first structural quantity we characterize is the mass density distribution across the slab, reported in Figure 3 at two different temperatures, i.e., T = 360 K, just above the experimental melting temperature T m = 351 K, and T = 300 K, with the systems in the glassy state. More precisely, Figure 3 reports the density divided by the average density of the homogeneous sample at the same temperature, reported in ref 19. Statistics could be improved by averaging with respect to the plane z = 0 that marks the center of the slab, but we decided to show the raw profile since the slight left−right asymmetry gives a measure of the statistical uncertainty. At both temperatures the profile shows a relatively sharp interface. Density oscillations of low amplitude are apparent in the condensed phase, of wavelength comparable to the size of the [DAEt][Tf] molecules. As often observed in computer simulations of liquid surfaces and interfaces, the amplitude of this oscillations decreases nearly exponentially with increasing depth from the surface. Because of the relatively thin width of the slab, the density (and charge) The Journal of Physical Chemistry B Article oscillations arising from two surfaces are not completely decoupled. As expected, sharpness and amplitude of the oscillations increase with decreasing temperature, the characteristic length of the oscillations remains nearly constant, and the penetration of oscillations inside the slab is enhanced at low T. Although these trends are clear from the results, the temperature dependence of the profile is mild from T = 360 K down to T = 300 K, consistently with the highly viscous state of the sample, becoming glassy (in our model) below T g = 310 K.
The picture changes somewhat below T = 300 K, since the density profile undergoes a nearly discontinuous change, with layering suddenly becoming much sharper at T = 280 K. Moreover, the amplitude of the density oscillations hardly decays inside the slab, pointing to a partially crystallized samples, or, more likely, to a nanocrystalline state. The comparison of the profiles at T = 300 K and at T = 280 K is shown in Figure S1 in the Supporting Information.
Because of their strong, ionic cohesion, RTILs are characterized by very low vapor pressure. Not surprisingly, we never observed ions or neutral ion pairs far from the condensed slab, unambiguously residing in the vapor phase.
The sharpness of the interface is quantified by computing the 90−10% separation, measuring the depth over which the density changes from 90% to 10% of the bulk density at the same conditions. Values at selected temperatures are reported in Table 2. The sizable increase of the error bar at T ∼ 300− 320 K is due to the dynamical coexistence of different surface structures that we have been unable to fully resolve, resulting in the broadening of the average density profile.
To the best of our knowledge, no accurate and accepted theory is available to predict the surface width of molecular liquids, and of organic ionic liquids in particular. For this reason, we are unable to interpret the trends that could be seen in the data of Table 2. The wide variety of RTILs, however, could provide an ideal testing ground for the development of such a theory.
The mass density distribution can be analyzed in terms of the composition considering distinct chemical species such as the F, O, N, C atoms, as indicated in Figure 4. The figure shows a layering much stronger than for the total density, in many ways representing a planar version of the mesophase formation well-known from bulk phases. Not surprisingly, the globally neutral and relatively inert −CF 3 group is exposed to the vacuum, apparently decreasing the surface tension. The enrichment of the surface layer in −CF 3 is apparent in the snapshot show in Figure 2. Then, the charged moieties represented by the O and N atoms form two adjacent layers apparently to optimize the system Coulomb energy. Inside the slab, the distributions of neutral and charged domains are closely superimposed. These features, present at all temperatures, become more apparent with decreasing temperature. The layering and the surface ordering of different species might explain the low (even negative) surface entropy of these samples. Below the first two surface layers, the radial distribution functions of all the species are close to those computed for the bulk case, 19 and they display a similar pattern of mesophases. III.A.3. The Electrostatic Double Layer at the Surface. The charge density profile ρ Q (z) is shown in panel a of Figure 5. To improve the signal-to-noise ratio, the charge density profile has [Tf] has been estimated by a semi-empirical van der Waals−Guggenheim approach, 39 using eq S1 of Supporting Information. Parentheses emphasise that the γ values reported in this table are an estimate of unknown uncertainty because of the tentative choice of the temperature T c (see Supporting Information).  The Journal of Physical Chemistry B Article been symmetrized with respect to the z = 0 plane, and the resolution has been slightly decreased with respect to the mass density profile. Despite the oscillations in the first few surface layers, the charge density is low everywhere, being less than 1 e/nm 3 , perhaps reflecting the close association of cations and anions in the system. Since the net average charge is the result of the partial compensation of a large number of atomic contributions, it is difficult to unambiguously assign each peak and valley of the charge density distribution to a unique atomic or ionic species. The corresponding electrostatic potential profile is computed according to where the value of the periodic potential has been set to zero at the simulation cell boundary: The result for ϕ(z) is shown in panel b of Figure 5. An important parameter that characterizes the charge separation at electrified interfaces is the electrostatic potential difference between the condensed and the vapor phases, that in our case is represented by Δϕ = ϕ(0) − ϕ(±L z /2) = ϕ(0). The Δϕ value for dry [DAEt][Tf] slabs is given in Table 3. The limited width of the slab results in oscillations in the electrostatic potential that extend up to the center of the slab. Therefore, the value of ϕ(0) used to compute the Δϕ reported in the table is the average of ϕ(z) over the interval −2.5 ≤ z ≤ 2.5 Å.
The computed Δϕ value is relatively small (see section IV), it increases slightly with decreasing T, and the rate of change with temperature decreases significantly below T = 340 K, probably reflecting the near vitrification of the sample.
The electrostatic potential drop Δϕ is directly related to the dipole moment per unit surface of the surface, according to At T = 300 K, the dipole moments turns out to be 0.70 D per nm 2 , with a sign that corresponds to a dipole vector pointing outward from the slab. The average number of H-bonds in the sample is significantly higher than in the homogeneous sample of the same size at the same temperature. The quantitative enhancement (50%) measured in the simulation reflects again the small size of the sample, since the relative change has to tend to zero with increasing thickness of the slab. Nevertheless, the sizable increase in the HB density at the surface could be revealed in thin films, or, more in general, measured by surface sensitive spectroscopy approaches. More quantitative details on this interesting feature are reported in the following section. In addition to anion−cation HBs, nearly every cation presents a distorted intraion HB. 19  The motion of ions is necessarily anisotropic (see Figure S2 in the Supporting Information) and only in the xy plane parallel to the surface one can observe genuine diffusion. Nevertheless, over the temperature range 280 K ≤ T ≤ 360 K, and over the ∼10 2 ns time scale of the simulations, the mean square displacement along z is one order of magnitude smaller than the square width of the slab, and the saturation in ⟨r 2 (t)⟩ due to finite size is not apparent yet.
While finite size has limited impact, layering of cations and anions is more important. Separate diffusion constants D xy and D z have been computed by dividing the time derivative of the mean square displacements along xy and along z, respectively, by twice the number of the corresponding dimensions. The results show that, on average, the diffusion coefficient along z is one-half of the value in the xy plane. Moreover, the value D xy along xy is enhanced and the value D z along z is reduced with respect to the bulk diffusion coefficient D reported in ref 19. To first approximation, the weighted average (2D xy + D z )/3 is close to the homogeneous bulk value D. Together with D z = D xy /2, this observation implies that in our samples , possibly reflecting the energy of the same −CF 3 layer exposed to the vacuum in the two cases. The density profile, shown in Figure S3 19 Hence, the somewhat reduced layering of the elements in[DABu] [Tf] implies that the analogy between mesophases formation and surface layering is not perfect, or, perhaps, it emphasizes that the balance of molecular packing, dispersion and Coulomb forces is not the same in the bulk and at surfaces.
As before, no single ion or ion pair is found at T = 360 K or T = 300 K in the vapor phase of [DABu][Tf] during the entire simulation.
As expected, the charge density distribution has lower amplitude in [DABu][Tf] than in [DAEt][Tf], as can be seen in Figure 6 comparing the charge density profile in the two cases. As a result, the size of |Δϕ| in [DABu][Tf] is reduced by ∼40% with respect to [DAEt][Tf], being |Δϕ| = 0.141 V at T = 300 K and |Δϕ| = 0.147 V at T = 360 K. Despite the   Before analyzing the equilibrium properties, the protocol to add water is briefly described, and the limited information that MD provides on the kinetics of this process is reported.
The sample was prepared by distributing 12 water molecules at rest in proximity of each surface of an equilibrated dry slab. The in-plane positions have been chosen at random with uniform probability, and the z coordinate has been selected in such a way that each water molecule is separated by ∼5 Å from the closest ion. In all cases, this last turns out to be an anion, since this species, or, more precisely, its −CF 3 moiety, is exposed to the vacuum side of the slab. At T = 360 K, it is virtually impossible to follow the time evolution of absorption, since water molecules cross the vacuum/[DAEt][Tf] interface without any delay, apparently progressing downhill on the potential and free energy surface. Each molecule reaches a subsurface position within a few picoseconds. The water distribution across the slab becomes stationary within 50−100 ps, but this short time obviously depends on the finite and small width of the slab.
The same addition procedure at T = 300 K takes longer, as expected. All water molecules reach a subsurface position within a nanosecond, and the time taken by water to diffuse across the entire slab is correspondingly longer. Also in this case, however, it is not possible to identify positions of metastable equilibrium, or to define stages in the absorption process. The analysis of animations of the simulation trajectories suggests that the important aspect in the downhill absorption and in the subsequent water diffusion across the slab is the close association of water with the [Tf] − anion, with a water−[TF] − average dissociation time within the slab of 150 ps at T = 360 K, growing to 1.9 ns at T = 300 K.
To ease a systematic comparison, the discussion of equilibrium properties follows the same scheme of section III.A on dry slabs. Once again, the surface energy as a function of temperature is computed as the difference of the average potential energy of the slab and of the homogeneous sample at the same number of [DAEt][Tf] and water molecules. The results are reported in Table 1. The surface energy of the wet slab is only slightly higher than that of dry slabs, and also increases slowly with decreasing temperature. These trends are consistent with the structural properties discussed below, and especially with the observation that the outermost layers are virtually devoid of water. Also in this case, however, the moderate temperature dependence implies that the surface energy is a fair approximation to the surface tension.
Since, by design, the total amount of water corresponds to about 1% weight of the entire sample, the mass density profile (see red line in Figure 7) is almost unchanged with respect to the dry slab case for 300 K ≤ T ≤ 360 K. The similarity of dry and wet slabs extends to the probability distribution of F, O, N, and C across the slab, as shown in Figure S5, Supporting Information. At T = 360 K the water density distribution (blue line in Figure 7) is nearly homogeneous over a width that on each side is ∼5 Å narrower than the total density distribution. Comparison with the density distribution of F, O, N, and C, shows that water avoids the surface layer occupies primarily by − CF 3 , while it covers the same range of the O (from the anion) and N (from cation) distributions. A pair of broad peaks in the water probability distribution correspond to peaks in the O (from [Tf] − ) distribution, apparently because of the close water-[Tf] − association. Since the −CF 3 layer exposed to the vacuum is nearly unchanged with respect to the dry case, the close agreement of surface energy between dry and wet slabs can easily be rationalized.
Also at T = 300 K, the mass density profile is nearly the same as in the dry case, the water profile, as expected, is more structured than at T = 360 K, and slightly asymmetric because of the glassy state of the sample. The T = 280 K simulation gives a somewhat different result from the dry-slab case. While the highly modulated density profile of the dry sample at T = 280 K qualitatively differs from the T = 300 K profile (See Figure S1 in Supporting Information), the temperature dependence of the total mass density is more continuous and mild in the case of water-contaminated samples. This could be due to the enhanced mobility of the wet samples, making them more fluid-like, 19 but it might also be the result of metastability and long relaxation times, preventing the system from reaching a nanocrystalline state. On the other hand, the water profile in [DAEt][Tf] at T = 280 K is much more structured than at T = 300 K. A detailed account of the temperature dependence of the water profile is given in Figure  S6 of the Supporting Information. Water molecules can be found in the portion of the simulation box that we can identify with the vapor phase. The average number of water molecules found beyond 5 Å from the average surface location (defined at 1% density) is low, hence the estimate of the vapor density ρ v ≡ ρ H 2 O (v) is affected by a large relative error bar. Nevertheless, the value of ρ H 2 O (v) as a function of T reported in Table 4 provides a first estimate of the excess free energy difference of water molecules within the slab and in the gas phase. To this aim, we express the chemical potential μ H 2 O of water molecules as From the equality of the chemical potential across the system we obtain 13) giving the values reported in Table 4. The assumption, supported by the plot in Figure 7, is that it is possible to define a region of nearly uniform water density in each of the two phases. This preliminary analysis is confirmed and extended in section III.C.
The water vapor density at T = 300 K reported in Table 4 corresponds to a pressure P = 3 × 10 −4 bar, which is right in the middle of the medium-vacuum range (from 3 × 10 −2 to 10 −6 bar), implying that it is not easy to decrease the water content of [DAEt][Tf] much below 1% by applying vacuum to the vapor side through a single stage pump.
The kinetics of water absorption and desorption through the surface of [DAEt][Tf] can be characterized by the following considerations. First, at equilibrium, the rate of water desorption equals the rate of absorption. In the dilute vapor, the motion of water molecules is virtually ballistic, and the flow of water toward the slab through a plane 5 Å beyond the surface (conventionally located at the coordinate of 1% bulk density) can be considered as the rate of absorption from the vapor. Analysis of trajectories shows that the rate of water crossings in each direction through this conventional surface is J = 0.075 ± 0.002 molecules/(nm 2 ns) at T = 360 K, and J = 0.0031 ± 0.0005 molecules/(nm 2 ns) at T = 300 K. The quantitative values, of course depend not only on T, but also on the chemical potential of water in the system. Hence, these two single values in themselves are not very relevant. Nevertheless, they are reported here to show how the combination of equilibrium MD, giving the ration ρ v /ρ l , free energy computations (discussed in section III.C), and analysis of trajectories provide a complete characterization of the kinetics of absorption and desorption over a wide range of thermodynamic conditions. III.B.2. The Electrostatic Double Layer at the Surface and Dielectric Properties of the Slab. The total charge distribution in wet [DAEt][Tf] slabs is shown in panel a of Figure 8, where it is also compared to the contribution arising from the position-dependent orientation of water molecules. The total electrostatic potential across the slab is shown in panel b of Figure 8, together with the contribution due to water, computed again from the charge distribution according to eq 10. The electrostatic potential profile closely follows that of the dry slabs, but the 1% weight contamination by water reduces Δϕ by nearly 10% at all temperatures but T = 280 K, as shown in Table 3. It is tempting to relate this observation to the dielectric properties of bulk water, but the results for [DABu][Tf] discussed below show that other factors enter the picture, since the water orientation not only depends on the local electric field but also is affected by the network of hydrogen bonds.
By definition ionic liquids are conductors. In the two directions parallel to the planar surface, dipole fluctuations and dielectric screening diverge in the static limit, and are difficult to determine also at low nonzero frequency. In the direction normal to a finite slab, however, dipole fluctuations are finite,  The time autocorrelation function for the z-component of the dipole moment M z (t) normalized to ⟨M z (t 0 )M Z (t 0 )⟩ t 0 , given by is shown in Figure 9 for the wet [DAEt][Tf] slab at T = 360 K. The autocorrelation ψ(t) has a purely relaxation character, without any oscillating component that could point to a weakly damped plasmon mode. At least one fast and one slow relaxation mode are apparent, and no fit with a single exponential can faithfully reproduce the simulation result. The superposition of two or three exponential terms provides a good fit to ψ(t), as can be seen in Figure 9 for the fit with three exponentials. The relaxation times determined by the fit (τ 1 = 0.361 ps, τ 2 = 16.85 ps, and τ 3 = 297.1 ps) span a broad range. No attempt has been made to identify the type of motion corresponding to each relaxation channel.
The fit by purely exponential relaxation terms implies that the frequency dependent dielectric function is the superposition of Debye dielectric relaxation terms: 40 j j j (15) where the {A i } are weight factors. Equivalent information is given by the direct computation of the dielectric function through the relation 41,42 (Fourier−Laplace transform): 16) where ϵ S is the static dielectric function. The result for the real and imaginary part are given in Figure 10.
A fit of comparable square deviation and marginally better χ 2 is obtained by writing ψ(t) as the superposition of a stretched exponential of the Kohlrausch−Williams−Watts 43 (KWW) form and a simple exponential relaxation: Since the two fits are virtually equivalent and both are very close to the simulation curve, the corresponding dielectric functions are nearly the same. The inclusion of the stretched exponential, however, points to a variety of relaxation processes beyond the single activated mechanism. Using two stretched exponentials instead of the expression in eq 17 does not improve the fit.
The equilibrium value of the square dipole fluctuation turns out to be ⟨M z (t 0 )M Z (t 0 )⟩ t 0 = 5.456 D/nm 2 , again for the wet [DAEt][Tf] slab at T = 360 K. This value enters the definition of the dielectric constant for the slab along the z coordinate (or the ϵ zz component of the static dielectric tensor), according to but the association of this single number to the inhomogeneous system represented by the slab is somewhat incomplete and ambiguous, as is apparent from the dependence of ϵ zz on the poorly defined volume of the slab. A similar analysis of dipole fluctuations and dielectric properties has been carried out also for the dry slabs. Since the two systems differ only by a slight water contamination, the  The Journal of Physical Chemistry B Article results are qualitatively and quantitatively similar. The relaxation times of the dipole autocorrelation function by three exponentials gives τ 1 = 0.380 ps, τ 2 = 18.22 ps, and τ 3 = 488.1 ps for the dry system. The equilibrium value of ⟨M z (t 0 ) M Z (t 0 )⟩ t 0 turns out to be 5.057 D/nm 2 in the dry [DAEt][Tf] slab, only slightly different from the value of the wet slab. Comparison of the results for the dry and for the wet slab shows that low concentration water contamination systematically decrease relaxation times, implying that water, while increasing the system fluidity and increasing diffusion and conductivity (see ref 19), it also increases the damping of modes, perhaps because it adds further relaxation channels.
III.B.3. Hydrogen Bonding and Diffusion Constants. As already found in the previous section for dry slabs, the number of hydrogen bonds in the wet slab is significantly higher than in the homogeneous sample of the same size, water content, and thermodynamic conditions. The quantitative enhancement with respect to the homogeneous case is particularly important (50%) for HBs connecting ions, and relatively less important (15−20%) for HBs involving water molecules, either among themselves or with ions. The difference might simply reflect the fact that the H-bonding capability of water molecules in [DAEt][Tf] is nearly saturated already in homogeneous samples. Because of the enhancement due to confinement, in fact, almost all HB sites are saturated both in the dry and in the wet slabs. In wet slabs, in particular, 95% of the −NH 3 + groups share two protons, in most cases with two different anions, their third proton being taken by a distorted intracation HB. 19 Moreover, 66% of the water oxygen atoms accept two protons from [DAEt] + , and virtually all water protons are shared with two (again different) [Tf] − . Water−water hydrogen bonding is rare at this low water concentration. As a result of these interactions, the entire slab represents a uniquely connected network of HBs. Somewhat surprisingly, the enhanced stability of the hydrogen bond network is not accompanied by a systematic decrease of the diffusion constant of cations and anions along the plane parallel to the surface.
The enhanced stability of hydrogen bonds seems to be due to geometrical ordering of ions at the surface, and, in particular, within the two packed surface layers enriched in −SO 3 and −NH 3 groups. This picture is supported by the plot of the spatial distribution of hydrogen bonds, shown in Figure  11, where the position of each HB is identified with the coordinate of the binding proton, and the separate populations of HB among ions and involving water molecules have been taken into account. The plot suggests also another reason for the lesser effect on H-bonding of water, since this species is less likely to reside in the highly ordered surface layers.
Analysis of the HB orientation, defined by the direction of the vector joining the N to the O atoms connected by a HB shows that in the outermost layers HB are highly oriented, laying in the plane parallel to the surface (see Figure S7, Supporting Information).
The addition of 1% water to the [DAEt][Tf] slab has the same qualitative effect on the diffusion of ions as for the homogeneous samples, increasing the cation and anion diffusion coefficient by up to 30%. The relative change is larger at low T than at high T, suggesting that the addition of water decreases the free energy barrier to diffusion. On the other hand, water addition does not qualitatively alter the effect of confinement along z on the diffusion of ions. Also in wet slabs, diffusion along xy is faster than in the homogeneous case, while the opposite is true of diffusion along z. As before, and similar to the bulk, cations diffuse faster than anions.
The effect of the slab geometry on the diffusion of water in [DAEt][Tf] is difficult to assess, since the few water molecules in the vapor phase dominate the mean square displacement of this volatile species, preventing a comparison with the homogeneous case. In other words, the continuous exchange between the liquid and the condensed phase hampers the independent evaluation of the water diffusion coefficient in the two phases. Perhaps, a simulation restraining water molecules to the inside of the RTIL slab, not carried out in the present study, could provide a fair estimate of water diffusion in wet The absorption of water once again is a downhill process in free energy, that occurs spontaneously and rapidly, with an obvious change of time scale between T = 300 K and T = 360 K. In our qualitative analysis, the two processes, in [ , addition of water increases the absolute value of the electrostatic potential across the slab (see Figure S8), and thus, it also increases the size of |Δϕ| at the surface, the change being ∼6% both at T = 300 K and at T = 360 K. In other terms, in [DABu][Tf] water does not screen the electric field The straight green dash lines give the total and the water-related density of HBs in the homogeneous sample at the same conditions. The distribution of hydrogen bonds across the dry slab at the same temperature is given as a comparison by the black dotted line.
The Journal of Physical Chemistry B Article due to the underlying ionic liquid, but it enhances it. This effect can be rationalized by recognizing that the surface orientation of water does not depend only on the local electric field, but also on packing and especially hydrogen bonding.
The water vapor outside the slab turns out to be slightly denser than in the [DAEt][Tf] simulations at equal conditions (see Table 5). Comparison with the density of the corresponding condensed phases shows that the difference in excess chemical potential of water molecules in the vapor and within the slab is 10% smaller (in absolute value) than in [DAEt][Tf], both at T = 300 K and at T = 360 K. These results are verified and discussed in the next section.
III.C. Free Energy Computations. The simulation results showing that water is readily incorporated into the slab, and the estimate of the excess free energy difference for water molecules inside the slab and in the vapor, although interesting, still provide a limited view of the kinetics of water sorption at the RTIL surface.
To complete the picture, and to confirm that no shallow subsurface intermediate state is present, we computed the free energy profile for moving one water molecule across the surface. The method briefly described in section II and implemented into Gromacs has been applied. In the case of the wet [DAEt][Tf], 31 umbrella sampling windows have been selected, each point has been equilibrated for 105 ns with a time step of 0.5 fs. Statistics has been collected during production runs of 17.5 ns, again at each of the 31 windows of the umbrella sampling.
The results for [DAEt][Tf] at T = 300 K and T = 360 K are shown in Figure 12. At both temperatures, the f exc (z) curve is smooth, monotonic and virtually featureless, apart from a broad and very shallow minimum at z ∼ 12 Å for T = 360 K. The alignment of the two curves is only approximate, having been set by shifting each curve in such a way that they approach zero far from the slab. The transition between the condensed and the vapor regions takes place on a range of about 10 Å somewhat wider than the geometric width of the interface measured by the 90−10% density width. The offset of the two curves of ∼7 kJ/mol measured at the center of the slab agrees fairly well with the estimate obtained from the water density in the condensed and in the vapor phase. The discrepancy is likely to be due to the approximate alignment of the zeros.
A similar computation for the surface of a wet [DABu][Tf] slab at T = 360 K gives the results reported in Figure 13. Also, in this case, 31 windows have been selected, distributed along the z coordinate at 0.1 nm separation. Equilibration  , but also in this case the free energy profile is virtually monotonic. Comparison of the two curves at the same temperature (see Figure 13) confirms the estimate of Δμ H2O exc reported in Table 4 and in Table 5, obtained from an independent route, based on the computation of the vapor and condensed water densities coexisting at each temperature.    The Journal of Physical Chemistry B Article and N atoms, providing the electrostatic signature of the interface. The marked layering of ionic and neutral functional groups can be seen as the 2D manifestation of mesophase formation in a planar inhomogeneous system. No ion or ion pair is detected in the vapor side of the interface, consistently with the low vapor pressure of RTILs. Somewhat unexpectedly, the layering at the surface promotes the formation of HBs, probably because of the orientational ordering and high density of adjacent O and H 3 N− surface layers. Because of the close association of cations and anions, the electrification of the interface is not very high, despite the clear ionic character of these compounds. The surface layering and the alternation of different atomic species cause the charge density distribution to alternate in sign, giving origin to oscillations that at low T persist deep into the liquid phase. The generally low charge density and the charge oscillations result in a surface dipole of less than 1 D/nm 2 at all temperatures 280 K ≤ T ≤ 360 K, corresponding to an electrostatic potential drop Δϕ ∼ 0.25 V between the center of the slab and the vacuum.

IV. CONCLUSIONS
The surface dipole, or, equivalently, the Δϕ computed here for dry and wet slabs of the order of 200 meV, is relatively small on the scale of the values computed for other ionic liquid surfaces. For instance, the |Δϕ| at the surface of molten alkali halides at their triple point ranges from the 33 meV of NaF to the 1415 meV of LiI. 44 Also for molten alkali halides, values are both positive and negative, reflecting the relative position of anions and cations with respect to the reference surface plane.
Over the 10 2 ns time scale of the simulations, the mean square displacement of cations and anions in [DAEt][Tf] is not greatly affected by the confinement along z. Surface layering, instead has a sizable effect, decreasing the diffusion constant along z and enhancing it in the xy plane. These statements are based on the computation of the average displacement of all ions, irrespective of their position with respect to the surface. At all temperatures down to T = 300 K, both cations and anions jump from layer to layer every few nanoseconds, and a quantitative determination of z-dependent coefficients D z and D xy is practically unfeasible.
The time autocorrelation function ψ(t) of the electrostatic dipole normal to the surface has a purely relaxation character, without any oscillation that could point to a weakly damped plasmon. This observation implies that bulk and surface plasmon modes are either strongly damped, or quenched by the sizable association of cations and anions, as found in ref 19. Damping, in particular, might be enhanced by the complex and flexible frame of the ions, that increases the number and variety of relaxation modes. As a result, the time autocorrelation of dipole fluctuations normal to the surface can be reproduced by a few decaying exponential terms. At the lowest T of our simulations, a good fit is provided by the superposition of a stretched exponential of the KWW form and a simple decaying exponential, consistently with the near glassy state of the system.
Results To investigate the changes of surface properties due to water absorption, 24  case. The absolute value of the electrostatic potential across the slab, and, in particular, the corresponding |Δϕ|, are enhanced by ∼6% by water addition at low concentration, emphasizing the role of packing and hydrogen bonding in determining the orientation of water molecules at the surface.
Water addition at low concentration changes dielectric properties only slightly, but in a systematic way. This is apparent from the relaxation times of the dipole timecorrelation function, that decrease substantially upon water absorption. The shorter relaxation times, indicating that the dipole loses the memory of its direction faster, are consistent with the increases fluidity of wet samples and point to an increase of dissipation mechanisms due to the presence of water.
Wet slabs are in equilibrium with a tenuous vapor phase, whose density can be estimated by simulation. The results, shown in Table 4 and in Table 5, provides a way to estimate the difference in excess chemical potential between the condensed and the vapor phases. This information is complemented by the computation of the free energy profile felt by water molecules moving along the normal to the surface, determined by the application of the umbrella sampling algorithm of statistical mechanics. This determination is an essential part of our study, and shows that the excess free energy of water molecules, due primarily to their interaction with the ionic liquid, decreases monotonically from the vapor phase to the inside of the liquid slab. The absence of secondary, local minima is consistent with the rapid absorption of water and with its equally fast redistribution throughout the slab.
The Journal of Physical Chemistry B

Article
The description of the kinetics of water absorption and desorption through the RTIL surface has been completed by the determination of the absolute rate of water molecules impinging to the RTIL surface from the vapor phase, which can be computed by analyzing the simulation trajectories.
In conclusion, molecular dynamics simulation of an atomistic model has been used to explore water sorption through the free surface of [DAEt][Tf] and [DABU][Tf], providing a wealth of data and details whose determination by experimental methods is made difficult by the tiny fraction of the total system that is represented by its surface.