The 1-ethyl-3-methylimidazolium bis(trifluoro-methylsulfonyl)-imide ionic liquid nanodroplets on solid surfaces and in electric field: A molecular dynamics simulation study

Room temperature ionic liquids (RTILs) have promising potential for applications in a variety of systems/devices, including supercapacitors, Li-ion batteries, and field-effect transistors.^1–3 Extensive research interest has been shown on the RTIL-based supercapacitors, usually with carbon-based electrodes, including graphene-like structures and carbon nanotubes. Despite a comprehensive study of bulk RTIL behavior at charged surfaces both by experiments and simulations,^1,4–8 the wetting and electrowetting states of RTIL nanodroplets remain barely investigated. Inspired by initial studies demonstrating that moving electrolytes on a conducting surface can generate an electrical current,^9,10 droplets of ionic liquids have been considered for potential applications in energy harvesting, sensors, and electronic devices.¹¹ Also, controlled manipulation of nanodroplets at interfaces allows for the design of new generation of optoelectronic devices.^12,13

Many potential applications of droplets are heavily related to their electrowetting states on charged surfaces or in the presence of external electrical fields.^14,15 Electrowetting states were reported to be independent of the direction of electrical field but strongly affected by the strength of the field, which has also been confirmed by experimental studies of macroscopic RTIL droplets.^15,16 At nanoscale, several molecular dynamics (MD) simulations showed that polar fluid nanodroplets can show symmetric and asymmetric (with respect to the polarity of the surface) electrowetting states depending on the strength and orientation of the field.^17,18 MD simulation studies of RTIL nanodroplets in the weak electric fields (i.e., at voltages typically used in batteries and supercapacitors and generating electric field on the order 0.5 V/nm in the considered simulation setup) observed almost no change in the electrowetting states as a function of voltage.¹⁹ These predictions are consistent with recent experimental observations for several RTILs on the positive Pt electrode, where almost no change in the droplet contact angle was observed as a function of voltage.²⁰ However, the same experiments showed that on the negative electrodes a reduction of contact angles to smaller values was observed for all RTILs in the electrode potential range between −0.2 and −1.0 V.²⁰ At stronger electric fields, pronounced effects on electrowetting states of RTIL droplets have been demonstrated both in experiments^21–23 and molecular simulations.²⁴ In the latter, the corresponding electric field between electrodes was around 2.0 V/nm (surface charge of 0.0029 e/atom in that study). At these conditions, the asymmetric electrowetting states of RTIL droplets have been observed, showing a strong dependence on the electrode polarity.²² A recent MD simulation study of RTIL nanodroplets in vacuum with an external electric field of 1.2 V/nm showed the charge separation and increasing dipole moment in the droplet.²⁵ 

In this paper, using MD simulations, we continue to investigate electrowetting states by examining 1-ethyl-3-methylimidazolium bis(trifluoro-methylsulfonyl)-imide (EMIM-TFSI) RTIL nanodroplets on various surfaces. In this study, an external electric field is introduced by adding a surface charge. The electrode surface was represented by a single graphene sheet, where the interaction of the graphene sheet with the ionic liquid was varied to mimic possible variations in electrode materials and substrates with graphene coatings. Due to its extraordinary mechanical, optical, and electrical properties, graphene is showing great potential as a material of choice for electrodes, graphene-based coatings, and microscopic fluidic devices. Graphene has demonstrated to be wetting transparent as a coating on copper, gold, and silicon surfaces, i.e., showing that the graphene layer has minimal effect on the droplets’ wetting states, which are determined by interactions between the substrate and the droplet. The wetting transparency of graphene is a result of long-range van der Waals (VDW) interactions between the substrate and the liquid droplet.²⁶ Therefore, the behavior of droplets on graphene sheets can be efficiently manipulated by changing the substrate below the graphene layer. However, there are several studies that showed that graphene is only partially wetting transparent and the transparency can break down on super-hydrophilic and super-hydrophobic substrates.^27,28 Despite this ongoing debate, there is a clear realization that diverse wetting states of droplets on graphene can be achieved by carefully manipulating the substrate/droplet materials combination or by changing the number of graphene layers at the surface.^26–29 To investigate the wetting states of RTIL nanodroplets on different graphene-coated surfaces, in our simulations the VDW interactions between the RTIL and graphene sheet were rescaled. The different surfaces were modeled by rescaling the strength of Lennard-Jones (LJ) potential between the droplet and graphene atoms (⁠$εε0$⁠), therefore allowing a series of VDW parameters representing surfaces ranging from complete wetting to non-wetting by the investigated RTIL.

