Molecular insights into an ether-functionalised ionic liquid electrolyte with hydrogen-modified anions at electrode interfaces

Abstract

The chemistry of cations and anions in an ionic liquid (IL) significantly affects its bulk phase and interfacial properties, requiring an in-depth understanding to advance its application in energy storage, including batteries and supercapacitors. In this work, a novel ionic liquid electrolyte featuring an ether-functionalised cation, N-methyl-N-methoxyethylpyrrolidinium [C_2O1mpyr], and a hydrogen-modified anion, (difluoromethanesulfonyl) (trifluoromethanesulfonyl)imide [DFTFSI], was studied computationally as a sodium battery electrolyte. The impact of the ether oxygen in the cation and hydrogen in the anion on the bulk phase structure, ion diffusion and interfacial chemistry was examined. The presence of the ether oxygen enhances ion diffusion and disrupts Na-anion aggregation through Na–O_cation coordination in both the bulk phase and at the electrode interface. The Na–O_cation coordination slightly affects the number of cations and anions in the innermost electrolyte layer but significantly influences the redox stability of the cation. The decomposition of the DFTFSI anion was explored, and two reduction routes and potential products were identified. This work provides new insights into the role of functionalised cations and anions in ionic liquid electrolytes, contributing to the development of advanced sodium batteries.

---

Introduction

Lithium-ion batteries (LIBs) are currently the dominant energy storage technology, extensively used in the automotive and electronics industries. However, concerns regarding the safety, cost, and finite supply of raw materials for LIB production have spurred interest in developing alternative battery technologies. Sodium metal batteries (SMBs) have attracted significant interest because of the abundance of sodium, its chemical similarity to lithium, and the potential for higher specific capacity compared to LIBs when using sodium metal as an anode.¹ However, the commercialisation of sodium batteries faces significant challenges, such as limited cycling stability, dendrite formation, and volumetric changes in the Na metal electrode.² Addressing these issues requires the development of compatible electrolytes with improved thermal and electrochemical stability, which is crucial for enabling high-voltage operation and enhancing the overall performance of SMBs.³

Ionic liquids (ILs) have emerged as promising electrolytes for lithium and sodium metal batteries due to their low volatility, good ionic conductivity, non-flammability, and high thermal stability.⁴ These liquid salts, composed of organic cations paired with organic or inorganic anions, can be tailored for specific properties by selecting different ion combinations or introducing chemical modifications, making them highly versatile electrolytes for energy storage devices.^5,6

Modifications in cations and anions affect the properties and bulk phase behaviour of ILs, changing ion dynamics, coordination, and various physical characteristics.^7–9 For instance, introducing an ether functional group into the cation side chain has been shown to enhance the conductivity of IL electrolytes.¹⁰ Chen et al. attributed this to the increased flexibility of the alkoxy side chain, hypothesizing that it reduces packing efficiency and increases the free volume.¹¹ However, Welton et al. later pointed out that the ether oxygen orients towards the organic cation charge centre, leading to an increase in density. This coordination also effectively shields the cationic charge and improves ion diffusion.¹² Other studies have shown that the alkoxy side chains can also interrupt ion aggregation, further improving ion diffusion.^9,13

Anion modification has also been explored to tune electrolyte properties. Zhang et al. introduced hydrogen atoms into the widely used bis(trifluoromethanesulfonyl)imide (TFSI) anion, adding functionalities and polarity, which in turn affect salt solubility, viscosity, and other electrolyte properties.^14,15 For example, H-substituents facilitate the formation of H-bonds in polymer electrolytes, thereby slowing down anion motion and increasing the Li transference number. However, increasing the number of H atoms in TFSI was found to raise the Li desolvation energy and reduce the anion's redox stability. The variant with a single hydrogen substitution on one side, named DFTFSI, struck the best balance by improving the Li transference number without significantly compromising conductivity and redox stability.¹⁶ Additionally, improved interfacial stability against metallic lithium (Li⁰) has also been observed with the DFTFSI anion, which was attributed to the possible formation of beneficial solid electrolyte interphase (SEI) products, such as LiH and LiF,¹⁴ although further investigation is needed.

The IL/electrode interfaces critically affect the electrochemical performance of IL electrolytes.^17,18 Cations, anions and salts influence the interfacial structure and chemistry of the IL, i.e., the number and arrangement of electrolyte species at the electrified interface.^19,20 This, in turn, directly impacts the early formation of the SEI layer, a key factor in electrochemical performance.

Classical molecular dynamics (MD) simulation is an effective tool for providing molecular-level insight into the ionic liquid interface near the charged electrode surface before SEI formation.²¹ Complementing computational analysis with experimental characterisation has become a powerful approach for studying the electrolyte interface and guiding interface design through a preconditioning process to enhance the electrochemical performance of IL electrolytes.^19,22 In addition, ab initio molecular dynamics (AIMD) has been used to study the decomposition of simple electrolyte components on an anode surface, primarily focusing on organic solvents and anions.^23,24 A few studies have also investigated IL electrolytes, revealing that cations remain stable within a few picoseconds (ps) while anions undergo decomposition.^25,26 However, research on electrolyte decomposition on sodium metal surfaces remains limited despite its importance for sodium battery applications.

