A molecular dynamics study of lithium-containing aprotic heterocyclic ionic liquid electrolytes

The steady rise in atmospheric CO₂ concentration and its relationship to climate change suggest the need to replace fossil fuels with renewable energy sources, such as wind and solar energy. These renewable energy sources require the development of cheap and reliable energy distribution and storage technologies.^1–3 

Marketed since the 1990s, lithium-ion batteries are now widely used for many different energy storage applications.^4,5 This is mainly due to their high density of charge compared with other commercial batteries. Lithium-ion batteries have been considered excellent candidates for large-scale applications, such as power plants and electrical vehicles. Current lithium batteries utilize organic electrolytes such as carbonates. This class of electrolyte presents two main limitations: the undesirable reactions with battery components that result in decreased battery lifetime and the risk of fire or explosion due to the high flammability of these solvents.^6–8 

One alternative to conventional electrolytes is the use of ionic liquids (ILs), which are salts that are liquids at low temperature, often room temperature or below. They are formed by the combination of typically asymmetric and/or bulky cations and anions. Many of them have favorable physicochemical properties,^9,10 such as good thermal and chemical stability, low vapor pressure, good conductivity, low flammability, and wide electrochemical window.¹¹ The combination of these physicochemical properties makes ILs excellent candidates for electrolytes in advanced batteries, but much is still unknown about their properties and performance.^6,12,13

Not surprisingly, there have been many studies over the past few years on IL electrolytes, and many of these have employed atomistic-level simulations. Such simulations enable researchers to compute several important macroscopic properties, such as conductivity and electrochemical window, while also providing molecular-level insight into the origins of these properties. The most commonly studied ILs are derived from the anions hexafluorophosphate ([PF₆]⁻), tetrafluoroborate ([BF₄]⁻), and bis(trifluoromethylsulfonyl)imide ([Tf2N]⁻) paired with imidazolium and pyrrolidinium cations.^14–19 

Zeiri and Dubnikova²⁰ also studied IL and Li⁺ mixtures. Based on density functional theory (DFT) calculations, they studied the influence of cations on different ILs containing [Tf2N]⁻ anions. They found that the small size and the large concentration of charge on the lithium ion create a strong Li⁺-anion interaction and a weakness in the IL cation-anion interaction and that the addition of a second [Tf2N]⁻ into the complex Li⁺-[Tf2N]⁻ leads to stabilization of the structure $Li+−([Tf2N]−)2$⁠.

Lesch and co-workers¹⁴ have studied the influence of IL cations on lithium ion coordination using molecular dynamics (MD) simulations with the APPLE&P (Atomistic Polarizable Potential for Liquids, Electrolytes, and Polymers) force field. Four mixtures of Li⁺ with two ILs, 1-ethyl,3-methylimidazolium bis(trifluoromethanesulfonyl) imide ([C₂mim][Tf2N]) and N-methyl-N-propylpyrrolidinium bis(trifluoromethanesulfonyl)imide ([pyr₁₃][Tf2N]), were studied. It was found that the strength of the cation-anion interaction has a significant influence on the lithium ion dynamics. In [pyr₁₃][Tf2N], the diffusivity of the lithium ion is low due to the high IL viscosity. However, this system has a relatively high Li⁺ transference number because of the weak anion-Li⁺ interaction and the subsequent short lifetime of the lithium-anion interaction.

Maginn and Liu studied the effect of the ion structure on Li⁺ mobility in IL/Li⁺ mixtures using MD simulations of 1-butyl,3-methylimidazolium bis(trifluoromethanesulfonyl) imide ([C₄mim][Tf2N]) and 1-butyl-4-methylpyridinium pyrrolide ([Bmpyr][pyl]) doped with Li⁺. It was found that the addition of lithium into the IL creates a decrease in ion mobility due the formation of a rigid structure Li⁺-(anion). Due to the planar shape of the [Bmpyr][pyl] and the parallel stacking between these ions, the decrease in dynamics of [Bmpyr][pyl] with increasing lithium concentration is less than that in [C₄mim][Tf2N].²¹ 

