The recent predictions of the self-consistent generalized Langevin equation theory, describing the existence of unusual partially arrested states in the context of ionic liquids, were probed using all-atom molecular dynamics simulations of a room-temperature ionic liquid. We have found a slower diffusion of the smaller anions compared with the large cations for a wide range of temperatures. The arrest mechanism consists on the formation of a strongly repulsive glass by the anions, stabilized by the long range electrostatic potential. The diffusion of the less repulsive cations occurs through the holes left by the small particles. All of our observations in the simulated system coincide with the theoretical picture.

Understanding all aspects of glass transition is one of the great challenges in modern science.1–3 Despite much research having been done in recent years,4 the possibility of observing mixed arrested states, in which only one of the components becomes glass while the rest of the system remains fluid, has received less attention.5 

Recently, the Self-Consistent Generalized Langevin Equation (SCGLE) theory was used to investigate the dynamically arrested states in the context of the primitive model (PM) of molten salts.6 The most important result of the cited work was the prediction of low-temperature regions where the larger ions are fluid and the smaller ions are arrested. These domains of mixed states are related with the size and charge asymmetries between ions and are accessible at very low temperatures.

Such kind of partially arrested states has been observed in materials known as superionic conductors,7 in which one type of the constitutive ions is able to move inside the solid matrix formed by the other components. Crystalline and glassy superionic conductors were extensively used in the industry of batteries since the 1980s.8 However, superionic conductors used to have a very complex composition, making difficult their direct comparison with our theoretical approach. Curiously, a novel class of ionic fluids, known as Room-Temperature Ionic Liquids (RTIL), possesses a glassy-behavior at room temperatures and resulted more adequate for probing the theoretical predictions.9,10 Clearly, due to that both systems are Coulombic fluids, most of the conclusions obtained in RTIL can be applied to superionic conductors.

RTIL are compounds usually formed by large organic cations combined with organic or inorganic anions which melt at temperatures less than 100°C.11 The study of these solvents has become very active due to their potential applicability towards industry.12 Despite the intense research around RTIL, scarce investigation has been done with respect to their glassy behavior.10,13–21 Many researchers have previously observed a slower diffusion of the smaller anions in RTIL, either by experiments21 or by molecular dynamics (MD) simulations.10,13–20 The most common explanation appeals to the dynamics along the direction of the carbon located between the two nitrogens in the imidazolium ring with a marginal dependence on the anion.22,23 However, the present work offers an alternative explanation for this observation on RTIL, one that is supported by an analysis of the simulation results in light of the SCGLE theoretical predictions. The existence of partially arrested states in RTIL may open the possibility for developing novel conducting materials near room temperature.

We studied the glassy behavior of the generic RTIL 1-ethyl-3-methyl-imidazolium tetrafluoroborate (EMIM+-BF4) using all-atom molecular dynamics (MD) simulations. For this particular compound, we have observed the onset in the formation of partially arrested states in which the small anion (BF4) becomes arrested while the large cation (EMIM+) is still fluid. In addition, a qualitative comparison between theory and simulation is presented in order to contrast with the general picture given by theoretical predictions. We can conclude that the partially arrested states in ionic liquids previously predicted by the SCGLE theory are consistent with our MD simulations of EMIM+-BF4.

The SCGLE theory of dynamic arrest leads to a remarkably simple equation for the asymptotic value γαlimt(ΔR(α))2 of the mean-squared displacement (MSD) of particles of species α. This property constitutes a convenient parameter with which to describe dynamic arrest in mixtures, since it indicates if the particles of certain species diffuse (infinite value of γα) or are arrested (finite value of γα). The existence of a partially arrested state implies the arrest of one component and the diffusion of the other one. The referred equation reads as24 

1γα=13(2π)3d3kk2{λ[λ+k2γ]1}αα×{cnSλ[Sλ+k2γ]1nh}αα,
(1)

where S is the matrix of partial structure factors and h and c are the Ornstein-Zernike (OZ) matrices of total and direct correlation functions, respectively. The matrix n is defined as [n]αβδαβnα, and λ(k) is a diagonal matrix given by λαβ(k)=δαβ[1+(k/kc(α))2]1, with kc(α)=8.17/σα.

We used a model of charged hard spheres, usually referred to as the primitive model (PM), for which the analytical solution of the OZ equation was given by Blum under the mean-sphere approximation (MSA).25 

