Reversible ionic aggregation kinetics in concentrated electrolytes

Abstract

Here we develop a formalism for reversible ionic aggregation kinetics in an example concentrated electrolyte, building on previous equilibrium work of McEldrew and co-workers, and thermoreversible polymers and patchy particle systems. This is achieved through solving a macroscopic rate equation of open/occupied association sites, shown to be a solution of the reversible Smoluchowski aggregation equation, which predicts how ionic associations in electrolytes change subject to a step-change in conditions. We tested the derived equations against atomistic molecular dynamics simulations of a salt-in-ionic liquid, where good qualitative agreement is obtained, but quantitative differences are found. This highlights the multiple time scales that exist in concentrated electrolytes, with a fast timescale preceding a longer timescale. We hope this formalism opens new avenues in understanding the dynamics and non-equilibrium behaviour in electrolytes. For example, the formalism can be developed further to investigate the non-Newtonian behaviour of concentrated electrolytes, double layer charging, and the slow dynamics of these electrolytes in confinement.

---

1. Introduction

Concentrated electrolytes, such as water-in-salt electrolytes (WiSES),^1–7 ionic liquids (ILs),^8–11 salt-in-ILs (SiILs),^12–15 amongst many other examples,^16–19 are increasingly of interest in various technologies.^{9–11,17,20–26} For example, ILs have potential applications as ‘green’ solvents for electrochemical reactions at electrified interfaces, and electrolytes for supercapacitors owing to their low vapour pressure and flammability, and because of their wide electrochemical stability window (ESW).^8–11 Similarly, WiSEs are promising alternative electrolytes for batteries which are non-flammable, owing to water being the solvent, but still have a wide ESW from the high salt concentration.^2,27–30 Finally, SiILs, or salt-doped ILs, are also finding applications in energy technologies, where the unique properties of ILs are being exploited as the ‘solvent’ to dissolve active ions.^31–33

These concentrated electrolytes also pose interesting fundamental questions.^24,34,35 For example, in ILs,^36–38 WiSEs^39,40 and SiILs¹⁵ (and also other concentrated electrolytes⁴¹), the observation of “anomalous underscreening” (unexpectedly extremely long-ranged force decay lengths in surface force balance/apparatus measurements) raised questions about the electrostatic correlations in bulk, equilibrium concentrated electrolytes.^{15,36–39,41,42} While recent research has indicated the non-equilibrium nature of these measurements cannot be neglected,⁴³ these measurements motivated significant work into understanding bulk, equilibrium concentrated electrolytes.^44–49 Furthermore, WiSEs behave as non-Newtonian fluids,⁵⁰ and have nano-channels of ionic aggregates interpenetrated by water domains that result in facile transport of the active cations.^27–30,51 In the context of SiILs, negative transference numbers of alkali metal cations at low mole fractions generated much interest in studying the transport properties of concentrated electrolytes.^12,13

These properties have been understood, at least in part, from the ionic aggregates and solvation environments in the concentrated electrolytes.^24,34,35,52 For example, the negative transference numbers of SiILs can be explained through the formation of a negative shell of anions around the alkali metal cations.^12,13 Moreover, more generally, ionic associations and solvation is known to be important in understanding interphase formation in batteries^53–55 and the kinetics of electrocatalysis reactions,²³ for example. Typically, these conclusions have been drawn from correlating the coordination environments of the electrolytes (and from the relative numbers of solvent separated ion pairs, contact ion pairs and “aggregates”), and these observed properties. Such approaches, however, provide little fundamental insight into the underlying laws which link these coordination environments to other properties of the electrolytes.

There are, however, microscopic theories which directly link the ionic associations and solvation to physiochemical properties of these electrolytes.^56,57 A notable example is the thermoreversible ionic aggregation model developed by McEldrew et al.,^57–62 which links coordination environments to activity and redox potentials,²⁶ transference numbers,^59–61 electrical double layers^63–66 and many more properties.⁵⁷ This framework is founded on the classical polymer theories of Flory, Stockmayer and Tanaka,^67–78 which describe the thermoreversible bonding of monomer units into a polymer, where McEldrew et al.⁵⁷ translated this to the thermoreversible association of ions into aggregates (i.e., an alternating copolymer of cations and anions, decorated by solvent).

While McEldrew et al. have mainly investigated equilibrium and linear-response transport properties,^57,59–62 in the context of polymers and patchy particle systems,^79–85 non-equilibrium properties have also been investigated. For example, Sciortino, Tartaglia and co-workers connected the kinetics of patchy particle aggregation to the distribution of clusters,^{79–81,86–90} and Tanaka^91,92 used this to understand the viscoelastic response of polymers with thermoreversible associations. This presents itself as a significant opportunity, where concentrated electrolytes could greatly benefit from drawing on these fields of research. For example, the dynamics of ionic associations are routinely studied (through residence times) of electrolytes,^12,13 and it is known that the desolvation time scale is important for intercalation rates.^93,94 However, more quantitatively linking these to measurable electrolyte properties through a microscopic theory remains lacking.⁹⁵

Here, we take the first steps to apply the non-equilibrium formalisms from polymer and patchy particle systems to concentrated electrolytes. Specifically, we investigate the dynamics of reversible ionic aggregation in bulk SiILs. Analytical expressions for how the coordination numbers of the ions vary subject to a step-change in conditions are derived. A single decay time is analytically obtained in the theory, and found to depend on the rate of forming a single association and the final coordination numbers. To validate these derived equations, we perform molecular dynamics simulations. Overall, qualitative agreement between theory and simulations is found. The quantitative differences between the theory and simulations reveal a single rate constant is not sufficient to describe the ionic aggregation over all time scales. The formalism provides insight into the kinetics of reversible ionic aggregation in SiILs, and indicates avenues for further developments of the theory to more quantitatively link the dynamics of associations to physiochemical properties of concentrated electrolytes.

This paper is structured as follows. First is the Theory Section. In this section, we initially cover the background and assumptions of the formalism of McEldrew and co-workers.^57,59–62 Next we go into a review of the bulk theory, providing a high-level overview initially. The remainder of the Theory Section introduces the non-equilibrium formalism from patchy particle systems, but casts it into the context of concentrated electrolytes, and then we solve for how the coordination numbers change in time from a step-change in conditions. Second, we provide the details of the atomistic MD simulations performed to test the derived equations. Third is the Results Section, which is self-contained and compares the non-equilibrium MD simulation results for how coordination numbers change in time against the derived, analytical equations. Finally, we discuss the comparison between theory and simulation, provide details of how the presented formalism can be compared against experiments, and further extensions of this formalism to transport properties, interfacial properties and rheological properties is outlined.

2. Theory

Here we study salt-in-ionic liquids (SiILs) as an example system to investigate the kinetics of reversible ionic aggregation. In particular, we study alkali metal doped ILs, where the anion is the same in each salt (3 species in total). Previously, McEldrew et al. developed an equilibrium theory of the bulk ionic aggregates,⁶⁰ quantified the effective charge of the alkali metal cation, and studied the electrical double layer.^15,66 We recommend that the reader familiarities themselves with these previous works,^15,60,66 also for other electrolytes.^26,59,61,62 Nonetheless, we provide a brief summary of the equilibrium theory before progressing onto the kinetics of aggregate formation. Before progressing onto the bulk and non-equilibrium properties, we first outline the assumptions of the theory and possible limitations.

2.1. Assumptions and definitions

The alkali metal cations and anions are assumed to form a polydisperse mixture of Cayley-tree aggregates through the associations in their fist coordination shell.^57,60 These aggregates have an alternating cation–anion structure, a reflection of overscreening, as a way of introducing correlations beyond mean-field electrostatics.^46,63,96 The Cayley tree approximation of aggregates is central to the analytical tractability of the theory, where long-range order can be constructed from short-range associations.^57,60 The IL cations are assumed not to participate in ionic associations directly, but can still influence the aggregates by interacting with the open anion association sites through a regular solution term.^15,60,97 This is schematically shown in Fig. 1. All associations and interactions are taken to be independent. We neglect all further interactions beyond these. In the absence of the IL cations, we treat the electrolyte as a ‘non-interacting’ mixture of these polydisperse aggregates.

Fig. 1 Schematic of non-equilibrium kinetics of ionic aggregation in the studied salt-in-ionic liquid (SIIL). The alkali metal cation and IL anion bind together to form aggregates. The alkali metal cation and IL anions are shown with 3 dangling bonds, to represent they can form associations, while the IL cations are shown without, since they do not associate strongly with the anions. The aggregates formed are denoted by dotted lines surrounding them, while the interactions between the IL cations and anions are shown with a bolt. From step-changing a property of the system (initially the interactions are "turned off", but then reinstated), we investigate how the aggregates evolve in time.

