Monte Carlo simulations of ion correlations and transport in highly concentrated liquid electrolytes

Abstract

Highly concentrated electrolytes (HCEs) exhibit promising properties for battery applications, such as low vapor pressure, high thermal stability and good compatibility with electrode materials. However, the charge and mass properties of HCEs in batteries are strongly influenced by volume conservation constraints as well as by interionic Coulomb interactions and cation–solvent interactions. In order to obtain a better understanding of the influence of volume conservation and different types of interactions on correlated movements of ions and solvent molecules as well as on the resulting transport properties, dynamic Monte Carlo simulations of highly concentrated electrolytes were carried out. Diffusion coefficients and Onsager transport coefficients were calculated for variable molar ratios of salt to solvent, variable partial volume ratios of solvent molecules to anions and for variable interaction strengths. The simulation results were compared to recent experimental results.

---

Introduction

In commercial lithium-ion batteries (LIB), an electrolyte consisting of 1 M LiPF₆ in a mixture of organic carbonates is employed. This electrolyte exhibits a high ionic conductivity, but its thermal and electrochemical stability is limited.¹ In addition, it is flammable leading to battery safety risks.¹ The formation of stable electrode/electrolyte interphases (SEI on the graphite particles of the negative electrode and CEI on the transition metal oxide particles of the positive electrode) allows the LIB to operate beyond the electrochemical stability window of the electrolyte.² In contrast, stable interphases do not generally form in contact to metallic lithium as the negative electrode^2,3 and in contact to high-voltage active materials in the positive electrode, such as LiNi_0.5Mn_1.5O₄.^4,5

Alternative electrolytes are highly concentrated electrolytes (HCEs) containing salt concentrations beyond 3 M.⁶ In these electrolytes, virtually all solvent molecules are bound within the solvation shell of the lithium ions, and virtually no free solvent molecules exist.^7,8 Due to the lack of Li⁺-solvating molecules, also anions take part in the Li⁺ solvation shell.^6,9 This structural difference of the solvation shell as compared to more diluted electrolytes alters the properties of HCEs.⁸ The vapor pressure of HCEs is very low, which significantly reduces their flammability.¹⁰ The thermal stability is significantly improved.¹¹ The SEI formed at the negative electrode is dominated by inorganic decomposition products of the anions, which leads to a denser and more robust SEI.^12–15 In the case of metallic lithium as negative electrode, such an SEI lowers the risk of Li dendrite formation.^16,17 Due to these properties, HCEs are promising electrolytes for high-voltage battery applications.^10,17,18

However, not all properties of HCEs are superior to those of more diluted electrolytes. The strong interactions in HCEs lead to a higher viscosities and to slower charge and mass transport.⁸ The charge and mass transport properties of an electrolyte can be characterized by four parameters, namely three Onsager conductivity coefficients σ₊₊, σ₋₋ and σ₊₋ and a thermodynamic factor Φ.^19–26 In the case of HCEs containing a monovalent salt, the Onsager conductivity coefficients σ_ij are linked to the Onsager transport coefficients L_ij = σ_ij/F², with F denoting the Faraday constant.^27,28 The four target parameters σ₊₊, σ₋₋, σ₊₋, and Φ can be determined by combining various experimental methods, such as electrochemical impedance spectroscopy, electrophoretic NMR and concentration cells with transference.^19,20 It has been shown that overdetermination of the four target parameters by measuring five experimental quantities reduces the uncertainties in the Onsager coefficients and the thermodynamic factor significantly.²⁹ When in addition, the self-diffusion coefficient of the cations D₊ and of the anions D₋ are determined by PFG-NMR, the Onsager coefficients σ₊₊ and σ₋₋ can be split up into their self parts, σ^self₊ and σ^self₋, respectively, and their distinct parts, σ^distinct₊₊ and σ^distinct₋₋, respectively.^{20,23–26,29–31} The coefficients σ^distinct₊₊, σ^distinct₋₋ and σ₊₋ yield information about cation–cation, anion–anion and cation–anion correlations, respectively.²³

Recent experimental results suggest that volume conservation plays an important role for the ion correlations in HCEs.³⁰ It has been demonstrated for several HCE systems that the laboratory frame can be considered as equivalent to a volume-fixed frame of ref. 29, 30, 32 and 33. Within the volume-fixed frame, the sum of the volume fluxes of all molecular species is zero (volume conservation), i.e., , with v_i and J_i denoting the partial molar volume and the molar flux of the molecular species i, respectively. Specifically, it was found experimentally that the relation is valid, if the partial molar volume of Li⁺ ions is neglected and the partial molar volume of the anions is identified with the partial molar volume of the salt.^29,30,32,33 Attractive cation–solvent interactions and interionic Coulomb interactions can also exert a significant influence on ion correlations and transport. Recent results suggest that cation–solvent interactions lead to positive cation–solvent correlations, indicative of a vehicular Li⁺–solvent transport mechanism.³⁰ Attractive cation–anion interactions result in cation and anion movements into the same direction and thus reduce the ionic conductivity.²³