We have performed all-atom MD simulations of ionic liquid EMIM+-BF4 using the MD package GROMACS.26 The BF4 anion has a tetrahedron structure and a mean size of 0.19 nm, while the EMIM+ cation is essentially planar and with a maximum size (taken from methyl to ethyl) of ∼0.69 nm (see Fig. 1). The force field of EMIM+ is based on the AMBER potential,27 and the BFBF4 force field was developed by de Andrade and colleagues.28 Periodic boundary conditions were applied to the three dimensions with a 1.6 nm cutoff for the atomic interactions. The long-range electrostatic interactions were handled with the particle mesh Ewald method.29 As an initial configuration for the EMIM+-BF4 RTIL, we have used a simulation box with 1024 ion pairs randomly positioned.

FIG. 1.

Schematic representation of the ions studied by MD simulation.

FIG. 1.

Schematic representation of the ions studied by MD simulation.

Close modal

The equilibration procedure is the same as described in Ref. 30. The equilibrated sample was probed from 700 K to 250 K in temperature steps of ΔT = 50 K. We simulated the system for a duration of 10 ns (10 000 configurations) at each temperature for the calculation of the mean-square displacement (MSD), the diffusion coefficient, and the radial distribution functions (RDFs). For obtaining reliable results at lower temperatures, we have used 100 ns durations for temperatures below 350 K.

Isobaric trajectories (P = 1 atm) were used in both simulation and theory. We employed the dimensionless quantitiescommonly used in the PM for: pressure PPϵ0σbs4/e2 (σbs(σbig+σsmall)/2, where σα is the diameter of the corresponding ion), temperature TkBTσbsϵ0/e2, and total volume fraction ϕπ(nbigσbig3+nsmallσsmall3)/6. For all the previous definitions, ϵ0 is the vacuum permittivity, kB is the Boltzmann constant, e is the electron charge, and nα corresponds to the number density of particles α. We have chosen a size asymmetry of δ=3.5 (δ=σbig/σsmall) in our theoretical calculations because it is close to the size ratio between EMIM+ and BF4 molecules. Clearly, the simulated ions are not charged hard spheres; therefore, we expect incorrect values of temperature and volume fraction for their dynamic arrest using the PM. However, let us suppose that the essential features of ionic liquids can be described by the PM.

In Fig. 2 we show the dynamic arrest line, calculated with Eq. (1), in terms of variables T and ϕ. The solid line corresponds to the transition to fully arrested states (G-G) in which both ions are arrested. The dashed line separates the fluid region from the partially arrested states (F-G) where only one of the components is arrested. The dotted line represents the isobaric trajectory used for our additional theoretical calculations. For reference, we have marked three representative points with large circles corresponding to T = 0.079, 0.025, and 0.012. In the inset, we plot the inverse of the localization length 1/γα as a function of T over the isobaric trajectory. The squares and the circles correspond to the data of the cations (large particles) and anions (small particles), respectively. We can observe that 1/γanion becomes different from zero at higher temperatures (T≈0.013) than compared with 1/γcation (T≈0.01), indicating that inside the F-G region the anion is arrested and the cations are liquid.

FIG. 2.

Arrest lines predicted by the SCGLE theory for the PM with a size asymmetry adequate for the simulated RTIL (1:3.5). The fluid region is labeled as F-F. The G-G region corresponds to fully arrested states. Inside the F-G region we found partially arrested. The dotted line is an isobaric trajectory equivalent to 1 atm of pressure. The points over the dotted line correspond to the following effective temperatures T = 0.079, 0.025, and 0.012. In the inset we plot the values of the parameter 1/γα following the isobaric trajectory. Circles correspond to anions (small particles) and squares to cations (large particles).

FIG. 2.

Arrest lines predicted by the SCGLE theory for the PM with a size asymmetry adequate for the simulated RTIL (1:3.5). The fluid region is labeled as F-F. The G-G region corresponds to fully arrested states. Inside the F-G region we found partially arrested. The dotted line is an isobaric trajectory equivalent to 1 atm of pressure. The points over the dotted line correspond to the following effective temperatures T = 0.079, 0.025, and 0.012. In the inset we plot the values of the parameter 1/γα following the isobaric trajectory. Circles correspond to anions (small particles) and squares to cations (large particles).

Close modal

The prediction of partially arrested states is confirmed by observing the inset of Fig. 3, where the mean square displacement (MSD) of the large cations (dashed lines) and small anions (solid lines) at temperatures T = 0.079, 0.025, and 0.012 are shown. We can observe that, even at temperatures above the dynamic arrest transition (T = 0.079), the slope of the MSD of the anion has a lower value compared with the slope of the cation. At lower temperatures (T = 0.025), the difference between the slopes increases. Finally, inside region F-G (T = 0.012) the anions become arrested while the cations remain mobile.