We treat the SiIL as an incompressible lattice fluid.^57,60 The alkali metal cations are denoted by +, and IL cations denoted by ⊕, and anions by −. The total volume fractions of ions, ϕ_j, are known, with 1 = ϕ₊ + ϕ₋ + ϕ_⊕ from the incompressibility constraint. The lattice site is assumed to be the volume of the alkali metal cation, v₊. The volumes occupied by the other ions are expressed relative to this volume through v₊ξ_j, where ξ_j = v_j/v₊. The (dimensionless) concentration of each ion per lattice site is given by c_j = ϕ_j/ξ_j.^57,60

The alkali metal cations can form a maximum of f₊ associations, and the anions a maximum of f₋ associations. This is referred to as the functionality.^57,60 These association sties, shown as dangling bonds in Fig. 1, represent associations between species (ions), which means that mono-dentate associations and bi-dentate associations, etc., are treated as a single association since they are between the same species.⁶⁰ As these functionalities are larger than 1, a polydisperse cluster distribution can form with clusters of rank lm, where l is the number of cations and m is the number of anions, with dimensionless concentration c_lm.^57,60 For functionalities equal to or larger than 2, a percolating ionic network can emerge.⁵⁷ This is referred to as the gel here. In the gel regime, we employ Flory's convention to determine the volume fractions of ions in the sol (ϕ^sol_±, i.e., not the gel phase) and gel phase (ϕ^gel_±), where ϕ_± = ϕ^sol_± + ϕ^gel_±. The total dimensionless concentration of alkali metal cations is given and anions . Note also helps convert between volume fractions and dimensionless concentrations.

2.2. Equilibrium theory overview

Before getting into the details of the theory, we outline the structure of the derivation and concepts introduced. The assumed free energy functional contains all the information necessary to derive the cluster distribution, c_lm, which is the central quantity in the theory. From taking functional derivative of this free energy, and establishing equilibrium between the resulting chemical potentials, we can derive an equilibrium equation for ϕ_lm. From simplifying this, we can arrive at c_lm. As c_lm is a function of the concentration of free species, which we wanted to predict, not be an input, we introduce the concept of association probabilities, and another set of mass action laws for these variables, which allows the full system of equations to be determined, just based on the association constant.

The free energy functional^15,57,60,66 is given by

where β⁻¹ = k_BT is thermal energy. The first term is the ideal entropy of the IL cations, and the second term is the ideal entropy from all clusters. The third term is the free energy of forming clusters, where Δ_lm is the free energy of forming each cluster of rank lm. The fourth term is the regular solution interaction between the IL cations and open anion binding sites,^15,60,66 with strength χ. The fifth term comes from the free energy of ions associating to the gel, Δ^gel_j.

The chemical potential of a rank lm cluster, µ_lm, can be determined by the functional derivative of the free energy with respect to c_lm. Establishing equilibrium between clusters of rank lm and their free counterparts, we have

lµ₁₀ + mµ₀₁ = µ_lm. (2)

Using the expressions for chemical potential, we obtain the cluster equilibrium

ϕ_lm = K_lmϕ₁₀^lϕ₀₁^m, (3)

where K_lm = exp{(l + m − 1)(1 + βχϕ_⊕) − βΔ_lm} is the equilibrium constant. To simplify this expression we have to specify Δ_lm, which is where our assumed Cayley-tree aggregates enter. We consider the free energy of forming a cluster of rank lm, Δ_lm,^57,60 to have three contributions

Δ_lm = Δ^comb_lm + Δ^conf_lm + Δ^bind_lm, (4)

where Δ^comb_lm is the combinatorial entropy (entropy related to the number of ways of arranging those ions in that Cayley-tree aggregate), Δ^conf_lm is the configurational entropy (related to those aggregates existing on the lattice) and Δ^bind_lm is the binding energy. The combinatorial entropy is given bywhere W_lm is related to the enumeration of the ways that a rank lm cluster can be formed. For Cayley trees, Stockmayer^68–70 determined the expression to beThe binding energy is approximated as the energy of a single association, Δu, multiplied by the total number of associations in a cluster (as they are assumed to be independent), which for Cayley tree clusters is l + m − 1. Thus, the binding energy is

Δ^bind_lm = (l + m − 1)Δu. (7)

The configurational entropy describes the entropy of placing a rank lm cluster on the lattice. We modify Flory's expression^{57,67,72,77,98} to bewhere Δs is a per-association entropy change.²⁶

From using Δ_lm, we obtain the thermodynamically consistent cluster distribution

where λ is the association constant given by

λ = exp{β(−Δu + TΔs + χϕ_⊕)}. (10)

The cluster distribution c_lm has direct utility, unlike eqn (3), since it essentially only depends on one parameter, the association constant, λ, which needs to be determined.

In eqn (9), c_lm is written in terms of the volume fraction of free alkali metal cations (ϕ₁₀) and anions (ϕ₀₁). These quantities are, however, not a priori known, and moreover, we aimed to predict them from the theory, not have them as inputs for it. To overcome this, we^{57,60,71–78} express ϕ_10/01 in terms of ϕ_± and introduce ion association probabilities, p_ij, which is the probability that an association site of species i is bound to species j. Simply put, these association probabilities are the coordination numbers over the functionalities, i.e., the actual number of associating species over the maximal number. Therefore, the volume fraction of free alkali cations can be written as ϕ₁₀ = ϕ₊(1 − p₊₋)^f+ and free anions as ϕ₀₁ = ϕ₋(1 − p₋₊)^f−, assuming independent associations.

The association probabilities can be determined through the mass action law

and conservation of associations

p₊₋ψ₊ = p₋₊ψ₋ = ζ, (12)

where ψ_i = f_iϕ_i/ξ_i is the total dimensionless concentration of association sites (per lattice).^57,60 Eqn (11) and (12) permit an explicit solution for the probabilitieswhere ψ_± = ψ₊ + ψ₋. It has been shown, from a bond percolation analysis, that the percolating ionic network onsets when p₊₋p₋₊ = [(f₊ − 1)(f₋ − 1)]⁻¹, which we also call the gel-point (see ref. 57 and 60). To determine sol and gel properties in the post-gel regime, we use Flory's treatment.⁵⁷

With the introduced association probabilities, the cluster distribution mainly depends on the value of the association constant, λ. This can be taken as a free parameter, but λ can also be constrained from MD simulations or Raman/IR measurements from using eqn (11) and (12) with the experimental coordination numbers.²⁶ For this to be possible, the functionalities of each species also needs to be known. Again, this can be directly inferred from MD simulations, and so these are also not free fitting parameters.^59–62 The only other parameters of the model are the volumes of species, but these can be determined from simple calculations.^59–62 If the free energy of forming a single association, Δf = Δu − TΔs, and the regular solution term, χ, need to be determined, the composition of the SiILs must be varied to extract λ′s composition dependence.⁶⁰ Therefore, while this is a phenomenological theory with parameters to tune for specific chemistries, none of these parameters are free, and so the predictions of the theory should be robust, provided its assumptions are being adhered to. The main assumption of the theory that can be broken is the Cayley-tree approximation. In SiILs, loop formation has been found at equilibrium, which causes disagreements between the theory and MD simulations.^15,60

2.3. Reversible aggregation kinetics

Next we outline the non-equilibrium theory for reversible ionic aggregation kinetics in bulk SiILs. This has not been explored in the context of concentrated electrolytes before. However, in patchy-particle systems, an ostensibly similar system, the kinetics of aggregation has been well studied.^{79,81,86–90} In this section, the presented formalism closely follows to the works of Sciortino, Tartaglia and co-workers,^{79,81,86–90} and Dongen and Ernst.⁹⁹ First we write down the considered aggregation equilibrium which are the reactions that go into the reversible Smoluchowski equation. The remainder of this section then specifies the rate kernels of these reversible aggregation and fragmentation processes.

The reversible Smoluchowski aggregation equation governs the time evolution of thermoreversible clustering. This master equation describes how the change in concentration of an aggregate of rank lm, in the two-component case studied here, with concentration c_lm(t), depends on the rates of formation and destruction of this aggregate

The first reversible aggregation process, eqn (14), describes the rate of formation of c_lm through the binding between and with a rate of ^90,99 Note we have left the rate constant to generally depend on the values of l, l′, m and m′, with exactly how the rate depends on these values being shown later. In this same equilibrium, c_lm can be removed through breaking an association to form and with rate . While the second reversible aggregation reaction, eqn (15), describes the removal of c_lm through the formation of an association with to form , with a rate of , again generally given. In addition, there is the formation of c_lm from the dissociation of , with rate ^90,99 These reversible aggregation reactions can be seen in the reversible Smoluchowski aggregation equationwhere the subscripts of and were dropped for clarity, and the explicit time-dependence of the concentrations suppressed.⁹⁰