In the present work, we have chosen to investigate the behavior of Li⁺ in a class of ionic liquids containing aprotic heterocyclic anions (AHAs)²² paired with imidazolium and phosphonium cations.^23,24 This class of ILs is rarely studied for electrochemical applications, but the nature of the anions suggests that they might have interesting properties for Li⁺ transport.²¹ 

Brennecke et al.^23,25–28 conducted measurements of different physical and electrochemical properties of pure AHA ionic liquids as well as those doped with lithium. The pure ionic liquids were found to have high thermal stabilities, high decomposition temperatures, low or nonexistent melting points, wide electrochemical windows, low viscosities, high diffusivities, and reasonable conductivities. For the mixtures with lithium, it was found that the effect of the Li⁺ addition on the dynamical properties was smaller than that observed with [Tf2N]⁻. Li⁺ had no impact on thermal stability. Brennecke et al. also show alternatives to decrease the viscosity and increase the conductivity in these ILs, choosing planar cations or small phosphonium cations with or without ether groups. Imidazolium cations are of interest due to their ubiquity, while phosphonium cations have received much attention recently due to their favorable reductive stability, reversible Li⁺ deposition in IL/Li⁺ mixtures, and good diffusivity and conductivity relative to their viscosity.^29–32 

Understanding the properties of lithium electrolytes is challenging due to the numerous processes that can contribute to their properties, such as ion solvation, ion association, and numerous dynamic properties in the system. MD is a powerful technique to examine electrolyte properties because a direct correlation between the atomistic level information and macroscopic properties can be obtained.^33–35 MD simulations were carried out on 12 ILs based on the combinations of three different cations, 1-ethyl-3-methylimidazolium [C₂mim]⁺, (methyloxymethyl)triethylphosphonium [P222mom], and triethylbutylphosphonium [P2224]⁺, and four anions, 2-(cyano)pyrrolide [2-CNpyr]⁻, 3-(cyano)pyrrolide [3-CNpyr]⁻, 1,2,3-triazolide [123Triaz]⁻, and 1,2,4-triazolide [124Triaz]⁻.

Four different mole fractions of lithium in the IL were considered at 403.15 K. Dynamical properties were analyzed, and different combinations of the structure functions were calculated to understand how the anion and cation structures influence lithium dynamics.

The MD simulations were performed utilizing the following potential function, consistent with the General Amber Force Field (GAFF):³⁶ 

in which the distance between atom i and atom j is r_ij, the partial charges are q_i and q_j, and the nominal bond angles are θ₀. In this work, an all-atom model was used. Unlike non-bonded interactions were modeled using a Lorentz-Berthelot combining rule. For interactions between atoms separated by three bonds, van der Waals and Coulombic interactions were scaled by a factor of f = 0.5 and f = 0.833, respectively, while intermolecular interactions between atoms separated by fewer than three bonds were neglected.

The parameters for bonds, angles, dihedrals, and Lennard-Jones interactions for the ILs were taken directly from GAFF,^36,37 while the partial charges were obtained using the Restrained Electrostatic Potential (RESP) method³⁸ at B3LYP³⁹/6-311++G(d,p) level of theory.⁴⁰ The parameters for Li⁺ were taken from the OPLS (optimized potential for liquid simulations)⁴¹ force field.

The structures of the three cations and four anions that make up the 12 different ILs studied in this work are shown in Fig. 1, while Table I lists the different compositions and system sizes examined. For the calculation of the lithium molar fraction, the IL pair was considered as one unit as well as the Li-anion salt.

Each simulation was carried out using the same six step process. First, initial structures were generated using the PACKMOL⁴² package at the compositions shown in Table I. Next, a steepest descent energy minimization was carried out to relax high energy states. A 7 ns canonical ensemble (constant NVT) simulation at 700 K was then run to randomize the positions. After this, an annealing simulation was performed in the isothermal-isobaric (constant NpT) ensemble, in which the temperature was decreased from 700 K to the temperature of interest at 1 bar in two cycles of 4 ns each. Following this, an equilibration NpT simulation was run for 10 ns at 403.15 K and 1 bar. After these procedures, the system was deemed to be equilibrated; production runs of 12 ns in the NpT ensemble were carried out and statistics collected.

