Mesoscopic Inhomogeneities in Concentrated Electrolytes.

A mesoscopic theory for water-in-salt electrolytes combining density functional and field-theoretic methods is developed in order to explain the unexpectedly large period of the oscillatory decay of the disjoining pressure observed in recent experiments for the lithium bis(trifluoromethylsulfonyl)-imide (LiTFSI) salt [T. S. Groves et al., J. Phys. Chem. Lett. 2021, 12, 1702]. We assumed spherical ions with different diameters and implicit solvent, inducing strong, short-range attraction between ions of the same sign. For this highly simplified model, we calculated correlation functions. Our results indicate that mesoscopic inhomogeneities can occur when the sum of the Coulomb and the water-mediated interactions between like ions is attractive at short and repulsive at large distances. We adjusted the attractive part of the potential to the water-in-LiTFSI electrolyte and obtained both the period and the decay rate of the correlations, in semiquantitative agreement with the experiment. In particular, the decay length of the correlations increases nearly linearly with the volume fraction of ions.

Introduction

For many years it was commonly assumed that dilute electrolytes, very well described by the Debye–Hückel (DH) theory, are more suitable for electrochemical devices than the concentrated ones. For this reason, neither experimentalists nor theoreticians paid much attention to the concentrated electrolytes. Recently, however, it was noted that the concentrated electrolytes have advantages such as a large electrochemical stability window, and they may find applications in electrochemical devices, for example, in lithium-ion batteries.[1−4] These observations motivated intensive experimental studies. On one hand, lithium salt/ionic liquid mixtures have been studied intensely because their unique properties of ionic liquids (IL) such as low volatility, low flammability, and high chemical and thermal stability make them attractive for technological applications as electrolytes. However, recent electrophoretic NMR measurements detected a negative alkali cation contribution to the total conductivity in a number of this type of mixture.[5,6] The experimental results were explained by the theory[7] and computer simulations.[8,9] On the other hand, water-in-salt electrolytes draw increasing attention as promising candidates for replacing current lithium-ion technology in state-of-the-art lithium-ion batteries.[1,2] It was found that when the concentration of ions, ρ, increases the deviation between experimental results and predictions of the DH theory becomes very large. Even qualitative trends well documented for the dilute electrolytes, such as decreasing screening length with increasing ρ, are opposite in concentrated electrolytes and IL solutions. Recent surface-force balance (SFB) experiments show that in the above systems the screening length λs is proportional to ρ, while the Debye screening length λD perfectly describing the dilute electrolytes is proportional to . The scaling behavior λs/λD ∼ (a/λD)3, where a is the average diameter of the ions, was found for a number of concentrated solutions of simple salts in water, IL solutions, and alkali halide solutions.[10−12] This universal behavior suggests that the observed decay of the disjoining pressure depends not on specific interactions but only on the Coulomb potential that is common for all the studied systems.

The strong disagreement of experimental results for concentrated electrolytes with DH theory predictions attracted the attention of theoreticians, but despite significant effort in theoretical and simulation studies, the experimental results are not fully explained yet.[13−22] In several theories and simulation studies, the scaling behavior λs/λD ∼ (a/λD)α was found, but the scaling exponent as well as λs in these studies were significantly smaller than in the experiments.[15−19,21,22] Correct scaling for the charge–charge correlation length (that should be equal to λs) was obtained in ref (23), where it was shown that the variance of the local charge density plays a significant role for large concentrations of ions. However, oscillatory decay obtained in theory is at variance with the asymptotic monotonic decay of the disjoining pressure observed in the experiments.[4,10−12]

Theoretical studies of ionic systems are very often based on the restricted primitive model (RPM) .[24,25] In the RPM, the ions are treated as charged hard spheres with the same diameter, and the Coulomb potential between the ions is assumed. The solvent, however, is treated as a dielectric continuum. Theoretical results for the RPM[26−28] show that when the density of ions is very low and the temperature is very high the charge–charge correlations exhibit exponential decay with the decay length λs = λD, signaling that the ions behave as point charges and the DH theory is valid. For increasing ρ and/or decreasing T, however, the sizes of the ions become more and more important, and λs deviates more and more from λD. Moreover, a second exponential term with a decay length smaller than λs and an amplitude with opposite sign becomes relevant. The two decay lengths merge at the Kirkwood line[29] on the (ρ, T) diagram. At the large-density, low-temperature side of the Kirkwood line, the charge–charge correlation function exhibits an oscillatory decay with the period ∼2a.[26,27] The neighborhood of opposite charges is more probable than the neighborhood of like charges, and at the Kirkwood line the density is large enough and the temperature is low enough to allow for formation of clusters with oppositely charged nearest neighbors, where the repeat unit is 2a. The predictions concerning the period of the damped oscillations of the correlation function were verified by simulations and experiments.[11,12,22,30−32] The decay length of the correlations in the RPM, however, depends on the approximation used in theoretical studies and remains a question of a debate.[23,26,28]