To solve the reversible Smoluchowski aggregation equation, we now must specify the rate kernels of the processes. First we can introduce the elementary rates for the formation of a single association, k, and the rate of a single dissociation, f, which is related to the association constant at equilibrium through λ = k/f.^90,99 Following ref. 90 and 99, the coalescence rate for two aggregates is related to the number of dangling association sites in each aggregate and the individual rate of forming associations. For a cluster of rank lm combining with a cluster of rank l′m′, this is

where, for example

σ_lm = f₊l + fm − 2(l + m − 1) = (f₊ − 2)l + (f₋ − 2)m + 2, (18)

is the number of open association sites in the aggregate of rank lm. For a free cation and free anion associating, the rate of this process is not simply k, but K_1,0,0,1^(r) = kf₊f₋

At equilibrium (t = ∞), the steady state solution of eqn (15) imposes⁹⁰

Therefore, from eqn (19), the full fragmentation rate is related to the rate of an individual fragmentation throughThe same procedure can be repeated for eqn (14), although not shown here.^90,99

In principle, this determines the reversible Smoluchowski aggregation equation to be solved.^90,99 Directly solving the reversible Smoluchowski aggregation equation is a formidable task, however. Even in the irreversible Smoluchowski equation, direct solutions are limited to simple initial conditions and certain choices of the rate kernels. However, Flory and Stockmayer, and later others, demonstrated that the association probabilities could be used as the central variable, greatly simplifying the problem.^90,99

Instead, the time dependence of the association probabilities is solved for, considering the addition and removal of single associations through this ensemble averaged quantity. It was shown that these time-dependent association probabilities used in the equilibrium cluster distribution, eqn (9), is a solution of the reversible Smoluchowski aggregation equation, in a quasi-equilibrium way.^90,99 In what follows, we solve the time-dependence of the association probabilities under different approximations. In the main text, we solve a simpler “symmetric case” where there is only one association probability,^90,99 and in the SI the more general “asymmetric” case of two association probabilities is shown.

2.4. Symmetric case

Typically, the functionalities of the alkali metal cation and anion are different in SiILs, which generally means their association probabilities are not equal.⁶⁰ For the association probabilities to be equal, we must have ψ₊ = ψ₋ (from the conservation of associations), which can (sometimes) be achieved at a specific composition in SiILs. Since f₊ = 4 and f₋ = 3,¹⁵ to achieve p₊₋ = p₋₊, there must be 3/4 times the number of alkali metal cations as anions. In this section, we consider p₊₋ = p₋₊ = p generally and do not specify exact values of f_i or ϕ_i explicitly, other than they satisfy ψ₊ = ψ₋ = ψ to ensure a single association probability.

As demonstrated in ref. 90 and 99, the time-dependence of p(t) is determined from the macroscopic rate equation

where [A^R] = ψp is the concentrated of bound association sites, and [A^F] = ψ(1 − p) is the concentration of unbound association sites. The rate of forming a single association, k, pre-multiplies [A^F]², to obtain the overall rate of forming associations. While the rate of breaking a single association, f, pre-multiplies [A^R], to obtain the overall rate of breaking associations.^90,99 Similar to the equilibrium theory, this macroscopic rate equation assumes independence and equal reactivities of each association. Thus, the first-order differential which describes the time-dependence of the association probability is given byIn the SI, we show how this equation can be derived from considering a limiting case of ion pair formation.⁶⁴ In the limit of long-time, dp(∞)/dt = 0, we findwhich is simply the mass action law, as in eqn (11), where the conservation of associations has been applied. Since the elemental rates, and therefore, the association constant are determined by the final state, the equilibrium state is naturally incorporated into the differential equation. Before there was a step-change in conditions that results in a variation of the associations, different rate constants and equilibrium constants applied. Upon the step-change, these rates are forgotten, but they persist through the initial values of the association probabilities.^90,99

Applying the initial condition of p(0) = p⁰, we find the solution to be^90,99

where p(∞) = p^∞ has been introduced, and we have defined the following parameterswhich is a statement of the direction that the association probabilities are changing (if they are increasing or decreasing) and by how much, andis the inverse time-scale that controls the evolution of the association probabilities.^90,99 We find that this inverse relaxation time depends on two main factors. Firstly, it is proportional to the rate of forming an association.^90,99 Therefore, different electrolytes could have different relaxation times based on this absolute scale, which is information that is removed in the association constant, since that depends on the ratio of the rates.^12,13,61 The second term is proportional to a function of the probabilities, which cause a divergent time scale when the association probabilities approach 1. Therefore, in principle, extremely long relaxation times of electrolytes could be obtained from reversible aggregation. Note that even though there are many values of the rate kernels from the different aggregates, these all collapse into a single exponential with one time constant, τ.

The solution of p(t), in eqn (24), then completely determines the changes in the aggregates in the system, which has been shown to be a solution of the reversible Smoluchowski aggregation equation. Therefore, the aggregates change in a quasi-equilibrium way, i.e., they always have a distribution that obeys an equilibrium cluster distribution [eqn (9) with the association constant which changes in time, calculated using the association probabilities], with their changes being completely described through p(t). This was noted by Sciortino and co-workers^79,81,88 as a connection between time and temperature through p, i.e., there is a time in the non-equilibrium aggregation when p(t) = p^∞(T), which can be thought as the temperature changing with time throughout the simulation.

In the post-gel regime, eqn (24) still applies, but the total association probabilities can be decomposed into sol and gel contributions, using Flory's post-gel treatment.^57,60 For more information on how the reversible Smoluchowski aggregation equation is modified, see ref. 99.

We can also study various limits of eqn (22). An example is irreversible aggregation, k ≫ f and λ = ∞, described by

Starting with an initial condition of p(0) = 0, the solution to this equation is p(t) = kψt/(1 + kψt). This can be seen to be a limit of eqn (24) from taking p(0) = 0 with λ ≈ ∞, i.e., p^∞ ≈ 1.

If, by contrast, we study irreversible fragmentation, where the initial aggregates can only break apart, we would have

The solution of this is simply p(t) = p⁰e^−ft. Therefore, the aggregates should exponentially decay, with the rate being solely determined by the rate of breaking a single association. This can be seen to be a limit of eqn (24)–(26), through taking p^∞ to be zero and λ = 0. Commonly, the lifetime of associations (residence times) are calculated from MD simulations, and this differential equation would describe it, if there was a single time-scale.^12,100,101

In this section and the previous one, it should be noted that the dynamical version of this theory introduces a new parameter, namely the rate of forming k (or breaking f) an association. Only a single parameter is introduced, since the rate of breaking (forming) an association is related through the association constant. This parameter needs to be determined in the theory from comparison against MD simulations or experiments. In principle, it could be extracted from residence times of equilibrium simulations, and applied to the non-equilibrium case, or directly fit from simulations/experiments.

3. Molecular dynamics

To test the equations derived in the Theory Section, we perform atomistic, classical molecular dynamics simulations of an example SiIL. We choose to study NaTFSI in EMIMTFSI. To perform these simulations, we use LAMMPS¹⁰² and the CL&P force field¹⁰³ with charges rescaled by 0.7 to mimic the effect of polarization.^12,13 The initial inputs were generated with fftool and packmol.¹⁰⁴ We studied mole fractions of x = 0.25, 0.5, 2/3, 0.75 in NaTFSI_xEMIMTFSI_1−x, with 400 TFSI⁻ anions in the simulation cell and varying numbers of each cation. A timestep of 1 fs was utilized, with all hydrogen atoms fixed by shake.

To equilibrate the system, we initially performed a temperature annealing cycle. Starting from a 1 ns run at 300 K (for example, as different temperatures are investigated, this is the final temperature is initiated here), we ramped the system to 450 K over 1 ns, and held at that temperature for another 1 ns, before the same procedure is repeated with a final temperature of 600 K (for higher temperature simulations, e.g., 400 K and 500 K, the same change in temperatures were performed for the annealing cycle). An equivalent cooling cycling was then applied. Then another 10 ns equilibration was performed at the target temperature. This was performed in NPT, where we use a temperature thermostat with 100 fs, and 1000 fs for the pressure. In NVT, at the equilibrium density, we then equilibrated for another 10 ns.