Temperature and pressure were maintained using the velocity-rescale⁴³ and Parrinello-Rahman⁴⁴ methods with a coupling time of 0.1 ps and 5 ps, respectively, except during the annealing simulations where the pressure was maintained by the Berendsen barostat⁴⁵ using a coupling time of 1.0 ps. Bond lengths were kept fixed at their nominal values using the LINCS algorithm,⁴⁶ and Coulombic interactions were handled with the smooth particle mesh Ewald method⁴⁷ with a direct space cutoff of 1.2 nm. A switching function was used to smoothly switch the Lennard-Jones potential to zero between 1.0 and 1.2 nm. The default Verlet neighbor list option was used to speed up the calculations. The leap-frog algorithm was used to integrate the equations of motion with a time step of 2 fs. All the simulations were performed using the Gromacs package, version 5.1.1.⁶⁵ All distribution functions were calculated using the Travis package,⁴⁸ while all other analyses were handled with Gromacs utility tools.

The calculated densities for the neat ionic liquids as a function of temperature are shown in Fig. 2. As expected, the densities all decrease linearly with increasing temperature. The densities of ILs with a common anion have the trend [C₂mim]⁺ > [P222mom]⁺ > [P2224]⁺. The ILs with triazolide anions have higher densities than those with cyanopyrrolide anions due to the smaller size of triazolide.²⁵ The computed densities of [C₂mim][123Triaz], [C₂mim][124Triaz], [C₂mim][2-CNpyr], and [P2224][2-CNpyr] are in good agreement with experimental data at 343 K,^23,24 differing by 2.37%, 2.79%, 3.20%, and 2.87%, respectively. No experimental data are available for the other ionic liquids, and no density data exist for any of the ILs with dissolved lithium.

Figure 3 shows the computed densities as a function of Li⁺ concentration at 403.15 K. The addition of lithium salt to these systems causes a nearly linear increase in density due to the strong interactions between Li⁺ and the anions.^19,25 All ILs doped with Li⁺ have a similar density increase when compared with the neat IL. This increase is less than 5% for the mixtures with a Li⁺ mole fraction of 0.25. This behavior is similar to that of other IL/Li⁺ mixtures reported in the literature.^16,49

The self-diffusivities of the ions were calculated using the Einstein relation^50,51

in which $r→i$ is the position vector of the center of mass for the ion i at a given time t.

The self-diffusion coefficients for neat [C₂mim][123Triaz], [C₂mim][124Triaz], [C₂mim][2-CNpyr], and [P2224][2-CNpyr] are consistent with experimental trends,^23,24 but the absolute values are lower than the experimental data by about a factor of ten. We are unaware of any experimental data for the other simulated ILs or for any of these ILs containing Li⁺.

The underestimation of dynamics in ILs is a well-known problem with fixed charged models, which have been shown to yield good structural and thermodynamic properties but give dynamics that are too slow due to a lack of charge transfer and polarization.^{14,15,52–54} Scaled charge models have been proposed to overcome this problem.^21,55,56 Because we are interested in liquid structure and relative trends in dynamics, we have opted to use a fixed charge force field with cation and anion charges of ±1 e. To test whether computed trends depended on the way charges were modeled, we carried out a series of simulations in which ion charges were uniformly scaled by a factor of 0.8. It was found that the main trends observed for the self-diffusivities were the same as those observed for the full charge models, although the overall dynamics were enhanced. Additional details may be found in the supplementary material.

Figure 4 shows the self-diffusion coefficients for the ILs containing the [P222mom]⁺ cation. As expected, the self-diffusion coefficients decrease with the increasing lithium concentration. It is well known that for IL/Li⁺ mixtures, there is a strong lithium-anion interaction, resulting in the formation of Li⁺-anion aggregates and a subsequent network between the ions, which is responsible for an increase in viscosity and decrease in ion mobility.^57–59 