For the charge density profile near a planar electrode, as well as for the disjoining pressure between parallel planar electrodes, the decay length and the period of oscillations are expected to be the same as the corresponding length in the charge–charge correlation function in the bulk electrolyte at the same thermodynamic state. In the SFB experiments, the oscillatory decay of the disjoining pressure was observed up to some distance between the crossed mica cylinders, but the asymptotic decay at larger distances was monotonic.[10−12]

Strong disagreement with the RPM predictions was observed in recent SFB experiments for concentrated lithium bis(trifluoromethylsulfonyl)-imide (LiTFSI) in water.[4] It was found that the period of oscillations of the disjoining pressure was twice as large as the sum of diameters of the cation and the anion, observed previously in many concentrated electrolytes and predicted by the RPM. The same length scale of inhomogeneities in the bulk was observed in scattering experiments for concentrated LiTFSI.[3] At large distances, the decay of the disjoining pressure changes from oscillatory to monotonic, as found before for the other systems. Importantly, the decay length increases with increasing salt concentration and is of the same order of magnitude as observed previously in various concentrated electrolytes.[4]

In the LiTFSI salt, the size and the chemical properties of the TFSI– and Li+ ions are significantly different.[3] The TFSI– ion is not spherical and is much larger than the Li+ ion. Moreover, the Li+ ions are very well solvated in water, in contrast to the hydrophobic TFSI– ions. Based on the above observations, we conclude that in the above water-in-salt electrolyte the size difference between the ions and/or the specific non-Coulombic interactions must play a significant role and cannot be neglected.

The effect of specific interactions was studied in the RPM supplemented with additional short-range (SR) interactions in ref (33). On the other hand, the size difference between the ions with neglected SR (primitive model (PM)) was studied in refs (34−37). It was found that the length scale of inhomogeneities depended on the strength of the SR interactions and on the size asymmetry. These general predictions were not verified by experiments, however. In the particular case of LiTFSI, the question of the origin of the scale of inhomogeneities and of the range of correlations is open.

In this work, we develop a minimal model for the water-in-salt electrolyte and fit the parameters to the particular case of the LiTFSI salt. The minimal model can allow us to see which details of the system properties can be neglected without changing the key features of the decay of the correlations. The important questions are: (i) To what extent can properties of the systems with large size asymmetry of ions and with strong specific interactions depend on details of the interactions and on the geometry of the ions? (ii) Can the solvent be treated as a dielectric continuum that mediates effective ion–ion interactions, or must it be taken into account explicitly? (iii) What types of approximations should be used to compute the correlation functions reproducing the qualitative trends found in experiments?

In our model, we assume that both the size difference and the effective SR interactions between the ions must be taken into account. The solvent, however, can be treated as a dielectric continuum. We further assume that water induces effective ion–ion interactions of a range much shorter than the range of the Coulomb potential. The model is introduced in Section . In Section , we develop an approximate form of the grand thermodynamic potential functional of local ionic densities that allows us to obtain correlation functions. In Section , we calculate the correlation functions first in mean-field approximation (MF) (Section ) and next beyond the MF, using the procedure developed in refs (38 and 39) (Section ). We obtain a semiquantitative agreement with experiment when we assume strong water-induced, short-range attraction between the cations. The last section contains the discussion and summary.

Construction of the Theoretical Model

In this section, we describe the assumptions and the approximations leading to the effective specific interactions between the ions in the LiTFSI salt dissolved in water. We take into account the ionic sizes reported in ref (4). In addition, we develop the approximation for the SR interactions based on the requirement that the model predicts inhomogeneities at the length scale ∼2(σ+ + σ–) = 4a, where σ± denotes the diameter of the corresponding ion. Once the form of the SR interactions is assumed, we calculate the correlation functions for different volume fractions of ions without further fitting of any parameters, using the method summarized in Section .

We first consider the excluded-volume interactions. The TFSI– ions are much bigger than the Li+ ions and have a shape of an ellipsoid. We assume that the size difference is more important than the nonspherical shape and assume that the Li+ and TFSI– ions can be modeled as charged hard spheres with the diameter σ+ = 0.2 nm and σ– = 0.6 nm, respectively.

From the previous studies,[3] we know that the Li+ ions are strongly hydrophilic, while the TFSI– ions are strongly hydrophobic. The Li+ ions are solvated by water, and the TFSI– ions are not. This means that the short-range non-Coulombic forces have a strong effect on the distribution of the ions and cannot be neglected. Even though the water–Li+ interactions play an important role, we assume that water can be treated as a dielectric continuum, as in the PM and the classical DH theory. We assume, however, that water molecules mediate effective interactions between the Li+ ions and that these solvent-induced effective interactions are short-ranged (compared to the Coulomb potential) but strong. In addition, we assume short-range interactions between the TFSI– ions but neglect the interactions between the Li+ and TFSI– ions other than the Coulomb potential. We do not attempt to determine the precise shape of these interactions because the correlations at large distances depend only on some gross features of the potentials, such as their range and strength. This is because in collective phenomena involving many ions and solvent molecules many details are washed out by averaging over all distributions of the ions and the solvent molecules.