To model the effect of the non-equilibrium kinetics of aggregation, we need to step-change a property/variable in the simulation. The most natural way of introducing non-equilibrium conditions is to step-change the temperature. Several simulations were performed, where, for example, we first equilibrated the electrolyte with a density that corresponds to 300 K (the final, target temperature) at a temperature of 500 K, which causes the associations to break from the higher temperature. Then we would step change the temperature to 300 K, i.e., the temperature of the density that the box is equilibrated to in NVT. This, however, only produces very modest changes in the coordination numbers. The results of these simulations are shown in the SI, but will be discussed in the main text.

To produce a more substantial effect, to investigate more thoroughly if the kinetic equations derived perform well, we investigated a more artificial test. In MD simulations, we have control over the interactions in the electrolytes. Therefore, we can turn off and on the electrostatic interactions, which is a main driving force for associations in the electrolyte. After the NVT equilibration at the target temperature with all interactions, as described before, we rescaled all charges by 10⁶, making them essentially 0, and equilibrated the electrolyte for 1 ns. Then the original charges (0.7 of the CL&P force field¹⁰³) of the atoms were reinstated, and the changes in coordination numbers tracked in time. Note that turning on and off the charges pumps/removes significant energy into the system. Therefore, while tracking how coordination numbers change in time, we more aggressively thermostated the electrolyte (2 fs for 1000 steps, then 5 fs for 1000 steps, and 10 fs for 1000 steps before the normal value was used for the remainder). This ensured the temperature remained approximately constant (if this was not performed, temperatures in the range of 10³ K were found). Note we did not find this aggressive thermostating qualitatively changed the results, with only a small quantitative effect observed.

The association probabilities defined in eqn (13) are directly linked to the average coordination number of alkali cations by anions (CN₊₋ = f₊p₊₋), as well the coordination number of anions by alkali cations (CN₋₊ = f₋p₋₊). From the MD simulations, we define if a Na⁺ and TFSI⁻ are associated if the O of TFSI⁻ is within 3.3 Å of Na⁺, or vice versa.^12,13,60 This was chosen based on the first time g(r) reaches 1 after the first peak,⁶¹ see SI for details. Other real-space cutoffs were investigated too. While the exact quantitative numbers change, the overall conclusions found here remain independent of this choice. Other methods for defining associations exist, such as kinetic criteria,¹⁰⁵ but as there are large energy and dynamical changes occurring in the non-equilibrium simulations, the real-space method appears to be the most robust. Note that if more than 1 O in TFSI⁻ is within the cutoff to the same Na⁺, this just counts as 1 association, i.e., bidendate associations only count as a single association here, as we are counting based on species. To quantify the ionic aggregates in the simulations, i.e., the number of ions in an aggregate and the number of associations within that aggregate, we construct an adjacency matrix of the ions.

To determine a persistence/residence time (shown in the SI), we use a Heaviside step function between Na⁺ and O in TFSI⁻ to count associations (so H(r, t) = 1 if r ≥ r_cut and H(r, t) = 0 if r < r_cut) to determine an autocorrelation function 〈H(r, t)H(r, 0)〉, which can be used to extract a persistence/residence time.^{12,100,101,105}

4. Results

In this section, the results from the MD simulations are shown and used to validate the derived equations from the Theory Section. This is a self-contained section and does not require the Theory Section to understand the main messages found. We describe how the coordination numbers change with time after the step-change is applied, then the percolation criteria is analyzed and the cluster distribution is used to confirm the existence of ionic networks for the symmetric case. In the SI additional results are shown for the asymmetric case, further supporting the observations shown in the main text.

4.1. Symmetric case

In Fig. 2 we show how the coordination numbers between Na⁺ and TFSI⁻ change with time from the non-equilibrium MD simulations, with a comparison against the developed theory. Here we study the symmetric case of the SiIL NaTFSI_0.75EMIMTFSI_0.25, since the functionality (maximum number of associations it can form) of Na⁺ is f₊ = 4 and TFSI⁻ is f₋ = 3, this electrolyte mixture has p₊₋(t) = p₋₊(t) (i.e., the probability of cations associating to anions, p₊₋, becomes equal to the probability of anions associating to cations, p₋₊). The coordination numbers are related to the association probabilities through CN_ij = f_ip_ij.^60–62 We choose to show the results in terms of coordination numbers, as these are readily extracted from MD simulations, and it gives distinct numerical values for Na⁺ and TFSI⁻ in this symmetric case. The results for the simulations in Fig. 2 correspond to the charge rescaling case (see Molecular dynamics Section for further description). In the SI, we show results for the temperature rescaling non-equilibrium conditions. The same conclusions are drawn from those, but where smaller changes in coordination numbers are observed.

Fig. 2 Coordination numbers between Na⁺ and TFSI⁻ in NaTFSI_0.75EMIMTFSI_0.25 as a function of time, for the charge rescaling case, at the indicated temperatures, (a) 300 K, (b) 400 K, (c) 500 K, from MD simulations and theory.

Fig. 2a shows the simulations at 300 K. Initially, even though the charges of all atoms are essentially 0, we find coordination numbers of ∼0.3. This can be attributed to the fact that a real-space cutoff has been used, and the volume has been constrained to the equilibrium value at 300 K with the normal values of the charge. We find that CN₊₋ (number of TFSI⁻ anions bound to Na⁺ cations) and CN₋₊ (number of Na⁺ cations bound to TFSI⁻ anions) initially increases rapidly over the first 2 ps of the simulation, after which they increase less rapidly. After 15 ps, these coordination numbers have approximately saturated at their equilibrium values.

To compare against the theory, we have to extract the initial and final association probabilities. Then the only unknown in the theory is the time-constant of the decay in the probabilities, τ, since the rates of forming and breaking associations (k and f, respectively) are not known. Using the MD data, we fitted CN_+−/−+ over the 17 ps displayed to extract τ, which we report in Table 1. We find that the fitted τ for each curve (CN₊₋ and CN₋₊) gives the same decay rate of ∼0.4 ps⁻¹. This corresponds to a decay time of ∼2.5 ps, which agrees well with the time-scale observed in Fig. 2a. Using this fitted decay constant, we plot the theory prediction in Fig. 2a, where reasonable agreement is found. Moreover, using eqn (26), we use this information to extract the rate constant of forming associations, k, which we also report in Table 1. We find a value of 3 ps⁻¹ for the rate of forming associations, which is larger than the fitted inverse time constant and can be attributed to the rate of forming associations needing to be faster than the observed inverse rate constant, since they are reversible associations.

Table 1 Summary of inverse time scales (τ⁻¹), rate of forming an association (k), and association constant (λ) for the studied temperatures of the symmetric case

T/K	τ^−1/ps^−1	k/ps^−1	λ
300	0.401	3.000	71.178
400	0.418	2.398	42.383
500	0.511	2.231	24.880

---

In Fig. 2b and c, we show analogous non-equilibrium coordination numbers as a function of time for 400 K and 500 K, respectively. Again fitting these (using the initial and final association probabilities from MD) provides the time constants ∼0.42 and ∼0.51 ps⁻¹, for 400 K and 500 K, respectively; and we find that the fitted decay constant for the theory results is a reasonable fit to the MD simulations. With increasing temperature, we observe the inverse time constant to increase, indicating the time constant decreases with temperature. Intuitively, this makes sense, since at higher temperatures the associations will be weaker, which should therefore take less time to form. In Table 1, we also report the association constant (λ) at the studied temperatures, where we find it decreases with temperature. Again, using eqn (26) we extract the rate constant of forming an association, k, and find ∼2.4 and ∼2.2 ps⁻¹, for 400 K and 500 K, respectively. These are summarized in Table 1. With increasing temperature, we find k decreases. If the process was activated, ∝e^−βE_a, it should increase with temperature. This suggests that the activation barrier, if present, is not large, and other factors are dominating the rate of forming an association.

While the theory also has an initial rapid increase over ∼2 ps, and slower increase after, the quantitative agreement with the MD simulations is not exact (same conclusion is obtained from the temperature rescaling simulations, see SI). We find that the theory under-predicts the MD simulations at short-times, but over-predicts the values of the coordination numbers at long-times. The only variable that could cause this disagreement is the decay constant (τ), or more specifically the rate constants of forming (k) and breaking (f) an association. Since a single value of this rate of forming an association does not fit the data well, it suggests that there is more than one time scale in the simulations.