In order to gain deeper insights into the influence of volume conservation, cation–solvent interactions and Coulomb interactions on the ion correlations in HCE, we have carried out Monte Carlo simulations of electrolytes consisting of monovalent cations, monovalent anions and solvent molecules carrying out jumps between the sites of a three-dimensional cubic lattice. We monitor the equilibrium displacement vectors of the three species and calculate diffusion coefficients D_i and Onsager coefficients L_ij. The molar ratio of salt to solvent, the partial molar volume ratio of anions to solvent and several interaction parameters are varied, and their influence on the ion correlations and transport is analyzed.

Methods

288 particles (cations, anions, and solvent molecules) carry out movements on a three-dimensional cubic lattice with spacing a = 1 and with L³ = 70³ lattice sites. The particles exhibit a hard-core radius of 5a, with the center of volume of particles being located on a lattice site. During a single Monte Carlo (MC) time step, the center of volume of a particle makes a jump attempt to a neighboring site. The small jump distance compared to the particle size implies quasi-continuous movements of the particles as expected for low-viscosity liquids. All lattice sites within the hard-core radius around the center-of-volume lattice site are considered as occupied, so that the minimum distance between two neighboring particles is 10a. This implies that the fraction of sites occupied by particles is 43.2%. A high particle density, yet lower than in crystals, is also expected for low-viscosity liquids. The particles move under periodic boundary conditions in all three dimensions. The number of univalent cations on the lattice, N₊, is identical to the number of univalent anions, N₋. The number of solvent molecules is given by N₀ = x·N₊, with x denoting the molar ratio of solvent to salt, see Table 1.

Table 1 Compositions of the studied electrolyte systems. x denotes the molar ratio of solvent to salt

	x	N+ = N−	N0
0.1	10	24	240
0.25	4	48	192
0.5	2	72	144
1.0	1	96	96

---

Fig. 1 illustrates the simulation algorithm. 500 simulation runs are carried out for each studied system. First, the system is brought to equilibrium through thermalization. To this end, 10⁶ Monte Carlo (MC) time steps, as illustrated on the left-hand side of Fig. 1, are performed, while monitoring the total energy of the system. The number of MC time steps is chosen such that the total energy of the system becomes constant during thermalization. Subsequently, an equilibrium simulation is carried out. During the equilibrium simulation, the number of MC time steps is chosen as 10⁵, which is sufficient for obtaining the long-time limits of diffusion coefficients and Onsager coefficients.

Fig. 1 Scheme of the simulation algorithm with the procedure for a single MC time step highlighted in brown and the procedure for a single jump attempt highlighted in grey.

One MC time step encompasses N₊ random selections of a cation, N₋ random selections of an anion, and N₀ random selections of a solvent molecule. Consequently, on the average, each particle is selected once per MC time step. In Fig. 1, the order of the different jumps performed during one MC time step is illustrated on the top right and highlighted in brown. After selecting a particle, the direction for a jump of the center of the particle from its original site (os) to a neighboring target site (ts) is randomly selected with identical probabilities for all six directions. If the jump leads to an overlap with another particle (hard-core interaction), the jump is rejected, and the next particle is selected. If there is no overlap, the energy difference between the particle on the os and on the ts, ΔE_i, is calculated, see Fig. 1. To this end, three different contributions to the system energy are taken into account:

(i) Since experimentally, volume conservation is found under the assumption that the partial molar volume of Li⁺ ions is negligible,^30,32 the center of volume of the particles is computed from the partial molar volumes and position vectors of the anions, v₋ and r⃑_i, respectively, and the partial molar volumes and position vectors of the solvent molecules, v₀ and r⃑_j, respectively:

On the average, the center of volume is located at the center of the lattice, r⃑_CV,avg. The mean-squared deviation of the center of volume from this average value:

R_CV² = (r⃑_CV − r⃑_CV,avg)² (2)

is used for calculating the system energy change caused by a center-of-volume shift, when a particle jumps from its os to its ts:

ΔE_CV = U_CV·ΔR_CV² = U_CV(R_CV,ts² − R_CV,os²) (3)