In this study, we conducted a comprehensive computational investigation of the IL [C_2O1mpyr] [DFTFSI] (Fig. 1) with a sodium salt. To assess the impact of the ether oxygen in the cation side chain, we compared this electrolyte with [C₄mpyr] [DFTFSI]. The neat ILs were first simulated, followed by the introduction of NaDFTFSI at a low concentration of 0.5 mol kg⁻¹, considering the possible low solubility of this sodium salt in these ILs. We compared their bulk phase structures and ion diffusion, and analysed their interfacial chemistry near different neutral and negatively charged metal surfaces. Finally, the decomposition of the DFTFSI anion was studied using the AIMD method to understand the effect of the introduced H atom. This work provides an in-depth understanding of this new IL electrolyte.

Fig. 1 The structures of the (a) C₄mpyr cation, (b) C_2O1mpyr cation, (c) DFTFSI anion and (d) Na employed in the study.

Results and discussion

Bulk phase study of IL electrolytes

The bulk phase properties of both neat ILs and their salt systems were first investigated. The C_2O1mpyr-containing electrolyte systems exhibit higher densities than their C₄mpyr counterparts (Table S1†), regardless of salt addition. Adding alkali salt further increases the density of both ILs. This finding aligns with the work of Welton et al.,¹² who reported an approximately 8% increase in the density of ether-functionalised ILs. In this study, we observed a 5.3% increase in density for neat ILs and a 4.5% increase for salt-containing ILs when comparing the C_2O1mpyr to the C₄mpyr system. Welton attributed this effect to a reduction in hydrodynamic radius, as the oxygen atom orients toward the central charged nitrogen, leading to more compact packing in ether-functionalized ILs. Here, we analyse the ion–ion coordination and packing through radial distribution functions (RDFs). The coordination number (CN) was determined at the first significant minimum of the RDF and is provided in Table S2.†

Fig. 2(a) and (b) present the RDFs of cation (N₄)–anion (NBT) (through nitrogen atoms) and Na-anion (NBT), which are nearly overlapped between the two systems, indicating a negligible difference in their coordination distances. The difference in the CN of cation–anion and Na–anion between the two ILs is minor, with C_2O1mpyr showing a slightly higher cation–anion CN (6.9 vs. 6.7) and C₄mpyr exhibiting a slightly higher Na-anion CN (4.1 vs. 4.0).

Fig. 2 RDFs calculated for (a) cation(N4)–anion(NBT) through nitrogen atoms, (b) Na–anion (NBT), (c) cation(N₄)–cation(N4), (d) N₄–O in C_2O1mpyr and N₄–C in C₄mpyr in which C is the second last carbon in butyl chain, (e) Na–C1 and (f) Na–C2 (C1 and C2 are carbon atoms in CF₂H and CF₃ groups in DFTFSI), (g) Na–O (ether oxygen in C_2O1mpyr), and (h) Na–Na.

However, a slight difference is observed in the first peak of the cation–cation RDF (Fig. 2(c)), which shifts to the left for C_2O1mpyr cations, suggesting their closer packing compared to the C₄mpyr cations. This observation is further supported by the spatial distribution function (SDF) analysis of cations surrounding a reference cation, as shown in Fig. S1.† Within the same cutoff distance, the SDF map for C_2O1mpyr displays a denser distribution compared to C₄mpyr. This closer arrangement between C_2O1mpyr cations is likely due to weaker electrostatic interactions between the ether oxygen and the cationic nitrogen centre, which are absent in C₄mpyr. This is supported by the N₄–O and N₄–C₈ RDFs in Fig. 2(d). A small N₄–O RDF peak appears near 4.5 Å, which is absent in the N₄–C₈ RDF, where C₈ refers to the second last carbon in the butyl side chain.

In addition, the first prominent peak in Fig. 2(d) indicates a different orientation of the ether oxygen in C_2O1mpyr compared to the equivalent C₈ site in C₄mpyr. The left shift of the N–O peak indicates that the oxygen is oriented more toward the charge centre, further supporting Welton's argument. Therefore, the higher density observed in the C_2O1mpyr system could be attributed to both the side chain orientation and the tighter cation packing.

Next, the Na coordination to the asymmetric anion was examined by calculating RDFs between Na and N and between Na and two C atoms in the CF₂H and CF₃ groups. Overall, the bi-dentate Na–DFTFSI coordination dominates as reflected by a sharp Na–N peak at around 4.5 Å in Fig. 2(b). Fig. 2(e) and (f) present the Na–C RDFs for the two ILs. In both cases, the Na–CF₂H coordination number is slightly higher, with values of 3.35 for C_2O1mpyr and 3.65 for C₄mpyr, compared to the Na–CF₃ coordination number, which is 3.22 for C_2O1mpyr and 3.29 for C₄mpyr within a cutoff distance of 6.5 Å. This is supported by the earlier NMR results for the same ILs, but with a LiDFTFSI salt.⁷

In addition, we investigated the coordination between Na and the ether oxygen in the cation side chain using the Na–O RDF (Fig. 2(g)). A prominent peak is observed near 2.6 Å, indicating their coordination, with a CN of 0.4 (Table S2†). The Na–Na RDF and CN profiles were generated, shown in Fig. 2(h), to provide insight into Na distribution. The multiple RDF peaks indicate structured ordering, suggesting Na-involved local nanostructures due to Na–anion aggregation. The higher CN in the C₄mpyr system reflects a denser Na distribution, possibly linked to the larger aggregate size.