For each IL in Fig. 4, the anion has a self-diffusivity that is higher than that of cation when no lithium is present. As the concentration of Li⁺ increases, however, there is a crossover in the self-diffusivities such that the [P222mom] cation has a larger self-diffusivity than the anions at highest Li⁺ concentration. This trend was also observed by Martins et al.²⁹ with similar phosphonium ILs. Using pulsed gradient spin-echo nuclear magnetic resonance (PGSE-NMR) and MD simulations, they found that the addition of lithium to (2-methoxyethyl)triethylphosphonium bis(trifluoromethanesulfonyl)imide ([P222(2O1)][Tf2N]) resulted in a similar crossover in self-diffusivity as was observed here. For triethylpentylphosphonium bis (trifluoromethanesulfonyl)imide ([P2225][Tf2N]), however, no such diffusivity crossover was observed. This difference in the behaviors of the two ILs was attributed to the presence of the alkyl ether chain in the [P222(2O1)]⁺ cation. The negative charge on the oxygen atom and the greater chain mobility created a weakening in the interactions between the ions, improving ion mobility. However, this weakening was not strong enough to break the Li⁺-anion interaction.²⁹ 

From an energy storage perspective, the dynamics of Li⁺ are particularly important. Li⁺ has the smallest self-diffusivity of all the ions in all mixtures due to its high charge density. Figure 4 shows that the addition of Li⁺ causes a decrease in Li⁺ self-diffusivity of 90% for the cyanopyrrolide-containing ILs with the [P222mom]⁺ cation. A similar decrease in dynamics is observed for cyanopyrrolide ILs having the [P2224]⁺ and [C₂mim]⁺ cation (see the supplementary material and Fig. 5). Interestingly, however, the dynamics of Li⁺ decrease much more modestly (about 50%) as a function of Li⁺ concentration for the triazolide-containing ILs with the [P222mom]⁺ cation. Roughly the same modest decrease in Li⁺ dynamics is observed for the triazolide ILs with the [P2224]⁺ cation, but the decrease is much larger for the triazolide ILs having the [C₂mim]⁺ cation. Although it appears that there is a plateau in the self-diffusion coefficient for the IL cation in the systems [P222mom][123Triaz] and [C₂mim][2-CNPyr] between Li⁺ mole fractions 0.10 and 0.18 (Figs. 4 and 5 respectively), there is actually a slight decrease with increasing Li⁺ concentration. This decrease is masked by the logarithmic scale and is greater than the computed uncertainty.

Figure 5 shows self-diffusivity results for the ILs containing the [C₂mim]⁺ cation and the four different anions. The [C₂mim]⁺ cation has a larger self-diffusivity than each anion, unlike the situation with the [P222mom]⁺ and [P2224]⁺ cations, where the cation and anion have approximately the same self-diffusivity. The fact that the imidazolium cation has a larger self-diffusivity than the anions is due to the planar shape of the cation, and it is consistent with previous observations.⁶⁰ Unlike the case with [P222mom]⁺, there is no crossover in the self-diffusivities of the cations and anions with the increasing lithium concentration.

The lithium ion transference number is a crucial property for a battery electrolyte and was calculated as

where the self-diffusivities were obtained from Eq. (2).

The lithium transference numbers in each IL are shown as a function of Li⁺ concentration in Fig. 6. The Li⁺ transference numbers generally increase with increasing lithium concentration over the entire concentration range. A few of the ILs, such as [C₂mim][123Triaz] and [P222mom][3-CNPyr], show a maximum in transference number at the intermediate Li⁺ concentrations. This is due to the competition between the number of Li⁺ ions in the solution and the increased viscosity associated with higher Li⁺ concentrations.⁴⁹ 

