Molecular-level heterogeneity in deep eutectic electrolytes

Abstract

Deep eutectic electrolytes (DEEs) have in recent years gained momentum as prospective compounds for next-generation (lithium) batteries. Yet, the connection between molecular-level properties and macroscopic performance of DEEs is poorly understood. We have here, by molecular dynamics simulations, studied in detail a few simple DEEs created using the hydrogen bond (HB) donor N-methyl-acetamide (NMA) and one of three different common lithium salts: LiBF₄, LiDFOB, and LiBOB, all in a 1 : 4 molar ratio composition. By altered anion size and symmetry, we observe varying molecular-level heterogeneity (MLH), reflecting differences in local structure and coordination as well as dynamic disorder, and especially the heterogeneity of the HB network. The HB network becomes more localized the more asymmetric and larger the anion is, which also affects the ion-pairing and aggregation of the DEEs. From a dynamic point of view, DEEs with larger MLH show slower ion self-diffusion, due to steric hindrances caused by the larger anions and the localized HB network. Overall, both the MLH and the HB network, inherent properties of DEEs, must be properly balanced to allow for adequate ion transport.

---

Introduction

The growing need for renewable energy integration and transport electrification requires efficient energy storage, i.e., battery technologies.¹ The electrolyte plays a crucial role for any battery technology, and while liquid electrolytes have been extensively researched for lithium-ion batteries (LIBs), they are limited by e.g., the thermal stabilities of the solvents and salts employed.^2,3 Deep eutectic electrolytes (DEEs) have emerged as a promising alternative to conventional liquid electrolytes and have been used both for LIBs and lithium metal batteries (LMB).^2,4–8 DEEs-in this context-are eutectic mixtures of pure compounds, often based on the sole addition of a hydrogen bond donor (HBD) to a lithium salt.^2,9,10 DEEs offer several advantages, such as ease of synthesis, tunability, wide electrochemical stability windows, and low vapor pressures.^2,11

DEEs have, for example, been applied to create LMBs using lithium bis(trifluoromethanesulfonyl) imide (LiTFSI) in combination with 2,2′-dipyridyl disulfide (DpyDS) as HBD, showing a strong ion–dipole interaction promoting enhanced ion conductivity.¹² Additionally, Hu et al.¹³ developed a DEE based on LiTFSI and acetamide, claiming exceptional physicochemical properties with an ionic conductivity exceeding 1 mS cm⁻¹. We here, however, focus on three different lithium salts, whereof lithium difluoro(oxalato)borate (LiDFOB) based electrolytes render improved ionic conductivities as compared to lithium tetrafluoroborate (LiBF₄) based ones, and in addition, a better-working low-temperature working range for LIBs.^14–16 Our third salt of choice, lithium bis(oxalato)borate (LiBOB), has gained interest mainly due to its fluorine-free design, lower cost, and higher thermal stability as compared to the for LIBs commonly used LiPF₆ salt.^17,18

However, the main reason as to why we here have chosen these three salts is that they represent a variation in anion symmetry and size (Table 1). The BF₄⁻ anion, the highest in symmetry, is pseudospherical, while formally belonging to the T_d point group, and is also the smallest. Moving to bulkier and lower in symmetry anions, BOB is elongated and has rings rich in oxygen sites, in stark contrast to BF₄⁻ with no such sites at all. DFOB naturally presents in-between characteristics.

Table 1 Anion structure, symmetry, and size

Anion	Structure	Symmetry	# Symmetry elements	Size^a
a Length of longest axis = maximum atomic distance.
BF4^−		T d	20	∼2.3 Å
DFOB		D 2h	8	∼4.3 Å
BOB		D 2d	7	∼6.6 Å

---

Overall, these differences in symmetry, size, and chemical composition render different charge distributions, as visualized by their electrostatic potentials (ESPs) (Fig. 1). This strategy thus allows us to systematically monitor effects on DEE properties.

Fig. 1 ESPs for: (a) BF₄⁻, (b) DFOB, (c) BOB, (d) NMA.

As the HBD, we have chosen N-methyl acetamide (NMA), for which, e.g. the work by Boisset et al.,⁵ shows that it acts as a complexing agent for both anions and cations and can thus potentially disrupt the cation–anion interactions. This functionality arises from its two polar groups—the NHR group, which coordinates with anions, and the CO group, which coordinates with cations–visible by its ESP (Fig. 1).

Bang and Yeu¹⁹ also highlighted the influence of the composition, i.e., the Li-salt : NMA molar ratio, on the ionic and electronic conductivity. Their DEEs exhibited increased ionic conductivity with decreased salt concentration when moving from molar ratios 1 : 1 to 1 : 5. While there is no fixed ratio that guarantees the eutectic point, the optimal molar ratio for ionic conductivity was between 1 : 4 to 1 : 5.