U_CV denotes an energy penalty parameter for center-of-volume shifts. The harmonic potential implies a backdriving force, which is proportional to the modulus of the center-of-volume displacement from the center of the simulation box.

(ii) Interionic Coulomb interactions are taken into account. The change in the Coulomb energy of the system, ΔE_Coul, caused by a jump of an ion with charge number z_b (z_b = 1 for cations and z_b = −1 for anions) from its os to its ts is given by:

with the parameter U_Coul characterizing the strength of the Coulomb interactions and d₀ = 10a denoting the minimum distance between the centers of two particles. d_i,ts and d_j,ts are the distances between the center of the considered ion on its target site to the other ions, while d_i,os and d_j,os are the same distances for the considered ion being located on its original site. The Madelung factors M account for the interactions of the ion making the jump attempt with ions in image lattices caused by the period boundary conditions. The Madelung factors were calculated prior to the simulations for all possible distance vectors d⃑ on the lattice using the Ewald summation method³⁴ and were used as tabulated values during the simulations.

(iii) We assume that the cation–solvent interactions are based on ion–dipole interactions, which decay as d⁻², with d denoting the distance between the center of a cation and the center of a solvent molecule. During the simulations, the interaction of a cation with all solvent molecules and of a solvent molecule with all cations, respectively, are summed up to a cutoff distance of . The system energy change due to these interactions after a jump of a particle from its os to its ts is then given by:

with N_b denoting the number of interacting particles of the other species b within the cutoff distance. d_i,ts and d_i,os are the distances between the selected particles and the particles of the other species at the target site and at the original site, respectively. U_cat–sol is the attractive energy between a cation and a solvent molecule at the minimum distance d₀ = 10a.

Since highly concentrated electrolytes typically contain large anions with highly delocalized negative charge, the anion–solvent interactions are typically much smaller than the interactions between solvent molecules and small alkali ions, such as Li⁺ ions. Consequently, anion–solvent interactions are neglected in our simulations. The total energy change of the system during a jump attempt of a cation, of an anion, and of a solvent molecule, respectively, is thus given by:

ΔE_cation = ΔE_Coul + ΔE_cat–sol (6)

ΔE_anion = ΔE_CV + ΔE_Coul (7)

ΔE_solvent = ΔE_CV + ΔE_cat–sol (8)

The jump probability w_i of the particle i is calculated according to the Metropolis-Hastings algorithm:³⁵

w_i = min{exp(−ΔE_i),1} (9)

and the jump attempt is accepted with this probability.

During the equilibrium simulation of a system, k independent simulation runs are carried out. During the m^th simulation run, the displacement vectors of the molecular species i and j, Δr⃑_m,i and Δr⃑_m,j, respectively, are monitored. During a single simulation run, s MC steps are performed. After completion of the k simulation runs, the diffusion coefficients of the species, D_i, and the Onsager coefficients L_ij are calculated in Monte Carlo units as follows:

with i,j = [+,−,0].

The diffusion coefficients are normalized in a way that they approach unity in the limit of a very low lattice occupation.

The transport coefficients L_ii can be divided into a self part L^self_i and a distinct part L^distinct_ii by means of the following equations:

L^self_i = N_i·D_i (12)

L^distinct_ii = L_ii − L^self_i (13)

For each system, three different partial volume ratios v₀/v₋ are chosen according to the following experimentally values: v₀/v₋ = 0.67 for sulfolane/LiTFSI, v₀/v₋ = 0.90 for DMC/LiFSI and v₀/v₋ = 1.05 for sulfolane/LiFSI.³⁰ Furthermore, we varied the energy parameters U_CV, U_Coul and U_cat–sol, see Table 2.

Table 2 Energy parameters of the studied sets of systems

Set of systems	UCV	UCoul	Ucat–sol
I	0	0	0
II	1000	0	0
III	1000	5	0
IVa	1000	5	−2.5
IVb	1000	5	−7.0

---

Results and discussion

Exemplary results for time-dependent diffusion coefficients D_i(s) and Onsager coefficients L_ij(s) of all systems are shown in Fig. S1–S5 of the SI. The long-time limits of these transport coefficients, D_i and L_ij, were obtained by averaging over the last 3 × 10⁴ MC time steps during the equilibrium simulation.