FIG. 3.

MSD of EMIM+-BF4 for the temperatures T = 0.0015 (700 K), 0.0011 (500 K), and 0.000 65 (300 K). Dashed lines represent our results for the cations and solid lines represent anions. At the lower temperature, the anion in MSD is almost bounded, indicating dynamic arrest. Inset: Theoretical calculations of the mean-square displacement for T = 0.079, 0.025, and 0.012 over the isobaric trajectory used in this work (big circles in Fig. 2). The first two temperatures are located in the F-F region and the last one inside the F-G region. Solid and dashed lines correspond to anions and cations, respectively.

FIG. 3.

MSD of EMIM+-BF4 for the temperatures T = 0.0015 (700 K), 0.0011 (500 K), and 0.000 65 (300 K). Dashed lines represent our results for the cations and solid lines represent anions. At the lower temperature, the anion in MSD is almost bounded, indicating dynamic arrest. Inset: Theoretical calculations of the mean-square displacement for T = 0.079, 0.025, and 0.012 over the isobaric trajectory used in this work (big circles in Fig. 2). The first two temperatures are located in the F-F region and the last one inside the F-G region. Solid and dashed lines correspond to anions and cations, respectively.

Close modal

The main panel of Fig. 3 contains our MD simulation results for the MSD. We show our calculations for three representative temperatures, i.e., T= 0.0015 (700 K), 0.0011 (500 K), and 0.000 65 (300 K). The solid and dashed lines correspond to the MSD of the BF4 anion (small particles) and the EMIM+ cation (large particles), respectively. At short times, we observe a faster mobility of BF4 compared with EMIM+ due to that, in the ballistic regime, the molecules perform a free flight dominated by the thermal velocity vT=KBT/M, which is scaled with the mass of the molecules. Hence, the massive EMIM+ (∼111 amu) moves slower than BF4 (∼87 amu). However, in the diffusive regime (long time) the molecular collisions dominate, and the diffusion of BF4 becomes slower than the EMIM+ diffusion. The differences between theory and simulations for short times are due to lacking of the ballistic regime in the SCGLE theory (originally developed for describing colloids). However, the dynamic arrest transition is a long-time effect. Thus, the behavior of the simulated RTIL at large times essentially reproduces the theoretical predictions for the PM.

The full squares in Fig. 4 correspond to our MD results of the diffusion coefficient of anions and the full circles to the cations, the lines are only guides for the eye. The intuitive idea that bigger particles have slower diffusion contrasts with the observed faster diffusion of the bigger EMIM+ in our MD simulations. In order to reinforce the validity of our results, we included experimental measurements (empty symbols) of diffusion coefficients (taken from Ref. 21), which confirms the tendency observed in the simulations. The comparison between simulation and experiments may be improved upon using polarizable force fields, for which the net effect is that the system becomes more mobile.31–33 However, the complexity of the force fields is surely not affecting the observation of anomalous diffusion of RTIL.

FIG. 4.

Diffusion coefficient as a function of the temperature of the EMIM+BF4 system. Lower diffusion coefficients were observed for the BF4 anion in a wide range of temperatures. Empty symbols represent experimental measurements showing the same tendency.21 Inset: Theoretical calculation of the diffusion coefficient as a function of T in the context of the primitive model. The solid line corresponds to the diffusion coefficient of anions and the dashed line to cations.

FIG. 4.

Diffusion coefficient as a function of the temperature of the EMIM+BF4 system. Lower diffusion coefficients were observed for the BF4 anion in a wide range of temperatures. Empty symbols represent experimental measurements showing the same tendency.21 Inset: Theoretical calculation of the diffusion coefficient as a function of T in the context of the primitive model. The solid line corresponds to the diffusion coefficient of anions and the dashed line to cations.

Close modal

As reference, the inset of Fig. 4 shows the SCGLE calculation of the self-diffusion coefficient as a function of T over the isobaric trajectory, which crosses through the F-C region. The solid line corresponds to the diffusion coefficient of cations and the dashed line to that of anions. We can observe that the diffusion coefficient vanishes at higher temperature for the anions than for the cations, implying the arrest of the anions and the diffusion of the larger cations. Despite the complexity of the constituent ions of RTIL and the fact that we have performed all-atom simulations, both simulations and experiments are in complete qualitative agreement with the predictions of the SCGLE theory. The evidence presented above shows the existence of partially arrested states in RTIL, where the small particles are arrested and the large particles are mobile.