The more general class of deep eutectic solvents (DESs), which form the basis of DEEs, have been extensively studied for a wide range of applications,^20,21 including as battery electrolytes, and their fundamental characteristics and behaviors carry more or less directly over into DEEs.^4–8 DESs have recently been reported to have a specific molecular ordering and a radially layered structure,^22,23 which is suggested to arise due to a combination of interactions, such as Coulombic, hydrogen bonding (HB), and van der Waals interactions, etc.,^22–24 and has also been described in terms of spatial heterogeneity on the molecular length scale.²⁰ Several of these insights have been obtained using MD simulations, highlighting the role of HB networks and structural organization in DESs systems.^20,22–25 One study described what was referred to as micro-heterogeneity (MH) as the separation of polar and nonpolar molecular parts, revealing that the HB network of the polar domains could be indirectly affected by the MH.²⁵ A complex HB network seems to be a main characteristic of DESs,^20,23–29 and we hypothesize here that this is the case also for DEEs. However, DEEs remain less studied from a molecular simulation perspective, and the specific role of anion symmetry and size in shaping heterogeneity has not yet been systematically addressed.

Thus far, however, the relationship between the DEE molecular-level structure and the local geometry/organization and symmetry and size of the ingoing DEE components is not well understood. We suggest that a deeper understanding of how, e.g., the anion symmetry and size influence the HB network, which in turn affects the structural organization and macroscopic behavior, could lead to DEEs with enhanced physicochemical and electrochemical properties. We therefore here introduce molecular-level heterogeneity (MLH), defined mainly by the heterogeneous distribution of the HB network arising from variations in vital molecular properties such as size, charge, shape, and manifests as coordination and solvation distribution differences, and also in dynamic properties, such as ion diffusion rates. Hence the MLH causes non-uniform behaviour that ultimately impacts the overall performance of the DEEs.

To investigate how primarily anion symmetry and size affect the HB network and thus, in turn, cause MLH, we as a proof-of-concept study use three DEEs, with components defined above, at a 1 : 4 molar ratio by molecular dynamics (MD) simulations.

Methodology

xTB MD simulations

xTB simulations^30,31 were implemented using CP2K³² with three-dimensional periodic boundary conditions using the geometry-, frequency-, noncovalent- and extended tight-binding (GFN-xTB) method.³¹ This semi-empirical quantum MD (xTB-MD) approach incorporates an augmented basis set with improved HB treatment and was specifically developed to yield accurate geometries, vibrational properties, and noncovalent interactions. It offers a favorable compromise between the efficiency of classical MD and the accuracy of DFT, enabling us to simulate the complex ion–solvent dynamics in DEEs.³¹

The initial configurations for the three DEEs were generated using CHAMPION,³³ which relaxes the system using the conjugate gradient method to avoid atom overlaps and refines the geometry output. The 1 : 4 LiBF₄/LiDFOB/LiBOB:NMA DEEs were composed of 486, 464, and 496 atoms, respectively (see SI for breakdown in molecules), and had cubic box side lengths of 16.7 Å, 17.3 Å and 16.7 Å, respectively, yielding initial densities within 0.98–1.1 g cm⁻³, consistent with previously reported values for similar systems. The total energy for the DEEs was obtained as:

E = E_el + E_rep + E_disp + E_XB (1)

where the total energy is represented by the electronic (E_el), atom-pairwise repulsion (E_rep), dispersion (E_disp), and halogen-bonding (E_XB) terms.^31,34 For each DEE, an initial MD equilibration run (40 ps) was performed using an NPT ensemble with 1 fs time-steps at 363 K. The temperature was chosen to ensure that the DEE system remains liquid and to facilitate configurational sampling. Equilibration was assessed by monitoring the convergence of potential energy and the stability of temperature, volume, and RDFs, indicating that structural and thermodynamic properties had stabilized. The barostat set the pressure to 1 atm and a CSVR was used as a thermostat. The xTB method was used in both equilibration and production runs and once the energy and the density had converged, and evaluating the total energy as a function of simulation time, MD production runs (60 ps) were made in NVT ensemble with 1 fs time-steps and a Nosé–Hoover thermostat³⁵ set to 363 K. Following this, the generated trajectories were analyzed using VMD,³⁶ and the (local) structure studied by calculating the g(r) to generate radial distribution functions (RDFs) as:where V represents the total volume, r is the distance between a pair of particles, p(r) is the average number of atom pairs found at a distance between r and r + dr and N_pairs is the number of unique pairs of atoms where one atom is from each of two sets selections.³⁷

The coordination number (CN) was subsequently calculated by integrating the g(r), over spherical coordinates, up to the first minimum of the RDF. Partial coordination numbers (pCNs) were similarly determined by integrating the partial RDFs, and also yield the cumulative distribution functions (CDFs) for specific pairs of atoms. The anion solvation numbers (SNs) were obtained by summing the pCN contributions of anion–solvent and anion–cation. In addition, we used Travis³⁸ to generate spatial distribution functions (SDFs), which provide a three-dimensional representation.