To investigate this further, we calculated the decay of the associations over time of, which can be used to extract residence/persistence times (see Molecular Dynamics Methods for the description). In other words, after an initial time in an equilibrium simulation, we counted the number of associations that remained after some given time. This is similar to the irreversible fragmentation solution discussed in the symmetric case. However, from the simulations, if those exact associations between a specific cation and anion reform, this causes the measure to increase, which would not be possible in the irreversible aggregation case. However, since species tend to diffuse and associate with other ions after dissociating, this measure generally decays with time. In the SI, we show this result for NaTFSI_0.75EMIMTFSI_0.25 at 300 K. What we find is that initially the associations decay by 10–20% very rapidly over several picoseconds. After which, the decay of the initial associations becomes more gradual, and only changes over a ∼ ns time scale. Similar results were found by Molinari et al. for SiILs,¹³ Self and Fong in the context of battery electrolytes,^100,101 and by Feng et al. in the context of ILs.¹⁰⁵ Therefore, there are clearly multiple time-scales for the kinetics of aggregate formation in SiILs, and the theory assuming there is one timescale is a limitation.

To test this further, we fitted the coordination numbers over time, but only for first 5 ps. In the SI we show these fits, where much better agreement is found, since it is described by a single, fast association-forming time scale (τ⁻¹ ≈ 1 ps⁻¹, which gives k ≈ 4.73 ps⁻¹). Moreover, we also show that the coordination numbers increase slightly from the values in Fig. 2, over 10 s of ps time scales. Therefore, in the MD simulations, the coordination numbers are initially controlled by the fast rate constant of forming associations of species that are already in proximity to each other, but at longer time scales this rate constant becomes limited by the diffusion/rearrangement of the electrolyte to allow for the equilibrium associations to be reached.^100,101 This explains why the theory underestimates the MD at short time scales, but overestimates the MD at long time scales. These observations pose a fundamental limitation of the theory, which assumes one rate constant, and motivates its development further. While the theory is not a perfect quantitative fit, it still provides insight into the underlying dynamics, and we further investigate it.

Thus far, we have systemically investigated the best possible agreement between the theory and simulations from fitting the decay constants. While we have found that to obtain a good quantitative fit, development of the theory to have multiple time scales is required, it is regardless unsatisfactory to have free fitting parameters in the theory. All that needs to be determined for a parameter-free comparison is the rate of forming/breaking an association k/f. Since we have found multiple time-scales, these rates need to be determined at short and long times. From the irreversible aggregation solution and the residence/persistence calculations, we can assume that an approximation to breaking an association (f) is given by the inverse residence/persistence time, which can be converted to the rate of forming an association through k = fλ. For short times, we find k_s ≈ 5.88 ps⁻¹ (given λ = 26.73), and at long-times we find k_l ≈ 0.03 ps⁻¹ (from λ = 163.02). At short-times, the extracted rate of forming associations is close to the fitted value, previously discussed. Moreover, in the SI, we show good agreement between the theory with this calculated rate (k_s) over short-times and the MD simulations. We also show the fit for the long-time extracted k_l, but much worse agreement is found overall, since most of the change occurs at short-times, and there is a memory effect in these non-equilibrium conditions. The long-time behavior is improved, however, with this extracted value k_l. Therefore, overall, it appears that with development of the rate-kernels, a quantitative agreement without free fitting parameters is within reach.

Using the time evolving association probabilities, we can determine the time a percolating ionic network, i.e., a gel, between Na⁺ and TFSI⁻ forms through tracking the product of the association probabilities, p₊₋p₋₊. When this value reaches [(f₊ − 1)(f₋ − 1)]⁻¹, that is the critical criteria for bond percolation on a Bethe lattice, i.e., a Cayley tree.^57,61 In SiILs, the percolating ionic network is know to occur at x ≈ 0.25 from simulations;⁶⁰ while the theory often predicts it at smaller mole fractions at equilibrium, owing to the formation of loops in the clusters.^15,60

In Fig. 3 we show p₊₋p₋₊ as a function of time for NaTFSI_0.75EMIMTFSI_0.25 at 300 K. From the MD simulations, we find that the critical criteria to form a gel is already satisfied after ∼0.5 ps. The theory predicts the percolating ionic network to form at ∼1 ps, but this is because the theory is underestimating the association probabilities at short times from the fitted decay constants.

Fig. 3 Gelation criteria in NaTFSI_0.75EMIMTFSI_0.25 as a function of time, for the charge rescaling case, at 300 K, for MD and theory. The percolation point is given by p₊₋p₋₊ = [(f₊ − 1)(f₋ − 1)]⁻¹.

To check if the critical gel criteria is actually being satisfied in the simulations, we can inspect the cluster distribution.^57,60,61 If we are in a pre-gel regime, we only expect small clusters (relative to the simulation size) to occur. As we approach the gel-point, the concentration of larger clusters continuously grows. When the gel-point is reached, there is simultaneously aggregates of all sizes (within the simulation size) and an aggregate that percolates throughout the simulation cell. In the post-gel regime, larger aggregates bind to the gel, increasing the gels volume fraction, leaving only small clusters. See ref. 57, 60 and 61 for more information.

In Fig. 4, we show the cluster distribution computed from MD simulations and theory at 3 times. At t = 0, we find that free anions and cations dominate, with a smaller number of ion pairs, and even fewer aggregates containing a handful of ions. At t = 0.4 ps, right at the critical gel time, we find that free anions still dominate the cluster distribution (as they usually do for SiILs^15,60), but there are aggregates which exist with over 30 ions. At this time, we also find a single aggregate with a significant portion of the ions in the simulation cell. At later times than t = 0.4 ps, the size of the aggregates decreases steadily. We show an example at t = 4 ps, where free anions still dominate, and there are at most 5 ions in an aggregate, which is smaller than when the charges were scaled to 0 at t = 0 ps. Therefore, the critical gel-time of the simulation predicted in Fig. 3 corresponds well to the behavior of the aggregates in the simulation, shown in Fig. 4, confirming the bond percolation criteria.^57,60,61 The theory predictions (calculated with the cluster distribution, and the association probabilities and association constant changing in time) match the simulated ones qualitatively very well. Previously, for SiILs, as discussed in ref. 60, the theory cluster distribution did not agree well with that of the MD simulations, owing to the presence of loops in the aggregates.¹⁵ In the SI, we computed the cluster bond density, which quantifies how many loops exist in the aggregates. Unlike the equilibrium case,¹⁵ we only find a few loops, which could be why good agreement is obtained here.

Fig. 4 Cluster distributions c_lm, with m TFSI⁻ anions and l Na⁺ cations in an aggregate, at various times in the MD simulation and theory, as indicated, at 300 K.

5. Discussion

To summarize, we developed the kinetics of aggregation formalism of Sciortino, Tartaglia and co-workers^{79–81,86–90} from patchy particle systems to work for concentrated electrolytes, extending the equilibrium formalism developed by McEldrew et al.^{26,57,59–62} Comparisons to molecular dynamics simulations found good qualitative agreement, but quantitative differences were observed. These were investigated further, and found to result from multiple time scales existing in the concentrated electrolyte.^13,100,101 Specifically, there is a fast time at short times, corresponding to associations forming between species that are already in close proximity. After which the time scale slows, being more limited by diffusion and rearrangements of ions.¹⁰⁰ While the theory is not quantitative, it does provide analytical expressions with insight into the overall mechanisms occurring and a framework to understand them. For example, we also demonstrated that the quasi-equilibrium solution applies to SiILs, where the cluster distribution changes in time from the changes in association probabilities, with this quasi-equilibrium solution connecting time and temperature through the association probabilities, since it can be thought that p(t) is changing through equilibrium states of p(T).^{79–81,86–90}

The analysis here has been completely theoretical and computational based, but this does not limit the applicability of the developed formalism. In fact, as analytical expression for coordination numbers changing in time after a step-change in conditions were derived, this directly lends itself to facile experimental and computational comparison. In battery electrolytes, Raman and IR spectroscopic methods are commonly employed to extract information about the coordination environments,²⁶ since there are often distinct features in the spectra for bound/unbound functional groups. If a step-change in temperature can be introduced into these electrolytes, and the response experimentally characterized, the developed equations can be used to obtain the rates of forming/breaking associations in these electrolytes directly from experiments. In addition, for computational work, the formalism presented here provides a framework to understand equilibrium and non-equilibrium changes in coordination environments and clusters, which can also be linked to other physiochemical properties.⁶⁰