For the sake of computational efficiency, in this work we employed a combined united-atom (UA)/explicit-atom (EA) force field that has been developed based on the polarizable, fully atomistic APPLE&P force field³⁰ and previously validated and utilized in simulations of supercapacitors.^31,32 This UA/EA force field successfully captured structural, dynamic, and thermodynamic properties of bulk RTILs predicted by the APPLE&P force field.^32–34 In this force field, the repulsion/dispersion (RD) interactions between two atoms i and j of types α and β are described by a four-parameter function,

where r_ij is the separation between atoms, A_αβ and B_αβ are the coefficients describing the repulsion term of α–β pairwise interaction, C_αβ is the dispersion interaction parameter, and D_αβ is an additional repulsive parameter which dominates the repulsion at r_ij less than 1 Å to avoid inflection of U^RD(r) at very short distances. To simplify the scaling of VDW interactions between the RTIL nanodroplet and the graphene surface, we have approximated the corresponding interaction types with the LJ potential by fitting Eq. (1) to reproduce the dispersion part of interaction between RTIL and graphene atoms while sacrificing the accuracy in the description of repulsive interactions at very close distances. In the LJ potential, the strength of the U^RD interaction can be modified by simply scaling the depth of the LJ interaction ε. Here we investigated systems with graphene-droplet ratios of $εε0$ (where ε₀ is the original potential parameter) varying from 1.0 to 0.2, which as we discuss below corresponds to completely wetting conditions (⁠$εε0$ = 1.0) and non-wetting conditions (⁠$εε0$ = 0.2). Parameters for the unscaled VDW interactions between the RTIL atoms and graphene are given in the supplementary material.

Simulations were conducted using the LAMMPS package.³⁵ The cut-off radius for the LJ interaction was set at 10 Å while for Coulomb interactions at 15 Å. The Particle-Particle-Particle-Mesh (PPPM) approach was employed to treat the long-range electrostatic interactions in the k-space with 10⁻⁴ force accuracy. Graphene atoms were fixed (i.e., forces acting on these atoms were set to zero) during simulation. System preparation was done following a three-step process. The bulk phase of RTIL was simulated in an NPT ensemble at 393 K and 1 atm for 10 ns. Then, a cubic-shape droplet of RTIL (unwrapped from periodic boundary conditions) was placed within 5 Å from the graphene sheet. Subsequently, MD simulations in the NVT ensemble at 393 K were conducted to allow relaxation of the RTIL to an equilibrium droplet shape. Production simulations were conducted once the shape of the droplet was stationary. The integration time step for simulations was 0.5 fs, and the total length of trajectory was longer than 40 ns for each system. To investigate the size effect on contact angles, we investigated systems with the number of IL ion pairs (N) of 100, 200, and 300. In previous MD simulations of wetting states of RTIL nanodroplets, it was shown that the contact angle predictions are converging when the number of ionic pairs is about 200 pairs.³⁶ 

A typical snapshot of a simulation cell is shown in Fig. 1. In order to keep the charge neutrality of the system during simulation with charged surfaces, the simulation cell contained both the negative and the positive surfaces each having an RTIL droplet as shown in Fig. 1(c). Charged surfaces with ±0.002 e or ±0.004 e on each carbon atom of graphene were investigated. These surface charge densities correspond to electric fields of 1.36 and 2.72 V/nm, respectively, for the cell geometry shown in Fig. 1(c) and assuming vacuum between two surfaces.

To calculate contact angles, we used the analysis of ions’ number density, following the procedure established in previous studies.^37,38 First, a three-dimensional grid was generated through a simulation cell, and for each point, the averaged number density was calculated every 1 ps from a 4 ns trajectory. Then, the density profile is calculated within a series of two-dimensional planes that are perpendicular to the surface and include the droplet center-of-mass. For each droplet, we have divided the trajectory into ten 4 ns intervals. For each interval, 1600 snapshots were analyzed, with each snapshot providing n = 4 planes rotated by an angle $πn$ (units in radian) between two neighboring planes, to generate the averaged number density profile. The resulting average 2D density profile of the droplet obtained from one of the 4 ns intervals is shown in Fig. 2. The contour line of the droplet is defined by an isosurface which has a value of the total density equal to half of the RTIL density in the middle of the droplet. A polynomial function is fitted to define the droplet contour line that subsequently allows defining the contact angle (θ) as illustrated in Fig. 2. We tested polynomial functions of different orders and found that there were no differences in the extracted contact angle upon increasing the polynomial order from 5th to 6th order. Therefore, all profiles were fitted using the 6th order polynomial function. This procedure was applied for each of the ten trajectory intervals considered as independent measurements of the contact angle. Standard error analysis was applied to estimate the error bars for the predicted contact angle for each system.