In Fig. 2, we show results for D₊, D₋, and D₀ obtained for the systems with a volume ratio v₀/v₋ = 0.67 and with different molar ratios of salt to solvent, 1/x. The diffusion coefficients obtained for the other volume ratios are shown in the SI (Fig. S8 and S12). In the systems with exclusively hard-core interactions (black data points in Fig. 2), the diffusion coefficients of all three species are identical, since particle transport is only limited by the availability of vacant space. When volume conservation with U_CV = 1000 is introduced as an additional constraint, a slight decrease of D₋ and D₀ is observed (blue data points in Fig. 2). This effect is mainly caused by a slightly decreasing jump rate of anions and solvent molecules due to the energy penalty for center-of-volume shifts, see Tables S3 and S6. A tiny decrease of D₊ is also noticeable. The jump rate of the cations is not directly affected by the volume conservation constraint (Tables S3 and S6), but the slight slowing down of anion and solvent molecule dynamics leads to slightly stronger single-particle backward correlations of the cations. This is schematically illustrated in Fig. 3. Without volume conservation, the particle dynamics is governed by the availability of vacant space. If a cation jumps into a neighboring vacant space (jump A), there are identical probabilities for a jump of a neighboring particle (a solvent molecule in Fig. 3) into the same direction (jump B) and for a backward jump of the cation (jump C). However, if the jump rate of anions and solvent molecules is reduced due to the volume conservation constraint, jump B becomes less likely and the backward jump C becomes more likely.

Fig. 2 Diffusion coefficients of (a) the cations D₊, (b) the anions D₋ and (c) the solvent D₀ vs. the molar ratio of salt to solvent, 1/x, for the different energy parameters and for a partial volume ratio of v₀/v₋ = 0.67.

Fig. 3 Schematic illustration of correlated particle movements. If only hard-core interactions are relevant, the particle dynamics is determined by the availability of vacant space (shaded in grey). If a cation jumps to the right into the vacant space (jump A), there are identical probabilities for a movement of a neighboring particle (here a solvent molecule in violet) into the same direction making use of the vacant space left behind by the cation (jump B) and for a backward jump of the cation (jump C). However, if the jump rate of anions and solvent molecules is reduced due to the volume conservation constraint, jump B becomes less likely and the backward jump C becomes more likely.

When Coulomb interactions between the ions are introduced in addition (U_Coul = 5, green data points in Fig. 2), the diffusion coefficients of both cations and anions, D₊ and D₋, decrease, since the attractive interactions between cations and anions hinder their transport. As the solvent movements are linked to the movements of the anions by the volume conservation constraint, D₀ decreases also slightly. The mean distance between the ions is between 12.1a (1/x = 1) and 19.3a (1/x = 0.1) and is thus significantly smaller than the Bjerrum length (50a), implying that the ions are subject to significant Coulomb interactions. On the other hand, higher values for U_Coul lead to ion clustering and even phase separation between ion clusters and solvent molecules and are thus not considered here.

Additional attractive cation–solvent interactions lead to a drop of the diffusion coefficients of cations and solvent molecules, D₊ and D₀, which becomes stronger with increasing interactions strength U_cat–sol and with increasing salt concentration. At high interaction strengths and high salt concentrations, virtually all solvent molecules are bound to cations and virtually no “free” solvent molecules exist in the system. The diffusion coefficient of the anions, D₋, is virtually unaffected by weak cation–solvent interactions (U_cat–sol = −2.5, pink data points in Fig. 2), but decreases significantly at stronger cation–solvent interactions (U_cat–sol = −7.0 red data points in Fig. 2). This decrease is caused by the interplay between cations–solvent interactions and volume conservations leading to interrelations between the movements of cations, anions and solvent molecules.

The results obtained for the Onsager coefficients L_ij at all volume ratios v₀/v₋ and molar ratios 1/x are shown in the SI (Fig. S6, S7, S10, S11, S14 and S15).

From the Onsager coefficients, the parameters α, β, and γ were calculated as:

The parameter α gives information about the relative mobility of cations and anions (α > 0.5: cations are more mobile; α < 0.5: anions are more mobile). The parameter β is a measure for cation–anion correlations (β = 1: strictly vehicular cation–anion transport and zero ionic conductivity; β = 0: no cation–anion correlations; β = −1: strong cation–anion anticorrelations). The parameter γ is a measure for cation–solvent correlations³⁰ (γ = 1: strictly vehicular cation–solvent transport; γ = 0: no cation–solvent correlations). In Fig. 4 we show plots of the parameters α, β, and γ versus the molar ratio 1/x at a volume ratio v₀/v₋ = 0.67. The results at the other volume ratios are shown in the SI (Fig. S9 and S13).

Fig. 4 Parameters α, β and γ plotted versus the molar ratio 1/x for the different energy parameters and for a volume ratio of v₀/v₋ = 0.67.