To analyze the HB distribution, hydrogen bonds were first identified using geometric cutoffs of HBA–HBD distance ≤3.5 Å and HBA–H–D angle ≥150°. We then employed a Python script that from the MD simulation trajectories calculates the probability of finding an HB within a specific radius, based on HBD–HBA distances. The probability distribution function (PDF) was generated by binning these distances and normalizing the data to reflect the probability of HB formation. Additionally, we characterized the HB strength and spatial distribution using a contour plot through Travis,³⁸ where the x-axis represents the HBD–HBA distances (e.g., O–H, F–H) and the y-axis the HB angles. To quantify “HB strength” and enable direct comparison across the three DEEs, we differentiated the cumulative distribution functions of these distance–angle bins to obtain a 2D probability density P(r,θ). All heat maps share a fixed color scale (v_min = 0; v_max = global maximum P over all DEEs), so identical color intensities correspond to identical absolute HB probabilities with warmer colors indicating more probable (stronger) HB geometries and cooler colors less frequent.

Classical MD simulations

MD simulations were performed to study dynamic properties and potential heterogeneity, as already roughly observed from the xTB simulations. The dynamic analysis is based on the calculation of self-diffusion coefficients from the mean-square displacements (MSDs).

The MD simulations were performed in LAMMPS,³⁹ using the non-polarizable OPLSAA force field.⁴⁰ To account for electronic screening effects and enhance the accuracy of inter-ionic interaction predictions, the partial charges of ions were scaled by a factor of 0.8.^41,42 The simulation boxes contained 200 cations, 200 anions, and 800 NMA molecules. For each box, topology files, bonded and Lennard-Jones parameters, and initial configurations were generated using Packmol⁴³ and fftool.⁴⁴ The simulations consisted of an energy minimization step employing conjugate gradients, followed by equilibration in the microcanonical ensemble (NVE) for 500 ps, succeeded by another 1 ns in the isothermal–isobaric ensemble (NPT) at 1 atm, and an additional 1 ns in the canonical ensemble (NVT), all carried out at 363 K. These initial runs were designated as equilibration runs and were excluded from the subsequent analysis. Production runs were performed in the NVT ensemble for 14 ns at 363 K and 1 atm, using the average simulation box size from the NPT equilibration run. A Nosé–Hoover thermostat was used with a temperature damping of 100 fs and a pressure damping of 1000 fs. Electrostatic interactions were computed using the particle–particle–particle–mesh scheme (PPPM), and periodic boundary conditions were applied in all directions. MSDs were calculated using the compute msd subroutine in LAMMPS, which tracks the displacement of atoms over time and averages the squared distance between their current and reference positions.⁴⁵ The self-diffusion coefficients were then calculated using a self-written script in MATLAB, determining the self-diffusion coefficient, D, via the Einstein relation:

where r_i(t) is the position of the i-th particle at the time t, and d is the dimensionality of the system (e.g., d = 3 for three-dimensional diffusion). The brackets 〈·〉 denote the ensemble average, which is often approximated by a time average over the duration of the simulation. In practice, reasonable averages require sufficiently long simulation durations to ensure that the individual ion movements become uncorrelated and that the dynamics as monitored by the MSD reach the diffusive regime.⁴⁶

Molecular-level heterogeneity

To quantify how primarily the HB network, but also other factors, causes MLH in the DEEs, an index is here proposed as a metric. This way, the MLH across the different DEEs can be compared, facilitating a deeper understanding. The index is designed to capture both structural and dynamic MLH, and is formulated as a weighted sum of selected features:where f_i represents individual molecular features, and w_i are their corresponding weights, constrained such that . The features incorporated include the HB localization factor, the diffusion ratio (D_anion/D_NMA), the anion size, the anion SN, and symmetry elements. All features are normalized by their maximum value prior to PCA to prevent any single feature from dominating due to scale. This ensures fair weighting regardless of units. Features inversely proportional to heterogeneity (e.g., diffusion ratio, symmetry) are transformed via,aligning their interpretation with MHI directionality. Testing using z-score standardization showed preserved trends, with minor shifts in absolute weights. The final correlation-based weight refinement further corrects for any over- or under-weighting caused by scaling. The HB localization factor, derived from the Gini coefficient, which is a statistical measure of dispersion,⁴⁷ is used here to quantify the HB distribution.

To compute the index, we employ a two-step hybrid approach combining Principal Component Analysis (PCA) and Correlation Analysis to determine the weights ω_i. PCA is first used to assess which features contribute most to the variance across the DEEs. However, since PCA emphasizes statistical variance rather than physical relevance, we refine the weights through correlation analysis by computing the correlation between each feature and the initial MLH estimate. This hybrid method ensures that the final weights reflect both statistical significance and physicochemical meaning, aligning the index with heterogeneity-driving behaviors in the system.

The general trend of MHI values was preserved under alternative weighting schemes (e.g., uniform or PCA-only weights), but the hybrid approach improved numerical resolution and better reflected qualitative HB and mobility observations.

Results and discussion

We here present a detailed analysis of the local structure using the RDFs, the CNs, the anion SNs, and HB networks, derived from our xTB–MD trajectories to study various ion–ion and ion–solvent interactions (Fig. 2). Next, we analyze the dynamic behavior through MSDs and self-diffusion coefficients using classical MD simulations. Finally, we conclude with the MLH and connect it back to the structure and dynamics of the DEEs.

Fig. 2 Schematic of the interactions: cation–anion, anion–solvent, and cation–solvent.

Local structure