Using the developed formalism and performed MD simulations, we extracted the rate for the formation of an association, which was found to change with time (from the short-time to long-time response) and, vary slightly with temperature and composition of the electrolyte (see SI). We found that the short-time rate of forming associations is of the order of 5 ps⁻¹, which corresponds to a typical time for ions to rattle in the cages that form in concentrated electrolytes.¹⁰⁵ Whereas, the long time scale corresponds to 10s of ps, which is closer to the scale of diffusive time scales in concentrated electrolytes. Therefore, the magnitudes of the rates appear to be physically sensible.^12,13 As discussed in detail in ref. 90 and 99, the rates of fragmentation (F) and coalescence (K) depend on the time-scale for an association to form and the diffusive time-scales of the species which associate, which we also found here through the multiple time scales in the electrolytes. This transition from “bond” forming time-scale to diffusion time-scale can be modeled with rate kernels which depend on the association probabilities, allowing for transitions from rates limited by “bond” forming to diffusion.⁸¹ These observations are interesting as multiple time scales in simulations have already been observed,^12,13 but if we have a distribution of clusters, it might be expected that there are a distribution of relaxation times. In the formalism developed here, however, this does not emerge, but one (or two) time scales emerge from the theory as this cluster distribution is accounted for in an ensemble averaged way through the association probabilities which change with time.

In the context of the temperature dependence, an exponential law would be expected for an activated process, meaning that the rate of forming associations would increase with increasing temperature. Since we do not observe this, it suggests this process is not an activated one. While in the context of different compositions (see SI), changes in dielectric environment and viscosity might also be expected to change these rates substantially. However, we also do not observe strong effects from these variables. Both of these observations could be a result of the short-time scale dominating the effect, where there is limited influence from temperature/composition on forming association between ions that are already in proximity. It is known that residence times of different SiILs,¹³ ILs,¹⁰⁵ WiSEs⁶¹ and battery electrolytes¹⁰⁰ can be quite different, which suggests different rates of forming associations, but this is likely dominated by the diffusive time-scale. Future work should investigate these observations in more depth.

While we have mainly extracted the rate constants from fitting the MD simulations, it was also demonstrated that approximate rate values can be obtained from bulk, equilibrium calculations of residence/persistence. The obtained rates at short-times agreed well with the MD simulation results at short times, demonstrating a parameter-free model could be obtained. However, the limitation of the multiple time scales prevented a full comparison here, with developments into the rate kernels changing throughout the processes being a promising direction for a parameter-free quantitative match.⁸¹

Another source of error between the theory and simulations could be a result of the employed force field. A simple non-polarizable force-field was utilized here, with scaled charges to model the effects of using a polarizable force-field.^12,13 If another force-field is utilized for the simulations, the non-equilibrium dynamics would almost certainly change, and this could change the agreement between the theory and simulations observed here. Moreover, comparisons to ab initio MD and machine learning interatomic potential simulations could also be readily achieved through step-changing temperature. It is not a priori clear, however, if changing the force-field will improve the agreement with the theory. It is clear, however, that the simulated time and length scales will limit the resulting rates of forming/breaking associations, while experimental measurements are macroscopic so the finite size effects should not play a role in bulk measurements.

Here we have investigated correlations between a select cation and anion (+−), using a real-space cutoff to define associations between species. This formalism is more general, however. It can be applied to understand like-charge correlations too, albeit being much weaker. In the SI, we show how the “coordination numbers” between like charges (++ and −−) change with time. We find that the ++ correlations only have a fast, sub-diffusive time scale and the theory is able to reproduce these dynamics well. However, the −− correlations have a non-monotonic dependence, initially reducing before returning to slightly above the initial value, and the theory is unable to describe its behaviour. Interestingly, this suggests that the anion plays a more major role in the multiple time scales of the electrolyte than the alkali metal cations. Beyond this, the requirement of associations between species (cation–anion) can be relaxed, so taking into account bidentate associations is possible, or higher-order encapsulating associations,¹⁴ for example. Additionally, other methods for defining the associations, such as kinetic-based cutoffs,¹⁰⁵ could be used to investigate the non-equilibrium effects studied here.

Moreover, from the electrolytes that have studied with this formalism so far, SiILs is one of the most challenging.^15,60 Its highly concentrated and correlated nature pushes the boundaries of the theory, where loops within the aggregates form (breaking the Cayley tree approximation), making it more difficult to apply the theory.^15,60 While SiILs are the most challenging example that have been applied the equilibrium theory to so far, we nonetheless found that the non-equilibrium formalism could capture well the dynamics of SiILs. In future work, applied to other electrolytes, we expect the developed formalism to perform even better.

It is hoped that this provides the first link between a microscopic theory of ionic associations and non-equilibrium changes to electrolyte structure, which is further built upon. In what follows, we outline possible extensions and connections to other areas in concentrated electrolytes, highlighting where experimental verification can enter. While this is not an exhaustive list, we hope it provides inspiration for what might now be possible.

In the context of ionic transport in concentrated electrolytes, there is simultaneously vehicular motion of the aggregates/clusters and structural diffusion.^100,101 To determine the dominant mechanism, Self and Fong et al.^100,101 suggested comparing length scales of the aggregates (L^s) to the length scales over which they move in a vehicular way. The later is determined from , where D_i is the self diffusion coefficient and τ_res is the residence time of the association. If L^c > L^s, there is vehicular dominated motion, while if L^c < L^s it is structural dominated motion. There are clearly multiple time scales for associations in electrolytes, as found here for SiILs, and by others.^{13,61,100,101,105} Self and Fong et al.^100,101 suggested using the long-time scale, and ignoring the short, sub-diffusive changes in associations. If the short time scale is used, which is attributed to the binding/breaking of associations in proximity, this would suggests structural diffusion dominates more. While if the long time scale is used, it suggest more of a vehicular motion could be dominant. The results here suggest that this short-time scale is important, indicating that transport in these electrolytes could be more dominated by structural diffusion. Using the formalism developed here, it might be possible to more quantitatively link vehicular and structural diffusion modes.

As discussed in the Results and Theory Sections, concentrated electrolytes have been predicted to form a percolating ionic network, i.e., a gel phase^57,59–61 (to date, only one experimental claim of a gel from associations has been reported¹⁰⁶). The terminology of the gel phase comes from the inherited language from polymer systems.^67–78 In thermoreversible polymers the time scales of the associations are much longer than electrolytes, since the associations (“bonds”) are much stronger. Upon gelation, this gives rise to a divergence in viscosity and the onset of an elastic response.⁵⁷ It is not clear if the elastic response will be measurable from the shorter time scales of the associations in concentrated electrolytes, however. In the context of experimental measurements, non-Newtonian behavior of WiSEs has been observed,⁵⁰ suggesting that in WiSE these associations could be long-lived enough. Significantly more work is required in this area, from theory, simulations and experiments, to understand how similar electrolytes are to polymers in this regard, i.e., how good the analogy between ionic associations and thermoreversible bonds. Using the developed formalism, it should be possible to include the finite time scales of associations, to more quantitatively predict the frequency dependence of the viscoelastic response of the fluid.^91,92

Here we have only studied bulk electrolytes, but the technologies which use these electrolytes often require contact with electrified interfaces.^9,23,24,107 How concentrated electrolytes accumulate at charged interfaces, i.e., the electrical double layer (EDL), has been studied at equilibrium in the context of McEldrew et al.’s formalism for ILs,^63,64 SiILs^15,66 and WiSEs.⁶⁵ Using the formalism developed here, the charging dynamics of the EDL could be studied with the relaxation of the aggregates. Experimentally, this could be achieved with time-resolved SIERAS or WAXS,³⁸ which would also allow for the changes in coordination environments to be probed. While the aggregation relaxation time scale appears to be faster than the electrolyte transport, studying the relaxation of aggregates in confinement could result in interesting interplay with the transport of ions. Building on this, the desolvation dynamics of battery electrolytes at electrodes would be another area to investigate, to more quantitatively determine the time scales of how the solvation shell is striped from the active cation as it intercalates into the electrode.⁹⁴ As a first step, the removal of one solvent from the diffuse EDL to Helmholtz layer could be investigated,¹⁰⁸ building all the way up to the complete removal and intercalation of the active ion.

Finally, in the context of “anomalous underscreening”,^36–42,109 it has been shown that with slower approaches of the interfaces the long-ranged interactions vanish.⁴³ Therefore, this effect is likely not from bulk electrostatic interactions, but a large part is hydrodynamic interactions.^15,43 It remains to be seen if it can be completely described by classical fluid mechanics, however. In confinement with these concentrated electrolytes, ionic aggregation effects could also play a role.¹⁵ In the context of SiILs, Zhang et al.¹⁵ found force separation, extremely long force decay lengths, and changes in refractive indices of the electrolyte between mica. This made them conclude there was confinement induced changes in electrolyte composition, possibly leading to the formation of a gel.¹⁵ Further developing the formalism here to confined geometries, also with changing composition, could provide insight into the role of ionic aggregation kinetic in these measurements. We have found that the time scale of ionic aggregation τ can diverge if the association probabilities approach 1, which could give rise to extremely long relaxation times. However, the short-time scale of association appears to be of the order of ∼ps (which could also arise from the finite size of the simulation cells). Further work should explore if aggregate relaxation could help explain these puzzling non-equilibrium measurements.