In systems with only hard-core interactions (black data points in Fig. 4), the collective transport of cations and anions is equally fast, manifesting in α = 0.5, independent of 1/x. The limited amount of vacant space leads to positive correlations between the movements of all particles, as sketched in Fig. 5.³⁶ With increasing molar ratio 1/x, more positively correlated cations and anions lead to higher values of L₊₋ and thus to higher values of β. The parameter γ increases also with increasing 1/x due to the increasing number of positively correlated cations and solvent molecules.

Fig. 5 Schematic illustration of correlated movements of cations (green), anions (red) and solvent molecules (violet) for different energy parameters and their impact on the Onsager cross coefficients L₊₋, L₋₀ and L₊₀.

In systems with strict volume conservation, the zero net volume flux together with a negligible partial volume of the cations leads to the following relations between different Onsager coefficients:³⁰

In the SI (Fig. S16) it is shown that an energy penalty parameter for center-of-volume shifts of U_CV = 1000 is sufficiently high for ensuring that the simulation results obey eqn (17)–(19). Overall, volume conservation leads to anion–anion, anion–solvent, and solvent–solvent anticorrelations. As shown in Tables S2 and S5, the Onsager coefficients L^distinct₋₋, L₋₀ and L^distinct₀₀, which are all positive for systems with only hard-core interactions, become negative for systems with strict volume conservation. A negative value of L^distinct₋₋ reduces the Onsager coefficient L₋₋ and increases the parameter α, see Fig. 4(a). The increase of α is strongest for the lowest volume ratio v₀/v₋ = 0.67, i.e. when the anions are significantly larger than the solvent molecules, see also Fig. S17.

As shown in the SI, strict volume conservation implies also that:

This equation can be used to analyze the relative contributions of anion–anion anticorrelations, solvent–solvent anticorrelations, and anion–solvent anticorrelations to volume conservation. Fig. 6 shows the relative weights of the three terms on the right-hand side of the equation for volume ratios of v₀/v₋ = 0.67 and v₀/v₋ = 1.05, respectively. The relative contribution of anion–anion correlations to volume conservation, as described by the first term on the right-hand side, increases strongly with increasing molar ratio of the salt, while the contribution of solvent–solvent anticorrelations, as described by the second term, shows the opposite trend. That implies that at low salt concentrations, solvent–solvent anticorrelations mainly ensure volume conservation, while with increasing salt concentration, anion–anion anticorrelations play an increasingly important role. The relative importance of anion–solvent anticorrelations increases also with increasing salt to solvent ratio. With increasing volume ratio v₀/v₋, the relative contribution of the first term becomes smaller, as becomes evident when comparing Fig. 6(a) and (b). Results for the volume ratio v₀/v₋ = 0.90 are shown in Fig. S18.

Fig. 6 Relative contributions of L^dist₋₋, L^dist₀₀ and L₋₀ to volume conservation, as described by eqn (20), for different molar ratios 1/x in systems with hard-core interactions between the particles and strict volume conservation (U_CV = 1000). Results for two different volume ratios (a) and (b) are shown.

Moreover, volume conservation diminishes the positive correlations between cations and anions and between cation and solvent molecules, which were found in the case of only hard-core interactions, resulting in β and γ values close to zero, see Fig. 4 and 5.

Next we consider systems with hard-core interactions, strict volume conservation (U_CV = 1000) and with Coulomb interactions between the ions characterized by the parameter U_Coul = 5. As mentioned above, the Coulomb interactions lower the diffusion coefficients of both anions and cations and thus the self parts L^self₊ and L^self₋ (eqn (12)). Furthermore, the anticorrelations between the anions are weakened, leading to more positive values of L^distinct₋₋, and the positive correlations between the cations are weakened, leading to more negative values of L^distinct₊₊ (see Tables S5 and S8 in the SI). In total, these effects lead to slightly reduced values of α, see Fig. 4.

At low molar ratios 1/x, the attractive cation–anion interactions result in positive correlations between the movements of cations and anions, manifesting in positive β values around 0.2, see Fig. 4. At higher molar ratios 1/x, the average distance between the ions becomes smaller, and the Coulomb barrier for a cation to escape from neighboring anions and for an anion to escape from neighboring cations, respectively, becomes lower. Consequently, cations and anions move more independently, manifesting in a drop of β. In Fig. S23, we show results for the dependence of the Onsager coefficient L₊₋ on the initial distance between cations and anions. Enhanced short-range positive cation–anion correlations under the influence of Coulomb interactions are observed, which can be interpreted as a dynamic ion pairing. The combined influence of Coulomb attraction between cations and anions and volume conservation leads to negative correlations between cations and solvent molecules, which are illustrated in Fig. 5 and which manifest in negative γ values.