We begin our analysis with the RDFs (Fig. 3), alongside the pCNs (Fig. 4) and SDFs (Fig. 5), to assess the local structure in the DEEs. We do this to scrutinize the roles of the different interactions present in the DEEs: cation–anion, cation–solvent, and anion–solvent interactions (Fig. 2).

Fig. 3 RDFs for cation–anion interactions (a)–(c), cation–solvent interactions (d), (e), and anion–solvent interactions (f) for all the DEEs.

Fig. 4 The three pCNs for each of the DEEs.

Fig. 5 Snapshots after the production run of the: (a) LiBF₄:NMA, (b) LiDFOB:NMA, and (c) LiBOB:NMA DEEs, respectively. Red is lithium cation, orange is anion, and blue is solvent.

From the RDFs (Fig. 3(a)–(c)), we observe that the lithium cation primarily interacts with the oxygen atoms of the DFOB and the BOB anions, with strong Li–O peaks at approximately 2 Å, indicating highly localized first cation coordination shells. The highest cation–anion pCN is observed in the LiBOB-based DEE (2.5), followed by the LiDFOB-based DEE and then the LiBF₄ one (1.4) (Fig. 4). This trend follows the anion size and available coordination sites: BOB > DFOB > BF₄⁻.

The increased ion–ion interactions reflect the size of the anion.⁴⁸

In more detail, in the LiDFOB-based DEE, the lithium cations tend to coordinate predominantly with the oxygen atoms (82%) rather than the fluorine atoms (18%), considering that a lack of preference would result in 66% and 33%. This selective interaction suggests that the presence of multiple coordination sites within DFOB facilitates more oxygen-centric solvation.

Turning to the cation–solvent interactions, the LiBF₄-based DEE shows the highest lithium cation solvation by NMA (pCN = 2.7), with a sharp Li–O(NMA) peak at 1.9 Å (Fig. 3(e)), reinforcing the role of solvent structure. The cation–solvent solvation decreases in the order: LiBF₄ > LiDFOB > LiBOB based DEEs, which hence inversely correlates with cation–anion coordination, and indeed the lower cation–solvent coordination in the LiBOB-based DEE we attribute to the stronger and dominant cation–anion interactions that limit the lithium cations’ ability to interact with NMA. Given that the total Li⁺ CN of ≈4.2 in the LiBF₄ DEE, we note its agreement with the tetrahedral coordination found experimentally by Seo et al.⁴⁹ in LiBF₄·(AN)₂ solvates comprising two acetonitrile nitrogen atoms and two BF₄ fluorine atoms, indicating that Li⁺ favors a tetrahedral arrangement of its nearest neighbours.

Finally, for anion–solvent interactions, the B–N RDF (Fig. 3(f)) indicates that BF₄⁻, with its smaller size and higher symmetry, has the highest anion–solvent coordination (pCN = 3.0), followed by DFOB and then BOB (pCN = 0.7). The isotropic nature of BF₄⁻ facilitates uniform solvation, whereas the steric hindrance of BOB restricts solvent interactions, making it less solvated in comparison. The solvation trends overall suggest that the LiBOB-based DEE exhibits stronger ionic association, forming a distinct heterogeneous lithium–anion domain with the solvent at a longer distance, i.e., mostly outside the first cation solvation shell, leading to enhanced ion-pairing.

The main features in the B–B anion–anion RDF are somewhat closer and stronger, and also binodal for the LiBOB-based DEE at 4.1 Å and 5.0 Å (Fig. 6(a)).

Fig. 6 RDFs of: (a) anion–anion, and (b) cation–cation correlations for the three DEEs.

This could suggest some kind of MH order, where these anions are localized together, something that is indeed visible in the snapshot (Fig. 5(c)).

Such binodal features are often associated with intermediate-range structuring or domain-like aggregation, consistent with MLH.

For the Li–Li cation–cation RDFs, there are no similar differences-all DEEs have the same main feature at approx. 3.5 Å (Fig. 6(b)). However, very noticeable is the much less ordered feature beyond this distance for the LiDFOB-based DEE, telling the story of the different cation solvation and cation–anion interactions in that particular DEE. These deviations in local ordering across salts support the manifestation of MLH, where local structural disorder varies based on anion symmetry and size.

And then, using the SDFs (Fig. 7) to visualize if there is a clear correlation between the choice of anion and the spatial distribution of ligands around them, we find that the LiBF₄-based DEE, which exhibits the highest anion and cation solvation, also has the SDF showing the most “closed” and spatially uniform isosurface, thus following anion size and symmetry. Specifically, this enables more NMA nitrogen atoms to approach and interact with the anion, creating a rather homogenously structured DEE altogether.

Fig. 7 SDFs for the: (a) LiBF4, (b) LiDFOB, and (c) LiBOB based DEEs showing lithium in purple and NMA in yellow around the anions.

In stark contrast, the SDFs for the DFOB and BOB-based DEEs display more localized and fragmented distributions, suggesting greater anisotropy and heterogeneity.

The cation coordination is more structured, especially for the LiBOB-based DEE, where distinct solvation “pockets” rather than a continuous shell are seen, due to the directional coordination preference introduced by the BOB anion.