With the above assumptions, we model the aqueous electrolyte solution as a binary mixture of oppositely charged hard spheres of different diameters (σ+ ≠ σ–) immersed in a structureless dielectric medium with the dielectric constant ε. We limit ourselves to the model with monovalent ions (z+ = z– = 1). The presence of the solvent is taken into account through the solvent-induced effective short-range attractive interactions between the ions of the same sign. Therefore, we assume that the pair interaction potentials between two ions for r > σαβ = (σα + σβ)/2, with α = +, – and β = +, – , can be presented in the formHere, UαβC(r) are the Coulomb potentials between the ions with the signs α and β. As a length unit we choose the sum of radii, a = (σ+ + σ–)/2. The Coulomb potentials for r* ≡ r/a arewhere a± = σ±/aand β = 1/kBT with kB and T denoting the Boltzmann constant and temperature, respectively. lB is the Bjerrum length in a units; EC is the electrostatic potential of the pair of oppositely charged ions at contact; and the reduced temperature T* is in units of EC. The unit step function θ(x) = 1 for x > 0 and θ(x) = 0 for x < 0 prevents contributions to the electrostatic energy of the pair of ions that would come from the forbidden overlap of the hard cores.

For UααA(r), we assume a short-range attractive potential. The form of the sum of the direct (van der Waals type) and solvent-induced interactions is unknown, but we assume that its detailed shape is not necessary for studies of the collective phenomena such as the long-distance correlations. For simplicity of calculations, we assume the attractive Yukawa potentials for UααA(r)where zα* is in the a–1 units, and we assume that z+* = z–* = z = 1.8 to ensure fast decay of these interactions. ϵαα* measures the strength of the effective non-Coulombic interactions in units of EC. Introducing the size-asymmetry parameterwe get a± = 1 ∓ δ.

We assume that for the TFSI– ions ϵ––* = 1 and treat ϵ++* as a fitting parameter. In order to find the best approximation for the Li+ ions in water, we calculate the length scale of inhomogeneities for arbitrary ϵ++* in Section . We find that ϵ++* = 5 leads to a satisfactory agreement of the length scale of inhomogeneities with the experimental results. Based on this observation, we assume ϵ++* = 5 for the considered system. Note that the second important length scale, λs, will be determined for several concentrations of ions without additional fitting.

In Figure (a) the potentials Uαβ(r) normalized by EC = e2/(aε) are shown for the model with ϵ++* = 5, ϵ––* = 1, and δ = 0.5. For the chosen parameters, the interaction potentials between like ions consist of short-range attraction (SA) and electrostatic long-range repulsion (LR). Competing interaction potentials of this kind are also known as SALR potentials. In one-component systems, the SALR-type interactions can lead to spontaneously formed stable aggregates of particles, such as spherical or elongated clusters, networks, or layers.[40−42] The Fourier transforms of the potentials Uαβ(r)/Ec for the above-mentioned model are shown in Figure (b), and the corresponding expressions are given in the Supporting Information. It is seen that Ũ++(k) and Ũ––(k) take minima at k ≠ 0 which are rather close to each other. Ũαα(k) < 0 means that the density wave of the α ions with the wavelength 2π/k leads to the energy lower than in the homogeneous state, signaling a high probability of such waves.

Figure 1

Interaction potentials U(r) (eqs –6) for the model with δ = 0.5, ϵ++* = 5, and ϵ––* = 1 in real space (a) and in Fourier representation (b). U are in units of EC. EC and δ are defined in eqs and 7. r and k are in a and a–1 units, respectively, where a = (σ+ + σ–)/2.

Theoretical Formalism: A Brief Summary