Na–anion aggregates form via anion-bridging coordination among Na ions and can therefore be visualized by connecting Na ions that fall within the cutoff distance of the first Na–Na RDF valley (∼10 Å). Fig. 3(c) and (d) illustrate these aggregates and clearly show the smaller aggregates in the C_2O1mpyr system. A complete distribution of cluster sizes was obtained by counting the frequency of each cluster size observed throughout the MD trajectory (Fig. S2†). Overall, smaller Na cluster sizes appear more frequently in the C_2O1mpyt system, whereas clusters with five or more Na ions are more often observed in the C₄mpyr system. The smaller aggregate size in the C_2O1mpyr system is attributed to Na–O_C2O1mpyr coordination, which involves approximately 18.8% of Na ions. This coordination interrupts Na–DFTFSI coordination and disrupts the formation of larger ion aggregates, which is expected to have a positive effect on ionic conductivity.

Fig. 3 MSD of different ions in the (a) neat and (b) salt IL systems. The snapshots illustrate the Na–DFTFSI aggregates by connecting Na ions in the same aggregate for the (c) C_2O1mpyr and (d) C₄mpyr systems.

Ion diffusion was studied using the mean square displacement (MSD), where a higher MSD value indicates faster ion diffusion. Fig. 3(a) and (b) suggest higher diffusion of all ions in the C_2O1mpyr system, both for the neat IL and the salt-containing systems. This is consistent with the reported lower viscosity and higher conductivity for [C_2O1mpyr][TFSI] compared to [C₄mpyr][TFSI] in the literature,¹¹ although the fully fluorinated TFSI anion was used in those studies. In both ILs, the diffusion of Na is lower than that of the organic cations and anions. This behaviour is commonly observed in ionic liquids and is primarily attributed to the stronger electrostatic interactions between metal ions and anions compared to those between IL cations and anions.²⁷

Interfacial study of IL electrolytes

The electrode/electrolyte interface plays a crucial role in determining the overall performance of the system. Here, we investigated IL electrolytes near the gold surfaces with different constant charge densities of 0, ±6, and ±16 μC cm⁻² in the fixed charge model (FCM). A comparison is also made with the results obtained using the more advanced constant potential electrode model (CPM) for the salt systems. The applied potentials are 0, 2.5 and 7.0 V, which correspond to average charge densities of 0, −5.9 and −16 μC cm⁻², respectively, for the negative electrode surface. These numbers are comparable with the charge densities used in the FCM. A more detailed comparison between the two electrode models in terms of derived average charged densities and potential differences is shown in Fig. S3 and Table S3.†

The number density profiles are generated for organic cations (via N4), anions (via NBT), and Na ions for both neat ILs and salt systems (Fig. 4). Multiple peaks are observed near uncharged or negatively charged gold surfaces, corresponding to multiple interfacial electrolyte layers, which have been reported previously through atomic force microscopy studies.^28–30 The neat ILs in Fig. 4(a) and (b) exhibit approximately three to four overlapping number density peaks of cations and anions near the uncharged surface (0C). Although the fourth peak is relatively minor, the overlap of the peaks indicates that these electrolyte layers are composed of both ionic species. When the surface is negatively charged, the first anion peak decreases with an increase in the height of the second peak due to the exclusion of anions from the innermost electrolyte layer. This results in more pronounced alternating cation/anion layers. The anion peak nearly vanishes at −16 μC cm⁻², accompanied by a significant growth in the cation number density peak due to the increased number of cations required to screen the more negative surface charge.

Fig. 4 The number density profiles of ions for neat C_2O1mpyrDFTFSI and C₄mpyrDFTFSI ILs (a) and (b) and their salt systems (c)–(f) near the negatively charged gold surface using the constant charge model (a)–(d) and the constant potential model (e) and (f). 1C = 1 μC cm⁻².

The addition of Na salt alters the composition of the innermost layer, especially near negative surfaces, as shown in Fig. 4(c) and (d). A small Na number density peak emerges at approximately 0.5 nm in the C₄mpyr system and grows with increasing negative surface charge densities for both ILs. The reduction in the first anion peak is less pronounced in the salt system, due to the enhanced stabilisation provided by Na–anion coordination. With the constant potential model, shown in Fig. 5(e) and (f), a similar trend was observed for the number density peak of each type of ions, although the small Na peak was also observed this time in the C_2O1mpyr electrolyte at 0 V.

Fig. 5 Average number of cations, anions, and Na ions in the inner-most electrolyte layer adsorbed on the different charged surfaces for (a) neat IL and salt-containing systems from the (b) FCM and (c) CPM. 1C = 1 μC cm⁻².