The wetting states of nanodroplets on graphene are strongly dependent on the strength of VDW potential between RTIL and graphene atoms. This effect is shown in Fig. 3 where a nanodroplet with N = 100 is totally wetting the surface at $εε0$ = 1.0, while a clear increase in the contact angle can be seen as $εε0$ reduces. The influence of the droplet size on the contact angle has been investigated by comparing droplets with N = 100, 200, and 300, as shown in Fig. 4. For all droplet sizes, the contact angle decreases almost linearly with the increase of $εε0$⁠, with perhaps some deviation from the linear behavior at $εε0>0.6$⁠. While there are some changes in the contact angle with increasing droplet size, those changes are not strong and qualitatively the dependence of the contact angle does not change. This is consistent with the recent work where the size effect for 1-ethyl-3-methylimidazolium tetrafluoroborate (EMIM-BF4) droplets has been thoroughly discussed, and properties were shown to converge for droplet sizes between 100 and 200 ionic pairs.³⁶ In another work, a linear dependence of contact angle for 1-ethyl-3-methylimidazolium glycine (EMIM-GLY) droplets was shown as a function of the droplet size.²⁹ Certainly the dependence of the contact angle on a nanodroplet size can vary depending on the chemical details/interactions within the droplet, strength and directionality of the field, and other factors. For the system investigated here, Fig. 4 shows that the difference between 200 and 300 ionic pairs is not substantial, and therefore further increase in the nanodroplet size is unlikely to affect the observed trends.

Our previous simulations of bulk RTILs at graphite surfaces have shown that EMIM cation has more favorable interaction with the neutral graphite (and graphene) surface compared to TFSI anion.^32,39 The ion number density profiles as a function of the distance to the graphene surface are shown in Fig. 5 for different interactions. From Fig. 5(a), it can be seen that nanodroplets with strong RTIL-graphene VDW potential (⁠$εε0$ = 1.0) are totally wetting on the graphene surface, and a single layer of ions was observed. As $εε0$ decreases [Figs. 5(b) and 5(c)], more ions are distributed further away from the graphene surface. For all levels of surface-RTIL interactions, cations are closer to the surface than the anions and ion layering is noticeable in all systems. Figs. 5(d)–5(f), which show fractions of cations and ions at a given separation from the surface, further illustrate the oscillatory density profile of the multi-layer structure inside nanodroplets. At lower $εε0$ ratio (less than 0.4), these oscillations are weaker, and at $εε0$ = 0.2, the nanodroplet approaches a bulk-like behavior (i.e., weak layering) as expected.

Next, we investigated electrowetting states of EMIM-TFSI droplets on the surfaces with different VDW interactions and, hence, different wetting properties. An electric field was introduced by placing point charges at the locations of graphene carbon atoms. Two values of surface charges ±0.002 e and ±0.004 e per carbon atom were investigated which correspond to levels of electric fields of 1.36 V/nm and 2.72 V/nm, respectively. The electrowetting states of RTIL nanodroplets in the electric field are shown in Fig. 6. The external electric field imposed by the surface causes asymmetrical wetting states of nanodroplets as illustrated in Fig. 6 by the side-by-side comparison in each panel of the half-droplet profiles from the positive and negative electrodes. At charged surface, the droplet is experiencing two competing effects. The ions near the surface respond to the surface charge leading to increased interaction between the surface and the layer of droplets close to the surface. The ions which are on the opposite side (top) of the droplet respond to the presence of electric field created by two charged surfaces, and hence here the droplet will tend to elongate along the field direction.²⁵ These competing effects lead to a non-trivial reshaping of the droplet on charged surfaces as illustrated in Fig. 6. The droplet asymmetry also shows strong dependence on the strength of VDW interaction between the surface and RTIL. For the $εε0$ = 1.0 and $εε0$ = 0.8 RTIL-surface interactions, EMIM-TFSI droplets form a single layer spreading on the graphene sheet. For weaker surface-RTIL interactions (⁠$εε0$ = 0.6–0.2) in 1.36 V/nm field (±0.002 e surface charge), the height of nanodroplets is increasing. For example, for $εε0$ = 0.6 and $εε0$ = 0.4, nanodroplets on the charged surface are about 8 Å and 11 Å higher than the corresponding droplets on a neutral surface. Also, a slight asymmetry in the droplet shape on the positive and negative electrodes can be seen for these systems. For $εε0$ = 0.6, the droplet on the positive electrode (left profile) is slightly higher than the droplet on the negative electrode (right profile). As $εε0$ reduced to 0.2 (non-wetting surface), the asymmetry is substantially increased, and the droplet on the positively charged surface is about 13 Å higher compared to the one on the negative surface. At the stronger electric field (2.72 V/nm or ±0.004 e), the interaction of ions with the surface becomes stronger and nanodroplets show an increased tendency for wettability. While the droplets still try to elongate along the field direction, the interaction with the surface dominates, and as a result, the droplet height decreases slightly compared to that on the neutral surface. At this stronger field and $εε0$ = 0.6, we also observe the shape asymmetry between droplets on the positively and negatively charged surfaces, but now the asymmetry is mostly in the lateral spreading of the droplet, rather than in the vertical height. As the strength of VDW interactions with the surface decreases [Fig. 6(f)], the asymmetry in height emerges but to a smaller extent compared to the corresponding system in the weaker electric field [compare Figs. 6(c) and 6(f)].