In this section, we present a brief description of the mesoscopic theory for inhomogeneous mixtures.[38,43] In this theory, we divide the system into subsystems having mesoscopic sizes and consider the local volume fraction of the anions and the cations in each subsystem. Thus, our approach with the focus on each mesoscopic part of the system and its state differs from the approach with the focus on the ions and their states. The linear size of the mesoscopic subsystems should not be larger than the scale of the inhomogeneities. The local volume fraction of the ions with the α sign in the mesoscopic subsystem with the center at r is denoted by ζα(r). By analogy with the macroscopic volume fraction, ζα(r) denotes the fraction of the volume of the mesoscopic region that is occupied by the α-sign ions. For a given macroscopic volume fraction, there are plenty of microscopic states with different positions of the ions. Likewise, for the given ζα(r), the positions of the ions can be different. The functions ζα(r) can be considered as constraints imposed on the microscopic states. This constraint imposed on a particular subsystem at r is analogous to the constraint of a fixed number of ions in the whole system in the canonical ensemble. The canonical and grand canonical ensembles are equivalent in the thermodynamic limit because the variance of the number of ions, σ = ⟨Nα2⟩–⟨Nα⟩2, is proportional to Nα, and for Nα → ∞. When mesoscopic inhomogeneities spontaneously appear in the fluid phase, however, the fluctuations of the local volume fraction, i.e., the local deviations of the volume fraction from its average value ζ̅α, cannot be neglected. The local volume fraction is larger or smaller than ζ̅α when the cluster enters or leaves the considered region, and the states with ζα(r) larger or smaller than ζ̅α give more important contributions to the grand potential of the whole system as inhomogeneity gets stronger. The strength of the inhomogeneity can be measured by the variance ⟨(ζα(r) – ζ̅α)2⟩ of the local volume fraction. The linear extent of the inhomogeneity is determined by the wavelength of the density waves that gives the largest decrease of the system energy compared to the homogeneous state and depends on the interaction potentials.

Let us first consider our system with a particular distribution of the local volume fractions and with frozen fluctuations. Such a system for a given ζα(r) represents a set of subsystems such that the exchange of the ions between different subsystems is forbidden, and the grand potential is Ωco = −kBT ln Ξco, with Ξco denoting the summation of exp(−βH), where H is the microscopic Hamiltonian, over the microscopic states compatible with ζα(r). Because of the frozen mesoscopic fluctuations, our “locally canonical” ensemble can be reasonably well described within the MF approximation.

The grand thermodynamic potential in the presence of the above mesoscopic constraints (frozen mesoscopic fluctuations) can be written in the formwhere Uco, S, and μα are the internal energy, the entropy, and the chemical potential of the species α, respectively. Hereafter, the summation convention for repeated indices is used. We make the approximation −TS = ∫drfh(ζ+(r), ζ–(r)), where fh(ζ+(r), ζ–(r)) is the free-energy density of the hard-core reference system in the local-density approximationwhere the first two terms come from the entropy of mixing, and fhs is the contribution to the free energy density associated with packing of hard spheres with two different diameters. The expression for fhs in the Carnahan–Starling approximation[44] is given in the Supporting Information. Since we are interested in correlations between local fluctuations in regions separated by large distances, more accurate forms of the contribution to the entropy associated with packing of hard spheres, such as in ref (45), are an unnecessary complication.

In MF, the internal energy Uco is given by the expressionBecause ζα = πρασα3/6 is used in the above definition, we have rescaled the interaction potential, Vαβ(r) = Uαβ(r)/(vαvβ), where vα = πσα3/6.

When the constraints imposed on the microscopic states by ζ̅α are released and the ions can move from one subsystem to the other, the microscopic states incompatible with ζ̅α can appear. The average volume fractions can remain equal to ζ̅α when the fluctuations Δζα(r) = ζα(r) – ζ̅α cancel one another, but the grand potential contains the fluctuation contribution and has the form[43]In eq , ∫DΔζα denotes the functional integral over the mesoscopic fluctuations Δζα that according to the Landau theory[46] appear with the probability proportional to the Boltzmann factor exp(−βHfluc), withdenoting the change of Ωco when the fluctuation Δζα appears. Thus, in our theory the contributions to the grand potential from the microscopic and the mesoscopic degrees of freedom are given in the first and second terms in eq , respectively. In MF, i.e., with frozen mesoscopic fluctuations, the second term on the RHS of eq is neglected.

We are interested in the correlation functions in the disordered phasewhere the averaging is with the probability proportional to exp(−βHfluc). The matrix G with the elements defined in eq satisfies the analogue of the Ornstein–Zernike equation, G = C–1, where the inverse correlation functions C̃αβ are the second functional derivatives of βΩ[ζ+, ζ–] with respect to ζα and ζβ.[43] In the lowest-order nontrivial approximation beyond MF[38,43]where f̃(k) denotes the function f in Fourier representation. In the above equationwith α = +, −. Note that in this approximation the dependence of C̃αβ(k) on k comes only from βṼαβ(k). The last term in eq is the fluctuation contribution and comes from the last term in eq . In the Brazovskii-type approximation[47]Note that is the local variance of ζα; i.e., is the standard deviation from the space-averaged value of the local volume fraction of the α-ions. The larger is, the stronger the mesoscopic inhomogeneity and the less accurate the MF approximation. Equations –12 have to be solved self-consistently. In general, it is a nontrivial task.