A quantitative analysis of ion adsorption on the gold surface at different surface charge densities is presented in Fig. 5. For both neat ILs (Fig. 6(a)), the number of adsorbed cations increases with increased negative charge density, while the number of anions decreases, becoming absent at −16 μC cm⁻². In the salt systems (Fig. 6(b) and (c)), a similar trend is observed, but the decrease in anions is less pronounced due to the concurrent increase in Na adsorption. The presence of Na helps stabilize anions at the interface. At −16 μC cm⁻² for the FCM or −7 V for the CPM, a slightly higher number of C_2O1mpyr and lower numbers of DFTFSI and Na ions are adsorbed on the gold surface compared to the C₄mpyr system. This consistency suggests that the more economical FCM approach reliably captures the main effects of increased current densities or voltages on the interfacial chemistry of these electrolytes, making it a cost-effective option for such studies. Since cations and anions contribute to the organic and inorganic SEI components, respectively, a moderate current density may help maintain a better balance between the two ion species compared to a high current density in this case.

Fig. 6 Snapshots of the innermost electrolyte layer adjacent to uncharged or negatively charged surfaces. (a) Fixed Charge Model (FCM): (i) C_2O1mpyrDFTFSI, (ii) C₄mpyrDFTFSI, (iii) C_2O1mpyrDFTFSI + NaDFTFSI, and (iv) C₄mpyrDFTFSI + NaDFTFSI. (b) Constant Potential Model (CPM): (i) C_2O1mpyrDFTFSI + NaDFTFSI and (ii) C₄mpyrDFTFSI + NaDFTFSI. The colour bar represents charge (e) fluctuations on the electrode surface.

Fig. 6 displays snapshots of the innermost electrolyte layer adsorbed on the Au surface under different applied surface charges or voltages. The neat ILs shown in Fig. 6(a)(i, ii) exhibit the same ion adsorption trends discussed earlier. In the salt-containing systems, a noticeable difference is seen in Na coordination, which occurs with both anions and cations through ether oxygens in C_2O1mpyr, but only with anions in the C₄mpyr system. This Na–C_2O1mpyr coordination enhances the number of cation-coordinated anions at low surface charge densities and interrupts the growth of Na–anion aggregates at high surface charge densities, as indicated by both FCM and CPM results. On the other hand, this difference in Na coordination is expected to affect the Na desolvation energy, which refers to the free energy required for Na to partially or completely shed its solvation shell in order to approach the anode surface and undergo reduction. Consequently, this may impact overpotentials during cycling and is worth further investigation in combination with experimental validation.

The Na–O (C_2O1mpyr) coordination also affects the redox behaviour. Table 1 suggests that the redox stability of the C_2O1mpyr cation is greater than that of C₄mpyr, with a less negative LUMO, a more negative HOMO and a wider electrochemical window. However, when Na is coordinated to C_2O1mpyr, the oxidation stability is further improved, but the reduction stability decreases as the LUMO orbital shifts to a more negative potential.

Table 1 HOMO, LUMO, and HOMO–LUMO gap calculated using Gaussian 09 at the MP2/6-311++G(d,p) level of theory

Cation	HOMO (eV)	LUMO (eV)	HOMO–LUMO gap (eV)
C2O1mpyr	−0.55951	−0.07123	0.48828
C2O1mpyrNa	−0.72036	−0.22665	0.49371
C4mpyr	−0.45438	−0.09721	0.35717
DFTFSI	−8.359	4.134	12.495
TFSI	−8.605	4.779	13.384

---

In terms of ion orientation on the surface, both cations and anions mainly lie on the surface, with cation side chains attached to the surface, regardless of the surface charge density. Cation rings align mostly parallel to the neutral and low-charged surfaces, as suggested by the angular analysis in Fig. S4,† which shows the probability distribution of the angle between the cation ring and the gold surface. A cation ring oriented parallel to the surface corresponds to an angle in the range of 0–30° and 150–180°. At −16 μC cm⁻², the percentage of tilted rings increases for the C_2O1mpyr cation but decreases for the C₄mpyr cation, especially in the salt system. The side chains of both cations align mostly parallel to the gold surface.

In addition, the snapshots in Fig. 6 indicate that the arrangement of the C_2O1mpyr cation on the gold surface is also affected by inter-cation interactions. The O–N interaction promotes alignment between adjacent C_2O1mpyr cations, in contrast to C₄mpyr where non-polar side-chain aggregation leads to microphase-separated polar and non-polar regions. This tighter cation–cation interaction also supports a slightly higher number of C_2O1mpyr adsorbed on the gold surface compared to the C₄mpyr.

AIMD study of anion decomposition on the Na electrode surface

The redox stability of DFTFSI was compared to that of TFSI through DFT calculations, as shown in Table 1. The introduction of H slightly decreases the redox stability of the anion, narrowing the electrochemical window from 13.384 eV to 12.495 eV, due to a decreased LUMO orbital and an increased HOMO orbital.

We further investigated the decomposition behaviour of the H-containing anion DFTFSI on a sodium surface through AIMD simulations. Given the high simulation cost, we conducted a 9 ps simulation to explore the early-stage decomposition processes. Four initial configurations of Na–DFTFSI ion pairs were placed on the Na (100) surface, as shown in Fig. S5.† These configurations were extracted from the MD interface simulations and include one bidentate Na–DFTFSI ion pair (Configuration A) and three monodentate Na–DFTFSI ion pairs: one where Na is present near the CF₃ group with the H-fragment attached (Configurations B), one where Na is unattached to the anode surface (Configuration C), and one where Na is present next to the CF₂H side (Configuration D). Snapshots illustrating the decomposition processes are presented in Fig. 7 alongside the corresponding reaction equations. The time required to reach the transition state varies among configurations. For each case, we first present the structure just 1 ps prior to the initial bond cleavage, followed by key snapshots capturing subsequent major bond-breaking events.