Two main trends can be seen in Fig. 6. First, Li⁺ has a higher transference number in ILs having a triazolide anion compared with those with a cyanopyrrolide anion. Second, the Li⁺ transference numbers in phosphonium triazolide ILs are higher than in [C₂mim]⁺ triazolide ILs. Both of these trends can be explained on the basis of the dynamics of the IL cations and anions. Triazolide anions have a stronger association with IL cations than do cyanopyrrolide anions, which leads to a higher viscosity and lower self-diffusivities of the IL ions. Similarly, the ILs with phosphonium cations have larger viscosities (and smaller ion self-diffusivities) than the [C₂mim] ILs because of the size and shape of the phosphonium cations. In both cases, the small self-diffusivities of the IL ions lead to larger Li⁺ transference numbers, consistent with Eq. (3). Similar results were observed by Lesch and co-workers,¹⁴ as described in the Introduction. They found a large transference number of Li⁺ in ILs with high viscosity. The strong interaction between cation and anion results in a larger number of free lithium ions in the mixture.¹⁴ However it is notable that the Li⁺ contribution to the overall mixture conductivity is small.⁵⁸ Overall we find that two main effects govern lithium diffusivity and transference number: the Li⁺-anion interaction and the IL viscosity.

The interaction between Li⁺ and the anions was analyzed in terms of partial radial distribution function g(r) between Li⁺ and the anion atoms. The partial g(r) between the Li⁺ and the anion nitrogens for the [P222mom]⁺ IL based systems is shown in Fig. 7. For [123Triaz]⁻, the first two peaks have similar distances at 0.23 nm because of the proximity between nitrogen atoms in the ring. The N1 atom has the highest intensity, however, due to the high concentration of negative charge on this atom (see Fig. S17 of the supplementary material). The coordination with [124Triaz]⁻ is established at 0.23 nm by the two ND atoms, while N1 atom has a small secondary contribution for the coordination at 0.42 nm.

Looking to the [2-CNPyr]⁻ g(r), we can see that the peaks for both nitrogens are localized at 0.23 nm. This suggests that due to the proximity between the nitrogens, lithium shows no preference for a single nitrogen but instead interacts equally with different nitrogens. For [3-CNPyr]⁻, on the other hand, there is a preference for Li⁺ to locate near the nitrogen of the ring and not the cyano nitrogen due to the positive charge on the carbon atom of the cyano group. These results are consistent with those observed in the center of mass g(r) shown in Fig. S4 of the supplementary material in which there is a sharp and intense peak for triazolide anions, while in the g(r) for cyanopyrrolide anions, there are multiple and broad peaks due to the presence of two different interaction sites in the anion structure (for more details, see the supplementary material). Looking at the g(r) for [C₂mim]⁺ and [P2224]⁺ ILs, we can see the same behavior observed for [P222mom]⁺ ILs, evidencing a minor influence of the IL cation in the lithium solvation. So, because of the similarity between them, all the structural analyses shown in this paper are for [P222mom]⁺ IL based systems; results for other ILs are provided in the supplementary material.

Due to the strong Li⁺-anion interaction, some anions are brought too close to each other by the lithium ions, and this effect can be seen in the anion-anion center of mass g(r). Figure 8 shows the g(r) of anion-anion interactions in the [P222mom]⁺-based IL systems. When no Li⁺ is present, there is a peak in each system at around 0.8 nm and a secondary peak at about 1.5 nm. These peaks are due to the normal charge ordering observed in ILs. When lithium is added to the ILs, a new peak appears at around 0.4-0.6 nm, depending on the IL. The peak intensity increases with increasing lithium concentration. This is due to the fact that Li⁺ is smaller than the IL cations and has the ability to coordinate to multiple anions. When it does so, it draws anions close to one another, which is why the peak at short distances appears. The intensity of the peak at around 0.8 nm decreases with increasing Li⁺ concentration because the number of anions coordinating with [P222mom]⁺ decreases. This is consistent with previous reports of the formation of Li⁺ |anion and anion|anion aggregates in the liquid.^16,17,21,49 For triazolide anions, the peak formed by the lithium addition has a structure much more well defined than cyanopyrrolide. Cation-cation and cation-anion g(r) plots show few differences between the anions and the Li⁺ concentrations, showing that the presence of Li⁺ affects just the anion; figures are provided in the supplementary material.