We focus on the disordered inhomogeneous phase and assume that the inhomogeneity occurs on a well-defined length scale corresponding to the largest decrease of the energy and is quite strong, as in the case of ref (4). Because of this assumption, our considerations from now on are limited to the thermodynamic states where strong inhomogeneity at well-defined length scales is present. For weaker inhomogeneity, the set of the integrals, eqs –12, and G = C–1 must be solved. In our case, the peak of G̃γδ(k) (proportional to the structure factor) is high and narrow. For functions with a high, narrow peak, the main contribution to the integral comes from the vicinity of the maximum. We assume that the maximum of all the integrands in eq is very close to the minimum at k = k0 of det C̃(k), and we make the approximationwhere [C̃αα(k)] = C̃ββ(k) and [C̃αβ(k)] = −C̃αβ(k) for α ≠ β andNear the minimum at k0, we have the approximationwhere D0 = det C̃(k0) and βW″(k0) is the second-order derivative of det C̃(k) with respect to the wavenumber k at k = k0. From the approximations 15 and 14, we obtain[43,48]With all the above assumptions, the problem reduces to determination of the minimum of det C̃(k) and to a solution of three algebraic equations for C̃αβ(k0) (see eqs 10 and 13) becausewhereThe last approximation is valid for k ≈ k0.

It should be noted that the results obtained within the framework of this theory for several models of inhomogeneos mixtures were verified by simulations.[38,39,49]

Results

MF Approximation

In MF, we neglect the last term in eq and easily obtain explicit expressions for the matrix C̃MF(k) inverse to the matrix of correlations. These expressions are shown in the Supporting Information. In MF, the disordered phase becomes unstable with respect to oscillatory modulations of the volume fractions of the ions at the so-called λ-line on the (ζ̅,T*) diagram. The λ-line marks the boundary of stability of the disordered phase with respect to mesoscopic fluctuations of the volume fractions and separates the phase space into regions corresponding to the homogeneous and inhomogeneous (on the mesoscopic length scale) phases. Thus, in MF the λ-line is interpreted as a continuous order–disorder transition. In order to get the λ-line, one should solve the system of equationsWhen the λ-line is crossed, the inhomogeneities in the distribution of ions occur on the length scale 2π/k0. In our model, k0 depends in particular on ϵ++* that we left as a free parameter. In order to fit ϵ++* to our water-in-salt system, we need to have 2π/k0 ≈ 4 in a units (k0 ≈ 1.6 in a–1 units), for the molarity M ∼ 3–5 for which the experimental data were obtained. The volume fraction of the spherical ions with σ+ = 0.2 nm and σ– = 0.6 nm is related to the molarity M by , where NA is the Avogadro number. We get ζ̅ ≈ 0.27 and ζ̅ ≈ 0.32 for the 3.8 M and 4.6 M systems, respectively. However, the above formula is a very rough estimation for ζ̅/M in view of the strong dependence of ζ̅ on the diameter of the ions and the ellipsoidal shape of the TFSI– anions, and it only gives the order of magnitude of M in the experimental system for the given ζ̅ in our theory. Thus, in our semiquantitative analysis, we will consider volume fractions up to ζ̅ = 0.55.

The plot of k0 as a function of ϵ++* for ϵ––* = 1, δ = 0.5, and ζ̅ = 0.45 is shown in Figure . We can see that for ϵ++* = 5 the length scale of inhomogeneities is 2π/k0 ≈ 3.9 (in a-units), which for a = 0.4 nm gives 1.56 nm, which is close to the experimental result of 1.4 nm. We thus choose ϵ++* = 5 in our further calculations.

Figure 2

Wavenumber of the density waves k0 (in a–1 units) as a function of ϵ++* for ϵ––*= 1, δ = 0.5, and the volume fraction of ions ζ̅ = 0.45. ϵ* describes the strength of the non-Coulombic interactions (see eq ), and the size asymmetry δ is defined in eq . For the considered system, a = (σ+ + σ–)/2 ≈ 0.4nm.

Figure shows the λ-line in the ζ̅—T* (panel a) and ζ̅—k0 (panel b) coordinates for the model with δ = 0.5, ϵ++* = 5, and ϵ––* = 1. For the thermodynamic states below the λ-line, the waves with the wavelength 2π/k0 are more probable than the constant volume fractions. For our model, the emergence of the inhomogeneous structure for det C̃MF(k0) < 0 may be associated with the formation of aggregates such as clusters or layers, rather than with a phase transition. For such thermodynamic states, the fluctuations dominating on the mesoscopic length scale should be taken into account in order to restore the stability of the disordered phase. As found in different systems with spontaneously appearing mesoscopic inhomogeneities, fluctuations induce a change of the continuous transition found in MF to the first-order crystallization that is also shifted to higher volume fractions. Because of the instability of the disordered phase for T* < Tλ* in MF, the asymptotic decay of the correlation functions GαβMF(r) can be analyzed only for T* > Tλ*.