In Fig. 5 we show the radial distribution functions (RDFs) for T = 0.0015 (red), 0.0011 (green), and 0.000 65 (blue), where the dashed, dashed-dotted, and solid lines correspond to g+−(r), g++(r), and g−−(r), respectively. Notice that the first peak of g+−(r) is located approximately at r ≈ 0.435 nm and the first peak of g++(r) can be found around 0.7 nm for all the temperatures probed. On the other hand, the first peak position of g−−(r) decreases from 0.73 nm to 0.63 nm as temperature lowers. From Fig. 5, we can see that while the g+−(r) peak is located at r ≈ 0.45 nm (σcation+σanion)/2 and g++(r) at r ≈ 0.68 nm σcation, the first peak of g−−(r) is located in between 0.6 nm and 0.7 nm which is larger than σanion. This behavior is typical of a Wigner glass in which the mean distance of the constituents of the glass is very large.34 

FIG. 5.

Radial distribution functions for EMIM+-BF4. The dashed, dashed-dotted, and solid lines represent g+−(r), g++(r), and g−−(r), respectively. We have plotted three representative temperatures, i.e., T = 0.0015 (red), 0.0011 (green), and 0.000 65 (blue).

FIG. 5.

Radial distribution functions for EMIM+-BF4. The dashed, dashed-dotted, and solid lines represent g+−(r), g++(r), and g−−(r), respectively. We have plotted three representative temperatures, i.e., T = 0.0015 (red), 0.0011 (green), and 0.000 65 (blue).

Close modal

What we are observing is that the cations have a low efficiency on screening the charge of the anions, due to their large volume, which implies that the anion-anion repulsion is greater than the cation-cation repulsion. In consequence, the anions become a glass due to the high repulsion between them, and the less repulsive cations have enough space for diffusion.

In Fig. 6 we show a snapshot of the system at T = 0.000 65. The BF4 anions are colored with green and the cations are represented as a continuous gray surface. We can observe that the mean distance of the anions is very large, showing the formation of a strongly repulsive glass. The space between the anions makes possible the diffusion of the less repulsive cations.

FIG. 6.

Snapshot of the system at T = 0.000 65. The green particles are the anions and the continuous gray surface represents the space occupied by the cations.

FIG. 6.

Snapshot of the system at T = 0.000 65. The green particles are the anions and the continuous gray surface represents the space occupied by the cations.

Close modal

In the present work we have probed the predictions of the SCGLE about the existence of partially arrested states at low temperatures in the context of ionic liquids using all-atom molecular dynamic simulations of EMIM+-BF4ionic liquid. The predicted arrested states are characterized by an atypical behavior in which the small particles, the anions in this case, are arrested while the large particles remain mobile. We have found that the essential features predicted by the SCGLE theory in the onset of glass transition were observed in our simulated RTIL. The most important common features are the following: (1) the diffusion of the anions, observed in the MSD and the diffusion coefficient, is slower for a wide range of temperatures and (2) the MSD of the anions becomes almost bounded at a higher temperature than the cations. In addition, we can propose the formation of a strongly repulsive glass by the anions as the arrest mechanism in these unusual states. The qualitative agreement between the theoretical prediction of these singular arrested states with simulation and experimental results gives a proof of the reliability of the theoretical predictions.

The authors appreciate the assistance of Monica I. Ihili and J. Limón Castillo. We acknowledge the financial support of National Basic Research Program of China (973 program, No. 2013CB932804), the National Natural Science Foundation of China (Nos. 11274319 and 11421063), and CONACYT through grants: Cátedras CONACyT No. 1631 and CB-2015-01 No. 257636.