In order to characterize the anion|anion aggregates and the lithium effect on the liquid structure, the combined distribution functions (CDFs) were computed, which is the combination of the radial distribution function [g(r)] with an angle distribution function (adf). The α is defined by the angle between the normal ring vector and the vector joining the two geometric ring centers in the anions, which characterizes the distance vector r. The angle and distance used in the CDF for the ILs formed by [P222mom]⁺ cation and triazolide anions are shown in Fig. 9. The left and right panels are the neat IL and the 0.25 Li⁺ fraction, respectively; both neat ILs have similar structures, with a weak tendency for planar stacking (high probability of α angles close to 0° and 180°) at 0.75 nm. The effect of lithium in the systems is quite different between the anions. The addition of Li⁺ in [123Triaz]⁻ induces a strong tendency for parallel stacking at 0.4 nm while the populations at 0.75 nm vanished. The parallel stacking is due to the presence of lithium in regions above and below the anion ring. For [124Triaz], there is a combination of parallel and perpendicular stacking at distances between 0.47 and 0.58 nm due to the positions of the nitrogens in the ring. The Li⁺-anion coordination is mainly established by the ND atoms with the lithium ion in the plane of the ring. If the anions adopt a combination of parallel and perpendicular stacking, however, it is possible to maximize the number of interactions between anions and lithium ions. For more information about the orientation of lithium around anions, see Figs. S17 and S22 and the discussion in the supplementary material.

Figure 10 shows the effect of Li⁺ on the cyanopyrrolide anions. For both neat ionic liquids, there are perpendicular and parallel stacking at 0.65 and 0.80 nm, respectively; when lithium is added into the [2-Cnpyr]⁻, there is a slight increase in the perpendicular stacking at 0.6 nm and the parallel stacking populations increase between 0.40 and 0.60 nm. The parallel stacking at 0.80 nm disappears due to the shift of the anion to close distances around Li⁺. The mixture Li⁺/[3-CNPyr]⁻ shows a much more simple structure than [2-Cnpyr]⁻. The populations in the neat IL are gone (as well as the g(r) peaks), and new broad parallel populations appear at 0.475 nm. The difference observed in different cyanopyrrolide anion systems is due to the different contributions of their nitrogens to the lithium coordination.

The formation of Li⁺–Li⁺ aggregates in IL has been reported in the literature via both experimental^61,62 and MD simulations.^15,16,49,63 These aggregates result in a reduction of Li⁺ mobility and solubility due to the stabilization of the structures and saturation in the coordination sites. Haskins and co-workers found the formation of Li⁺–Li⁺ aggregates in three different ionic liquids, [Pyr14][Tf2N], [Pyr13][FSI] (1-propylmethylpyrrolidinium bis(fluorosulfonyl)imide), and [C₂mim][BF₄] (1-ethyl-3-methylimidazolium tetrafluoroborate). They claim that the aggregates are the result of the small distance between the Li⁺ ions due to the anion size and the fact that lithium ions share anions between them. They claim that the presence of well defined peaks under 1 nm in the Li⁺–Li⁺ g(r) are the evidence of Li⁺ aggregates in the liquid structure.⁴⁹ 

In the current work, lithium aggregates were analyzed in terms of Li⁺–Li⁺ g(r) following the definition of Haskins et al.⁴⁹ For all g(r) shown in Fig. 11, there are well defined peaks below 1 nm. However, the structure of Li⁺–Li⁺ aggregates and the influence of lithium concentration are quite different between the two classes of anions. For triazolide ILs, there are intense and sharp peaks below 0.5 nm (two peaks for [123Triaz]⁻ and one peak for [124Triaz]⁻). For cyanopyrrolide ILs, however, there are several split peaks below 1 nm. The sharp peaks close to 0.25 nm in the triazolide anions correspond to a close and homogeneous interaction between the metal ions due to the small size of the anion, while the broad and multiple peaks in the cyanopyrrolide ILs are due to the presence of two interaction sites in the anion and relatively large size of the cyanopyrrolide anions. It is interesting to highlight here that the presence of sharp and intense peaks at distances below 0.34 nm is unusual for the most common anions studied in IL electrolytes.