Fig. 7 The snapshots from AIMD simulations show the decomposition process of the DFTFSI anion on the Na(100) surface with four initial configurations of A, B, C, and D. Atom types are indicated.

At the transition state, all four configurations exhibit oxygen atoms attached to the sodium surface, with both the CF₃ and CF₂H groups oriented upward. Configurations A and B, originating from different bi-dentate and monodentate initial structures, evolve into similar transition state geometries. Configuration A has four oxygen atoms and Configuration B has three oxygen atoms attached to the surface, with the Na cation lying on the surface, coordinated to the SO₂CF₃ group through one oxygen in A and two oxygen atoms in B. Configuration C also has four-oxygen attached to the sodium surface, similar to A, but the Na cation is positioned above the surface and coordinated to the CF₃ group. In Configuration D, two oxygen atoms are attached to the surface, and the Na ion is also located above the surface but coordinated to the CF₂H group.

Interestingly, both Configurations A and C have four oxygen atoms attached to the surface and follow the same decomposition pathway, despite the different positions of the Na cation. This observation may suggest that the decomposition of DFTFSI is affected mainly by its initial binding structure on the surface. The S–C bond cleavage occurs first, forming –SO₂CF₃ and –NSO₂CF₂H fragments at 1.27 ps and 2.49 ps for Configurations A and C, respectively. The nitrogen-containing NSO₂CF₂H fragment remains stable throughout the remainder of the simulation. In contrast, the –SO₂CF₃ fragment undergoes further decomposition, breaking down into SO₂ and CF₃. The SO₂ molecule remains intact, while the CF₃ group subsequently decomposes by gradually losing its fluorine (F) atoms one by one. The final decomposition products of Configuration A are NaF, SO₂, and –NSO₂CF₂H. For Configuration B, only the first F-cleavage event is observed within the simulation timeframe.

Configurations D results in the same decomposition products as Configuration A, but through a slightly different pathway in the first and second decomposition steps. This divergence is likely associated with its different initial structure, where only two oxygen atoms are attached to the surface before decomposition. The decomposition begins with the S–C bond cleavage, forming a CF₃ group and a larger –NSO₂HCF₃SO₂ fragment, both of which are unstable. The SO₂ cleavage occurs almost immediately (after 0.01 ps), producing the more stable –NSO₂CF₂H fragment. The CF₃ group then undergoes the same decomposition pathway, losing F atoms one by one.

Finally, a different decomposition pathway and products are observed for Configuration B. The S–C bond cleavage happens first at 4.72 ps on the SO₂ group, which has both of its oxygen attached to the surface. This produces –NSO₂CF₃ and –SO₂CF₂H fragments, with the N attaching to the SO₂CF₃ moiety in this case. This nitrogen-containing fragment remains stable throughout the remaining time. The CF₂HSO₂ fragment undergoes further decomposition, losing two F atoms at 5.36 ps and 6.00 ps, respectively. Next, one oxygen atom detaches from the –SO₂CH fragment at 6.06 ps, and the remaining –SOCH fragment remains intact. Different decomposition products were observed here, including NaF, NaO, SOCH and NSO₂CF₃ fragments.

Notably, the stability of N-containing fragments like NSO₂CF₃ has also been reported in other studies investigating the decomposition of the TFSI anion on the Li electrode surface, where AIMD simulations were run up to 230 ps.^23,26 Additionally, we did not observe further decomposition of the SOCH fragment or the formation of NaH, which may require further investigation.

Conclusions

This research investigated an ether-functionalised ionic liquid, C_2O1mpyr DFTFSI, with a hydrogen-modified anion to understand how these chemical modifications affect its electrolyte behaviour with a sodium salt. Both bulk phase and interfacial structures were studied through molecular dynamics simulations and compared with those of its counterpart C₄mpyrDFTFSI.

Bulk-phase analysis reveals that while the C_2O1mpyr IL and its salt system exhibit higher densities, they also show greater ion diffusivities compared to C₄mpyrDFTFSI. The higher density could be attributed to the orientation of the ether oxygen in the side chain towards its nitrogen charge centre. Additionally, the weaker inter-cationic interaction between the ether oxygen and nitrogen charge centre also contributes to a more tightly packed cation–cation distance. The coordination between Na⁺ and C_2O1mpyr via the ether oxygen leads to a reduction in NaDFTFSI aggregates. Na⁺ coordinates with DFTFSI, an asymmetrical anion, in both bidentate and monodentate modes, with monodentate coordination to the CF₂H side being slightly more prevalent than to the CF₃ side.

Interfacial studies were performed near a gold electrode with increasing negative surface charge using two theoretical electrode models. Both ILs show similar trends of increased cations and Na⁺, along with decreased anions, as the surface charge becomes more negative. While the quantitative differences in the numbers of cations, anions, and Na⁺ on the negative gold surface between the two ILs are minor, their coordination structures differ significantly when sodium salt is present. In the C_2O1mpyr system, Na⁺ interacts with the ether oxygen, leading to fewer anions and smaller Na–DFTFSI aggregates on the surface compared to the C₄mpyr system. This Na–O coordination also reduces the reductive stability of the C_2O1mpyr cation while enhancing its oxidative stability.