Figure 3

λ-line in ζ̅—T* (panel a) and ζ̅—k0 (panel b) coordinates for the model with δ = 0.5, ϵ++*=5.0, and ϵ––*=1.0. T* and δ are defined in eq and eq , respectively. ζ̅ = ζ̅++ζ̅– is the total volume fraction of ions, ζ̅α = πρασα3/6, k0 is in the a–1 units, a = (σ+ + σ–)/2 ≈ 0.4nm.

In general, the long-range behavior (r ≫ 1) of Gαβ(r) is described by the function[50]In eq , α0 = 1/λs and α1 = 2π/λ are the imaginary and real parts of the leading order pole of G̃αβ(q) in the complex q-plane, which is determined as the complex root q = iα0 ± α1 of the equation det C̃(q) = 0 having the smallest imaginary part. Since all G̃αβ(q) have a common denominator det C̃(q), they exhibit the same pole structure and have the same exponential contributions. Only the amplitudes and the phases θαβ differ for different αβ combinations.

We calculate α0 and α1 for our model in MF from the equation det C̃MF(q) = 0. The ζ̅-dependence of both the decay length λs = α0–1 and the period of oscillations λ = 2π/α1 of the correlation functions Gαβ(r) is presented in Figure for two values of the reduced temperature, T* = 1.5 and T* = 1.6 (T* > Tλ*). Note that the decay length α0–1 tends to ∞ (α0 → 0), and simultaneously the period of oscillations λ tends to 2π/k0 ≈ 4a when T* → Tλ*. More precisely, we get λ = 3.96 a for T* = 1.5, ζ̅ = 0.5 and λ = 3.98 a for T* = 1.6, ζ̅ = 0.55.

Figure 4

Model with δ = 0.5, ϵ++* = 5.0, and ϵ––* = 1.0. The decay length λs = α0–1 and the period of oscillations λ = 2π/α1 of the correlation functions Gαβ(r) as a function of the total volume fraction of ions ζ̅ for T* = 1.5 (solid line) and T* = 1.6 (dash-dotted line) in the MF approximation. α0 and α1 are in the a–1 units, with a = (σ+ + σ–)/2 ≈ 0.4 nm.

There is a very small difference between the values of λ obtained for the two temperatures, and the difference decreases with an increase of ζ̅. Moreover, for ζ̅ > 0.3 the dependence of λ on ζ̅ is weak, as in ref (4).

The values of the reduced temperature T* > Tλ* for large ζ̅, however, are too high when compared to room temperature. A rough estimate of the reduced temperature that corresponds to the conditions for the LiTFSI salt in water at room temperature (T = 300 °C and ε = 80) is about 0.5. Assuming that the dielectric constant of bulk water is decreased proportionally to the ion concentration, we should consider T* < 0.5.

Beyond the MF Approximation

In order to calculate the fluctuation contribution to the inverse correlation functions C̃αβ(k) in the Brazovskii-type approximation, we take into account the last term in eq 10 and solve the closed set of four equations for the unknowns k0 and C̃αβ(k0). The explicit forms of these equations are given in the Supporting Information. Once k0 and C̃αβ(k0) are determined, the inverse correlation functions C̃αβ(k) can be obtained from eqs and 17.

From G = C–1, one can calculate the correlation functions in Fourier representation. In Figure , we show the reduced correlation functions in Fourier representation G̃αβ*(k) = G̃αβ(k)/ζ̅αζ̅β for the fixed total volume fraction ζ̅ = 0.55 and for two temperatures, T* = 0.25 (panel a) and T* = 0.2 (panel b). The three correlation functions G̃αβ*(k) have sharp maxima for k = k0 ≃ 1.58. For both temperatures, the height of the G̃++*(k) maximum is about 100 times higher than the maximum of G̃––*(k). It should be noted that the dependence of k0 on T* for the fixed ζ̅ is negligible, especially in the range T* = 0.2–0.3 (see Table ).

Figure 5

Correlation functions in Fourier representation with the effect of fluctuations taken into account for T* = 0.25, ζ̅ = 0.55 (panel a) and for T* = 0.2, ζ̅ = 0.55 (panel b). G̃αβ* = G̃αβ/ζ̅αζ̅β, ζ̅α = πρασα3/6, and the wavenumber k is in a–1 units with a = (σ+ + σ–)/2. The results are for the model with δ = 0.5, ϵ++* = 5, and ϵ++* = 1. In the insets, we show sharp peaks of G̃––*(k).

Table 1

Wave Number k0, the Decay Length λs = α0–1, and the Period of Oscillations λ = 2π/α1 of the Pair Correlation Functions Gαβ(r) Depending on the Total Number Density ζ̅ for Fixed Values of Temperature T*a