Author contributions

Z. A. H. G. conceptualized, curated, analysed, wrote and edited all parts of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this paper, including the coordination numbers as a function of time, adjacency matrices, and cluster distributions at selected times are available on the zenodo repository at https://doi.org/10.5281/zenodo.20207432. The analysis scripts for determining the coordination numbers, cluster distributions and cluster bond density can be found at https://github.com/zg1413/kinetics_aggregation.git.

Supplementary information (SI): additional derivations and results. See DOI: https://doi.org/10.1039/d6sc01986b.

Acknowledgements

I am extremely greatful for the support and stimulating discussion from M. McEldrew, A. Kornyshev, C. Phelan and R. Espinosa-Marzal. I also acknowledge discussions with J. P. de Souza, R. Weatherup, N. Molinari, M. Markiewitz and M. Bazant. Z.A.H.G acknowledges support through the Glasstone Research Fellowship in Materials and The Queen's College, University of Oxford.

References

1. L. Suo, Y.-S. Hu, H. Li, M. Armand and L. Chen, Nat. Commun., 2013, 4, 1481.
2. L. Suo, O. Borodin, T. Gao, M. Olguin, J. Ho, X. Fan, C. Luo, C. Wang and K. Xu, Science, 2015, 350, 938.
3. J. Wang, Y. Yamada, K. Sodeyama, C. H. Chiang, Y. Tateyama and A. Yamada, Nat. Commun., 2016, 7, 12032.
4. F. Wang, O. Borodin, T. Gao, X. Fan, W. Sun, F. Han, A. Faraone, J. A. Dura, K. Xu and C. Wang, Nat. Mater., 2018, 17, 543.
5. Y. Yamada, K. Usui, K. Sodeyama, S. Ko, Y. Tateyama and A. Yamada, Nat. Energy, 2016, 1, 16129.
6. L. Suo, O. Borodin, Y. Wang, X. Rong, W. Sun, X. Fan, S. Xu, M. A. Schroeder, A. V. Cresce and F. Wang, et al., Adv. Energy Mater., 2017, 7, 1701189.
7. Q. Dou, S. Lei, D.-W. Wang, Q. Zhang, D. Xiao, H. Guo, A. Wang, H. Yang, Y. Li, S. Shi and X. Yan, Energy Environ. Sci., 2018, 11, 3212.
8. A. A. Kornyshev, J. Phys. Chem. B, 2007, 111, 5545.
9. M. V. Fedorov and A. A. Kornyshev, Chem. Rev., 2014, 114, 2978.
10. T. Welton, Chem. Rev., 1999, 99, 2071–2084.
11. J. P. Hallett and T. Welton, Chem. Rev., 2011, 111, 3508–3576.
12. N. Molinari, J. P. Mailoa and B. Kozinsky, J. Phys. Chem. Lett., 2019, 10, 2313.
13. N. Molinari, J. P. Mailoa, N. Craig, J. Christensen and B. Kozinsky, J. Power Sources, 2019, 428, 27.
14. N. Molinari and B. Kozinsky, J. Phys. Chem. B, 2020, 124, 2676.
15. X. Zhang, Z. A. Goodwin, A. G. Hoane, A. Deptula, D. M. Markiewitz, N. Molinari, Q. Zheng, H. Li, M. McEldrew and B. Kozinsky, et al., ACS Nano, 2024, 18, 34007.
16. C. M. Efaw, Q. Wu, N. Gao, Y. Zhang, H. Zhu, K. Gering, M. F. Hurley, H. Xiong, E. Hu and X. Cao, et al., Nat. Mater., 2023, 22, 1531.
17. K. Xu, Chem. Rev., 2014, 114, 11503.
18. S. C. Kim, J. Wang, R. Xu, P. Zhang, Y. Chen, Z. Huang, Y. Yang, Z. Yu, S. T. Oyakhire, W. Zhang, L. C. Greenburg, M. S. Kim, D. T. Boyle, P. Sayavong, Y. Ye, J. Qin, Z. Bao and Y. Cui, Nat. Energy, 2023, 8, 814.
19. Y. Xie, J. Wang, B. H. Savitzky, Z. Chen, Y. Wang, S. Betzler, K. Bustillo, K. Persson, Y. Cui, L.-W. Wang, C. Ophus, P. Ercius and H. Zheng, Sci. Adv., 2023, 9, eadc9721.
20. K. Xu, Chem. Rev., 2004, 104, 4303.
21. J. B. Goodenough and K. S. Park, The Li-ion rechargeable battery: A perspective, 2013.
22. M. Li, C. Wang, Z. Chen, K. Xu and J. Lu, Chem. Rev., 2020, 120, 6783–6819.
23. M. A. Gebbie, B. Liu, W. Guo, S. R. Anderson and S. G. Johnstone, ACS Catal., 2023, 13, 16222.
24. J. Zheng, J. A. Lochala, A. Kwok, Z. D. Deng and J. Xiao, Adv. Sci., 2017, 4, 1700032.
25. C. M. Phelan, E. Björklund, J. Singh, M. Fraser, P. N. Didwal, G. J. Rees, Z. Ruff, P. Ferrer, D. C. Grinter and C. P. Grey, et al., Chem. Mater., 2024, 36, 3334.
26. C. Phelan, A. Bhandari, J. Singh, M. F. nad Erik Björklund, J. Swallow, P. Bencok, C. Grey, Z. Goodwin, C. Skylaris, and R. Weatherup, Chemrxiv, 2025,  DOI:10.26434/chemrxiv-2025-k903g.
27. Z. Yu, L. A. Curtiss, R. E. Winans, Y. Zhang, T. Li and L. Cheng, J. Phys. Chem. Lett., 2020, 11, 1276.
28. O. Borodin, L. Suo, M. Gobet, X. Ren, F. Wang, A. Faraone, J. Peng, M. Olguin, M. Schroeder and M. S. Ding, et al., ACS Nano, 2017, 11, 10462.
29. P. Lannelongue, R. Bouchal, E. Mourad, C. Bodin, M. Olarte, S. le Vot, F. Favier and O. Fontaine, J. Electrochem. Soc., 2018, 165, A657.
30. J. Vatamanu and O. Borodin, J. Phys. Chem. Lett., 2017, 8, 4362.
31. L. G. Chagas, S. Jeong, I. Hasa and S. Passerini, ACS Appl. Mater. Interfaces, 2019, 11, 22278.
32. T. C. Lourenco, L. G. Dias and J. L. Da Silva, ACS Appl. Energy Mater., 2021, 4, 4444.
33. R.-S. Kühnel, D. Reber and C. Battaglia, ACS Energy Lett., 2017, 2, 2005.
34. O. Borodin, J. Self, K. A. Persson, C. Wang and K. Xu, Joule, 2020, 4, 69.
35. Z. Yu, N. P. Balsara, O. Borodin, A. A. Gewirth, N. T. Hahn, E. J. Maginn, K. A. Persson, V. Srinivasan, M. F. Toney, K. Xu, K. R. Zavadil, L. A. Curtiss and L. Cheng, ACS Energy Lett., 2022, 7, 461.
36. M. A. Gebbie, M. Valtiner, X. Banquy, E. T. Fox, W. A. Henderson and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 9674.
37. M. A. Gebbie, H. A. Dobes, M. Valtiner and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 7432–7437.
38. M. Han, H. Kim, C. Leal, M. Negrito, J. D. Batteas and R. M. Espinosa-Marzal, Adv. Mater. Interfaces, 2020, 7, 2001313.
39. M. Han, R. Zhang, A. A. Gewirth and R. M. Espinosa-Marzal, Nano Lett., 2021, 21, 2304.
40. T. S. Groves, C. S. Perez-Martinez, R. Lhermerout and S. Perkin, J. Phys. Chem. Lett., 2021, 12, 1702.
41. A. Lee, C. S. Perez-Martinez, A. M. Smith and S. Perkin, et al., Faraday Discuss., 2017, 199, 239.
42. A. M. Smith, A. A. Lee and S. Perkin, J. Phys. Chem. Lett., 2016, 7, 2157.
43. B. Cross, L. Garcia, E. Charlaix and P. Kekicheff, Proc. Natl. Acad. Sci. U. S. A., 2026, 123, e2517939123.
44. A. Härtel, M. Bültmann and F. Coupette, Phys. Rev. Lett., 2023, 130, 108202.
45. R. Kjellander, Phys. Chem. Chem. Phys., 2016, 18, 18985.
46. Y. Avni, R. M. Adar and D. Andelman, Phys. Rev. E, 2020, 101, 010601.
47. R. M. Adar, S. A. Safran, H. Diamant and D. Andelman, Phys. Rev. E, 2019, 100, 042615.
48. J. P. de Souza, Z. A. Goodwin, M. McEldrew, A. A. Kornyshev and M. Z. Bazant, Phys. Rev. Lett., 2020, 125, 116001.
49. E. Krucker-Velasquez and J. W. Swan, J. Chem. Phys., 2021, 155, 134903.
50. T. Yamaguchi, A. Dukhin, Y.-J. Ryu, D. Zhang, O. Borodin, M. A. Gonzàlez, O. Yamamuro, D. L. Price and M.-L. Saboungi, J. Phys. Chem. Lett., 2024, 15, 76.
51. J. Lim, K. Park, H. Lee, J. Kim, K. Kwak and M. Cho, J. Am. Chem. Soc., 2018, 140, 15661.
52. Y. S. Meng, V. Srinivasan and K. Xu, Science, 2022, 378, 1065.
53. O. Borodin and D. Bedrov, J. Phys. Chem. C, 2014, 118, 18362.
54. O. Borodin, X. Ren, J. Vatamanu, A. von Wald Cresce, J. Knap and K. Xu, Acc. Chem. Res., 2017, 50, 2886.
55. Q. Wu, M. T. McDowell and Y. Qi, JACS, 2023, 145, 2473.
56. A. V. Tkachenko, C. Cao, A. C. Marschilok, and D. Lu, arXiv, 2025, preprint, arXiv:2512.18167,  DOI:10.48550/arXiv.2512.18167.
57. M. McEldrew, Z. A. Goodwin, S. Bi, M. Z. Bazant and A. A. Kornyshev, J. Chem. Phys., 2020, 152, 234506.
58. M. McEldrew, Ion Aggregation, Correlated Ion Transport and the Double Layer in Super-Concentrated Electrolytes, Ph.D. thesis, MIT, 2021.
59. M. McEldrew, Z. A. H. Goodwin, H. Zhao, M. Z. Bazant and A. A. Kornyshev, J. Phys. Chem. B, 2021, 125, 2677–2689.
60. M. McEldrew, Z. A. H. Goodwin, N. Molinari, B. Kozinsky, A. A. Kornyshev and M. Z. Bazant, J. Phys. Chem. B, 2021, 125, 13752.
61. M. McEldrew, Z. A. H. Goodwin, S. Bi, A. Kornyshev and M. Z. Bazant, J. Electrochem. Soc., 2021, 168, 050514.
62. Z. A. H. Goodwin, M. McEldrew, B. Kozinsky and M. Z. Bazant, PRX Energy, 2023, 2, 013007.
63. Z. A. H. Goodwin, M. McEldrew, J. P. de Souza, M. Z. Bazant and A. A. Kornyshev, J. Chem. Phys., 2022, 157, 094106.
64. Z. A. H. Goodwin and A. A. Kornyshev, Electrochim. Acta, 2022, 434, 141163.
65. D. M. Markiewitz, Z. A. H. Goodwin, Q. Zheng, M. McEldrew, R. M. Espinosa-Marzal and M. Z. Bazant, ACS Appl. Mater. Interfaces, 2025, 17, 29515.
66. D. M. Markiewitz, Z. A. H. Goodwin, M. McEldrew, J. P. de Souza, X. Zhang, R. M. Espinosa-Marzal and M. Z. Bazant, Faraday Discuss., 2024, 365.
67. P. J. Flory, J. Chem. Phys., 1942, 10, 51.
68. W. H. Stockmayer, J. Chem. Phys., 1943, 11, 45.
69. W. H. Stockmayer, J. Chem. Phys., 1944, 12, 125.
70. W. H. Stockmayer, J. Polym. Sci., 1952, 9, 69.
71. M. Ishida and F. Tanaka, Macromolecules, 1997, 30, 3900.
72. F. Tanaka, Macromolecules, 1989, 22, 1988.
73. F. Tanaka, Macromolecules, 1990, 23, 3784.
74. F. Tanaka and W. H. Stockmayer, Macromolecules, 1994, 27, 3943.
75. F. Tanaka and M. Ishida, J. Chem. Soc. Faraday Trans., 1995, 91, 2663.
76. F. Tanaka, Phys. A, 1998, 257, 245.
77. F. Tanaka and M. Ishida, Macromolecules, 1999, 32, 1271.
78. F. Tanaka, Polym. J., 2002, 34, 479.
79. F. Sciortino and E. Zaccarelli, Curr. Opin. Solid State Mater. Sci., 2011, 15, 246.
80. F. Sciortino, C. D. Michele, S. Corezzi, J. Russo, E. Zaccarelli and P. Tartaglia, Soft Matter, 2009, 5, 2571.
81. S. Corezzi, D. Fiorettoac and F. Sciortino, Soft Matter, 2012, 8, 11207.
82. A. Gröschel, A. Walther, T. I. Löbling, F. H. Schacher, H. Schmalz and A. H. E. Müller, Nature, 2013, 503, 247.
83. A. Umar and S.-M. Choi, J. Phys. Chem. C, 2013, 117, 11738.
84. A. Zaccone, J. J. Crassous and M. Ballauff, J. Chem. Phys., 2013, 138, 104908.
85. A. Zaccone, J. J. Crassous, B. Béri and M. Ballauff, Phys. Rev. Lett., 2011, 107, 168303.
86. S. Corezzi, D. Fioretto, C. D. Michele, E. Zaccarelli and F. Sciortino, J. Phys. Chem. B, 2010, 114, 3769.
87. E. Bianchi, P. Tartaglia, E. L. Nave and F. Sciortino, J. Phys. Chem. B, 2007, 111, 11765.
88. S. Corezzi, C. D. Michele, E. Zaccarelli, P. Tartaglia and F. Sciortino, J. Phys. Chem. B, 2009, 113, 1233.
89. F. Sciortino, C. D. Michele and J. F. Douglas, J. Phys.: Condens. Matter, 2008, 20, 155101.
90. P. Tartaglia, J. Phys.: Condens. Matter, 2009, 21, 504109.
91. F. Tanaka, Macromolecules, 2024, 57, 10600.
92. F. Tanaka, Gels, 2023, 9, 379.
93. K. Xu and A. von Wald Cresce, J. Mater. Res., 2012, 2327.
94. I. Rubinstein, M. Bazant, N. Kononenko, V. Nikonenko and B. Zaltzman, Phys. Rev. E, 2026, 113, 045404.
95. Z. Tian, Y. Zou, G. Liu, Y. Wang, J. Yin, J. Ming and H. N. Alshareef, Adv. Sci., 2022, 9, 2201207.
96. M. Z. Bazant, B. D. Storey and A. A. Kornyshev, Phys. Rev. Lett., 2011, 106, 046102.
97. T. Welton, Chem. Sci., 2025, 16, 18976.
98. P. J. Flory, Principles of polymer chemistry, Cornell University Press, 1953.
99. P. G. J. van Dongen and M. H. Ernst, J. Stat. Phys., 1984, 37, 301.
100. J. Self, K. D. Fong and K. A. Persson, ACS Energy Lett., 2019, 4, 2843.
101. K. D. Fong, J. Self, K. M. Diederichsen, B. M. Wood, B. D. McCloskey and K. A. Persson, ACS Cent. Sci., 2019, 5, 1250–1260.
102. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comp. Phys. Comm., 2022, 271, 10817.
103. J. N. C. Lopes and A. A. Pádua, Theor. Chem. Acc., 2012, 131, 1129.
104. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157.
105. G. Feng, M. Chen, S. Bi, Z. A. Goodwin, E. B. Postnikov, N. Brilliantov, M. Urbakh and A. A. Kornyshev, Phys. Rev.: X, 2019, 9, 021024.
106. D. Reber, R. Grissa, M. Becker, R.-S. Kühnel and C. Battaglia, Adv. Energy Mater., 2021, 11, 2002913.
107. H. Cheng, Q. Sun, L. Li, Y. Zou, Y. Wang, T. Cai, F. Zhao, G. Liu, Z. Ma, W. Wahyudi, Q. Li and J. Ming, ACS Energy Lett., 2022, 7, 490.
108. Z. A. Goodwin, D. M. Markiewitz, Q. Wu, Y. Qi and M. Z. Bazant, ACS Appl. Energy Mater., 2025, 8, 8376.
109. R. M. Espinosa-Marzal, Z. A. Goodwin, X. Zhang, and Q. Zheng, in One Hundred Years of Colloid Symposia: Looking Back and Looking Forward, ACS Publications, 2023, pp. 123–148.

---