AIMD simulations were conducted to explore the early-stage decomposition of the DFTFSI anion on a Na metal surface. Two decomposition pathways were identified: one initiated by S–N bond cleavage and the other initiated by S–C bond cleavage. The nitrogen-containing fragments remain stable throughout the 9 ps simulation, while the remaining part undergoes a quick decomposition. One route results in the formation of NaF, SO₂, and NSO₂CF₂H, and the second route produces NaF, NaO, SOCH and NSO₂CF₃ fragments. No NaH was observed within the simulation timeframe, likely due to the stability of C–H bond. This work provides molecular-level insight into this novel IL with modified cation and anion chemistries, increasing our understanding of how these chemistries can be used for IL electrolyte development in Na batteries.

Computation methodology

Molecular dynamics force field and validation

The force field parameters (bond, angle, torsion angle, Lennard-Jones parameters, and charges) of the C_2O1mpyr and C₄mpyr cations, sodium and the DFTFSI anion were generated using the OPLS_AA force field via the online generator LigParGen.³¹ A charge scaling factor of 0.7 was adopted to account for polarisation effects.³² The force field was validated by comparing the simulated densities and ion diffusivities with experimental results,⁸ as shown in Table S1.† A small error of 1.4% was observed for density, which is considered acceptable. Ion diffusivities were calculated using the mean square displacement (MSD) through the Einstein relation. They also show good agreement with the pulsed field gradient nuclear magnetic resonance (NMR) measurements for both C_2O1mpyr and DFTFSI. Although the MD model slightly underestimates diffusivity, which is a known limitation of classical MD simulations using a non-polarisable force field, it accurately predicts the D(C_2O1mpyr)/D(DFTFSI) ratio of 1.25, consistent with the NMR results. This suggests that the MD model effectively captures the relative diffusivity difference between the cation and anion, thereby indicating the reliability of the model. The force field parameters for the gold electrode were adopted from the literature.

MD simulation details

The bulk phase investigations were carried out using GROMACS 2022.3 version.³³ PACKMOL³⁴ was utilized to generate the initial system, consisting of 216 randomly packed ion pairs. Coulombic and van der Waals interactions were calculated within a 1.2 nm cutoff. Long-range electrostatic interactions were calculated using Particle Mesh Ewald (PME) and long-range van der Waals interactions were accounted for by using energy dispersion correction. Bonded hydrogens were constrained using LINCS, and a time step of 1 fs was used, and the integrator was leapfrog. The initial system was subjected to energy minimization using the steepest descent minimization (STEEP) algorithm. Following minimization, an annealing process was performed by raising the temperature to 700 K and then cooling it down to 393 K in steps, using the NPT ensemble. The annealing process lasted 18 ns, with at least 9 ns of equilibration at 700 K. V-scale temperature coupling and Berendsen pressure coupling were applied to ensure proper mixing of cations and anions while enabling efficient system optimization. Once the annealing was completed, the system underwent further equilibration sequentially at a target temperature of 293 K, using the NPT ensemble for 30 ns. Nose–Hoover temperature coupling and Parrinello–Rahman pressure coupling were adopted for NPT equilibrium calculations. Lastly, the NVT simulation was conducted for 20 ns for ion diffusion investigation.

For the interface model, the electrolyte was confined between two face-centred Au (111) electrodes, with the z-distance determined by the bulk phase electrolyte density. The Au (111) electrode dimensions were 4.48 × 4.48 nm. A vacuum region, twice the size of the electrode separation, was added to separate the electrodes and eliminate artifact interactions arising from the opposite sides of the charged electrodes due to the use of periodic boundary conditions (PBCs).³⁵

The fixed charge model (FCM) was adopted in GROMACS to prepare the electrified electrodes, with positive or negative charge evenly distributed among the gold atoms on the surfaces adjacent to the electrolyte. The electrodes were fixed in the X, Y and Zdirections throughout the simulation using freezegrps and freezedims options. Electrostatics were treated using the Particle Mesh Ewald (PME) method with Ewald-geometry = 3dc setting, which applies a slab correction to reduce spurious interactions along the Z direction. The initial system was subjected to a STEEP minimization, followed by an annealing–cooling procedure similar to that used in bulk simulations, using the NVT ensemble. A final NVT production run was performed for 30 ns using Nose–Hoover temperature coupling. Bonded hydrogens and both short-range and long-range interactions were considered in a similar manner to the bulk simulations. A velocity-Verlet integrator was used here.

The interface simulation with the Constant Potential Model (CPM), which allows for a more realistic charge oscillation in response to the electrolyte, was conducted using the LAMMPS 28 March 2023 version with the ELECTRODE package.^36,37 Short coulombic and van der Waals interactions were calculated within a 1.2 nm cut-off, similar to the previous setup. Long-range electrostatic and van der Waals interactions were accounted for with the PPPM (Particle–Particle Particle-Mesh) with slab 3.0 and tail methods, respectively. Bonded hydrogens were constrained using SHAKE, and the electrodes were fixed by neutralizing the force at each time step using the fix setforce option. A time step of 1 fs was used. Initially, the system was minimized using the conjugate gradient method. Subsequently, the system was annealed–cooled using the same protocol described previously for the bulk and FCM models. Finally, an NVT production run was performed for 8 ns. Both Gromacs and Lammps output trajectories were saved in XTC format, enabling the use of Gromacs analysis tools, such as RDF and MSD.