T*	ζ̅	k0	α0	α1	α0–1	2π/α1
0.4	0.45	1.597	0.225	1.613	4.444	3.895
0.4	0.5	1.591	0.163	1.599	6.151	3.930
0.4	0.55	1.582	0.119	1.587	8.424	3.959
0.3	0.45	1.602	0.124	1.607	8.042	3.910
0.3	0.5	1.593	0.088	1.595	11.426	3.939
0.3	0.55	1.584	0.063	1.585	15.796	3.964
0.25	0.45	1.605	0.084	1.607	11.852	3.911
0.25	0.5	1.594	0.059	1.595	16.920	3.939
0.25	0.55	1.584	0.043	1.585	23.411	3.964
0.2	0.45	1.607	0.053	1.608	19.016	3.907
0.2	0.5	1.595	0.037	1.596	27.191	3.937
0.2	0.55	1.585	0.027	1.585	37.610	3.963

T* is defined in eq , ζ̅ = ζ̅+ + ζ̅–, ζ̅α = πρασα3/6, and k0, α0, and α1 are in a–1 units.

The reduced correlation functions in real-space representation, Gαβ*(r) = Gαβ(r)/ζ̅αζ̅β, are obtained from the inverse Fourier transformation of G̃αβ*(k). They are shown in Figure for T* = 0.25, ζ̅ = 0.55 (panel a) and for T* = 0.2, ζ̅ = 0.55 (panel b). As seen, Gαβ*(r) shows exponentially damped oscillatory behavior, as described by eq . The period of damped oscillations is about 4a. We study the asymptotic decay of the correlation functions Gαβ(r) using the pole analysis. The results of this numerical analysis for T* = 0.2, 0.25, 0.3, and 0.4 and for ζ̅ = 0.45, 0.5, and 0.55 are summarized in Table . For the fixed volume fraction, the period λ = 2π/α1 coincides with 2π/k0 and is rather kept constant for T* ≤ 0.3. For the fixed temperature, λ is a weakly increasing function of ζ. It should be noted that the period of damped oscillations obtained with the effect of fluctuations taken into account is very close to the period obtained in MF for the higher temperature. By contrast, the decay length α0–1 noticeably increases with an increase of ζ̅ for the fixed temperature, and this increase is more rapid for lower temperatures. In Figure , we present the decay length α0–1 as a function of the total volume fraction of ions ζ̅ for fixed temperatures (panel a) and as a function of the Bjerrum length lB for fixed volume fractions (panel b). One can observe that α0–1 has a nearly linear dependence on the volume fraction for fixed T*, with the slope decreasing with T* and a nearly linear dependence on lB for fixed ζ̅, with the slope increasing with ζ̅.

Figure 6

Correlation functions in real space with the effect of fluctuations taken into account for ζ̅ = 0.55, T* = 0.25 (panel a) and T* = 0.2 (panel b). Gαβ* = Gαβ/ζ̅αζ̅β, ζ̅α = πρασα3/6, and r is in a units with a = (σ+ + σ–)/2. The results are for the model with δ = 0.5, ϵ++* = 5, and ϵ++* = 1.

Figure 7

Decay length λs = α0–1 of the correlation functions Gαβ(r) as a function of the total volume fraction of ions ζ̅ for T* = 0.2, 0.25, 0.3, and 0.4 (from the top to the bottom line) (panel a) and as a function of the Bjerrum length lB for ζ̅ = 0.45, 0.5, and 0.55 (from the bottom to the top line) (panel b). 1/α0 and lb are in a units, a = (σ+ + σ–)/2 ≈ 0.4 nm. The results are for the model with δ = 0.5, ϵ++* = 5, and ϵ++* = 1.

Discussion

We have developed a highly simplified model for water-in-salt electrolytes and focused on the salt LiTFSI that was a subject of recent experiments.[3,4] In our model, the ions are treated as charged hard spheres with different diameters, and we assumed additional, water-mediated specific interactions between the ions of the same sign. We assumed that the solvent influences the distribution of the ions mainly by inducing effective interactions between them, and otherwise the solvent can be neglected. Next, we assumed that the detailed shape of the specific interactions is not important, as long as these interactions are strongly attractive but of a short range. We chose the Yukawa potentials for the specific interactions and adjusted the parameters to the LiTFSI by requiring the same scale of inhomogeneities as found experimentally. Importantly, the obtained sum of the Coulomb and specific interactions for both the anions and the cations is attractive at short distances and repulsive at large distances.