HB analysis

In general, fluorine atoms tend to play a dominant role in HB formation in the DEEs, with F–H interactions being more prevalent than O–H interactions, due to their stronger electronegativity and spatial freedom. And indeed, in the LiDFOB-based DEE exhibits the most HBs due to its combination of fluorine and oxygen HBA sites, with pCN values of F–H = 1.3 and O–H = 0.8 (Fig. 8).

Fig. 8 RDFs of the HBs in the three DEEs.

The distributions of the HB networks (Fig. 9) show that while both the LiDFOB-based and the LiBF₄-based DEEs form similarly strong HBs, the latter shows slightly stronger HBs at shorter distances, corresponding to F–H distances of 1.5–2.0 Å (Fig. 9). Overall, the strongest HBs are found at shorter distances and larger angles (160–180°). In contrast, the BOB anion, being the largest anion with also directional preferences, sterically hinders HB formation and creates a localized HB network, as seen in the HB distribution, which shows a concentration at 6 Å (Fig. 9) and the HB strength, which shows the least intense heat map for its O–H HBs (Fig. 9). In line with these resuults, the HB lifetime analysis (Table S6) shows longer persistence for LiBOB and LiDFOB (O) compared to LiBF₄ and LiDFOB (F), supporting the same trend.

Fig. 9 Top: The HB strength heat map as a function of HBA–HBD distances r in pm and angles in degrees, with warmer colors indicating more probable (stronger) HB geometries. Bottom: The HB distributions of the three DEEs.

Dynamics

Turning swiftly to the dynamics of the DEEs, the total MSD curves (Fig. 10(a)) show an initial non-linear behavior, indicative of a ballistic regime followed by a diffusive motion.

Fig. 10 (a) Total MSDs and (b) self-diffusion coefficients for the ions and the solvent.

The LiBF₄-based DEE exhibits the steepest MSD slope (Fig. 10(a)), indicating the highest ion mobility, followed by LiDFOB and then LiBOB. This trend is also quantified in the inset of Fig. 11(a), where the mobility difference spans nearly a factor of two between LiBF₄ and LiBOB electrolytes. This order is consistent with the calculated self-diffusion coefficients in Fig. 10(b), supporting that LiBF₄ enables the fastest transport behavior.

Fig. 11 (a) Radar plot comparing key molecular features: diffusion ratio, HB localization, anion size, and number of symmetry elements, and (b) bar chart of MHI for the LiBF₄, LiDFOB, and LiBOB-based DEEs.

Fig. 10(b) further reveals distinct differences in species mobility not only across salts but also within each DEE. In the LiBF₄-based DEE, the anion diffuses faster than the cation, with the solvent diffusing fastest overall. In contrast, the cation and anion in the LiDFOB-based and LiBOB-based DEE exhibit nearly identical self-diffusion coefficients, indicating stronger ion pairing and increasingly correlated motion. This narrowing of the cation–anion mobility gap from LiBF₄ → LiDFOB → LiBOB DEEs reflects enhanced ionic association in the latter electrolytes. Similar lithium diffusivity has been reported for LiBF₄ and LiTFSI in ADN-based electrolytes, closely matching the LiBF₄-based DEE.⁵⁰ Meanwhile, ethaline- and glycerol-based DEEs show slower lithium diffusion, comparable to or lower than those of LiDFOB and LiBOB in our study.⁵¹ In contrast, LiPF₆-based electrolytes generally exhibit higher self-diffusion coefficients.^52,53

Cross-comparing species within each electrolyte: in the LiBF₄-based DEE, the trend observed is solvent > anion > cation; while in the LiDFOB-based and LiBOB-based DEEs, solvent > cation ≈ anion. The solvent mobility drops sharply between the LiBF₄-based and the LiDFOB-based DEEs, after which it remains similarly suppressed in the LiBOB-based DEE. This shows that solvent motion becomes increasingly hindered as HB network localization and cation–anion ordering increase, suggesting a dynamically arrested environment in the most structured DEE.

The inverse correlation between localized HB networks (Fig. 8) and ion mobility (Fig. 10(b)), emphasizes the role of intermolecular interactions in determining not only extended structure but also transport properties.^54,55

Molecular-level heterogeneity

Trying to get a comprehensive view of the MLH, we first gather some molecular features: HB localization factor, diffusion ratio, anion size, and number of symmetry elements (Fig. 11(a)), which do show quite some variability across the DEEs but do not allow for any clear discrimination between them.

Quantifying the structural and dynamic disorder within the DEEs by calculating an MHI (Fig. 11(b)), shows how the LiBOB- based DEE exhibits the largest MLH, due to its large anion size, strong cation–anion interactions, and localized HBs.

In contrast, the LiBF₄-based DEE, with its higher ion solvation and more homogeneous HB distribution, shows the lowest MLH. Finally, the LiDFOB-based DEE, with a more moderate MHI suggests both spatial and dynamic heterogeneity, though to a lesser degree than the LiBOB-based DEE.