Ab initio molecular dynamics (AIMD)

AIMD simulations were conducted to investigate the reduction reactions of the DFTFSI anion on the sodium surface, using the CP2K software package.³⁸ A Na–DFTFSI ion pair was placed on the Na(100) metal surface, with different initial geometries suggested by MD simulations. The simulation box dimensions were set to 11.5 × 11.5 × 31.0 Å to ensure sufficient vacuum space and avoid interactions. The Becke-Lee-Yang-Parr (BLYP) DFT functional^39,40 was adopted for geometry optimization, with Grimme's DFT-D3 dispersion correction to account for van der Waals interactions.⁴¹ For the AIMD simulations, the DZVP-MOLOPT-SR-GTH basis set⁴² was used, along with Goedecker−Teter–Hutter (GTH) pseudopotentials⁴³ to represent the core electrons. A real-space grid cutoff of 800 Ry was applied for the electron density representation within the GPW (Gaussian and Plane Wave) framework, ensuring convergence of the total energy during geometry optimization and AIMD simulations. For system preparation, a geometry optimization (GO) was conducted first, where all atoms were allowed to relax, and the Na atoms were not frozen during the GO step to allow them to fully relax. After completing the GO step, the system was subjected to AIMD simulations. During the AIMD simulations, the bottom three layers of Na slabs were frozen. The system was treated under the NVT ensemble, with a Nose–Hoover thermostat maintaining the temperature at 350 K. A time step of 1 fs was used, and the total simulation time was extended to 9 ps. The system was simulated under periodic boundary conditions in the xy-plane, with a vacuum region in the z-direction to prevent interactions between periodic images. The simulations were conducted under the NVT ensemble, with a Nose–Hoover thermostat maintaining the system at 350 K. Lastly, the orbital energy gaps were calculated using the Gaussian 09 (G09) software package. Molecular structures were optimized using density functional theory (DFT) with the B3LYP functional and the 6-311++G(d,p) basis set.⁴⁴ The HOMO and LUMO energies was then calculated using Møller–Plesset second-order perturbation theory (MP2) with the 6-311++G(d,p) basis set in the gas phase.

Data availability

All the data that support the findings of this study are available from the corresponding authors upon reasonable request.

Author contributions

F. C. conceived the idea, and led and supervised the research. S. B. R. conducted the research and drafted the manuscript with F. C. J. S. provided guidance on simulations. A. P., L. A. O., M. A., P. H., and M. F. contributed to scientific discussions on the project and the results. All authors contributed to the manuscript revisions.

Conflicts of interest

The authors declare no conflict of interest.

Acknowledgements

All authors acknowledge the Australian Research Council for financial support through DP210101172. This work was supported by the Australasian Leadership Computing Grants scheme, with computational resources provided by NCI Australia, an NCRIS-enabled capability supported by the Australian Government.

References