When additional attractive interactions between cations and solvent molecules are taken into account, cations are slowed down by the interactions with both solvent molecules and anions. This effect is strongest at low molar ratios 1/x leading to α values close to 0, see Fig. 4(a). With increasing molar ratio 1/x, the barriers for a cation to escape from neighboring anions and neighboring solvent molecules decreases, and consequently the transport of the cations becomes faster. At high molar ratios 1/x, the α values are significantly larger than 0.5, since the anion transport is hindered by volume conservation, while the Coulomb barriers for the cation transport are low. Regarding the cation–anion correlation parameter β, there is a competition between cation–solvent and cation–anion interactions. As mentioned above, the combined action of volume conservation and attractive cation–anion interactions leads to positive values of β and to negative values of γ. Strong cation–solvent interactions, on the other hand, lead to positive cation–solvent correlations, as sketched in Fig. 5, and thus to positive values of γ, which in combination with volume conservation implies negative values of β. This effect becomes maximal at 1/x = 1.0 with β = −0.27 and γ = +0.36. Values in this range were also found experimentally for high-concentration electrolytes, like LiFSI or LiTFSI in sulfolane and LiFSI in dimethlycarbonate.^29,30 This good agreement between simulations results and experimental results confirms the delicate interplay between volume conservation, Coulomb interactions and cation–solvent interactions with regard to ion correlations and transport in HCEs. Remarkably, hydrodynamic effects, as observed by Diddens and Heuer in MD simulations of liquid electrolytes,³⁷ were also found in our MC simulations when analyzing the dependence of Onsager cross coefficients on the initial distance between the particles, see Fig. S24 in the SI. If the initial distance between an anion and a solvent molecule is short, there are positive anion–solvent correlations due to the limited availability of vacant sites. Macroscopically, this effect is overcompensated by anticorrelated movements of anions and solvent molecules on larger distances (quasi-hydrodynamic flow). Overall, this leads to negative values of L₋₀.

Finally, we consider the charge and salt transport properties of the simulated electrolytes, which can be described by the transport coefficients:

L_ion = L₊₊ + L₋₋ − 2L₊₋ (21)

For enabling fast cation transport in a battery under anion-blocking conditions, both transport coefficients should be high.^21,22 Fig. 7 illustrates the influence of various energy parameters on these two coefficients for the most concentrated electrolytes with 1/x = 1.0. As shown in the SI, the ratio of the transport coefficients is related to the parameters α and β via:

Fig. 7 Influence of different energy parameters on the salt transport coefficient L_salt (orange data points) and on the charge transport coefficient L_ion (blue data points) of concentrated electrolytes with 1/x = 1.0 and . (I) Only hard-core interactions, (II) U_CV = 1000, (III) U_CV = 1000 and U_Coul = 5 and (IV) U_CV = 1000, U_Coul = 5 and U_cat–sol = −7.0.

In the case of only hard-core interactions, salt and charge transport are almost equally fast. As compared to an ideal electrolyte with α = 0.5, β = 0, and L_salt/L_ion = 0.25, the hard-core interactions lead to positive cation–anion correlations, which reduce charge transport and enhance salt transport. Volume conservation reduces the positive cation–anion correlations and slows down neutral salt transport. In addition, volume conservation leads to faster cation transport as compared to anion transport (α ≈ 0.8). This reinforces the slowing down of neutral salt transport, since salt transport is limited by the transport of the slower ion. Additional Coulomb interactions lead to positive correlations between cations and anions (β ≈ 0.06), which reduces charge transport, while neutral salt transport remains virtually unchanged. Additional strong cation–solvent interactions slow down the dynamics of all species in the electrolyte and reduce their transport coefficients. Furthermore, in combination with volume conservation, the interactions lead to significant anticorrelations between cations and anions, which slow down neutral salt transport and reduce the ratio L_salt/L_ion to values around 0.1. These values are in good agreement with experimental values of L_salt/L_ion between 0.06 and 0.13 obtained for sulfolane/LiTFSI (1/x = 0.50), sulfolane/LiFSI (1/x = 0.42) and DMC/LiFSI (1/x = 0.91) electrolytes.^29,30 This good agreement suggests, on the one hand, that in these HCEs, cation–solvent interaction strengths are at least comparable or higher than Coulomb interactions strengths. In this sense, the solvents in these HCEs can be considered as strongly solvating. On the other hand, our simulation results show that very strong attractive cation–solvent interactions in a highly concentrated electrolyte are detrimental for the application in a battery, since neutral salt transport is slowed down strongly.