For this model, we calculated correlation functions for concentrated electrolytes for reduced temperatures close to room temperature, using our theory for binary mixtures with competing interactions.[38,39,43] There is no unique way of relating the volume fraction in the approximate theory to experimental molarity since we assumed spherical rather than ellipsoidal anions, and the values of the ion diameters in the theory are not precise. We considered volume fractions 0.45 ≤ ζ̅ ≤ 0.55 somewhat larger than in the experiment but of the same order of magnitude (the molarity 0.38–0.46 considered in ref (4) for spherical ions of our sizes gives 0.27 ≤ ζ̅ ≤ 0.32). We obtained exponentially damped oscillations for the correlation functions with the period λ ≈ 4 in a ≈ 0.4 nm units (in the experiment,[3,4] λ ≈ 1.4 nm), very weakly depending on ζ̅ in the considered range of concentrations. In MF, λ decreases a bit. When the fluctuations are taken into account, λ increases with ζ̅. Interestingly, λ slightly decreases, increases, or does not change with increasing ζ̅, depending on the method of determining it in the experiment.[4]

We found that the decay length of the correlations, λs, increases almost linearly with ζ̅ for fixed reduced temperature T*. It increases also with the Bjerrum length for temperatures of the order of room temperature. In particular, for the reduced temperature T* = 0.3 (Bjerrum length lB = 3.3 in a ≈ 0.4 nm units), we obtain λs ≈ 3.2–6 nm, and for T* = 0.2 (Bjerrum length lB = 5) we obtain λs ≈ 8–15 nm for ζ̅ = 0.45–0.55. These values are in semiquantitative agreement with the experimental decay lengths, λs = 8.3–11.5 nm, for the molarity 3.8–4.6 M (ζ̅ = 0.27–0.32).

The approximate relation λs ∼ lBρ was discovered earlier for concentrated simple salts in water, ILs, and alkali halide solutions[12] and confirmed theoretically for the RPM (charged hard spheres with equal diameters and positive and negative charges of the same magnitude).[23] Our present results show that this relation is valid when ions with significantly different sizes interact with additional strong short-range potentials as well. It means that the dependence of λs on lBρ does not depend on the length scale of inhomogeneity, at least in the range 2a – 4a.

The very weak dependence of λ on ζ̅ and the semiquantitative agreement of the decay length with the experimental results indicate that the properties of the correlation functions do not depend sensitively on the details of the interactions. The ellipsoidal anions are approximated by the spherical ones; implicit solvent-inducing effective anion–anion and cation–cation interactions are assumed; and we neglected fluctuations of the dielectric constant induced by the concentration fluctuations. The dependence of the reduced temperature on ζ̅ (through the dependence of ϵ on ζ̅) was disregarded as well. The latter dependence for the rather narrow range of ζ̅ is not very strong, however. Finally, we rather arbitrarily assumed the shape of the specific interactions. With the above simplifications, we got semiquantitative agreement with experiments. It means that the above features are important only on the quantitative level.

We conclude that the key property determining the inhomogeneities on the mesoscopic length scale is the shape of the sum of the Coulomb and the solvent-induced specific interactions. In order to induce mesoscopic inhomogeneities, this sum should be attractive at short and repulsive at large distances, with the ranges and strengths of the attractive and repulsive parts determined by the properties of the ions and the solvent. If these anion–anion and cation–cation potentials have a negative minimum followed by a positive maximum, then layers of ions of the same sign and of the thickness determined by the width and depth of the attractive well can be formed. We believe that this conclusion is not restricted to the particular case of the water-in-LiTFSI but applies to concentrated solutions of salts with significantly different solubility of cations and anions in charge-neutral solvents, such that the effective interactions have the form similar to the one shown in Figure . Hydrophilic cations and hydrophobic anions are present in the so-called “antagonistic salts”, such as for example sodium tetraphenylborate (NaBPh4). This salt induces strong mesoscopic inhomogenity on the length scale 10 nm in the D2O/3-methylpyridine (3MP) mixture.[51−54] In concentrated solutions in water of this and some other antagonistic salts, mesoscopic charged regions correlated over large distances may be present. On the other hand, we expect that for a more complex anion structure than the structure of TFSI– the approximation of the anion shape by a sphere cannot be sufficient that it, in turn, can worsen the agreement between theoretical and experimental results. In particular, more elongated shape may lead to orientational ordering of the anions.

It is worth noting that our theory is formulated for the case when the valences of the cation and the anion are equal in magnitude but not necessarily equal to 1. Our theory can be extended to the case of unequal valences by changing the parameters in the Coulomb potentials 2–4. In principle, the theory could be extended to lithium salt/IL mixtures, but in a minimal model at least two types of cations and one type of anion, i.e., a three-component mixture, should be considered.

Finally, we have shown that the self-consistent theory with the local fluctuations of the volume fractions taken into account[38,39] can predict the structure with local inhomogeneities on a semiquantitative level.

There remains one unsolved problem—namely, the experimental disjoining pressure between crossed mica cylinders decays monotonically at large distances,[4] whereas our theory predicts the oscillatory decay. The same problem concerns simple salts and some other ILs modeled by the RPM.