Figure 7 shows the dependence of contact angle on the strength of electric field for surfaces with different strengths of VDW interaction. For the weak and moderate VDW interactions (⁠$εε0$ ≤ 0.6), a change in contact angle upon switching from neutral to the surface with 0.002 e charge (1.36 V/nm electric field) is not significant (less than 10°), indicating that this electric field is not strong enough to noticeably influence the nanodroplet shape. The spreading of RTIL droplets was observed at $εε0$ = 0.8, showing a complete wetting at this charged surface. The asymmetry in the droplet structure/shape on the positive and negative surfaces is also manifested in the contact angle. The contact angle on the negative surface is somewhat lower than that on the positively charged surface. The contact angle on the positively charged surface was calculated to be 121.9°, 100.2°, and 70.0° for $εε0$ = 0.2, 0.4, and 0.6, respectively, while on the negatively charged surface, values of 114.5°, 96.6°, and 64.7° were found for corresponding VDW interactions. In the recent MD simulations with RTIL bridging two charged graphene surfaces,²⁴ the asymmetric contact angle induced by the surface charge was also reported. At a charge density of ±0.001 962 e/atom, a 6.6° difference in the contact angle has been observed. Similar differences between contact angles on the positive and negative surface are observed in our simulations at 1.36 V/nm electric field (±0.002 e/atom).

As the electric field increases to 2.72 V/nm, the dependence of contact angle on charge density is significant, showing a pronounced drop in the contact angle. The nanodroplets begin to show better wettability even for surfaces with weak VDW interactions $εε0$ = 0.4 and $εε0$ = 0.2, where more than 50° decrease in the contact angle has been observed compared to the neutral surface. At this field, the interaction between ions and charges on the surface is a dominant factor. Also, the asymmetry in the contact angle becomes more prominent. For the positively charged surface, the contact angles for $εε0$ = 0.6 and $εε0$ = 0.4 are 33.8° and 51.0°, respectively, while 24.7° and 45.2° were obtained for the corresponding nanodroplets on the negatively charged surface.

The direction of the external field has been shown to result in a non-trivial effect on electrowetting states of nanodroplets.^17,18,40 Several MD simulations have been conducted on insulating surfaces with an external electric field parallel and perpendicular to the surface.^17,40–42 In all cases, the nanodroplets have been observed to expand in the field direction.¹⁷ In our simulation, for nanodroplets on a positively charged surface, the effects of Coulombic interaction and VDW interaction are competitive. EMIM is favored by the VDW interaction, yet, “disliked” by positive charges. On the negatively charged surface, both the VDW interaction and Coulombic potential lead to a favorable interaction with EMIM. Therefore, the wetting states are expected to be asymmetric on surfaces with different polarity.