Conclusions

Dynamic Monte Carlo simulations of binary highly concentrated liquid electrolytes were carried out at four molar ratios of salt to solvent, 1/x = 0.10, 0.25, 0.5 and 1.0, and three partial volume ratios of solvent to anions, v₀/v₋ = 0.67, 0.90 and 1.05. The influence of hard-core interactions, volume conservation, interionic Coulomb interactions and attractive cation–solvent interactions on dynamics and transport of ions and solvent molecules was analysed. In systems with exclusively hard-core interactions, positive correlations exist between the movements of all particles due to the limited amount of vacant space. Volume conservation reduces these positive correlations significantly and causes anticorrelated movements of anions and solvent molecules. Additional Coulomb interactions result in positive cation–anion correlations reducing charge transport. Additional attractive cation–solvent interactions compete with cation–anion interactions, slow down the overall dynamics and promote a vehicular-type cation–solvent transport. In combination with volume conservation, the vehicular transport goes along with anticorrelated cation–anion movements. Consequently, strong cation–solvent interactions hinder neutral salt transport strongly and are thus detrimental for cation transport under anion-blocking conditions in a battery.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: Fig. S1–S5: time-dependent diffusion coefficients and Onsager coefficients; Tables S1–S15: self and distinct parts of Onsager coefficients and jump rates of particles; Fig. S6 and S7: composition-dependent Onsager coefficients for systems with ; Fig. S8–S15: composition-dependent diffusion coefficients and correlation parameters for systems with and with ; Fig. S16: simulation results for different energy penalty parameter for center-of-volume shifts U_CV; Fig. S17: influence of the volume ratio v₀/v₋ on the parameter α; derivation of eqn (20) for analyzing the relative contributions of anions and solvent molecules to volume conservation; Fig. S18: relative contributions to the volume conservation for systems with ; derivation of eqn (23) for the ratio L_salt/L_ion. Fig. S19-S22: Composition-dependent Onsager coefficients for systems with and two different cutoffs for the attractive cation-solvent interaction. Fig. S 23: Dependence of the Onsager coefficient L₊₋ on the initial cation-anion distance d for two systems with a molar ratio of 1/x = 1, and with strict volume conservation (U_CV = 1000) and Coulomb interactions (U_Coul = 5). Fig. S24: Dependence of the Onsager coefficients L^distinct₊₊, L^distinct₋₋, L^distinct₀₀, L₊₋, L₊₀, and L₋₀, on the initial distance between the respective particles for systems with , 1/x = 1 and 1/x = 0.25, U_CV = 1000, U_Coul = 5, and U_cat–sol = −7. See DOI: https://doi.org/10.1039/d6cp00330c.

Acknowledgements

We would like to thank the German Science Foundation (DFG) for financial support of this work under the project number Ro1213/19-1. Furthermore, we are grateful to the Competence Center for High Performance Computing in the State of Hesse (MarC3a Cluster, University of Marburg) for providing computational resources.

References