It is known that the dynamical properties may often be related to free volume (FV) effects. The free volume (FV) was obtained for all simulated systems following the methodology described elsewhere.⁶⁴ The FV was computed by placing a point probe at random locations in the simulation box. This was repeated 10⁶ times, and the present FV was taken as the fraction of times the probe did not overlap with any atoms in the system. Each atom was assumed to be a hard sphere whose size was given by its Lennard-Jones diameter, σ. Figure 12 shows the percentage of free volume for [P222mom]⁺ based systems as a function of the Li⁺ mole fraction. The FV decreases in all cases as Li⁺ is added due to the structuring induced by Li⁺. The change in the FV is most pronounced for the cyanopyrrolide anions. These results are consistent with those observed in self-diffusivities and transference numbers.

In this work, the effect of lithium addition on thermodynamic properties and solvation structure in different ionic liquids was studied. Molecular dynamics simulations for ILs based on the combinations of three different cations and four aprotic heterocyclic anions (AHAs) doped with lithium at different molar fractions were performed at 403.15 K. Diffusion coefficients of all ions were calculated, and diverse structural analyses were performed to characterize the effect of lithium addition into the IL. The influence of IL ion system on the liquid dynamics and packing structure in the mixture was also discussed.

The structural analysis shows that the Li⁺-anion interaction leads to a new structure in the ILs. Due to the small size of lithium, multiple anions can have close contact with a Li⁺, resulting in the formation of various |anion|Li⁺|anion aggregates. By comparing the liquid structure in different ILs studied in this work, it was found that the lithium solvation environment depends just on the anion structures. The IL cation has negligible influence on the lithium solvation structure or the aggregate structure. On the other hand, lithium does not exert much influence on the cation-anion structure packing.

The addition of Li⁺ to the ILs also affects the dynamics of the ions. There is a decrease in the mobility of all species with the increase of Li⁺ molar fraction due to the strong Li⁺-anion interaction and the formation of aggregates in the liquid structure. For the [C₂mim]⁺ based ILs, it was found that the cations have faster diffusivity than the anions at all Li⁺ mole fractions, while in the phosphonium based ILs, some ILs exhibit a crossover between cation and anion diffusivity due to the significant reduction in the mobility of the anions because of Li⁺ anion association. In all simulated systems, the Li⁺ was found to have the lowest diffusion coefficient due to the high charge density.

The dynamics in phosphonium cyanopyrrolide ILs was affected more by the increase in Li⁺ molar fraction than that in phosphonium triazolide ILs. This difference is also seen in the Li⁺ transference number. A larger transference number was found for the triazolide systems, evidencing a larger number of free lithium ions in these systems than in the cyanopyrrolide-based ILs.

This study has shown that the dynamics in lithium containing ILs is complex and hard to predict. The tendency of Li⁺ to complex with the anions and reduce IL cation-anion interactions is key to determining self-diffusivities and transference numbers. The addition of lithium to the ILs lowers the free volume, which correlates with the reduction in dynamics.

The choice of an ionic liquid electrolyte requires the investigation of many physicochemical properties, and it is usually necessary to find a balance between them. In this work, the phosphonium triazolide based ionic liquids showed the best potential for electrolyte applications due to their good lithium transference number and the small influence of Li⁺ concentration on the system dynamics.

See supplementary material for the diffusion coefficients obtained from the scaled charge models; the center of mass g(r) for the interactions Li⁺-anion, anion-anion, cation-cation, and cation-anion for all systems simulated in this work at 403.15 K; the partial g(r) for the interaction Li⁺-anion in the [C₂mim]⁺ and [P2224]⁺ based ILs; the CDF for the interaction Li⁺-anion in the [C₂mim]⁺ ILs; and the CDFs for the interaction anion-anion in the [C₂mim]⁺ and [P2224]⁺ based ILs.

We thank FAPERJ, especially by the process 201.995/2016, and CAPES for the fellowship given to TCL, LAME-UFF, and Notre Dame Center for Research Computing for the computing support. Y.Z. and E.J.M. were supported by the U.S. Department of Energy, Basic Energy Science, Joint Center for Energy Storage Research under Contract No. DE-AC0206CH11357. L.T.C. acknowledges the CNPq fellowship and FAPERJ, JCNE No. 214996/E06/2015.