1.
F. H.
Stillinger
and
P. G.
Debenedetti
,
Annu. Rev. Condens. Matter Phys.
4
,
263
(
2013
).
2.
P. G.
Debenedetti
and
F. H.
Stillinger
,
Nature
410
,
259
(
2001
).
3.
C. A.
Angell
,
K. L.
Ngai
,
G. B.
McKenna
,
P. F.
McMillan
, and
S. W.
Martin
,
J. Appl. Phys.
88
,
3113
(
2000
).
4.
L.
Berthier
and
G.
Biroli
,
Rev. Mod. Phys.
83
,
587
(
2011
).
5.
J.
Hendricks
,
R.
Capellmann
,
A. B.
Schofield
,
S. U.
Egelhaaf
, and
M.
Laurati
Phys. Rev. E
91
,
032308
(
2015
), and references therein.
6.
L. E.
Sanchez-Díaz
,
A.
Vizcarra-Rendón
, and
R.
Juárez-Maldonado
,
Phys. Rev. Lett.
103
,
035701
(
2009
).
7.
A.
Laskar
,
Superionic Solids and Solid Electrolytes Recent Trends
(
Elsevier
,
2012
).
8.
A.
Varzi
,
R.
Raccichini
,
S.
Passerini
, and
B.
Scrosati
, “
Challenges and prospects of the role of solid electrolytes in the revitalization of lithium metal batteries
,”
J. Mater. Chem. A
(published online).
9.
M. G.
Del Popolo
and
G. A.
Voth
,
J. Phys. Chem. B
108
,
1744
(
2004
).
10.
J.
Habasaki
and
K. L.
Ngai
,
J. Chem. Phys.
142
,
164501
(
2015
).
11.
W.
Freyland
,
Coulombic Fluids: Bluk and Interfaces
(
Springer
,
2011
).
12.
N. V.
Plechkova
and
K. R.
Seddon
,
Chem. Soc. Rev.
37
,
123
(
2008
).
13.
A.
Triolo
 et al,
J. Phys. Chem. B
110
,
21357
(
2006
);
[PubMed]
K. R.
Harris
,
L. A.
Woolf
, and
M.
Kanakubo
,
J. Chem. Eng. Data
50
,
1777
(
2005
).
14.
J.
Habasaki
and
K. L.
Ngai
,
Anal. Sci.
24
,
1321
(
2008
).
15.
J.
Habasaki
and
K. L.
Ngai
,
J. Chem. Phys.
129
,
194501
(
2008
).
16.
J.
Habasaki
and
K. L.
Ngai
,
J. Non-Cryst. Sol.
357
,
446
(
2011
).
17.
N. C.
Forero-Martinez
,
R.
Cortes-Huerto
, and
P.
Ballone
,
J. Chem. Phys.
136
,
204510
(
2012
).
18.
O.
Andreussi
and
N.
Marzari
,
J. Chem. Phys.
137
,
044508
(
2012
).
19.
M. H.
Kowsari
,
S.
Alavi
,
M.
Ashrafizaadeh
, and
B.
Najafi
,
J. Chem. Phys.
129
,
224508
(
2008
).
20.
H. A.
Every
,
A. G.
Bishop
,
D. R.
MacFarlane
,
G.
Oradd
, and
M.
Forsyth
,
Phys. Chem. Chem. Phys.
6
,
1758
(
2004
).
21.
A.
Noda
,
K.
Hayamizu
, and
M.
Watanabe
,
J. Phys. Chem. B
105
,
4603
(
2001
).
22.
S. M.
Urahata
and
M. C. C.
Ribeiro
,
J. Chem. Phys.
120
,
1855
(
2004
).
23.
S. M.
Urahata
and
M. C. C.
Ribeiro
,
J. Chem. Phys.
122
,
024511
(
2005
).
24.
R.
Juarez-Maldonado
and
M.
Medina-Noyola
,
Phys. Rev. E
77
,
051503
(
2008
).
26.
H. J. C.
Berendsen
,
D.
Vanderspoel
, and
R.
Vandrunen
,
Comput. Phys. Commun.
91
,
43
(
1995
).
27.
W. D.
Cornell
 et al,
J. Am. Chem. Soc.
117
,
5179
(
1995
).
28.
J.
de Andrade
,
E. S.
Boes
, and
H.
Stassen
,
J. Phys. Chem. B
106
,
13344
(
2002
).
29.
U.
Essmann
,
L.
Perera
,
M. L.
Berkowitz
,
T.
Darden
,
H.
Lee
, and
L. G.
Pedersen
,
J. Chem. Phys.
103
,
8577
(
1995
).
30.
P. E.
Ramírez-González
,
G.
Ren
,
G.
Saielli
, and
Y.
Wang
,
J. Phys. Chem. B
120
,
5678
(
2016
).
31.
O.
Borodin
,
J. Phys. Chem. B
113
,
11463
(
2009
).
32.
T.
Yan
,
Y.
Wang
, and
C.
Knox
,
J. Phys, Chem. B
114
,
6886
(
2010
).
33.
T.
Yan
,
Y.
Wang
, and
C.
Knox
,
J. Phys. Chem. B
114
,
6905
(
2010
).
34.
H.
Tanaka
,
J.
Meunier
, and
D.
Bonn
,
Phys. Rev. E
69
,
031404
(
2004
).