The ion number density profiles for nanodroplets in 1.36 V/nm electric field are shown in Figs. 8(a)–8(c), and the corresponding ion fractions are shown in Figs. 8(d)–8(f). The positively and negatively charged surfaces are indicated by P and N, respectively. As can be seen from Figs. 8(a)–8(c), nanodroplets on the positively and negatively charged surface show a different packing structure and this disparity increases with the weakening of the strength of the VDW interaction between the droplet and surface. At VDW interaction $εε0$ = 1.0 (wetting case), on both electrodes the cation and anion layers are merged into one layer within the distance less than 6.2 Å from the surface [Figs. 8(a) and (d)]. At the surface with $εε0$ = 0.6 VDW interaction, the height of the droplet [defined by the extent of density profile in the Z-direction, Fig. 8(b)] and the oscillatory layering of cations and anions are very similar for the positive and negative electrodes [Fig. 8(e)]. It is also not substantially different from profiles at the neutral surface, and therefore the contact angles on electrodes with this charge state are only slightly lower than those on the neutral surface. For these systems, we can conclude that the VDW interaction between the RTIL droplet and the surface is the dominant interaction defining droplet shape. As the VDW interactions reduced to $εε0$ = 0.2, the polarity of the surface significantly influences the ion distribution and droplet shape. As can be seen from Fig. 8(c), on the positively charged surface, the total density of ions near the surface is substantially lower than that on the negative electrode. As a result, the droplet extension in the Z-direction (height) is noticeably different. At these conditions (VDW interaction and charge density), the energy of electrostatic interaction of ions with the charged surface becomes comparable to the difference in the VDW interaction of cations and anions with the surface.

The ionic number distributions for nanodroplets in 2.72 V/nm electric field are shown in Fig. 9. The preference of cations on a positively charged surface holds true only at $εε0$ = 1.0. For VDW interaction weaker than $εε0$ = 1.0, the electrostatic interaction with the surface is the overwhelming factor that dominates the RTIL-graphene interaction and cation-anion separation. As can be seen from Figs. 9(b) and 9(c), anions are noticeably favored by the positively charged surface, i.e., almost no cations are located near the positively charged surface. An enhanced cation-anion separation has also been observed for nanodroplets on the negatively charged surface.

The layered structure of ions in the droplet and how it changes with increasing of surface charge are further illustrated in Fig. 10 where the average charge density profiles are shown for the system with the intermediate strength of VDW interaction. Consistent with the density profiles shown in Figs. 8 and 9, we can see a change in the ordering of RTIL charge layering as a function of electrode surface charge.

To investigate transitions from poor to good wetting states, a series of simulations with different electric fields were conducted for the system with $εε0=0.4$⁠. Figure 11 shows the contact angle of the droplet in this system on both electrodes as a function of surface charge. Both on positively and negatively charged surfaces, a transition of wetting states takes place at charge density between 0.0025 and 0.003 e/atom, which corresponds to 2.04 V/nm and 2.38 V/nm electric fields. For surface charges between 0.0020 e/atom and 0.0025 e/atom, only a slight decrease in contact angle (by 8.0° and 4.4°) has been observed for droplets on both charged surfaces. As the surface charge increases to 0.0030 e/atom, a significant drop in the contact angle can be seen, demonstrating a transition from poor to good wettability.

Wetting and electrowetting states of EMIM-TFSI room temperature ionic liquid on graphene sheets have been studied via MD simulations. By changing ε of the Lennard-Jones potential between ionic liquid atoms and graphene atoms, the surface-nanodroplet interactions were adjusted to represent different types of surfaces that can be obtained with graphene coating. Nano-sized droplets, with various wetting states ranging from a contact angle larger than 90° (non-wetting) at $εε0$ = 0.4 and $εε0$ = 0.2 VDW interaction to complete wetting at $εε0$ = 1.0 were investigated. The electrowetting behavior of RTIL nanodroplets was studied by introducing an external electric field imposed by surface charges on two graphene electrodes. Variation in electrowetting states on surfaces with different charge densities was investigated at 1.36 V/nm and 2.72 V/nm electric fields. Asymmetric shape of nanodroplets was observed on the positively and negatively charged surface, and the extent of this asymmetry depends on the strength of the VDW interaction between the RTIL and the surface as well as the magnitude of the electric field. At 1.36 V/nm field, the nanodroplets expanded in the direction perpendicular to the surface (along the electric field) compared to nanodroplets on a neutral surface. At this field, the asymmetry in droplet height is the most pronounced for non-wetting surfaces (⁠$εε0$ = 0.2), where a 13 Å difference in droplet height was observed between positive and negative surfaces. At the higher electric field of 2.72 V/nm, “good” wettability (showing contact angle less than 90°) has been observed for all cases including those with the weak VDW interaction between the RTIL and surface. At this electric field, the asymmetry of wetting states on positive and negative electrodes was also observed in the direction parallel to the surface, therefore mostly affecting the contact angle rather than the height of the droplet.

See supplementary material for force field parameters used to describe van der Waals interactions in simulated systems.

Authors gratefully acknowledge the support from the project sponsored by the Army Research Laboratory under Cooperative Agreement No. W911NF-12-2-0023. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of ARL or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein. We also would like to acknowledge the Center of High Performance Computing at the University of Utah for technical support and generous allocation of computing resources.