1. A. Mauger and C. M. Julien, Ionics, 2017, 23(8), 1933.
2. K. Xu, Chem. Rev., 2004, 104(10), 4303.
3. W. Xu, J. Wang, F. Ding, X. Chen, E. Nasybulin, Y. Zhang and J.-G. Zhang, Energy Environ. Sci., 2014, 7(2), 513.
4. A. Manthiram, ACS Cent. Sci., 2017, 3(10), 1063.
5. N. Nitta, F. Wu, J. T. Lee and G. Yushin, Mater. Today, 2015, 18(5), 252.
6. O. Borodin, J. Self, K. A. Persson, C. Wang and K. Xu, Joule, 2020, 4(1), 69.
7. Y. Yamada and A. Yamada, J. Electrochem. Soc., 2015, 162(14), A2406–A2423.
8. D. W. McOwen, D. M. Seo, O. Borodin, J. Vatamanu, P. D. Boyle and W. A. Henderson, Energy Environ. Sci., 2014, 7(1), 416.
9. X. Cao, H. Jia, W. Xu and J.-G. Zhang, J. Electrochem. Soc., 2021, 168(1), 10522.
10. J. Wang, Y. Yamada, K. Sodeyama, C. H. Chiang, Y. Tateyama and A. Yamada, Nat. Commun., 2016, 7, 12032.
11. S. Kim, B. Seo, H. V. Ramasamy, Z. Shang, H. Wang, B. M. Savoie and V. G. Pol, ACS Appl. Mater. Interfaces, 2022, 14(37), 41934.
12. T. T. Hagos, B. Thirumalraj, C.-J. Huang, L. H. Abrha, T. M. Hagos, G. B. Berhe, H. K. Bezabh, J. Cherng, S.-F. Chiu, W.-N. Su and B.-J. Hwang, ACS Appl. Mater. Interfaces, 2019, 11(10), 9955.
13. G. Jiang, F. Li, H. Wang, M. Wu, S. Qi, X. Liu, S. Yang and J. Ma, Small Struct., 2021, 2(5), 2000122.
14. Y. Jie, X. Ren, R. Cao, W. Cai and S. Jiao, Adv. Funct. Mater., 2020, 30(25), 1910777.
15. T. Li, X.-Q. Zhang, N. Yao, Y.-X. Yao, L.-P. Hou, X. Chen, M.-Y. Zhou, J.-Q. Huang and Q. Zhang, Angew. Chem., Int. Ed., 2021, 60(42), 22683.
16. C. Tian, K. Qin and L. Suo, Mater. Futures, 2023, 2(1), 12101.
17. Y. Li, M. Liu, K. Wang, C. Li, Y. Lu, A. Choudhary, T. Ottley, D. Bedrov, L. Xing and W. Li, Adv. Energy Mater., 2023, 13(30), 2300918.
18. X. Fan, L. Chen, X. Ji, T. Deng, S. Hou, J. Chen, J. Zheng, F. Wang, J. Jiang, K. Xu and C. Wang, Chem, 2018, 4(1), 174.
19. D. Dong, F. Sälzer, B. Roling and D. Bedrov, Phys. Chem. Chem. Phys., 2018, 20(46), 29174.
20. S. Pfeifer, F. Ackermann, F. Sälzer, M. Schönhoff and B. Roling, Phys. Chem. Chem. Phys., 2021, 23(1), 628.
21. B. Roling, J. Kettner and V. Miß, Energy Environ. Mater., 2022, 5(1), 6.
22. B. Roling, V. Miß and J. Kettner, Energy Environ. Mater., 2024, 7(1), e12533.
23. N. M. Vargas-Barbosa and B. Roling, ChemElectroChem, 2020, 7(2), 367.
24. S. Ikeda, S. Tsuzuki, T. Sudoh, K. Shigenobu, K. Ueno, K. Dokko, M. Watanabe and W. Shinoda, J. Phys. Chem. C, 2023, 127(28), 13837.
25. K. Shigenobu, K. Dokko, M. Watanabe and K. Ueno, Phys. Chem. Chem. Phys., 2020, 22(27), 15214.
26. K. Shigenobu, S. Tsuzuki, F. Philippi, T. Sudoh, Y. Ugata, K. Dokko, M. Watanabe, K. Ueno and W. Shinoda, J. Phys. Chem. B, 2023, 127(48), 10422.
27. K. D. Fong, H. K. Bergstrom, B. D. McCloskey and K. K. Mandadapu, AIChE J., 2020, 66(12), e17091.
28. S. Kjelstrup, A. F. Gunnarshaug, Ø. Gullbrekken, S. K. Schnell and A. Lervik, J. Chem. Phys., 2023, 159(3), 34104.
29. T. Pothmann, M. Middendorf, C. Gerken, P. Nürnberg, M. Schönhoff and B. Roling, Faraday Discuss., 2024, 253, 100.
30. H. Kilian, T. Pothmann, M. Lorenz, M. Middendorf, S. Seus, M. Schönhoff and B. Roling, Phys. Chem. Chem. Phys., 2025, 27(3), 1593.
31. Z. Li, R. Bouchal, T. Mendez-Morales, A.-L. Rollet, C. Rizzi, S. Le Vot, F. Favier, B. Rotenberg, O. Borodin, O. Fontaine and M. Salanne, J. Phys. Chem. B, 2019, 123(49), 10514.
32. M. Lorenz, F. Kilchert, P. Nürnberg, M. Schammer, A. Latz, B. Horstmann and M. Schönhoff, J. Phys. Chem. Lett., 2022, 13(37), 8761.
33. F. Kilchert, M. Lorenz, M. Schammer, P. Nürnberg, M. Schönhoff, A. Latz and B. Horstmann, Phys. Chem. Chem. Phys., 2023, 25(38), 25965.
34. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford Science Publications, Clarendon Pr, Oxford, 1987.
35. W. K. Hastings, Biometrika, 1970, 57(1), 97.
36. B. J. Morgan, R. Soc. Open Sci., 2017, 4(11), 170824.
37. D. Diddens and A. Heuer, J. Chem. Phys., 2023, 158(15), 154112.

---