The physical meaning of the MHIs extend beyond relative comparison. Low MHIs (e.g., for LiBF₄) reflect a more homogeneous molecular environment, associated with delocalized HB, relatively uniform diffusion behavior, and higher symmetry. In contrast, high MHIs (e.g., for LiBOB) indicate heterogeneous structuring and mobility, often arising from large, asymmetric anions and localized HB. Intermediate MHIs (e.g., for LiDFOB) capture partial heterogeneity in either structure or mobility, such as symmetry-breaking without full diffusion suppression. Thus, MHI serves as a semi-quantitative scale for molecular-level disorder, offering insights into dynamic uniformity or solvation–shell asymmetry, which is not apparent from single descriptors.

While the current MHI captures key aspects of MLH, it is worth noting that disorder can manifest itself in various forms. The index relies on five selected descriptors for interpretability, though future refinements could incorporate dynamic properties like HB lifetimes or local viscosity. PCA was limited to two components due to the small number of systems; expanding this would improve variance resolution. Lastly, although normalization by maxima reduces scaling bias, alternative standardization approaches may yield different weight distributions and deserve further investigation.

Conclusions

This research highlights how systematic variation in anion size and symmetry modulates solvation behavior, molecular interactions, transport dynamics, and molecular ordering in DEE-based electrolytes. The small, symmetric LiBF₄-based DEE promotes isotropic coordination and uniform hydrogen bonding, enabling greater ion dissociation and solvent mobility. The LiBOB-based DEE, on the other hand, demonstrates the strongest cation–anion interactions, a localized HB network, and the highest MLH, with its ionic domains more segregated, possibly impacting both ion conduction pathways and perhaps even the electrochemical performance. The LiDFOB-based DEE falls, as expected, in between these two extremes, displaying moderate solvation, HB localization, and ionic mobility. Ultimately, these findings underscore the influence of anion structure on electrolyte properties, hydrogen-bond network arrangements, and MLH, providing valuable insights into the molecular-level interplay affecting electrolyte functionality. Such understanding can inform the rational design of future electrolytes tailored for improved ionic conduction and optimal electrochemical performance.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting the findings of this study are provided in the SI. Supplementary information: This includes representative input files for both xTB (for production runs) and LAMMPS (for molecular dynamics equilibration and production runs). These files outline the key parameters and settings used to generate the simulation data discussed in the manuscript. All additional parameters, such as box sizes, system composition, etc, are described in the main text. Simulations were performed using the open-source software packages xTB (https://xtb-docs.readthedocs.io/) and LAMMPS (https://lammps.org/), and the generated trajectory data were used to compute structural and dynamic properties. The methods used for data analysis and post-processing are detailed in the main text, along with relevant equations and literature references, facilitating the simulation and analysis setup to be reproducible. Further instructions for running the input files and interpreting outputs. See DOI: https://doi.org/10.1039/d5cp01828e

Additional data or clarifications are available from the authors upon reasonable request.

Acknowledgements

The authors would like to acknowledge the financial support from the Swedish Research Council (grants #2020 – 03988 and #2021 – 00613) and VINNOVA/Batteries Sweden (BASE) (grant 2019-00064).

Notes and references

1. J. Cabana, L. Monconduit, D. Larcher and M. R. Palacín, Beyond Intercalation-Based Li-Ion Batteries: The State of the Art and Challenges of Electrode Materials Reacting Through Conversion Reactions, Adv. Mater., 2010, 22, E170–E192.
2. A. Navarro- Suárez and P. Johansson, Perspective—Semi-Solid Electrolytes Based on Deep Eutectic Solvents: Opportunities and Future Directions, J. Electrochem. Soc., 2020, 167, 070511.
3. A. Manthiram, X. Yu and S. Wang, Lithium battery chemistries enabled by solid-state electrolytes, Nat. Rev. Mater., 2017, 2, 1–16.
4. M. E. Di Pietro and A. Mele, Deep eutectics and analogues as electrolytes in batteries, J. Mol. Liq., 2021, 338, 116597.
5. A. Boisset, S. Menne, J. Jacquemin, A. Balducci and M. Anouti, Deep eutectic solvents based on N-methylacetamide and a lithium salt as suitable electrolytes for lithium-ion batteries, Phys. Chem. Chem. Phys., 2013, 15, 20054–20063.
6. H. Liang, H. Li, Z. Wang, F. Wu, L. Chen and X. Huang, New binary room-temperature molten salt electrolyte based on urea and LiTFSI. The, J. Phys. Chem. B, 2001, 105, 9966–9969.
7. O. Geiculescu, D. DesMarteau, S. Creager, O. Haik, D. Hirshberg, Y. Shilina, E. Zinigrad, M. Levi, D. Aurbach and I. Halalay, Novel binary deep eutectic electrolytes for rechargeable Li-ion batteries based on mixtures of alkyl sulfonamides and lithium perfluoroalkylsulfonimide salts, J. Power Sources, 2016, 307, 519–525.
8. Z. Hu, F. Xian, Z. Guo, C. Lu, X. Du, X. Cheng, S. Zhang, S. Dong, G. Cui and L. Chen, Nonflammable nitrile deep eutectic electrolyte enables high-voltage lithium metal batteries, Chem. Mater., 2020, 32, 3405–3413.
9. L. J. B. M. Kollau, M. Vis, A. van den Bruinhorst, A. C. C. Esteves and R. Tuinier, Quantification of the liquid window of deep eutectic solvents, Chem. Commun., 2018, 54, 13351–13354.
10. M. A. R. Martins, S. P. Pinho and J. A. P. Coutinho, Insights into the Nature of Eutectic and Deep Eutectic Mixtures, J. Solution Chem., 2019, 48, 962–982.
11. E. L. Smith, A. P. Abbott and K. S. Ryder, Deep Eutectic Solvents (DESs) and Their Applications, Chem. Rev., 2014, 114, 11060–11082.
12. J. Song, Y. Si, W. Guo, D. Wang and Y. Fu, Organosulfide-Based Deep Eutectic Electrolyte for Lithium Batteries, Angew. Chem., Int. Ed., 2021, 133, 9969–9973.
13. Y. Hu, H. Li, X. Huang and L. Chen, Novel room temperature molten salt electrolyte based on LiTFSI and acetamide for lithium batteries, Electrochem. Commun., 2004, 6, 28–32.
14. S.-D. Han, J. L. Allen, E. Jónsson, P. Johansson, D. W. McOwen, P. D. Boyle and W. A. Henderson, Solvate structures and computational/spectroscopic characterization of lithium difluoro (oxalato) borate (LiDFOB) electrolytes, J. Phys. Chem. C, 2013, 117, 5521–5531.
15. C. W. Yang, Y. H. Ren, B. R. Wu and F. Wu, Formulation of a New Type of Electrolytes for LiNi1/3Co1/3Mn1/3O2 Cathodes Working in an Ultra-Low Temperature Range, Adv. Mater. Res., 2012, 455–456, 258–264.
16. V. Aravindan, J. Gnanaraj, S. Madhavi and H. Liu, Lithium-Ion Conducting Electrolyte Salts for Lithium Batteries, Chem. – Eur. J., 2011, 17, 14326–14346.
17. J. Huang, L.-Z. Fan, B. Yu, T. Xing and W. Qiu, Density functional theory studies on the B-containing lithium salts, Ionics, 2010, 16, 509–513.
18. V. Aravindan and P. Vickraman, Effect of aging on the ionic conductivity of polyvinylidenefluoride–hexafluoropropylene (PVdF–HFP) membrane impregnated with different lithium salts, Indian J. Phys., 2012, 86, 341–344.
19. B. Bang and T. Yeu, Characterization of lithium tetrafluoroborate and N-methylacetamide complex as electric and ionic conductors, Korean J. Chem. Eng., 2016, 33, 1441–1446.
20. S. Kaur, A. Gupta and H. K. Kashyap, Nanoscale spatial heterogeneity in deep eutectic solvents, J. Phys. Chem. B, 2016, 120, 6712–6720.
21. Q. Li and H. Ardebili, Flexible thin-film battery based on solid-like ionic liquid-polymer electrolyte, J. Power Sources, 2016, 303, 17–21.
22. O. S. Hammond, D. T. Bowron and K. J. Edler, Liquid structure of the choline chlorideurea deep eutectic solvent (reline) from neutron diffraction and atomistic modelling, Green Chem., 2016, 18, 2736–2744.
23. S. J. Bryant, A. J. Christofferson, T. L. Greaves, C. F. McConville, G. Bryant and A. Elbourne, Bulk and interfacial nanostructure and properties in deep eutectic solvents: Current perspectives and future directions, J. Colloid Interface Sci., 2022, 608, 2430–2454.
24. O. S. Hammond, D. T. Bowron and K. J. Edler, The effect of water upon deep eutectic solvent nanostructure: an unusual transition from ionic mixture to aqueous solution, Angew. Chem., Int. Ed., 2017, 129, 9914–9917.
25. V. Alizadeh, D. Geller, F. Malberg, P. B. Sánchez, A. Padua and B. Kirchner, Strong microheterogeneity in novel deep eutectic solvents, ChemPhysChem, 2019, 20, 1786–1792.
26. D. O. Abranches and J. A. Coutinho, Everything you wanted to know about deep eutectic solvents but were afraid to be told, Annu. Rev. Chem. Biomol. Eng., 2023, 14, 141–163.
27. L. Percevault, A. Jani, T. Sohier, L. Noirez, L. Paquin, F. Gauffre and D. Morineau, Do deep eutectic solvents form uniform mixtures beyond molecular microheterogeneities?, J. Phys. Chem. B, 2020, 124, 9126–9135.
28. Y. Cui and D. G. Kuroda, Evidence of molecular heterogeneities in amide-based deep eutectic solvents, J. Phys. Chem. A, 2018, 122, 1185–1193.
29. S. Kaur, A. Malik and H. K. Kashyap, Anatomy of microscopic structure of ethaline deep eutectic solvent decoded through molecular dynamics simulations, J. Phys. Chem. B, 2019, 123, 8291–8299.
30. C. Bannwarth, S. Ehlert and S. Grimme, GFN2-xTB—An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions, J. Chem. Theory Comput., 2019, 15, 1652–1671.
31. S. Grimme, C. Bannwarth and P. Shushkov, A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements (Z = 1–86), J. Chem. Theory Comput., 2017, 13, 1989–2009.
32. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, cp2k atomistic simulations of condensed matter systems, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 15–25.
33. R. Andersson, F. Årén, A. A. Franco and P. Johansson, CHAMPION: Chalmers hierarchical atomic, molecular, polymeric and ionic analysis toolkit, J. Comput. Chem., 2021, 42, 1632–1642.
34. S. Spicher and S. Grimme, Robust Atomistic Modeling of Materials, Organometallic, and Biochemical Systems, Angew. Chem., Int. Ed., 2020, 59, 15665–15673.
35. D. J. Evans and B. L. Holian, The Nose–Hoover thermostat, J. Chem. Phys., 1985, 83, 4069–4074.
36. W. Humphrey, A. Dalke and K. Schulten, VMD: Visual molecular dynamics, J. Mol. Graphics, 1996, 14, 33–38.
37. B. G. Levine, J. E. Stone and A. Kohlmeyer, Fast analysis of molecular dynamics trajectories with graphics processing units—Radial distribution function histogramming, J. Comput. Phys., 2011, 230, 3556–3569.
38. M. Brehm and B. Kirchner, TRAVIS-a free analyzer and visualizer for Monte Carlo and molecular dynamics trajectories, J. Chem. Inf. Model., 2011, 51(8), 2007–2023.
39. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1–19.
40. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids, J. Am. Chem. Soc., 1996, 118, 11225–11236.
41. E. S. C. Ferreira, I. V. Voroshylova, C. M. Pereira, M. Natália and D. S. Cordeiro, Improved Force Field Model for the Deep Eutectic Solvent Ethaline: Reliable Physicochemical Properties, J. Phys. Chem. B, 2016, 120(38), 10124–10137.
42. Z. Yu, P. E. Rudnicki, Z. Zhang, Z. Huang, H. Celik, S. T. Oyakhire, Y. Chen, X. Kong, S. C. Kim, X. Xiao, H. Wang, Y. Zheng, G. A. Kamat, M. S. Kim, S. F. Bent, J. Qin, Y. Cui and Z. Bao, Rational solvent molecule tuning for high-performance lithium metal battery electrolytes, Nat. Energy, 2022, 7, 94–106.
43. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, PACKMOL: A package for building initial configurations for molecular dynamics simulations, J. Comput. Chem., 2009, 30, 2157–2164.
44. A. A. H. Pádua, Resolving dispersion and induction components for polarisable molecular simulations of ionic liquids, J. Chem. Phys., 2017, 146(20), 204501.
45. S. Plimpton LAMMPS Documentation. https://www.lammps.org, 2024.
46. D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Elsevier, 2nd edn, 2002.
47. F. Cowell, Measuring Inequality, Oxford University Press, 2011.
48. M. Ue, A. Murakami and S. Nakamura, A Convenient Method to Estimate Ion Size for Electrolyte Materials Design, J. Electrochem. Soc., 2002, 149(10), A1385.
49. D. M. Seo, P. D. Boyle, J. L. Allen, S.-D. Han, E. Jónsson, P. Johansson and W. A. Henderson, Solvate Structures and Computational/Spectroscopic Characterization of LiBF₄ Electrolytes, J. Phys. Chem. C, 2014, 118(32), 18377–18386.
50. O. Borodin and G. D. Smith, Mechanism of Ion Transport in Adiponitrile-Based Electrolytes: Molecular Dynamics Study of LiBF₄ and LiTFSI in ADN, Phys. Chem. Chem. Phys., 2018, 20, 28294–28306.
51. Y. Ding, D. Sun and J. Xu, Lithium Ion Transport Properties in Deep Eutectic Solvents: A Molecular Dynamics Study of Ethaline and Glyceline-Based Electrolytes, Mol. Simul., 2021, 47, 1177–1186.
52. J. Landesfeind and H. A. Gasteiger, Temperature and Concentration Dependence of the Ionic Transport Properties of Lithium-Ion Battery Electrolytes, J. Electrochem. Soc., 2019, 166, A3079–A3097.
53. A. B. Gunnarsdóttir, S. Vema, S. Menkin, L. E. Marbella and C. P. Grey, Investigating the Effect of a Fluoroethylene Carbonate Additive on Lithium Deposition and the Solid Electrolyte Interphase in Lithium Metal Batteries Using In Situ NMR Spectroscopy, J. Mater. Chem. A, 2020, 8, 14975–14992.
54. U. L. Abbas, Q. Qiao, M. T. Nguyen, J. Shi and Q. Shao, Molecular dynamics simulations of heterogeneous hydrogen bond environment in hydrophobic deep eutectic solvents, AIChE J., 2022, 68, e17382.
55. M. T. Ong, O. Verners, E. W. Draeger, A. C. Van Duin, V. Lordi and J. E. Pask, Lithium ion solvation and diffusion in bulk organic electrolytes from first-principles and classical reactive molecular dynamics, J. Phys. Chem. B, 2015, 119, 1535–1545.

---