1. J.-Y. Hwang, S.-T. Myung and Y.-K. Sun, Chem. Soc. Rev., 2017, 46, 3529–3614.
2. S. Wang, B. Peng, J. Lu, Y. Jie, X. Li, Y. Pan, Y. Han, R. Cao, D. Xu and S. Jiao, Chem.–Eur. J., 2023, 29, e202202380.
3. Q. Li, J. Chen, L. Fan, X. Kong and Y. Lu, Green Energy Environ., 2016, 1, 18–42.
4. A. Balducci, in Ionic Liquids II, ed. B. Kirchner and E. Perlt, Springer International Publishing, Cham, 2018, pp. 1–27,  DOI:10.1007/978-3-319-89794-3_1.
5. M. Palluzzi, A. Tsurumaki, H. Adenusi, M. A. Navarra and S. Passerini, Energy Mater., 2023, 3, 300049.
6. A. Lewandowski and A. Świderska-Mocek, J. Power Sources, 2009, 194, 601–609.
7. D. Gyabeng, L. Qiao, H. Zhang, U. Oteo, M. Armand, M. Forsyth, F. Chen and L. A. O'Dell, J. Mol. Liq., 2021, 327, 114879.
8. D. Gyabeng, L. Qiao, H. Zhang, U. Oteo, M. Armand, M. Forsyth, F. Chen and L. A. O'Dell, J. Phys. Chem. C, 2021, 125, 14818–14826.
9. M. Kunze, E. Paillard, S. Jeong, G. B. Appetecchi, M. Schönhoff, M. Winter and S. Passerini, J. Phys. Chem. C, 2011, 115, 19431–19436.
10. L. J. A. Siqueira and M. C. C. Ribeiro, J. Phys. Chem. B, 2009, 113, 1074–1079.
11. Z. J. Chen, T. Xue and J.-M. Lee, RSC Adv., 2012, 2, 10564–10574.
12. F. Philippi, D. Rauber, B. Kuttich, T. Kraus, C. W. Kay, R. Hempelmann, P. A. Hunt and T. Welton, Phys. Chem. Chem. Phys., 2020, 22, 23038–23056.
13. L. C. Branco, J. N. Rosa, J. J. Moura Ramos and C. A. Afonso, Chem.–Eur. J., 2002, 8, 3671–3677.
14. H. Zhang, U. Oteo, X. Judez, G. G. Eshetu, M. Martinez-Ibanez, J. Carrasco, C. Li and M. Armand, Joule, 2019, 3, 1689–1702.
15. H. Zhang, U. Oteo, H. Zhu, X. Judez, M. Martinez-Ibanez, I. Aldalur, E. Sanchez-Diez, C. Li, J. Carrasco and M. Forsyth, Angew. Chem., 2019, 131, 7911–7916.
16. U. Oteo, M. Martinez-Ibanez, I. Aldalur, E. Sanchez-Diez, J. Carrasco, M. Armand and H. Zhang, Chemelectrochem, 2019, 6, 1019–1022.
17. D. A. Rakov, F. Chen, S. A. Ferdousi, H. Li, T. Pathirana, A. N. Simonov, P. C. Howlett, R. Atkin and M. Forsyth, Nat. Mater., 2020, 19, 1096–1101.
18. M. V. Fedorov and A. A. Kornyshev, Chem. Rev., 2014, 114, 2978–3036.
19. D. Rakov, M. Hasanpoor, A. Baskin, J. W. Lawson, F. Chen, P. V. Cherepanov, A. N. Simonov, P. C. Howlett and M. Forsyth, Chem. Mater., 2021, 34, 165–177.
20. D. A. Rakov, J. Sun, S. A. Ferdousi, P. C. Howlett, A. N. Simonov, F. Chen and M. Forsyth, ACS Mater. Lett., 2022, 4, 1984–1990.
21. S. Begic, E. Jonsson, F. Chen and M. Forsyth, Phys. Chem. Chem. Phys., 2017, 19, 30010–30020.
22. U. Pal, D. Rakov, B. Lu, B. Sayahpour, F. Chen, B. Roy, D. R. MacFarlane, M. Armand, P. C. Howlett and Y. S. Meng, Energy Environ. Sci., 2022, 15, 1907–1919.
23. H. S. Dhattarwal, Y.-W. Chen, J.-L. Kuo and H. K. Kashyap, J. Phys. Chem. C, 2020, 124, 27495–27502.
24. D. Stottmeister and A. Grob, Batteries Supercaps, 2023, 6, e202300156.
25. J. Clarke-Hannaford, M. Breedon, T. Rüther and M. J. S. Spencer, ACS Appl. Energy Mater., 2020, 3, 5497–5509.
26. H. Yildirim, J. B. Haskins, C. W. Bauschlicher Jr and J. W. Lawson, J. Phys. Chem. C, 2017, 121, 28214–28234.
27. O. Borodin, G. A. Giffin, A. Moretti, J. B. Haskins, J. W. Lawson, W. A. Henderson and S. Passerini, J. Phys. Chem. C, 2018, 122, 20108–20121.
28. D. S. Silvester, R. Jamil, S. Doblinger, Y. Zhang, R. Atkin and H. Li, J. Phys. Chem. C, 2021, 125, 13707–13720.
29. X. Mao, P. Brown, C. Červinka, G. Hazell, H. Li, Y. Ren, D. Chen, R. Atkin, J. Eastoe, I. Grillo, A. A. H. Padua, M. F. Costa Gomes and T. A. Hatton, Nat. Mater., 2019, 18, 1350–1357.
30. J. M. Black, M. Zhu, P. Zhang, R. R. Unocic, D. Guo, M. B. Okatan, S. Dai, P. T. Cummings, S. V. Kalinin, G. Feng and N. Balke, Sci. Rep., 2016, 6, 32389.
31. L. S. Dodda, I. Cabeza de Vaca, J. Tirado-Rives and W. L. Jorgensen, Nucleic Acids Res., 2017, 45, W331–W336.
32. C. Schröder, Phys. Chem. Chem. Phys., 2012, 14, 3089–3102.
33. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
34. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
35. I.-C. Yeh and M. L. Berkowitz, J. Chem. Phys., 1999, 111, 3155–3162.
36. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
37. L. J. Ahrens-Iwers, M. Janssen, S. R. Tee and R. H. Meibner, J. Chem. Phys., 2022, 157.
38. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 15–25.
39. A. D. Becke, Phys. Rev., 1988, 38, 3098.
40. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B:Condens. Matter Mater. Phys., 1988, 37, 785.
41. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132.
42. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127.
43. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B, 1996, 54, 1703.
44. M. Frisch, F. Clemente, Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino and G. Zhe, Gaussian, 2009, vol. 9.

---

Footnote

† Electronic supplementary information (ESI) available: Density and diffusivity comparison between experimental and simulation data, density comparison between neat and salt systems, coordination numbers (CNs) and cutoff used for CN calculations between cations and anions, potential difference vs. potential of zero charge, cation arrangement in the bulk phase and its spatial distribution function, potential plots, angular distribution of cations on the electrode surface, initial structures of AIMD calculations and AIMD videos in the videos folder. See DOI: https://doi.org/10.1039/d5ta02471d

---